Fig. 1

Overview of the experimental design and analysis pipeline. a Summary of the analysis pipeline. b The 16-channel EEG configuration, with Fz as the ground electrode and A2 as the reference electrode. c Diagram illustrating the screen layout during the bimanual visuomotor tracking task; arrow colors represent the direction of the joysticks and cursors, \(c\) denotes cursor position, and \(g\) denotes track position. d Timing of the experimental trial

Participants

Twenty participants (18 males, 2 females; \(M_ = 24.7\), \(SD_ = 2.3\)) without any known neurological disorders were recruited for the study. All participants had normal or corrected-to-normal vision and were confirmed as right-handed using the Edinburgh Handedness Inventory (Oldfield 1971).

The experimental protocol was approved by the Ethics Committee of the Graduate School of Engineering, Tohoku University (Approval Number: 24A-5). Prior to participation, all individuals provided informed consent following the Declaration of Helsinki. Participants were informed that identifiable data and images may be included in publications, but these would be anonymized or pseudonymized to ensure privacy and confidentiality.

EEG signal acquisition and preprocessing

The overview of the signal processing and analysis pipeline is summarized in Fig. 1a. EEG signals were acquired and amplified using a 16-channel g.USBamp amplifier (g.tec Medical Engineering GMBH) at a sampling rate of 512 Hz. Wet active Ag/AgCl electrodes were placed at AF3, AF4, FC3, FCz, FC4, C3, Cz, C4, T7, T8, CP3, CPz, CP4, Pz, O1, and O2. The ground electrode was located at Fz, while the reference electrode was attached to the right earlobe (A2), as depicted in Fig. 1b.

An eighth-order Butterworth band-pass filter, with cutoff frequencies set between 0.5 and 100 Hz, was applied to the EEG signals. Furthermore, a fourth-order Butterworth notch filter was employed to reduce power line interference at 50 Hz. Ocular, muscle, and movement artifacts were identified through visual inspection and subsequently removed using independent component analysis. Additionally, the signals were re-referenced using a surface Laplacian filter to enhance localized activity and suppress diffuse signals.

Experimental design

Participants were positioned approximately 1 m from a display screen, which featured two tracks on the left and right sides for bimanual tracking using joysticks. The goal of the experiment was to assess participants’ adaptability in bimanual visuomotor tracking tasks. Two joystick-to-cursor mappings were used: Up-down, in which joystick movements were rotated 90 degrees counterclockwise relative to cursor movements, and Inverted, where joystick movements were mirrored relative to cursor movements, as shown in Fig. 1c.

At the start of each trial, the joystick cursors appeared in red, positioned at the center of each track. After a 1-second delay, the cursors changed to yellow, indicating the preparatory phase, and after another second, they turned green, signaling the commencement of the trial (GO signal). Each trial lasted 60 s, during which the tracks moved downward. Participants were required to align the cursors with the tracks. Any deviation or contact with the track boundaries caused the cursors to turn red, prompting participants to realign their tracking, as shown in Fig. 1d.

The experiment consisted of three sessions: pre-training, training, and post-training. In the pre-training and post-training sessions, participants completed four trials–two for each task presented in a random sequence–at the hardest difficulty level. In the training session, six trials were conducted, with three trials for each task in a random sequence. The difficulty gradually increased across these trials, corresponding to easy, moderate, and hard levels. For each difficulty level, the frequency parameter \(f\) was drawn from different ranges. Specifically, for the easy level, \(f \in \mathcal (0.01, 0.5)\); for the moderate level, \(f \in \mathcal (0.01, 0.75)\); and for the hard level, \(f \in \mathcal (0.01, 1)\).

The tracks used in the task were generated pseudorandomly by combining 20 sine waves, as described by the equation:

$$\begin g = \sum _^ A_i \sin \left( f_i \left( y + \frac \right) + \phi _i \right) \end$$

(1)

where \(g\) is the x-coordinate of the track, \(y\) is the y-coordinate, and \(t\) represents time. The parameters \(A_i \sim \mathcal (1, 10)\) are amplitude parameters drawn from a uniform distribution, and \(\phi _i \sim \mathcal (-180, 180)\) are phase parameters also drawn from a uniform distribution. Before starting the experimental trials, participants completed several familiarization trials. A brief break was provided after each trial to minimize fatigue. Behavioral data were sampled at 50 Hz, in line with the monitor’s refresh rate.

