Voxel-Wise Motion Artifacts in Population-Level Whole-Brain Connectivity Analysis of Resting-State fMRI

Abstract

Functional Magnetic Resonance Imaging (fMRI) based brain connectivity analysis maps the functional networks of the brain by estimating the degree of synchronous neuronal activity between brain regions. Recent studies have demonstrated that “resting-state” fMRI-based brain connectivity conclusions may be erroneous when motion artifacts have a differential effect on fMRI BOLD signals for between group comparisons. A potential explanation could be that in-scanner displacement, due to rotational components, is not spatially constant in the whole brain. However, this localized nature of motion artifacts is poorly understood and is rarely considered in brain connectivity studies. In this study, we initially demonstrate the local correspondence between head displacement and the changes in the resting-state fMRI BOLD signal. Than, we investigate how connectivity strength is affected by the population-level variation in the spatial pattern of regional displacement. We introduce Regional Displacement Interaction (RDI), a new covariate parameter set for second-level connectivity analysis and demonstrate its effectiveness in reducing motion related confounds in comparisons of groups with different voxel-vise displacement pattern and preprocessed using various nuisance regression methods. The effect of using RDI as second-level covariate is than demonstrated in autism-related group comparisons. The relationship between the proposed method and some of the prevailing subject-level nuisance regression techniques is evaluated. Our results show that, depending on experimental design, treating in-scanner head motion as a global confound may not be appropriate. The degree of displacement is highly variable among various brain regions, both within and between subjects. These regional differences bias correlation-based measures of brain connectivity. The inclusion of the proposed second-level covariate into the analysis successfully reduces artifactual motion-related group differences and preserves real neuronal differences, as demonstrated by the autism-related comparisons.

About the Authors:

Tamás Spisák

* E-mail: [email protected]

Affiliation: Department of Nuclear Medicine, Medical and Health Science Centre, University of Debrecen, Debrecen, Hungary

András Jakab

Affiliation: Department of Radiology, Medical University of Vienna, Vienna, Austria

Sándor A. Kis

Affiliation: Department of Nuclear Medicine, Medical and Health Science Centre, University of Debrecen, Debrecen, Hungary

Gábor Opposits

Affiliation: Department of Nuclear Medicine, Medical and Health Science Centre, University of Debrecen, Debrecen, Hungary

Csaba Aranyi

Affiliation: Department of Nuclear Medicine, Medical and Health Science Centre, University of Debrecen, Debrecen, Hungary

Ervin Berényi

Affiliation: Department of Biomedical Laboratory and Imaging Science, Medical and Health Science Center, University of Debrecen, Debrecen, Hungary

Miklós Emri

Affiliation: Department of Nuclear Medicine, Medical and Health Science Centre, University of Debrecen, Debrecen, Hungary

Introduction

Background

Motion artifacts are problematic for all types of MRI including resting-state functional MRI (fMRI) and therefore motion correction is a vital step in every work-flow during fMRI analysis. According to population-level analysis and group comparisons, retrospective motion related artifact removal strategies can be performed at five different stages of the data processing pipeline: (i) motion correction of fMRI time-series by realignment to a reference image (using automatic co-registration approaches) [1]; (ii) censoring data to exclude periods of high motion (scrubbing, de-spiking) [2], [3]; (iii) modeling the effect of motion-related nuisance parameters on blood oxygen level dependent (BOLD) signal [4]–[7]; (iv) temporal filtering of BOLD timecourses to discard frequencies encumbered by motion artifacts and (v) correct for subject-specific motion effects on population-level (descriptive summary statistics of subject-specific motion as second-level model regressors) [8]–[10].

Traditional realignment-based correction approaches ensure that different time-points of the BOLD signal correspond to the same location. However, such methods do not handle voxel-level intensity confounds originating from the establishment of magnetic gradients and subsequent readout of the BOLD signal [2], [11]. Furthermore, automatic co-registration approaches may introduce spurious displacements in intervals of relatively low motion [12]. Nonetheless, subject movement is often measured with parameters based upon the resulting image realignment transformations.

Large BOLD intensity confounds (spikes) in time-frames with extreme, abrupt movement can be eliminated from the analysis by simply dropping the corrupted data (“scrubbing”) [2] or by spike-regression [13]. However, the reduction in time points is associated with an increase in the likelihood of high correlation scores [8]; moreover, recent findings [8], [14], [15] suggest that in population-level functional connectivity studies, scrubbing can be omitted from the analysis when using proper second-level correction.

Intensity confounds originating from micro-movements (as small as 0.1 mm from one time point to the next) can also disrupt results, especially in case of correlation based connectivity analysis methods, where such small but temporally concordant noise leads to spurious increase in connectivity strength [2], [9], [16]. Nuisance signal regression approaches aim to eliminate the signal components of non-neuronal origin from the raw BOLD data utilizing linear regression. These confounder signals can be defined by dedicated physiological monitoring devices during the scan, calculated from motion parameters extracted during spatial realignment based motion correction or derived directly form the data itself, using a “noise ROI” [6], [7], [17].

Temporal filtering is also a crucial step in the reduction of fMRI brain connectivity artifacts, including motion confounds. Most connectivity studies apply a band-pass filter with a high-pass cutoff of 0.008–0.01 Hz and a low-pass threshold of 0.08–0.1 Hz [18]. While there is evidence that resting-state networks are present at a relatively broad band in the frequency spectra [19], slight modifications in the frequency band have been suggested [3].

Some of these techniques can effectively reduce not only motion-related effects, but also physiological noise (e.g. cardiac or respiratory confounds) or hardware drifts and instabilities. However, recent studies [3], [8], [20] report that clear artifacts remain in the data even after such regression and filtering approaches, and that these artifacts have systematic effects upon resting-state functional MRI connectivity (rs-fcMRI) patterns.

When performing group comparisons in functional connectivity studies, one can account for these motion-related residual artifacts during the second-level analysis by inclusion of motion-related, subject-specific covariates into the population-level model. A common choice is to include a measure of the average patient movement [8]–[10]. Alternatively, the value of global voxel-to-voxel correlation (GCOR) can be utilized as confounder, as in [20], although the latter quantity can also carry valuable neurological information.

Due to in-scanner head rotation, the effect of patient movements on the BOLD signal is not spatially constant in the whole brain; however, this local relationship is poorly understood and is rarely considered in brain connectivity studies. The pattern of voxel-wise motion not only varies among different loci of the same subject, but also among subjects. According to Satterthwaite et al. [9], between-subject differences of motion are stable and hence, in-scanner head motion should be considered as a trait. Thus, the effect of the subject-specific spatio-temporal motion pattern on the BOLD signal could bias group analysis when different groups have different tendencies in their spatio-temporal motion patterns. This is particularly problematic in studies when regional connectivity deficits are associated with a pathological condition, and thus, limits the usability of functional connectivity as a biomarker of disease. These biases in the functional connectivity pattern can lead to invalid conclusions regarding biomedical hypotheses, as denoted by [21] and demonstrated by [20], especially in patholological conditions associated with hyperkinetic patients (epilepsy, attention deficit hyperactiviy disorder, some forms of autism). Group-wise inconsistencies in motion patterns can arise from different patient positioning and multi-center studies are also challenging in this regard.

