Twenty-three volunteers (7 females and 16 males) aged from 21 to 36 (26 ± 4) participated in the experiments. The volunteers were right-handed according to the Edinburgh Handedness Inventory and had no reported neurological or other disorders. All of them gave written informed consent.

None of the participants had an experience of controlling a BCI. The ability of the participants to perform motor imagery was not assessed prior to the experiments.

The original EEG dataset for the work was obtained using the same procedure as described in earlier works (Frolov et al., 2012, 2017a). Each subject had to sit relaxed, looking at the center of the screen or had to imagine flexion of either his left or right hand according to visual cue presented on the screen. Each task had to be performed for 10 s, with a 2 s preparation period. There were 20 cue presentations for each hand motor image, which were split into blocks containing 2 cues for each hand. In each block, the cues were given in a random order. Each motor imagery task was preceded by the relaxation period.

There were from 10 to 20 experimental days (14 ± 3) for each subject, with 1–3 day gaps possible between the sessions.

The EEG data were recorded using the NVX52 acquisition device (Medical Computer Systems, Russia) with 32 electrodes placed according to the 10–10 system. The data were band-pass filtered with a 5–30 Hz Butterworth filter as well as a 50 Hz notch filter to suppress power line interference.

The EEG data were classified with a Bayesian classifier under the assumption that the signal has a multivariate Gaussian distribution with a zero mean and a different covariance matrix for each task classified:

where C_i is a covariance matrix of the EEG corresponding to the ith task, and x is an EEG sample. According to the Bayes rule, the probability that a given sample will correspond to the ith class is proportional to P (x|i)P(i), where P(i) is the probability of the ith task being cued. Thus, for a signal epoch logarithm of the probability of all its samples to correspond to the ith class is proportional to 1nt(−ΣtxtTCi−1xt)−ln det Ci+ln P(i)+const, where summation is performed over all the epoch samples, n_t is the epoch length in samples, and the constant term is independent of both the signal and the class number. Since 1nt(−ΣtxtTCi−1xt) can be substituted with trace(cov(X)Ci-1), where cov(X) is the epoch covariance estimate under the assumption of a zero mean, the class number is

We used a window of 1 s length sliding with a 100 ms shift to present the feedback during the sessions. The probabilities P (i) were estimated based on the cue durations. The classifier was updated after each block of cues, with no feedback during the first block.

All of the blind and non-blind source separation methods used can be, in general, represented in the form of multivariate linear decomposition with an optional noise term.

where X is the decomposed signal in matrix form, a_i, i = 1, …, n, are columns of A that define the source (component) weights representing the source topographic maps, ξ_i, i = 1, …, n, are row vectors of the source activities, n is the number of EEG channels, and is the optional noise term. Matrix A is called the mixing matrix, and its inverse, W, is called the unmixing matrix. The invertibility of A is based on the assumptions that there are as many components as there are signal channels and that the signal correlation matrix has full rank. These assumptions are in general not required but are common for linear source separation methods (Cichocki and Amari, 2002; Hyvärinen et al., 2004), and all of the decomposition methods used in the paper were derived under these assumptions.

Usually, the signal is whitened before searching for decomposition (3) so that the covariance matrix of X becomes an identity matrix. Also, it is well-known that the component activity variances are not defined, thus requiring the application of some constraints when the decomposition is computed. That is why we shall further consider the decomposed signal to be whitened, i.e.

where V is a matrix such that Vcov(X)V^T = I. If we also suppose that all ξ_i, i = 1, …, n, have unit variance, then the mixing matrix for Z is orthogonal, Z = UΞ, UU^T = I, and the EEG mixing matrix is A equals VU.

We used 16 source separation methods to obtain the decomposition (3).

PCA, or Principle Component Analysis, is a classic example of a technique for obtaining the decomposition (3). It uses Singular Value Decomposition to diagonalize the signal covariance matrix. Usually, the components with too much or too little variance are discarded. If one rather chooses to keep all of the components under the constraints that their variances are equal to 1, this will become signal whitening. In this case, U = I, and the unmixing matrix is not orthonormal.

KURT is an ICA method based on maximizing the difference of the component distributions from normal. The difference is measured by the absolute value of the component excess kurtosis, which serves as a cost function for iterative component search (Delfosse and Loubaton, 1995; Hyvärinen et al., 2004), i.e., the columns of unmixing matrix W are estimated sequentially by maximizing the excess kurtosis for wiTZ under the assumption that W is orthogonal.

CUMUL is an ICA method using the signal non-stationarity as the criterion of statistical independence (Hyvarinen, 2001). The non-stationarity is measured by the fourth-order cumulant:

