Global Cortical Network Distinguishes Motor Imagination of the Left and Right Foot

Conventional passive lower limb rehabilitation is suboptimal since the brain is not actively involved in the training. An autonomous motor imagery brain-computer interface (MI-BCI) could potentially improve rehabilitation outcomes. However, motor cortex regions associated with the individual feet are anatomically close to each other. This presents a difficulty in distinguishing the left and right foot MI during rehabilitation therapy. To overcome this difficulty, we extracted functional connectivity to measure the global cortical network via electroencephalography (EEG) signals. Fourteen spatial connections (P3-Fp1, P3-F3, P3-F7, P3-C3, T5-F7, T5-C3, T5-T3, Fp2-T5, Fp2-P3, T6-Fp2, T6-T4, Cz-Fp1, Cz-F7 and Fp2-F7) found across twelve subjects significantly differed between the left and right foot MI, evidencing nonlocalized brain activity during MI. Foot MI were distinguished using machine learning algorithms in terms of the time- and frequency-domain connectivities extracted from Pearson’s correlation, multivariate autoregression (MVAR), bandpass correlation, and partial directed coherence (PDC) models. The results showed that connectivity extracted by pairwise Pearson’s correlation could be distinguished with 86.26 ± 9.95%, while in the frequency-domain, the gamma band presented the best classification accuracy of 73.55 ± 17.11%. We attempted to simulate asynchronous real-time classification paradigms in order to evaluate the classification performance of connectivity features compared to common spatial pattern (CSP) and band power (BP). The results indicate correlation-connectivity has the best outcome, attaining an accuracy of 80.75 ± 9.51% in asynchronous classification.


