Removal of artifacts from resting-state fMRI data in stroke.

We examined the effect of lesion on the resting-state functional connectivity in chronic post-stroke patients. We found many instances of strong correlations in BOLD signal measured at different locations within the lesion, making it hard to distinguish from the connectivity between intact and strongly connected regions. Regression of the mean cerebro-spinal fluid signal did not alleviate this problem. The connectomes computed by exclusion of lesioned voxels were not good predictors of the behavioral measures. We came up with a novel method that utilizes Independent Component Analysis (as implemented in FSL MELODIC) to identify the sources of variance in the resting-state fMRI data that are driven by the lesion, and to remove this variance. The resulting functional connectomes show better correlations with the behavioral measures of speech and language, and improve the out-of-sample prediction accuracy of multivariate analysis. We therefore advocate this preprocessing method for studies of post-stroke functional connectivity, particularly in samples with large lesions.

Introduction

Functional magnetic resonance imaging (fMRI) provides a noninvasive indirect measure of brain activity (Logothetis and Wandell, 2004). Traditionally, each voxel of the brain is independently tested to see if the observed signal correlates with behavioral events (e.g. observed stimuli or responses; e.g., Worsley and Friston, 1995). However, one can also measure how the signals measured at different voxels correlate with each other. These analyses of functional connectivity (Friston, 1994) have become a popular approach to study the activity of cortical networks. Functional connectivity can be measured while the participant is engaged in a particular task (e.g., Friston et al., 1993), or when the participant is not given any particular instruction other than to stay still, yielding a measure of resting-state functional connectivity (Biswal et al., 1995). The latter procedure may be particularly useful for clinical studies involving participants who may not be able to carry out a given behavioral task while in the scanner or where some or all of the participants respond with such low accuracy that the functional MRI may not reflect task performance. Analysis of resting-state functional connectivity has been performed on data obtained from various clinical conditions such as Alzheimer's disease, attention deficit disorder, autism, schizophrenia, stroke as well as several other disorders (see review by Lee et al., 2013). In theory, this neuroimaging modality could have a major clinical impact for the diagnosis, prognosis and individualized treatment for neurological disorders.

Functional connectivity studies of stroke are common (see reviews by Grefkes and Fink, 2011, and by Thiel and Vahdat, 2015); however, brain lesions caused by stroke present some special challenges for functional connectivity analyses. The tissue at the site of the infarction, being deprived of blood supply, is at various stages of necrosis and gliosis; therefore, the local BOLD signal may not reflect meaningful neuronal activity. Currently, there is no consensus on if or how the signal measured in the frank lesion should be preprocessed before carrying out connectivity analyses. In healthy individuals, it is common to regress out signal variability in either white matter or cerebro-spinal fluid (Weissenbacher et al., 2009). In some stroke studies, no specific steps are reported to account for potentially aberrant signal in the stroke lesion (e.g., Nair et al., 2015, Marangolo et al., 2016, Nijboer et al., 2017, Zhang et al., 2017). In chronic participants, the site of the lesion is composed mostly of cerebro-spinal fluid; so, regressing out the CSF signal (e.g. see Tuladhar et al., 2013, van Hees et al., 2014, Xu et al., 2014, Bannister et al., 2015, Gili et al., 2016, Kielar et al., 2016, Sebastian et al., 2016, Siegel et al., 2016, Tang et al., 2016, Adhikari et al., 2017, Hu et al., 2017, Yang et al., 2017) might remove the variance in signal unrelated to brain activity. Yet another approach is masking out the site of the lesion during the computation of functional connectivity (e.g., Balaev et al., 2016, New et al., 2015).

In the present study, we investigated resting-state functional connectivity of 74 chronic stroke participants. Our preliminary analyses revealed that signal intensities measured at different locations within the site of infarction tended to be correlated. This yielded high values of functional connectivity between brain areas destroyed by the lesion, making it impossible to distinguish from functional connectivity measured between two intact areas that were part of the same cortical network. This, in turn, provided some paradoxical results when we analyzed the associations between functional connectivity and behavioral impairment. Regressing out the CSF signal did not alleviate the problem, yielding further paradoxical results.