Performance was quantified using the Pearson correlation coefficient (\(r\)) to assess the linear relationship between the cursor positions (\(c\)) and the track positions (\(g\)). The formula for \(r\) is given by:

$$\begin r = \frac^n (c_i - \bar)(g_i - \bar)}^n (c_i - \bar)^2} \sqrt^n (g_i - \bar)^2}}, & \text r \ge 0 \\ 0, & \text r < 0 \end\right. } \end$$

(2)

where \(c_i\) and \(g_i\) represent the \(i\)-th data points for the cursor and track positions, respectively, and \(\bar\) and \(\bar\) are their mean values. \(n\) denotes the total number of data points. This constraint ensures the value of \(r\) lies between 0 and 1, where values closer to 1 indicate better performance.

The Pearson correlation coefficient was chosen over metrics such as absolute error (AE) or root mean square error (RMSE) because it captures the strength of the linear relationship between cursor and track positions, regardless of scale or offset. This is well-suited for tasks where participants are allowed some deviation and aim to stay within track boundaries rather than precisely follow the centerline. Unlike AE and RMSE, which are sensitive to outliers, \(r\) emphasizes co-variation and is less affected by occasional large deviations, making it a more robust indicator of overall trajectory tracking performance.

Performance for each session was calculated by averaging the values of all trials within a session, irrespective of task conditions, for each participant. Moreover, principal component analysis (PCA) was employed to reduce left- and right-hand performance measures into a single dimension, referred to as bimanual performance. This approach is favored over simply averaging the two performances, as PCA maximizes the variance across both hands, accounting for both individual and combined contributions. Importantly, the performance of the left and right hands may differ significantly; thus, averaging the two values could obscure meaningful asymmetries and lead to a loss of critical information regarding hand-specific performance characteristics.

Functional connectivity analysisFunctional connectivity estimation

Prior to functional connectivity estimation, each EEG trial was segmented into epochs using sliding windows. A window length of 4 s and a stride of 2 s were applied, resulting in 120 windows per trial. These windows were subsequently used for functional connectivity analysis.

Brain functional connectivity was estimated using the corrected imaginary phase-locking value (ciPLV) (Bruña et al. 2018), an extended version of the original phase-locking value (PLV) (Lachaux et al. 1999). This method was selected for its robustness against volume conduction and source leakage, its ability to disregard zero-lag connectivity, its sensitivity to non-linear dynamics, and its low computational complexity. The PLV was calculated as follows:

$$\begin \textrm = \left| \mathbb \left[ S_ /\left| S_\right| \right] \right| \end$$

(3)

where \(\mathbb [\cdot ]\) represents the expected value over windows, and \(S_\) is the cross-spectral density between EEG signals \(x\) and \(y\).

However, the PLV measure is sensitive to volume conduction and source leakage effects. To address these limitations, the measure was extended to ciPLV, which reduces the influence of these confounding factors by isolating the imaginary component of phase differences and applying normalization. It is formulated as follows:

$$\begin \textrm = \frac\left[ \Im \left( }/|}\right) \right] \right| }\left[ \Re \left( }/|}\right) \right] \right| ^2}} \end$$

(4)

where \(\Im (\cdot )\) represents the imaginary part, while \(\Re (\cdot )\) represents the real part. The ciPLV was computed for all possible pairs of the 16 electrodes, producing a 16-by-16 connectivity matrix with 120 unique connections, as ciPLV is an undirected measure. The ciPLV values range from 0 to 1, where 0 indicates no phase synchronization and 1 indicates perfect synchronization. The calculations were conducted separately for four EEG frequency bands: theta (3–7 Hz), alpha (7–13 Hz), beta (13–35 Hz), and gamma (35–50 Hz). The resulting values were averaged across sessions, yielding two connectivity matrices (i.e., pre- and post-training sessions) for each participant and frequency band. Function connectivity estimation was carried out using the MNE-Connectivity Python package (Gramfort et al. 2014).

The resulting connectivity matrices may contain weak or spurious connections, which may introduce noise in subsequent analyses; therefore, a percolation analysis was performed to eliminate these connections. The analysis iteratively removes the weakest connections until the graph becomes disconnected, preserving the last connected state. This is analogous to identifying the maximum threshold that removes weak connections while maintaining the full connectivity of the graph.