I. INTRODUCTION
Rehabilitation after stroke aims to promote neuron plasticity in the affected motor cortex. Conventional lower limb rehabilitation usually requires multiple therapists to manually assist patients in moving their foot. This passive therapy might fail to produce the desired outcome since the brain is not actively involved in the training. Advancements in mechatronics have allowed for the implementation of robotic exoskeleton support systems for facilitating lower limb rehabilitation [1]- [5].
The associate editor coordinating the review of this manuscript and approving it for publication was Bin Liu . Active training with autonomous involvement has been made possible by monitoring electromyography (EMG) signals in the affected foot. Recent interest in upper limb rehabilitation has included electroencephalography (EEG) monitoring, including the proposal of several motor imagery (MI) based rehabilitation systems [6]- [8]. A clinical study showed substantial improvement of arm motor function in an MI rehabilitation group, compared to a passive control group [9]. However, the recently proposed lower limb rehabilitation systems were binary (rest vs. imagery) [10], [11]. Bipedal MI training was near to impossible because the foot motor regions are close to each other.
Brain activity during upper limb movement has been well studied throughout the years. However, it has been difficult to distinguish the MI of individual feet from the associated regions of the cortex. The motor cortex areas responsible for the activity in each foot are anatomically close to each other [12]. Functional magnetic resonance imaging (fMRI) has shown that the spatial distributions of blood-oxygen-level-dependent (BOLD) activity are similar between left and right foot MI tasks [13]. A functional near-infrared spectroscopy (fNIRS) study was conducted to differentiate between the left and right foot during both motor execution and motor imagery [14]. The team found highly similar fNIRS patterns between left and right lower limb activity and concluded there was difficulty in classifying foot-specific brain activity. Besides, an EEG study that investigated alpha band power from the Cz channel found no significant differences between unilateral foot MI [15]. It remains difficult to localize foot-associated areas using EEG due to its low spatial resolution.
Brain-computer interfaces (BCIs) are systems developed to allow direct communication between a human brain and a computer. Brain activity patterns acquired through fNIRS, fMRI, EEG, or electrocorticography (ECoG) are decoded into an input signal to manipulate external devices such as cursors [16]- [18], wheelchairs [19], [20] or exoskeletons [21]- [23]. BCIs have also been used to detect physiological states such as attention [24], drowsiness [25], [26] and various emotions [27]- [29]. To date, the three main techniques to acquire an interpretable signal for BCIs are steady-state visually evoked potential (SSVEP), P300 oddball stimulus, and MI somatosensory rhythms. Although both SSVEP-and P300based BCIs have achieved promising results, these systems require an external stimulus in order to operate. MI is the overt mental representation of a motor task without any muscle movement [30]. One great advantage of MI as a BCI is that an external stimulus is unnecessary. Due to the spatially distinct motor areas associated with the upper limbs, hand MI is widely used in MI-based BCI. Several studies have also proposed tasks such as tongue or both feet MI (for review refer to [31]). However, there are limited studies regarding the classification of left versus the right foot MI; it remains a challenge in MI-based BCI, due to the closeness of the left and right foot motor areas. Recent research has shown that individual foot MI, similar to hand MI, produces a rebound in the beta frequency band, which could be potentially used as a feature to classify foot laterality [15]. Theta band event related desynchronization (ERD) was found to significantly differ between the MI of left and right foot stepping [32]. Correspondingly, beta rebound in bipolar C1-Cz and C2-Cz EEG was computed from 9 subjects to classify left-right foot MI [33]; an average accuracy of 69.3% was reported for a synchronous BCI design and 57.4% for an asynchronous system. Empirical mode decomposition (EMD) extracted from 3 s of EEG recording yielded an accuracy of 83.82% in discriminating foot MI [34]. An effort to classify left-right foot laterality utilizing steady-state somatosensory evoked potentials (SSSEPs) reported an average accuracy of 72.4% [35]. However, the abovementioned paradigms often require extensive algorithm optimization and yielded lower than optimum distinguishability.
Brain connectomics is an emerging technique to study the global cortical networks between various brain regions. Generally, there are two main types of brain networks: anatomical networks and functional networks. The anatomical networks contain the mappings of the white matter tracts, where axons and dendrites in the brain are bridged by synapses. Functional networks, on the other hand, are estimated from statistical models applied to brain signals recorded using functional imaging modalities such as fMRI or EEG. Functional brain connectivity has revealed differences in the default mode network (DMN) in the brain for various neuropsychiatric disorders, such as Alzheimer's disease [36], Parkinson's disease [37], epilepsy [38], attention deficit hyperactivity disorder (ADHD) [39], mood disorders [40], and schizophrenia [41]. Recent findings on EEG functional networks derived from MI tasks were discussed in Hamedi et al. [42]. Statistical models such as phase locking value (PLV), magnitude squared coherence (MSC), Pearson's correlation, multivariate autoregression (MVAR), partial directed coherence (PDC), directed transfer function (DTF), and transfer entropy have been applied to compute functional connectivity from hand, foot, and tongue MI. However, to our knowledge, the global cortical network between the left and right foot MI remains unexplored. Global cortical networks could potentially describe a more complete picture of EEG functional interactions during foot MI and improve lower limb MI-BCI performance.
The main objectives of this study are as follows: 1) To examine the global cortical network between left and right foot motor imagery. 2) To distinguish the global cortical networks of left and right foot MI by utilizing machine learning algorithms. 3) To evaluate the classification performance of network features in comparison to benchmark MI-BCI features, such as band power (BP) and common spatial pattern (CSP), in a simulated asynchronous experiment.

II. GLOBAL CORTICAL NETWORK EXTRACTION AND CLASSIFICATION
Global cortical networks are measured by the functional connectivity among the brain's cortical sites. We extracted functional connectivity from 12 healthy subjects undertaking foot MI tasks. The connectivity estimates derived from Pearson's correlation and MVAR were used as features to train classifiers. Furthermore, we simulated an asynchronous experiment by a sliding window approach to compare the classification performance between connectivity, BP and CSP features.