The decomposition seeks to maximize the sum of the source cumulants. The proof of non-stationarity maximization being a criterion of component independence can be found in Hyvärinen et al. (2004), chapter 18. The time lag τ serves as the method parameter. In our experiment, it was set to 100 ms to extract the components with dominant alpha-band frequencies.

FastICA is an ICA method in which mutual information reduction between the components serves as the statistical independence criterion. The mutual information between the components ξ_i, i = 1, …, n, is defined as

where H (·) denotes entropy. The entropy of linear transformation is H (WZ) = H (Z) + |det(W)|. Under the assumption of Z being whitened and the components having unit variance, W is orthogonal, which implies that det (W) = 1 and H (WZ) does not depend on det (W). As a result, H(wiTZ) can be used as a cost function for sequential search of the unmixing matrix columns w_i.

Since computing the entropy requires knowledge on the signal distribution, different parametric approximations are used. FastICA uses approximation in the form

where g (x) is a non-linear function and v is a random normally distributed variable with unit variance. FastICA can be computed using different non-linearities (Bingham and Hyvärinen, 2000; Hyvärinen et al., 2004). In this work, we used two functions, namely, tanh(x) and xexp(−x²/2). The corresponding methods are denoted as FastICAT and FastICAG.

RunICA, or extended infomax, is an ICA method that is similar to FastICA. The method takes into account both super-and sub-Gaussian distributions by applying two different non-linearities, g₊ (x) = −2 log cosh (x) and g-(x)=log cosh (x)-x2/2, and switching between them. Originally, the infomax algorithm was introduced by Bell and Sejnowski (1995) based on a neural network approach to maximize the entropy of the neural network output. It was later reformulated in terms of parametric non-linearities for source entropy estimation (Lee et al., 1999).

AMICA, which stands for Adaptive Mixture Independent Component Analysis, is an ICA method based on the assumption that all EEG components have generalized Gaussian mixture distributions

AMICA can also segment the signal into several parts corresponding to different models (8) when the number of the models is set to higher than 1. All distribution parameters, unmixing matrix coefficients, and model scores (in the case of several models) are estimated using the expectation-maximization algorithm (Palmer et al., 2006, 2008, 2012). In our work, we compared AMICA with the default settings (2 models, 2 distributions in mixture) to AMICA with two models and Gaussian distributions with means fixed to zero (2 models, 2 distributions in mixture, ρi,jl=2, μi,jl=0, where l indicates the model, l = 1, 2, and i = 1, …, n, j = 1, 2), and to AMICA with a single model and Gaussian distributions with means fixed to zero (1 model, 2 distributions in mixture, ρ_ij = 2, μ_ij = 0, i = 1, …, n, j = 1, 2). The two latter cases will be denoted as AMICA1 and AMICA2, respectively.

PWCICA is an ICA method mapping the signal into a complex domain using the first derivative (rate) of the signal as its imaginary part. The algorithm is described in Ball et al. (2016). It uses the complex variant of FastICA in order to find decomposition (3) in the complex domain (Bingham and Hyvärinen, 2000). After that, the real-valued unmixing matrix is computed using the algorithm provided in Ball et al. (2016).

SOBI, or Second Order Blind Identification, is a blind source separation method that uses joint diagonalization to remove signal cross-covariance for a set of defined time lags (Belouchrani et al., 1997). Particularly, a set of time lags δ_k is specified, and corresponding cross-covariance matrices are estimated:

Next, the unmixing matrix W is found by joint diagonalization of the estimated cross-covariance matrices:

where ||·||_off denotes the Euclidean norm of the off-diagonal part of a matrix. The norm of the unmixing matrix column is constrained to be equal to 1. There exist several algorithms for solving the problem (11). We used one proposed in Ziehe et al. (2004) to obtain the solution.

CSP, or Common Spatial Patterns (Ramoser et al., 2000; Bashashati et al., 2007, 2015), is a non-blind decomposition method, in contrast to the source separation techniques described above. This method utilizes the information on signal segmentation with respect to the experimental tasks in order to find the decomposition (3). The method is used for feature extraction in two-class separation problems. The decomposition is sought that maximizes the component variance ratio for the states classified under the assumption that total component variance is fixed. This is equivalent to finding matrix W such that

where C_i is the i-th class covariance matrix and D₁ is a diagonal matrix. The equations (12) have an explicit solution where W and D₁ result from singular-value decomposition (SVD) of the matrix D-1/2 UTC1UD-1/2 with orthogonal matrix U and diagonal matrix D given by SVD of C₁ + C₂. We used several CSP decompositions in this work, namely, CSP12 comparing relaxation and left hand motor imagery, CSP13 comparing relaxation and right hand motor imagery, CSP23 comparing left with right hand motor imagery, and CSP1X comparing relaxation and both left and right hand motor imagery. Although equations (12) are solved using signal epochs corresponding to the selected tasks only, the unmixing matrix was used to estimate the activity of the components at all other time moments.

