The present study is a collaboration between the Center for Applied Genomics (CAG) at Children's Hospital of Philadelphia (CHOP) and the Brain Behavior Laboratory at the University of Pennsylvania (Penn). Study procedures were reviewed and approved by the Institutional Review Board of both CHOP and Penn. The target population-based sample is of 10,000 youths who presented to the CHOP network for a pediatric visit and volunteered to participate in genomic studies of complex pediatric disorders. All research participants undergo assessment with a structured neuropsychiatric interview, a neuroscience based computerized neurocognitive battery (Gur et al., 2010
), and review of electronic medical records. A subsample of 1,000 subjects is selected for neuroimaging. This report represents an interim analysis of the initial consecutively acquired subjects that underwent neuroimaging, including 456 individuals aged 8-23 years old (mean age 15.6 years; standard deviation 3.4 years); 199 were male. All imaged subjects were recruited as part of the larger study examining neurodevelopmental genomics. Inclusion criteria for the overall study
: 1. Able to provide signed informed consent. For participants under age 18 assent and parental consent were required. 2. English proficiency. 3. Physically and cognitively able to participate in computerized neurocognitive testing. Notably, the imaged sample was not recruited on the basis of any specific disorder or symptoms. All imaged subjects were free from active, ongoing medical disorders that would impact brain function; all had normal or corrected-to-normal vision and met no exclusionary criteria for a MRI study.
2.2 Image acquisition
All subject data were acquired on the same scanner (Siemens Tim Trio 3 Tesla, Erlangen, Germany; 32 channel head coil) using the same imaging sequences. Blood oxygen level dependent (BOLD) fMRI was acquired using a whole-brain, single-shot, multi-slice, gradient-echo (GE) echoplanar (EPI) sequence of 124 volumes with the following parameters: TR/TE=3000/32 ms, flip=90°, FOV=192×192 mm, matrix=64×64, slice thickness/gap =3mm/0mm. The resulting nominal voxel size was 3.0×3.0×3.0mm. A fixation cross was displayed as images were acquired. Subjects were instructed to stay awake, keep their eyes open, fixate on the displayed crosshair, and remain still. Prior to time-series acquisition, a 5-minute magnetization-prepared, rapid acquisition gradient-echo T1-weighted (MPRAGE) image (TR 1810ms, TE 3.51 ms, FOV 180×240 mm, matrix 256×192, effective voxel resolution of 1 × 1 × 1mm) was acquired to aid spatial normalization to standard atlas space. Prior to scanning, in order to acclimate subjects to the MRI environment, a mock scanning session was conducted for each individual using a decommissioned MRI scanner and head coil. Mock-scanning was accompanied by acoustic recordings of the noise produced by gradient coils for each scanning pulse sequence. During these sessions, feedback regarding head movement was provided using the MoTrack (Psychology Software Tools, Inc, Sharpsburg, PA) motion tracking system. In order to further minimize motion, subjects' heads were stabilized in the head coil using one foam pad over each ear and a third over the top of the head.
2.3 Assessment of in-scanner head motion
As in Van Dijk et al., our primary measurement of in-scanner head motion was mean relative volume-to-volume displacement. This standard measure summarizes total volume-to-volume translation and rotation across all three axes. While this summary measure was used in all of our analyses, we also examined the correlation between mean relative displacement and other commonly used movement parameters, including maximum relative displacement, standard deviation of relative displacements, maximum absolute displacement, mean absolute displacement, and standard deviation of the absolute displacement. Here, absolute displacement refers to the displacement versus a median reference volume, rather than displacement when compared to the previously acquired volume as in measures of relative displacement.
2.4 Subject sub-samples
As described below (see section 3.5), subject age and subject motion were highly related in the complete sample of 456 subjects, and this sample included subjects with gross motion that would be excluded from any typical analysis. Therefore, in order to examine the effects of motion on connectivity without age being a confounding factor, we selected an age/motion-unrelated subsample of 348 subjects (age range 8-23 years; mean age 16.6 y, SD 3.0 y; 146 male). Following the exclusions of subjects with gross motion (see below), age and motion were still significantly correlated in this age/motion-related subsample. This relationship was driven primarily by young subjects who had high levels of motion (see ). By excluding subjects with relatively high levels of motion (i.e., +2SD MRD), this relationship was attenuated (as high motion subjects tended to be younger), but a significant correlation between motion and age remained. Therefore, we next sequentially removed younger subjects with higher levels of motion until the absolute value of the correlation between age and motion was less than or equal to r=0.03. In this sample mean relative displacement was not significantly different among age groups when binned by quartile. We used this age/motion-unrelated subsample for both replication of the within-network effects described by Van Dijk et al., (see section 2.7) as well as the subsequent analyses that investigated approaches beyond within-network connectivity (e.g., modularity, dual-regression ICA, ALFF/fALFF).
Figure 5 Relationship between age and in-scanner head motion. (A) Relative mean displacement by age of the complete original sample of 456 subjects. (B) Relationship between exclusion threshold and % of sample retained (blue) and the correlation between age and (more ...)
As our final aim of this study was to investigate how in-scanner motion might bias estimates of how connectivity changes with age, we compared the effect of age in this age/motion-unrelated
sample to a sample where age and motion were related
. This age/motion-related
sample of 421 subjects (age range 8-23 years; mean age 15.9 y, SD 3.3 y; 181 male) was a subset of the complete sample of 456 subjects, created by simply excluding subjects with gross in-scanner motion, without any attempt to match age and motion. Gross motion was defined as relative mean displacement >0.55mm. While the selection of any exclusion threshold is ultimately arbitrary, this threshold was in line with prior work, and produced a sample that had a similar amount of motion to other large-studies of developmental connectivity (Fair et al., 2008
; Fair et al., 2007
). Note that these subjects with gross motion (>0.55mm) were excluded from all analyses, and also were not part of the age/motion-unrelated
subsample. Summary motion statistics for the inclusive sample and each subsample are presented in .
2.5 General image preprocessing
While each analytic approach required specialized preprocessing, some steps were common to all. All fMRI data processing was conducted using tools that are part of FSL (FMRIB's Software Library, www.fmrib.ox.ac.uk/fsl
). BET was used to remove non-brain areas (Smith, 2002
). The first four volumes were removed to allow BOLD signal stabilization. All functional timeseries were slice-time corrected, motion corrected to the median image using a tri-linear interpolation with six degrees of freedom (Jenkinson et al., 2002
), spatially smoothed (6mm FWHM, isotropic), and grand-mean scaled using mean-based intensity normalization.
2.6 Preprocessing for seed analyses
In additional to the general preprocessing described above (section 2.5), for both within-network and whole-brain seed analyses several additional steps were performed. During preprocessing, data were band-pass filtered to retain frequencies between 0.01Hz and 0.1 Hz (van den Heuvel et al., 2008
). As in Van Dijk et al. and multiple prior studies (Andrews-Hanna et al., 2007
; Fox et al., 2006
; Vincent et al., 2006
), we removed the influence of several confounding signals using linear regression. These included: 1) six motion parameters, 2) mean white matter timecourse, 3) mean CSF timecourse, and 4) mean whole-brain timecourse. CSF and white matter timecourses were extracted from subject-specific tissue-segments of the T1-weighted image created using FAST (Zhang et al., 2001
). Following regression of confound timecourses, the residual functional timeseries for each subject was co-registered with the anatomical image and transformed to standard anatomical space (T1 Montreal Neurological Institute template, voxel dimensions of 2×2×2 mm) using a non-linear normalization algorithm (FNIRT).
2.7 Replication of Van Dijk et al.: within-network connectivity
We conducted seed-based analyses in order to independently replicate the results reported by Van Dijk et al. (2011)
. We examined the effect of motion on pairwise connectivity relationships in two of the networks investigated in that study: the default mode network and the frontoparietal control network. In order to replicate their results, we used the same nodes for the default mode network (originally based on Raichle et al., 2001
), which included the posterior cingulate cortex (MNI coordinates: 0, −53, 26), medial prefrontal cortex (0, 54, −4), and left and right inferior parietal lobule (−48, −62, 36 and 50, −62, 32).The frontoparietal control network (originally based on Vincent et al., 2008
) was also defined as per Van Dijk et al. (2011)
, and included the anterior prefrontal cortex (−36, 57, 3 and 36, 57, 3) and inferior parietal lobule (−44,−52, 54 and 48, −50, 52). Data were preprocessed as described in sections 2.5 and 2.6. Subsequently, at each ROI, we extracted the mean timecourse of the standard-space residual timeseries. This ROI timeseries was correlated with the timeseries for every other ROI within each network using Pearson's correlations, generating a connectivity matrix describing each network for every subject. In order to improve normality, a Fisher's r-to-z transformation was applied prior to statistical testing. As in Van Dijk et al., averaging across all within-network connections created a composite of total within-network connectivity; this average value was correlated with mean relative displacement. Nonlinear effects were investigated using quadratic regression. All statistical analyses of non-voxelwise data were conducted in Matlab 7.11 (Math Works, Natick, MA, USA) using the Statistical Toolbox. Unless specified otherwise, alpha for all comparisons was set to p=0.05.
2.8 Whole-brain network definition
Next, in order to further examine the effect of motion on measures of whole-brain network topology, we calculated the pairwise connectivity relationship among a much larger set of ROIs that covered the entire brain. Recent evidence suggests that anatomic ROIs may contain functionally heterogeneous sub-regions, leading to the mixing of disparate signals, with resultant negative consequences for network estimation (Smith et al., 2011
). Accordingly, we used 160 ROIs (5mm radius spheres) that were derived from a meta-analysis of a large sample of task-based fMRI studies (Dosenbach et al., 2010
). Subsequent processing was otherwise the same as above.
2.9 Modulation of motion effects by inter-node distance
Van Dijk et al. demonstrated that motion causes increased nonspecific coupling at ranges <12mm compared to longer ranges, but did not quantify the effect of motion on connectivity over a full range of node distances. Using the 160×160 connectivity matrix, we explored the degree to which inter-node distance modulates the effect of motion on a continuous basis.Mean relative displacement was correlated across subjects with the connectivity measure for every node pair, generating a second-level 160×160 correlation matrix containing 12,720 unique cells. This matrix described how motion impacted pairwise connectivity. In order to establish how distance influenced this relationship, a linear mixed effect regression model was fit using inter-node distance as a fixed factor. Node variability was modeled as a random factor. To ascertain if non-linear effects were present, a second model was run that included both linear and quadratic effects of distance. In addition, this effect was examined separately in positive and negative connections. Lastly, to examine the influence of subject age, we conducted the analyses separately in younger and older subjects as defined by a median split of subject age.
2.10 Modularity analysis
Within the whole-brain network of 160 ROIs, we evaluated the effect of motion on network modularity. Modularity provides a measure of the degree to which a given set of nodes can be parsed into coherent sub-networks (Rubinov and Sporns, 2010
). Prior studies of network topology have used thresholded correlation matrices (Achard et al., 2006
; Bassett et al., 2008
). However, threshold choice may be arbitrary, and significantly impacts estimates of network topology (Rubinov and Sporns, 2011
). We therefore used the threshold-free modularity analysis recently introduced by Rubinov and Sporns (2011)
, part of the Brain Connectivity Toolbox (https://sites.google.com/a/brain-connectivity-toolbox.net/
). This measure allows estimation of modularity (using the Louivan algorithm) for fully connected networks with positive and negative weights. Modularity was linearly related to relative mean displacement using Pearson's correlations. Nonlinear effects were evaluated with quadratic regression.
2.11 Independent component analyses
Temporal concatenation group independent component analysis (TC-GICA) was conducted on the age/motion related sample using MELODIC (Beckmann et al., 2005
). Prior to TC-GICA, in addition to the general preprocessing steps described above (Section 2.5), data were high pass filtered (0.01Hz) to remove drift, spatially smoothed at 6mm FWHM, and resampled to 4mm isotropic voxels in MNI template space. Component dimensionality was estimated automatically. All components resulting from TC-GICA were submitted for subsequent dual-regression analysis. We focused our analysis on three commonly-observed components, corresponding to the default mode network, the right-lateralized frontoparietal control network, and the dorsal attention network (Biswal et al., 2010
; Damoiseaux et al., 2006
; Smith et al., 2009
; Zuo et al., 2010b
; see section 3.4). As part of the dual-regression procedure (Biswal et al., 2010
; Filippini et al., 2009
; Zuo et al., 2010b
), all group-level components were first used as a spatial regressor for the participant's functional data, resulting in a timeseries for each component reflecting how closely the functional data spatially resembled the component at each timepoint. Next, these timeseries was together fed into a second regression as a temporal regressor for the functional data (effectively functioning as a complex seed), producing subject-level spatial maps of component connectivity containing both positive and negative values. These subject-level maps were subjected to a final group-level regression analysis examining the effect of motion (see section 2.13).
2.12 ALFF and fALFF
In contrast to seed-based or pairwise approaches, ALFF and fALFF provide a measure of low frequency oscillations that contribute to connectivity but are fundamentally different from between-region measures of connectivity (Zuo et al., 2010a
). ALFF measures the amplitude of low frequency oscillations, while fALFF measures the relative predominance of low frequency amplitude to the amplitude of all oscillations across the entire power spectrum (Zou et al., 2008
).fALFF has greater specificity than ALFF in that it is less susceptible to artifact in the ventricles and near large blood vessels (Zou et al., 2008
; Zuo et al., 2010a
); however, both measures have been found to have high test-retest reliability (Zuo et al., 2010a
). Prior to calculation of ALFF and fALFF, data underwent general preprocessing as described in section 2.5. fALFF was calculated as described by Zuo et al. (2010a)
, and was implemented using scripts written by Biswal et al. (2010)
as part of the 1000 functional connectomes project (www.nitrc.org/projects/fcon_1000/
). Briefly, at each voxel, ALFF is the sum of amplitudes within a specified low-frequency range of oscillations (0.01-0.1 Hz). In contrast, fALFF is calculated as the sum of amplitudes within this range compared to the total sum of amplitudes across the entire frequency range. Prior to group level analyses, all subject-level ALFF and fALFF maps were Z-transformed to improve normality, and registered to a 2mm MNI template as above. In order to provide an estimate of the effect size of motion on both ALFF and fALFF, we averaged ALFF or fALFF signal across all gray matter voxels in template space (using the Harvard-Oxford atlas, thresholded at probability >0.25) and correlated this average cortical ALFF or fALFF value with mean relative displacement. As above, nonlinear effects were investigated using quadratic regression.
2.13 Voxelwise group level analyses of motion
Subject-level ICA dual regression maps and ALFF/fALFF maps were subjected to identical group-level analysis procedures. Group-level analyses were performed in the age/motion-related sample using a GLM where three regressors were included: age, sex, and mean relative displacement. Mean relative displacement was the effect of interest. All group-level voxelwise analyses were conducted using non-parametric permutation testing methods implemented using RANDOMISE (Nichols and Holmes, 2002
). Type I error was controlled with threshold-free cluster enhancement (TFCE); clusters of 100 voxels above a corrected p value <0.01 were considered significant and reported in the text. For display, images were rendered on a cortical map of the right hemisphere using Caret V5.6 (Van Essen et al., 2001
). All figures display uncorrected, unmasked t-statistic maps.
2.14 Relevance of motion to estimation of relationship between age and connectivity
While the above analyses used an age/motion related sample to investigate how in-scanner motion impacted multiple different analytic approaches, the final aim of the paper was to examine how in-scanner motion was related to age, and might impact estimates of how connectivity changes with age. Therefore, the relationship between relative mean displacement and age was assessed in the initial inclusive sample of 456 subjects, and each sub-sample (related and unrelated) using Pearson's correlations. In order to establish the relationship between inclusion threshold and sample size, we plotted the percent of the original sample retained using exclusion thresholds ranging from 0mm to 1mm mean relative displacement. In addition, we examined the effect of exclusion threshold on the correlation between age and motion over the same range of exclusion thresholds.
Next, we examined the potential confounding effects motion might have on studies of the relationship between age and connectivity. Prior investigations have highlighted the association of age with increased long-range connectivity of cortical hubs involved in higher-order mental processes (Fair et al., 2007
; Fair et al., 2008
), and more generally demonstrated that with age short-range connectivity decreases whereas long-range connectivity is enhanced (Dosenbach et al., 2010
).We examined age effects on a voxelwise basis from the seed region in the posterior cingulate cortex described above. We chose to examine the effects of age using this specific analysis as such a voxelwise approach is sensitive to the modulation of the effects of age by distance. The voxelwise approach chosen here may be more prone to type I error, but this was controlled using the same statistical thresholding as used for the other voxelwise analyses (TFCE corrected p<0.01; see section 2.13). Voxelwise PCC time-series analysis was carried out using FILM with local autocorrelation correction as implemented in FEAT (Woolrich et al., 2001
). As in the pairwise seed analyses, we included confound regressors in the GLM in addition to the PCC timecourse (e.g., motion parameters, mean white matter timecourse, mean CSF timecourse, and mean whole-brain timecourse). To investigate motion's impact on the estimate of the effect of age, we conducted three group-level analyses. First, we examined the effect of age without controlling for motion in the age/motion-related
sample (n=421).We compared this to two analyses where motion was accounted for: an analysis of age in the age/motion-unrelated
sample (n=348), and an analysis of age in the age/motion-related
sample with each subject's mean relative displacement included in the group-level regression as a confound variable (“motion regressed
”). Group level statistics were conducted with RANDOMISE using the procedures described above (see section 2.13). For further illustration, we also report correlation and partial correlation values between the PCC seed and the MPFC ROI for each sample. It should be noted that as these sub-samples were overlapping, in this descriptive analysis direct statistical testing between subsamples was not conducted.