Participants

We recruited 16 healthy participants (sub1–sub16, 25.9 ± 2.7 years old, 6 males, 10 females). All subjects are right-handed. Fourteen of them completed the Edinburgh Handedness Inventory (https://www.brainmapping.org/shared/Edinburgh.php) and obtained average scores of 91.9 ± 8.5 on the augmented index and 90.3 ± 10.5 on the laterality index. Before the experiment, all subjects were informed about the study details and gave their consent. Six of them (sub11–sub16, 4 females) underwent a multi-session experiment (Supplementary Table s1), with 3, 3, 3, 2, 2, and 5 sessions, respectively. Therefore, there are a total of 28 sessions of the experiment. Participants were remunerated per session. This study was approved by the Ethical Committee of the University Hospital of KU Leuven (UZ Leuven) under reference number S6254.

Experiment setup

During the experiment, subjects needed to follow the instructions shown on the screen (ViewPixx, Canada) in front of them, while their brain signals and finger trajectories from their right hand were simultaneously recorded. We used high-density EEG, 58 active electrodes covering frontal, central, and parietal areas with positions following the 5% electrode system [34], and a Neuroscan SynAmps RT device (Compumedics, Australia) for recording. The ground electrode was set at AFz, and the reference electrode was at FCz. All electrode impedances were kept below 5 kΩ before recording. The sampling rate was set to 1000 Hz. Right-hand finger flexions and extensions were tracked using a digital data glove (5 Ultra MRI, 5DT, Irvine CA, USA). We designed the experimental paradigm by relying on Psychtoolbox-3 (www.psychtoolbox.net) to synchronize the EEG and glove data per trial.

Finger flex-maintain-extend paradigm

We designed a finger flex-maintain-extend paradigm including both individual (5 fingers) and coordinated (4 gestures) finger movements, as shown in Fig. 1b. The ‘no movement’ class was designed as the baseline for comparison. A single-session experiment comprises 30 blocks, with each block consisting of a single round of the 10 finger movement scenarios, namely, 30 trials per scenario for each session. Before the experiment, the subjects were told to relax and keep their right hand naturally open with the palm facing upwards on the table (considered the rest position). Figure 1a shows the timing of an exemplary Thumb trial. At the beginning of the trial, a picture is displayed on the screen for 2 s indicating which movement the subjects need to perform in that trial. When this movement scenario cue disappears, a grey circle shows up and starts shrinking at a fixed speed. On top of the circle, there is a cross with scales. When the circle reaches the outer scale (3 s), the subjects must immediately flex their finger(s) and maintain the action for 4 s. Until the circle reaches the inner scale, the subjects need to immediately extend their finger(s) back to the rest position. They were given 2-s rest between trials. We opt for this shrinking-circle design as it diminishes the effect of visual cues [27]. When a trial was started, the subjects were required to only move the indicated finger(s) according to the scenario cue. For the None class, the subjects had to keep their hand at rest while the circle shrinks.

Fig. 1

Paradigm details. a Timing of a trial. The cross is with scales that indicate when the subjects need to flex or extend the corresponding finger(s). b Different finger movement scenarios. During the experiment, the subject’s hand was positioned on the table with the palm facing upwards

Finger trajectory processing

The kinematic data from the data glove were used to precisely detect the onset of movement (finger flexion and extension). We first obtained the finger trajectories based on normalized bending sensor output. Then, we smoothed the trajectories and calculated the trajectory velocity for each movement’s representative finger. For individual finger movements, the representative one was the cued finger, and for coordinated finger movements, we selected the index finger for Pinch, ThumbsUP, and Fist, and the middle finger for Point. Next, the onset of movement was determined by the time when velocity exceeded a threshold of 0.2 times the maximal value (minimal value for finger extension). A graphical explanation is shown in Supplementary Fig. s1.

EEG data preprocessing

Data preprocessing was done by customized scripts and Fieldtrip functions [35]. Raw EEG data were first downsampled to 250 Hz for ease of computation. An antialiasing filtering was applied during this process. Then, the power line noise at 50 Hz was removed by a 3rd-order two-pass band-stop Butterworth IIR filter. Using the same type of band-pass filter, the EEG data were filtered between 0.1 and 70 Hz. We visually inspected faulty channels and excluded them for further preprocessing. Next, Independent Component Analysis was used, and components related to eye movements and abnormal artifacts were identified and removed. Last, the cleaned data went through common average referencing (CAR). We epoched the recordings according to trial markers once the continuous EEG data were preprocessed. We used 4 criteria to find bad channels in each trial. Specifically, a channel was considered bad when any of its kurtosis, mean value, and variance exceeded three times the standard deviation of the mean for all electrodes, or its peak-to-peak amplitude exceeded 200 microvolts. Bad trials, either noisy or containing undesired finger movements, were determined by visual inspection. The bad channel and trial information were kept. For later analysis, bad trials were excluded, and faulty channels and bad channels were interpolated with the average value of neighboring ones. We used the triangulation method in Fieldtrip to calculate each electrode’s neighbors. Ultimately, we obtained an average of 28.4 ± 2.0 clean trials per movement scenario across subjects and sessions.

EEG correlates

Each subject’s single-session data were analyzed to investigate EEG correlates, in which we focused on low-frequency band signals and ERD/ERS.

Low-frequency band signals

We obtained cleaned epochs in the low-frequency band (0.3–3 Hz). For each epoch, we looked into 2-s pre-movement and 2-s post-movement by indexing finger movement (flexion and extension) onset according to kinematic information. For the None case, we extracted epochs according to the corresponding trial marker. We averaged all epochs per finger movement which resulted in a low-frequency EEG template of dimensions 58 × 1000 × 10 (channels × time points × finger movements) for each subject. Two aspects of low-frequency EEG correlates were analyzed, i.e., their time series and MRCPs. The first aspect was examined by showing the temporal evolution of the amplitude topoplots between each movement and None. The second aspect was to analyze, for selected representative channels, their MRCPs.

ERD/ERS

We first segmented the cleaned data based on trial markers and then extracted the epochs. Then, we implemented the Morlet wavelet time–frequency transformation (ft_freqanalysis() function in Fieldtrip) on each epoch. The frequency of interest was set to 1–50 Hz with a resolution of 1 Hz. The time resolution was set to 0.01 s. The resulting power spectra of all trials from the same movement were averaged per subject. Finally, ERD/ERS for each movement was derived as:

$$ERD/ERS\left(f,t,c\right)=\frac_}\sum_^_}\frac_\left(f,t,c\right)-_^\left(f, c\right)}_^\left(f, c\right)}\times 100\%$$