Algorithm 1

Network-based statistic (NBS) (Adapted from Zalesky et al. (2010))

Node-wise analysis

Closeness centrality (\(C\)) was used to assess the importance or influence of a node within a network. This metric is based on the assumption that a node with a short average path length can interact with many other nodes through only a few intermediate links, indicating its topological centrality (Beauchamp 1965). In other words, nodes with high closeness centrality can quickly reach other nodes in the network, making them crucial for efficient communication or information flow. Closeness centrality is calculated as:

$$\begin C(i) = \frac l(i, j)} \end$$

(5)

where \(N\) is the number of nodes and \(l(i, j)\) is the shortest path length between nodes \(i\) and \(j\). Note that the path length here is defined as the reciprocal of the weight, which is calculated by ciPLV, because higher connectivity (represented by stronger connections) corresponds to shorter path lengths, reflecting more efficient communication between nodes.

Edge-wise analysis

The network-based statistics (NBS) toolbox (Zalesky et al. 2010) was used to detect subnetworks with edges showing significant differences between sessions while controlling for the family-wise error rate. The null hypothesis assumes that the edge weights or connectivity strengths have equal means across the two sessions. The key steps of the NBS procedure are outlined in Algorithm 1. Changes in network components were assessed using a t-value threshold (\(T\)) of 2.75. Permutation testing (\( k = 100,000 \)) was subsequently applied to calculate the p-value. Statistical significance was determined with a threshold of \( \alpha = 0.05 \).

Correlation analysis and linear regression model

For each frequency band, the goal was to determine how well the node metric predicted behavioral performance. As the metrics may introduce multicollinearity in linear regression, principal component regression (PCR) was applied. Initially, performance data from both hands were reduced to a single component using principal component analysis (PCA). As a feature selection step, Spearman correlations between the bimanual performance and node closeness centrality were computed. However, if participants showed improved performance between sessions, this training effect may be associated with changes in connectivity. Hence, correlations between performance and connectivity could be confounded by these session-related or participant-specific effects. To address this potential confounding, partial correlations were computed while controlling for both session and participant. If the partial Spearman correlations, after Bonferroni correction with a significance threshold of \( \alpha = 0.05 \), were found to be significant, partial ordinary least squares regression was then conducted using the following formula:

$$\begin y' = \beta _0 x' + \beta _1 \end$$

(6)

where \(y'\) represents the residual of the dependent variable \(y\) (bimanual performance) after removing the effect of the control variables. Similarly, \(x'\) is the residual of the predictor variable after accounting for the influence of the control variable. \(\beta _0\) and \(\beta _1\) are the slope and intercept of the model, respectively. The model was trained and evaluated using leave-one-out cross-validation to evaluate the predictability.

Connectivity strength was not included in the analysis, as the percolation process led to a high number of zero connections, which could diminish the statistical power of both correlation and regression analyses.

Statistical analysis

Statistical analyses were conducted using R (version 4.4.3), focusing on group-level effects. Due to the non-normality and skew of the data, a generalized estimating equations (GEE) approach with a gamma log link function was employed to evaluate differences. This was implemented using the geepack package (Højsgaard et al. 2006). GEE accounts for within-subject correlations by incorporating a dependence structure, with an exchangeable correlation structure assumed in this analysis.

To ensure adequate statistical power for detecting a within-subject effect, a sample size estimation was performed using the lmmpower function from the longpower package (Iddi and Donohue 2022; Diggle 2002). Assuming an expected slope of .6, a residual standard deviation of .5, an intra-subject correlation of .6, and two repeated measurements per subject, the analysis indicated that a sample size of N = 18 would be required to achieve.8 power at a significance level of .05.

For the node-wise analysis, closeness centrality was modeled as a function of session, node, and their interaction. Participants were treated as the clustering variable to account for repeated measures at the group level.

To further interpret significant effects, estimated marginal means (EMMs) were computed using the emmeans package (Lenth 2025). Pairwise comparisons were performed for the session \(\times \) node interaction (node-wise), with a Bonferroni correction applied for multiple comparisons. Statistical significance was set at \( \alpha <.05 \).

A similar GEE model was applied to the behavioral analysis, examining the effects of session, hand (left vs. right), and their interaction on behavioral performance. Additionally, edge-wise statistical analysis was performed using NBS as described in Sect. 2.4.