Large-Scale Cortical Network Analysis and Classification of MI-BCI Tasks Based on Bayesian Nonnegative Matrix Factorization

Motor imagery (MI) is a high-level cognitive process that has been widely applied to clinical rehabilitation and brain-computer interfaces (BCIs). However, the decoding of MI tasks still faces challenges, and the neural mechanisms underlying its application are unclear, which seriously hinders the development of MI-based clinical applications and BCIs. Here, we combined EEG source reconstruction and Bayesian nonnegative matrix factorization (NMF) methods to construct large-scale cortical networks of left-hand and right-hand MI tasks. Compared to right-hand MI, the results showed that the significantly increased functional network connectivities (FNCs) mainly located among the visual network (VN), sensorimotor network (SMN), right temporal network, right central executive network, and right parietal network in the left-hand MI at the $\beta $ (13-30Hz) and all (8-30Hz) frequency bands. For the network properties analysis, we found that the clustering coefficient, global efficiency, and local efficiency were significantly increased and characteristic path length was significantly decreased in left-hand MI compared to right-hand MI at the $\beta $ and all frequency bands. These network pattern differences indicated that the left-hand MI may need more modulation of multiple large-scale networks (i.e., VN and SMN) mainly located in the right hemisphere. Finally, based on the spatial pattern network of FNC and network properties, we propose a classification model. The proposed model achieves a top classification accuracy of 78.2% in cross-subject two-class MI-BCI tasks. Overall, our findings provide new insights into the neural mechanisms of MI and a potential network biomarker to identify MI-BCI tasks.


I. INTRODUCTION
M OTOR imagery (MI) is a multidimensional high-level cognitive simulation process [1].It is a mental rehearsal of specific motor acts without any motor output [2].Over the past two decades, MI has been widely used to control intelligent devices, provide a quality life for disabled patients, and improve brain-computer interface (BCI) systems [3].Besides, MI plays a key role in studying motor-related brain cognitive functions.However, due to individual differences, the decoding of MI tasks still faces challenges, and the large-scale brain network neural mechanisms underlying its clinical application are unclear, which seriously hinders the development of MI-based clinical applications and BCI systems.Thus, it is necessary to profoundly understand the brain network neural mechanism of MI and explore new network biomarkers to identify MI tasks.
In recent years, large-scale brain network analysis methods [4], [5] based on electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) neuroimaging data have provided a new way to understand the neural mechanism of MI.For example, Zhang et al. [6] investigated the reconfiguration patterns of large-scale networks in distinct MI tasks using fMRI and revealed that somatomotor network and dorsal attention network were important in switching from resting-state to task state.Based on EEG analysis, Li et al. [7] probed the large-scale functional network connectivity (FNC) from resting to MI task states, and they reported that the connections between bilateral motor areas (i.e., premotor cortex and supplemental motor area, SMA) were increased in MI compared to resting-state.Leeuwis et al. [8] found that MI-BCI participants performed better when the network functional connectivity in the right hemisphere was enhanced.These existing studies consistently suggest that network connectivity analysis may provide a valuable characteristic feature in understanding MI illiteracy and MI classification.Thus, large-scale network analysis based on EEG and fMRI data can help us better understand the neural mechanism of MI.
In practical applications, EEG-based MI can be more conveniently and economically used in clinical rehabilitation and BCI control [9], [10].However, EEG signal has a drawback, that is, its low spatial resolution and is easily contaminated [11], making it difficult to construct largescale networks.To overcome these shortcomings and perform large-scale network analysis, some brain source localization methods have been proposed to improve spatial resolution, especially the standardized low-resolution electromagnetic tomography (sLORETA) method which offers better localization results [12].In our study, the sLORETA method is used for EEG source reconstruction.Based on the high spatial resolution EEG source localization data, the group independent component analysis (ICA) and Bayesian nonnegative matrix factorization (NMF) methods were commonly used for temporal-spatial decomposition to recognize large-scale networks.ICA is widely used as a powerful and data-driven tool for evaluating the underlying spatiotemporal structure using brain imaging data [13].Bayesian NMF is an unsupervised learning method for improving the interpretability of latent variable decompositions on brain imaging data [14].According to a previous study, Sockeel et al. [15] combined source localization and spatial ICA to identify large-scale networks in resting-state and indicated that group ICA could characterize temporally and spatially brain networks of EEG due to its high temporal resolution.Moreover, Bayesian NMF has the advantages of excellent identification capabilities and soft-partitioning in overlapping community detection [16].However, ICA employs a firm restriction of temporal and spatial independence on the components and the negative components provided by ICA are not physically interpretable [17].Thus, we combined source reconstruction and Bayesian NMF to investigate the large-scale network patterns of MI.
Furthermore, there are still many challenges in EEG-based MI-BCI identification, especially cross-subject classification.Therefore, deep learning methods are widely used in the classification of MI.For example, Huang et al. [18] introduced the local reparameterization trick into convolutional neural networks (CNN) to decode four-class MI-BCI tasks resulting in 92.41% classification accuracy.Although the above deep learning study achieves relatively high classification accuracy, it mainly concentrates on the subject-specific strategy and may have high time costs.It is known that discriminative features and appropriate classifiers can help improve model classification accuracy and reduce calculation time.There has been growing evidence suggesting that large-scale brain network analysis methods can provide a new way to explore discriminative network biomarkers to decode MI [6].For example, Gu et al. [19] extracted features from FNCs by phase locking value (PLV) and graph theory indices using support vector machine (SVM) achieving 75% accuracy for discerning left-hand MI and right-hand MI tasks.Ai et al. [20] extracted features from network properties based on FNCs to distinguish four-class MI tasks achieving 79.7% accuracy.Thus, we used the large-scale network patterns as features to construct the classification model of MI tasks, which may contribute to identifying new biomarkers to distinguish multiple MI tasks and developing efficient MI-BCI systems for engineering and clinical applications.
In this study, we combined EEG data source reconstruction and Bayesian NMF methods to construct large-scale networks of MI-BCI tasks at different frequency bands.We first compared the FNC and network property differences between the left-hand MI and right-hand MI at the α (8-13Hz), β (13-30Hz), and all (8-30Hz) frequency bands.Then, we assessed the relationships between large-scale network patterns and MI-BCI performance.Finally, based on the largescale networks, the spatial pattern of networks (SPN) was used to extract the combined spatial patterns as features to construct MI-BCI classification models, where three widely used classifiers were selected, such as SVM, linear discriminant analysis (LDA), and light gradient boosting machine (LightGBM).The proposed method captures brain information with large-scale network features, aiming to reduce the impact of individual differences and to be available and effective for practical BCI applications.

