The overarching goal of this study was to identify potential brain plasticity effects of BCI. Motivated by the above‐mentioned previous results indicating rapid signs of brain plasticity in MRI measures, we hypothesized that such effects might be seen after only 1 h of purely mental BCI. Particularly for the engaging motor imagery task, where a performance increase during MI‐BCI is usually observed, we expected to find MRI signatures of plasticity. We further hypothesized that the two types of BCI (MI‐ and ERP‐BCI) differentially impact on brain structure and function by modulating the respective targeted areas, namely sensorimotor areas associated with motor imagery (MI‐BCI) and areas associated with visual target detection (ERP‐BCI).

There is an abundance of applications for BCI, including assistive devices, computer gaming (Blankertz et al . 2016 ), and also clinical aids and interventions (Millán et al . 2010 ). While neurofeedback BCI tasks are associated with transient changes in brain activity during a task, it is not known whether they also alter neural structure and function beyond performance of the task. Such BCI‐induced neural plasticity could be especially relevant for potential therapeutic applications of BCI. Using non‐invasive MR‐imaging, it has become possible to identify short‐term and long‐term brain plasticity, usually associated with different types of learning in human subjects. For example, task‐based fMRI has been applied successfully to show brain plasticity effects for motor and other types of tasks (Karni et al . 1995 ; Poldrack et al . 1998 ; Fletcher et al . 1999 ; Pleger et al . 2003 ; Limanowski et al . 2017 ). Since task‐based fMRI tends to be confounded by changes in task performance, several studies have successfully used task‐free ‘resting‐state functional connectivity fMRI’ to identify neural connectivity changes (Taubert et al . 2011 ; Zhang et al . 2014 ; Ge et al . 2015 ; Amad et al . 2017 ; Mawase et al . 2017 ). Noteably, ‘structural’ MRI measures (T1‐weighted and diffusion‐MRI) have also been shown to identify signs of brain plasticity. After the seminal work of Draganski and colleagues, who showed changes in grey matter density (GMD) after 3 months of juggling training, several other studies have also demonstrated changes in GMD after different types of learning within several days or weeks (Draganski et al . 2004 ; Taubert et al . 2010 ; Zatorre et al . 2012 ). Recently, several MRI studies successfully used ‘structural MRI measures’ to detect cerebral modifications after only a few hours of training or stimulation (Sagi et al . 2012 ; Hofstetter et al . 2013 ; Taubert et al . 2016 ; Schmidt‐Wilcke et al . 2018 ), consistent with animal studies that showed signs of brain plasticity already within a few minutes of training (Xu et al . 2009 ; Moczulska et al . 2013 ; Kuhlman et al . 2014 ).

Common types of BCI are based on sensorimotor rhythms modulated by a motor‐imagery task (MI‐BCI), or event‐related‐potentials (ERP) elicited by visual paradigms. The modulation of the visual‐evoked‐potential (ERP‐BCI) is usually induced by observing and counting flashing letters on a screen. Subjects rapidly achieve good performance that typically does not increase over time (Guger et al . 2009 ), but rather is modulated by factors such as fatigue, habituation, or fluctuating attention (Polich, 2007 ). In contrast, the modulation of sensorimotor rhythms by imagining limb movements (MI‐BCI) (Pfurtscheller & Lopes da Silva, 1999 ) is a more active and complex task where, if the system is adequately tuned to the subject, increased performance is often observed (Blankertz et al . 2010; Lorenz et al . 2014 ; Sannelli et al . 2016 ). While there is a large fraction (about 30%) of MI‐BCI users that are unable to gain control (Guger et al . 2003 ), recent progress has been made to drastically reduce this percentage by using co‐adaptive approaches, i.e. not only do the MI‐BCI users modulate their sensorimotor rhythms to gain BCI control, but the algorithm also adapts to the user to increase the system's accuracy (Vidaurre et al . 2011a , Acqualagna et al . 2016 ).

Then the EEG data from the first and last runs in the selected channel were again filtered between 8 and 12 Hz and the sgn‐r 2 values of both runs in that channel were computed. Finally, this resulted in one difference value of sgn‐r 2 from the first to the last run for each subject, and these values were correlated (on the second level) with the change in the MRI data from pre‐ to post‐BCI training using Spearman's correlation as implemented in Matlab (corr.m). The correlation between the sgn‐r 2 difference and the MRI data difference was calculated for each voxel, giving their respective correlation coefficients and respective P (or t ) values. To correct the statistical T‐maps for multiple comparison to the alpha‐level, P corr < 0.05, AlphaSim correction as implemented in AFNI (Version 16.2.12; Cox, 1996 ) was applied within the bilateral somatosensory‐motor mask to identify clusters of significant correlation between change in sgn‐r 2 values and change in MRI parameters. As MRI parameters, T1‐weighted signal intensity (derived from the structural scans) and eigenvector centrality were analysed separately. For the correlation with ECM, we used the difference of the two scans that were acquired directly before and after the BCI training (RestPost – RestPre 3 ).