(1)

where \(_\left(f,t,c\right)\) denotes the power spectrum of the i-th subject at frequency f, time t, and EEG channel c. \(_^(f,c)\) is the baseline, selected from the middle 1-s power spectra of the None movement, and Ns the total number of subjects. A negative value indicates ERD and vice-versa ERS.

Similarity analysis and clustering

We relied on Riemannian distance as the dissimilarity metric to assess the finger representations obtained from EEG [36, 37]. First, we obtained cleaned epochs during finger flexion and extension within the 0.3–70 Hz frequency band. One epoch contains 1-s pre-movement and 0.5-s post-movement. Then, each epoch was transformed into a covariance matrix that lies in the Riemannian manifold [38]. For finger flexion or extension, we obtained the centroids of each type of finger movement’s covariance matrices and calculated the pairwise Riemannian distance between them. A larger distance indicates a larger dissimilarity. Finally, we could get a symmetric representational dissimilarity matrix (RDM) that reflects the structure of the broadband EEG responses for different finger movements during flexion or extension. Hierarchical clustering was done based on this matrix. The resulting dendrogram was analyzed by looking into clustered finger movements.

Decoding models and implementation detailsFeature extraction

We tested mainstream feature extraction methods in the literature related to hand and upper-limb movement classification tasks. First, we obtained the cleaned epochs from multiple frequency bands, including the low-frequency (0.3–3 Hz), delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and low gamma (30–70 Hz) bands. Then, for the low-frequency band, feature extractors including time-domain amplitude [27, 29, 30, 39], discriminative spatial patterns (DSP) [40], and discriminative canonical pattern matching (DCPM) [41] were implemented. For the other frequency bands, we extracted band power, common spatial pattern (CSP) [42], and Riemannian geometry tangent space (RGT) [38] features.

