BCI Training Effects on Chronic Stroke Correlate with Functional Reorganization in Motor-Related Regions: A Concurrent EEG and fMRI Study

Brain–computer interface (BCI)-guided robot-assisted training strategy has been increasingly applied to stroke rehabilitation, while few studies have investigated the neuroplasticity change and functional reorganization after intervention from multimodality neuroimaging perspective. The present study aims to investigate the hemodynamic and electrophysical changes induced by BCI training using functional magnetic resonance imaging (fMRI) and electroencephalography (EEG) respectively, as well as the relationship between the neurological changes and motor function improvement. Fourteen chronic stroke subjects received 20 sessions of BCI-guided robot hand training. Simultaneous EEG and fMRI data were acquired before and immediately after the intervention. Seed-based functional connectivity for resting-state fMRI data and effective connectivity analysis for EEG were processed to reveal the neuroplasticity changes and interaction between different brain regions. Moreover, the relationship among motor function improvement, hemodynamic changes, and electrophysical changes derived from the two neuroimaging modalities was also investigated. This work suggested that (a) significant motor function improvement could be obtained after BCI training therapy, (b) training effect significantly correlated with functional connectivity change between ipsilesional M1 (iM1) and contralesional Brodmann area 6 (including premotor area (cPMA) and supplementary motor area (SMA)) derived from fMRI, (c) training effect significantly correlated with information flow change from cPMA to iM1 and strongly correlated with information flow change from SMA to iM1 derived from EEG, and (d) consistency of fMRI and EEG results illustrated by the correlation between functional connectivity change and information flow change. Our study showed changes in the brain after the BCI training therapy from chronic stroke survivors and provided a better understanding of neural mechanisms, especially the interaction among motor-related brain regions during stroke recovery. Besides, our finding demonstrated the feasibility and consistency of combining multiple neuroimaging modalities to investigate the neuroplasticity change.

The first 10 volumes were discarded to assure the remaining volumes of fMRI data were at magnetization steady state. The remaining volumes were corrected with slice timing and realigned for head motion correction. Nuisance variables were then regressed out including white matter, cerebrospinal fluid (CSF), global mean signal and Fristion 24 head motion parameters [1]. To further control for head motion, scrubbing process were done for the volumes with framewise displacement (FD) value execeed 0.3 [2]. Then the anatomical dataset was aligned to the functional dataset. Detrending and temporal band-pass filtering (0.01 Hz -0.1 Hz) [3] were performed to remove higher frequency physiological noise and lower frequency scanner drift. Subsequently, the functional images were spatially normalized to the Montreal Neurological Institute (MNI) template (MNI152: average T1 brain image constructed from 152 normal subjects), resliced to 2 mm × 2 mm × 2 mm voxels, and smoothed with a Gaussian kernel with a full-width at half-maximum (FWHM) of 6 mm.

Defining the Seed Locations
When acquiring task-based fMRI, subjects were asked to try to do grasp and open using their paretic hand when a mark of "L" or "R" (decided by each subject's paretic hand) appeared on the screen and were also asked to maintain 6 seconds until the mark disappeared from the screen. An event-related design was adopted and randomized time intervals from 14 to 20 seconds were assigned between every two tasks. Two 6-minute task-based fMRI runs were performed for each subject.
The task-based fMRI data were also preprocessed using DPARSF toolbox. Similar preprocessing steps were performed on task-based fMRI data except that the threshold for FD value was set to 0.7 in the motion scrubbing step [4] and no band-pass filter was used. Subjects with left-hemispheric lesions were also flipped along the midsagittal plane. The preprocessed task-based fMRI data were fitted into a GLM for subject-level analysis where each of the events were convolved using a canonical hemodynamic response function (HRF) and used as a regressor. Besides, six motion regressors were also included in the design matrix to regress out motion-related fluctuations in the BOLD signal. As a result, each subject acquired a t-map at each session. The statistical maps from all sessions were used to do the group level analysis (one sample t-test). Voxels were identified as significantly active if they surpassed a threshold of z > 2.7 and corrected using Gaussian random field theory at a threshold of p < 0.05 at a cluster level. The seed locations for iM1 and cM1 were defined based on the group-level activation map.

EEG Data Preprocessing
Under the condition where MRI was acquired simultaneously, the switching of magnetic field gradients would pollute and overwhelm EEG signal which led to low signal to noise ratio (SNR). A principle component analysis (PCA)-based optimal basis set (OBS) algorithm [5] was adopted to remove the MRI gradient artifact and the onset markers indicating the beginning of each fMRI volume, generated by MRI scanner were also provided for better extraction and selection of artifactual features. The output EEG signal were double-checked visually to ensure that the amplitude was not grandiosely large. The time course of heartbeat artifact was determined with a R-peak detection algorithm [6]. The final ECG artifact was eliminated channel-wisely using the strategy which combined independent component analysis, OBS and an information-theoretic rejection criterion developed by Liu and colleagues [6].
After that, EEG signal was band-pass filtered from 2 to 40 Hz using a Butterworth non-causal filter. Subsequently, bad channels were removed and reconstructed using spherical spline interpolation with neighbor electrodes. Following that, all data were common average referenced. According to the fMRI trigger markers, these data were segmented into non-overlapping two-second epochs where the first and last several data segments were removed due to signal instability. Bad epochs were rejected based on statistical measurement metrics (e.g. z-score, variance, min, and max etc.) with remaining ones further inspected visually to guarantee the signal quality. We utilized adaptive mixture independent component analysis (AMICA) algorithm [7] to separate EEG signals into spatially static and maximally temporally independent components [8]. The components related to residual artifact induced by MRI scanning, Electrooculogram (EOG) artifact and muscular artifact were rejected. Processes of remaining components were then projected back to all original channels. Finally, we applied a surface Laplace filter with the spherical spline method [9] to increase the topographical selectivity, eliminate the volume conduction and highlight the high-spatial-frequency components while attenuating low ones [10].

Correlation of Information Flow When iM1 Was the Source Region
Then we explored the relationship between training effect and information flow changes when iM1 was treated as a source region (i.e. information flow from iM1 to cPMA or SMA). However, neither information flow change from iM1 to cPMA (r = −0.405, p = 0.7660, Bonferroni corrected) nor from iM1 to SMA (r = 0.1189, p = 1, Bonferroni corrected) correlated significantly or strongly with FMA score change. Figure S1. The correlation between functional connectivity change and information flow change from iM1 to cPMA (left), iM1 to SMA (right).