A. EEG DATASET
The EEG dataset used in this study was published by Kaya et al. [43], with the following EEG setup parameters. VOLUME 8, 2020 EEG signals were recorded using an EEG-1200 JE-921A EEG system (Nihon Kohden, Japan) based on standard 10/20 international configuration. The sampling frequency was 200 Hz and the impedance was below 10 k throughout the recording sessions. Two earbud electrodes A1 and A2 were used as the ground, while 'System 0 V' was used as the virtual reference point and 'System 0 V' was described in the EEG-1200's technical manual as 0.55 × (C3 + C4)V .
No custom filtering was applied to the recorded signal, except a 50 Hz notch filter, to reduce electrical noise. The dataset consists of five paradigms recorded from healthy individuals, within which was the hand/leg/tongue (HaLT) MI paradigm. A total of six mental tasks including left hand, right hand, left foot, right foot, tongue, and rest were carried out. During the recording, a stimulus appeared on screen for 1 s, which prompted subjects to imagine the corresponding motor movement within the simulation period. After each trial, subjects were asked to fix their gaze on a fixation point for 1.5-2.5 s. Each subject implemented a total of 143-156 trials per task.
From this dataset, we utilized the HaLT MI paradigm recording for our analysis. We omitted 2 occipital signals (out of 19 EEG channels) due to minimal relevancy to the motor functions of interest. Since the original virtual reference was at the sensorimotor cortex, we rereferenced the signals in order to retain EEG signals during MI. The unipolar physical reference was usually selected from the EEG channels located at the midline (Fz, Cz and Pz) [44]. However, both motor areas (Cz) and the prefrontal cortex (Fz) were found to contribute to the MI [45]- [47]. Hence, we rereferenced the EEG signals corresponding to the Pz channel, resulting in 16 remaining channels (Fp1, F3, F7, C3, P3, T3, T5, Fz, Cz, Fp2, F4, F8, C4, P4, T4, and T6). Since the EEG signals were validated technically in the original paper, no additional preprocessing such as source localization was carried out as it could disrupt the informative connections presence. For our in-depth analysis of foot MI, we selected only the left foot and right foot MI trials from the 12 subjects for downstream analysis.

B. EXTRACTION OF FUNCTIONAL CONNECTIVITY
Pearson's correlation, multivariate autoregressive (MVAR), and partial directed coherence (PDC) models were used to extract the brain-wide weighted connectivity networks between signals recorded from EEG electrodes. Pearson's correlation, r, measures the pairwise correlation of EEG signals between pairs of electrodes (Eq. 1). The connectivity measured by Pearson's correlation is structured to be weighted and undirected. Correlation can estimate the strength of connections (correlation coefficient), but not the direction of information flow; it can only measure the similarity between signal x and y by dividing their covariance by the product of their standard deviations.
wherex andȳ are the means of signals from channels x and y, respectively, while N t is the number of time samples.
On the other hand, MVAR estimates the weighted, directed connectivity in terms of Granger causality. MVAR computes the regression between signal y at time t with signal x at time t − , where x t− is signal x with lag (Eq. 2). A MVAR model can be described as the influence of lagged signal x t− to signal y t . MVAR with order 1 was utilized in this study in order to preserve the same number of features across different connectivity estimates.
where t represents the inclusion of white noise, L is the model lag-order, and |P| > 0 represents the connectivity exerted by signal x t− to signal y t . In multichannel applications, the MVAR model could be written as a multivariate linear regression, where β is a N ×N matrix consisting of the MVAR coefficient between EEG signals recorded from N numbers of channels. The connectivity coefficient β is commonly estimated from the conditional least square [48] as Unlike correlation and MVAR, which measure brain connectivity in the temporal domain, PDC calculates the functional connectivity in the frequency domain. PDC, as proposed by Sameshima and Baccala [49], estimates signal dependency in a specific frequency, f (Eq. 5).
where P(f ) is the Fourier transform of the MVAR coefficient.
, is a normalized ratio between the outflow of connections from channel x j to y i , and the total outflows from x j to all channels at frequency f .

C. COMPLEX NETWORK ANALYSIS
Recently, a complex network analysis based on graph theory has been proposed to characterize complex, high-dimensional brain connectivity [50], [51]. Complex network analysis reduces the dimension of features by computing meaningful information from brain connectivity networks. We extracted four complex network measures to characterize local and global network integration, as well as network segregation. Degree, k w i , measures the local network integration in individual channels (Eq. 6), while global efficiency, E w , summarizes brain-wide global network integration (Eq. 7).
where w ij is the weighted connectivity between EEG channels i and j, and d w ij is the shortest weighted path length between nodes i and j. Clustering coefficient, C w , and transitivity, T w , were computed to characterize the local (Eq. 8) and global (Eq. 9) network segregations, respectively.
where k i is the weighted degree as calculated in Eq. 6, and t i is the weighted mean of the triangle around node i, defined as t w i = 0.5 j,h∈N w ij w ih w jh 1/3 .