Purposes

Here, we hypothesize that at least some of the above-mentioned artifactual effects may originate from the complex voxel-wise spatio-temporal nature of head displacement, and can be modeled more efficiently using this information.

The possibilities for utilizing the voxel-wise nature of in-scanner motion in artifact removal approaches has not been intensively investigated, as yet. As recently reported by two studies [3], [14] and confirmed by our preliminary analysis, including voxel-wise displacement parameters as voxel- or region-wise covariates in the appropriate nuisance signal regression model does not significantly improve motion artifact removal, compared to the usual technique, the regression of spatially averaged global displacement parameters. However, Yan et al. in [14] also brings up the possibility that an appropriate correction technique may have greater success in using the rich information encapsulated by voxel-specific indices.

Our study was designed to characterize the impact of voxel-wise head motion artifacts in population-level rs-fcMRI brain connectivity studies and investigate how this local information on displacement can be utilized for artifact removal.

We initially demonstrate the local correspondence between head displacement and the changes in the rs-fcMRI BOLD signal. We then aim to investigate how functional connectivity strength is affected by the deviations in the average regional spatial displacements on the population-level. We propose Regional Displacement Interaction (RDI), a novel modeling approach for second-level brain connectivity analysis, which provides the opportunity to incorporate voxel-wise motion information into the population-level model and to account for corresponding artifactual effects. The effectiveness of this motion artifact reduction technique is evaluated by investigating the variance explained by the proposed confound modeling covariates in the model. The method is than demonstrated in group comparisons of cohorts with differing average voxel-wise displacement patterns. Due to the disagreement [3], [14], [20], [22] about the optimal first-level nuisance signal regression technique, we perform a comparison of prevailing first-level nuisance signal regression approaches and characterize their interference with the proposed method. Finally, to test whether the proposed method preserves group differences of neuronal origin, a comparison of autistic and control groups is performed.

Materials and Methods

Image acquisition

Analysis was performed on the resting-state fMRI data of 184 patients obtained from the Autism Brain Imaging Dataset Exchange database [15], [23], [24] (ABIDE). All of the images were acquired at the NYU Langone Medical Center using a 3 Tesla Siemens Magnetom Allegra syngo MR 2004A. A T1-weighted sagittal MP-RAGE structural image was obtained (TE = 3.25 ms, TR = 2530 ms, TI = 1100 ms, flip angle = 7, 256 slices with 1.3×1×1.3 mm voxels). Functional images were obtained using a BOLD contrast sensitive gradient echo echo-planar sequence (TE = 15 ms, flip angle = 90, in-plane resolution = 3×3 mm; volume TR = 2000 ms). Whole-brain coverage for the functional data was obtained using 33 contiguous interleaved 4 mm axial slices. During the resting-state fMRI scan, most participants were asked to relax with their eyes open, while a white cross-hair against a black background was projected on a screen. However, data were also included for some individuals who were asked to keep their eyes closed; and, in a few cases, participants closed their eyes regardless of instructions to keep them open.

