Decoding Coordinated Directions of Bimanual Movements From EEG Signals

Bimanual coordination is common in human daily life, whereas current research focused mainly on decoding unimanual movement from electroencephalogram (EEG) signals. Here we developed a brain-computer interface (BCI) paradigm of task-oriented bimanual movements to decode coordinated directions from movement-related cortical potentials (MRCPs) of EEG. Eight healthy subjects participated in the target-reaching task, including (1) performing leftward, midward, and rightward bimanual movements, and (2) performing leftward and rightward unimanual movements. A combined deep learning model of convolution neural network and bidirectional long short-term memory network was proposed to classify movement directions from EEG. Results showed that the average peak classification accuracy for three coordinated directions of bimanual movements reached <inline-formula> <tex-math notation="LaTeX">$73.39~\pm ~6.35$ </tex-math></inline-formula>%. The binary classification accuracies achieved <inline-formula> <tex-math notation="LaTeX">$80.24~\pm ~6.25$ </tex-math></inline-formula>, <inline-formula> <tex-math notation="LaTeX">$82.62~\pm ~7.82$ </tex-math></inline-formula>, and <inline-formula> <tex-math notation="LaTeX">$86.28~\pm ~5.50$ </tex-math></inline-formula>% for leftward versus midward, rightward versus midward and leftward versus rightward, respectively. We also compared the binary classification (leftward versus rightward) of bimanual, left-hand, and right-hand movements, and accuracies achieved <inline-formula> <tex-math notation="LaTeX">$86.28~\pm ~5.50$ </tex-math></inline-formula>%, <inline-formula> <tex-math notation="LaTeX">$75.67~\pm ~7.18$ </tex-math></inline-formula>%, and <inline-formula> <tex-math notation="LaTeX">$77.79~\pm ~5.65$ </tex-math></inline-formula>%, respectively. The results indicated the feasibility of decoding human coordinated directions of task-oriented bimanual movements from EEG.


