Functional disconnection of associative cortical areas predicts performance during BCI training.

Brain-computer interfaces (BCIs) have been largely developed to allow communication, control, and neurofeedback in human beings. Despite their great potential, BCIs perform inconsistently across individuals and the neural processes that enable humans to achieve good control remain poorly understood. To address this question, we performed simultaneous high-density electroencephalographic (EEG) and magnetoencephalographic (MEG) recordings in a motor imagery-based BCI training involving a group of healthy subjects. After reconstructing the signals at the cortical level, we showed that the reinforcement of motor-related activity during the BCI skill acquisition is paralleled by a progressive disconnection of associative areas which were not directly targeted during the experiments. Notably, these network connectivity changes reflected growing automaticity associated with BCI performance and predicted future learning rate. Altogether, our findings provide new insights into the large-scale cortical organizational mechanisms underlying BCI learning, which have implications for the improvement of this technology in a broad range of real-life applications.

Introduction

Voluntarily modulating brain activity is a skill that can be learned by capitalizing on the feedback presented to the user. Such an ability is typically used in neurofeedback control to self-regulate putative neural substrates underlying a specific behavior, as well as in brain-machine interfaces, or brain-computer interfaces (BCIs) (McFarland and Wolpaw, 2018), to directly regulate external devices. Despite the potential impact, from elucidating brain-behavior relationships (Shibata et al., 2011) to identifying new therapeutics for psychiatric (Kim and Birbaumer, 2014) and neurological disorders (Pichiorri et al., 2015; King et al., 2013), both neurofeedback and BCIs face several challenges that affect their usability. This includes inter-subject variability, uncertain long-term effects, and the apparent failure of some individuals to achieve self-regulation (Vidaurre and Blankertz, 2010). To tackle these issues, investigators have searched for better decoders of neural activity (Lotte et al., 2018) as well as for psychological factors (Kleih et al., 2010) and appropriate training regimens (Jeunet et al., 2016) that can influence the user’s performance. On the other hand, neuroplasticity is thought to be crucial for achieving effective control and this has motivated a deeper understanding of the neurophysiological mechanisms of neurofeedback and BCI learning (Sitaram et al., 2016). At small spatial scales, the role of cortico-striatal loops with the associated dopaminergic and glutamatergic synaptic organization has been demonstrated in human and animal studies suggesting the procedural nature of neurofeedback learning (Koralek et al., 2012). At larger spatial scales, evidence supporting the involvement of distributed brain areas related to control, learning, and reward processing has been provided in fMRI-based neurofeedback experiments (Emmert et al., 2016). Recently, a motor imagery (MI) BCI study based on ECoG recordings showed that successful learning was associated with a decreased activity in the dorsal premotor, prefrontal, and posterior parietal cortices (Wander et al., 2013). To date, however, the evolution of large-scale cortical network changes during BCI training has not been tested directly.

On the above-mentioned grounds, we hypothesized that BCI learning would result from a large-scale dynamic brain connectivity reorganization. More specifically, based on previous evidence documenting user’s transition from a deliberate mental strategy to nearly automatic execution (Sitaram et al., 2016), we expected that the reinforcement of activity in the cortical areas targeted by the BCI would be accompanied by a progressive decrease of functional integration in regions associated with the cognitive processes of human learning (Seger and Miller, 2010). Furthermore, we hypothesized that the characteristics of such network changes would contain valuable information for the prediction of the BCI learning rate.

To test these predictions, we simultaneously recorded high-density EEG and MEG signals in a group of naive healthy subjects during a simple MI-based BCI training consisting of 4 sessions over 2 weeks. For both EEG and MEG, we derived cortical activity signals by performing source-reconstruction and we studied the longitudinal task-modulated changes in different frequency bands. Specifically, we evaluated the spatial extension of the activated cortical areas as well as the regional connectivity strength over time. Finally, we tested their relationships with learning as measured by the BCI performance.

Materials and methods

Participants and experiment

Twenty healthy subjects (aged 27.5 ± 4.0 years, 12 men), all righthanded, participated in the study. Subjects were enrolled in a longitudinal EEG-based BCI training (twice a week for two weeks). All subjects were BCI-naive and none presented with medical or psychological disorders. According to the declaration of Helsinki, written informed consent was obtained from subjects after explanation of the study, which was approved by the ethical committee CPP-IDF-VI of Paris. All participants received financial compensation at the end of their participation (around 30 euros per hour).