D. MACHINE LEARNING ALGORITHMS
In offline analysis, support vector machine (SVM) [52] was used to distinguish the N × N connectivity matrices and 34 × 1 complex network features extracted from left and right foot MI. SVM is a machine learning algorithm, widely used for classification due to its ability to maximize the margins between the classifier's separating hyperplane and support vectors. No thresholding approach was implemented in the connectivity matrices, to allow classifiers to learn from all features. We compared the classification accuracy of both linear and polynomial (order 2 & 3) SVMs, as well as linear discriminant analysis (LDA) and k-nearest neighbor (KNN; K = 3) for distinguishing connectivity features during imagery foot movement. We also evaluated the SVM performance in classifying connectivity features derived from different numbers of EEG channels, as well as EEG signals from various frequency bands. CSP and BP are two commonly used features to distinguish brain patterns of motor imagery. We designed a simulated paradigm to evaluate the asynchronous real-time classification of CSP, BP, and connectivity features. For the test set, we randomly stratified raw EEG signals from 20 MI trials (10 from each foot), and the remaining data were used to train the SVM. We utilized a sliding window, with length 1 s and stride of 0.1 s, across the test set. Features were extracted from each window and were classified in real-time; classification accuracy, correlation, and loss were computed as performance evaluation metrics.
In all classifications, except intersubject classifications, the machine learning algorithms were trained with ten-fold cross-validation (90% training data and 10% testing data). The random guessing threshold for two-class classification was 50%. The reported subject-specific accuracies are the average accuracies across ten folds, with corresponding standard deviation. In other to minimize subject-bias and foldbias, the mean standard deviations of accuracies, A were calculated by where σĀ is the standard deviation of mean accuracy, σ d s represents the standard deviation of each subject s across d-fold, and N s is the number of subjects. The results comparisons were conducted using Bonferroni corrected [53] pairwise nonparametric Wilcoxon signed rank tests [54].

III. OFFLINE ANALYSIS RESULTS
Functional connectivity measures were used to extract the interactions between various brain regions. Connectivity features were estimated not only from the somatosensory area, but from brain regions including frontal, parietal, and temporal areas. Machine learning algorithms were trained to classify the functional connectivity of the left and right foot MI. We compared the classification performance of time-domain (correlation & MVAR), frequency-domain (bandpassed correlation & PDC), and feature reduced (complex network) connectivity features. The paradigm of offline analysis is depicted in Figure 1.

A. CONNECTIVITY FEATURES EXTRACTION
Pearson's correlation was applied to estimate pairwise connectivity between EEG signals recorded from various scalp regions. The correlation coefficient, r xy ∈ [−1, 1], is the connection strength between brain areas, where a positive correlation, r xy > 0, represents mutual potential fluctuation in time-aligned EEG activity, while a negative correlation, r xy < 0, indicates the potential of signals x and y are inversely proportionate. A 16 × 16 symmetric connectivity matrix was constructed, with each element representing the correlation between a pair of signals. We performed a connectionwise, two-sample Student's t-test, with α = 0.05, to measure significant differences of the individual connections between the left and right foot MI. The α was adjusted with the Holm-Bonferroni [55] correction for multiple comparisons, as it is less conservative than the classical Bonferroni adjustment when the number of comparisons is high. Thresholding was not performed because spurious connections from neighboring electrodes might unlikely pass Student's t-test, as these connections were inherently present during left or right foot MI tasks, given that both tasks were carried out in a stratified random manner. Moreover, additional analysis revealed that classification performance deteriorated after applying a threshold to filter weak connections. This indicates that weak connections could contribute to the distinguishability of foot MI. Figure 2 visualizes the statistically significance of trial-averaged connections from individual subjects; the insignificant connections were set to zero. It is interesting to note that there are numerous significant positive (red) and negative (blue) connections observed across the subjects, indicating that foot MI could potentially be discriminated via the global cortical network. We noticed that all subjects exhibited significant connections to the fronto-central regions, particularly nine subjects who showed significant connections between the left motor cortex and the P3 region.
To quantify the subject independent connections, we investigated the inherent intersubject connections (ITSCs) that significantly differentiate left and right foot MI. Figure 3 shows fourteen inherent connections present in more than half of the subjects. We found stronger P3-Fp1, P3-F3, P3-F7, P3-C3, T5-F7, T5-C3, T5-T3, Fp2-T5, Fp2-P3, T6-Fp2 and T6-T4 connections in the left foot MI. The right foot MI exhibit stronger Cz-Fp1, Cz-F7 and Fp2-F7 connections compared to the left foot MI. Our results suggest an intercooperation between brain regions during MI tasks. Contrary to the popular belief that motor-related tasks only produce changes in motor areas, we found interdependency among the frontal, temporal, parietal, and motor regions during foot MI.

