Decoding Motor Imagery through Common Spatial Pattern Filters at the EEG Source Space

Brain-Computer Interface (BCI) is a rapidly developing technology that aims to support individuals suffering from various disabilities and, ultimately, improve everyday quality of life. Sensorimotor rhythm-based BCIs have demonstrated remarkable results in controlling virtual or physical external devices but they still face a number of challenges and limitations. Main challenges include multiple degrees-of-freedom control, accuracy, and robustness. In this work, we develop a multiclass BCI decoding algorithm that uses electroencephalography (EEG) source imaging, a technique that maps scalp potentials to cortical activations, to compensate for low spatial resolution of EEG. Spatial features were extracted using Common Spatial Pattern (CSP) filters in the cortical source space from a number of selected Regions of Interest (ROIs). Classification was performed through an ensemble model, based on individual ROI classification models. The evaluation was performed on the BCI Competition IV dataset 2a, which features 4 motor imagery classes from 9 participants. Our results revealed a mean accuracy increase of 5.6% with respect to the conventional application method of CSP on sensors. Neuroanatomical constraints and prior neurophysiological knowledge play an important role in developing source space-based BCI algorithms. Feature selection and classifier characteristics of our implementation will be explored to raise performance to current state-of-the-art.


Introduction
Brain-Computer Interface (BCI) is emerging as a promising rehabilitation technology, that aims to establish a connection between brain activity and external devices. Recent advances in invasive BCIs have demonstrated the feasibility of performing complex motor tasks using brain signals by people with disability such as severe spinal cord injury and quadriplegia [1]. As invasive BCIs use intracranial electrodes to measure electrical activity of the cerebral cortex, either implanted or directly lying on the cortical surface such as electrocorticography (ECoG), their usage is limited due to ethical, medical, and physiological issues [2]. ese limitations are not present with noninvasive BCIs, and the most widely used noninvasive modality, electroencephalography (EEG), uses electrodes over the scalp to measure inferred cerebral cortical activity.
A variety of brain signal types and features have been used to decode user intent in noninvasive EEG-based BCIs, such as visual evoked potential (VEP), P300 response, slow cortical potentials (SCP), and sensorimotor rhythm (SMR), to name but a few [3]. SMR or mu (μ) rhythm, typically measured at the alpha band of [8][9][10][11][12][13] Hz over the scalp area overlying the sensorimotor cortex, can be modulated during motor execution or motor imagery (MI) tasks, and the BCIs decoding this type of signal are referred to as SMR-BCIs. Motor imagery displays similar patterns of brain activation and communication to motor execution [4,5] while research and development in the domain of SMR-BCIs has brought some remarkable applications ranging from the accurate control of a cursor in 2-D space [6], control of a quad-copter in 3D space [7], control a robotic arm for reach and grasp tasks [8], and control of a wheelchair [9] proving the potential of this technology.
Nonetheless, noninvasive BCIs also feature a number of limitations with regards to reliability, speed, and accuracy and have many challenges to overcome to meet both research and casual everyday use needs. Key features for the success of SMR-BCIs involve the classification accuracy, performance robustness, and asynchronous and intuitive control that requires the decoding of multiple motor imagery tasks. Control of an external complex device with multiple degrees of freedom, such as a robotic arm or an artificial limb, can be better achieved by utilizing motor imagery classes that are related to the intended end effector movement [10,11], making control more intuitive and thus requiring less time for training.
Moreover, intrinsic drawbacks of EEG include low signal-to-noise ratio (SNR), low spatial resolution, and imprecise and indirect measuring of brain activity mainly attributed to the volume conduction effect. is effect describes the spread of the brain's electrical field while it is transmitted from the source space through the cerebrospinal fluid, skull, and scalp to reach the scalp surface where the electrodes lay, known as the sensor space [12]. To reduce the volume conduction effect and study the brain activity on the cortex, source imaging techniques are commonly used that map the scalp potentials measured by EEG sensors to cortical activations on the cortical mantle [13,14]. Low SNR led the researches to search for spatial filters that extract the EEG components that reflect user intention. In this context, Common Spatial Pattern (CSP) method was proposed to extract spatial features of event related de/synchronization during motor imagery [15]. CSP filters are spatial filters designed to maximize the power difference on their outputs given different EEG classes [15]. CSP filters are considered as an effective way to discriminate classes and are one of the most popular feature extraction methods in the BCI field [16], which also have multiple extensions [17][18][19][20]. Remarkable classification results have already been reported by studies that implemented the CSP algorithm or its variants [21,22].
In the current study, we describe the development of a BCI algorithm, aiming to decode multiple (4) MI tasks. In order to overcome the issues associated with low spatial resolution, we use source imaging and extract features in the cortical source space from selected Regions of Interest (ROIs), using Common Spatial Pattern filters. Finally, the classification is performed with an ensemble classification model that synergistically uses the classification models of selected ROIs, in order to increase classification accuracy.