We speculate that the presence of the lesion introduces a particular artefact into the fMRI data which is not removed by standard preprocessing techniques (e.g. CSF regression). We obtained notably better results when these lesion-related artifacts were removed from the fMRI data using Independent Component Analysis (ICA; Beckmann and Smith, 2004). This method separates the fMRI signal into different sources of variance, which are mutually independent in the spatial domain. ICA has been advocated for removing artifacts related to physiological noise (Perlbarg et al., 2007, Tohka et al., 2008, Griffanti et al., 2014). We noticed that in our sample of chronic stroke survivors, ICA identified the sources of variance driven by the lesion: the spatial maps of the corresponding independent components overlapped with manually-drawn lesion masks. When these “lesion-driven” components were filtered out of the fMRI data, we saw consistent improvement: (a) the BOLD correlations between lesioned areas decreased; (b) a larger number of functional connections were associated with behavioral impairment; and (c) predictive power of the functional connectome, as measured by multivariate analysis within a cross-validation framework, improved. The first observation makes the results appear more plausible. The second observation supports the notion that the ICA improves the relative signal-to-noise ratio. Finally, the third observation provides objective evidence for the benefit of this approach, such that this novel pre-processing step improves our ability to predict behavioral status. A common challenge with many resting-state analyses is determining whether a particular pre-processing procedure improves or hinders performance, as we do not know if detecting more correlations (or more anti-correlations) is better or worse, or whether simulated data accurately captures the noise observed in real data (Weissenbacher et al., 2009). However, using the resting state data to predict individual variability on independent out-of-sample behavioral measures provides an objective measure to assess the influence of different pre-processing steps. Therefore, we suggest that ICA is effective in identifying and isolating the artifacts related to the lesion, and advocate the removal of the corresponding independent components in preprocessing of fMRI for functional connectivity studies.

Methods

Participants

Participants were recruited from the local community as part of a larger study of aphasic impairment associated with left hemisphere stroke. The research was approved by the Institutional Review Board at the University of South Carolina. Only individuals with aphasia resulting from a single ischemic or hemorrhagic stroke involving the left hemisphere were included. Participants with lacunar infarcts or with damage that only involved the brainstem or cerebellum were excluded. The behavioral assessment of the participants took place between May 2007 and March 2015, and 74 individuals were included in the data analyses. The mean ± SD age was 60.3 ± 9.4 years (range, 38–81 years), and 23 were women. All participants were at least 6 months post-stroke, and the mean ± SD time since stroke onset was 38.6 ± 48.7 months (range, 6–276 months).

Behavioral evaluation

Aphasic impairment was assessed using the Western Aphasia Battery (WAB; Kertesz, 1982). For this study, we used six composite measures of impairment:

information content of the participant's speech;

speech fluency;

auditory comprehension;

speech repetition;

oral naming;

aphasia quotient, which is an overall score of aphasia severity.

Data acquisition

MRI scanning was performed within two days of behavioral testing of language abilities. Images were acquired on a Siemens Trio 3 T scanner equipped with a 12-element head coil located at the University of South Carolina. Three images relevant for this work were acquired for each participant:

T1-weighted image utilizing an MP-RAGE sequence with 1 mm isotropic voxels, a 256 × 256 matrix size, and a 9-degree flip angle. For 69 individuals we used a 160 slice sequence with TR = 2250 ms, TI = 900 ms, TE = 4.52 ms. For the remaining 6 individuals we used a 192 slice sequence with TR = 2250 ms, TI = 925 ms, TE = 4.15 with parallel imaging (GRAPPA = 2, 80 reference lines). Each of these scans required approximately 7 min to acquire.

T2-weighted image using a sampling perfection with application optimized contrasts using a different flip angle evolution (3D-SPACE) sequence. This 3D TSE scan uses a TR = 2800 ms, a TE of 402 ms, variable flip angle, 256 × 256 matrix scan with 192 slices (1 mm thick), using parallel imaging (GRAPPA =, 120 reference lines). This used the same slice center and angulation as the T1 scan.