The population sample consisted of 79 patients with autism spectrum disorders (7.1–39.1 years) (53 Autistic Disorder, 21 Asperger's Disorder, 5 Pervasive Developmental Disorder-Not Otherwise Specified) and an age and gender-matched group of 105 typical control subjects (6.5–31.8 years).

Data collection for the ABIDE dataset was approved by the institutional review boards of the New York University School of Medicine and New York University. Prior to participation, written informed consent and assent (for participants>18 years) were obtained from all participants and their parents/legal guardians (for participants<18 years). Participants received monetary compensation for completing the study. In this study, the patient data were analyzed anonymously.

Preprocessing

FMRI time series were co-registered and frame-wise estimation of head displacement was performed using FSL McFlirt [25], [26]. Matrices defining the rigid-body (three translation and three rotation parameter) transformation that fit each frame to the reference frame (at the middle time-point) were saved for further use. The first five volumes of each dataset were discarded from further analysis to allow for T1 equilibration effects. BET was used to remove non-brain areas [27]. The resulting pre-processed fMRI data were nonlinearly co-registered to the brain-extracted anatomical image, and then, spatially standardized to the MNI152 space using the FLIRT and FNIRT utilities [28] of the FSL package, to achieve spatial correspondences for group analysis. Since further processing steps utilized averaged regional time courses, no smoothing was applied on the images.

Calculation of voxel-wise displacement

With an in-house-developed utility based on the m3i software library system [29], transformation matrices outputted by McFlirt were converted to world coordinate origin. The respective inverse transformations were applied to each frame of the fMRI time-series and the root mean squared voxel position change in world coordinates was calculated for each voxel of each frame. The first derivate of the resulting local displacement time-series was saved in NIfTI format dynamic images in the same space as the fMRI time-series (see Fig. 1 for demonstrative images), and then realigned to standard space.

[Figure omitted. See PDF.]

Figure 1. Voxel-wise characteristics of head motion during an fMRI scan.

Examples of different patterns of voxel-wise displacement within the time frames of one illustrative subject (A) and temporally averaged voxel-wise displacement of six illustrative subjects (B).

https://doi.org/10.1371/journal.pone.0104947.g001

ROI definition

In order to improve the signal-to-noise ratio and reduce the amount of data to analyze, all regional timecourses (regional BOLD signal, temporal derivate of its root mean squared variance, regional displacement) and corresponding correlation coefficients presented in this paper were drawn from a set of ROIs (M = 88) that were defined based on the Harvard-Oxford Cortical and Subcortical brain atlases [30]. Probability maps for all regions were accessed and region borders were delineated by retaining voxels with a probability greater than 25%. Voxels associated with multiple regions (in case of overlapping regions) were assigned to the region in which the underlying probability was higher. To avoid very small regions with poor signal-to-noise ratio, ROIs having a volume less than 30 were merged into neighboring ROIs. A complete list of the brain regions and the modifications are summarized in Table 1. Fig. 2 presents the axial projection of brain regions in the glass-brain plot used to demonstrate results.

[Figure omitted. See PDF.]

Figure 2. Brain atlas regions.

Short names of atlas regions in the glass brain plot used to demonstrate results. Full names and additional information about regions can be seen in Table 1. Red spheres imply the axial projection of the center of mass of brain regions. Note that this plot does not indicate the axial depth of the regions.

https://doi.org/10.1371/journal.pone.0104947.g002

[Figure omitted. See PDF.]

Table 1. Brain atlas regions, left hemisphere.

https://doi.org/10.1371/journal.pone.0104947.t001

Calculation of regional and frame-wise displacement

We defined two metrics of displacement: regional and frame-wise displacement (RD and FD). RD time-courses were calculated as averaged voxel-wise displacements over ROIs, while FD is the analogous measurement for the entire brain. This method of calculating FD and RD is analogous to the parameter described in [14].

Quantification of global and regional BOLD intensity change

While DVARS (D referring to temporal derivative of timecourses and VARS referring to RMS variance over voxels) [31] indexes the rate of change of BOLD signal across the entire brain at each time-point of the data, Regional DVARS (RDVARS) shows the same rate for each ROI:(1)where is the image intensity at locus on frame , angle brackets denote the spatial average over all voxels within ROI () and where is the number of time-frames. In order to effectively relate RDVARS to RD, RDVARS was calculated on the timeseries following re-alignment, but prior to confound regression and filtering.

Investigating the effect of regional displacement on DVARS

We reproduced results showing global motion-related BOLD changes [2], [14], [16] by computing the correlation coefficient between FD and DVARS. Then, to distinguish global from local effects, we defined two measures, residual RD and residual RDVARS (denoted with and ), as follows:(2)and(3)where , is the number of time-frames and () identifies the region. After computing these measures for every subject and every ROI, we computed their correlation coefficient and investigated whether it depends on the degree of the global motion-BOLD relationship among subjects.

Functional connectivity processing and graph formation

For rs-fcMRI analysis, additional preprocessing steps were utilized on regional BOLD timecourses to reduce spurious variance that was unlikely to reflect neural activity. These steps included: (i) multiple regression of nuisance variables; and (ii) a temporal band-pass filter on residual data using a standard fourth-order Butterworth band-pass filter, retaining frequencies between 0.01 and 0.1 Hz.

The detailed data processing steps involved in the following strategies are discussed in the relevant works. Here, we only summarize the protocols based on basic criteria, such as the sources of nuisance signal, the number of such signal time-courses and whether global signal regression was performed.

1. NOREG: No nuisance regression,

2. WMCSF: average BOLD signals of eroded white matter and cerebrospinal fluid ROI-s, segmented using FAST [27],

3. GSREG: regression of whole-brain global signal as a covariate,

4. COMPCOR: Nuisance regression of five principal components of a noise ROI, defined based on the temporal signal-to-noise ratio, as proposed in [7],

5. NOREG+M6: NOREG + six motion parameters,

6. WMCSF+M6: WM+CSF regression + six motion parameters,

7. GSREG+M6: GSREG + six motion parameters,

8. COMPCOR+M6: COMPCOR + six motion parameters,

9. SAT36: A 36-parameter nuisance regression model proposed in [3] (incorporating global signal regression).

To ensure that neighboring regional BOLD time courses do not show spurious increase in connectivity strength, no smoothing was applied. Connectivity strength between pairs of preprocessed regional time series were calculated as the Fisher-Z transformed Pearson product moment (), and ordered into 88×88 correlation matrices for each subject (hitherto referred to as for subject ) and for each nuisance regression technique, resulting in a total of 184*9 = 1656 matrices.

Group formation

One subject with extreme in-scanner motion (average FD greater than 0.7 mm) was excluded from further analysis. The remaining 183 subjects were arranged into various groups:

* Autistic and normal control groups. Based on the clinical neuropsychological diagnostic tests for autism, as detailed in the original study description of the ABIDE dataset [15]. A group of 49 autism subjects and a neurotypical control group (n = 105) were defined.

* Groups of healthy control patents with different average voxel-wise displacement patterns, based on the group mean voxel-wise displacement maps. The healthy control group was divided into two sub-groups randomly. Temporally averaged standard-space voxel-wise displacement maps were averaged across subjects of both sub-groups separately and the spatial Pearson correlation coefficients for the group-mean voxel-wise displacement maps were calculated between the sub-groups (hitherto referred to as ). This random group formation was repeated 5000 times to generate groups-pairs with different between-group voxel-wise displacement correlation. The histogram of observed values can be seen in Fig. 3. Group-pairs were chosen so that the corresponding values are relatively low, meaning that the subjects in one sub-group tend to have different voxel-vise displacement patterns than in the other. We chose eight pairs of groups with coefficients between 0.89–0.95 (Table 2). was also computed for the autism-related group-pair.

[Figure omitted. See PDF.]

Figure 3. Histogram of group averaged voxel-wise displacement correlations based on 5000 random permutations.

Temporally averaged standard-space voxel-wise displacement maps were averaged across the subjects of two randomly assigned groups. Spatial Pearson correlation was calculated between these group-mean voxel-wise displacement maps ). The histogram of this inter-group voxel-wise displacement correlation was computed based on 5000 random group formulation. Group-pairs with extrem inter-group differences in voxel-wise displacement were chosen for further analysis.

https://doi.org/10.1371/journal.pone.0104947.g003

[Figure omitted. See PDF.]

Table 2. Voxel-vise displacement-dependent groups.

https://doi.org/10.1371/journal.pone.0104947.t002

Group comparisons

Group comparisons were performed by fitting Generalized Linear Models (GLM) [32] and arranging statistical parameters into differential statistical parametric networks (SPNs) [33]. We investigated the effect of the grouping variable on functional connectivity strength. Additional covariates that might significantly influence functional connectivity, and thus, disturb comparison, were also included in the models. These are phenotypic covariates describing age, full-scale Wechsler Abbreviated Scale of Intelligence (full IQ), gender, and subject-specific mean FD. Connectivity strength is therefore modeled as:(4)where is the measured connectivity strength between regions and for subject (element at the row and column of the matrix), is the dummy variable coding groups to compare, and , , and are the aforementioned subject-specific confounder variables, , and s are the coefficients to estimate, and is the the independent identically distributed normal error. Models were fitted utilizing an iteratively reweighted least squares (IWLS) algorithm. T-scores and p-values of the effect of interest were obtained by dividing the coefficient of interest by the estimated standard error.

Regional Displacement Interaction – RDI.

According to our hypothesis, differences in the mean RDs of region-pairs and may have an important effect on the population-level distribution of correlation coefficients corresponding to the given connection. This can be tested by adding new terms to the linear model: the temporally averaged RDs of regions and and, furthermore, the interaction term between these two covariates. Since RDs are strongly correlated with each other and with the global FD, we introduced which is an alternative to the average RD, but with FD subtracted, in order to avoid multicollinearity in the model. Accordingly, can be defined as follows:(5)or using the terminology of equation (2):(6)where denotes the subject, the time frame, the number of time frames in the fMRI time series of subject (after exclusion of first five volumes).

Thus, in case of a GLM model for the connection between regions A and B, equation (4) extends to:(7)

Henceforward, we refer to the model specified in equation (4) as a standard (STD) model and, to the terms as regional displacement interaction (RDI) and, to the model defined by equation (7) as an STD+RDI model.

Characterization of the RDI effect.