The BCI task consisted of a standard 1D, two-target box task (Wolpaw et al., 2003) in which the subjects modulated their α [8–12 Hz] and/or β [14–29 Hz] activity to control the vertical position of a cursor moving with constant velocity from the left to the right side of the screen. To hit the up-target, the subjects performed a sustained motor imagery of right-hand grasping (MI condition) and to hit the down-target they remained at rest (Rest condition). Each run consisted of 32 trials with up and down targets, consisting in a grey vertical bar displayed on the right part of the screen, equally and randomly distributed across trials. Each session was divided into two phases (see Supplementary Materials Fig. S1):

The training phase consisted of five consecutive runs without any feedback. For a given trial, the first second corresponded to the inter-stimulus interval (ISI), while the target was presented during the subsequent 5 s. From the data obtained during this phase, contrast maps between conditions were computed for each bin comprised between 5 and 40 Hz to elicit the features, i.e. (channel; frequency) couples of interest that best discriminate the subjects’ mental state over the left motor area and within the mu-beta frequency ranges (see Supplementary Materials Fig. S2). For that purpose, we used the R-square as a metric of such a discrimination between the conditions. To avoid redundant information in the classification, we only selected the couples that show the highest R-square values (typically two channels associated with two frequency bins each). As a rule of thumb, from four to six features were manually chosen (Schalk et al., 2004)).

The testing phase consisted of six runs with a cursor feedback. For a given trial, the first second corresponded to the ISI, while the target was presented throughout the subsequent 5 s, just as in the training phase. The visual feedback, displayed from 3 s to 6 s, consists of a cursor that starts from the left-middle part of the screen and moves with a fixed velocity to the right part of the screen. The subjects were asked to control the vertical position by modulating their brain activity. During our experiments, the online features (i.e. power spectra estimated by an autoregressive model based on the Maximum Entropy Method (Kay, 1988) every 28 ms on time window of 0.5 s with a model order of 28) were classified by using the Linear Discriminant Analysis (LDA) method. The classifier was retrained in each session based on the (channel; frequency bin) couples preselected during the on-going session. Its output enabled the control of the vertical position of the moving cursor through the linear combination of power spectra computed at the pre-selected (channel; frequency bin) couples by applying the moving average method (Ramoser et al., 2009). The present work relies only on the data obtained from the testing phase.

Experiments were conducted with a 74 EEG-channel system, with Ag/AgCl passive sensors (Easycap, Germany) placed according to the standard 10–10 montage. EEG signals were referenced to mastoid signals, with the ground electrode located at the left scapula, and impedances were kept lower than 20kOhms. A system composed by 102 magnetometers and 204 gradiometers collected MEG data (Elekta Neuromag TRIUX MEG system). EEG and MEG signals were simultaneously recorded in a magnetic shielded room with a sampling frequency of 1 kHz and a bandwidth of 0.01–300 Hz. The amplification of the EEG signals was directly performed by the Elekta acquisition system. To digitize the head positions, we used the Polhemus Fastrak digitizer (Polhemus, Colchester, VT). Nasion, left and right pre-auricular points were used as landmark points to provide co-registration with the anatomical MRI. Four Head Position Indicator (HPI) coils were attached to the EEG cap. The subjects were seated in front of a screen at a distance of 90 cm. To ensure the stability of the position of the hands, the subjects rested their arms on a comfortable support, with palms facing upward. We also recorded electromyogram (EMG) signals from the left and right arm of the subjects. Expert bioengineers visually inspected EMG activity to ensure that subjects were not moving their forearms during the recording sessions. We carried out BCI sessions with EEG signals transmitted to the BCI2000 toolbox (Schalk et al., 2004) via the Fieldtrip buffer (Oostenveld et al., 2010). To obtain accurate head models using surface-based alignment (Gross et al., 2013), individual T1 sequences (256 sagittal slices, TR = 2.40 ms, TE = 2.22 ms, 0.80 mm isotropic voxels, 300 × 320 matrix; flip angle = 9°) have been obtained by using a 3T Siemens Magnetom PRISMA after the fourth session. The experiments consisted of a 15 min-resting-state task. Images were preprocessed via the FreeSurfer toolbox (Fischl, 2012) and directly imported (15002 vertices) to the Brainstorm toolbox. In this work, we used the Destrieux atlas (Destrieux et al., 2010). We digitized the location of the EEG electrodes using the FastTrak 3D digitizer (Polhemus, Inc., VT, USA), the landmarks (nasion, left and right preauricular points), and at the scalp. We aligned these locations with the MRI using the Brainstorm toolbox (Tadel et al., 2011).