I. INTRODUCTION
B RAIN-COMPUTER interface (BCI) can provide a communication pathway between the brain and external devices [1], [2]. To date, a large number of studies have demonstrated that kinematic parameters of hand movements are strong relative to non-invasive electroencephalogram (EEG) signals [3], [4], [5]. Parameters such as position [6], [7], speed [7], [8], and direction [9], [10], [11] have been successfully decoded from these recordings. It is believed that such investigations focusing on hand kinematics decoding can potentially lead to intelligent manipulation of upper limb robots [12] and hand-related rehabilitation [13].
Bimanual coordination is very important in daily life and necessary in arm rehabilitation to achieve complete functional recovery. Two hands are mostly required to be used simultaneously to manipulate an object or perform a task [14]. Bilateral coordination training plays an irreplaceable role compared with unilateral training in robot-assisted rehabilitation and thus has been attracting lots of attention [15], [16]. Furthermore, direction decoding is one of the most popular BCI topics to achieve independent multidimensional control, which can be used in many rehabilitation applications [17]. These motivated us to verify the feasibility of decoding coordinated directions of task-oriented bimanual movements [18] from EEG, which has been reported in very few studies.
Currently, most existing studies in decoding hand directions from EEG signals are limited to a single hand. Some studies focused on classifying any two directions of four orthogonal directions (a total of six binary classifications) [19], [20], [21]. For example, Robinson et al. decoded orthogonal directions of a single hand using a waveletcommon spatial pattern algorithm with an average accuracy of 87.85% across six binary classifications [19]. In a similar paradigm, Chouhan et al. achieved an average accuracy of 76.85% via wavelet phase-locking values [21]. There is also a study in which a four-class classification of decoding centerout single-hand movement directions was performed with the phase-locking value method, achieving an average decoding performance of 67% [9]. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Compared with the majority of studies focusing on singlehand kinematics, to the best of our knowledge, there is only one study that attempts to classify bimanual directions from EEG signals in 2021 [11]. The authors conducted a multiclass classification, including two unimanual and four bimanual movements, by means of linear discriminant analysis (LDA) and support vector machine (SVM) classifiers based on temporal and spectral features. The experimental paradigm was designed as follows: (1) regarding the right hand for a dominant role and left hand for secondary movement, (2) moving the right hand in horizontal directions (right or left) and left hand in vertical directions (up or down). The peak accuracy reached 70.29% ± 10.85%. One point to be noted is that this study used orthogonal manipulation of both hands in asymmetric space (right hand along with horizontal directions and left hand for vertical directions). Such bimanual mode is uncommon in daily life. In contrast, some practical life tasks require the parallelly coordinated use of both hands in symmetric space [22], [23], such as moving a heavy box and pulling a chest expander.
Decoding natural bimanual movements has also been investigated, but mostly on non-human primates utilizing invasive brain signals. Two typical examples, Ifft et al. [24] developed a BCI that enables monkeys to control two avatar arms simultaneously in two-dimensional reaching tasks via decoding large-scale neuronal signals recorded by implanted multielectrode arrays. Choi et al. [25] decoded monkeys' bimanual movements of three-dimensional reaching tasks, but with a different type of signal (epidural electrocorticography). Although these techniques have decoded bimanual arm movements from brain signals successfully, problems related to security and the long-term usability of invasive BCI can lead to limited clinical applications [26], [27].
To summarize, while bimanual movements are indispensable to activities of daily life, effective decoding of them from non-invasive EEG signals is still in a very primitive stage for potential robotic control or rehabilitation purposes. Hence, this study proposes a BCI to decode coordinated directions of task-oriented bimanual movement from EEG recordings. Compared with decoding bimanual movements in asymmetric spatial orientations [11], the proposed BCI paradigm focuses on task-oriented coordinated movements in spatial symmetry, which is more natural to match real-life tasks. The main contributions of this study include: (1) a new BCI paradigm of task-oriented bimanual movements for decoding coordinated directions was proposed for the first time; (2) its feasibility was experimentally demonstrated with eight human subjects.
The paper is organized as follows: Section II presents the experimental paradigm, human participants, and data processing, together with a hybrid classification model by combining convolution neural network (CNN) and bidirectional long short-term memory (BiLSTM) neural network [28]; Section III presents the experimental results in terms of brain activity analysis and classification performance; Section IV presents the discussion, followed by the conclusion in the end.

A. Subjects and Experimental Paradigm
Eight healthy subjects (5 males, and 3 females), aged between 21 and 28 years, participated in the experiment. This study was approved by the Ethics Committee of the Southern University of Science and Technology (20190004), Shenzhen, China. All participants signed informed consent and were confirmed to be right-handed by the Edinburgh Handedness Inventory.
The experimental setup is illustrated in Fig. 1 (a). The experiments were carried out in an electromagnetically shielded room. The subjects were required to sit on a chair in front of a desk. A 43-inch display monitor, showing the rectangle target and the dots representing hand position, was placed at a comfortable viewing distance of approximately 2 m from the subject. A guide rail was used to direct subjects to move hands in one-dimensional space. The tasks performed in the experiment involved bimanual and unimanual targetreaching movements. For bimanual movements, the target randomly appeared on the left, middle, or right of both hands, as indicated in Fig. 1 (b), respectively, corresponding to leftward, midward, and rightward bimanual movements. Differently, the target randomly appeared on either left or right for unimanual movements as shown in Fig. 1 (c) and (d). Both hands needed to reach the appeared target during bimanual tasks simultaneously. In the unimanual version of tasks, subjects only moved a single hand to the appearing target but remained the other hand still.
Subjects were pre-trained on the tasks for about 15 mins before the formal experiment and were asked to avoid blinking and body movements during experiments. The formal experiments were divided into three sessions: both-hand, lefthand, and right-hand movements. Each session contained six runs, with a rest period of 2 mins between two runs. Each run of bimanual session consisted of 15 trials that are five leftward, five rightward, and five midward movements. Each run of unimanual sessions included 10 trials that are five leftward and five rightward movements. Thus, there were a total number of 210 trials for each subject.
The timeline of a single trial was consistent across all sessions. An example of bimanual movements is presented in Fig. 1 (e). To be specific, each trial started with the appearance of one red rectangle target on the monitor. After a 2 s preparation, two dots appeared at home positions indicated by virtual hands in Fig. 1 (b). The home position was defined as the starting position of movement and was fixed across all trials. Both hands needed to move back to the home positions before the next trial. During movement, the red dot position was controlled directly by the real-time position of human hands along the guide rail. The task required human participants to reach the target within 5 s and then hold on for 1 s. Only when the center of the red dot entered the red rectangle within 5 s, a trial was considered successful. Otherwise, trials were aborted and considered unsuccessful if the dot had not reached the target after 5 s. A successful trial was indicated by the target changing from red to yellow. After that, all visual cues disappeared on the monitor and subjects had a rest of 5 s before the next trial.