BCI Competition Dataset.
e BCI Competition IV 2a dataset was used to develop and test the BCI decoding algorithm. e dataset contains recordings from 9 healthy subjects that perform 4 motor imagery tasks, right arm, left arm, feet, and tongue [23]. e data of a subject consist of 2 sessions, one intended for training and the other for evaluation. Each session is comprised of 72 trials for each MI task, 288 trials in total, recorded with 22 EEG channels and 3 monopolar electrooculogram (EOG) channels (with left mastoid serving as reference). In our study, only data from the training session are used.
At the beginning of each trial (t � 0 s), a white cross on black background appeared, and after 2 s, an arrow oriented right, left up, or down informed the subject to perform the corresponding MI task ( Figure 1). e arrow appeared for 1.25 s, and the subject was asked to keep on performing the MI task until the white cross disappeared (t � 6 s).

Signal Preprocessing.
Signal analysis was performed solely on the EEG electrodes, and the EOG channels were excluded. Average reference was used, and the data were band-pass filtered at 7-15 Hz using a zero-phase FIR filter in order to capture the event related desynchronization and synchronization (ERD/ERS) activity [24]. Subsequently, data were down-sampled at 100 Hz and epoched for 500 msec after the visual cue with epoch duration of 3000 msec. Data were visually inspected for bad channels but none was excluded. All preprocessing was performed using a custom Fieldtrip script [25].

Inverse Problem Solution.
EEG source imaging was deployed to mitigate low spatial resolution and low SNR caused by volume conduction. EEG source imaging maps sensor activity to brain neural current distribution at fixed positions over the cortex. e source activity is defined in terms of current dipoles, at a grid of vertices on the MNI cortical surface template, that model electrical activity of neuronal groups firing synchronously [26,27]. e estimation of the sources from the EEG recordings constitutes the solution for the inverse problem, while the forward problem is described by the following equation (assuming zero noise) [27]: where M is the N c × T matrix of the EEG data, G is the leadfield matrix (also referred as gain matrix) that maps the source data to sensor data (N c × N d ), and D is the dipole current density (N d × T). N d is the number of current dipoles, N c is the number of EEG channels, and T is the number of measurements. Solving the forward problem   [28] was used as the default subject anatomy to compute a three-compartment (scalp, skull, and cortex) head model with symmetric boundary element method (BEM) using OpenMEEG [29]. Default Brainstorm Colin 27 cortex was down-sampled using iso2mesh [30] to 5023 vertices, and the relative conductivity values of Scalp/Skull/Brain was assumed to be 1 : 1/15 : 1, with σ brain � σ scalp � 0.33 S/m and σ skull � 0.0042 S/m [31]. All 5023 dipoles are assumed constrained to the cortical surface with an orientation perpendicular to the surface, based on the assumption that EEG primary signal sources are local groups of pyramidal neurons firing synchronously, located on the cortex and arranged perpendicular to its surface [26,31].
Given the lead-field matrix, the inverse EEG problem consists of finding the dipole current density D in (1). is is a highly underdetermined problem since the number of dipoles (sources) is at the order of thousands and the number of EEG channels is at most at the order of hundreds, which in practice means that different current distributions (brain activity) can lead to exact EEG sensor values. Among different methods for solving the inverse problem, here it was solved with the weighted minimum norm estimate (wMNE) method using the Brainstorm toolbox [32][33][34]. Sensor noise covariance matrix, required for the computation of the solution, was calculated on the resting state period at the start of the session.