Several neurocognitive questionnaires were proposed to subjects to assess their specific traits such as the self-esteem (Rosenberg, 1965), and the global motivation (Guay et al., 2003), as well as the ability to perform a motor imagery task (Roberts et al., 2008). Before each session, the subjects’ anxiety was also measured (Spielberger et al., 1983). To reinforce the learning and to improve the subjects’ autonomy (Lotte et al., 2013), the subjects were asked to train at home and alone with a 10-min video, which corresponded to 3 training runs without any provided feedback, each day between two sessions. Between two sessions, a new video was sent to the subjects.

M/EEG preprocessing and source reconstruction

MEG signals were first preprocessed by applying the temporal extension of the Signal Space Separation (tSSS) to remove environmental noise with MaxFilter (Taulu and Simola, 2006). MEG and EEG signals were downsampled to 250 Hz before performing an ICA (Independent Components Analysis) with the Infomax approach using the Fieldtrip toolbox (Bell and Sejnowski, 1995; Oostenveld et al., 2010). The number of computed components corresponds to the number of channels, i.e. 72 for EEG (T9 and T10 were removed). Only the independent components (ICs) that contain ocular or cardiac artifacts were removed. The selection of the components was performed via a visual inspection of the signals (from both time series and topographies). On average, 2 ICs were removed. Data were then segmented into epochs of 7 s corresponding to the trial period. Our quality check was based on the variance and the visual inspection of the signals. For each channel and each trial, we plotted the associated variance values. We kept a ratio below 3 between the noisiest and the cleanest trials. The percentage of removed trials was kept below 10% of the total number of trials (Gross et al., 2013).

After having average referenced the signals, we performed source reconstruction by computing the individual head model with the Boundary Element Method (BEM) (Fuchs et al., 2001; Gramfort et al., 2010). BEM surfaces were obtained from three layers associated with the subject’s MRI (scalp, inner skull, outer skull) that contain 1922 vertices each. Then, we estimated the sources with the weighted Minimum Norm Estimate (wMNE) (Fuchs et al., 1999; Lin et al., 2006; Gramfort et al., 2014) using the Brainstorm toolbox (Tadel et al., 2011). Here, we used the identity matrix as the noise covariance matrix. The minimum norm estimate corresponds in our case to the current density map. We constrained the dipole orientations normal to the cortex. To perform the group analysis, we projected the sources estimated on each subject, and each session, onto the common template anatomy MNI-ICBM152 (Mazziotta et al., 2001) via Shepard’s interpolation. From these estimated signals, we computed the associated power spectra. To identify the anatomical structures associated with the obtained clusters without restricting our work on motor or sensorimotor areas, we used the Destrieux atlas (Destrieux et al., 2010).

Band power and connectivity metrics. Statistical analysis

For each subject, session, and trial in the source space, we computed the power spectra. We used the Welch method with a window length of 1 s and a window overlap ratio of 50% applied during the feedback period that ranged from t = 3 s to t = 6 s. In the case of the group analysis presented in Fig. 2A, we worked within the ICBM152-MNI template. Elsewhere, we used the individual anatomical space. To take into account the subjects specificity when defining frequency ranges (Klimesch, 1999), we used a definition of the α and β bands that relies on the Individual Alpha Frequency (IAF) obtained from a resting-state recording that lasted 3 min (with the subjects eyes open). Similarly to (Pichiorri et al., 2015), the IAF corresponds to the first peak of the power spectrum between 6 and 12 Hz. The α1 ranges from IAF - 2 Hz to IAF, α2 from IAF to IAF +2 Hz, β1 from IAF +2 Hz to IAF +11 Hz and β2 from IAF +11 Hz to IAF +20 Hz. Eventually, we averaged the power spectra values across the frequency bins within each predefined band.

Fig. 2.

Cortical activity changes during BCI training. (A) Task-related activity maps obtained with EEG-source reconstructed power spectra in the α2 and β1 frequency band. The colors code the statistical difference obtained by contrasting motor-imagery and rest conditions through cluster-based permutation t-tests performed at the group level. For illustrative purposes, we show the obtained p-values multiplied by the sign of the t-values. (B) Normalized power spectra across the sessions for the motor-imagery (blue) and rest (red) condition from EEG signals. The significant clusters of activity in each individual were obtained with respect to the inter-stimulus intervals (ISI). The group results for the α2 frequency band are illustrated by the boxplots on the left side of the panel, while results from the β1 frequency band are shown on the right. For illustrative purposes, we plotted the log transformed values. Similar results were obtained with MEG signals.

