EEG sensorimotor rhythms’ variation and functional connectivity measures during motor imagery: linear relations and classification approaches

Dataset and number of subjects

Data of left and right hands MI from a 64 channel EEG were obtained from the Physionet’s open database (Goldberger et al., 2000; Schalk et al., 2004), from which acquisitions from 10 subjects were analyzed in this study.

Data were acquired at a 160 Hz sampling rate. Experimental protocol consisted of randomly alternating blocks of task (right or left hand MI) or rest periods. Each block lasted approximately 4 s (for further details of the acquisition protocol, please refer to Schalk et al., 2004). Three runs were recorded for each subject, and each run contained between seven and eight blocks of each MI tasks, depending on the recording. To increase the quantity of available samples for our classifier, we analyzed them all as an ensemble, with each sample consisting of the whole 4 s block. For each subject, a total of 24 samples were available for training and testing the classifier.

Data pre-processing

Data were filtered in the two frequency bands of interest most prominent in MI studies: µ(7–13 Hz) and β (13–30 Hz), using the standard FIR (finite impulse response) filter (Rabiner & Gold, 1975) of EEGLab (Delorme & Makeig, 2004), a MATLAB suite. This filtering operation would not have been necessary if we were to focus our analyses only on the signals PSD. However, since the signal time domain was considered to build the graphs, its component frequencies could influence calculations of the corresponding adjacency matrices. Thus, the filtering process enabled us to build specific graphs for each frequency band.

To attenuate common artifacts arising at all channels at the same time, data were spatially filtered using a CAR (common average removal) filter (Ludwig et al., 2009). Our reason for choosing this specific filter was three-folded: (1) CAR has been related to decreasing volume conduction effects, which is of particular interest for scalp EEG functional connectivity studies (Thakor & Tong, 2009; Yuksel & Olmez, 2016; Brunner et al., 2016). (2) Although more complex spatial filters have also been suggested (Thakor & Tong, 2009; Yuksel & Olmez, 2016; Brunner et al., 2016), a previous work has shown that the CAR filter provided similar MI classification results to more sophisticated techniques under specific feature selection situations (Uribe et al., 2016). Since in this work we were also interested in this issue, using the CAR filterin as a first approach was suggested by this finding. (3) Studies of the electrode referencing effect on the functional connectivity outcome have shown that using the average of the electrodes’ signal is one of the least distortive approacheswhen compared to other strategies, such as using the Cz electrode or the mastoids (Chella et al., 2016) (note that applying the CAR filter re-references the electrodes’ signals to their average).

Graphs construction

Graphs connectivity matrices were calculated using the motifs’ synchronization method (Rosario et al., 2015). The general principle is to divide the original time-series into smaller ensembles of data points (here we used three data points, as suggested by (Rosario et al., 2015)) and to label these new patterns according to Fig. 1A. This translates the original recorded EEG series into a new, labeled one—the motifs series (Fig. 1B). See, for example, the sample signal of Fig. 1B (blue curve with red dots corresponding to the sampled data points). The time-series is analyzed within a three-points sliding window, and the segment in this window is labeled accordingly. This process is performed for the time series of all electrodes. Due to the highly noisy nature of the EEG signal, we did not consider the possibility of two adjacent signal samples having the same value (see Fig. 1A).

The next step is to define a similarity coefficient between electrodes i and j (s_ij) (i.e., between their time series). We calculated this by counting the number of coincidences (N) amongst their motif series’ elements and normalizing the result by the series’ length (L); that is: (1)${\displaystyle s_{ij}=\frac{N}{L}.}$ Under this scheme, s_ij is normalized within the range [0,1]. To avoid losing important information by arbitrarily thresholding our graph’s edges for binarizing them, we chose to work with weighted graphs.

In summary, s_ij represents the i, j-th element of the functional connectivity matrix (FCM) related to the electrophysiological brain activity measured by EEG.