B. INTRASUBJECT CLASSIFICATION
Machine learning algorithms were utilized to distinguish foot MI from functional connectivity. Ten-fold cross-validation was performed. Subject-specific foot MI trials were randomly partitioned into ten equal portions containing the same amount of data from each class. Nine out of ten portions were used to train classifiers and the remaining one was used to evaluate the trained model. This process was repeated ten times until a trained classifier was evaluated exactly once using each of the ten portions. The mean accuracy and standard deviation across ten folds were reported for each individual subject.

1) TIME-DOMAIN CONNECTIVITY
Linear support vector machine was trained on both MVAR-based and correlation-based connectivity. The accuracies reported are the averages of ten-fold cross-validations. The results in Table 1 show that better performance was achieved by correlation-connectivity than MVAR-connectivity (p < 0.01), with average accuracies of 86.26 ± 9.95% and 81.48 ± 12.39%, respectively. Seven out of twelve subjects demonstrated a classification accuracy of more than 80% using correlation-connectivity as a feature, with only one outlier (Subject 7) performing lower than 70%. The outlier, for which the accuracy was more than two standard deviations of the mean, was excluded from downstream analysis. The results in Table 2 indicate linear SVM is more robust than nonlinear SVM in classifying connectivity features with high dimensional space; SVM tends to overfit easily on high dimensional training data with an increasing polynomial order. We also compared the performance of SVM with LDA and KNN, which showed suboptimal performance compared to linear SVM. However, the performances of linear SVM and LDA were not significantly difference, suggesting that linear classifiers were better in discriminating foot MI connectivity. We further investigated the classification performance of the N × N correlation-connectivities, measured from different sets of EEG electrodes (namely, 3E, 8E 16E and ITSCs). • ITSCs (14 connections): Fourteen inherent intersubject connections as described in Section III-A. The first electrode set covers only the sensorimotor areas while the frontal and parietal areas were taken into account in 8E. Based on these results reported in Table 3, functional connections solely from the motor cortex (3E) were insufficient to discriminate left and right foot MI. The performance   drastically improved when frontal and parietal connections were incorporated. The intersubject connections showed the best accuracy-to-feature ratio, achieving a mean classification accuracy of 79.77 ± 11.92% with only 14 connections, with five out of eleven subjects achieving classification accuracies of more than 80%.

2) FREQUENCY-DOMAIN CONNECTIVITY
We extracted frequency-specific connectivity features using two approaches: bandpass correlation and partial direction coherence (PDC). We applied a minimum-order, linearphase, finite impulse response (FIR) filter to bandpass EEG signals into five frequency bands, delta (1-3 Hz), theta (4-7 Hz), alpha (8)(9)(10)(11)(12)(13), beta (14-29 Hz) and gamma . Band-pass connectivity was extracted from bandpassed signals using pairwise Pearson's correlations, via eq. 1. On the other hand, PDC was estimated from the MVAR coefficients, as shown in eq. 5. Table 4 and Table 5 summarize the performances of frequency-domain connectivity in classifying left and right foot MI. It was surprising to note that the gamma band achieved the best performance among the frequency ranges, despite major studies that have classified MI using band power from low frequency ranges [56]- [58]. Differences of > 4.86% and > 1.71% were noted between the performance of the gamma band and the remaining frequency bands, in bandpass correlation and PDC connectivity, respectively. However, the differences did not pass the Wilcoxon signed rank test. Frequency-domain connectivity performed poorer than time-domain connectivity (p < 0.05), suggesting the MI discriminating features are highly represented in the time-domain connectivity.

3) COMPLEX NETWORK FEATURES
Complex network analysis was used to extract high hierarchical network topology from the connectivity matrices. Connectivity networks often carry complex, high-dimensional features that are difficult to interpret. In the machine learning domain, model overfitting can occur during training with high-dimensional features. We applied complex network analysis to the 16 × 16 correlation-connectivity network, to summarize the 34 high-order topology features, comprising of 16 channelwise degrees, 16 channelwise clustering coefficients, 1 brain-wide global efficiency metric, and 1 brain-wide transitivity metric. Figure 4 shows the SVM classification performance using these complex network features from the eleven subjects. An average accuracy of 81.94 ± 11.32% was achieved. Although the accuracy dropped slightly in comparison to the correlation-connectivity (86.28 ± 9.95%), complex network features performed better than the ITSCs features selected manually from the EEG electrodes (79.77 ± 11.93%).