To perform the analysis presented in Fig. 2B, we computed statistical differences among activations recorded in the MI and the rest conditions at the group level or at the subject level via a paired t-test. Since we expected a desynchronization between the two conditions, we applied a one-tailed t-test. Statistics were corrected for multiple comparisons using the cluster approach (Oostenveld et al., 2010; Tadel et al., 2011). We fixed the statistical threshold to 0.05, a minimum number of neighbors of 2 and a number of randomization of 500, which lead to stable results (Pernet et al., 2015). Clustering was performed on the basis of spatial adjacency. Cluster-level statistics were obtained by using the sum of the t-values within every cluster. To obtain the relative power Δ; we computed the relative difference, in terms of power spectra, between the two conditions, as follows: , where P and P correspond, respectively to the averaged power calculated across the cluster from MI and Rest trials. The cluster size C was obtained by estimating the number of cortical vertices that belong to the cluster that presented the best discrimination between the conditions. To perform the study for each condition separately (Fig. 2), we normalized the power spectra with respect to the inter-stimulus interval (ISI) with the Hilbert transform, similar to the approach reported in (Wander et al., 2013).

The connectivity analysis (Fig. 3) was based on the cross-spectral estimation computed with the Welch method. To reduce dimensionality, we extracted the first principal component obtained from the power spectra calculated across the dipoles within each ROI. Then, we computed the imaginary coherence between each pair of ROIs based on the definition proposed in (Sekihara et al., 2011). From the resulting connectivity matrix, we next computed the relative node strength Δ similarly to what we did for the relative power. The strength of the i-th node was here calculated by summing the values of the i-th row of the connectivity matrix.

Fig. 3.

Cortical connectivity changes during BCI training. Task-related connectivity networks obtained with MEG-source reconstructed signals in the α2 and β1 band are represented on a circular graph. The nodes correspond to different regions of interest (ROIs) and the links code the statistical values resulting from a paired t-test performed between the motor-imagery and rest conditions performed at the group level. Only significant links (p < 0.005) are illustrated for the sake of simplicity. The color of each node, corresponds to a specific macro-area as provided by the Destrieux atlas. Similar results were obtained with EEG signals.

To evaluate the session effect on BCI scores, Δ, C, and Δ one-way repeated non-parametric ANOVAs were applied with the session number as the intra-subject factor. The computation of the p-values was based on the bootstrapping approach with 200 repetitions (Delorme and Makeig, 2004). In the specific case of Δ, the ANOVA was performed separately for each ROI.

To estimate the correlations between BCI scores and, respectively, Δ, C, and Δ, we performed repeated-measures correlations (Bakdash and Marusich, 2017) which control for non-independence of observations obtained within each subject without averaging or aggregating data. In the specific case of Δ, the correlation analysis was performed separately for each ROI. Results presented in sections 3.3 and 3.4 referred to a statistical threshold of 0.05 corrected for multiple comparisons by adopting a false discovery rate (FDR) criterion (Benjamini and Yekutieli, 2001), which is a method extensively used in biological studies (McAuley et al., 2009; Sanders et al., 2012; Matthews and Farewell, 2015) (the associated p values are noted p).

Results

Behavioral performance and BCI controlling features

At the beginning of each experimental session, we identified the controlling EEG features among the electrodes over the contra-lateral motor area and within the standard α and β frequency ranges during a calibration phase (SI Fig. S2). We found that the ability to control the BCI significantly increased across sessions (days) but not within sessions (hours) (SI Fig. S3). The session effect was also present when we averaged the BCI accuracy scores across the runs of each session (one-way ANOVA, F3;57 = 13:9; p = 6:56:10−7). Despite the expected high inter-subject variability (> 8:95%), 16 subjects out of 20 learned to control the BCI by the end of the training, with accuracy scores above the chance level of 57% (Müller-Putz et al., 2008) (Fig. 1A and SI Table S1).

Fig. 1.

(A) Evolution of BCI performance over sessions. Individual performance is measured by considering the average BCI accuracy score (i.e. percentages of correctly hit targets) of the 96 trials in each session. In the violin plots, the black line corresponds to the group-averaged BCI score and the outer shape represents its distribution. The horizontal dashed grey line shows the chance level (57%), which is here considered as learning threshold. (B) Representation of the selected EEG controlling features across all subjects. On the left, we show occurrences obtained across subjects and sessions in terms of pre-selected channels; on the right, we show occurrences in terms of frequency bins selected over the sessions.