A. Participants
The public Motor Movement/Imagery EEG dataset can be downloaded at https://physionet.org.In the experiment, 109 healthy subjects participated in the multiple motor/ imagery tasks.S088, S092, S100, and S104 were excluded due to their damaged recordings, so 105 subjects were used in our study.The EEG dataset was recorded by 64-channel (10-10 international system) with a sampling rate of 160 Hz and the predefined bandpass filter was 0-80 Hz.

B. MI-BCI Experiment Procedure
The EEG experimental paradigm contains 14 runs, including two one-minute baseline runs individually (eyes open and eyes closed), and three two-minute runs individually (open and close left or right fist, imagine opening and closing left or right fist, open and close both fists or both feet, and imagine opening and closing both fists or both feet).Our study used the tasks of imaging opening and closing left or right fist tasks (corresponding to runs 4, 8, and 12).For the MI tasks, each run contains 15 trials, resulting in 45 trials for each subject.During the experiment, a cue appeared and stayed on the left or right side of the screen, and the subject imagined the corresponding fist opening and closing.Then the subject could rest until the next trial.In the experimental paradigm, -2-0s corresponds to the pre-MI resting period, 0-4.1s corresponds to the MI tasks period, and 4.1-6.1scorresponds to the after-MI resting period.Fig. 1 shows the detailed experimental paradigm.

C. MI-EEG Data Preprocessing
The MI-EEG dataset of 105 subjects was processed with a standard procedure.In this paper, we extracted frequency bands the α (8-13Hz), β (13-30Hz), and all (8-30 Hz) from the EEG signal using band-pass filtering in the finite impulse response filter.Then, the average reference method was used to eliminate the effects caused by variations in the original reference electrodes.EEG signals were segmented into epochs of 6s (from -2s to 4s) duration, and baseline correction was performed by subtracting the average value on the pre-MI 200ms period (from -0.2s to 0s).Next, the artifact trail removal set ±100µV as the threshold, where 18 subjects were removed leaving 87 subjects.Eye movement and blink artifacts were removed by ICA.Finally, we selected the 4s (from 0s to 4s) data and averaged trials for each subject in each task.The above processings are based on the MNE of Python software.