B. Data Acquisition
EEG signals of 64 channels were collected from eight subjects by the g.HIAMP 256 biosignal amplifier (g.tec medical engineering GmbH, Austria) at a sampling rate of 256 Hz. Electrode impedances were calibrated to be less than 30 K . The electrode locations were referred to as the international 10-10 system with a forehead ground electrode at AFz. A notch filter of 50 Hz was used to remove the powerline interference. The movement position and velocity of both hands were captured by an optical hand tracking module called Leap Motion Controller (Ultraleap Inc., California, USA).

C. Data Preprocessing
EEG signals were first band-pass filtered from 0.1 to 30 Hz by using a zero-phase filter and then down-sampled to 100 Hz. Generally, the bad channels appear noisier than the other surrounding channels and their spectra are outliers. Therefore, the bad channels were first identified visually by browsing the raw data and confirmed by channels' spectra. Then bad channel repair was implemented using spherical spline interpolation [29]. In this method, all channel locations were projected onto a unit sphere and the signal at the bad locations was interpolated based on the signals at the good locations. Next, the EEG data were re-referenced to the averaged earlobes (i.e., A1 and A2) and the common average of all channels successively. The motion-related artifacts were removed by the artifact subspace reconstruction (ASR) algorithm [30]. Additionally, eye blinks and eye movement artifacts were visually identified and removed from EEG using the independent component analysis (ICA) method [31]. The EEG signal generally consists of a mixture of various brain and noise sources. ICA is a widely-used blind source separation technique and is particularly useful for the removal of artifacts embedded in EEG. First, the raw multichannel EEG was decomposed into a sum of temporally independent components (ICs) by the ICA algorithm. Then, ICs related to eye movements and blinks were identified and removed by their tempo-spectral features and scalp distribution. Finally, the artifact-corrected EEG was reconstructed using the inverse ICA transformation on the remaining ICs.
In general, there is a delay from movement cue onset to actual movement onset because of reaction time, and the delay is subject-specific across trials and individuals according to personal response-ability. Hence, we calibrated the time of actual movement onset by Leap Motion Controller for each participant, defined as subjects-specific actual movement onset. The 2 s epoch between −0.5 and 1.5 s from the actual hand movement onset during each trial was extracted from EEG for analysis. Each epoch was then baseline corrected by subtracting the average voltage of data between −0.5 and 0 s of movement onset. Finally, we performed a bandpass filter between 0.1 and 4 Hz to epochs before feeding them into the decoding model. Unsuccessful trials of approximately 5.6% were deleted, which may cause a side effect on classification accuracy. We implemented all data preprocessing in opensource MNE-Python software [32].