Implementation

Denote the i-th EEG trial as \(}\left(i\right)\in }^\), where C and P indicate the number of channels and sampling points, respectively. In this study, we fixed the time window to be 1.5 s, and thus P = 375. The time-domain amplitude was taken every 0.12 s, resulting in a 754-dimensional feature vector for each trial with C = 58 channels. Considering the scarcity of training data, we removed redundant features using Lasso regularization with a regularization coefficient of 0.05 [43]. For DSP, we selected the top 10 eigenvectors as spatial filters, hence the trial channel dimension was reduced from 58 to 10. Then, the average value of each channel was extracted as a feature. For DCPM, the top 10 eigenvectors were selected as spatial filters during computation, and finally, the model outputs a 3-dimensional feature vector for each trial. To extract power features, we took the average square value of each trial’s channel, resulting in a 58-dimensional feature vector. For CSP, we selected the paired first and last 3 spatial filters and generated a 6-dimensional feature vector. Last, for RGT, the mapped features in Riemannian Tangent space have an original dimension of C × (C + 1)/2 but were reduced again using Lasso regularization. Note that the above description of feature dimensions is the theoretical value of a trial in one single frequency band. When multiple frequency bands’ information was fused, the dimensions changed accordingly.

Classification task

The shrinkage linear discriminant analysis (sLDA) model was used as the classifier for its excellent performance in single-trial EEG classification [27, 30, 39, 44]. The regularization parameter was set to 0.8 according to a trial–error test on one subject. Based on this model, we aim to investigate: (I) which feature extractor and frequency band contributes the most to finger movement detection and pairwise classification, and (II) whether we could build a model based on those contributing features and obtain an overall performance improvement. In order to address task I, we trained and tested several sLDA classifiers based on each feature type on the low-frequency, delta, theta, alpha, beta, and low gamma band, individually. The contribution of features and frequency bands was analyzed. For task II, we gathered all classifiers trained on the selected features and frequency bands and used majority voting for prediction. Based on this ensemble model, we also investigated the impact of data sizes, time window choices, and EEG electrode layouts on model performance. We tested the impact of data sizes on decoding performance with increasing sub11-sub16’s multi-session data. For time window choices, we considered three primary time windows [−1.6, −0.1]s, [−1, 0.5]s, and [0, 1.5]s to the movement onset (0 s). Likewise, different EEG electrode layouts (Supplementary Fig. s2) were selected from the 58 electrodes in total and compared based on the ensemble model. All data (28 sessions) were used for each task, except for the data size testing one. We performed tenfold cross-validation on all tasks. All trained models are subject-specific. The chance levels were estimated following [45], and we obtained 0.6225 for single-session and 0.5722 (estimated based on the averaged trial numbers across sub11–sub16) for multi-session pairwise classification (alpha = 0.05).

Statistical analyses

All statistical tests were conducted using MATLAB with a significance level of 0.05. For multiple within-factor conditions, such as finger movement scenarios, time windows, and electrode layouts, we relied on one-way repeated measures ANOVA. Then, a post hoc multiple comparisons test with Bonferroni correction was used to identify significant pairwise differences. For pairwise conditions, we relied on the Wilcoxon signed-rank test. We marked associated p-values using asterisks in figures (N.S.: not significant; *: p < 0.05, **: p < 0.01; *** p < 0.001).