D. Analysis Procedure
Fig. 2 shows the analysis procedure.First, the raw EEG data were preprocessed and concatenated.Then, we reconstructed cortical sources from EEG scalp using the sLORETA software.Bayesian NMF was utilized to identify large-scale networks related to MI at group level.Furthermore, we construct FNC and calculate the network properties for the left-hand and right-hand MI tasks.We also compared the FNC and network property differences between the two tasks and revealed the relationships between the network features and BCI performance.Finally, the obtained network properties and large-scale FNC were used to construct a classification model to discern the two tasks.The specifics of these steps are described in the following sections.

E. Source Reconstruction of MI-EEG Data
Reconstructing deep cortical source activity of EEG accurately is widely used to measure neural activity of scalp potentials and decode cognitive processes [21].The key to the process is a solution to inverse problem, which can provide significant information on the time course and brain function localization [22].In our study, we have 87 subjects (2 trials per subject) and we concatenated the left-hand and right-hand MI data along the time points to 87 × 2×640 = 111360, resulting in a 2D matrix of 64 × 111360.Then, the LORETA software (v20210701) was used to reconstruct brain sources from EEG scalp.Three-shell spherical head model was used to estimate brain activities in the sLORETA.The standardized boundary element method was used for head modeling.Reconstructed cortical brain activations matched to Montreal Neurological Institute space.Finally, the 3D space contained 6239 voxels at a 5mm spatial resolution.For the 64-electrode system, the dimension of the lead field matrix was 64 × 3×6239.After source reconstruction, the data was transformed into the matrix of 6239 × 111360.Furthermore, we transposed the matrix to 111360 × 6239 before Bayesian NMF estimation.

F. Bayesian NMF
Bayesian NMF derived Gibbs sampler to approximate the posterior density of the NMF factors, based on the normal likelihood and exponential priors.The NMF problem [23] can be formulated as X = AB + E, where X ∈ R I ×J needs to be decomposed into two nonnegative matrices A and B. A ∈ R I ×N + and B ∈ R N ×J + with R + indicating the nonnegative real number.E ∈ R I ×J represents the residual matrix.Based on the Bayesian procedure, we assumed that the residuals E i, j are i.i.d.zero mean normal with variance σ 2 , which yields the likelihood as follows: where )) is the normal density, and θ = {A,B,σ 2 } are all parameters.
Matrices A and B are assumed to be exponentially distributed respectively, with scales α i,n and β n, j : where ε(x; λ ) = λ exp(-λ x)u(x) is the exponential density, and u(x) is the unit step function.An inverse gamma density with shape k and scale θ is selected as the prior for the noise variance: Based on Gibbs sampling, an efficient Markov chain Monte Carlo method was used to evaluate its posterior density.The conditional densities were denoted by R(x; µ, σ 2 , λ ) ∝ N(x;µ, σ 2 )ε(x; λ ).The conditional density of A i,n is as follows: where A i,n indicates all elements of A without A i,n , and the expression for B follows from the symmetry.The conditional density of σ 2 is an inverse-gamma: Then, the posterior is approximately estimated by sampling in order from above conditional densities.In the study, it is important to choose the number of the components N .Chib's method is used to estimate the marginal likelihood, based on the relation p (X ) = p(X |θ ) p(θ) p(θ |X ) .Thus, the marginal likelihood can be obtained by estimating the 2N runs of the Gibbs sampler, which works well with few components for NMF models.