To investigate the RDI interaction effect on connectivity strength we utilized the following models:(8)and(9)which are alternative versions of models (7) and (4), respectively, with GRP variable excluded. Since models (8) and (9) are nested, we can compare the reduction in deviance to residuals utilizing an F-test under the null hypothesis that none of the additional RDI covariates in the STD+RDI model is related to the measured connectivity strength. The resulting statistical parameters for each connection were ordered into nine differential SPNs for each nuisance regression model, showing connections significantly related to RDI. Model (8) was also used to demonstrate RDI in case of a single, representative connection.

We furthermore fitted a model defined as:(10)where so that and so that and so that , or, in other words, the dependent and independent variables of model (8) are concatenated following the A and B () atlas regions. In contrast to the previously introduced models, here, we utilize one model for all connections of all subjects. The effect of interest is the RDI interaction term. This model was used to analyze the overall nature of the interaction term.

To avoid the disturbing effect of autistic differences when characterizing the interaction term, only the healthy control population was involved when applying models defined by Eqs. (8), (9) and (10).

Voxel-wise motion-related group comparison.

Differences between low and high motion groups were investigated by two GLM models for every connection: STD (Eq. (4)) and STD+RDI (Eq. (7)) models were applied where GRP was a dummy variable that defined the motion-related groups as listed in Table 2. Results were ordered into 8*9*2 t-score SPNs, summarizing motion-related group differences (eight pairs of motion-related groups, nine first-level nuisance regression methods, and two second-level regression models [STD and STD+RDI]). Since, in these comparisons, the variable of interest (grouping factor) and the motion-related covariates (FD and the RDI covariates) are potentially related, we computed the variance inflation factor (VIF) [34] for each model. VIF values suggest that these modeling approaches are free from multicollinearity issues. (Maximal observed VIF values for the variable of interest are reported in Table 2.)

Autistic-control comparison.

We tested the proposed RDI second-level interaction covariate set by comparisons of autistic and control groups defined by phenotypic information that was provided with the Autism Brain Imaging Dataset Exchange database. Results were arranged into 9*2 SPNs for nine nuisance regression methods and two second-level regression methods (STD and STD+RDI).

Computations and network visualization

Computations in this study, when not specified otherwise, were performed using R statistical programming language [35], using the packages “glm” [36], “fdrtool” [37], [38], “HH” [39] and “visreg” [40]. Differential SPNs were thresholded and visualized with the in-house developed software BrainCON (www.minipetct.com/braincon) [29].

Results and Discussion

fMRI motion artifacts have spatial predisposition

As noted by Power et al. [2], [41], the effect of motion appears to scale with the amplitude of the displacement over the whole brain: frames with greater amplitude displacements are associated with a greater change in BOLD signal. As a first step, we reproduced these results by observing a (p<0.000001) correlation between frame-wise displacement FD and DVARS.

Our results also show that this effect is not spatially constant. Although as also reported in [3], RD time courses a show high correlation in the entire brain of one subject, if the global effect is subtracted from the regional measures (FD from RD and DVRAS from RDVARS, respectively) the resulting residual measures and (Eqs. (2) and (3)) still show a significant correlation of (p<0.0001). In addition, as presented in Fig. 4, the correlation of these residual regional measures increases with the global motion-BOLD relationship throughout subjects.

[Figure omitted. See PDF.]

Figure 4. Regional motion-BOLD relationship scales with the global motion-BOLD relationship.

Within-subject correlation of and (horizontal axis) plotted against the correlation of FD and DVARS (vertical axis) for all brain regions and time frames of N = 183 subjects. The regional effect of motion on BOLD seems to scale with the global relationship which is correlated ( (p<0.000001) with the mean FD.

https://doi.org/10.1371/journal.pone.0104947.g004

In-scanner head motion consists not only of translations. Rotational components make the displacement diverse in distinct locations: it becomes greater when moving further from the center of rotations. A plausible explanation of the reported regional relationship is that this complex spatio-temporal pattern of displacement implicitly affects local BOLD signal changes.

These results can also yield a potential explanation of the phenomenon that motion tends to increase connectivity for locally adjacent nodes, but reduces connectivity between distant nodes [2], [8], [9], [14], [16]: neighboring regions having more similar RD will share more similar motion artifacts than regions being far from each other and this effect biases the distance dependence of connectivity strength.

However, we should point out that all the discussed voxel-wise and regional displacement measures are only “apparent displacements”: their estimated values may have been affected by phenomena other than head motion, including physiological noise, magnetic field inhomogeneities, instrumental instabilities, as well as BOLD activity of neuronal origin. This effect should be more pronounced for regional displacements than for the frame-wise displacements. Furthermore, the computation of RD implicitly include integrating of motion effects within each fMRI volume. Rapid head movements occurring on time scale shorter than the fMRI repetition time TR may affect different slices within an fMRI volume differently [42]. Effects of such rapid movements cause slice-specific image distortions that cannot be accurately taken into account by the volume realignment-based procedure, but can still affect fMRI functional connectivity results. Averaging over all time frames and within time frames is a simplification in modeling the spatial predisposition of head motion. However, as suggested by the significant correlation between and , this simplification seems to be reasonable. Nevertheless, sub-TR frequency components of in-scanner motion deserve more attention and their voxel-wise effect should be investigated in more detail in future publications.

The interaction of regional displacements affects measured connectivity strength

Evidences of a spatially non-constant motion artifact in brain connectivity analysis, like an increase in short-range and a decrease in long-range connectivity, or the special pattern of related changes in connectivity strength reported in [2], [8], [9], [14], [16], suggest that the reported local relationship between motion and BOLD signal changes should be considered when performing correction techniques. This is especially true for correlation-based functional connectivity analysis, where the similarity of two regional BOLD time courses can be increasingly affected by these small but systematic variations.

However, including voxel-wise motion parameters in nuisance signal regression does not seem to be efficient (as reported by [3], [14] and also found in our preliminary analysis). Yan et al. [14] used voxel-wise displacement as a reference to evaluate the differential region-specific impact of motion on the BOLD signal. Although these authors presented significant correlations with a spatial pattern similar to that previously reported, it is still not clear whether those patterns can be explained only by locally differential BOLD answers to a global motion effect, or, alternatively, by a real local relationship with the spatio-temporal motion pattern.

To investigate this question, we defined RDI, a set of second-level regression covariates that models the interaction effect between the temporally averaged regional displacements of the regions involved in the connection. In regression analysis, an interaction effect is said to exist when the effect of the focal independent variable on the dependent variable differs depending on the value of a third variable [43], called the moderator variable. (Statistically, the choice of which of the two independent variables should be the moderator variable is unimportant.) The proposed RDI interaction term in the second-level GLM model defined by Eq. (10) was found to be significant (p<0.000001) for all first-level nuisance regression methods), which means that the effect of the average RD of one region on the dependent variable (connectivity strength) changes when the average RD of the other region changes. The effect is visualized in Fig. 5 by the filled contour plots. The predicted connectivity strength changes depending on the simultaneously varying values of and , in case of all investigated first-level nuisance regression methods. These results imply that throughout the population, connectivity strength between two brain regions tends to increase if the average RD of the regions is similar (eg. both are larger or smaller than the average FD) and tends to decrease otherwise.

