Neural Correlates of Motor Recovery after Robot-Assisted Training in Chronic Stroke: A Multimodal Neuroimaging Study

Stroke is a leading cause of motor disability worldwide, and robot-assisted therapies have been increasingly applied to facilitate the recovery process. However, the underlying mechanism and induced neuroplasticity change remain partially understood, and few studies have investigated this from a multimodality neuroimaging perspective. The current study adopted BCI-guided robot hand therapy as the training intervention and combined multiple neuroimaging modalities to comprehensively understand the potential association between motor function alteration and various neural correlates. We adopted EEG-informed fMRI technique to understand the functional regions sensitive to training intervention. Additionally, correlation analysis among training effects, nonlinear property change quantified by fractal dimension (FD), and integrity of M1-M1 (M1: primary motor cortex) anatomical connection were performed. EEG-informed fMRI analysis indicated that for iM1 (iM1: ipsilesional M1) regressors, regions with significantly increased partial correlation were mainly located in contralesional parietal, prefrontal, and sensorimotor areas and regions with significantly decreased partial correlation were mainly observed in the ipsilesional supramarginal gyrus and superior temporal gyrus. Pearson's correlations revealed that the interhemispheric asymmetry change significantly correlated with the training effect as well as the integrity of M1-M1 anatomical connection. In summary, our study suggested that multiple functional brain regions not limited to motor areas were involved during the recovery process from multimodality perspective. The correlation analyses suggested the essential role of interhemispheric interaction in motor rehabilitation. Besides, the underlying structural substrate of the bilateral M1-M1 connection might relate to the interhemispheric change. This study might give some insights in understanding the neuroplasticity induced by the integrated BCI-guided robot hand training intervention and further facilitate the design of therapies for chronic stroke patients.