Resting state fMRI: For 59 participants we acquired an EPI sequence with a 208 × 208 mm field of view, a 64 × 64 matrix size, and a 75-degree flip angle, 34 axial slices (3 mm thick with 20% gap yielding 3.6 mm between slice centers), TR = 1850 ms, TE = 30 ms, GRAPPA = 2, 32 reference lines, sequential descending acquisition, 196 volumes acquired. For the remaining 16 participants were acquired using a multiband sequence (× 2) with a 216 × 216 mm field of view, a 90 × 90 matrix size, and a 72-degree flip angle, 50 axial slices (2 mm thick with 20% gap yielding 2.4 mm between slice centers), TR = 1650 ms, TE = 35 ms, GRAPPA = 2, 44 reference lines, interleaved ascending slice order, 370 volumes acquired. One functional run was acquired (6 min for the single-band sequence, and 10 min for the multi-band sequence). The participants were instructed to stay still during the acquisition.

Preprocessing of structural scans

Lesions were manually drawn on the T2 weighted image by a neurologist, who was blinded to the participant's language scores at the time of the lesion drawing. The T2 image was co-registered to the T1 image, and these parameters were used to reslice the lesion into the native T1 space. The resliced lesion maps were smoothed with a 3 mm full-width half maximum (FWHM) Gaussian kernel to remove jagged edges associated with manual drawing. We then performed enantiomorphic normalization (Nachev et al., 2008) using SPM12 and Matlab scripts we developed (Rorden et al., 2012) as follows: First, a mirrored image of the T1 scan (reflected around the midline) was created, and this mirrored image was coregistered to the native T1 image. We then created a chimeric image based on the native T1 scan with the lesioned tissue replaced by tissue from the mirrored scan (using the smoothed lesion map to modulate this blending, feathering the lesion edge). SPM12’s unified segmentation-normalization (Ashburner and Friston, 2005) was used to warp this chimeric image to standard space, with the resulting spatial transform applied to the actual T1 scan as well as the lesion map. The normalized lesion map was then binarized, using a 50% probability threshold. Fig. 1 shows the overlap of lesions for our sample of participants.

Fig. 1

Overlap of lesions in the 74 participants form our study. The color scale shows the number of individuals with damage at each voxel. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

To assess regional brain damage, we used the brain parcellation developed by Joliot et al. (2015), which segments the grey matter into 384 regions of interest (ROIs). The anatomical brain atlas containing the parcellation was aligned with each individual's T1-weighted image. The T1-weighted image was segmented into probabilistic grey and white matter maps, and the grey matter map was divided into regions according to the atlas. Then, regional damage was computed as the proportion of intact (i.e. not lesioned) voxels per each ROI.

Preprocessing of fMRI

Functional MRI data were preprocessed using four different pipelines, each of them having a different approach to removal of lesion-driven artifacts. All pipelines started with the same initial steps, commonly used to preprocess fMRI data for neurologically intact participants as well as for clinical populations. A description of this initial preprocessing pipeline is provided below, followed by the description of the procedures specific to the four pipelines.

Initial procedure

The fMRI data were corrected for motion using the SPM12 “realign and unwarp” procedure with default settings. After that, we performed brain extraction using the SPM12 script pm_brain_mask with default settings. Slice time correction was also done using SPM12 (for the 16 participants acquired with a multiband sequence, this step was skipped). After that, the mean fMRI volume for each participant was aligned to the corresponding T2-weighted image to compute the spatial transformation between the fMRI data and the lesion mask. The fMRI data were then spatially smoothed with a Gaussian kernel with FWHM = 6 mm.

Pipeline A

This pipeline did not perform any specific steps to remove the lesion-driven artifacts. The voxel-wise fMRI timecourses were detrended using the following regressors: mean signal from the white matter, obtained from the chimeric T1-weighted image; time courses of the six motion parameters estimated at the motion correction step; linear, quadratic and cubic trends. Then, the timecourses were bandpass-filtered using the 0.01–0.1 Hz frequency band. Global mean signal has not been regressed out.

Pipeline B

In this pipeline, removal of the lesion-driven artifacts was performed by regressing out the mean signal from the CSF maps (which included the site of the lesion) obtained from the T1-weighted images. This signal was added to the list of regressors used in the detrending as described in Pipeline A.

Pipeline C