[Figure omitted. See PDF.]

Figure 5. Filled contour plots visualizing the Regional Displacement Interaction (RDI) effect: how the predicted connectivity strength (color-coded) changes depending on the simultaneously varying values of and , in case of no nuisance regression (A∶NOREG) and all investigated first-level nuisance regression methods, i.e., NOREG+M6 (E), WMSCF (B), COMPCORR (C), GSREG (D), WMSCF+M6 (F), COMPCORR+M6 (G), GSREG+M6 (H), and SAT36 (I).

Vertical and horizontal axes of plots B-I are the same as those of plot A. Gray bars next to the legends indicate the (−1,1) interval to ease interpretation of color-coded Z-score values.

https://doi.org/10.1371/journal.pone.0104947.g005

This effect is demonstrated for a representative connection (occipital fusiform gyrus - prefrontal gyrus) in Fig. 6. While the presented partial residual plot reveals no significant relationship between connectivity strength and , the interaction effect is significant implying that the effect of on connectivity strength differs depending on the value of . This is demonstrated by dividing the data into four groups based on the value of and visualizing the corresponding cross-sectional CCPR (component and component-plus-residual) plots, which reveal that in each group the relationship between and partial residual connectivity strength is significant in all cases but the regression lines have different slopes. Thus, this latent relationship is not observable without accounting for the interaction of the regional displacement covariates.

[Figure omitted. See PDF.]

Figure 6. RDI effect in case of a demonstrative connection: occipital fusiform gyrus (A = 20) - prefrontal gyrus (B = 12).

On the left, the partial residuals for from the model defined by Eq. (8) (dependent variable: connectivity data of the healthy control population (N = 105), without nuisance signal regression (NOREG)) are plotted against . Although the model reveals no significant relationship (t = −1.46,p = 0.15) between connectivity strength and , the interaction effect is significant (t = 3.31, p = 0.0013), implying that the effect of on connectivity strength differs depending on the value of . This is demonstrated by dividing the data into four groups based on the value of (color-coded on the left plot) and visualizing the corresponding cross-sectional CCPR (component and component-plus-residual) plots (on the left). Mean value of corresponding to the cross-section is indicated. Partial residuals are plotted with colored dots corresponding to the cross-section group. The corresponding regression line estimated from the full model fit and the corresponding 95% confidence interval is displayed in black and gray, respectively. The horizontal and vertical axes of the cross-sectional CCPR plots are the same as those of the partial residual plot on the left. Cross-sectional CCPR plots imply that in each group the relationship between and partial residual connectivity strength is significant but the regression lines have different slopes, which makes this latent relationship not observable without accounting for the interaction of the regional displacement covariates.

https://doi.org/10.1371/journal.pone.0104947.g006

To characterize how this phenomenon is related to the spatial patterns of measured functional connectivity and to what extent it is present in case of the applied first-level nuisance signal regression techniques, we performed a comparison of the STD (Eq. (8)) and STD+RDI (Eq. (9)) models. The model comparison was realized by F-tests between the models and resulted in the SPNs shown in Fig. 7. In this figure, connections are visualized, where the STD+RDI model explains significantly more variance than the STD model (the null hypothesis of the F-test can be rejected) with a false discovery rate of q<0.05. The proposed method, RDI proved to be the most efficient with nuisance signal regression methods NOREG, WMCSF, WMCSF+M6, GSREG, and COMPCOR+M6. The explanatory power added by RDI is most pronounced typically in case of middle- and long-range connections of the temporal poles.

[Figure omitted. See PDF.]

Figure 7. Network pattern of connections where utilizing RDI significantly improves second-level modeling.

Statistical parametric networks presenting the model comparison performed by F-tests between the STD and STD+RDI models (Eq. (9) and (8)). Connections are only visualized, when the STD+RDI model explains significantly more variance than the STD model (the null hypothesis of the F-test can be rejected) with a false discovery rate of q<0.05. The proposed STD+RDI method proves to be most efficient with the nuisance signal regression methods NOREG, WMCSF, WMCSF+M6, GSREG, and COMPCOR+M6, and seems to demonstrate no significant improvement in case of SAT36, after correction for multiple comparison.

https://doi.org/10.1371/journal.pone.0104947.g007

The interaction of regional displacements may bias functional connectivity group comparisons

The correlation coefficient of two regional BOLD time courses can be sensitive to the regional displacement time course (RD) of both regions. Even if the motion-related artifactual component is small, it still can significantly affect correlations, depending on the degree to which it is shared between the time courses.

This phenomenon becomes even more problematic on the population-level. Satterthwaite et al. [9] reported that between-subject differences in head motion are stable: that subjects who tend to move on one occasion tend to move on another occasion. This means that analyses of functional connectivity needs to consider the possibility that certain aspects of head motion behave as a trait. Accordingly, even if the above-mentioned effect is otherwise small, it can disturb group comparisons and lead to erroneous conclusions, since it is of non-neural origin.

This assumption can be admitted easily by considering two patient cohorts where region A and B have similar RD within each subject of one group (e.g., due to relatively smaller rotations), and different RD in the subjects of the other (e.g., more prevalent rotations with a center to which A and B are located asymmetrically). The significant RDI interaction effect means that the corresponding correlation coefficient will be biased, and tends to be larger in the first group, even in the absence of a real functional difference. However, it is still not clear how different grouping conditions interfere with the tendency of motion patterns.

One can hypothesize that the spatio-temporal pattern of motion can be significantly different between groups that were defined by a factor in relation to motion. This could be the case in group comparisons in several physiological and pathological conditions co-occurring with hypo- or hyperkinetic signs.

Our results show that the spatio-temporal pattern of head motion biases measured connectivity strength on the population level, and, practically speaking, the proposed second-level covariates (RDI) can be utilized as a method to incorporate these individual regional differences of in-scanner head motion into the model, and thus, reduce artifactual variance in the data.

Including RDI as a covariate in second-level regression efficiently reduces group differences caused by differences in voxel-wise motion

To demonstrate the reported confounding effect of voxel-wise motion on population-level analysis, we performed eight group comparisons (See Table 2) where the groups to compare have different average voxel-wise displacement patterns.

Results are summarized in Table 3. In Fig. 8, the number of significantly (p<0.01) differing connections is plotted against the group-defining mean FD threshold for each nuisance signal regression method. Results show that, for less extensive nuisance regression methods (NOREG, NOREG+M6, WMCSF, WMCSF+M6), when the voxel-wise displacement maps between groups are more different (lower between-group mean voxel-wise displacement map correlations), more group differences appear. These findings seem to confirm that group comparisons may be biased when groups show different tendencies in voxel-wise motion. Results, in conjunction with the results of the F-test based model comparisons Fig. 7, also show that the inclusion of RDI covariates seems to decrease these artifactual group differences, especially by moderate nuisance regression methods.