D. Decoding Model
The changes of spatial and temporal dynamics in different EEG frequency bands are related to directions of hand movements. Therefore, effective extraction of corresponding features is the key to high decoding performance. In this work, we proposed a hybrid deep learning model based on CNN and BiLSTM to extract complex features automatically from EEG and decode movement directions. The visualization and full description of the network architecture are shown in Fig. 2 and Table I, respectively. The structure and parameters are the same for multiple or binary classification used in this study except for the last full connection layer and activation function.
1) Network Architecture: The proposed architecture was built on temporal CNN block, spatial CNN block, and temporal BiLSTM block, enabling spatial-temporal frequency features extraction. All related features were integrated within the final fully connected layer for calculating output classification probability. The CNN-based layers were configured without a bias and followed by a batch normalization (BatchNorm), exponential linear unit (ELU) activation function, dropout, and average pooling layers except for the first full convolution. The temporal and spatial CNN block was inspired by EEGNet [33].  EEGNet has proven to be an effective model for EEG feature extraction, which has compact architecture and generalization ability across different paradigms. These abilities are very important, especially for our small EEG dataset. Hence, the same structure of EEGNet was used in this study to learn EEG features. Whereas, the parameters configuration was adjusted and optimized according to our paradigm. For example, the kernel size of the first convolution layer was set to be (1, 50) due to the lower sampling rate of EEG. The dropout rate was 0.3, and the kernel size of the pooling layer was (1, 2). More details are described as follows.
In the temporal CNN block, we fitted eight convolution filters of size (1,50). The output of eight-channel feature map contained the EEG signal at different frequency bands. No activation function was used after the first convolutional layer because it had no significant effect on the final accuracy.
In the spatial CNN block, we performed two CNN layers in sequence, called the DepthwiseConv2D layer and SeparableConv2D layer, respectively. Sixteen-channel feature map were outputted by applying a depthwise spatial convolution with a kernel size of (N c , 1) and depth multiplier D = 2, where N c indicates the number of EEG channels.
The depthwise convolution acted on each channel of feature map separately. To combine the information across the feature map channels, the SeparableCon2D layer was then used, which consisted of first performing a depthwise convolution followed by a pointwise convolution. The pointwise convolution performed 1 × 1 convolution to capture new features by mixing the resulting output channels of feature map.
It is well known that BiLSTM can effectively increase the amount of information available to the network by learning the temporal dependence in both forward and backward directions of time series [28], [34]. Considering the information in EEG signals before and after actual movement, which may contribute to coordinated directions classification, BiLSTM was used to further extract temporal features of EEG following CNN blocks.
In BiLSTM block, input data was first fed into the forward LSTM layer, and then reversed and fed into the other backward LSTM layer as shown in block C of Fig. 2. Each layer of LSTM contains 25 LSTM cells, each of which consists of 16 hidden units. A single LSTM cell was implemented and updated as follows: where W f , W i , W o denote the weighted matrices and b i , b f , b o denote the biases in LSTM cell. σ and are the sigmoid function and element-wise multiplication, respectively. h t indicates a hidden state. The output of BiLSTM at the step t is represented as follows: where f h t and bh t correspond to forward and backward LSTM, respectively.
2) Training Strategy: Before training the network, input data were preprocessed as follows: (1) data normalization, (2) data reshape. By subtracting the mean and then dividing by the standard deviation, data normalization first transformed data to be zero mean and unit variance for improving model accuracy. Then we reshaped three-dimensional sample data with format (N s , L, N c ) to four-dimension with format (N s , L, N c , 1), where N s is the number of input samples, L is the data length of EEG, and N c is the number of EEG channel. Models were trained using an adaptive moment estimation (Adam) optimizer. For multi-classification, binary cross-entropy and softmax were selected as loss function and activation function in the final layer, respectively. In contrast, cross-entropy loss function and sigmoid activation function were used for binary classification.
To obtain the optimal hyperparameters such as epoch,  determined by grid search techniques for the highest crossvalidation score. First, 10-fold cross-validation was used to split the entire dataset into 10 non-overlapping groups. Then, nine folds were used as a training set and the remaining one fold was used as a hold-out set for validation. This process was repeated 10 times and each time a different fold was used for testing. Therefore, a total of 10 models were trained to cover all combinations of parameters. The combination that yielded the best average performance over 10 models was identified as the optimal parameter. Furthermore, to avoid data leakage, data preparation such as normalization was prepared only on the training set instead of on the entire dataset, and then applied to the training and test sets within each fold of the crossvalidation. That is, we first split data into training and test sets, then fitted the data preparation on the training set, and applied the transform to the training and test sets. Finally, the hyperparameters were set to epoch = 150, batch size = 16, and number of neurons = 16.
Furthermore, there are a total of 6307 and 6241 trainable parameters for multiple and binary classification used in this study as shown in the 5 th column of Table I. To avoid overfitting during training these complex models, several measures were taken: (1) the 10-fold cross-validation was used to assess the overall performance of the model, (2) dropout was added after each convolutional layer to reduce the complexity of the model, (3) we would also early stop the learning process when noticed the accuracy on the validation set remained unchanged or even decreased over a series of epochs. (4) learning decay with a factor of 0.2 were implemented to adjust the learning rate only when the optimizer cannot improve the results over a patience number of epochs.