This pipeline used Independent Component Analysis (ICA) to identify and remove the lesion-driven artifacts. First, the fMRI data were processed with Pipeline A. Then, we applied FSL MELODIC package to compute the Z-scored spatial maps of the independent components. These spatial maps were thresholded at p < 0.05, and compared with the lesion mask for that participant. Since both the lesion mask and the thresholded IC map were binary images, we were able to utilize the Jaccard index (the number of voxels in the intersection divided by the number of voxels in the union) to quantify the amount of spatial overlap. If the Jaccard index was > 5%, the corresponding component was deemed to be significantly overlapping with the lesion mask (the threshold of 5% was chosen to resemble the threshold of 0.05 that is frequently used for testing the significance of a p value). All such components were then regressed out of the fMRI data using the fsl_regfilt script from the FSL package. The matlab script that identifies the lesion-driven independent components was written in MATLAB using both SPM12 and FSL packages. The script, called nii_filter_lesion_ICs.m, can be downloaded from https://github.com/neurolabusc/nii_preprocess.

Pipeline D

This pipeline excluded the lesion voxels from the connectome computation. For each ROI, the voxels that were included in the lesion map were masked out and did not contribute to the mean fMRI timecourse for that ROI. If the ROI was completely covered by lesion, its connectivity with other ROIs was treated as a missing value.

Construction of functional connectomes

After processing the data with the four pipelines described above, we computed the functional connectomes. Functional connectivity (FC) was computed as Pearson's correlation coefficient between the mean timecourses of each of the 384 grey-matter ROIs in our parcellation. Therefore, we obtained four 384 × 384 matrices for each participant, corresponding to four preprocessing pipelines. The connectomes created with Pipeline D had missing values corresponding to the ROIs completely covered by the lesion.

Univariate analysis of connectome-behavior relationships

To investigate the relationship between functional connectivity and behavioral impairment, we calculated the relationship between functional connectomes and WAB scores. For each pair of ROIs and each WAB score, we computed the General Linear Model with the score as the outcome and the functional connectivity as the predictor. The t values for the predictor were Z-transformed using SPM's spm_t2z function. The connectomes created with Pipeline D contained missing values, therefore, for this particular pipeline, the number of participants (and the associated degrees of freedom) varied across functional connections. Because each functional connection was analyzed independently of the others, this is an example of univariate analysis. The analysis was limited to the cortical ROIs in the left hemisphere (basal ganglia and thalamus were excluded because the ventricles were greatly enlarged in participants with big lesions, which presented severe difficulties in across-subject alignment of subcortical ROIs); there are 173 cortical grey-matter ROIs in the left hemisphere. Bonferroni correction was used to correct for multiple comparisons to ensure that the same multiple comparison-corrected threshold was applied in every analysis. The surviving functional connections were visualized using the Surfice software (https://www.nitrc.org/projects/surfice).

Multivariate analysis of connectome-behavior relationships

As a complementary analysis, we performed support vector regression (SVR) within a leave-one-out cross-validation framework to predict the WAB score from the functional connectome. SVR was performed using a linear kernel and a C value of 0.01; scripts from the LIBSVM library (Chang and Lin, 2011) were used to train and test the SVR model. As in univariate analyses, the inputs to the SVR consisted of functional connections between left-hemisphere cortical regions. This framework has been described in our previous publication (Yourganov et al., 2016). At each iteration, one participant was left out, and the SVR model was trained on the remaining 73 participants; this model was then used to estimate the WAB score for the left-out participant. The correlation between actual and predicted WAB scores served as a measure of accuracy of the SVR prediction. Both univariate and multivariate analyses were performed using our NiiStat package that can be downloaded from https://github.com/neurolabusc/NiiStat.

Results

Identification of lesion-driven independent components

Many participants had substantial lesions covering multiple ROIs. Fig. 2 shows the distribution of damaged regions (where lesion covers at least 95% of the ROI volume) in our sample of participants. For the majority of our participants (54 out of 75), Pipeline C identified at least one lesion-driven independent component. The number of such components was strongly dependent on the size of the lesion: participants with large lesions had a larger number of lesion-driven independent components. Fig. 3 shows the scatter plot (each circle representing a participant) that shows this dependency. Here, the lesion size was computed on the brains that were spatially normalized to the standard template.