We next investigated the characteristics of the EEG controlling features. From a spatial perspective, the electrodes above the primary motor area of the right hand (C3 and CP3) tended to better discriminate the MI and Rest mental states (Fig. 1B). The most discriminant frequencies occurred between high- α and low-β ranges (Fig. 1B). These results are in line with previous studies (Neuper and Pfurtscheller, 2001). Notably, we observed a progressive focus over CP3 and low-β ranges throughout the sessions.

Among the demographical and psychological items that we measured before the experiment, only the kinesthetic imagery score (Roberts et al., 2008) was moderately correlated with the ability of subjects to control the BCI (Spearman, r = 0.45, p = 0.045, not significant after FDR correction). However, these items were not correlated with the evolution of BCI accuracy over time (SI Table S2), and it is possible that a more complete psychological assessment could determine BCI performance (Jeunet et al., 2015).

Reinforcement of sensorimotor cortical activity during BCI training

We evaluated the spatiotemporal cortical changes associated with the BCI training by performing a statistical analysis of the MEG and EEG signals at the source-space level. Separately for each neuroimaging modality, we computed the associated task-related brain activity by statistically comparing the power spectra of the MI versus the Rest condition in each session, across subjects. In both α and β frequency ranges, we found a progressive involvement of EEG sources in the cortical hemisphere contralateral to the movement (Fig. 2). The involved regions exhibited a significant power decrease (p < 0.025), a phenomenon known as event-related desynchronization (ERD) which reflects sensorimotor brain activity (Lopes da Silva, 2013).

In session 3, ERDs were particularly significant in the α2 and β1 frequency bands, and mainly spanned the primary sensorimotor cortex (pre- and postcentral gyri, central sulcus, inferior and superior parts of the precentral sulcus) and secondary higher-order premotor and somatosensory areas (Fig. 2 and SI Figs. S5 and S8). At the end of training, ERDs were more localized in the contralateral paracentral lobule, precentral gyrus, and superior parietal lobule (Fig. 2 and SI Figs. S5 and S9), which are typically involved in hand motor tasks (Yousry et al., 1997) as well as in motor imagery (Solodkin et al., 2004; Lotze and Halsband, 2006) and motor learning (McDougle et al., 2016). No other comparable significant differences were observed in the other frequency bands (SI Fig. S5).

To quantify these changes at the individual level, we calculated in each subject the size C and the relative power Δ of the most significant ERD cluster. These quantities exhibited a significant session effect only in the α and β frequency ranges (p < 0.03, SI Table S3). Notably, BCI training was accompanied by a reinforcement of task-related activity (C and Δ), including areas that are mainly within the sensorimotor territory. These longitudinal changes were explained by the significant decrease of relative power in the MI condition (F3;57 = 4:82; p = 0:003; F3;57 = 3:09; p = 0:024 respectively in the α2 and β1 bands), while the Rest condition did not vary across sessions (α2) or varied with lower extent compared to MI (β1) (Fig. 2). These findings confirm that during the training the subjects actually focused on how to perform MI of the hand.

Similar results were obtained when we considered MEG source-reconstructed signals (SI Figs. S4, S6 and S7). Given the higher spatial resolution of MEG, we will focus on the results obtained with this modality while providing detailed EEG analysis in the supplementary materials.

Progressive functional disconnection of associative cortical areas

To evaluate the cortical changes at the network level, we considered functional connectivity (FC) patterns that have been previously shown to be sensitive to BCI-related tasks (Mottaz et al., 2018) as well as to learning processes (Bassett and Khambhati, 2017). For this purpose, we calculated the imaginary coherence between the source reconstructed signals of each pair of regions of interest (ROIs) corresponding to the Destrieux atlas. Imaginary coherence is a spectral measure of coherence weakly affected by volume conduction and spatial leakage (Nolte et al., 2004; Sekihara et al., 2011).

By statistically comparing the MEG-based FC values between MI and Rest conditions across subjects, we found a progressive decrease of task-related connectivity in both α and β frequency ranges along sessions (SI Figs. S11 and S12). In α frequency ranges, the strongest decreases involved fronto-occipital (α1, α2) and parieto-occipital (α2) interactions. In β frequency ranges, significant decrements involved parieto-occipital (β1) but also fronto-central and bilateral temporal interactions (β1, β2).