[Figure omitted. See PDF.]

Figure 8. The effect of RDI on motion-related group differences.

The number of significantly (p<0.01) differing connections is plotted against the spatial correlation coefficient of group-mean voxel-wise displacement maps for RD-based group-pairs for each nuisance signal regression method. The number of significant group differences is plotted on logarithmic axis. Improvement in the reduction of motion-related group differences was most pronounced for the NOREG, WMCSF and COMPCOR methods. Nuisance regression methods incorporating GSREG seem not to be sensible for differences in group-mean voxel-wise displacement patterns however, they show relatively high number of group differences by all group comparisons.

https://doi.org/10.1371/journal.pone.0104947.g008

[Figure omitted. See PDF.]

Table 3. Motion-related group differences.

https://doi.org/10.1371/journal.pone.0104947.t003

The change in the corresponding connectivity pattern is visualized in Fig. 8 for the nuisance signal regression methods NOREG, COMPCOR, and GSREG, and for each motion-related group comparison defined by . Results with the STD and STD+RDI second-level regression models are presented in the upper and lower rows on each panel, respectively. Group differences with probability (p> = 0.01) are not visualized.

As predicted by the F-test-based model comparisons (Fig 7), the reduction of false group differences caused by voxel-wise motion most markedly improved in the NOREG, WMCSF and COMPCOR methods. Regressing out six motion parameters did not seem to be highly efficient but may be beneficial when no other nuisance signal regression covariate is applied (NOREG+M6).

However, the small number of differences surviving the q<0.05 false discovery rate criterion implies that - when utilizing a proper second-level model - all the investigated methods are able to reduce motion-related group comparison artifacts to a decent extent. With the RDI correction, no FDR significant group differences remained in any of the comparison cases.

We note, that, in contrast to most prior studies, here, we performed more than a single pair of group comparisons, thus avoiding that results reflect only the random effects of the grouping condition.

Including RDI as a covariate in second-level regression preserves autism-related group differences

To test the efficiency of the proposed correction method, we performed group comparisons where both motion-related artifacts [20], [21] and real neuronal differences [10], [15], [44] were expected to be present. We compared the functional networks of autistic and control patients.

The correlation between the group-mean voxel-wise correlation map between the autistic and the normal control groups was found to be 0.98. It suggests that this grouping condition should be only slightly biased by the effect of regional displacement on the measured connectivity strength. Thus, as predicted by voxel-wise displacement pattern-related group comparisons, utilizing the RDI correction method should introduce only minor changes in the differential connectivity patterns.

Results are presented in Fig. 9 and Tables 4 and 5.

[Figure omitted. See PDF.]

Figure 9. Autism related group comparisons.

Autism related group differences for the investigated correction strategies. Colors denote significance levels as detailed in the legend.

https://doi.org/10.1371/journal.pone.0104947.g009

[Figure omitted. See PDF.]

Table 4. Autism-related group differences (without nuisance signal regression of six motion parameters).

https://doi.org/10.1371/journal.pone.0104947.t004

[Figure omitted. See PDF.]

Table 5. Autism-related group differences (with nuisance signal regression of 6 motion parameters).

https://doi.org/10.1371/journal.pone.0104947.t005

More than 200 connectivity differences survived the q<0.05 false discovery rate criterion in NOREG and COMPCOR (both with and without RDI) and none survived in GSREG+RDI, GSREG+M6, GSREG+M6+RDI, SAT36 and SAT36+RDI. As predicted above, the inclusion of RDI introduces only slight changes in the pattern of autism-related group differences.

All the evaluated signal regression approaches revealed presumably autism-linked impairments of functional connectivity. Autism was mainly characterized by decreased synchronicity, i.e., under-connectivity. This finding is in line with the majority of intrinsic functional connectivity studies. The spatial predisposition, along with the 30 most significant differences, is presented in Fig. 9 and Table 4 and 5.

While including RDI significantly reduced (presumably artifactual) differences between voxel-wise displacement-related subject cohorts, differences in the autism-related comparison were more or less preserved. These results suggest that the proposed correction method, while effectively reducing motion artifacts in group comparisons, preserves the sensitivity to neural differences.

A critical interpretation of our autism-linked findings in this study is not directly possible. This mainly stems from the lack of ground truth information about the basic neuropathology of the disease. Furthermore, larger-scale, multi-centric comparisons would be optimal to test the reproducibility of any finding.

Contrasting the impact of head motion and nuisance signal regression strategies

The efficiency of nuisance signal regression techniques in the context of rs-fcMRI motion artifacts analysis has been intensively investigated in the last few years; however, a significant part of these studies did not apply second-level regression covariates [2], [3], [9], [16] in population-level analysis. Their overall conclusion was that global signal regression, high-parameter nuisance signal regression, scrubbing, and de-spiking are potentially beneficial. However, recent studies have questioned some of these conclusions. In the following, we summarize the latest findings in the literature in contrast to our results.

Effect of global signal regression.

In our analysis, patterns of group differences become extremely different when regressing out the whole-brain signal from regional BOLD time courses (GSREG, GSREG+M6 and SAT36). Most of the above-mentioned studies, which did not utilize second-level correction, applied GSREG in their analysis pipeline. Recent studies using mean FD as a second-level regressor [14], [45] also concluded that GSREG mitigates the effects of motion-related differences among subjects, but warns that investigators must weigh up the pros and cons of GSREG when deciding whether to employ it in the context of testing specific hypotheses.

Our results show that, after the correction of group-wise differences in head movement in the autism-control comparison, over-connectivity and under-connectivity can simultaneously appear. However, global signal regression in the processing pipeline appears to bias results toward over-connected differential networks. In addition to under-connectivity, many authors suggested short-distance over-connectivity (in the frontal lobe or globally) as a possible finding in autism spectrum disorders [46]. This theory was then questioned on the grounds that head-movement could induce similar effects [21]. As reported in [47], correlation estimates obtained after GSREG are more susceptible to the presence of motion and exacerbate distance-dependent bias. Moreover, as reported in [20], [22], correlation patterns and group differences may become distorted after GSREG (depending on, e.g., region size or the underlying true connectivity structure). According to Müller et al. [48], pre-processing strategies greatly affect the spatial patterns of autism-linked connectivity traits, although under-connectivity is the most prevalent across studies. It is, therefore, safe to conclude that GSREG not only introduces anti-correlations in the functional connectome [49], but can also confound case-control comparisons in autism.

Effect of including motion parameters.