Fig. 2

Distribution of damaged regions (where at least 95% of the region is covered by lesion) in our participant sample. The total number of ROIs is 173.

Fig. 3

The relation between lesion volume and the number of independent components that represent lesion-driven artifacts. Each circle represents a participant.

Correlations of BOLD within the lesion

Fig. 4 illustrates the importance of removing the lesion-driven artifacts from the fMRI data. We inspected the lesions for each participant, identified the ROIs that were at least 95% covered by a lesion, and inspected the functional connectivity between these ROIs. Fig. 4 shows the histogram of these correlations (pooled across participants). If the lesion artifacts are not removed (pipeline A), we see a large number of high FC values (the histogram peaks at 0.6, and the median FC value is 0.5). The reason for such high FC values between the two areas that are lesioned is that the BOLD signal at these two areas is sampling the physiological noise from two parts of the lesion, which is very likely to be highly correlated. Regression of mean CSF signal (pipeline B) decreases the FC between lesioned areas, bringing the median down to 0.28 and the histogram peak down to 0.5. However, it also introduces correlations that are strongly negative (4% of the observed correlations are below −0.4, compared to 1.6% in pipeline A). Pipeline C, where the lesion-driven artifacts are identified with ICA, is more effective than Pipeline B in bringing down the FC values (median FC: 0.21; histogram peak: 0). In addition, the distribution of negative FC values in Pipeline C (1.9% of FC values are below − 0.4) is very similar to Pipeline A. We cannot conclude that Pipeline C does not introduce strongly negative correlations, but it doesn't change their distribution (i.e., introduction of strongly negative correlations is balanced by removal of other strongly negative correlations). Also, the inter-quartile range (IQR) of the FC values was smaller in Pipeline C compared to pipelines A and B (Pipeline A: IQR = 0.51; Pipeline B: IQR = 0.54; Pipeline C: IQR = 0.47). In the connectomes created with Pipeline D, the connections stemming from fully lesioned ROIs are not computed; for this pipeline, Fig. 4 shows the distribution of FC values for the ROIs with damage from 95% up to (but not including) 100%; the histogram of these values has a sharp peak at zero.

Fig. 4

Histogram showing the distribution of functional connectivity between pairs of lesioned areas, as computed using four separate pipelines (Pipeline A: no removal of lesion-driven artifacts; Pipeline B: regression of mean CSF signal; Pipeline C: removal of independent components associated with the lesion; Pipeline D: masking out lesioned voxels).

We analyzed the correlations between the speech and language scores and the functional connectomes created by the four pipelines. Fig. 4, Fig. 5 display the connections that relate significantly to each of the six behavioral scores (and survive the Bonferroni correction for multiple comparisons). Fig. S1, Fig. S2 show the corresponding results obtained for the 59 participants scanned with a single-band sequence. The connections displayed in yellow have a positive association with the score, that is, increased functional connectivity is associated with better behavioral performance. The connections displayed in green have a negative association: increased functional connectivity is associated with worse behavioral performance. Table 1 lists the number of positively and negatively associated connections for each pipeline and behavioral score.

Fig. 5

Functional connections that were significantly correlated with behavioral scores (information content, speech fluency, and auditory comprehension). Yellow lines correspond to positive associations and green lines correspond to paradoxical negative associations (e.g. where increased connectivity is correlated with poorer behavioral (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)performance).

Fig. S1

Functional connections that were significantly correlated with behavioral scores (information content, speech fluency, and auditory comprehension), for the 59 participants scanned with a single-band sequence. Yellow lines correspond to positive associations and green lines correspond to paradoxical negative associations (e.g. where increased connectivity is correlated with poorer behavioral performance).

Fig. S2

Functional connections that were significantly correlated with behavioral scores (repetition, naming, and aphasia quotient), for the 59 participants scanned with a single-band sequence. Yellow lines correspond to positive associations and green lines correspond to paradoxical negative associations (e.g. where increased connectivity is correlated with poorer behavioral performance).

Table 1

Number of functional connections significantly correlated with behavioral scores (either positively or negatively).