For each subject, we quantified the regional connectivity changes by computing the relative node strength Δ in the α and β frequency ranges. Significant across-session declines were spatially distributed involving bilaterally primary visual areas and associative regions (p < 0.025, SI Fig. S13). Specifically, in the α2 band, the Δ values of the ROIs typically associated with visuo-spatial attentional tasks (Silver et al., 2007) (e.g. cuneus) decreased significantly with the training (SI Table S5). In the β1 band, we observed a significant reduction for the orbital part of the inferior frontal gyrus, which is involved in mental rotation (Milivojevic et al., 2009) and working memory (Wilson et al., 2014; Christophel et al., 2017) (SI Table S5).

Collectively, the results indicate that BCI training is associated with a progressive reduction of integration among cortical systems that are specialized for different functions such as motor imagery and learning, visual attention, and working memory. Similar results from the analysis of EEG-based FC networks are reported in SI Figs. S10, S12, S13, SI Table S4.

Predictive markers of BCI performance and learning rate

To better understand how the observed cortical changes in α and β frequency ranges were associated with performance, we next performed a repeated-measures correlation analysis which takes into account the longitudinal nature of our data (Bakdash and Marusich, 2017). In general, we found strong correlations in the α2 and β1 frequency bands where both activity and connectivity metrics significantly correlated with the BCI performance within the same session (SI Tables S6 and SI Fig. S15). Higher BCI scores were associated with a higher number of task-related cortical sources C (Fig. 4A) and with a stronger decrease of relative power Δ (Fig. 4B). From a network perspective, better performance was associated with larger reduction of relative node strength Δ in associative areas and, to a minor extent, in primary visual regions (Fig. 4C and SI Fig. S15).

Fig. 4.

Correlation between activity/connectivity changes and BCI performance. The first row shows the results obtained in the α2 band. The second row shows the results obtained in the β1 band. (A) Scatter plots with cluster sizes C and the BCI scores of all the subjects. Colors identify the values obtained for the same individual across sessions. (B) Correlation values between the relative power Δ and the BCI scores of all the subjects. Same color conventions as in (A). (C) Correlation values between the relative node strengths Δ and the BCI scores in the same session. All correlations values (r) are calculated through a repeated-measures correlation coefficient, with a statistical threshold (p < 0.025). For a detailed account of these results, see SI Table S6 and Fig. S15. Similar results were obtained with EEG signals.

Specifically, we observed significant correlations in regions known to be involved in cognitive aspects of human learning, such as decision making and memory consolidation (middle-anterior part of the cingulate gyrus) (Kolling et al., 2016), change detection and shifts in behavior (posterior-ventral part of the cingulate gyrus) (Pearson et al., 2011), as well as motion detection and tracking (lingual gyrus) (Waberski et al., 2008; Kamitani and Tong, 2006). Furthermore, we observed that areas known to be involved in both mental rotation and working memory (e.g. orbital part of the inferior frontal gyrus) (Milivojevic et al., 2009; Wilson et al., 2014; Christophel et al., 2017), were also correlated with BCI scores. No other comparable significant differences were observed in the other frequency bands (SI Table S6).

In terms of future prediction, we found that only the relative node strength (Δ) was significantly correlated with the learning rate, defined as the relative difference of BCI accuracy between consecutive sessions (p < 0.025, SI Figs. S15 and SI Table S6). Strongest correlations were found in α2 and β1, where higher values of Δ were associated with the a larger learning amount in the following session.

In these bands, the most predictive cortical areas were the anterior part of the cingulate gyrus and the orbital part of the inferior frontal gyrus (Fig. 5), both known to be involved in human learning (Euston et al., 2012). Significant predictions were also reported for the fronto-marginal gyrus in the β1 band and for the superior parietal lobule in the α2 band (SI Fig. S15), which are typically associated with learning and motor imagery tasks (Stephan et al., 1995; Johnson et al., 2002; Solodkin et al., 2004).

Fig. 5.

Prediction of BCI learning rate from regional connectivity strengths. The first row shows the results obtained in the α2 band. The second row shows the results obtained in the β1 band. (A) Colors show the correlation values for the ROIs with a significant effect (p < 0.025). (B) Scatter plots show the values of relative node strengths Δ and the learning rates of all the subjects for the most significant ROIs (p < 0.002). Colors identify the values obtained for the same individual across sessions. The r values correspond to the repeated measures correlation coefficients. For a detailed account of these results, see SI Fig. S15. Similar results were obtained with EEG signals.

We obtained similar results from the analysis of EEG source-reconstructed activity and connectivity SI Tables S6 and SI Fig. S14. Altogether, these findings demonstrate that the observed dynamic cortical changes at the network level were intrinsically associated with successful BCI learning.

Discussion

Neuroplasticity and motor learning