In the case of voxel-wise motion pattern-related group comparisons, inclusion of the six motion parameters in first-level nuisance signal regression showed no obvious improvement in motion-artifact reduction. However, the inclusion of these parameters had a pronounced effect in autism-related group comparisons, especially in NOREG, WMCSF and COMPCOR techniques: the number of significant group differences decreased. Whether this phenomenon is due to improved motion-artifact reduction in the special case of autistic group comparisons, or due to overfitting and the removal of real brain signal, remains a question. However, evidence that spatial realignment-based estimation of motion parameters may yield poor results in periods of small movements [12] point toward the conclusion that motion parameter estimates should be applied carefully. As a possible solution, we suggest utilizing thresholded motion estimates and avoiding the use of motion parameters in periods of relatively small movements.

Optimal choice of individual-level motion-correction technique.

In some cases, it is not clear whether various high-parameter nuisance regression techniques eliminate group differences due to increased specificity or decreased sensitivity. Considering most significant autism-related group differences, minimal nuisance regression techniques (NOREG, NOREG+M6) show a high consensus with moderate (WMCSF, WMSCF+M6) and more complex (COMPCOR, COMPCOR+M6) methods. This points to the conclusion that when analyzing a sufficiently large sample and utilizing an appropriate second-level model, the choice of individual-level signal nuisance regression technique becomes less crucial. However, when analyzing small samples and also with individual analysis, the role of these techniques is unquestionably important.

The role of confounds not related to subject motion.

This article focused on motion-related artifacts, which are only one, although a conspicuous source of confounding effects in functional MRI. In thes context, the performance of NOREG and NOREG+M6 methods in motion-related comparisons deserves attention. One possible explanation is that, although artifacts of other sources are obviously present in the data, their individual spatio-temporal pattern is more constant and their population-level distribution is more similar among the investigated subpopulations compared to motion-related confounds. Thus, these artifacts may have only a moderate disturbing effect in large-sample group comparisons. However, by grouping conditions in relation to physiological conditions, like blood pressure and blood oxygenation (or artifacts of scanner-related sources in multi-center studies), non-motion originated artifacts can appreciably affect the results. Thus in such experimental designs, nuisance regression methods may have a more important role.

This is also suggested by the results of autism-related group comparisons, where pronounced differences were experienced among various nuisance signal regression techniques. A potential explanation for these deviations is that, although the proposed motion-correction technique successfully reduces motion-related erroneous group differences, artifacts of other physiological and scanner-related sources affect autism-related group comparisons and these phenomenon is handled differently by various signal regression techniques.

Conclusion

In this study, we demonstrated that small movements during scanning can cause different displacements in various locations of the brain, and, accordingly, motion-related BOLD signal changes also depend on location. We characterized the effect of this spatio-temporally complex BOLD artifact pattern on functional connectivity. We proposed RDI, a set of regression covariates for the population-level correction of motion artifacts arising from local head motion. As shown with comparisons of groups with differing average voxel-wise motion pattern, the proposed correction technique efficiently reduces artifacts caused by differences in voxel-wise motion patterns in population-based connectivity analysis; and meanwhile, as demonstrated by comparing autistic and control groups, preserves differences corresponding to neural origin. Our findings suggest that, especially by moderate nuisance correction methods, the inclusion of RDI as second-level nuisance covariates is generally appropriate and may become increasingly necessary when the variable of interest is interrelated with altered subject kinetics.

A limitation of the proposed method is that it cannot be effectively applied in case of individual studies or small sample sizes. Nevertheless, one should note that, in situations where the variable of interest is correlated with motion, second-level regression-based, motion-correction approaches can be conservative, as they remove common variation among regressors. Furthermore, the proposed method is based only on a simplified measure of motion and does not handle rapid sub-TR displacements, which may play an important role in regional motion artifact interactions.

The question of what is the optimal individual-level signal regression technique for motion correction remains open, but, seems less crucial for large-sample, group-level studies using a proper second-level correction method.

This article focused on motion-related artifacts, which are only one, although a conspicuous source, of confounding effects in functional MRI. In-scanner head motion is relatively easy to measure, and thus, corresponding artifacts are actively investigated. However, as suggested by the observed differences among nuisance signal regression techniques, physiological and scanner-related artifacts may also have an essential impact on fc-MRI studies.

Author Contributions

Conceived and designed the experiments: TS ME AJ EB. Analyzed the data: TS AJ CSA. Contributed reagents/materials/analysis tools: GO SK CSA. Wrote the paper: TS AJ ME.

Citation: Spisák T, Jakab A, Kis SA, Opposits G, Aranyi C, Berényi E, et al. (2014) Voxel-Wise Motion Artifacts in Population-Level Whole-Brain Connectivity Analysis of Resting-State fMRI. PLoS ONE 9(9): e104947. https://doi.org/10.1371/journal.pone.0104947

References

1. Kim B, Boes JL, Bland PH, Chenevert TL, Meyer CR (1999) Motion correction in fmri via registration of individual slices into an anatomical volume. Magnetic Resonance in Medicine 41: 964–972.

2. Power JD, Barnes KA, Snyder AZ, Schlaggar BL, Petersen SE (2012) Spurious but systematic correlations in functional connectivity mri networks arise from subject motion. Neuroimage 59: 2142–2154.

3. Satterthwaite TD, Elliott MA, Gerraty RT, Ruparel K, Loughead J, et al. (2013) An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data. Neuroimage 64: 240–256.

4. Friston KJ, Holmes AP, Worsley KJ, Poline JP, Frith CD, et al. (1994) Statistical parametric maps in functional imaging: a general linear approach. Human brain mapping 2: 189–210.

5. Fox MD, Zhang D, Snyder AZ, Raichle ME (2009) The global signal and observed anticorrelated resting state brain networks. Journal of neurophysiology 101: 3270–3283.

6. Glover GH, Li TQ, Ress D (2000) Image-based method for retrospective correction of physiological motion effects in fmri: Retroicor. Magnetic Resonance in Medicine 44: 162–167.

7. Behzadi Y, Restom K, Liau J, Liu TT (2007) A component based noise correction method (compcor) for bold and perfusion based fmri. Neuroimage 37: 90–101.

8. Fair DA, Nigg JT, Iyer S, Bathula D, Mills KL, et al. (2012) Distinct neural signatures detected for adhd subtypes after controlling for micro-movements in resting state functional connectivity mri data. Frontiers in systems neuroscience 6: 80.

9. Satterthwaite TD, Wolf DH, Loughead J, Ruparel K, Elliott MA, et al. (2012) Impact of in-scanner head motion on multiple measures of functional connectivity: relevance for studies of neurodevelopment in youth. Neuroimage 60: 623–632.

10. Jakab A, Beres M, Berenyi E, Emri M, Spisak T, et al. (2013) Autistic traits in neurotypical adults: Correlates of graph theoretical functional network topology and white matter anisotropy patterns. PloS one 8: e60982.