G. Functional Subnetworks Identification With Bayesian NMF
As mentioned above, the Bayesian NMF method can decompose signals into nonnegative independent components (ICs).Each IC contains an independent spatial activation map (subnetwork) and corresponding time course (TC).In our study, subject-specific spatial activation maps and TCs were obtained for the left-hand MI and right-hand MI using the Bayesian NMF and source reconstruction.For the parameters of the model, we chose the flat prior α i,n =0, β n,i = 0, and the prior of the noise variance was k = θ = 0. We took 20000 points from the posterior and discarded the first 10000 points to allow the sampler to burn in [24].After taking the parameters and data into the model, matrices A and B were obtained.Each row of matrix B indicates the functional activation IC (subnetwork), and each column of matrix A indicates the TCs corresponding to each subnetwork.Then, matrix B was matched to 3D MNI space to visualize activation distributions of subnetworks via source localization by source reconstruction.Before visualizing the subnetwork, the z-score transform was performed on matrix B, where the transformed values smaller than the µ + σ were zeroed.All MI-related ICs were visually inspected to exclude the artificial components.In addition, the selection of the IC number is a key parameter.Inspired by previous studies [24], we set N from 20 to 50 at intervals of 5, and the process was repeated 5 times to find steady components.Finally, we set the number of the IC to 35.We performed the NMF to construct large-scale cortical networks at the α (8-13 Hz), β (13-30 Hz), and all (8-30Hz) frequency bands, respectively.We identified 8, 9, and 10 subnetworks individually.
In order to verify the validity of Bayesian NMF, we also obtained the subnetworks using Gift software (GroupI-CATv4.0b;https://www.nitrc.org/projects/gift).The information maximization algorithm was used to estimate the ICs.Regular was selected as the stability analysis type.In order to be unified with NMF method, we also determined the number of components to be 35.Furthermore, we reconstructed all components by sLORETA and selected the best-fit networks.

H. Functional Networks Connectivity (FNC)
After the IC selection, we identified 8, 9, and 10 MI-related subnetworks at the α(8-13 Hz), β(13-30 Hz), and all (8-30Hz) frequency bands.Each subnetwork contains multiple synchronously activated cortical regions.We extracted the TC of the subnetwork for each subject across states and three frequency bands.Then, the PLV was calculated to construct the FNC matrix for the left-hand MI and right-hand MI.PLV is formulated as follows: PLV indicates the phase synchronization value of signals x and y. ϕ n (t) = ϕ x (t) − ϕ y (t) is the phase difference corresponding to the signals x and y at the time point t.T is the total time points of signals.

I. Networks Properties
To quantitatively characterize the large-scale brain networks, we calculated four network properties based on FNC, including clustering coefficient (CC), local efficiency (LE), global efficiency (GE), and characteristic path length (CPL).The brain connectivity toolbox was used to calculate the network Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
properties (https://sites.google.com/site/bctnet/).In the definitions, m is the number of nodes in the networks, K is the set of all nodes, and (i, j) represents the link between the node i and j(i,j∈ K ).w i j is the connection weight of (i, j) with the range of [0,1].d i j is the shortest weighted path length of (i, j).
Clustering coefficient is formulated as follows: Local efficiency is formulated as follows: Global efficiency is formulated as follows: Characteristic path length is formulated as follows: J. Common Spatial Pattern of Large-Scale Networks (SPN) The network properties may not contain all information to discern the MI tasks.Thus, we used the SPN approach to extract implicit information contained in the spatial topology of the large-scale networks.The key of SPN is to identify spatial filters which minimize the variance of one class and maximize another, extracting the features with significant differences [25].Let A 1 and A 2 be the adjacency matrices for the FNC of left-hand and right-hand MI respectively.In essence, the spatial filters are the projections, equal to maximize the formulation: where M 1 and M 2 are the covariance matrices of A 1 and A 2 for left-hand and right-hand MI individually.Since the object value cannot be affected by the scaling of the projection, the previous formulation can be translated into the following restricted optimization problem: By adding the Lagrange multiplier, the goal function can be expressed as follows: By using the derivative of (15) about λ while satisfying the condition of ∂ L/∂γ = 0, the goal value is estimated as follows: where γ and λ represent the eigenvector and eigenvalue.Then, ( 16) can be solved with multiple m SPN filters as follows: where V includes all the eigenvectors of M −1 2 M 1 .For example, we used three pairs of SPN filters, and then we get six-dimensional SPN features.In theory, the three pairs of SPN filters form a matrix V S P N = [Filter1, Filter2, . . ., Filter6].In this study, for a 10 × 10 matrix P, each SPN filter is a 10-length vector, so the dimension of V S P N is 10 × 6.Additionally, the SPN features are calculated by F S P N = log(var(V T S P N P)), which results in a vector of length 6.

K. Classification of Left-Hand MI and Right-Hand MI
In this study, we first defined the SPN features based on the FNC matrices.Given that the number of SPN filters could influence the classification performance, we set several pairs of SPN filters, such as 2, 4, and 6.Then, we extracted the SPN features for each SPN filter, resulting in 4-, 8-, and 12-dimensional SPN features.In addition to SPN features, we also consider the joint features of SPN and network properties.For classifiers, SVM (SVM with linear and radial basis function (RBF) kernels), LDA, and LightGBM were used.To evaluate the classification performance, leave-onesubject-out cross-validation (LOSOCV) method was utilized.In the study, we have 87 subjects and each subject includes one left-hand MI trial and one right-hand MI trial.In LOSOCV, one subject is used for evaluation and the remaining subjects are used for model training.Finally, 87 times were contributed to this process until all subjects had performed for model evaluation and the results were averaged over all subjects.
Additionally, we used accuracy (ACC), sensitivity (SEN), and specificity (SPE) to evaluate the classification performance.
where n le f t and n right are the correctly predicted number of left-hand MI and right-hand MI tasks, and N le f t and N right are the real number of left-hand MI and right-hand MI tasks.
L. BCI Performance BCI performance represents the evaluation of individual MI ability, which is associated with large-scale network patterns (FNC and network property).Inspired by a recent study conducted by Kang et al. [26], where they calculated individual MI performance using four different deep learning methods.In our study, we chose the results calculated by the DeepConvNet method as the evaluation of individual BCI performance.Due to the missing information of 4 subjects (subjects 34, 43, 78, and 102), only the information of 83 subjects was used for correlation analysis.Detailed information about the BCI performance is provided in Supplementary TABLE SI.

M. Statistical Analysis
We used paired t-tests (two-tailed) to compare the FNCs and network properties differences between left-hand MI and righthand MI tasks.For FNC differences, the significant difference level was set at p < 0.01 (FDR corrected, p < 0.01).For network properties differences, the significant difference level was set at p < 0.01.We used Pearson's correlation analysis to evaluate the relationships between the MI-related indexes and BCI performance at three different frequency bands, where the significant correlation was set at p < 0.05.It is worth noting that the BCI performance represents the individual overall ability to identify left-hand and right-hand MI, so the sum of the network pattern (FNC and network property) of left-hand and right-hand MI was used to calculate the correlation.

A. Identification of Large-Scale Functional Networks
Using the Bayesian NMF method, we identified 8 ICs/ subnetworks at the α  3 shows the spatial maps of ICs at three frequency bands, where the activations were transformed to z-scores and the values below µ + σ were zeroed.
To verify the utility of the Bayesian NMF method, group ICA method was also performed to identify the large-scale networks based on the same MI-EEG data.Compared to the group ICA method, we found that the NMF method can well deal with the negative activation generated by source localization and can effectively characterize the patterns of the networks.This is in line with previous studies [24].The details of these networks using the Bayesian NMF and group ICA methods including activated regions, main Broadman areas, MNI coordinates, peak activation values, and the number of clustering voxels, are summarized in the supplementary material (see Supplementary TABLE SII).

B. Differences in FNCs Between Left-Hand and Right-Hand MI
Fig. 4 shows the mean FNC matrices across all subjects and the FNC differences between left-hand MI and right-hand MI tasks at the α (8-13Hz), β (13-30Hz), and all (8-30Hz) frequency bands.The paired t-tests were used to assess the significant difference (FDR corrected, p < 0.01) in FNCs.For the α (8-13 Hz) frequency band, we found that the FNC was significantly increased in SMN-VN in the left-hand MI compared to the right-hand MI task.For the β (13-30Hz) frequency band, we found that the FNCs were significantly

C. Differences in Network Properties Between Left-Hand MI and Right-Hand MI
Compared to right-hand MI, we found that the CC, GE, and LE were significantly increased ( p < 0.01) and the CPL was significantly decreased ( p < 0.01) in left-hand MI at the β (13-30Hz) and all (8-30Hz) frequency bands.There is no significant difference at the α (8-13Hz) frequency band.Fig. 5 shows the differences in network properties between the left-hand MI and right-hand MI tasks at the α (8-13Hz), β (13-30Hz), and all (8-30Hz) frequency bands.

E. Overall Classification Results and Comparison of Left-Hand MI and Right-Hand MI
We take the SPN features and network properties to construct classification models based on various classifiers, such as LDA, SVM, and LightGBM at three frequency bands.TABLE I shows the classification results of combined network properties and SPN features using 2, 4, and 6 spatial filters in different classifiers at three frequency bands.As shown in TABLE II, we compared the classification accuracy results of our methods with some studies including iterative multi-objective optimization for channel selection with filter band CSP (FBCSP) [27], meta-learning [28], EEG-Net [29], bilinear neural network with 3-D attention [30], and CNN [31].The results showed that our proposed method got a high classification accuracy for cross-subject two-class classification.The highest classification accuracy was obtained by combining two pairs SPN features and network properties at the all (8-30 Hz) frequency band using the LDA classifier (ACC=78.2%,SEN=79.3%,SPE=77%).

IV. DISCUSSION
Exploring large-scale network patterns and finding meaningful network markers are important for both neural mechanisms understanding and practical application of MI.In the current study, we combined source reconstruction and the Bayesian NMF method to identify the large-scale networks of left-hand and right-hand MI tasks at different frequency bands.Based on the identified large-scale networks, we compared the FNC and network properties differences between left-hand and right-hand MI tasks at the α (8-13Hz), β (13-30Hz), and all (8-30Hz) frequency bands.The results showed that the left-hand MI task produced stronger FNC (i.e., increased VN-SMN connectivity) and modulated more network resources (i.e., higher CC and GE), specifically at the β and all (8-30Hz) frequency bands.Moreover, the FNC and network properties could be used as important features to construct classification models to discern left-hand and right-hand MI tasks.These findings highlight the utility of the NMF method in construing robust EEG FNC and provide fresh insights into the neural mechanisms of MI.In terms of methodology, group ICA [13] and Bayesian NMF [32] are both commonly used source separation methods in EEG and fMRI analysis.Group ICA can decompose the aggregated signals into multiple ICs (subnetworks), each IC contains positive and negative activations.However, the negative activation of IC is usually not considered in previous studies, which may ignore a part of important information and violate the physiological source activations [24].Compared with group ICA, the Bayesian NMF method only captures the non-negative spatial distribution of the networks, thus handling the negative activation of components well.Besides, it has unique ability to infer features compared with ICA [33].For these reasons, the Bayesian NMF method has been used in source separation, latent factor extraction, and signal processing [32].For example, Yi et al. [24] applied Bayesian NMF with EEG source imaging to construct a large-scale network and investigated the neural mechanisms of decisionmaking tasks.Thus, we utilized the Bayesian NMF with source reconstruction to construct MI-related large-scale cortical brain networks.Furthermore, we made a comparison with the group ICA method as shown in Supplementary Fig. S1 to verify the validity of the Bayesian NMF method.Our results suggest that the Bayesian NMF method is available to identify MI-related large-scale networks.
In our study, using the Bayesian NMF method, we identified 8, 9, and 10 components at three different frequency bands, respectively.The RPN is absent at the α frequency band and the LPN is absent at the β frequency band, while all components are identified at the all (8-30 Hz) frequency band.Although we identified different ICs, these large-scale networks are highly correlated with MI.For example, the SMN is a motor-related network including premotor cortex (PMC) (see Supplementary TABLE SII), which is associated with MI and motor execution (ME) [34].The hand MI evoked strong activation in the sensorimotor areas [35].Additionally, MI and ME consist of similar neural networks, including sensory and motor-related networks [1].The VN (i.e., IC11 and IC14 at the α frequency band) is a sensory network, involving visual processing, motor-related functions, and high-level cognition including motor planning [36].The CEN specializes in controlling interference and properly allocating cognitive resources, playing an important role in MI [37].The FN areas contribute to motor and cognitive functions.The TN areas were activated during action observation [38].Moreover, the PN areas could affect the motor performance, locating between the sensory and motor cortex [39].In addition, we found that the identified networks were mostly related to motor performance and visual information processing, which were highly correlated with MI.Thus, our results suggest that identification-based large-scale networks are meaningful for studying MI.
Based on these large-scale cortical networks, we further investigate the large-scale FNC differences at three different frequency bands.For the α frequency band, we observed that the FNC between the SMN and VN was significantly increased in the left-hand MI compared to the right-hand MI.The neural underpinnings of MI is related to mu rhythm [40].The mu rhythm in the α frequency range was measured over sensorimotor and occipital areas during MI [41].Ding et al. [42] found that the reduced SMN-VN functional connectivity was highly correlated with the decline of motor functions in patients with white matter lesions.In addition, the SMN-VN connection was also significantly increased at the β and all frequency bands.We also found that the SMN-VN connection was significantly positively correlated with BCI performance at the α and all (8-30Hz) frequency bands.These results implied that the SMN-VN connection played a key role in Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the performance of MI.In the β and all (8-30Hz) frequency bands, we found that the FNCs were significantly increased for left-hand MI compared to right-hand MI, and these increased network interactions were mainly located in the right hemisphere, such as SMN-VN, RPN-VN, RPN-RTN, and RTN-RCEN.These results are consistent with previous MI-related reports [43], [44].For instance, Dentico et al. [43] revealed that the VN-PN (from occipital to parietal) connection was increased in mental imagery compared to visual perception, using dynamic causal modeling and Granger causality analysis.Moreover, the correlation analysis showed that the SMN-VN, LCEN-RPN, and LCEN-RTN were significantly positively correlated with BCI performance.These indicated that the increased FNCs might promote individual MI performance.By comprehensively comparing the different frequency bands, we found that the number of FNC differences was larger at the β and all frequency bands, and the FNC differences mainly concentrated on the right hemisphere and the connections between VN and other networks in the left-hand MI compared to right-hand MI.Previous studies have shown that left-hand MI exhibits right hemispheric laterality, while right-hand MI exhibits left hemispheric laterality [45], [46].Thus, we also suggest that the left-hand MI produces stronger FNC and requires modulation of more right hemispheric network resources than right-hand MI at the large-scale level.
In addition to the FNC, the network property was another important index to investigate the pattern of large-scale networks in MI.Consistent with previous studies [7], we found no significant differences in network properties between left-hand MI and right-hand MI at the α frequency band.To our best knowledge, few studies have explored differences in large-scale network patterns between left-hand MI and righthand MI at the high frequency bands.Previous studies have indicated that the MI-related information at the β frequency band served as a coordinative role during the sensorimotor process [47], while the information at the α frequency band involved in local information processing at higher connection costs [48].In our study, the results showed that CC, GE, and LE were significantly increased and CPL was significantly decreased in the left-hand MI compared to the right-hand MI at the β and all (8-30 Hz) frequency bands.Considering that the FNC difference is also mainly shown in the high frequency bands, we speculate that the coupling of network information at the α frequency band is the basis of MI, whereas the coupling pattern of network information at the high frequency band is an important specificity indicator for distinguishing left-hand and right-hand MI tasks.Moreover, the network properties reflect the transfer efficiency of information globally and locally [49].For example, the LE and GE indicate the efficiency of the information transmission via the networks, and the higher values indicate the stronger ability of information integration [50].Our findings show that the left-hand MI has better network information transfer efficiency than the right-hand MI in the high frequency band.Similar results were also reported that the CC, GE, and LE during left-hand MI was higher than right-hand MI [51].We also found that the GE was significantly positively and CPL was significantly negatively correlated with BCI performance at the all (8-30 Hz) frequency band, which were consistent with Zhang et al. [52].These results indicated that an efficient network facilitates MI performance.Indeed, for right-handed subjects, we have reason to believe that imagining right-handed movements is easier and requires fewer brain resources than imagining left-handed movements.In addition, Li et al. [7] revealed that the left-hand MI needs more brain resources of information processing than the right-hand MI.Similar to that of Zhang et al. [6], who revealed that functional connectivities were significantly increased during left-hand MI compared to right-hand MI.These are consistent with our findings.Thus, we suggest that left-hand MI provided more network efficiency of information transfer and modulated more network resources compared to right-hand MI.
In addition, large-scale brain FNCs and network properties not only help us understand MI but can also be used as discriminative features to build MI-BCI classification models [53].Cross-subject and within-subject MI-BCI classification are the two most concerning issues in current MI-BCI applications.Among them, cross-subject MI classification is challenging and it is difficult to achieve a high classification accuracy.The classification model we proposed based on FNC and network properties achieved a high classification accuracy (see TABLE II and TABLE SIII).We also compare with several state-of-the-art methods that work on the same MI-BCI dataset.The results show that the classification accuracy of our proposed method is 0.2%-14.6%higher than other methods, indicating that the large-scale network feature is also effective for MI-BCI classification.The proposed method has some advantages for classification.In the study, we extracted spatial features by converting EEG signals from sensor domain to source domain which may be possible to reduce the effect caused by individual differences and provide more valuable information.Then, our proposed NMF-based model can deal well with the community detection problem in complex networks and may reduce the loss of useful information caused by negative activation.Moreover, we employed the SPN filters to retain the effective features.The proposed technique can improve discriminant ability, tackling the multi-classification in MI-BCI tasks well.The results indicated that our proposed model based on large-scale network features can effectively recognize MI-BCI tasks.

V. LIMITATIONS
There are some considerations and methodological limitations in this study.First, the ICs are manually selected by prior MI-related information, which may be subjective.Future work needs to find a valid way to select MI-related ICs automatically.Although large-scale FNC analysis can provide fresh sight into the cortical network patterns of MI, such static patterns of FNC cannot capture the dynamic network patterns of interactions among the networks.In future work, the dynamic FNC should be considered to investigate the dynamic patterns of large-scale networks during MI.Moreover, we only explored the neural mechanism in the task state.We should combine resting and task data to comprehensively consider the potential neural mechanism of MI in future work.

VI. CONCLUSION
In this study, we combined the EEG source localization and Bayesian NMF methods to examine the FNC patterns of MI and identify discriminative network biomarkers for MI-BCI control.Source localization method can improve the spatial resolution of EEG signals, and the NMF method can well identify MI-related brain network components, which provides us with good technical support for studying largescale networks.The results show that the FNCs and network efficiency of information transfer were significantly increased in the left-hand MI compared to the right-hand MI, specifically at the β (13-30Hz) and all (8-30Hz) frequency bands.In particular, we also found that the stronger the FNC of VN-SMN, LCEN-RPN, and LCEN-RTN, the better the participants' MI-BCI performance.In the study, our large-scale FNC pattern analysis is conducive to understanding the underlying neural mechanism of MI illiteracy and provides new theoretical support for clinical rehabilitation applications.In addition, the cross-subject classification model based on the FNCs and network properties proposed in our study can find meaningful network connections suitable for recognizing the individual's left-hand and right-hand MI process, which is conducive to the development of personalized training methods for different clinical rehabilitation patients and BCI operators.They can improve the understanding of the network mechanism of MI-based rehabilitation and the development of MI-BCI.
(8-13 Hz) frequency band, 9 ICs at the β (13-30 Hz) frequency band, and 10 ICs at the all (8-30Hz) frequency band.These large-scale networks are the sensorimotor network (SMN), left central executive network (LCEN), right central executive network (RCEN), visual network (VN), primary visual network (PVN), left frontal network (LFN), right temporal network (RTN), left parietal network (LPN), and right parietal network (RPN).The VN and PVN are the same type of networks with similar functions.When defining the name of ICs, we define ICs as the same name with the same functions.In the study, we refer to VN and PVN collectively as VN.At the α (8-13 Hz) frequency band, IC 11 and IC 14 are categorized into VN.At the β (13-30 Hz) frequency band, IC 1, IC 10, and IC 21 are categorized into VN.At the all (8-30Hz) frequency band, IC 18, IC 19, and IC 25 are categorized into VN.The rest networks have one component.Fig.

Fig. 4 .
Fig. 4. The mean FNC across all subjects for left-hand and right-hand MI and FNC differences between the two groups at the α, β, and all (8-30Hz) frequency bands.The star denotes a significant FNC difference in the left-hand compared to the right-hand MI.The red lines represent that the FNCs are significantly increased in the left-hand compared to the right-hand MI.

Fig. 5 .
Fig. 5.The differences in four network properties (i.e., CC, CPL, GE, and LE) between the left-hand MI and right-hand MI at the α frequency band (a), β frequency band (b), and all (8-30 Hz) frequency band (c).In the figures, the asterisk represented the significant difference between the two MI tasks (p < 0.01), and the colored circles represented the subjects.

Fig. 6 .
Fig. 6.Significant correlations between BCI performance and FNCs and network properties at different frequency bands (p < 0.05).The black lines denoted the fitted curve, and the colored circles represented the subjects.

TABLE I CLASSIFICATION
RESULTS BASED ON SPN FEATURES AND NETWORK PROPERTIES