Identifying the large-scale neural mechanisms underlying plasticity is fundamental to understand human learning (McDougle et al., 2016; Hétu et al., 2013; Hardwick et al., 2018; Dayan and Cohen, 2011). The ability to voluntarily modulate neural activity to control a BCI appears to be a learned skill. Investigators have repeatedly documented that task performance typically increases over the course of practice (Ganguly and Carmena, 2009; Moritz et al., 2008), while BCI users often report transitioning from a deliberate cognitive strategy (e.g. MI) to a nearly automatic goal-directed approach focused directly on effector control. This evidence is indicative of a learning process taking place in the brain that is consistent with procedural motor learning. Efforts in understanding the neural dynamics underlying BCI skill acquisition have been made by using neuroimaging techniques in primates (Ganguly and Carmena, 2009; Carmena et al., 2003) as well as in humans (Pichiorri et al., 2011; Wander et al., 2013). These works suggest that even if the use of a BCI only requires the modulation of activity in a motor-related brain area, a dynamic and distributed network of remote cortical areas is involved in the early acquisition of BCI proficiency. However, how such network reorganizes over time and which are the connectivity mechanisms subserving BCI learning is still largely unknown.

Here, we showed that a MI-based BCI learning is associated with a progressive decrease of functional integration of associative cortical regions and with the reinforcement of sensorimotor activity targeted by the experiment. Together with sensorimotor areas, associative areas play a crucial role in motor sequence learning as well as in abstract task learning (McDougle et al., 2016; Hétu et al., 2013; Hardwick et al., 2018; Dayan and Cohen, 2011). In particular, we showed that MI-based BCI learning is accompanied by the disconnection of specific areas related to working memory and decision-making (Fig. 3, SI Figs. S12–S13). Both these cognitive processes are known to be involved in the supervisory attentional system (van Zomeren and Brouwer, 1994; Wolpert and Landy, 2012), which is an important prerequisite for successful motor learning (Wulf, 2007; Lohse et al., 2014; Dayan et al., 2000; Gottlieb, 2012) and also present during MI tasks (Guillot and Collet, 2010). Altogether, we speculate that the observed progressive disconnection of associative areas would mainly reflect the attentional effort fading related to the optimization of the MI strategy to control the BCI and to the gradual skill acquisition.

Brain network predictors of BCI learning

Forecasting behavior from brain functioning is one of the main challenges in human neuroscience. In BCI contexts, the identification of neural features explaining BCI performance will allow to better design adaptive BCIs. Investigators have recognized the need for adaptive BCI architectures that accommodate the dynamic nature of the neural features used as inputs (Vidaurre et al., 2011). Initial efforts have taken into account psychological (e.g. anxiety) or demographical items to explain BCI performance (Benaroch et al., 2019). Other studies have focused on the use of power spectra within theta, alpha and gamma bands as potential predictors of BCI scores (Ahn et al., 2013a, 2013b; Bamdadian et al., 2014; Jeunet et al., 2015). In the specific case of the gamma band, we actually observed a reduced desynchronization effect, more circumscribed within motor areas (SI Figs. S4 and S5). However, a significant correlation with BCI scores has been observed only with the cluster size (and not with the relative power) within the gamma-band (r = 0.47 and p = 10−4 for EEG; r = 0.32 and p = 0.012 for MEG), without being possible to predict the future amount of learning. The effects measured within the gamma band were less significant and harder to interpret than those associated to alpha/beta sub-bands. We also notice that in (Jeunet et al., 2015) the authors failed to replicate results obtained in (Grosse--Wentrup and Schölkopf, 2012), highlighting the difficulty to compare and to reproduce previous studies. Thus, the role of gamma band in long-term training appears still controversial.

Only recently, FC-based metrics have been shown to correlate with the user’s performance suggesting potential strategies for improving MI-based BCI accuracy (Sugata et al., 2014; Pichiorri et al., 2015; De Vico Fallani et al., 2013). However, these findings only referred to the same experimental session and did not inform on the prediction in the future sessions. In this field, there is a critical need for biologically informed computational approaches to identify the neural mechanisms of BCI learning that predict future performance, thereby enabling the generalization of these results across subject cohorts, and the optimization of BCI architectures for individual users (Perdikis et al., 2014; De Vico Fallani and Bassett, 2019). Here, we showed for the first time that the regional connectivity strength of specific associative cortical areas is not only able to explain the BCI performance in the same session (Fig. 4C) but can also predict the learning rate in the subsequent session (Fig. 5). Notably, higher values of relative node strength Δ were associated with larger learning rates, indicating that the potential to improve performance is higher when the functional disconnection of the associative areas has not yet started. These findings could be used to inform future decisions on how to train individuals depending on the current properties of the functional brain network organization. For example, new adaptive BCIs could be designed to target the associative ROIs with highest node strength so to promote the disconnection process associated with learning.