Based on the results of the previous spectral analysis, the EEG data was preprocessed as follows in order to obtain values to correlate with MRI data: the EEG data was filtered in the mu band between 8 and 12 Hz using a Butterworth filter of order 5. Then, the logarithmic band power of the signal in all channels was computed using the data from the last run. From those values, the sgn‐r 2 coefficient was calculated and the channel on the left hemisphere exhibiting the most positive or negative sgn‐r 2 value was selected. The sgn‐r 2 is a measure of discriminability of univariate features, in this case, the ability of the mu band power to discriminate between the two MI‐BCI tasks, and thus it is related to the classification performance.

In order to find group effects on the power of brain signals after MI‐BCI training that would be candidates to perform a correlation analysis with MRI data, a spectral analysis employing subject‐specific information was performed. First, the inspection of all CSP patterns that were used to deliver feedback in run 4 revealed that all subjects had filters corresponding to the right‐hand (left hemisphere). Then, the trial‐wise power on all channels of the left hemisphere was extracted in the subject‐specific frequency band and time interval of run 4, and the signed square of the point‐biserial correlation coefficient ( sgn‐r 2 ) of each channel was computed. The sgn‐r 2 expresses to what extent the tasks (right hand versus feet motor imagery) are discriminable, and thus we could select the channel with highest absolute r 2 . Then, a spectral analysis over the complete frequency range of interest (5–30 Hz, thus including mu and beta bands) on this selected channel was performed for the first and last run and the result averaged over all subjects. The discriminability over all frequency bins was computed based again on sgn‐r 2 values.

Users conducted a 1 h BCI session of an ERP (even‐related potential)‐based spelling task derived from the Hex‐O‐Speller as in (Treder et al . 2011 ). They performed one calibration run followed by three feedback runs. The first two feedback runs consisted of copy spelling the phrases ‘let your brain talk’ and ‘the summer comes’, respectively. The last run was a free spelling recording in which the users spelled a sentence of their choice, which was communicated to the experimenter for the later computation of the performance. The design of the Hex‐O‐Speller consists of six discs arranged in a circle that allows users to choose one out of 30 different symbols (letters of the English alphabet, punctuation marks ‘.’ and ‘,’, a space symbol and a backspace symbol). It is a two‐stage speller, where a symbol group (e.g., ‘ABCDE’) is selected at the first stage (Treder et al . 2011 ). The symbols of the chosen group are then expanded on the other discs so that, at the second stage, the target symbol can be selected. For selection of a disc, the six discs are subsequently intensified, that is, increased in size and at the same time, they flash in a specific colour. The time between subsequent onsets of intensifications was 200 ms and each disc was intensified for 100 ms. Users were instructed to silently count the number of times their selected disc was intensified. The users’ focus on the target led to a noticeable increase in the P300 amplitude that could be extracted from the EEG (Treder et al . 2011 ). For online and offline ERP classification, EEG data were downsampled to 100 Hz and divided into overlapping epochs ranging from −200 to 800 ms relative to the onset of the stimulus (intensification). Baseline correction was done on the 200 ms pre‐stimulus interval and no artifact rejection was performed. LDA with shrinkage of covariance matrix was applied as classifier (Treder et al . 2011 ). As spatial features, all 62 electrodes except for 12 frontal and temporal electrodes (i.e. Fp1,2; AF3,4,7,8; FT7,8; T7,8; TP7,8) were included. The signed square of the point‐biserial correlation coefficient ( sgn‐r 2 ) was used as a measure for the discriminability of targets vs . non‐targets. For temporal feature selection, between four and seven temporal windows were automatically extracted using a heuristic search for peaks in the sgn‐r 2 between targets and non‐targets. This yielded a feature vector with 50 spatial features times four to seven temporal features, hence a total of between 200 and 350 spatio‐temporal features.

For the second level, the features were again logarithmic band power, but they were spatially optimized using CSP filters and changing Laplacian derivations, i.e. it included spatial flexibility of the features. The classifier was LDA with shrinkage of the covariance matrix and was recomputed after each trial (Vidaurre et al . 2011a ). Level three consisted of logarithmic band power of fixed CSP filters classified with LDA adapted in an unsupervised manner (Vidaurre et al . 2011a ), to be able to assess the BCI control of the user.

The first level consisted of a simple BCI system where the features were the logarithmic power values of fixed mu and beta bands computed from small Laplacian derivations over C3, Cz, and C4. The classifier was a linear discriminant and was updated after every trial (Vidaurre et al . 2011a ). This simple system allowed fast and supervised adaptation over the two frequency bands of interest for BCI (mu and beta bands). However, due to the lack of data from the BCI subjects in this study, the initial classifier of run 1 was subject‐unspecific, i.e. did not contain specific information from the subjects performing the BCI experiment. Instead it was computed from previous datasets of 50 users whose performance was above 70% of accuracy in the class pair of interest, i.e. right hand and feet motor imagery (Sannelli et al. 2019 ).