E. Statistical Analysis
To estimate the differences between the groups and conditions, a one-way analysis of variance (ANOVA) was used. Classification results were tested for normality by using the Shapiro-Wilk test. Leven's test was then applied for the homogeneity test. Results with significance were calculated for all pairwise comparisons with Tukey's honest significant difference method. A significance threshold was set to p = 0.05 in all tests. The classification accuracies were presented as mean ± standard deviation (std) in all tables.

III. RESULTS
In this section, we first analyzed brain activity changes related to bimanual and unimanual movements in time and frequency domains, including event-related desynchronization/synchronization (ERD/ERS), topographical maps, and movement-related cortical potentials (MRCPs). Then, three coordinated directions (i.e., leftward, midward, and rightward) were classified to validate the feasibility of decoding bimanual movements. Finally, to compare the decoding performance of bimanual and unimanual movements, we calculated binary classification accuracy in the leftward and rightward directions of both bimanual and unimanual movements.
A. Movement-Related Changes in EEG Activities 1) Time-Frequency Representation: Task-related imagery and execution of movements have been found to lead to amplitude (or power) decrease in the alpha band of EEG, known as ERD. We computed time-frequency representation to observe this physiological pattern using the Morlet wavelet transform. Fig. 3 showed the averaged time-frequency representation at C3 and C4 channels across subjects. The ERDs associated with movement execution were observed on the alpha band under all movement conditions. Strong ERDs were presented in the contralateral hemisphere during the unimanual tasks, which reflected sensorimotor activations at C4 and C3 during left-hand and right-hand movement, respectively. These contralateral ERDs were also accompanied by weaker ipsilateral ERDs during unimanual tasks. In contrast, symmetrical bilateral ERDs were evident during bimanual movement, which indicated the cooperative activations of two cerebral hemispheres.
2) Topological Maps: Fig. 4 showed the topographical maps in the time window from −0.5 to 0.9 s of actual movement onset for all EEG channels. The obvious changes of evoked potential occurred around 0 s. Under unimanual movement conditions, noticeable evoked potential changes were observed over the contralateral brain hemisphere. However, bilateral activations were apparent during bimanual movements. These patterns were broadly consistent with the findings in time-frequency representation analysis mentioned above. Furthermore, most of these patterns were located on C3, C4, and the surrounding areas.
3) Movement-Related Cortex Potentials: To further study the differences under different movement conditions, we also investigated the averaged MRCPs on C1, Cz, and C2 channels as shown in Fig. 5. The first row in Fig. 5 presented the MRCPs of bimanual movement along directions of leftward, midward, and rightward. A positive peak appeared at about −0.25 s before the actual movement onset. Thereafter, the amplitude pronouncedly decreased to the negative peak which occurred between 0.25 and 0.75 s after movement onset. To compare the difference of MRCP patterns between bimanual and unimanual movements, the curves of MRCP under leftward and rightward movement conditions were presented in the second row of Fig. 5. It is noted that MRCP during bimanual midward movement was not provided here due to no so-called unimanual midward movements. Similar results were found between bimanual and unimanual movements in terms of the occurrence of positive and negative peaks. In addition, to further compare their difference in amplitude, a statistical analysis was performed on the negative and positive peaks of MRCP, as shown in Figs. 6 and 7, respectively. Bimanual movements generated a lower negative peak of MRCP compared with unimanual movements in spite of no significance in some cases. However, no similar results were found for positive peaks.