Behavioral score	Pipeline A		Pipeline B		Pipeline C		Pipeline D
	Positive	Negative	Positive	Negative	Positive	Negative	Positive	Negative
Information content	3	0	0	4	12	0	7	0
Speech fluency	29	2	0	8	34	0	10	0
Auditory comprehension	2	3	1	1	12	0	4	0
Repetition	6	0	0	0	14	0	2	0
Naming	7	3	4	4	23	0	9	0
Aphasia quotient	13	3	2	4	22	0	9	0

Pipeline A (top row on Fig. 5, Fig. 6), where the lesion-driven artifacts are not removed, identified some positively associated functional connections for each behavioral score; most of them involved regions well-known for their importance in processing speech and language, such as inferior frontal, posterior temporal, inferior parietal, and insular regions. However, we revealed some negative associations as well (for speech fluency, auditory comprehension, and naming scores, as well as for the aphasia quotient). These negative associations were in the inferior parietal, posterior temporal, and insular parts of the brain, and these negative associations are paradoxical: why would higher functional connectivity between them be associated with impaired performance? The answer is explained by Fig. 4, where we see that Pipeline A produces high positive correlations between lesioned areas. All the areas involving negative associations are perfused by the middle cerebral artery and are frequently lesioned together in participants with large lesions; therefore, the functional correlation between them could be high whereas performance could be severely impaired. In people with smaller lesions (and more intact behavior), some of these areas might be preserved, and a correlation between a preserved area (reflecting neural activity) and a lesioned area (reflecting physiological noise) might be quite low. Therefore, across participants, Pearson's correlation between these functional connections and behavioral score would manifest as a negative association.

Fig. 6

Functional connections that were significantly correlated (yellow) or anticorrelated (green) with behavioral scores (repetition, naming, and aphasia quotient). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Pipeline B (second row on Fig. 5, Fig. 6), instead of solving this problem of negative associations, greatly exacerbates it. As can be seen in these figures and in Table 1, the majority of significant connections is negatively associated with behavioral scores. For information content and fluency scores, no positive associations have been found. On the other hand, Pipeline C (third row on Fig. 5, Fig. 6) gets rid of negative associations completely. Moreover, relative to Pipeline A, a much larger number of positive associations is identified for every score. Pipeline D (bottom row of Fig. 5, Fig. 6) also gets rid of negative associations, but the number of positive associations is markedly less than the number observed for Pipeline C (Table 1). Importantly, Fig. 5, Fig. 6 show that these positive associations are largely driven by lesion frequency: most of them stem from the insula, which is the most frequently lesioned part of the brain in our sample (Fig. 1), and tend to be short-range (i.e. terminate in the vicinity of the insula, which is likely to be lesioned as well).