MCSP is a generalization of the CSP method for multiple-class problems. There are several methods of generalization (Lee et al., 2005; Bashashati et al., 2015). The most frequent approaches are to either compute all pairwise CSP and combine the corresponding components into a new signal to be classified or to use classifier voting or a tree of some sort when only two classes are compared at each step (e.g., relaxation vs. motor imagery with subsequent discrimination between the hands if the motor imagery was recognized). However, these techniques would not give us any additional information in terms of sources, as we have already considered all pairwise CSP decompositions for our experimental tasks, as well as the CSP1X comparing relaxation and motor imagery. That is why we have decided to use a different approach. Problem (12) was generalized so as to find W such that

The equations are written for the case of three tasks, but the generalization for the case of an arbitrary number of classes is straightforward. The problem in general has no explicit solution and can be solved numerically using any joint diagonalization technique. The joint diagonalization algorithm was the same as was used to perform the SOBI decomposition (Ziehe et al., 2004). It should be noted that unlike the two-class CSP, component variance ratio maximization (or minimization) for different classes is not guaranteed, but if any component has high (small) variance for a certain class, it will have small (high) variance for other classes.

The mutual information reduction (MIR) provided by the decomposition methods is defined as

where the mutual signal information I is given by equation (6). Using (6), we get

Since X = AΞ,

and thus

The entropy of each individual channel and component was estimated and bias-corrected as in Delorme et al. (2012).

The decomposition (3) shows that the contribution of each component to the signal registered by the EEG electrodes is given by the weight vector, which can be treated as potential distribution over the scalp surface. This allows the source of electrical activity represented by the component to be localized. As has been shown for the case of visual memory tasks, the physiologically meaningful components tend to have dipolar distributions of their weights, i.e., the distributions that could be adequately approximated by a single current dipole inside a head model (Delorme et al., 2012). We had neither anatomical magnetic resonance imagery (MRI) scans nor digitized electrode positions for most of our subjects. This is why a standardized head model based on the ICBM MNI atlas (Fonov et al., 2009, 2011) was used for estimating the dipolarity of the component weights. The dipolarity was computed as a residual variance of the best single dipole fit for the corresponding distribution. The fit was considered acceptable if the residual variance did not exceed 10%.

The number of the same components found by different decomposition methods can be used as a measure of similarity between the methods. On the other hand, the more methods reveal the component, the less likely it is that the component is an artifact of a certain mathematical procedure underlying the decomposition algorithm. That is why we focused our attention on the components found by different methods.

Two component similarity measures were used to determine which components were shared between the methods. The first measure was the absolute value of the cosine between normalized component weight vectors given by the corresponding mixing matrix columns. The second one was the absolute value of the Pearson correlation between the component activities. Two thresholds were chosen for the measures. The components were considered the same when both measures exceeded the corresponding thresholds. We have chosen 0.9 as the threshold for the topographic map similarity and 0.8 as the threshold for activity correlation.

When the number of shared components is obtained, the similarity between the decomposition methods can be measured as

where n_s is the number of shared components, and n_i and n_j are the numbers of components obtained by the i-th and j-th methods, respectively. The similarity measure varies from 0, when there are no shared components, to 1, when the methods have provided identical decomposition. Also, only the dipolar components can be considered when the similarity (17) is calculated, allowing non-dipolar components that are unlikely to have physiological meaning to be discarded. The number of methods that found a component minus one will be called the component rank. The component rank equals 0 when it is found by only one method and 15 when it is found by all of the 16 methods considered.

Given a set of the EEG components, their activity can be used for estimating the accuracy of discriminating between the experimental tasks. This estimate can serve as a measure of components' specificity for the performed tasks, allowing the task-specific EEG patterns to be extracted and both artifacts and irrelevant EEG activity to be removed.

We used a greedy algorithm to find the task-specific components as follows. Only the dipolar components were considered. At the first step, classification accuracy was estimated for each set containing one or two components. The best set was then selected. The remaining components were added to the set one by one so that the new set would provide the highest classification accuracy. The components from the set for which the classification accuracy was maximal (over all of the steps) were considered as the task-specific components. The classification was tested by 120 trials of cross-validation where 7 blocks were chosen for the training set and the 3 remaining blocks were used as the testing set.

The accuracy of task classification was measured by Cohen's kappa index, κ (Vieira et al., 2010), which was calculated from the elements of the confusion matrix (g_ij) resulting from the cross-validation. The element g_ij counts how many times the i-th task was recognized by the classifier under the condition that the classified EEG epoch corresponded to the j-th cue. Given g_ij, the κ index is