fMRI Preprocessing
The resting-state fMRI data were preprocessed using the Data Processing Assistant for Resting-State fMRI (DPARSF) toolbox [1] based on Statistical Parametric Mapping (SPM12) (http://www.fil.ion.ucl.ac.uk/spm). 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 Friston 24 head motion parameters [2]. To further control for head motion, the scrubbing process was performed for the volumes with frame-wise displacement value exceeding 0.3 [3]. If over 25% of all the volumes exceed the frame-wise displacement threshold, the data for this subject would be discarded, and no subject was discarded in resting-state fMRI analysis. Then the functional dataset was aligned to the anatomical dataset. Detrending and temporal band-pass filtering (0.01-0.1 Hz) [4,5] 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, resliced to 2×2×2 mm 3 voxels, and smoothed with a Gaussian kernel with a fullwidth at half-maximum (FWHM) of 6 mm. In the statistical analysis, subjects with left-hemispheric lesions were flipped along the midsagittal plane using MRIcron (www.mccauslandcenter.sc.edu/mricro/mricron). This procedure was also adopted in previous studies [6,7], so that the lesions of all subjects were in the right hemisphere.

EEG Preprocessing
EEG data were mainly preprocessed with EEGLAB [8], Fieldtrip toolbox [9], Matlab signal processing toolbox and custom-made codes (Mathworks Natick, MA, USA). The EEG signal had a relatively low signal-to-noise ratio (SNR) when simultaneously collected with fMRI due to the switching of magnetic field gradients. Therefore, a principal component analysis (PCA)-based optimal basis set (OBS) algorithm [10] was utilized to remove the MRI gradient artifact. The triggers marking the onset of each fMRI volume were also provided to better select and extract the artifactual features. In our study, the missing triggers for two subjects were repaired manually by an expert radiographer. The resulting EEG signal was inspected again to ensure no massive residual gradient artifact was observed. Ballistocardiographic (BCG) artifact related to heart pulse was another source distorting the EEG. To remove the artifact, the QRS complexes of each heart pulse were initially automatically determined with an R-peak detection algorithm. Then, a strategy developed by Liu and colleagues, which combined independent component analysis (ICA), OBS and an information-theoretic rejection criterion, was adopted to eliminate the BCG artifact channel-wisely [11]. After gradient and BCG artifact removal, the EEG signals were resampled to 250 Hz and further band-pass filtered with a 2-40 Hz Butterworth non-causal filter. Bad channels with extensive artifacts were removed and reconstructed through spherical spline interpolation with neighbor electrodes. All data were then common average referenced. Subsequently, the EEG data were segmented into non-overlapping twosecond epochs according to the triggers from MRI scanner. The first and last four epochs were dropped, considering the signal instability. The remaining epochs were inspected visually accompanied by statistical metrics (e.g., z-score, variance, min and max, etc.), and bad epochs were rejected. Finally, an adaptive mixture independent component analysis (AMICA) algorithm [12] was utilized to wipe off the components potentially related to residual artifacts induced by MRI scanning, eye movements, and muscle contraction. The time courses of remaining components were projected back to all channels. To run appropriately the group-level analysis, EEG data for patients with left-hemispheric lesions were left-right flipped along the middle line before signal processing procedures. (Corresponding to step 5 in Figure 1.)

EEG distributed source estimation
The preprocessed EEG data were reconstructed into the distributed source space using Brainstorm software [13]. First, EEG electrode positions were manually co-registered to the individual T1 images based on the standard 10-20 EEG electrode locations and polygon models of the cortical surface were constructed. The three realistic layers (scalp, inner and outer skull) were extracted from individual MRI data to generate the boundary element model (BEM) surface. The gain matrix, called the lead field, representing the linear relationship projecting the source time courses on the cortical surface to the measurements on EEG channels, was computed using OpenMEEG [14]. The Tikhonovregularized minimum-norm algorithm was adopted to perform EEG source localization [15], with Tikhonov parameter which regularized the smoothness of the images set to 10% of maximum singular value of the lead field. Herein, one constraint was imposed that only a single current dipole was assumed at each vertex point and oriented normally to the cortical surface based on the anatomical observation. The noise covariance matrix was calculated over a long period of recordings, as recommended by Tadel and colleagues [13]. (Corresponding to step 6 in Figure 1.) To investigate the neural activities from primary motor areas, two spherical seeds with a radius of 5 mm at ipsilesional M1 (iM1) and contralesional M1 (cM1) were determined. The iM1 and cM1 seeds were set at (38,-22,56) and (-38,-22,56) in the MNI space, respectively. The standard seeds in MNI space were wrapped to subject individual cortex space and the time series of vertices within the seed were averaged. The extracted time series from iM1 and cM1 seeds were mainly explored in further analysis. (Corresponding to step 7 in Figure 1.)

EEG-informed fMRI analysis
The extracted EEG source time courses in iM1 seed for each subject were further used in EEG-informed fMRI analysis [16,17]. Time-frequency power spectrum was estimated for each extracted EEG source time course using short-time Fourier transform and the length of the sliding window was 2 seconds without overlap between adjacent windows. The EEG bands were defined as theta (4-7 Hz), alpha (8)(9)(10)(11)(12), and beta (13-30 Hz). Band power signals were calculated by averaging the power spectrum within each EEG band. The band power signals were further convolved with the hemodynamic response function (HRF) in SPM12, to obtain the EEG regressors of interest in partial correlation analysis [18,19]. (Corresponding to step 12 in Figure 1.) For each EEG band and each subject, general linear models (GLM) were fitted for every voxel in the brain with the corresponding EEG regressor and six motion parameters in fMRI analysis as covariables to control the effects of head motion. The GLM model was formed as follows [20]: where Y represented the measured fMRI signal of a single voxel. X was the design matrix. In this study, the first column of X was EEG regressor and the remaining ones were composed of six motion parameters. β was the weighting vector which was to be estimated and the entry corresponding to the EEG regressor would be treated as the partial correlation for this voxel. was the random error term. The map of partial correlation would be produced after traversing all voxels. (Corresponding to step 11 in Figure 1.)

Conventional seed-based fMRI connectivity analysis
Except for EEG-informed fMRI analysis, conventional seed-based fMRI functional whole-brain connectivity analysis was also performed. The same seeds of iM1 in MNI space were used, and the average time course of the BOLD signal within the seeds was calculated as the regressors. Following that, functional connectivity maps were generated.

Fractal Analysis
The fractal analysis on the iM1 and cM1 source time series was characterized by fractal dimension (FD) [21], a measure which has been shown to be able to capture the non-linear information on scale-free properties and to extend the understanding of brain mechanisms at system levels [22,23]. Usually, brain damage could result in less complexity of the underlying neuronal activity, which is indicated in the decrease of FD value. Plenty of algorithms are available to calculate the FD and box-counting method [24] is a classical one among them. In our study, we adopted a strategy similar to Raghavendra and Dutt [25] to increase the box size from the finest time resolution to the maximum coarse time resolution with the increasing factor equal to two. The FD values were calculated in the 10-second sliding time windows with 50% overlap and averaged over all windows. The balance of homologous motor regions is of importance in rehabilitation for stroke patients with motor impairment. In order to investigate the complexity unbalance, F D asymmetry was computed to characterize the interhemispheric asymmetry of motor areas as follows: where F D cM 1 and F D iM 1 are FD values for cM1 and iM1 respectively. In this study, we utilized the pre-post change of FD asymmetry to characterize the process of interhemispheric rebalance after the training intervention. (Corresponding to step 8 in Figure 1.)

DTI analysis
Diffusion-weighted images before the training were used for further analysis and first preprocessed using FMRIB Software Library (FSL), including correction for eddycurrent effect and brain extraction. Diffusion tensors were then estimated from the preprocessed diffusion images. Fractional anisotropy (FA) [26], a measure quantifying the degree of anisotropic diffusion, was calculated from the diffusion tensors of each voxel. All the diffusion images were registered to MNI space. To assess the transcallosal connections between two hemispheres, especially the reserved connection between primary motor areas, a recently published standard template TCATT (Transcallosal Tract Template) [27] was used to map the structural integrity. Average FA values within the M1-M1 tract were computed as the structural integrity between the bilateral motor cortex. (Corresponding to step 1, 2 and 3 in Figure 1.)

Extraction of contralesional frontal-parietal network and sensorimotor network
The contralesional frontal-parietal network and sensorimotor network were extracted using independent component analysis (ICA). All fMRI scans from all sessions were combined and run in a group ICA using the temporal concatenation method implemented in MELODIC [28], constrained to identify 30 components. An initial analysis of the population using model order estimation was used to calculate the appropriate number of components. The contralesional frontal-parietal network and sensorimotor network were chosen based on previous studies [29].
4. Interhemispheric asymmetry before and after training Figure S2: The interhemispheric asymmetry values before and after training for each individual subject.

EEG-informed fMRI analysis using surface-based method
MRI data were preprocessed using the fMRIPrep software version 1.5.0 [30]. T1weighted images were corrected for intensity non-uniformity [31] and skull-stripped using ANTS 2.2.0 [32]. High resolution cortical surfaces were reconstructed with FreeSurfer 6.0.1 [33] based on T1-weighted images, and then registered to the fsaverage5 template [34]. Functional data was slice-time corrected using 3dTshift [35], motion corrected using MCFLIRT [36], distortion corrected, and then resampled to the fsaverage5 template based on boundary-based registration [37]. Nuisance variables were then regressed out, including white matter, cerebrospinal fluid (CSF), global mean signal, and Friston 24 head motion parameters [2]. The surface data were further smoothed within 6 vertices. Similar EEG-informed fMRI analysis as the volume-based analysis was performed on the preprocessed surface data. The two-way ANOVA with time (pre and post) and frequency bands as factors were performed. Vertex-level threshold was set at p < 0.005. No significant interaction result was found for iM1 seed. To compare with the volume-based results, paired t-tests were performed between pre and post training for 3 frequency bands (theta, alpha, and beta) and two hemispheres. Vertexlevel threshold was set at p < 0.005. Significant clusters were only observed for alphaband and theta-band for the iM1 seed. Consistent with the volume-based results, significant increased partial correlation was observed in the contralesional superior frontal area for alpha-band, and significant decreased partial correlation was observed in the ipsilesional supramarginal gyrus for theta-band. Besides, significant increased partial correlation was also observed in the contralesional superior frontal area for the theta-band. (illurstrated in Figure S3) Figure S3: The brain regions which showed significant partial correlation change before and after training using surface-based method, given the regressor embedding spectral information (including theta, alpha frequency bands) derived from EEG source time course of iM1.