C. INTERSUBJECT CLASSIFICATION
Major efforts in BCI research concentrate on model generalizability. An ideal BCI system should be able to generalize its performance to different individuals. However, to date, there is no single model that is capable of fitting all subjects. We evaluated the generalization of the correlation-connectivity model using a leave-one-subject-out (LOSO) paradigm. An SVM classifier was trained with features from all subjects except one. The excluded subject was then used as the testing set to assess the model. This process was repeated until the feature from each subject was left out exactly once; the classification results are shown in Figure 5. The average accuracy was 67.02 ± 11.91%, with four out of eleven subjects having a performance better than 70%.
Our results indicate that intersubject variability is present in the functional connectivity, which could potentially hinder the performance of a BCI.

IV. SIMULATED REAL-TIME ANALYSIS
Asynchronous BCI could provide a more natural training paradigm for poststroke patients. We evaluated the classification performance of connectivity features compared to benchmark MI-BCI features, such as band power (BP) and common spatial pattern (CSP), in a simulated asynchronous experiment. BP from five main EEG bands was extracted from all recording channels with ranges: delta (1-3 Hz), theta (4-7 Hz), alpha (8)(9)(10)(11)(12)(13), beta (14-29 Hz) and gamma . On the other hand, the two most discriminating CSP features were computed for this study. The asynchronous real-time paradigm was simulated by randomly aligning the raw EEG signals of 20 MI trials (10 from each foot). The features extracted from the remaining trials were used to train the SVM classifier. A 1 s sliding window, with a 0.1 s stride, was applied on the aligned raw signals, and the MI features were extracted from each window for classification, as depicted in Figure 6. These procedures were repeated ten fold, by which in each fold, 20 random nonrepeated trials were aligned to evaluate the SVM model. We compared the performances among CSP, BP, and the connectivity features, in terms of classification accuracy, correlation, and loss, as shown using the best performing subject (Subject 9) in Figure 7. The accuracy and correlation were calculated between the true and predicted class labels, while the prediction loss was the root mean square error (RMSE) between the true label and posterior probability of the SVM classifier. The individual performance is tabulated in Table 6. Both MVARand correlation-connectivity features performed better than CSP and BP in asynchronous foot MI classification. For the two connectivity measures, correlation-connectivity reported the best results, with a mean accuracy of 80.75 ± 9.51%, correlation of 0.62 ± 0.19, and loss of 0.36 ± 0.08. Nine subjects successfully achieved real-time classification accuracy above 70%, for which Subject 9 performed best with an accuracy of 92.46%, correlation of 0.85 and loss of 0.24. These results suggest that foot MI could be classified in the majority of subjects by utilizing an EEG-based global cortical network as features.

V. DISCUSSION
Based on our findings, left foot MI exhibits stronger P3-Fp1, P3-F3, P3-F7, P3-C3, T5-F7, T5-C3, T5-T3, Fp2-T5, Fp2-P3, T6-Fp2 and T6-T4 connections. On the other side, right foot MI shows higher connectivity strength in Cz-Fp1, Cz-F7 and Fp2-F7. Classification of these fourteen intersubject connections shows promising accuracy, indicating brain activity associated with foot-based motor tasks is not limited solely to the responsible motor cortices. We found numerous frontal, parietal and temporal connections significantly differ between left and right foot MI. The contribution of frontal brain region activity during foot tasks has been studied using fMRI. In Sahyoun et al. [59], higher BOLD activity was found in both the frontal and motor cortex during the preparation of foot movement; it was proposed that anterior prefrontal regions are involved in decision making during movement preparation. The parietal cortex has been found to be involved in motor planning [60], [61]. Upper limb motor imagery tasks were successfully classified from the parietal cortex of a tetraplegic patient using implanted electrodes [62]. Our findings suggest that temporal cortex connectivities are involved in foot MI as well. Temporal lobe functions encompass memory and learning, and their activity during MI might indicate subjects practicing to perform MI tasks. By integrating the connections among various brain regions, we developed a more complete description of the functional interactions of the brain's electrical activity during left and right foot MI tasks.
Our results showed better performance in undirected correlation-connectivity than directed MVAR-connectivity. MVAR calculates directed connectivity by Granger causality. By using the least-square method, the 16 × 16 connections were estimated from 1 s (200 time points) segments of EEG signals. With a large number of estimation parameters and a short time series, the model suffers from a low degree of freedom (DOF), which might not produce a reliable estimation. MVAR models are robust for EEG signals that have a    long recording span, but are suboptimal in rapid MI tasks. On the other hand, Pearson's correlation estimates pairwise connectivity. A single correlation parameter was estimated from each pair of 1 s EEG signals, resulting in a model with higher DOF and better estimation. This could explain the performance variation between the two connectivity features. In addition, we showed features from the sensorimotor cortex are insufficient to discriminate left and right foot MI; inclusion of features from the frontal, parietal and temporal cortex improved the discriminating power. We observed the gamma frequency band as the most discriminatory frequency range, even though conventionally BCI research has rarely adopted gamma activity as features. More recently, studies have found strong relationships between gamma activity and motor functions [63]- [65], which our results suggest their potential extensibility to motor imagery of the foot.
However, the global cortical network model developed still suffers from low generalization across subjects, which could potentially be improved in the future with implementation