Regions of Interest. Cortical Regions of Interest (ROIs)
were defined on the sensorimotor cortex to reduce the dimension of the source data derived from the inverse problem solution, having anatomical constraints and aiming at extracting valuable information related to MI tasks [11,35]. 24 ROIs were defined based on neuroanatomical landmarks and Broadman areas and are depicted in Figure 2. Defined ROIs include bilaterally presupplementary motor area (pSMA), supplementary motor area (SMA), cingulate motor area (CMA), dorsal premotor cortex (PMd) and ventral premotor cortex (Pmv), primary foot motor area (M1F), primary hand motor area (M1H) and primary lip motor area (M1L), primary foot somatosensory area (S1F), primary hand somatosensory area (S1H), secondary somatosensory area (S2), and somatosensory association cortex (SAC). During ROI analysis, source times-series that lay only on the defined ROIs on the mantle are analyzed, excluding from analysis all the other sources.

Feature Extraction.
Feature extraction was performed at the source level, on ROIs data in particular. Common Spatial Pattern (CSP) filters are one of the most used feature extraction methods in BCI domain [16]. Assuming data of two classes, for example, the motor imagery of right and left, CSP algorithm calculates spatial filters that maximize the ratio of variance of data stemming from the two classes. Consequently, the extracted signals are optimally discriminating two different EEG classes while they are revealing spatial patterns of different classes [15,17]. e spatially filtered signal S of an EEG trial is given by where M is a N c × T matrix representing the EEG measurement of data for the given trial and W is L × N c matrix referred as CSP projection, whose rows are the spatial filters designed to output signals whose ratio of variances are maximally discriminating input data of two different classes.
Original CSP algorithm has been developed for two class problems, though there exist multiclass extensions [17,37]. Since the classification problem of this work is multiclass, a multiclass extension of CSP was deployed using the Onevs-Rest scheme, with L � 8 filters, the last and the first eigenvectors of each class were selected [15]. CSP filters were calculated during training phase, on the mean covariance matrices of the data conditioned to the four classes. In this work, CSP filters were applied to the source data, and they were calculated on every ROI current dipole timeseries. Assuming D q the current dipole times-series of ROI q , that resulted from the solution of the inverse problem (1), and W q the CSP filter computed on the data of ROI q , the output of the ROI-CSP filters is e feature vector of ROI q , v ROI q ∈ R L is extracted from CSP filters output, and each of its components, v p , p � 1, . . . , L, is given as where Z p is the p row of the matrix Z, that is, the output signal of the pth CSP filter output. Repeating this procedure for each selected ROI results to Q feature vectors of L � 8 elements, where Q is the number of selected ROIs.

2.6.
Classification. An ensemble classification model was used for the prediction of the MI task [38], illustrated in Figure 3. It is supported that an ensemble method using multiple independent classification models can increase the classification performance [39,40]. An independent classification model was built for each of the selected ROIs, and the final classification outcome was selected by an inference (fusion) mechanism. K-nearest neighbors (kNN), Naive Bayes, Decision Tree, and Linear Discriminant Analysis (LDA) classifier were tested with LDA having superior performance as it is demonstrated in the Results section. e ROI classification model was based on LDA, and the inference mechanism was the majority vote of the selected ROI classification models. e selected ROIs were the Q most accurate ROIs according to a selection procedure that is presented in the next section.

