For this study, we collected EEG and EMG data from 14 participants performing MI and created an MI-EEMG dataset. The group consisted of six males (average age, 21.5 ± 3.2 years) and eight females (average age, 22.4 ± 4.6 years). All participants were right-handed, had normal vision, no motor impairments, and signed consent forms before the experiment. Participants were informed regarding the experimental procedures and precautions.

The MI-EEMG dataset included three types of MI actions: left hand, right hand, and resting state. Each subject underwent 30 trials of 10-s MI experiments, as shown in Figure 1. EEG electrodes were placed according to the international 10–20 system (Klem et al., 1999), with 32 channels and a sampling frequency of 500 Hz; four EMG electrodes were located on the left biceps (LB), left triceps (LT), right biceps (RB), and right triceps (RT), all with a sampling frequency of 500 Hz.

The WAY-EEG-GAL dataset included data from 12 participants who were right-handed, had normal vision, and no motor impairments. During the EEG and EMG data collection process, the participants were instructed to reach out and grasp an object when prompted, lift it using their thumb and index finger, hold it for a few seconds, place it back on the support surface, release it, and finally return their hand to a specified resting position. EEG signals were captured using 32 channels from the international 10–20 system, and the EMG electrodes were positioned on the hands and arms. The experiment involved a total of 3,936 EEG and EMG data collections.

The ESEMIT dataset included data from 26 participants, including 10 males (average age, 25.4 ± 7.4 years) and 16 females (average age, 23 ± 3.2 years). All the participants were right-handed, had normal vision, and no motor impairments. The ESEMIT dataset includes three types of MI actions: left hand, right hand, and resting state. Each participant underwent a 4-min resting state and 15 min of MI induced by VR for collecting EEG and EMG data. EEG electrodes were placed according to the international 10–20 system, with 32 channels and a sampling frequency of 500 Hz; two EMG electrodes were located on the LB and RB, with a sampling frequency of 500 Hz.

The original EEG signals were recorded in this study as shown in Equation 1: X = x 1 , x 2 , , x C 1 ∈ R C 1 × L 1 (1) Where C 1 represents the number of channels for EEG acquisition ( C 1 = 32 in this study), L 1 denotes the number of samples, with L 1 = T × f E E G , T is the sampling time, and f E E G is the EEG sampling frequency.

EEG data are very weak bioelectrical signals and factors such as ECG signals, electromagnetic waves generated by power components, and inherent noise from the acquisition equipment can cause interference during data collection (Ferracuti et al., 2021). Therefore, it is necessary to preprocess EEG signals to eliminate interference and improve the signal-to-noise ratio. In this study, EEG preprocessing (Sun and Mou, 2023) included baseline correction, 50 Hz notch filtering, 4–50 Hz bandpass filtering, independent component analysis (Zhukov et al., 2000), and removal of EOG artifacts (ElSayed et al., 2021).

The original EMG signals were recorded in this study as shown in Equation 2: Y = y 1 , y 2 , , y C 2 ∈ R C 2 × L 2 (2) Where C 2 represents the number of channels for EMG acquisition ( C 2 = 4 in this study), L 2 denotes the number of samples ( L 2 = T × f E M G ), T denotes the sampling time, and f E M G is the EMG sampling frequency.

Similar to EEG signals, EMG signals also require preprocessing to remove noise. This study employed a bandpass filter from 30 to 150 Hz to eliminate the influence of ECG signals, along with a 50 Hz notch filter to remove specific frequency noise.

As shown in Figure 5, this process consists of two steps: encoding and feature extraction of the bioelectrical signals. First, the collected bioelectric signals were encoded according to Equation 3 and Equation 4, to obtain X ψ , ψ ∈ E E G , E M G . Then, X ψ was input into the VAE network for feature extraction, resulting in relevant deep features Z ψ , ψ ∈ E E G , E M G .

Below is a detailed description of the process for extracting deep features from EEG and EMG signals using VAE. Let the input EEG and EMG signals be denoted as X , with its true distribution probability represented as P θ X . The relevant latent features Z are obtained from X using a VAE, and the reconstructed data are denoted as X ′ . Then the inferential network can be expressed as Q ϕ Z | X , and the generative network can be expressed as P θ Z P θ X ′ Z .