FIGURE 7.
Simulated asynchronous real-time classification using different EEG features. The figure shows the results from the best performing Subject 9 in a single fold. Accuracy and correlation were calculated from the true and predicted labels. RMSE between the classifier posterior probability and true label was computed as the model loss.
of transfer learning techniques. Recent work has employed transfer learning to confront intersubject variability in BCI [66]. Additionally, deep learning-based adversarial variational autoencoders transfer learning has been implemented with promising results [67]. Such transfer learning techniques could be applicable to foot MI BCI and utilized to improve cross-subjects classification performance.

VI. CONCLUSION
We proposed the use of global cortical network measures to study distinct brain activity between left and right foot MI. We found significant variations in the frontal, parietal, temporal, and motor cortex connections between the bilateral tasks, where left foot MI exhibits stronger P3-Fp1, P3-F3, P3-F7, P3-C3, T5-F7, T5-C3, T5-T3, Fp2-T5, Fp2-P3, T6-Fp2 and T6-T4 connections and right foot MI exhibits stronger Cz-Fp1, Cz-F7 and Fp2-F7 connections. We achieved a noteworthy accuracy of 86.26±9.95% for classifying correlationfunctional-connectivity of foot MI using linear SVM. Gamma frequency bands outperformed other EEG bands in distinguishing foot MI, suggesting the potential application of gamma connectivity in BCI systems. Complex network measures also showed promising discriminative power, reaching a mean classification accuracy of 81.94 ± 11.32%. Correlation-connectivity performed best in simulated asynchronous real-time classification paradigms, in comparison to MVAR-connectivity, CSP, and BP features. These findings inform which global cortical network features could provide an effective basis for an asynchronous BCI system for lower limb rehabilitation. A recently proposed mixed-reality (MR) music rehabilitation system introduced a novel scheme for poststroke lower limb training [68]. Our framework has the potential to integrate with the MR system to provide an intrinsic rehabilitation paradigm.
CHUN-REN PHANG is currently pursuing the Ph.D. degree with the International Ph.D. Program in Interdisciplinary Neuroscience, Neural Engineering Laboratory (NEL), National Chiao Tung University, Taiwan. His primary research interests include neural engineering, brain-computer interface, statistical analysis of biosignals, machine learning, and neural networks.
LI-WEI KO (Member, IEEE) is currently a Professor with the Institute for Bioinformatics and Systems Biology, National Chiao Tung University, Taiwan. He is also the Director of the Interdisciplinary Neuroscience Ph.D. Program, National Chiao Tung University. His primary research interests include neural engineering, including brain-computer interfaces, neural rehabilitation, and computational neuroscience, and especially for developing the mobile and wireless brain-machine interface for daily life applications. He is the Chair of IEEE Computational Intelligence Society Taipei Chapter. He has served as an Associate Editor for the IEEE TRANSACTIONS ON NEURAL NETWORKS and LEARNING SYSTEMS, from 2010 to 2015. He has been serving as an Associate Editor for the IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, since 2019. VOLUME 8, 2020