where g₀ is the sum of all of the elements of the confusion matrix, ∑igii is the sum of all diagonal elements of the confusion matrix, and g_j is the sum of all of the elements of the j-th column of the confusion matrix. The index varies from −1 to 1 (perfect classification) and equals 0 in case of random classification.

The decomposition (3) gives at least as many components as there are electrodes for each of the methods used. Consequently, the number of components computed for all experimental recordings is too high to check the components manually. In order to investigate EEG patterns most typical for our experimental tasks, we used cluster analysis to group the components according to their similarity. The clustering was performed using the Attractor Neural Network with Increasing Activity (ANNIA), which was originally proposed for Boolean factor analysis (Frolov et al., 2007). Adaptation of the technique for cluster analysis is described in Bobrov et al. (2014). When the ANNIA is used for clustering, the strength of synaptic connections between the neurons of the network is determined by the similarity between the corresponding elements rather than by Hebbian learning. The stopping criterion is based on a threshold that specifies the minimal average similarity between the elements of a cluster. In Bobrov et al. (2014), the ANNIA was used to group the components with respect to their topographic map similarity. In this paper, we have also accounted for the component activity similarity using the correlation coefficient between the components' activity power spectral densities estimated for each of the experimental tasks. Spectral analysis was used since the components obtained from different recordings were compared, in contrast with the shared components search described earlier. Both the component similarity measures considered vary from 0 to 1, and their mean was taken as the similarity measure for ANNIA.

The methods were compared according to the criteria described in the previous section. The percentage of dipolar components found by each method is presented in Table 1 (median values and quartiles are presented). The table also contains average values of mutual information reduction values for each of the methods. The p-values resulting from pairwise comparison of the percentage of the dipolar components found are presented in the upper diagonal part of Table 2. The values were obtained using a Wilcoxon test of the values pooled from all of the participants and sessions with subsequent Benjamini-Hochberg correction (initial significance level: 0.05; resulting critical value: 0.0429). The results indicate that the PWCICA and AMICA methods provide more dipolar components than the other methods. The difference is significant for the AMICA methods and insignificant when PWCICA is compared to other BSS methods. The PCA and CSP methods extracted significantly fewer dipolar components than the other methods.

Pairwise Wilcoxon comparison of the MIR values with Benjamini-Hochberg correction (initial significance level: 0.05; resulting critical value: 0.0275) shows that the ICA methods provide significantly higher mutual information reduction than the PCA and CSP methods. The difference between ICA and other methods in terms of MIR was insignificant, as was the difference between the PCA and CSP methods.

The method similarity was estimated by the fraction of shared components according to (17), providing a method similarity matrix. The similarity matrices for the cases when either all or only dipolar components were considered were used for mapping the method onto a 2D plane with a multidimensional scaling technique. The results of the mapping are shown in Figure 1. Apparently, the methods tend to group according to the decomposition criteria, one group containing the ICA methods except PWCICA, and another group containing the CSP methods. SOBI and PWCICA are distant from both groups, while the PCA is as an outlier.

The results for classification accuracy, obtained after searching for the task-specific components, are shown in Figure 2. Mean κ values are shown for each subject and each method. The methods are sorted according to the κ value averaged over all the subjects. The subjects are sorted according to the κ value averaged over all the methods. The lower diagonal part of Table 2 contains p-values resulting from the pairwise comparison of the κ values for different methods. The values were obtained using Wilcoxon test of the values pooled from all of the participants and sessions with Benjamini-Hochberg correction (initial significance level: 0.05; resulting critical value: 0.320). The CSP methods provide patterns that are significantly better classified than those found by the other methods. The PCA patterns are significantly worse-classified than the patterns found by the other methods. The task-specific component search yielded accuracy values significantly exceeding those obtained during on-line BCI operation, which is likely due to discarding irrelevant, artifact, and noisy components.

The clusters obtained using the ANNIA were sorted according to their component occurrence, i.e., the percentage of the experimental session in which the components of the cluster were found. Figure 3 shows topographic maps and activity power spectral densities (PSD) for the components of the first 12 clusters. The maps and PSDs were obtained by averaging over all of the cluster components regardless of the subject or session. Table 3 presents the occurrence, average dipolarity, average rank, and specificity of the components of each cluster. As stated above, the component dipolarity was measured by the residual variance of the topographic map fit with a potential distribution resulting from a single current dipole, and the component rank is the number of methods that found the component minus one. The specificity was calculated as the percentage of cases when the component was included in the best set of components (the set search is described in the Methods section) among the cases when the component was found.

It should be mentioned that AMICA restricted to only Gaussian distributions in mixtures performed as well as AMICA with the generalized super-Gaussian mixtures, suggesting that dipolar components obtained from band-pass filtered EEG signal have activity the distribution of which can be adequately approximated with a mixture of two Gaussians.