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

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.


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][5][6][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][9][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 spikeregression [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][9][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.

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.36161.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 = 363 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 timeseries 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 timeseries (see Fig. 1 for demonstrative images), and then realigned to standard space.

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 cm 3 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.

Calculation of regional and frame-wise displacement
We defined two metrics of displacement: regional and framewise 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 FD vox 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:  Lateralization, and long and short names of brain atlas-based ROIs used for estimating regional BOLD and motion related measures. The sources of brain regions are the Harvard-Oxford Cortical and Subcortical brain atlases [30]. Probability maps for all regions were accessed and region borders were delineated by kretaining 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. ROIs having a volume less than 30 cm 3 were merged into neighboring ROIs, as indicated by column five. The table lists only regions in the left hemisphere. The naming conventions and the region merging procedure was analogous for their contra-lateral pairs. doi:10.1371/journal.pone.0104947.t001 Voxel-Wise Motion Artifacts in fMRI Brain Connectivity Analysis where I i (x) is the image intensity at locus x on frame i, angle brackets S.T (R) denote the spatial average over all voxels x within ROI R (R~1 . . . 88) and i~2 . . . N where N 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 DRD and DRDVARS), as follows: and where i~2 . . . N, N is the number of time-frames and R (R~1 . . . M) 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 fourthorder 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.
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 (z~atanh(r)), and ordered into 88688 correlation matrices for each subject (hitherto referred to as C i for subject i ) 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: N 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. N Groups of healthy control patents with different average voxelwise displacement patterns, based on the group mean voxelwise 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 r W D ). This random group formation was repeated 5000 times to generate groups-pairs with different between-group voxel-wise displacement correlation. The histogram of observed r W D values can be seen in Fig. 3. Group-pairs were chosen so that the corresponding r W D values are relatively low, meaning that the subjects in one subgroup tend to have different voxel-vise displacement patterns than in the other. We chose eight pairs of groups with r W D 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. doi:10.1371/journal.pone.0104947.g002 coefficients between 0.89-0.95 (Table 2). r W D was also computed for the autism-related group-pair.

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: where y (A,B) i is the measured connectivity strength between regions A and B for subject i (element at the A th row and B th column of the C i matrix), GRP i is the dummy variable coding groups to compare, and AGE i , FIQ i , SEX i and FD i are the aforementioned subject-specific confounder variables, a, b and i s are the coefficients to estimate, and E i is the the i th 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 b coefficient of interest by the estimated standard error.
Regional Displacement Interaction -RDI. According to our hypothesis, differences in the mean RDs of region-pairs A and B 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 A and B and, furthermore, the interaction term between these two covariates. Since RDs are strongly correlated with each other and with the global FD, we introduced DRD which is an alternative to the average RD, but with FD subtracted, in order to avoid multicollinearity in the model. Accordingly, DRD can be defined as follows: or using the terminology of equation (2): 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 (r W D )). 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. doi:10.1371/journal.pone.0104947.g003 where i denotes the subject, t the time frame, T i the number of time frames in the fMRI time series of subject i (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: Henceforward, we refer to the model specified in equation (4) 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: 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: where Y j~Y 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 motionrelated 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 denotes the obtained maximal variance inflation factor (VIF) (throughout all connections) corresponding to the grouping factor, in models STD and STD+RDI, respectively. None of the groups introduce multicollinearity in the models. (however there is a slight difference between the FDs of the group-pairs with a spatial voxel-wise displacement correlation of 0.98). doi:10.1371/journal.pone.0104947.t002 database. Results were arranged into 9*2 SPNs for nine nuisance regression methods and two second-level regression methods (STD and STD+RDI).

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 r~0:485 (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 DRD (R) i and DRDVARS (R) i (Eqs. (2) and (3)) still show a significant correlation of r~0:278 (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.
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 framewise 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 DRD (R) i and DRDVARS (R) i , 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 secondlevel 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. i , 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. 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 DRD A , the DRD A DRD B interaction effect is significant implying that the effect of DRD A on connectivity strength differs depending on the value of DRD B . This is demonstrated by dividing the data into four groups based on the value of DRD B and visualizing the corresponding crosssectional CCPR (component and component-plus-residual) plots, which reveal that in each group the relationship between DRD A 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.
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.

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 betweensubject 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.
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 motionrelated group comparison defined by r W D . 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 (r W D ) 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.
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 . Mean value of DRD B corresponding to the cross-section is indicated. Partial residuals are plotted with colored dots corresponding to the crosssection 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 DRD A 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. doi:10.1371/journal.pone.0104947.g006 Voxel-Wise Motion Artifacts in fMRI Brain Connectivity Analysis [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 abovementioned 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- 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. doi:10.1371/journal.pone.0104947.g007 Table 3. Motion-related group differences. 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 highparameter nuisance regression techniques eliminate group differ-ences 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 largesample 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 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. doi:10.1371/journal.pone.0104947.g008 Figure 9. Autism related group comparisons. Autism related group differences for the investigated correction strategies. Colors denote significance levels as detailed in the legend. doi:10.1371/journal.pone.0104947.g009 Table 4. Autism-related group differences (without nuisance signal regression of six motion parameters).  Table 5. Autism-related group differences (with nuisance signal regression of 6 motion parameters). Voxel-Wise Motion Artifacts in fMRI Brain Connectivity Analysis 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 spatiotemporally 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 scannerrelated artifacts may also have an essential impact on fc-MRI studies.