B. Decoding Bimanual Movements Performance
To verify the feasibility of decoding coordinated directions of bimanual movements, a multi-classification by the proposed model was performed. Table II presented the accuracies of classifying three coordinated directions (i.e., leftward, midward, and rightward) from EEG. The mean and std in the last row indicated the grand average and standard deviation of mean values on all subjects. The grand average peak accuracy of decoding bimanual movement directions reached 70.10% ± 6.88% in the 2 nd column of Table II. This grand average peak accuracy was calculated using a general time window (0 s, 1 s) corresponding to the start time point of 0 s (green dot) as shown in Fig. 8. The subject-specific model performances in Fig. 8 were calculated using a shifted window with a fixed length of 1 s and a step size of 0.1 s. The start time points of shifted windows ranged from −1.4 to 1.4 s of actual movement onset. For each shifted window, we trained and evaluated the model based on 10-fold cross-validation for each subject. Each value of accuracies for each subject was the average from 10fold cross-validation. Consequently, each subject had a specific accuracy peak (red dots in Fig. 8) at a different starting time point. The mean value of subject-specific peak accuracies reached 73.39 ± 6.35% as shown in the 3 rd column of Table II. It was noted that all subject-specific peak accuracies occurred after the visual cue onset, more specifically, between visual cue and actual movement onset. This time interval was defined as visual-motor reaction time, which was presented by the green arrow in Fig. 9. The visual-motor reaction time was subject-related, and even not equal between trials within one experimental run. Obtaining such visual-motor reaction time  in each trial, all trials were aligned to the actual movement onset which is marked as t = 0 s in Fig. 8.
Furthermore, we performed binary classification for directions of leftward versus midward, rightward versus midward and leftward versus rightward during bimanual movements. The results were presented in Table III. The mean and std in the last row indicate the grand average and standard deviation of mean values on all subjects. The time window was set to (0 s, 1 s). The highest accuracy was obtained for the leftward versus rightward classification with a mean accuracy of 86.28 ± 5.50%. The average classification accuracies

C. Decoding Performance Between Bimanual and Unimanual Movements
We also compared the decoding performance between unimanual and bimanual movements for classifying directions of leftward and rightward. The results were shown in Table IV. Two traditional feature-based methods (i.e., LDA and SVM) used in the previous work [11] and three deep learning-based methods including BiLSTM, EEGNet, and the combination of ShallowConvNet [35] and BiLSTM were also used for comparison. All models were trained in the same way. The mean and std indicate the grand average and standard deviation of mean values on all subjects. Results showed that our proposed method yielded the highest average accuracies across eight subjects under all movement conditions. For example, compared with SVM, 14.61%, 17.38%, and 12.56% average increase was found in decoding accuracies under left, right, and bimanual movements, respectively. The average increase of 2.67%, 1.89% and 4.28% was achieved compared to ShallowConvNet+BiLSTM model. In addition, all models obtained higher classification accuracies on bimanual movement task than unimanual task as shown in the last row of  To assess the statistical differences of classifications under three movement conditions using different methods, an ANOVA was performed for accuracies of eight subjects (Table IV). The statistical results were presented in Fig. 10. It can be seen from Fig. 10 (a) that the proposed method significantly surpassed LDA and SVM methods under all movement conditions. Our proposed model also achieved higher accuracies compared with the other three deep learningbased models under unimanual movement conditions, but there was no significant difference. Especially, the proposed method significantly performed better than BiLSTM and EEGNet under bimanual movement conditions. In Fig. 10 (b), the decoding performance was compared between unimanual and bimanual movement conditions. The proposed method, ShallowConvNet+BiLSTM, LDA, and SVM performed significantly better on decoding bimanual than unimanual movements. It also should be mentioned that no significant difference was found between left-and right-hand movement decoding based on all methods.