BCI “illiteracy” and performance assessment

BCIs are increasingly used for control and communication as well as for the treatment of neurological diseases (Daly and Wolpaw, 2008). However, a non-negligible portion of users (between 15% and 30%) exhibit an inability to interact accurately with a BCI (Allison et al., 2010). This is a well-known phenomenon that is informally referred as to “BCI illiteracy”. Critically, it affects the usability of BCIs in the user’s daily life but the reasons (Blankertz et al., 2010; Jeunet et al., 2015) and even the definition for such inability (Thompson, 2018) is still under debate. It has been suggested that typical accuracy metrics, based on behavioral performance, could be affected by decoder recalibration (Perdikis et al., 2016), re-parameterizations of the BCI, and the adoption of different mental strategies (Kober et al., 2013; Perdikis et al., 2014).

However, our results indicate that the behavioral performance (i.e. number of hit targets) was significantly correlated with the modulation of sensorimotor power spectra (Fig. 4). This separation between mental states has been previously adopted as a potentially more appropriate indicator of performance in MI-based BCIs (Wander et al., 2013; Perdikis et al., 2018). Thus, at least in our study, we believe that the behavioral performance can be used to measure subject’s performance.

It is important to state that our findings do not represent in any way a general and definitive answer to the underlying causes of BCI illiteracy. Instead, the reported spatiotemporal cortical and functional connectivity changes supporting BCI skill acquisition, may give complementary insights into the understanding of the neural processes behind the BCI illiteracy. Eventually, these results could open the way to the definition of more robust perfomance metrics, integrating both behavior and functional brain changes associated with individual learning (Perdikis et al., 2018).

Caveats and limitations

The temporal window of two weeks considered in our experiment prevents us from observing behavioral and neural changes over longer timescales and therefore, might not be sufficient to observe the full learning process (Yin et al., 2009). Here, BCI skill acquisition was paralleled by a progressive focused activity over the sensorimotor areas together with a loss of large-scale connectivity, which altogether indicate the initiation of an automaticity process typical of procedural motor learning (Seger and Miller, 2010). Future studies are necessary to assess whether and how the observed cortical patterns will evolve with longer BCI training.

While the BCI accuracy was highly variable across individuals, the group-averaged performance at the end of the training is in line with the state-of-the-art (Vidaurre and Blankertz, 2010). It is important to mention that the main goal of the present work was not to maximize the performance but to study the neural mechanisms underlying BCI learning. In this respect, all our experimental subjects were BCI-naive and exhibited on average an increase of performance reflecting a successful BCI skill acquisition.

While EEG and MEG results presented the same trends, some differences can be appreciated as well. MEG-ERD tended to be more localized in parietal and occipital areas while central and frontal areas were more involved in EEG-ERD, especially within the α band (see SI Figs. S4 and S5). Also, the number of ROIs significantly correlated with BCI scores (or with the learning rate) was larger in MEG than EEG (see Figs. S14 and S15). This could be explained by the fact that EEG signals, attenuated by tissues that present a difference in terms of conductivity (Puce and Hämäläinen, 2017), are sensitive to both tangential and radial components of dipolar sources whereas MEG presents a better spatial resolution and its signals are mainly sensitive to the tangential components (Ahlfors et al., 2010). Thus, EEG and MEG do not necessarily capture the same information. In addition, the subjects were instructed to lean their head back within the MEG helmet. As a result, MEG sensors closer to the parietal/occipital areas could better capture the posterior activity than those located above the frontal areas.

From this work alone, we are unable to determine whether or not learning is the only possible modulator of the observed cortical changes. While no correlation has been found with behavioral factors (i.e. anxiety), complementary experiments could be designed to test whether the observed cortical changes are also modulated by fatigue or exogenous stimulants to increase motor excitability (i.e., transcranial direct current stimulation, tDCS (Nitsche et al., 2003)).

Conclusion

Consistent with our hypothesis, we have identified specific cortical network changes that characterize dynamic brain reorganization during BCI training. We found that the progressive functional disconnection of associative areas is crucial for the BCI skill acquisition process. These network signatures varied over individuals and, more importantly, were significant predictors of the BCI learning rate. Taken together, our results offer new insights into the crucial role of brain network reconfiguration in the prediction of human learning.