ROI Selection.
e defined ROIs extend all over the motor cortex, while the cortical activity related to the performed motor imagery tasks is derived only from a subset of the defined ROIs. ROIs were selected based on their classification model accuracy. In order to select the most accurate ROIs, 10-fold cross-validation using the LDA classifier was performed on ROI level, and this was repeated 10 times to ensure more robust results (in every run, CSP filters are calculated on different data  Figure 3: Outline of the implemented decoding algorithm. EEG sensor time series are transformed to current dipole time series. Data from the Regions of Interest (ROIs) are spatially filtered by ROI-CSP filters, to extract features to be classified by independent ROI classification models. Predicted class is the most voted class of the ROI classification models. On the classification model scheme, the predicted class is the outcome of an inference mechanism (majority vote). e inference mechanism takes as input the predicted class from the individual ROI classification models.
accuracy. Parametric analysis was run for different subjects, and the number of selected ROIs was set to Q � 8. e performance of the classification scheme on the source space was further compared to the performance on the sensor space using the same setting (10-fold crossvalidation of the LDA classification repeated 10 times). For the sensor space, the CSP filters are computed on the preprocessed EEG data. Moreover, to better assess the developed method, performance in terms of Cohen's kappa statistic, a useful metric for multiclass prediction problems, was compared to the winner of BCI Competition IV of dataset 2a [23,33]. e winner of the competition deploys CSP on multiple frequency bands (FBCSP) as feature extraction method, Mutual Information-based Best Individual Feature (MIBIF) algorithm for feature selection, and Naïve Bayesian Parzen Window (NBPW) classifier [21,41].

Classification Accuracy.
Four different classifiers were tested to select the classifier to make the predictions. LDA, kNN, Naive Bayesian, and Decision Tree were tested by performing 10-fold cross-validation, 10 times on all subjects. LDA had superior performance with the highest prediction accuracy among all subjects, with mean accuracy 54.1%. Naive Bayesian was second with 46.9%, followed by Decision Tree and kNN with 45.5% and 44.5%, respectively (Figure 4). e source method of classification achieved consistently higher accuracy rates across all subjects (43.7% to 74.5%), when compared to the sensor method (37.7% to 73.4%), as illustrated in Figure 5 and displayed in Table 1 below. Comparison of the developed method's performance to the winner of BCI Competition IV of dataset 2a in terms of Cohen's kappa statistic [42,43] (multiclass prediction) is presented in Table 2.
Classification sensitivity and specificity, also referred to as true positive and true negative rate respectively, between the source and sensor method are demonstrated in Tables 3-6. Among the subjects, the source method has mean 11.1% higher true positive rate for the left arm, 5.2% higher for right arm, and 3.3% and 1.9% better rate for foot and tongue imagery, respectively. e mean differences of sensor to source true negative rate metric are low for all classes, −1.2%, 1.9%, 2.9%, and 3.6% for left, right arm, foot, and tongue imagery, respectively.

Selected ROIs.
e ROI selection procedure was performed for all the subjects, exhibiting interesting intersubject properties. As illustrated in Figure 6, the symmetrical left and right S1H, M1L, M1H, and CMA ROIs were the, Q � 8, most selected among all the subjects, with the left and right S1H, and left M1L, M1H, and CMA being selected for all 9 subjects. For the subjects A02T, A06T, A09T the pMd_R, SAC_L, and S2_R were selected instead of M1L_R, M1H_R, and CMA_R, respectively, with still 7 out of 8 selected ROIs being on the most frequent ones. Most frequently selected ROIs are illustrated on Figure 7 on the cortical mantle model.

Discussion.
Noninvasive BCI systems emerge as a promising and safe solution for rehabilitation purposes in contrast with invasive BCIs that are associated to health risks and ethical issues [44], but their commercial use is still hindered by low performance and instability. Despite a number of already demonstrated SMR-BCI applications [8,9,45], noninvasive BCIs still suffer from low SNR. In our work, we investigate the use of source imaging and Computational Intelligence and Neuroscience subsequent application of CSP on the source space to compensate for the head volume conduction by mapping scalp potentials to cortical activations [46]. ere are several studies supporting that BCI algorithms based on source space features are superior to the sensor ones [11,47], an observation that is confirmed in our study, comparing the classification results on the sensor and source space. Our BCI algorithm uses, in particular, sources belonging to select Regions of Interest (ROIs) on motor cortex for feature extraction and an ensemble classification model to take advantage of ROI data. Despite the fact that our algorithm did not reach the accuracy levels of the wining method of the BCI Competition, during the ROI selection procedure, common ROIs emerged among all subjects. e emerged ROIs are anatomically and neurophysiologically related to the MI tasks, linking the method results with neurological data. Given that the motor tasks of the competition involved motor imagery of both arms, tongue, and feet, consistent selection of primary hand motor areas and primary lip motor area (cortical representations of hands and face on the primary motor cortex) seems very promising. Cingulate motor areas are also considered very important nodes of the sensorimotor network, having been demonstrated to drive the sensorimotor process [48,49]. It is our conviction that selected ROIs, as produced by the developed algorithm, validate our results since there is a clear neuroanatomical and neurophysiological link between these ROIs and the Motor Imagery tasks performed in the dataset.
Our method appeared to improve mean accuracy by 5.6% and by 0.07 Cohen's kappa value among all subjects, with respect to sensor method. When our algorithm is compared      with the winner of the BCI Competition (FBCSP), the mean accuracy is considerably lower by 0. 19 Cohen's kappa value. e performance of our algorithm based on kappa value is considered moderate while that of the winning implementation is considered substantial [40]. We believe this difference is attributed to feature selection and classifier used by the winner. FBCSP generates CSP features in different frequency bands resulting to multiple features, while feature selection procedure is a vital component to detect the most discriminable features [41]. On the other hand, our sourcebased algorithm seems to increase the classification accuracy of subjects with the worst performance, namely the A04, A05, & A06, as it can be illustrated in Figure 5, outperforming sensor algorithm by a mean accuracy rate of 7.5%. Nevertheless, the trend identified cannot lead to safe conclusions yet, since we cannot infer statistical significance of the results, requiring further investigation on data with larger population of subjects. Moreover, ensemble classification was used, in an effort to increase classification accuracy, by synergistically deploying the independent ROI classification models. Majority vote of ROI classification models was used as final classification outcome, although a weighted vote taking into account the ROI-MI task relation could be considered in the future.

Limitations and Future
Steps. In this study, a generic template three-compartment BEM head model was utilized to solve the forward problem. Forward problem solution induces an important error in the source estimation, as has been explored extensively in previous studies [50,51]. Main forward problem error inducing factors are (a) the use of the MNI template MRI data rather than the subjects' individual neuroanatomy and (b) the absence of cerebrospinal fluid (CSF) in the forward modeling. Template anatomy was used for all subjects, missing important geometrical information for every subject, producing a lead-field matrix that transforms the EEG sensor data into a template cortical manifold different from the real one. CSF compartment has big influence on both signal topography and magnitude, resulting in strong signal attenuation for superficial sources on gyral crowns [52]. is effect is termed to the high increase of conductivity between the sensors and sources. In a future effort to address this problems, a 4-or 5-compartment head model including CSF and skull anisotropy will be used, modeled with finite elements [53].
ere are two main shortcomings in the use of CSP method that were not dealt in this work. e first is that the CSP filters are prone to noise and overfitting, and the second is that the CSP performance is highly dependent on the input signal frequency band the individual subject BCI performance is dependable on individual frequency band used [20,21].
ere are many variants of conventional CSP algorithm designed to overcome the limitations, with popular variants being the RCSP that tackle noise and overfitting with regularization and the FBCSP that better capture the individual subject multiple frequency band filters feature selection, respectively, while a newly introduced method combines aforementioned methods [20,21,39].
Future work will focus on better CSP filters extraction and use feature selection and more sophisticated ensemble models, in an effort to increase the performance of the algorithm. Since the anatomy used for the forward model is common among all subjects and the selected ROIs are common among all subjects, we would like to check the potential of the algorithm in transfer learning between subjects. ere is a study supporting that transfer learning between different subjects by means of source space can achieve higher average single-trial classification accuracy than with a conventional method [54]. Beyond the BCIC IV 2a dataset that is a common ground for the evaluation of methods decoding multiple MI, we aim to evaluate the improved method on dataset we compiled for the

Conclusions
Source estimation and application of CSP filters at the source space constitute a promising solution to increasing classification accuracy of noninvasive BCIs. Our method has demonstrated capability in decoding multiple motor imagery tasks with better accuracy than the equivalent sensor method. While our implementation still is not superior to the state of the art of BCI algorithms, feature selection and classifier characteristics can improve performance. Neuroanatomical constraints and prior neurophysiological knowledge has been shown to play an important role in developing source spacebased BCI algorithms. Our results indicate that the selected ROIs are common among all subjects, which worth further investigation probably in the context of transfer learning between different subjects.

Data Availability
e BCI Competition IV dataset is available at http://www. bbci.de/competition/iv/. e data from the hereby described analysis can be made available from the authors upon request.

Ethical Approval
is study describes a novel analysis of a publicly available dataset. It does not describe new experiments on human subjects.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this article.