The motifs’ synchronization method should partially compensate for the volume conduction problem, since only the direction of change of the signal is considered, not contemplating effects that directly influence the signal’s magnitude. Moreover, since the construction of the functional networks was performed after the CAR filter was applied, the volume conduction issue should already have been at least partially accounted for.

Features extraction, selection and data classification

Features of two distinct methods were tested, obtained either from the signal PSD or the connectivity approach, represented by the FCM.

The PSDs were estimated via Welch’s transform (Welch, 1967), and values from multiple frequencies were gathered for each electrode for further use by the classifier algorithm. More specifically, all frequencies of each band’s range were contemplated, collecting PSD values at unity steps. Therefore, each electrode yielded 7 and 18 frequencies (features) in the mu and beta band, respectively. In the case of the FCM, each one of its elements ( s_ij) could be selected as a feature.

To study the effect of FS methods on the outcome classification results, we analyzed three approaches: (1) we designed a wrapper that should find the optimal set of electrodes, which would maximize the classifier’s outcome response. The wrapper’s search ended when no improvement occurred after four consecutive iterations. The other two FS strategies consisted of using a Pearson’s (2) or a Fisher’s (3) filter (Duch, 2006), combined with the wrapper described in (1). The basic idea of these filters is to rank the best attributes to be selected for classification according to a specific criterion. Pearson’s filter estimates how strongly a feature and its labeled-class are correlated by using Pearson’s linear correlation (Campbell & Swinscow, 1996). On the other hand, Fisher’s filter ranks attributes based on the criterion of the Fisher’s discriminant (McLachlan, 2004); that is, on the ratio between the difference of means and variances across data classes for that specific feature.

A least squares-based linear discriminant analysis (LDA) was used for MI data classification. This method was chosen due to its simplicity and robustness. Moreover, it has been commonly employed in BCI research, with results as good as those of other more complex classifiers (see, e.g., Carvalho et al., 2015).

All classification tests were performed using the leave-one-out scheme.

Correlations between PSD and FCM

To investigate if any linear relationship between ERDs due to the hands MI tasks and variations on the FCM (when compared to the rest condition) could be observed, we performed analyses using Pearson’s correlation (Campbell & Swinscow, 1996).

Relative variations (either as ERDs or ERSs) of the signal’s PSD at a given frequency f on electrode i at a specific task block t (ΔPSD(f)_i,t) were estimated as: (2)${\displaystyle \Delta PSD{\left(f\right)}_{i,t}=\frac{PSD{\left(f\right)}_{i,r}-PSD{\left(f\right)}_{i,t}}{PSD{\left(f\right)}_{i,r}}.}$ In Eq. (2), PSD(f)_i,r represents the PSD value for electrode i and frequency f during an average of rest blocks. Note that a positive value for ΔPSD(f)_i,t indicates an ERD occured, whereas a negative one expresses an ERS.

Since each one of the 64 graph nodes contained information about other 63 connections, analyzing the s_ij values directly would require an elevated computational cost. Therefore, our correlation analyses were estimated indirectly by studying how the weighted degree of each node varied according to the ΔPSD magnitude, rather than the s_ij themselves. In doing so, we hypothesized that alterations in the functional networks (represented by s_ij) should yield variations in the weighted degree for the graph nodes.

For a weighted graph, such as the ones built here, the weighted degree for a node i (W_i) can be calculated as (Zhang et al., 2012): (3)${\displaystyle W_i=\sum _j{s}_{ij}.}$ Variations relative to rest periods were then estimated for the weighted degree similarly to Eq. (2): (4)${\displaystyle \Delta W_{i,t}=\frac{W_{i,r}-W_{i,t}}{W_{i,r}}.}$ In Eq. (4), indices t and r refer to a task block and the average of rest blocks, respectively.

Therefore, our correlation analyses considered the values of ΔPSD and ΔW for all trials of each electrode.