To see the effect of frank damage on the relationship between functional connectivity and behavioral impairment, we looked closely at the connectivity between two ROIs which are highly involved in speech and language (see e.g. Fridriksson et al., 2016): pars triangularis of the inferior temporal gyrus (which is part of Broca's area), and the posterior part of the superior temporal gyrus (also known as area Spt, and is part of Wernicke's area). Functional connectivity between these two ROIs was found to be significantly correlated with all behavioral scores, regardless of the pipeline used to create the connectomes. Table 2 shows these correlations (as Pearson's correlation coefficient; p < 0.05 for all reported correlations); the corresponding scatter plots are provided in the Supplementary material. One can see that the strongest correlations are consistently observed in the connectomes created with Pipeline C, and the weakest with Pipeline D (to ensure the fairness of comparison across pipelines, we excluded the four subjects who had complete damage in at least one of the two ROIs, and, therefore, the functional connectivity was not computed by Pipeline D). Moreover, these correlations are significant even if the damage to these two regions is controlled for: Table 3 displays the corresponding partial correlations, and, again, Pipeline C yields the strongest correlations, which are significant even after controlling for the summary damage to the two ROIs, and Pipeline D yields the weakest correlations.

Table 2

Correlation between the functional connectivity of the two illustrative ROIs (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and behavioral scores.

	Pipeline A	Pipeline B	Pipeline C	Pipeline D
Information content	0.459	0.466	0.625	0.322
Speech fluency	0.435	0.468	0.592	0.377
Auditory comprehension	0.307	0.296	0.550	0.322
Repetition	0.453	0.451	0.581	0.381
Naming	0.409	0.384	0.615	0.309
Aphasia quotient	0.439	0.445	0.626	0.375

Table 3

Partial correlation between the functional connectivity of the two illustrative ROIs (see Table 2) and behavioral scores, controlling for the total amount of damage in the two ROIs. Asterisks denote the correlation coefficients with p < 0.05.

	Pipeline A	Pipeline B	Pipeline C	Pipeline D
Information content	0.340*	0.368*	0.441*	0.169
Speech fluency	0.302*	0.368*	0.375*	0.237*
Auditory comprehension	0.107	0.116	0.289*	0.150
Repetition	0.329*	0.345*	0.365*	0.246*
Naming	0.223	0.190	0.346*	0.108
Aphasia quotient	0.305*	0.336*	0.417*	0.228

Univariate analysis identifies the functional connections that are, by themselves, correlated with behavioral measures. Multivariate analysis, on the other hand, uses all functional connections to predict the behavioral scores. To compare the predictive power of the four pipelines, we predicted the six scores within a leave-one-participant-out procedure and computed Pearson's correlation coefficients between actual and predicted scores. These correlations are listed in Table 4; pipeline C consistently yields more accurate predictions, as indicated by higher correlation coefficients, than the other pipelines. It should be noted, however, that the connectomes created with Pipeline D contain missing values and therefore present a problem for multivariate analysis such as SVR; LIBSVM discards the functional connections which have any missing values, which results in poor prediction evident in Table 4.

Table 4

Correlation between actual and predicted behavioral scores (as estimated by multivariate support vector regression).

	Pipeline A	Pipeline B	Pipeline C	Pipeline D
Information content	0.426	0.542	0.563	0.223
Speech fluency	0.615	0.691	0.700	0.435
Auditory comprehension	0.570	0.519	0.652	0.455
Repetition	0.505	0.535	0.552	0.352
Naming	0.534	0.449	0.627	0.394
Aphasia quotient	0.567	0.644	0.671	0.418
Mean	0.536	0.563	0.628	0.380

The Z scores that correspond to the correlations (according to Fisher's r-to-Z transform) are not significantly different across pipelines A, B, and C for any given score (p > 0.06). Pipeline D, however, yields significantly worse actual-versus-predicted correlations than pipelines B and C. The Z-transformed correlation coefficients are significantly different (p < 0.05) between pipelines B and D for the prediction of Information Content, Fluency, and Aphasia Quotient scores. Between pipelines C and D, there is a significant difference for five out of six scores (p = 0.066 for the repetition score, and p < 0.05 for the remaining five scores). The differences between pipelines A and D are not significant (p > 0.06).

Discussion

We have demonstrated that functional connectivity analysis in stroke survivors has some particular challenges due to the presence of the lesion. The image signal sampled from two locations within a lesion is likely to be highly correlated, because it is sampling physiological noise from the same pool of necrotic tissue. Therefore, participants with large lesions are likely to have high values of functional connectivity between lesioned areas. These participants are also likely to have significant behavioral impairment, which leads to paradoxical results when poor behavioral performance is associated with seemingly high functional connectivity (“negative associations”). Regression of mean CSF signal, rather than solving this problem, further exacerbates it: it increases the number of functional links negatively associated with the behavioral scores, and reduces the number of positively associated functional links. This might be related to the theoretical result, reported by Murphy et al. (2009), that regression of global mean signal introduces negative correlations in the data. Here, we observe a similar effect in the lesioned areas (which consist mostly of CSF), where regression of mean CSF signal introduces negative correlations.

We propose a solution to this problem by employing Independent Component Analysis to isolate the sources of variance in fMRI data that are driven by lesion, which can be then easily filtered out using the scripts from the FSL MELODIC package. The superiority of this approach is demonstrated in several ways:

This approach reduces the correlations within lesioned areas and does not introduce negative correlations, in effect getting rid of the negative associations between functional connections and behavioral scores;

It also improves sensitivity of the univariate analysis of FC-behavior relationship, increasing the number of functional connections that are correlated with improvement in behavior;

It improves the predictive power of the functional connectome by increasing the accuracy of multivariate prediction of behavioral scores.

An alternative approach, where the lesion voxels are masked out during the computation of the connectome, was found inferior to the ICA-based lesion-artefact removal. Masking leads to imbalanced number of voxels (and, therefore, a difference in statistical power) in the same ROIs across subjects, which might be an additional source of noise in already noisy data. Also, the functional connections involving the completely lesioned ROIs are treated as missing values, which introduces a further difference in statistical power across the pairs of ROIs. These missing values present a particular problem for multivariate analysis. On the other hand, univariate analysis computed on these connectomes is heavily influenced by lesion location: the functional connections that survive the multiple comparison correction are located in the insular region and its immediate neighborhood, which is the most frequent locus of lesions in our sample (as is common in middle cerebral artery occlusions).

Our study only included participants at the chronic stage of post-stroke recovery. However, our preprocessing pipeline could also be used to process data from acute patients. The only change would be the creation of lesion masks, which would utilize diffusion-weighted images rather than T1-weighted images as typically done in chronic participants.

Currently, the most stable finding of functional connectivity analysis in stroke studies is the association between behavioral impairment and a decrease in interhemispheric connectivity, in particular between homotopic regions (for various types of post-stroke behavioral impairment, this was reported by Carter et al., 2010, Xu et al., 2014, New et al., 2015, Bannister et al., 2015, Siegel et al., 2016, Tang et al., 2016, Adhikari et al., 2017, Sandberg, 2017, Yang et al., 2017; see also review by Grefkes and Fink, 2011). It is possible that this significant and stable association might be explained by frank lesion effect, rather than by change in neuronal synchrony within the relatively intact regions of a network. In the healthy brain, the strongest correlations are often found between homotopic regions (Mišić et al., 2014). Cerebral infarctions are often limited to one hemisphere, so if the brain area that is heavily involved in a particular behavioral function is lesioned, the connectivity with its contralesional homologue (and, generally, with the remaining network in the contralesional hemisphere) decreases. This association between inter-hemispheric connectivity and behavioral impairment might be so strong that it could obscure the relationship between behavioral impairment and connectivity within a single hemisphere. We speculate that our suggested approach to lesion-artefact removal, by increasing the predictive power of the functional connectome, could clarify the relationship between the FC within the ipsilesional hemisphere and post-stroke impairment. This is particularly relevant to behavioral functions largely involving brain areas within a hemisphere, such as different aspects of speech and language processing.

The following are the supplementary data related to this article.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the aphasia quotient. The functional connectome was obtained with pipeline A.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the auditory comprehension scores. The functional connectome was obtained with pipeline A.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the fluency scores. The functional connectome was obtained with pipeline A.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the information content scores. The functional connectome was obtained with pipeline A.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the naming scores. The functional connectome was obtained with pipeline A.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the repetition scores. The functional connectome was obtained with pipeline A.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the aphasia quotient scores. The functional connectome was obtained with pipeline B.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the auditory comprehension scores. The functional connectome was obtained with pipeline B.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the fluency scores. The functional connectome was obtained with pipeline B.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the information content scores. The functional connectome was obtained with pipeline B.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the auditory comprehension scores. The functional connectome was obtained with pipeline B.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the repetition scores. The functional connectome was obtained with pipeline B.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the aphasia quotient scores. The functional connectome was obtained with pipeline C.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the auditory comprehension scores. The functional connectome was obtained with pipeline C.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the fluency scores. The functional connectome was obtained with pipeline C.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the information content scores. The functional connectome was obtained with pipeline C.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the naming scores. The functional connectome was obtained with pipeline C.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the repetition scores. The functional connectome was obtained with pipeline C.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the aphasia quotient scores. The functional connectome was obtained with pipeline D.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the auditory comprehension scores. The functional connectome was obtained with pipeline D.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the fluency scores. The functional connectome was obtained with pipeline D.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the information content scores. The functional connectome was obtained with pipeline D.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the naming scores. The functional connectome was obtained with pipeline D.

Scatter plot of correlations between two language areas (pars triangularis of the inferior frontal gyrus, and posterior part of the superior temporal gyrus) and the repetition scores. The functional connectome was obtained with pipeline D.