IV. DISCUSSION
Bimanual coordination is an important and complex daily skill [14]. In this work, we presented the decoding of coordinated bimanual movements in leftward, midward, and rightward directions for the first time. The experimental results showed that the proposed deep learning model based on the combination of CNN and BiLSTM achieved the accuracy of 70.10 ± 6.88% for classifying three-directional coordinated movements from EEG signals. We also compared the decoding performances of bimanual and unimanual movements in leftward and rightward directions. The binary classification accuracies reached 86.28 ± 5.50%, 75.67 ± 7.18%, and 77.79 ± 5.65% for bimanual, left-and right-hand, respectively, based on our proposed method. In addition, the performance of the proposed decoding model surpassed not only traditional LDA and SVM which are highly dependent on handcrafted features, but also deep learning-based models including BiLSTM, EEGNet, and ShallowConvNet+BiLSTM.
Five different methods were used to decode coordinated bimanual movements. BiLSTM and EEGNet have proven to be great potential in EEG-based BCIs [33], [34]. Effective feature extraction by learning enables them to perform better than traditional methods depending on hand-crafted features. EEGNet constructed by convolutional layers mainly learns temporal and spatial features of EEG but is unable to capture long temporal dynamics. In contrast, LSTM can identify time dependencies on a time series by introducing memory cells [28], which is important in temporal information classification. Consequently, a combined architecture including both EEGNet and BiLSTM was used for the classification of coordinated directions. The results in Fig. 10 (a) indicated the benefit of the combined model with higher accuracy than a single network architecture. Furthermore, a similar approach which consists of ShallowConvNet and BiLSTM was utilized Fig. 10. Statistical comparisons of decoding performance for eight subjects using ANOVA with significant analysis: (a) across classification models, (b) across movement conditions. The significance threshold was * p < 0.05. "ns" within groups denoted non-significant.
for further comparison. ShallowConvNet originally reported by [35] is another popular used model for EEG decoding. It was found that ShallowConvNet tended to perform worse than EEGNet in the case of MRCP classifications [33]. This is the possible explanation for the lower accuracy of ShallowConvNet+BiLSTM than EEGNet+BiLSTM.
Differences in EEG activities between bimanual and unimanual movements were found in this study. Unimanual movement execution led to a significant contralateral ERD in the alpha band. Interestingly, this ERD was accompanied by a weaker ipsilateral ERD during unimanual tasks, which was consistent with findings reported by [36] in unilateral wrist extension tasks. Similar to previous work [37], a clear bilateral ERD was observed during coordinated bimanual movements. Our results also showed that the MRCPs began with an upward deflection at around 250 ms before the actual movement onset. Afterward, a slow negative shift was observed, followed by a negative peak after movement execution, which was identified in the work [38]. We also found that bimanual movements generated a lower negative peak of MRCP compared with unimanual movements as shown in Fig.6. This is consistent with the studies of local field potential based on invasive techniques, in which the authors concluded that the amplitude of movement-evoked potentials in both supplementary motor area (SMA) and primary motor cortex (M1) is larger during bimanual than during unimanual movements [39]. In fact, we also found many previous studies have explored the cortical mechanisms of bimanual coordination based on functional imaging or EEG approaches. Three definite special brain activations of the coordination can be concluded: (1) interhemispheric interactions contribute to bimanual synchronization, and they are significantly higher during bimanual tasks compared with unimanual conditions [40], (2) brain activation levels during bimanual coordination exceed the sum of the activation during single-limb task [14], [24], [43], (3) the primary motor cortex, premotor cortex, and supplementary motor area are more involved during coordination tasks [14], [44]. These special characteristics of the coordination also support that bilateral activation of the brain during bimanual movement may provide more discriminable information in EEG for direction classifications than unimanual movement. This was implied by the higher classification accuracy of bimanual movement ( Fig. 10 (b)) compared with unimanual movement.
Bimanual coordination requires two sides of the body to actively and cooperatively act in executing tasks. There are two basic models in bimanual coordination (i.e., spatial symmetry and asymmetry). The former requires both limbs to move in parallel. This can be further subdivided into isodirectional and opposite movement patterns. In contrast, the latter refers to that the limbs move orthogonally. Compared with the orthogonal mode in [11], the parallel movement pattern in our paradigm is less influenced by spatial interference and exhibits tighter inter-limb synchrony [22], [23]. In general, the parallel mode can be more natural to match task-oriented movements.
One limitation of this study is that the MRCP associated with only actual movements was processed for direction decoding. Motor execution (ME) and motor imagery (MI) are two popular kind of tasks in EEG-based BCI [5]. Many studies have confirmed their consistency in brain activity modulation [41], [42], [43], such as the ERD/ERS [3] and MRCP [36]. For example, the MRCP occurs consistently between ME and MI [44], [45]. Therefore, although we did not take into account the MI task, it is reasonable that the proposed method for decoding ME task can be transferred to MI tasks, as has been done in [46]. This suggests the potential of the proposed BCI paradigm be integrated with our previously developed bilateral robot [12]. Under this circumstance, the patients will be able to control a bilateral robot by simply imagining hand movement without ME. These have been confirmed in unimanual MI tasks on stroke patients [17] and bimanual tasks on monkeys [24]. However, it is also important to note that ME has been shown to provide a higher decoding performance compared to MI due to more pronounced MRCPs [46], [47].
It is well known that task-oriented bimanual movements can happen along one, two, or three dimensions in our daily life. In contrast, the proposed paradigm was implemented in only the horizontal dimension, which can be considered as the other limitation of this study. As for the cases with multiple dimensions, they will be investigated in the future. However, to our best knowledge, this is the first study attempting to decode task-oriented coordinated movements. This study focuses on validating the feasibility of the proposed BCI paradigm as a preliminary study. Furthermore, future research also may consider the perspectives that participants can try MI tasks [48]. In this way, stroke patients are likely able to control the bilateral robot by simply imagining the actions of the bimanual movement guided by the target. We will further validate the rehabilitation outcome of the patient due to MI-based bilateral physical training.

V. CONCLUSION
In this paper, we aimed to decode the coordinated directions of natural bimanual movements from EEG. For that, a bimanual task-oriented BCI paradigm of performing target-reaching tasks was designed. The experimental results mainly showed that (1) peak classification accuracy of 73.39 ± 6.35% was obtained in three coordinated directions (i.e., leftward, midward, and rightward) through a hybrid CNN-BiLSTM model; (2) the binary classification accuracies achieved 80.24 ± 6.25, 82.62 ± 7.82, and 86.28 ± 5.50% in coordinated directions of leftward versus midward, rightward versus midward and leftward versus rightward, respectively; (3) the binary classification accuracies (leftward versus rightward) of bimanual, left-hand, and righthand movements achieved 86.28 ± 5.50%, 75.67 ± 7.18%, and 77.79 ± 5.65%, respectively. These findings support the feasibility of decoding coordinated directions of task-oriented bimanual movements from EEG.