11. Friston KJ, Williams S, Howard R, Frackowiak RS, Turner R (1996) Movement-related effects in fmri time-series. Magnetic resonance in medicine 35: 346–355.

12. Freire L, Mangin JF (2001) Motion correction algorithms may create spurious brain activations in the absence of subject motion. NeuroImage 14: 709–722.

13. Lemieux L, Salek-Haddadi A, Lund TE, Laufs H, Carmichael D (2007) Modelling large motion events in fmri studies of patients with epilepsy. Magnetic resonance imaging 25: 894–901.

14. Yan CG, Cheung B, Kelly C, Colcombe S, Craddock RC, et al. (2013) A comprehensive assessment of regional variation in the impact of head micromovements on functional connectomics. Neuroimage 76: 183–201.

15. Di Martino A, Yan C, Li Q, Denio E, Castellanos F, et al. (2013) The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism. Molecular psychiatry 19: 659–667.

16. Van Dijk KR, Sabuncu MR, Buckner RL (2012) The influence of head motion on intrinsic functional connectivity mri. Neuroimage 59: 431–438.

17. Jo HJ, Saad ZS, Simmons WK, Milbury LA, Cox RW (2010) Mapping sources of correlation in resting state fmri, with artifact detection and removal. Neuroimage 52: 571–582.

18. Cordes D, Haughton VM, Arfanakis K, Carew JD, Turski PA, et al. (2001) Frequencies contributing to functional connectivity in the cerebral cortex in resting-state data. American Journal of Neuroradiology 22: 1326–1333.

19. Van Someren E, Van Der Werf YD, Roelfsema PR, Mansvelder HD, Lopes da Silva FH (2011) Spectral characteristics of resting state networks. Slow Brain Oscillations of Sleep, Resting State and Vigilance 193: 259.

20. Gotts SJ, Saad ZS, Jo HJ, Wallace GL, Cox RW, et al. (2013) The perils of global signal regression for group comparisons: a case study of autism spectrum disorders. Frontiers in human neuroscience 7: 356.

21. Deen B, Pelphrey K (2012) Perspective: brain scans need a rethink. Nature 491: S20–S20.

22. Saad ZS, Gotts SJ, Murphy K, Chen G, Jo HJ, et al. (2012) Trouble at rest: how correlation patterns and group differences become distorted after global signal regression. Brain connectivity 2: 25–32.

23. Imperati D, Colcombe S, Kelly C, Di Martino A, Zhou J, et al. (2011) Differential development of human brain white matter tracts. PloS one 6: e23437.

24. Koyama MS, Di Martino A, Zuo XN, Kelly C, Mennes M, et al. (2011) Resting-state functional connectivity indexes reading competence in children and adults. The Journal of Neuroscience 31: 8617–8624.

25. Jenkinson M, Beckmann CF, Behrens TE, Woolrich MW, Smith SM (2012) Fsl. NeuroImage 62: 782–790.

26. Jenkinson M, Bannister P, Brady M, Smith S (2002) Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage 17: 825–841.

27. Smith SM (2002) Fast robust automated brain extraction. Human brain mapping 17: 143–155.

28. Jenkinson M, Smith S (2001) A global optimisation method for robust affine registration of brain images. Medical image analysis 5: 143–156.

29. Spisák T, Opposits G, Kis SA, Clemens B, Emri M (2012) Braincon: Software tool for graph theory based multimodal brain connectivity analysis and visualization. Magn Reson Mat Phys Biol Med 25: 616.

30. Desikan RS, Ségonne F, Fischl B, Quinn BT, Dickerson BC, et al. (2006) An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest. Neuroimage 31: 968–980.

31. Smyser CD, Inder TE, Shimony JS, Hill JE, Degnan AJ, et al. (2010) Longitudinal analysis of neural network development in preterm infants. Cerebral Cortex 20: 2852–2862.

32. Dobson AJ (2010) An introduction to generalized linear models. CRC press.

33. Ginestet CE, Simmons A (2011) Statistical parametric network analysis of functional connectivity dynamics during a working memory task. Neuroimage 55: 688–704.

34. Stine RA (1995) Graphical interpretation of variance inflation factors. The American Statistician 49: 53–56.

35. R Core Team (2013) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/. R package version 3.1.0.

36. Venables WN, Ripley BD, Venables W (1994) Modern applied statistics with S-PLUS, volume 250 . Springer-verlag New York.

37. Klaus B, Strimmer K (2013) fdrtool: Estimation of (Local) False Discovery Rates and Higher Criticism. URL http://CRAN.R-project.org/package=fdrtool. R package version 3.1.0.

38. Strimmer K (2008) fdrtool: a versatile r package for estimating local and tail area-based false discovery rates. Bioinformatics 24: 1461–1462.

39. Heiberger RM, Holland B (2004) Statistical Analysis and Data Display: An Intermediate Course with Examples in S-Plus, R, and SAS. Springer-Verlag, New York, first edition. URL http://springeronline.com/0-387-40270-5.

40. Breheny P, Burchett W (2012) Visualization of Regression Models Using visreg. URL http://myweb.uiowa.edu/pbreheny/publications/visreg.pdf. R package version 3.1.0.

41. Power JD, Barnes KA, Snyder AZ, Schlaggar BL, Petersen SE (2012) Steps toward optimizing motion artifact removal in functional connectivity mri; a reply to carp. Neuroimage 76: 439441.

42. Zotev V, Yuan H, Phillips R, Bodurka J (2012) Eeg-assisted retrospective motion correction for fmri: E-remcor. NeuroImage 63: 698–712.

43. Allison PD (1977) Testing for interaction in multiple regression. American Journal of Sociology 83: 144.

44. Belmonte MK, Allen G, Beckel-Mitchener A, Boulanger LM, Carper RA, et al. (2004) Autism and abnormal development of brain connectivity. The Journal of Neuroscience 24: 9228–9231.

45. Yan CG, Craddock RC, He Y, Milham MP (2013) Addressing head motion dependencies for small-world topologies in functional connectomics. Frontiers in human neuroscience 7: 910.

46. Courchesne E, Pierce K (2005) Why the frontal cortex in autism might be talking only to itself: local over-connectivity but long-distance disconnection. Current opinion in neurobiology 15: 225–230.

47. Jo HJ, Gotts SJ, Reynolds RC, Bandettini PA, Martin A, et al. (2013) Effective preprocessing procedures virtually eliminate distance-dependent motion artifacts in resting state fmri. Journal of applied mathematics 2013: 9.

48. Müller RA, Shih P, Keehn B, Deyoe JR, Leyden KM, et al. (2011) Underconnected, but how? a survey of functional connectivity mri studies in autism spectrum disorders. Cerebral Cortex 21: 2233–2243.

49. Murphy K, Birn RM, Handwerker DA, Jones TB, Bandettini PA (2009) The impact of global signal regression on resting state correlations: are anti-correlated networks introduced? Neuroimage 44: 893–905.