The MI‐BCI experiment consisted of a 1 h BCI session with four runs of 100 trials each. All users imagined right hand versus feet motor imagery. The subjects were asked to perform two different movements with the right hand (feet) separately: the rotation of the wrist (ankles) or the flexion/extension of the wrist (ankles). Then, they were asked to focus on the sensations that the movements caused on their hand, wrist and arm (feet, ankles and legs). After that, they were asked to choose the movement that caused more intense feelings on their body. Then, they were asked to perform the selected movement several times more. In order to quickly train the imagination of movements, after each performance the participants were asked to stay still and feel the movement as if they were executing it, but without actually moving. Then, the subjects were instructed to fixate on a cross in the middle of the screen. Trials started with an arrow behind the fixation cross, indicating the motor imagery task, that is, an arrow pointing right or down for right hand or feet motor imagery. For online classification, every 40 ms features were calculated as the log‐variance of the band‐pass and spatially filtered data (Vidaurre et al . 2011a ) of the last 750 ms and classified by linear discriminant analysis (LDA) classifier. Feedback for the subjects started at t = 2s until t = 6s, that is, it lasted 4 s. Depending on the result of the LDA, the fixation cross moved to the right or down, updated every 40 ms (this is approximately the speed at which the human eye perceives continuous movement) (Vidaurre et al . 2011a ). Performance in a trial was ‘correct’, if the position of the fixation cross at the end of a trial ( t = 6s) corresponded to the direction of the arrow. Between runs, users were asked to explain the motor imagery strategy that they had followed in the previous run to identify participants not complying with the task. Two subjects reported trying to use eye movements instead of imagining motor movements of their limbs and were eliminated from further analysis. The experimental paradigm consisted of three different methodological adaptation levels as in (Vidaurre et al . 2011a , b ) with run 1 belonging to the first level, runs 2 and 3 to the second, and run 4 to the third level. In those experiments, we observed a significant increase in BCI control in the group of users who found it difficult to achieve BCI accuracy. Therefore, we chose the same architecture for the present study, in which the subjects undertook a BCI session for the first time, to induce fast BCI learning.

EEG was recorded by 62 Ag/AgCl electrodes concentrated on the central areas using a commercial EEG system (BrainAmp, Brain Products, Germany). Mastoid electrodes were used as reference and the ground electrode was placed on the frontal part of the EEG cap. Electrode impedances were less than 5 kΩ and the data were downsampled to 100 Hz before processing. EEG was measured in a noise‐protected and electrically shielded room while the users were sitting in an armchair in front of a computer screen.

For each subject of both groups, the first‐level general linear model (GLM) contained the three different motor imagination conditions of both 10‐min scans (i.e., ‘left hand’, ‘right hand’, ‘feet’). For each condition, one ‘imagination’ regressor was included in the GLM. For each regressor, the imagination was modelled with the boxcar function covering the 20 s cue duration. These boxcar functions were convolved with the standard haemodynamic response function (HRF) as implemented in SPM 8. The long inter‐trial intervals of 20 s were not explicitly modelled with the first level GLMs and hence represented an implicit baseline measure. For each of the three imagination conditions we computed individual β‐maps. Functional sessions pre‐ and post‐BCI training were analysed separately. On the second level (group) analyses, the main effect of each of the three imagination conditions (i.e., ‘left hand’, ‘right hand’, ‘feet’) was visualized by applying the individual β‐maps (pre‐BCI‐training sessions of both groups together) to one‐sample t tests to compare them against the null hypothesis. For the analysis of the pre, post and group comparisons, we applied a flexible factorial design with the factors ‘group’ (MI‐BCI vs . ERP‐BCI) and ‘time’ (pre vs . post) within SPM8.

The two results of the between‐group comparisons of (1) the difference T1‐weighted maps and (2) the corresponding difference ECM maps (RestPost3 – RestPre1) showed overlapping effects in the precuneus. This region of interest (ROI) was used to perform a whole‐brain seed‐based functional connectivity analysis for the ‘RestPost3’ and ‘RestPre1’ block. Custom‐built MATLAB scripts were used to calculate spatial correlation maps between the seed region and all other voxels in the GM mask. Paired t tests were used for within‐group comparison of the resulting spatial correlation maps. The between‐group comparison was performed by contrasting the difference between correlation maps (RestPost3 – RestPre1) of the two groups in a two‐sample t test. Monte Carlo simulations with 1000 iterations were used to correct all resulting T‐maps for multiple comparisons.

For analysing functional connectivity, we first used eigenvector centrality mapping (ECM) as a data‐driven approach that characterizes each voxel's connectivity within the whole‐brain GM mask (Lohmann et al . 2010 ; Nierhaus et al . 2012 ). One centrality map is calculated for each resting‐state block and Z ‐standard transform (i.e. for each voxel, subtract the mean value of the whole brain then divide by the standard deviation of the whole brain) was performed on the individual ECM maps (Zuo et al . 2012 ). Within‐group comparisons of the resulting centrality maps were performed using paired t tests (e.g. contrasting RestPost2 and RestPost1 for MI‐BCI and ERP‐BCI group, respectively). The between‐group comparisons were performed by contrasting the difference between ECM images (e.g. RestPost3 – RestPre1) of the two groups in a two‐sample t test. All resulting T‐maps underwent a multiple comparison procedure for significance thresholding based on Monte Carlo simulations with 1000 iterations (vmulticomp, LIPSIA toolbox; Lohmann et al . 2001 ; Nierhaus et al . 2015 ).