Because the potential feature Z is unobservable and cannot be solved directly, Z is estimated by calculating the lower bound. The true posterior distribution that the inference network could not determine is denoted as P θ Z | X , and it is assumed that the inferred network Q ϕ Z | X is a known distribution, and Q ϕ Z | X is used to approximate the replacement of P θ Z | X , such that both the VAE inferred and generated networks are known distributions. To make Q ϕ Z | X as close to P θ Z | X as possible, KL divergence is used to measure the degree of similarity between the two, and is minimized by optimizing the constraint parameters θ and ϕ . It can be expressed as shown in Equation 5: ϕ , θ = arg   min   D K L Q ϕ Z | X P θ Z | X = E Q ϕ Z | X 1 b Q ϕ Z | X − 1 b P θ Z | X + 1 b P θ X (5)

Noted as Equation 6: L θ , ϕ ; X = E Q ϕ Z | X 1 b Q ϕ Z | X − 1 b P θ Z | X (6)

It can be obtained as Equation 7: 1 b P θ X = D K L / ( Q ϕ Z | X P θ Z | X + L θ , ϕ ; X (7)

Since the KL divergence D K L * ≥ 0 is always true, therefore 1 b P θ X ≥ L θ , ϕ ; X is always true; thus, the variational lower bound function of the marginal likelihood 1 b P θ X of the input electrical signal X ψ can be obtained as L θ , ϕ ; X . Assuming Q ϕ Z | X follows a normal distribution, the optimization objective for the inferred network is Equation 8: ϕ , θ = arg   max ϕ , θ L θ , ϕ ; X (8)

To generate samples, the conditional distribution is generally a Bernoulli or Gaussian distribution whose probability density function is obtained using neural network computation. The optimization objective for the generative network is given in Equation 8.

The optimization objective of both the inferential and generative networks is to maximize the variational lower-bound function; therefore, the optimization objective of the VAE is to maximize the variational lower-bound function. Subsequently, the auxiliary parameter ε is introduced to transform the distribution of Q ϕ Z | X to obtain G ϕ ε , X with Z ψ = G ϕ ε , X , where ε ∼ p ε , and p ε have a known marginal likelihood distribution. Assuming that Q ϕ Z ψ | X ψ follows a normal distribution, the sampling of Z ψ can be described as Equation 9: Z i = μ i + σ i ε i (9) where μ represents the mean of the posterior probability and σ 2 represents the variance of the posterior probability. Subsequently, the variational lower bound can be simplified as Equation 10: L θ , ϕ ; X = − D K L Q ϕ Z | X P θ Z | X + 1 L ∑ l = 1 L 1 b P θ X i ′ | Z i , j = 1 2 ∑ 1 b σ i 2 − μ i 2 − σ i 2 + 1 + 1 b P θ X i ′ | Z i (10)

The introduction of auxiliary parameters changes the relationship between Z and σ and μ from sampling calculation to numerical calculation, which can be optimized using the stochastic gradient descent method (Dai and Wipf, 2019). The mean and standard deviation of P θ X i ′ | Z i can be calculated using neural network, enabling the calculation of each term in the lower bound of the variational inference. Consequently, the structure of the VAE model can be established.

To ensure that the classification results were close to the actual situation in the MI classification task, cross-entropy was adopted as the loss function. Cross-entropy loss was first introduced by Rubinstein in 1999 (Rubinstein, 1999), and it measures the discrepancy between the predicted and true distributions through simple calculations. The calculation is shown in Equation 15: L c r o s s − e n t r o p y y , y ^ = − ∑ k = 1 k y k ⁡ log ⁡ y ^ k , k = 0 , 1 (15) where y and y ^ represent the true label and classification model prediction probabilities, respectively. The result of the MI classification task is the type with the highest probability.

In feature metric learning as the auxiliary task, the triplet-center loss (TCL) (He et al., 2018) was used as the loss function. TCL combines the advantages of triplet loss for interclass relationship processing with those of center loss for intraclass relationship processing and reduces certain computational complexities to effectively improve the accuracy of feature metric learning, which is calculated as in Equation 16 L T C L = ∑ i = 1 M ⁡ max   D f i , c y i + m − min j ≠ y i D f x i , c j , 0 (16) Where, y i represents the label of the i -th type, f · represents the feature vector extracted in the neural network, f i represents the feature vector extracted by f from the i -th class, c j represents the center of the j -th class, and D · is the euclidean distance function is expressed as Equation 17: D f x i , c y i = 1 2 f x i − c y i 2 2 (17)

In this study, the abovementioned loss functions L c r o s s − e n t r o p y and L T C L and the loss functions L V A E − E E G and L V A E − E M G of the VAE network were combined to obtain the final loss function: L t o t a l = λ 1 L c r o s s − e n t r o p y + λ 2 L V A E − E E G + λ 3 L V A E − E M G + 1 − λ 1 − λ 2 − λ 3 L T C L (18) where, λ 1 , λ 2 , and λ 3 represent the hyperparameters for the contributions of models L c r o s s − e n t r o p y , L V A E − E E G , and L V A E − E M G , respectively, satisfying the condition λ 1 + λ 2 + λ 3 ≤ 1 . To determine the hyperparameters in Equation 18, a preliminary study was conducted using a grid-search algorithm. The sets of parameters λ 1 , λ 2 , and λ 3 in the grid search algorithm were defined as {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7}. The experimental results indicate that λ 1 , λ 2 , and λ 3 should be set to 0.7, 0.1, and 0.1, respectively. In addition, to complete the multitask ablation experiments, the parameters for two-task and three-task collaborative training were investigated. The parameters for the two-task model 2M-hBCINet-CEEGR were λ 1 = 0.7 and λ 2 = 0.3 , while the parameters for the model 2M-hBCINet-CEMGR were λ 1 = 0.7 and λ 3 = 0.3 . The parameter for the model 2M-hBCINet-CF was λ 1 = 0.7 . For the three-task model 2M-hBCINet-CEEMGR, the parameters were λ 1 = 0.6 , λ 2 = 0.2 , and λ 3 = 0.2 . The parameters for the model 2M-hBCINet-CEEGRF were λ 1 = 0.6 and λ 2 = 0.2 , and for the model 2M-hBCINet-CEMGRF, the parameters were λ 1 = 0.6 and λ 3 = 0.2 .

Common machine learning evaluation metrics were selected for this study to serve as evaluation indicators, with the performance of the 2M-hBCINet model described from various perspectives. The evaluation indicators included accuracy, precision, recall, and F1-score. The calculation methods are shown in Equations 20–23: Model performance was assessed by plotting the receiver operating characteristic (ROC) curve and calculating the area under the curve (AUC), known as the ROC-AUC. A c c u r a c y = T P + T N / T P + T N + F P + F N (20) P r e c i s i o n = T P / T P + F P (21) R e c a l l = T P / T P + F N (22) F 1 − s c o r e = 2 × P r e c i s i o n × R e c a l l / P r e c i s i o n + R e c a l l (23) where T P , F P , T N , and F N represent true positive, false positive, true negative, and false negative, respectively.

The LOOCV method (Walter et al., 2013) was employed to validate the performance of the 2M-hBCINet model. In the LOOCV method, one sample is extracted from the dataset to serve as the validation set for testing at each iteration until all samples are used as the validation set. The number of validations correspond to the number of samples, and the validation results are expressed as the average of all experimental outcomes, as indicated by Equation 24. The advantages of the LOOCV method include its ability to fully utilize data in small sample sizes for MI and to effectively prevent overfitting, thereby enabling a comprehensive assessment of the model’s generalization ability. e M = 1 N ∑ n = 1 N e r r M 2 M − h B C I N e t x n , y n M ∈ A c c u r a c y , P r e c i s i o n , R e c a l l . F 1 − s c o r e , R O C − A U C (24) Where N represents the number of samples, n represents the index of the selected sample, x n represents the samples, y n represents the corresponding labels of the samples, M represents the selected evaluation indices, e r r M · denotes the performance assessment results of the 2M-hBCINet model using different evaluation indices, and e M represents the validation results obtained for the corresponding evaluation metrics using the LOOCV method.

This study used parameter settings for three models: V A E E E G , V A E E M G , and total loss function. A symmetrical structure consisting of four convolutional layers was utilized for encoding and decoding in both V A E E E G and V A E E M G . Each layer was composed of convolutional kernels of sizes 3 × 3 × 128 , 3 × 3 × 256 , 3 × 3 × 256 , and 3 × 3 × 512 , with a stride of 1. Hyperparameters λ 1 , λ 2 , and λ 3 were set to 0.7, 0.1, and 0.1, respectively, for the total loss function.

To validate the effectiveness of each module of the proposed 2M-hBCINet model, an ablation study was conducted using the MI-EEMG, WAY-EEG-GAL, and ESEMIT datasets. LOOCV was used as the validation method. The validated models included model V A E E E G − C B A without the module V A E E M G , model V A E E M G − C B A without the module V A E E E G , model V A E E E G − V A E E M G without the module C B A , and the complete 2M-hBCINet model.

In the MI-EEMG dataset, 70 MI classification experiments were conducted for each model (14 participants × 5 trials). Figure 7A illustrates the performance of each participant in this experiment, and Figure 8A presents the significance analysis results for each model. In the WAY-EEG-GAL dataset, 60 MI classification experiments were performed for each model (12 participants × 5 trials). Figure 7B shows the performance of each participant in the experiment, and Figure 8B shows the significance analysis results for each model. In the ESEMIT dataset, 130 MI classification experiments were performed for each model (26 participants × 5 trials). Figure 7C shows the performance of each participants in this experiment, and Figure 8C displays the significance analysis results for each model. The statistical results of the experiments for each model across all the datasets are presented in Table 1.

Table 1 shows that the accuracy of the models missing various modules decreased by 6.1%–10.3%, 5.7%–8.0%, and 7.1%–10.5% compared to the 2M-hBCINet model. The precision, recall, and F1-score also decreased, demonstrating the necessity of each submodule. Furthermore, in both the MI-EEMG and ESEMIT datasets, the model V A E E M G − C B A achieved better performance than V A E E E G − V A E E M G and V A E E E G − C B A . This indicates that the features obtained from EMG signals provide greater expressiveness for MI tasks and allow for better representation than features obtained from EEG signals.

In the significance analysis presented in Figure 8, the p-values are categorized into five levels: p ≥ 0.05 , p < 0.05 , p < 0.01 , p < 0.005 , and p < 0.001 ; statistical significance was set at p < 0.05 , with smaller p-values indicating greater significance and p ≥ 0.05 indicating that the difference between the two models was not statistically significant. A statistically significant improvement in the accuracy of the 2M-hBCINet model compared to other models missing various modules indicated a notable difference, and demonstrated the significant superiority of the 2M-hBCINet model.

To validate the effectiveness of the MTL approach in the proposed 2M-hBCINet model, comparative experiments were performed using the MI-EEMG, WAY-EEG-GAL, and ESEMIT datasets, with LOOCV as the validation method. The validated models included the single-task model 2M-hBCINet-C, the two-task models 2M-hBCINet-CEEGR, 2M-hBCINet-CEMGR, and 2M-hBCINet-CF, as well as the three-task models 2M-hBCINet-CEEMGR, 2M-hBCINet-CEEGRF, and 2M-hBCINet-CEMGRF, and the four-task model 2M-hBCINet. In this context, C, EEGR, EMGR, and F denote MI classification (primary task), EEG reconstruction, EMG reconstruction, and feature metric learning (auxiliary tasks), respectively. The dataset and number of experiments used in this study were consistent with those used in the ablation study. Figure 9 illustrates the experimental performance of the participants across different datasets, and Figure 10 presents the significance analysis results for each model, with A, B, and C corresponding to the three types of datasets. The experimental results are presented in Table 2.

Table 2 shows that the 2M-hBCINet model, which employs all tasks, achieved an average accuracy improvement of 1.2%–11.0%, 1.5%–10.5%, and 1.7%–10.0% with the MI-EEMG, WAY-EEG-GAL, and ESEMIT datasets, respectively, compared to other models. Additionally, the precision, recall, F1-score, and ROC-AUC evaluation metrics showed improvements, demonstrating the effectiveness of each task during the model training process. This enhancement is attributed to the strong correlation between the selected EEG reconstruction, EMG reconstruction, and feature metric learning tasks, and the primary task of MI classification. The collaborative training process allows for mutual reinforcement and complementary learning among these tasks. Overall, MTL effectively enhanced the generalization ability and accuracy of the 2M-hBCINet model.

To further validate the superiority of the 2M-hBCINet model across different EEG frequency bands, we categorized EEG signals into four frequency bands: theta (four to eight Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–50 Hz), and conducted comparative experiments for each band. The comparison models included the single EEG model V A E E E G − C B A , single EMG model V A E E M G − C B A , and reference EEG models Reference-Leeb (Leeb et al., 2011), Reference-Barachant (Barachant et al., 2012), and Reference-Hutchison (Hutchison et al., 2010).

Altogether, 420 experiments were conducted with the MI-EEMG dataset (6 models × 5 frequency bands × 14 participants), and the results are shown in Figure 11A. Figures 12A–D illustrate the performance of each participant across different frequency bands. For the WAY-EEG-GAL dataset, 360 MI classification experiments were performed for each model (6 models × 5 frequency bands × 12 participants), and the results are shown in Figure 11B. Figures 12E–H depict the performance of each participant across different frequency bands in this dataset. A total of 780 MI classification experiments were performed with the ESEMIT dataset for each learning method (6 models × 5 frequency bands × 26 participants), and the results are summarized in Figure 11C. Figures 12I–L illustrate the performance of each participant across different frequency bands.

From Figure 11, it can be observed that in the theta (four to eight Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–50 Hz) frequency bands, the 2M-hBCINet model achieved average accuracy improvements of 4.5%–13.6%, 1.9%–12.9%, 1.7%–11.2%, and 3.7%–12.8%, respectively, compared to the models V A E E E G − C B A , V A E E M G − C B A , Reference-Leeb, Reference-Barachant, and Reference-Hutchison for the MI-EEMG dataset. Considering the WAY-EEG-GAL dataset, the improvements were 4.6%–12.1%, 2.0%–11.4%, 0.5%–10.9%, and 0.6%–9.7%. For the ESEMIT dataset, the model demonstrated accuracy increases of 7.1%–14.3%, 3.5%–11.2%, 2.4%–11.7%, and 3.6%–10.5%. Overall, the 2M-hBCINet model exhibited the best performance across all datasets, confirming its superiority. The reason for the optimal performance of the 2M-hBCINet model across different frequency bands is its ability to extract more discriminative EEG features using unsupervised methods. In addition, unlike other models that rely solely on EEG, such as Reference-Leeb, Reference-Barachant, and Reference-Hutchison, the 2M-hBCINet model effectively integrates deep EMG features through the CAM, further enhancing its performance across different frequency bands. This resulted in superior MI classification capabilities across all EEG frequency bands.

In practical MI classification tasks, issues such as muscle fatigue and decreased muscle strength can occur in participants. In this study, a muscle fatigue experiment was designed to evaluate the performance of the 2M-hBCINet model under muscle fatigue conditions. Referencing a previously described method (Xi et al., 2021), we simulated the muscle strength state of participants experiencing fatigue by reducing the amplitude of the EMG signals (with an amplitude reduction resolution of 5%). The comparative models included a single EEG model ( V A E E E G − C B A ), a single EMG model ( V A E E M G − C B A ), and EMG classification models that used the same feature selection but different classifiers, namely, Reference-Leonardis (Leonardis et al., 2015) and Reference-Tang (Tang et al., 2014). A total of 525 experiments were conducted with the MI-EEMG, WAY-EEG-GAL, and ESEMIT datasets (5 models × 21 muscle strength states × 5 repetitions), and the visualization results are presented in Figure 13.

Figure 13 shows the trends in accuracy for each model across the MI-EEMG, WAY-EEG-GAL, and ESEMIT datasets. For models relying solely on EMG, such as V A E E M G − C B A , reference-Leonardis, and reference-Tang, the accuracy dropped to 71.6%–75.7%, 64.6%–77.7%, and 61.2%–69.9% when the muscle strength was 50%. At a muscle strength of 25%, the accuracy further declined to 62.3%–67.8%, 54.8%–60.2%, and 54.5%–57.9%. As the muscle strength continued to decrease, the accuracy of these EMG-only models stabilized at approximately 50%. In contrast, for the proposed EEG-EMG classification model, 2M-hBCINet, the accuracy was 86.1%–87.6% at a muscle strength of 50% and 82.6%–85.6% at 25%. When the muscle strength was further reduced, the accuracy of 2M-hBCINet stabilized at 80.1%–81.3%.