Pre‐processing of the resting‐state and functional MRI data was performed using the LIPSIA toolbox (Lohmann et al . 2001 ) including head motion correction, slice time correction, high pass filtering (at 1/100 s), and spatial smoothing (7 mm kernel). All images were co‐registered to the individual T1‐weighted structural image. For normalization to MNI space we used the same non‐linear deformation matrix which resulted from the DARTEL analysis of the T1 weighted images. CompCor analysis was done using the DPABI toolbox (toolbox for Data Processing & Analysis of Brain Imaging, http://rfmri.org/dpabi ) within the CSF/white matter mask (Behzadi et al . 2007 ), and the first four principal components together with six head motion parameters were used as nuisance signals to regress out associated variance.

T1‐weighted images were analysed using the DARTEL (Diffeomorphic Anatomical Registration Through Exponentiated Lie algebra) approach (Ashburner, 2007 ) as implemented in the Statistical Parametric Mapping software version 8 (SPM 8, Wellcome Trust Centre for Neuroimaging, University College London, London, UK), running in Matlab 8.2 (Mathworks Inc., Natick, MA, USA). All images were bias‐field corrected and segmented according to tissue probability maps incorporated in SPM8. Using DARTEL, all images were registered to an average size template generated from their own mean. The template was also used for registration to standard Montreal Neurological Institute (MNI) space. Finally, grey matter (GM) images were modulated with the Jacobian determinants to preserve the total amount of signal and smoothed using a Gaussian kernel of 8 mm full width at half‐maximum (FWHM). For each T1‐scan, this preprocessing produces one map containing the probability (values between 0 and 1) of each voxel belonging to GM. To ensure that further analysis was restricted to GM, threshold masking was applied to the GM images by excluding all voxels with probability below 0.2. For statistical analysis, we first performed a within group paired t test contrasting the GM images from pre‐ and post‐MRI sessions for the MI‐BCI and ERP‐BCI groups, respectively. Family‐wise error correction at cluster size level ( P < 0.05) was used to control for multiple comparisons. The between‐group comparison was performed by contrasting the difference between T1‐weighted GM maps (post‐pre) of the two groups in a two‐sample t test. The resulting T‐maps underwent a multiple comparison procedure for significance thresholding based on Monte Carlo simulations with 1000 iterations (vmulticomp, LIPSIA toolbox (Lohmann et al . 2001 ; Nierhaus et al . 2015 )).

The experiment consisted of three parts: ‘pre‐MRI’, ‘EEG‐BCI’ and ‘post‐MRI’. The experimental procedure was the same for both groups, only differing in the type of BCI task (MI‐BCI group vs . ERP‐BCI group). An MRI session including structural and functional scans was performed before and after conducting the EEG‐BCI session. The ‘pre‐MRI’ session started with a T1‐weighted anatomical scan followed by resting‐state and functional EPI scans. In the ‘post‐MRI’ session, the T1‐weighted anatomical scan was measured at the end (see Fig. 1 ). During MRI scanning, subjects were instructed to relax with eyes open. Resting‐state scans lasted 7 min each and were collected to investigate the effect of the BCI training on functional connectivity. ‘Functional blocks’ had durations of 2 × 10 min and consisted of a simple motor imagery task without feedback. This task was the same for both groups, as we wanted to compare the impact of the two different BCI trainings on evoked BOLD activity. Subjects were instructed to imagine movement of the right hand, left hand, or feet, indicated by a visually displayed arrow pointing right, left, or down, respectively. The three conditions were randomly presented in a 20‐s block design, intermingled with 20 s of rest. In a 10‐min block, each condition was presented five times. After the ‘pre‐MRI’ session, subjects were taken to the EEG room and fitted with the EEG cap for the BCI experiment. Overall, the experiment lasted approximately 4 h.

The study was designed and conducted according to the Declaration of Helsinki (except for registration in a database) and was approved by the local ethics committee of the University of Leipzig (Reference no. 953). Twenty‐three healthy subjects (11 females, mean age 26.7 years, SD = 4.1) participated in the motor imagery BCI experiment (MI‐BCI group). Two participants were eliminated from the later analysis, as they failed at performing the motor imagery task. Thus, 21 subjects were included for further analysis. Nineteen healthy subjects (9 females, mean age 27.2 years, SD = 3.1) performed a visual speller BCI experiment (ERP‐BCI group). All participants were right‐handed and gave written informed consent to participate in the experiment. Prior to the experiment, all subjects underwent a neurological examination and confirmed they were not taking any medication.

A and B , clusters of significant differences in correlation to the seed region, resulting from a paired t test contrasting correlation maps of ‘RestPost3 vs . RestPre1’ block within the MI‐BCI ( A ) and ERP‐BCI ( B ) groups, respectively. C , significant difference when comparing MI‐BCI and ERP‐BCI groups (two‐sample t test, RestPost3 – RestPre1). Warm colours indicate an increase and cold colours a decrease in correlation to the seed region. L, left hemisphere. Please note: this is an explanatory analysis given that the seed area showed statistically significant changes in the ECM analysis. Thus, the statistical values here are not valid but rather are meant to highlight a qualitative difference.

While the ECM comparison between two resting‐state runs identifies the voxel for which the whole‐brain correlation is altered, it does not show which specific connections are affected (Nierhaus et al . 2015 ; Antonenko et al . 2018 ). In order to identify potentially differential specific correlation changes between the MI‐BCI versus the ERP‐BCI group, we performed a seed‐based analysis of the region that showed the strongest consistent contrast in ECM and T1‐weighted analysis between the two groups. The overlap of the two ‘between‐group‐comparisons’ (MI‐BCI vs . ERP‐BCI group) for MPRAGE scans and the corresponding ECM comparison (Fig. 5 ) reveals a ROI in the precuneus that shows increased T1‐weighted signal and increased ECM for the MI‐BCI group. This ROI was used for a seed‐based correlation analysis of the resting‐state data. The comparison of the resulting spatial correlation maps for RestPost3 – RestPre1 shows increased correlated activity between the precuneus and the bilateral sensorimotor system for the MI‐BCI group, but decreased correlation between these areas for the ERP‐BCI group (Fig. 9 ).

In another analysis, we correlated the MRI data with EEG data that contained more individual BCI‐feedback information. In particular, we calculated band power features on the left hemisphere for BCI runs 1 and 4 in the same way as was done during the experiment, i.e. we used the subject‐specific spatial information and frequency band, and thus the data were no longer restricted to the mu‐band. Similar to the previous analysis, we found a significant correlation with the T1‐weighted signal in the right sensorimotor cortex, but no significant correlation for the ECM data (data not shown).

For the MI‐BCI group, the change in the sgn‐r 2 values (which expresses the discriminability of the two tasks in the EEG) from start (run 1) to end (run 4) of the BCI session shows a correlation with the change in EC from pre to post BCI (RestPre3 vs . RestPost1, Fig. 8 , upper plot, red‐yellow), and a correlation with the change in T1‐weighted MR intensity from start to end (Fig. 8 , upper plot, blue‐green) in the sensorimotor cortex. The sgn‐r 2 value is based on the difference in mu‐power between the two classes. The change in the mu‐power difference from run 1 to run 4 in the BCI also shows a significant correlation with ECM and T1‐weighted signal in the sensorimotor cortex (Fig. 8 , lower plot). This correlation is directly related to the change in mu power over the left hemisphere during the BCI session (compare Fig. 3 ).

When comparing the functional sessions pre vs . post BCI, we found increased BOLD activation only for the MI‐BCI group in the left sensorimotor cortex (right hand area) for the imagination condition ‘right hand’ (Fig. 4 , lower plots). The imagination conditions ‘left hand’ and ‘feet’ showed no significant BOLD activation and we did not find significant differences for the ERP‐BCI group. Furthermore, there was no interaction effect of the factors ‘group’ × ‘time’, namely the difference in the functional sessions from pre to post BCI was not significant for the comparison between MI‐ and ERP‐BCI group. The main effects (BOLD activation) for the three imagination conditions are shown in Fig. 7 .

In addition, for the MI‐BCI group we found increased centrality in the sensorimotor areas after performing the motor imagery task in the scanner, but not for the ERP‐BCI group (Fig. 6 ). While the comparison of RestPre2 – RestPre1 (before BCI session, Fig. 6 , left) shows only moderate changes in eigenvector centrality for both groups, the group comparison of RestPost2 – RestPost1 (after BCI session, Fig. 6 , right) shows increased EC for the MI‐BCI group in the left inferior parietal/superior temporal gyrus, left inferior frontal gyrus and left SII/operculum (the latter also overlaps with the T1‐weighted effect found for the MI‐BCI group, Fig. 4 , upper left plot).

We used eigenvector centrality mapping (ECM) as a data‐driven analysis tool to describe whole‐brain functional connectivity changes between two resting‐state runs without any prior assumption. Comparing the resting‐state scans acquired directly before and after the BCI session (RestPost1 vs . RestPre3, Fig. 4 , 2nd row) shows increased EC in left and medial sensorimotor cortices (right hand and feet areas), as well as decreased EC in the precuneus and prefrontal cortex for the MI‐BCI group, and increased EC in the left inferior frontal gyrus, cerebellum and middle temporal gyrus, as well as decreased EC in the occipital lobe and anterior cingulate cortex for the ERP‐BCI group. Comparison of the resting‐state scans performed alongside the MPRAGE‐anatomical scans (RestPost3 vs . RestPre1, Fig. 4 , 3rd row) shows increased EC in the precuneus for the MI‐BCI group. For the ERP‐BCI group this contrast showed decreased EC in the precuneus and medial sensorimotor cortex, as well as increased EC in the left middle temporal gyrus and bilateral inferior frontal gyrus. The corresponding between‐group comparison (MI‐BCI vs . ERP‐BCI group) for the ‘difference ECM images’ (RestPost3 – RestPre1) shows a significant cluster of larger centrality for the MI‐BCI group bilaterally extending over the precuneus, cuneus and calcarine gyrus (Fig. 5 , lower plot). This cluster overlaps with the precuneus cluster from the between‐group comparison of the T1‐weighted analysis (Fig. 5 , upper plot). This overlap was later used as a ROI for an exploratory seed‐based analysis of the resting‐state data.

Upper plots: clusters of significantly increased T1‐weighted MR intensity when comparing post‐MPRAGE. pre‐MPRAGE scans (pairedtest, glass brain illustration, FWE cluster corrected0.05). See Fig. 5 for between‐group comparisons. Second row: significant clusters of modulated eigenvector centrality observed by comparing the resting‐state scans directly before and after BCI (RestPost1RestPre3). Third row: significant clusters of modulated eigenvector centrality observed for the resting‐state scans performed alongside the MPRAGE‐anatomical scans (RestPost3. RestPre1). Warm colours indicate increased EC and cold colours indicate decreased EC from pre to post BCI. See Fig. 6 for modulation of ECM by performance of the MI task in the scanner. Lower plots: increased BOLD response in the functional sessions for the motor imagery condition ‘right hand’ (pre. post BCI training, FWE cluster corrected0.05). See Fig. 7 for the BOLD activation of the three imagination conditions. L, left hemisphere.

The majority of the participants in the MI‐BCI group showed a significant performance increase, and all users except one performed above the criterion of control (70%) in the last run. This ‘accuracy improvement’ is also a consequence of the algorithms used by the BCI system and their ability to increase the separability of the signals during the ongoing experiment. The use of co‐adaptive BCI systems implies that the algorithms continuously change during the experimental session and thus, in principle, an improvement could be entirely due to a more efficient algorithm, that is, the performance increase would only be due to the machine, not to the user. Therefore, we tested for changes related to the user by analysing feature shifts between the first and last runs. For that, the features (individual logarithmic frequency band power spatially filtered with subject‐optimized CSP) were projected on a subspace where the angles of the original space were preserved. In this manner, we could analyse the translation of the features in the space between first and last runs. More details of these procedures are available in Shenoy et al . ( 2006 ). In 19 of the 21 subjects, a feature shift was observed. Additionally, an analysis of the discriminability of EEG signals without using classifier information was also carried out. This procedure (see the Methods section) is a broad‐band (5–30 Hz) spectral analysis performed on subject‐specific optimal channels for discrimination of the two classes. The result of grand‐averaging the spectra of all participants is visible in Fig. 3 (left). We could observe an improvement of discriminability in the mu‐band from first to last run (Fig. 3 , left panel). That is, the sgn‐r 2 value is larger in run 4 compared to run 1, indicating that the quality of the signals, in terms of discriminability, had improved. The mu band was therefore selected for subsequent correlation analyses with MRI data. For more information, on the right panel of Fig. 3 , the number of selections of each channel in the left hemisphere for the MI‐BCI group is depicted.

The online BCI performance averaged across all users shows an accuracy increase for the MI‐BCI group over time (left panel, Fig. 2 ), whereas the ERP‐BCI group rapidly achieved good performance, even in the first run (right panel, Fig. 2 ). Each horizontal bar represents the mean value for one run. For the MI‐BCI group each point is the mean value across 20 trials. Users are considered to have control if they reach an average performance of 70% in the last run, in accordance with the results of Kübler et al . ( 2004 ) for a two‐class BCI. For each user, binomial tests on the hit/miss distribution of the last run were performed and compared to the accuracy reached in the first run. For the MI‐BCI group there was a statistically significant increase for fourteen users. Six participants performed similarly well in the first and last run, and one user performed worse in the last run. Only one user performed significantly worse than 70% (criterion of control) in the last run. For the ERP‐BCI group, during the last run six users significantly outperformed their expected performance in run 1, nine performed similarly well, and four performed significantly worse.

Discussion

In two different types of BCI, we demonstrate changes in functional and structural MRI after only 1 h of BCI. For the visual ERP‐BCI group, structural MRI shows increased T1‐weighted intensity of grey‐matter in occipital/parietal areas, for the MI‐BCI group in left sensorimotor, medial parietal and left occipital areas. Additionally, we find increased eigenvector centrality (EC) in left sensorimotor areas and precuneus for the MI‐BCI group, and increased EC in the inferior frontal gyrus as well as decreased EC in precuneus, occipital cortex and medial sensorimotor cortex for the ERP‐BCI group. For the MI‐BCI group, we observed increased central background (mu) rhythm activity, and the mu‐rhythm change as well as the MI‐BCI task discriminability (sgn‐r2 value measured in the EEG) shows a strong correlation with T1‐weighted and ECM change in sensorimotor cortices. Furthermore, when using the subject specific features of run 4 and run 1, we observed a similar correlation with T1‐weighted MRI data in sensorimotor cortex, but no correlation for ECM.

We regard our findings as indicators of BCI‐induced brain plasticity and, as subsequently outlined, even the very early time point of the observed changes seems to be consistent with this notion: Confirming invasive animal studies, an abundance of non‐invasive MRI studies have provided evidence that the adult human brain is capable of functional and structural reorganization (Draganski & May,2008). Cross‐sectional studies have indicated changes in grey‐matter‐density and/or volume in certain brain regions presumably engaged in long‐term activities, such as increased hippocampus size in London taxi drivers (Maguire et al. 2000), alterations in planum temporale and corpus callosum of musicians (Schlaug et al. 1995a,b), and in the striatum of professional athletes (Taubert et al. 2015). The correlational evidence from cross‐sectional studies has been supported by causal evidence from interventional studies using functional MRI (Karni et al. 1995; Poldrack et al. 1998; Fletcher et al. 1999; Pleger et al. 2003; Limanowski et al. 2017), but also structural MRI (Draganski et al. 2004; Taubert et al. 2010, 2012; Zatorre et al. 2012). While early intervention studies employing structural MRI (T1‐weighted) reported ‘brain plasticity’ within weeks/months of training onset (Draganski et al. 2004; Taubert et al. 2010), more recent studies have reported changes within few hours after task performance or sensory stimulation (Sagi et al. 2012; Hofstetter et al. 2013; Taubert et al. 2016; Schmidt‐Wilcke et al. 2018). Interestingly, the locations of these early changes are not necessarily the same as for long‐term effects (Taubert et al. 2016). It seems that different brain areas and networks are involved during different phases of learning (multiphase‐model of brain plasticity; Taubert et al. 2012). Furthermore, (T1‐weighted) MRI findings at different time points after training onset probably do not reflect the same underlying neurophysiological/cellular processes. Long‐term plasticity changes are characterized by persistent structural changes, i.e. increased number of axons, dendrites, synapses, glial cells and vessels, as well as (perhaps) neurogenesis in some areas (e.g. hippocampus, striatum). These long‐term neuroanatomical changes most likely correspond to changes in T1‐weighted MRI; however, they occur only within days/weeks (Anderson et al. 2002; Kleim et al. 2004) and thus are unlikely to underlie the very early MR‐signal changes. For rapid structural plasticity, the direct link to the MR‐signal is still poorly understood. It is known that within the first hour of training, new spines are already being generated (Xu et al. 2009; Moczulska et al. 2013; Kuhlman et al. 2014). But due to their extremely small size it is rather unlikely that this can explain the ‘macroscopic’ MR‐signal changes. However, microglial invasion in the surrounding of new synapses has been shown to be a very rapid early process (Sagi et al. 2012; Kettenmann et al. 2013; Tremblay et al. 2010). Also, the motility of synapse‐associated astrocytic processes has been shown to be enhanced as early as 5 min after stimulus‐induced long‐term potentiation (Perez‐Alvarez et al. 2014). As neuroglial transmission has been shown to be fundamental in maintaining learning‐induced cortical stability (long‐term potentiation) (Jammal et al. 2018), such rapid neuroglial processes, together with local increases in metabolism, seem likely to explain early MR‐signal changes.

Based on these considerations, we regard it as plausible that our findings reflect this early phase of brain plasticity, especially because potential confounders for the interpretation of ‘structural MRI’ (i.e. influences of blood flow, blood volume and/or tissue oxygenation on T1‐relaxation time) cannot explain our findings. We have previously shown that rapid changes in T1‐weighted MR‐intensity (induced after 1 h of balance training) are not accompanied by major blood flow changes (Taubert et al. 2016). But what about (minor) local changes in blood flow or oxygenation? Recently, it has been shown that a breathing‐induced increase in the concentration of oxygen decreases T1‐relaxation time, the so‐called tissue oxygenation level‐dependent (TOLD) contrast (Haddock et al. 2013; Tardif et al. 2017), and likewise hypercapnia decreases T1 (Tardif et al. 2017). These effects may be due to changes in blood flow/volume and/or tissue oxygenation. When performing straightforward analyses of cortical thickness and grey‐matter density, both hypercapnia and ‘breathing oxygen’ lead to ‘apparent increases’ (Tardif et al. 2017). Could such an effect explain our findings? The post T1‐weighted scan was performed more than 1 h after the BCI‐session and about 15 min after the last functional MRI run. Since task‐induced changes in blood flow and oxygenation normalize after a few minutes (Bandettini et al. 1997; Howseman et al. 1998), they are unlikely to explain the persistent increase in T1‐weighted signal. Furthermore – keeping in mind that the ‘last task 15 min before the structural MRI’ was a sensorimotor task in all subjects – even if task‐induced blood flow/volume increases lasted longer than 15 min, there should be sensorimotor related changes in both groups (MI‐BCI and ERP‐BCI). However, the observed changes were noted in visual areas for the ERP‐BCI group and sensorimotor areas for the MI‐group, i.e. areas which are specifically related to the BCI‐task (more than 1 h before the structural scan) but NOT related to the functional task in the MR‐scanner. It is therefore much more likely that the increased T1‐weighted MR‐intensity, which we find for the respective groups in different areas, reflects early events of brain plasticity such as local accumulation of microglia or long‐lasting plasticity‐associated increase in metabolism induced by the BCI tasks.

For the MI‐BCI group, we provide evidence that this brain plasticity is related to modulation of background rhythms during the BCI session. The performance increase during MI‐BCI is related to ‘the adaptation of the algorithm to the subject’ as well as ‘the subject's ability to modulate their sensorimotor rhythms’. Overall, the co‐adaptive set‐up leads to a performance increase. Additionally, subjects’ EEG data show increased background rhythm power at the end of the BCI session, especially in the mu‐frequency band over left central areas. This synchronized oscillatory activity may reflect BCI‐induced recruitment of additional neurons that are coherently involved in generating the oscillatory pattern and, after some time, would strengthen their connections (Ros et al. 2014). Furthermore, our finding of increased BOLD activity in the left hemisphere after the MI‐BCI session hints at the recruitment of ‘more neurons’ during imagination of right hand movements, in contrast to the ‘left hand’ condition, which was not trained in the BCI and thus does not show a change in BOLD activity. Moreover, the modulation of sensory background rhythms during BCI is directly related to BCI performance. Thus, the correlation between changes in EEG‐BCI data (mu‐rhythm increase as well as BCI discriminability improvement) and changes in MRI data (T1‐weighted MR‐intensity as well as ECM), which we found in sensorimotor areas, is consistent with a role of the underlying background rhythm modulations in inducing brain plasticity. These findings are in line with previous evidence of homeostatic plasticity found after alpha suppression trained by means of neuro‐feedback (Ros et al. 2014, 2017). Among the two MR‐correlates of BCI‐induced neural changes (T1, ECM), ECM changes seem to be more closely related to mu‐rhythm, as no significant correlation was found when the EEG data are not restricted to the mu‐frequency band, whereas T1‐weighted MRI changes also showed for the subject‐specific rhythms a strong correlation in the sensorimotor cortex, similar to those related to changes in mu‐rhythm and BCI discriminability.

Our findings of increased T1‐weighted intensity and EC in the precuneus, as well as the increased temporal correlation of the BOLD signal time courses between the precuneus and sensorimotor system for the MI‐BCI group, are in line with previous studies showing that the precuneus plays an important role for visuospatial integration, e.g. hand‐eye coordination (Ferraina et al. 1997) and reaching (Caminiti et al. 1999). Also, increased precuneus activation during motor imagery (vs. real motor execution) was shown for finger movements (Gerard et al. 2000) and sequential finger tapping (Hanakawan et al. 2003). An MEG study revealed precuneus dipole activity for motor imagery that suggested a signal transfer to supplementary motor areas (Ogiso et al. 2000). It seems that the information exchange between precuneus and sensorimotor areas is necessary for motor imagery, independent of perceptual motor feedback. This might explain our finding of enhanced functional connectivity between these brain areas after MI‐BCI, whereas visual ERP‐BCI seems not to depend on this connection.

Recently, modulation of resting‐state functional connectivity was demonstrated after 2 weeks of motor imagery training, in terms of modulated default‐mode network connectivity (Ge et al. 2015), as well as increased connectivity strength of the sensory‐motor network in the precuneus and fusiform gyrus (Zhang et al. 2014). Interestingly, significantly increased connectivity in the fusiform gyrus does not appear in our data, but we found increased T1‐weighted signal for the MI‐BCI group in this area. It seems that the fusiform gyrus is engaged in the MI‐BCI task, and thus the T1‐weighted signal is modulated in our data presumably due to persistent changes in metabolism.

The increased T1‐weighted MR‐intensity in occipital areas and increased EC in the inferior frontal gyrus (IFG) that we found for the ERP‐BCI group are probably explained by the performance of the visual spelling task. The middle occipital gyrus (BA18 – secondary visual cortex) not only processes visual properties of shape (Hegdè & Essen,2000) or orientation (Anzai et al. 2007), but also activates with sustained (visual) attention (Le et al. 1998) and, together with IFG, is involved in word retrieval processes (Abrahams et al. 2003). The left IFG especially is related to linguistic processes such as syntax (Tyler et al. 2011), lexical search (Heim et al. 2005), or word generation (Friedman et al. 1998). The right IFG seems to be involved in target detection (Hampshire et al. 2009; de Haan et al. 2015), but also contributes to language function (Hartwigsen et al. 2013). The increased T1‐weighted MR‐intensity that we find in left occipital areas for the MI‐BCI group might be due to the lateralization of the MI‐task (training right hand and feet associated with cursor movement to the right and down, respectively), i.e. ‘successful’ imagination of right hand movement results in a right‐directed cursor movement on the screen (visual stimulation in the right hemifield), thus involving left occipital areas.

In conclusion, in young and healthy volunteers, our study demonstrates rapid effects of purely mental BCI tasks on MRI parameters that are thought to reflect structural and functional brain plasticity. Future studies need to identify how these early changes develop into long‐term neural reorganization, and how BCI‐induced brain plasticity may be therapeutically useful in patients with brain disorders, as suggested by interesting pilot studies using BCI in stroke patients (Varkuti et al. 2013; Young et al. 2014). The spatial specificity of the BCI effects that we observe in this study creates an opportunity for tailoring BCI‐based therapeutic approaches individually to, for example, stroke patients, according to the individual patient's lesion location.