Head motion is a significant source of artifacts in resting-state fMRI (rsfMRI) studies and has been shown to affect resting-state functional connectivity (rsFC) measurements. In many rsfMRI studies, motion parameters estimated from volume registration are used to characterize head motion and to mitigate motion artifacts in rsfMRI data. While prior task-based fMRI studies have shown that task-evoked brain activations may induce temporally correlated bias in the motion estimates, resulting in artificial activations after registration, relatively little is known about neural-related bias in rsfMRI motion parameter. In this study, we demonstrate that neural-related bias exists in rsfMRI motion estimates and characterize the potential effects of the bias on rsFC estimates. Using a public multi-echo rsfMRI dataset, we use the differences between motion estimates from the first echo and second echo data as a measure of neural-induced bias. We show that the resting-state global activity of the brain, as characterized with the global signal (GS), induces bias in the motion estimates in the y- and z-translational axes. Furthermore, we demonstrate that the GS-related bias reflects superior-inferior and anterior-posterior asymmetries in the GS beta coefficient map. Finally, we demonstrate that regression with biased motion estimates can negatively bias rsFC estimates and also reduce rsFC differences between young and old subjects.

Head motion is a major source of artifacts in fMRI data, and motion estimates from data registration are commonly used to characterize and mitigate motion-related artifacts (Dijk et al., 2012; Friston et al., 1996; Power et al., 2012, 2014, 2015; Satterthwaite et al., 2012; Yan et al., 2013). Prior task-based fMRI studies have demonstrated that task-evoked brain activations can induce bias in motion estimates that is temporally correlated with brain activations (Freire & Mangin, 2001; Freire et al., 2002). In contrast, relatively little is known about neural-related bias in motion estimates from resting-state fMRI (rsfMRI) data.

Detecting such bias is challenging due to the absence of ground truth head motion and neural signals in rsfMRI studies. To first order, fMRI signal changes resulting from head motion and neural activity can be distinguished based on their dependence on echo time (TE): neural activity causes TE-dependent blood-oxygen-level-dependent (BOLD) fluctuations, whereas head motion largely contributes to TE-independent non-BOLD signal changes. Taking advantage of these differences, prior studies have demonstrated that BOLD and non-BOLD signals can be distinguished using multi-echo fMRI (MEfMRI) that acquires data at different TEs (Bright & Murphy, 2013; Buur et al., 2009; Kundu et al., 2012, 2017). Using short TE first echo data as a regressor to remove motion-related artifacts from data acquired at a longer TE with significant BOLD-weighting, Buur et al. (2009) demonstrated improvements in task-related functional activation maps (computed using the longer TE data) for individual subjects. The primary rationale behind their work was that the short TE data primarily reflected initial signal intensity S0 fluctuations (due to confounds such as motion) that have minimal BOLD weighting. Bright and Murphy (2013) later demonstrated that Burr et al.’s approach could be applied to improve the accuracy of resting-state functional connectivity maps.

Burr’s findings imply that it may be feasible to identify neural-related bias by comparing motion estimates obtained from the first and second echo data. Since head motion primarily results in TE-independent signal changes, its effects should be captured in the motion estimates derived from both the first and second echo data. In contrast, there is typically minimal BOLD weighting in the first echo data but strong BOLD weighting in the second echo data. As prior work has shown that a higher level of brain activation can lead to greater bias in motion estimates (Freire & Mangin, 2001; Freire et al., 2002), potential neural-related bias should exhibit greater magnitude in the motion estimates from the second echo data as compared to those from the first echo data. Taken together, potential neural-related bias may be isolated from head motion by examining the difference between the motion estimates obtained from the first and second echo data.

Prior task-based fMRI studies have shown that brain activations with a larger spatial extent can lead to a higher level of bias in the motion estimates (Freire & Mangin, 2001; Freire et al., 2002). This finding motivates us to examine whether global brain activity in rsfMRI leads to bias in the motion estimates. In this study, we used the global signal (GS), calculated as the mean fMRI signal over the voxels within the brain, as a proxy for global brain activity. While the interpretation of the GS is still controversial (reviewed in T. T. Liu et al., 2017), there is growing evidence suggesting that the GS is linked to global neural activity (Falahpour et al., 2016, 2018; Gu et al., 2021; X. Liu et al., 2018; Wong et al., 2013, 2016), in addition to other components such as low frequency confounds and spin-history effects that are not related to neural activity.

In this work, we provide a mathematical rationale for how the difference in motion estimates can isolate the effects of head motion from GS-related effects. We then characterize the GS-related bias in rsfMRI motion parameters as estimated by AFNI 3dvolreg with a public multi-echo fMRI dataset (Setton et al., 2022; Spreng et al., 2022). To reflect the GS-related effects that would be typically observed in single-echo fMRI studies, we use the second echo data, which were acquired at a commonly used TE (30 ms), for computation of the GS. Furthermore, we investigate the consequences of using biased estimates as regressors in resting-state functional connectivity (rsFC) analyses.

2.1 Subjects and MEfMRI data acquisition

We used a public dataset (denoted as the Cornell-York dataset and described in Setton et al. (2022) and Spreng et al. (2022)) downloaded from OpenNeuro (dataset ds003592). The Cornell-York dataset includes multi-echo fMRI data collected from 301 healthy subjects (181 younger: mean age = 22.59 y, SD = 3.27 y, 57% female; and 120 older adults: mean age = 68.63 y, SD = 6.44 y, 54% female). The data from 238 subjects were acquired on a 3 T GE Discovery MR750 MRI scanner with a 32-channel head coil. The data from the remaining 63 subjects were collected on a 3 T Siemens Trio MRI scanner with a 32-channel head coil. For each subject, two 10-min resting-state runs were acquired using an ME EPI sequence on the GE scanner (204 volumes; TR = 3000 ms; TE = 13.7,30,47 ms; flip angle = 83°; FOV = 210 mm; voxel size = 3 × 3 × 3 mm3; matrix size = 72 × 72 × 46; SENSE acceleration factor = 2.5; phase encoding direction: A-P) or on the Siemens scanner (200 volumes; TR = 3000 ms; TE = 14,29.96,45.92 ms; flip angle = 83°; FOV = 216 mm; voxel size = 3.4 × 3.4 × 3 mm3; matrix size = 64 × 64 × 43; GRAPPA acceleration factor = 3; phase encoding direction: A-P). For both scan protocols, there was no acceleration in the slice direction. The subjects were instructed to stay awake and lie still with their eyes open during the scans.

2.2 Data preprocessing

AFNI and MATLAB were used for data preprocessing (Cox, 1996). As described in detail below and summarized in Supplementary Material Figure S1, we used standard AFNI algorithms for resampling (3dResample), registration (3dVolreg), mask estimation (3dAutomask), and transformation to standard space (align_epi_anat.py, auto_warp.py, and 3dNwarpApply). The rationale behind our approach was to use the minimal processing needed to demonstrate the relationship between the various metrics, such as the global signal and motion estimates. The fMRI data from the first and second echoes were used and denoted as e1 and e2, respectively. The data were first reoriented to Right-Anterior-Inferior (RAI) orientation (AFNI 3dresample -orient RAI), and the first 6 TRs of the data were discarded to allow magnetization to reach a steady state.

2.3 Calculation of the motion estimates and global signal

In this study, we examined the motion estimates from AFNI 3dvolreg (Cox & Jesmanowicz, 1999) using the first volume as the reference volume. The default weights used in registration were disabled by feeding an all-ones image to the -weight option. This approach weights all voxels equally during registration and simplifies the theory. As shown in Supplementary Material (Fig. S2), similar results were obtained when using the default weights. For each run, the motion was estimated separately for the e1 and e2 data, resulting in two sets of estimates. The motion parameters were multiplied by —1 to represent the movement of the volumes as compared to the reference volume (see Appendix A).

For each run, the global signal (GS) was calculated from the registered e2 data. As shown in Supplementary Material (Fig. S2), nearly identical results were obtained when the GS was calculated from the unregistered e2 data. Before calculating the GS, each voxel’s signal was normalized to represent the percent BOLD signal change from the voxel’s mean signal value (computed over the duration of the scan), as the percent signal change has a meaningful connection with brain physiology (T. T. Liu et al., 2017). Then, the GS was computed by averaging over all voxels within the brain (brain masks were formed by AFNI 3dAutomask). Finally, the mean and the linear and quadratic trends were regressed out from the motion estimates and GS, as these are typically considered to represent confounds not related to neural activity (T. T. Liu et al., 2017).

2.4 Identifying BOLD-weighted GS bias in the motion estimates

To characterize BOLD-weighted GS bias in the motion estimates, we considered a simple signal model for the e1 and e2 motion estimates. Let me1K×1 and me2K×1 represent the motion estimates of one motion axis each from e1 and e2, respectively, where K is the number of volumes. We model the motion estimates as the weighted sum of head motion, BOLD-weighted GS bias, and estimation error,

me1=αe1mh+βe1TE1gs+ϵe1me2=αe2mh+βe2TE2gs+ϵe2
(1)

where mhK×1 represents head motion, αei is the regression weight corresponding to head motion for the ith echo, gsK×1 is the GS, TE1 and TE2 are the first and second echo times, respectively, βei is the regression weight of the BOLD-weighted GS bias for the ith echo, and ϵeiK×1 is the estimation error for the ith echo. With the above model, subtracting me1 from me2 yields

Δm=me2me1=(αe2αe1)mh+(βe2TE2βe1TE1)gs+ϵe2ϵe1
(2)

Since prior findings suggest that head motion can be accurately estimated from both e1 and e2 data (Buur et al., 2009; Speck & Hennig, 2001), we assume that αe1αe2, yielding

Δm(βe2TE2βe1TE1)gs+ϵe2ϵe1
(3)

Note that Δm isolates BOLD-weighted GS bias from head motion. Therefore, we can examine the presence of potential bias by assessing the significance of the correlation between Δm and the GS. For each run and motion axis, we calculated the correlation between Δmand the GS (denoted as r(Δm,GS)) and assessed the statistical significance of r(Δm,GS) values on a per-run basis using empirical null distributions. For each motion axis, a null distribution was formed by calculating r(Δm,GS) values using all possible permutations across runs, that is, pairing the GS from one run to Δm from other runs and looping over all runs. The resulting null distributions consisted of 361,802 samples. Then, for each motion axis, we used the null distribution to compute the two-sided p-value associated with the r(Δm,GS) value calculated from each run’s measured data. A p-value threshold of 0.05 was divided by 602 (the number of runs) to correct for multiple comparisons, and the Bonferroni-corrected threshold was used to determine the significant correlations. Finally, we calculated the percentage of runs that show significant r(Δm,GS) values for each motion axis.

In this study, we used the GS to represent the global activity of the brain. To reduce the potential effect of motion artifacts in the GS, we repeated the above analysis after regressing out the e1 motion regressors from both the GS and Δm. The e1 motion regressors included the six motion parameters estimated from the e1 data and their first derivatives.

2.5 Spatial maps underlying r(Δm,GS)

To provide insight into the mechanisms underlying BOLD-weighted GS bias, we derived an empirical approximation for r(Δm,GS). As shown in Appendices B and C,

r(Δm,GS)gsΔm(βgsT,e2de2de22βgsT,e1de1de12)
(4)

where βgs,eiN×1 is the GS beta coefficient map for the ith echo, deiN×1 is the spatial derivative image with respect to (w.r.t.) one motion axis of the ith echo, N is the number of voxels, and denotes the L2-norm. For each run, the GS beta coefficient map βgs,ei=Yeigsgs2 was calculated from the linear fit of the GS to the unregistered functional data YeiN×K of the ith echo.

We calculated the spatial derivative images w.r.t. the motion axes following the algorithm implemented in AFNI 3dvolreg. For each run and each echo, the spatial derivative images were calculated based on the reference volume used in motion estimation. Denoting the reference image of the ith echo as yr,eiN×1, the spatial derivative image of the ith echo w.r.t the jth motion axis was calculated as

dei,j=Tj(yr,ei,ϵj)Tj(yr,ei,ϵj)2ϵj,i=1,2,j{Rx,Ry,Rz,Tx,Ty,Tz}
(5)

where Tj:N×1N×1 is the function that transforms the reference volume along the jth motion axis; ϵj is the transformation parameter; Rx, Ry, and Rz represent x-, y-, and z-rotation, respectively; and Tx, Ty, and Tz represent x-, y-, and z-translation, respectively. When calculating the spatial derivative images, the transformation was performed with AFNI 3drotate -heptic. For the translational motion axes, ϵ was set to 2.1 mm. For the rotational motion axes, ϵ was set to 0.4. These ϵ values were determined by AFNI 3dvolreg based on the spatial resolution of the functional data.

2.6 Effect of the bias on ROI-ROI FC via motion regression

In this work, we investigated the effect of regression with biased motion estimates on rsFC estimates between nodes within the default mode network (DMN) and dorsal attention network (DAN). The regions of interest (ROIs) within these networks were defined in Dijk et al. (2010), with four ROIs in the DMN (posterior cingulate cortex (PCC), lateral parietal cortex (LatPar), medial prefrontal cortex (mPFC), and Hippocampal formation (HF)) and three ROIs in the DAN (frontal eye field (FEF), intraparietal cortex (IPS), and middle temporal area (MT+)). Seed ROIs were created using a sphere with a diameter of 12 mm centered about each seed coordinate (Wong et al., 2013). The left and right ROIs were combined to form bilateral ROIs (See Supplementary Material Fig. S3 for depiction of ROIs). Prior to averaging signals within each ROI, the e2 data after volume registration were transferred to MNI space using AFNI programs align_epi_anat.py, auto_warp.py, and 3dNwarpApply. As supplementary analyses, we also considered (1) transformation to standard space using a custom young-old template (Laurita et al., 2020) and (2) rsFC analysis using a 200-region parcellation of the brain (Schaefer et al., 2018).

To reduce the confounding effects of head motion, motion censoring was performed before motion regression. For each run, framewise displacement (FD) was calculated based on six motion parameters estimated from the e1 data (the calculation of FD follows the description in Power et al. (2012)). Volumes with FD values larger than 0.2 mm were censored. We show in Supplementary Material (Figs. S7 and S8) that motion censoring has little impact on the effect of regression with biased motion estimates on rsFC.

After motion censoring, we calculated and compared the ROI-ROI FC after e2 and e1 motion regression to assess the effect of potential bias on FC analysis. For each ROI, the ROI-based seed signal was calculated by averaging the percent change BOLD signal over the voxels in that ROI. Then, the motion regressors, including the six motion parameters and their first derivatives, were regressed out from the ROI signals. For each pair of ROIs, the ROI-ROI FC was computed as the Pearson’s correlation coefficient between the ROI average signals. Correlation values were converted to z-scores using the Fisher-z transformation. Note that all statistical comparisons were performed on the z-scores. In some plots, the r-values are also presented to support the interpretation of the results. The differences in the ROI-ROI FC calculated after e2 and e1 motion regression were calculated to assess the effect of the bias on FC. The e2-e1 differences in r-values and z-scores were denoted as Δr and Δz, respectively.

For each run, the level of head motion was measured by the mean FD value calculated by averaging over the FD values across time. The level of the GS was measured by the GS amplitude (aGS) computed as the standard deviation of the GS after motion censoring. We used 2D scaled histograms to examine the distribution of the Δz values (averaged over all ROI pairs) as a function of aGS and motion (mean FD). Based on this examination, we then arranged the data into four groups based on aGS and mean FD levels and characterized the mean Δz values for each ROI pair averaged over runs within each group. All the runs were first divided into two groups based on their mean FD values. Runs with mean FD values larger than the group median were classified as high motion runs, and the remaining runs were classified as low motion runs. Then, within each FD-based group, the runs were further divided into two groups based on aGS. Runs with aGS values larger than the group median were classified as high aGS runs, and the remaining runs were classified as low aGS runs. Consequently, we formed four groups of runs: (1) low motion and high aGS runs, (2) low motion and low aGS runs, (3) high motion and high aGS runs, and (4) high motion and low aGS runs. A one-way ANOVA was calculated on mean Δz over all ROI pairs to assess whether there was a group effect. Post-hoc two-sample t-tests were calculated to characterize the differences between pairs of groups.

Additionally, to verify that the effect of the bias on FC was dominated by the motion estimates from the axes where we found GS-induced bias, we evaluated the effect of two subsets of motion regressors. One of the subsets included the motion estimates and their first derivatives from the Ty and Tz axes (where we found GS-induced bias), while the other subset included the motion regressors from the other four motion axes, where minimal GS-induced was observed.

2.7 Effect of the bias on young vs. old group-level FC analysis

We first examined whether the bias in motion estimates affects the ROI-ROI FC of the young and old subjects differently by comparing Δr and Δz between the young and old runs. Furthermore, we investigated if the bias alters the group-level FC analysis between the young and old subjects. For each pair of ROIs, the significance of the FC differences between the young and old subjects was assessed by a permutation test with 1×107 random permutations to allow us to apply p-value thresholds of 0.01, 1×103 and 1×106. Also, the effect size of the differences was measured with Cohen’s d. The FC differences and the significance and effect size of the differences calculated after e1 and e2 motion regression were compared.

3.1 Examples of the GS and motion estimates

As discussed below, we found a significant association between the GS and Δm in the Tz and Ty axes. To provide a qualitative view, Figure 1 (a, b) shows motion estimates in the Tz axis, including me1 (blue), me2 (green), and Δm (red) from two example runs with (a) low and (b) high levels of motion. Note that for the high motion run, me1 and me2 are plotted at 1/20th scale to facilitate comparison with me1 and me2 from the low motion run. As shown in these subfigures, Δm estimates from both runs fluctuate in a similar range from 0.05 to 0.05 mm (std(Δm)=0.017 and 0.025 for the low and high motion runs, respectively). In contrast, the standard deviations of me1 and me2 from the high motion run (std(me1)=0.583, std(me2)=0.602) are an order of magnitude larger than the standard deviations of me1 and me2 from the low motion run (std(me1)=0.027, std(me2)=0.034).

Fig. 1.

Motion estimates in the Tz axis, including me1 (blue), me2 (green), and Δm (red) from (a) low and (b) high motion runs. Panels (c) and (d) show the GS (blue) and Δm (red) from the low and high motion runs, respectively. Panels (e) and (f) show the GS (blue) and Δm (red) after motion regression (denoted as MotReg; with me1 as regressor) from the low and high motion runs, respectively.

Fig. 1.

Motion estimates in the Tz axis, including me1 (blue), me2 (green), and Δm (red) from (a) low and (b) high motion runs. Panels (c) and (d) show the GS (blue) and Δm (red) from the low and high motion runs, respectively. Panels (e) and (f) show the GS (blue) and Δm (red) after motion regression (denoted as MotReg; with me1 as regressor) from the low and high motion runs, respectively.

Close modal

Figure 1 (c) and (d) show the GS and Δm from the low motion and high motion runs, respectively. For the low motion run, Δm covaries with the GS throughout the run, leading to a strong r(Δm,GS) of 0.93. The high motion run shows a weaker r(Δm,GS) of 0.54 as compared to the low motion run, and reflects motion artifacts in the GS. After motion regression (with me1; panels e and f), the r(Δm,GS) for the high motion run increases to 0.62, while the r(Δm,GS) of the low motion run remains at a high value of 0.92, suggesting that the relation between the GS and Δm is enhanced when motion artifacts are minimized.

3.2 Significance testing for r(Δm,GS) values

To examine the presence of BOLD-weighted GS bias over runs and motion axes, we assessed the significance of r(Δm,GS) values on a per-run and per-axis basis using permutation-based empirical null distributions. Figure 2 shows two-sided violin plots of the distributions of r(Δm,GS) values (blue) and the empirical null distributions (green) for all six motion axes (a) before and (b) after e1 motion regression. The blue solid lines and circles represent the median values for each distribution of measured r(Δm,GS) values. The dashed lines represent the r(Δm,GS) values corresponding to a Bonferroni-corrected p-value threshold of 0.05 (two-sided) assessed from the empirical null distributions. The dark red square markers represent the percentages of runs showing r(Δm,GS) values that are significantly different from zero.

Fig. 2.

Two-sided violin plots showing the distributions of r(Δm,GS) values (blue) and the empirical null distributions (green) for all six motion axes (a) before and (b) after e1 motion regression. The blue solid lines and circles represent the median values for each data distribution. For each motion axis, a Bonferroni-corrected run-wise p-value threshold of 0.05/602 was used, where 602 is the number of runs. The black dashed lines show r(Δm,GS) values that correspond to the p-value thresholds (two-sided) assessed from the empirical null distributions. The dark red dashed lines and square markers represent the percent of the runs with significant r(Δm,GS) values.

Fig. 2.

Two-sided violin plots showing the distributions of r(Δm,GS) values (blue) and the empirical null distributions (green) for all six motion axes (a) before and (b) after e1 motion regression. The blue solid lines and circles represent the median values for each data distribution. For each motion axis, a Bonferroni-corrected run-wise p-value threshold of 0.05/602 was used, where 602 is the number of runs. The black dashed lines show r(Δm,GS) values that correspond to the p-value thresholds (two-sided) assessed from the empirical null distributions. The dark red dashed lines and square markers represent the percent of the runs with significant r(Δm,GS) values.

Close modal

In the Tz axis, 63.8% and 69.3% of the runs show significant positive r(Δm,GS) values before and after e1 motion regression, respectively. The group median r(Δm,GS) value increases from 0.57 to 0.59 after e1 motion regression, with 97 out of 602 total runs showing r(Δm,GS) values larger than 0.8 after e1 motion regression. In the Ty axis, 30.7% and 40.0% of the runs show significant r(Δm,GS) values before and after e1 motion regression, respectively. The group median r(Δm,GS) value increases from 0.31 to 0.35 after e1 motion regression. Figures S4 and S5 show examples of the GS and the motion estimates, including me1, me2, and Δm in the Tz and Ty axes, respectively. Together, these findings indicate the presence of BOLD-weighted GS bias in the Tz and Ty motion estimates.

For the other motion axes (Rx, Ry, Rz, and Tx), the percent of significant r(Δm,GS) values fluctuates around 5%, ranging from 1.3% to 7.1%, indicating minimal BOLD-weighted GS bias for these axes. As described in Supplementary Material, a preliminary examination shows that bias is also observed when using other software packages, such as SPM (Friston et al., 1995), FSL (Jenkinson & Smith, 2001; Jenkinson et al., 2002), and ANTS (Avants et al., 2011, 2014). The bias observed with SPM is similar to that observed with AFNI, while FSL and ANTS show relatively lower levels of bias.

3.3 Visualizing the spatial maps underlying r(Δm,GS) values

In the previous section, BOLD-weighted GS bias in the motion estimates was identified by examining the temporal correlations between the GS and Δm. Furthermore, as described in Eq. 4, r(Δm,GS) can be approximated as the product of the difference of βgsTdd2 values from e2 and e1 and a scaling factor gsΔm. The approximation provides a unique angle to interpret the GS bias by looking at the relation of the βgs and d spatial maps. Figure 3 visualizes βgs, dd2, βgsdd2, and the βgsdd2 difference maps in the Tz and Tx axes from an example run, where represents element-wise multiplication. Note that larger magnitudes in the GS beta coefficient βgs maps reflect a greater presence of the GS in a voxel’s time series; larger magnitudes in the normalized spatial derivative dd2 maps indicate voxels with relatively higher spatial derivatives; larger magnitudes in the element-wise multiplication βgsdd2 maps indicate voxels for which both the GS beta coefficients and normalized spatial derivatives have relatively high values; and larger magnitudes in the element-wise difference maps typically indicate voxels for which the values in the element-wise map for the second echo βgs,e2de2de22 are greater than the corresponding values in the element-wise map for the first echo βgs,e1de1de12, where the differences are driven for the most part by the greater βgs values observed for the second echo.

Fig. 3.

Visualization of the spatial maps underlying r(Δm,GS) in the Tz (rows 2-6 on the left) and Tx (rows 2-6 on the right) axes from an example run. For each map, three representative slices (one axial, one sagittal, and one coronal) are plotted. From top to bottom, rows 1 through 3 show the βgs, dd2, and βgsdd2 maps from e1 and e2, where represents element-wise multiplication. Row 4 shows βgsdd2 difference (e2–e1) maps. From left to right in rows 2 and 3, the first and second columns show the e1 and e2 maps in the Tz axis, and the third and fourth columns show the e1 and e2 maps in the Tx axis. The red and blue arrows point to brain regions showing high positive and negative values in the maps, respectively.

Fig. 3.

Visualization of the spatial maps underlying r(Δm,GS) in the Tz (rows 2-6 on the left) and Tx (rows 2-6 on the right) axes from an example run. For each map, three representative slices (one axial, one sagittal, and one coronal) are plotted. From top to bottom, rows 1 through 3 show the βgs, dd2, and βgsdd2 maps from e1 and e2, where represents element-wise multiplication. Row 4 shows βgsdd2 difference (e2–e1) maps. From left to right in rows 2 and 3, the first and second columns show the e1 and e2 maps in the Tz axis, and the third and fourth columns show the e1 and e2 maps in the Tx axis. The red and blue arrows point to brain regions showing high positive and negative values in the maps, respectively.

Close modal

In the Tz axis, dd2 exhibits high positive and negative values at the superior and inferior edges of the brain, respectively, whereas βgs exhibits positive values across the cortex, including the superior edge of the brain, but exhibits relatively low values (approaching zero) along the inferior edge. As a result, most of the high values in βgsdd2 are positive with greater amplitudes for e2, resulting in high positive values in the βgsdd2 difference map. Consequently, summing up the βgsdd2 difference map leads to a high r(Δm,GS) value. In contrast, for the Tx axis, dd2 shows high positive and negative values on the left and right edges, respectively, whereas βgs shows high positive values across the cortex, including both the left and right edges. The resulting high positive and negative values observed in the βgsdd2 difference map tend to cancel out, resulting in a lower r(Δm,GS) value. Together, these observations suggest that the GS bias in the Tz motion estimates results from βgs and d sharing a relatively similar top-bottom asymmetric spatial pattern.

In the Ty axis, as shown in Figure S6, dd2 exhibits negative and positive values at the anterior and posterior edges of the brain, respectively, and βgs exhibits positive values across the cortex, including the posterior edge of the brain, but exhibits relatively low values (approaching zero) along the anterior edge. As a result, the positive values in the βgsdd2 difference map show relatively greater amplitudes as compared to the negative values, resulting in a high r(Δm,GS) value.

Taken together, the GS-induced bias in the Ty and Tz motion estimates arises because of the interaction between the roughly odd symmetry in the respective spatial derivative maps, which exhibit high positive and negative values at the edges, and the asymmetries in the respective βgs maps. If the βgs maps exhibited even symmetric values at the edges, then the odd symmetry in the spatial derivative maps would be largely preserved when multiplying by the βgs maps, such that the high positive and negative contributions would cancel out when computing the sum over the βgsdd2 difference maps. However, because the βgs maps show higher values at the superior and posterior parts of the brain as compared to the inferior and anterior parts, the odd symmetry in the spatial derivative maps is not preserved and the high positive and negative contributions do not cancel out when computing the sum of the βgsdd2 difference maps.

As shown in Supplementary Material, the superior-inferior and posterior-anterior asymmetric spatial patterns in the βgs maps show a significant spatial correlation with the patterns observed in the vigilance template, suggesting that the asymmetric patterns may partly reflect spatial variations in the effect of vigilance variations on the fMRI signal. Furthermore, as described in Supplementary Material, image intensity variations due to factors such as inhomogeneities in the receive coil sensitivity patterns (Belaroussi et al., 2006) are another potential source of the posterior-anterior asymmetric pattern in βgs.

3.4 Effect of the bias on ROI-ROI FC via motion regression

In this section, we investigate how the presence of the GS-induced bias in the motion estimates affects the ROI-ROI FC by examining the differences in the ROI-ROI FC calculated after e2 and e1 motion regression. The e2–e1 differences in r-values and z-scores are denoted as Δr and Δz, respectively.

Figure 4 shows the distribution of mean Δz values (averaged over all ROI pairs) as a function of aGS and mean FD. When all runs are considered together (panel a), the most negative (dark blue) Δz values occur for runs with high aGS and low mean FD. This behavior is largely driven by the distribution of values for the younger subjects (panel b), which is in contrast to the distribution of values for the older subjects (panel c), where there are relatively few samples at high aGS and low mean FD. The age-related difference in the distributions gives rise to the group analysis effects addressed in Section 3.4.1. As shown in Figure S9, similar distributions are observed when considering only those runs acquired on the GE scanner.

Fig. 4.

2D scaled histograms showing differences (e2–e1) in mean ROI-ROI FC (z-scores) as a function of aGS and mean FD (mm), where the heights of the bars indicate the count and the color indicates the mean Δz values.

Fig. 4.

2D scaled histograms showing differences (e2–e1) in mean ROI-ROI FC (z-scores) as a function of aGS and mean FD (mm), where the heights of the bars indicate the count and the color indicates the mean Δz values.

Close modal

Figure 5 shows Δr and Δz for all pairs of ROIs averaged over four groups of runs: (a) low motion and high aGS runs, (b) low motion and low aGS runs, (c) high motion and high aGS runs, and (d) high motion and low aGS runs. All four groups demonstrate negative Δr and Δz values for all pairs of ROIs, suggesting that the bias in the motion estimates can reduce FC estimates via motion regression. Comparing groups, we observed that the runs with low motion and high aGS exhibit the most negative Δr and Δz values (Δr ranging from -0.02 to -0.07 and Δz ranging from -0.50 to -1.10), while the runs with high motion and low aGS show Δr and Δz values close to zero (Δr ranging from -0.00 to -0.02 and Δz ranging from -0.02 to -0.26).

Fig. 5.

Differences (e2–e1) in ROI-ROI FC calculated after e2 and e1 motion regression. The differences were averaged across runs within each of four groups: (a) low motion and high aGS runs, (b) low motion and low aGS runs, (c) high motion and high aGS runs, and (d) high motion and low aGS runs. Each subplot is divided into an upper right triangle showing the average differences in r-values and a lower left triangle showing the average differences in z-scores.

Fig. 5.

Differences (e2–e1) in ROI-ROI FC calculated after e2 and e1 motion regression. The differences were averaged across runs within each of four groups: (a) low motion and high aGS runs, (b) low motion and low aGS runs, (c) high motion and high aGS runs, and (d) high motion and low aGS runs. Each subplot is divided into an upper right triangle showing the average differences in r-values and a lower left triangle showing the average differences in z-scores.

Close modal

Furthermore, a one-way ANOVA on the mean Δz values calculated over ROI pairs showed that there is a significant (p<1×106, F3,598=46.27) difference among groups. The results from the post-hoc t-tests indicate that the low motion and high aGS runs exhibit significantly (p<1×104) more negative mean Δz values as compared to the other three groups, whereas the high motion and low aGS runs exhibit significantly (p<1×103) less negative mean Δz values as compared to the other three groups. Also, the group of low motion and low aGS runs shows a significantly (p=6.5×106) less negative mean Δz value as compared to the group of high motion and high aGS runs. The distributions of the mean Δz values over ROI pairs for the four groups are shown in Figure S8.

Figures S10 and S11 demonstrate that the observed effect of the bias is dominated by the motion estimates in the Ty and Tz axes. In Figure S10, when including only the Ty and Tz motion regressors, the results are similar to those shown in Figure 5. In contrast, when excluding the Ty and Tz motion regressors (Fig. S11), we observed minimal differences between the FC estimates obtained with e1 and e2 motion regression. Together, these results indicate that regressing out the motion estimates with GS-induced bias may lead to reductions in FC estimates, with a stronger effect for runs with higher aGS and lower head motion.

3.4.1 Effect of the bias on young vs. old group-level analysis

As shown in Figure 6 (a) and (b), young subjects show significantly (p<1×106) higher aGS and lower mean FD values as compared to the old subjects, with the higher aGS values in the young consistent with the higher correlation values seen in the young GS topography maps presented in Supplementary Material Figure S12. Consequently, as shown in Figure 6 (c), the young subjects show significantly (p<1×106) more negative mean Δz (e2-e1) values over ROI pairs as compared to the old subjects. The mean Δz values over runs and ROI pairs are -0.53 and -0.11 for the young and old subjects, respectively. Furthermore, Figure 6 (d) and (e) visualize the mean Δr and Δz values over the young and old subjects for all ROI pairs, with the underlying z-score and r-value maps shown in Supplementary Material Figure S13 (in addition, maps with global signal regression are shown in Fig. S14). The Δz values were only weakly correlated with the z-score differences for e1 motion regression versus no motion regression, accounting for 6.8% and 2.0% of the variance in the young and old groups, respectively. We observed that for all ROI pairs, the young subjects show more negative Δr and Δz values as compared to the old subjects. Since Δz and Δr reflect the effect of the bias when performing motion regression, these results suggest a greater impact of the bias on the young subjects as compared to the old subjects. Similar age-dependent reductions in Δz and Δr are observed across multiple brain networks, as shown in Supplementary Material Figure S15 for a 200-region parcellation of the brain.

Fig. 6.

(a-c) Two-sided violin plots showing the distributions of (a) aGS, (b) mean FD, and (c) mean Δz (e2-e1) over ROI pairs for the young (green) and old (purple) subjects. Two-sided permutation tests were calculated to assess the significance thresholds for the differences between the young and old subjects. Young subjects show significantly (p<1×106) larger aGS, smaller mean FD, and more negative mean Δz as compared to the old subjects. Panels (d-e) show the mean Δz and Δr over (d) young and (e) old subjects for all ROI pairs. Each subplot is divided into an upper right triangle showing the average differences in r-values and a lower left triangle showing the average differences in z-scores.

Fig. 6.

(a-c) Two-sided violin plots showing the distributions of (a) aGS, (b) mean FD, and (c) mean Δz (e2-e1) over ROI pairs for the young (green) and old (purple) subjects. Two-sided permutation tests were calculated to assess the significance thresholds for the differences between the young and old subjects. Young subjects show significantly (p<1×106) larger aGS, smaller mean FD, and more negative mean Δz as compared to the old subjects. Panels (d-e) show the mean Δz and Δr over (d) young and (e) old subjects for all ROI pairs. Each subplot is divided into an upper right triangle showing the average differences in r-values and a lower left triangle showing the average differences in z-scores.

Close modal

Moreover, we examined whether the bias affects the group-level analysis between the young and the old subjects. Figure 7 (a) shows the differences in ROI-ROI FC between the young and old subjects calculated after motion censoring and e1 motion regression. We observed that the young subjects show significantly (p<0.01) higher FC than the old subjects for all ROI pairs, except for the pairs between LatPar and the DAN ROIs (i.e., FEF, IPS, and MT+). However, after e2 motion regression, as shown in Figure 7 (b), the significant differences for four pairs of ROIs become insignificant (p>0.01), including the pairs between (1) PCC and FEF, (2) PCC and IPS, (3) PCC and MT+, and (4) mPFC and MT+. Additionally, Figure 7 (c) and (d) show that e2 motion regression reduces the FC differences between young and old for all pairs of ROIs as compared to e1 motion regression. A supplementary analysis using a custom young-old template for transformation to standard space (Laurita et al., 2020) yielded similar results, as shown in Figure S16. Furthermore, as shown in Figures S17 and S18, we verified that the observed effect is dominated by regressing out the Ty and Tz motion regressors. In addition, a supplementary analysis using an expanded model with FD as a confound yielded similar reductions in young-old FC differences with e2 motion regression, as shown in Figure S19. Together, these results demonstrate that motion regression with the e2 motion estimates can lead to a reduction in FC differences between the young and old subjects as compared to regression with the e1 motion estimates. Because both the e1 and e2 motion estimates carry information related to actual motion, the differences most likely reflect the effect of additional information carried in the e2 motion estimates.

Fig. 7.

(a-b) ROI-ROI FC differences (young-old) between young and old subjects calculated after (a) e1 and (b) e2 motion regression. Each subplot is divided into an upper right triangle showing the averaged differences in z-scores between young and old subjects and a lower left triangle showing the effect size of the differences. The black asterisks indicate the statistical significance of the differences in z-scores assessed by permutation tests (*: p<0.01, **: p<1×103, ***: p<1×106). A positive value (red color) indicates that the young subjects show higher connectivity as compared to the old subjects. (c) Scatter plot comparing the differences in ROI-ROI FC between the young and the old subjects calculated after e1 and e2 motion regression (denoted as e1 MotReg and e2 MotReg, respectively) for all pairs of ROIs. (d) Scatter plot comparing the effect sizes of the differences in ROI-ROI FC between the young and the old subjects calculated after e1 and e2 motion regression for all pairs of ROIs.

Fig. 7.

(a-b) ROI-ROI FC differences (young-old) between young and old subjects calculated after (a) e1 and (b) e2 motion regression. Each subplot is divided into an upper right triangle showing the averaged differences in z-scores between young and old subjects and a lower left triangle showing the effect size of the differences. The black asterisks indicate the statistical significance of the differences in z-scores assessed by permutation tests (*: p<0.01, **: p<1×103, ***: p<1×106). A positive value (red color) indicates that the young subjects show higher connectivity as compared to the old subjects. (c) Scatter plot comparing the differences in ROI-ROI FC between the young and the old subjects calculated after e1 and e2 motion regression (denoted as e1 MotReg and e2 MotReg, respectively) for all pairs of ROIs. (d) Scatter plot comparing the effect sizes of the differences in ROI-ROI FC between the young and the old subjects calculated after e1 and e2 motion regression for all pairs of ROIs.

Close modal

In this study, we used a public resting-state MEfMRI dataset to examine whether global brain activity can lead to bias in rsfMRI motion estimates and to characterize the potential impact on rsFC estimates. By examining the correlation between the GS and the difference in the motion estimates from the first and second echoes, we found evidence for GS-related bias in the Tz and Ty motion estimates. We also demonstrated that the GS-induced bias can lead to underestimation of rsFC estimates when using motion regression, with low motion and high aGS runs exhibiting the greatest reductions in FC due to the bias and high motion and low aGS runs showing minimal reductions in FC. Finally, we showed that regression with biased motion estimates can reduce rsFC differences between groups of young and old subjects, due in part to different levels of aGS and head motion between the groups.

Extending prior studies that examined bias in task-based fMRI motion estimates, we identified BOLD-weighted bias in rsfMRI motion estimates and investigated its effect on rsFC. Moreover, utilizing multi-echo fMRI data, we proposed a novel method to detect the BOLD-weighted bias in real motion estimates over a large sample of runs with rigorous statistical tests. Furthermore, to investigate the underlying cause of the observed bias, we proposed an empirical approximation to r(Δm,GS) and used this approximation to demonstrate how the presence of the GS-induced bias in the motion estimates may be attributed to the superior-inferior and posterior-anterior asymmetric spatial patterns in the GS beta coefficient map βgs. In prior work (Freire & Mangin, 2001; Freire et al., 2002), the existence of the bias in task-based fMRI was primarily illustrated using simulations and the use of experimental data was limited to one scan.

While the interpretation of the GS is still controversial (reviewed in T. T. Liu et al., 2017), there is growing evidence suggesting that the GS is linked to vigilance (also known as arousal level) (Falahpour et al., 2016, 2018; X. Liu et al., 2018; Wong et al., 2013, 2016) and that arousal-related GS peaks are related to changes in head motion parameters (Gu et al., 2020). For example, Falahpour et al. (2016) found that the GS is negatively correlated with EEG measures of vigilance. Furthermore, studies have shown that a vigilance template, calculated from voxel-wise correlations between the EEG vigilance measures and the fMRI signal, can be used to estimate vigilance fluctuations in fMRI scans (Chang et al., 2016; Falahpour et al., 2018; Goodale et al., 2021). As noted in Section 3.3, the presence of a significant spatial correlation between the vigilance template and the βgs maps suggests that the asymmetric spatial patterns in the βgs maps may partly reflect vigilance effects.

In this work, we performed a detailed analysis of the bias in motion estimates obtained with AFNI 3dvolreg, which is widely used for preprocessing of fMRI data. As noted in Section 3.2, a preliminary examination indicates that lower levels of bias may be achieved with other software packages such as FSL and ANTS. Future work is needed to thoroughly examine how the choice of algorithm (and associated cost functions) affects the level of the bias in the motion estimates and to determine whether migrating to alternate algorithms constitutes an effective mitigation strategy. For example, such an effort could examine performance across the range of cost functions in the AFNI 3dAllineate program. In addition, future work focused on developing registration algorithms with reduced bias would be of interest.

In this study, we demonstrated that regression with biased motion estimates can reduce differences in rsFC between the young and old subjects, as assessed with ROIs in the DMN and DAN. Importantly, we showed that the significant rsFC differences for four pairs of ROIs became insignificant (p>0.01) due to regression with the biased motion estimates, with all four pairs corresponding to connections between the DMN and DAN. These findings suggest that regression with biased motion estimates may impede the detection of rsFC differences between the young and old groups. Because rsFC of the DMN plays a key role in understanding the aging brain (reviewed in Ferreira & Busatto (2013)), it can be important for future aging studies to minimize the effect of the bias on rsFC analyses. Based on our results, investigators who are analyzing single-echo fMRI studies may want to compare findings obtained with and without Ty and Tz motion regressors. In addition, as discussed below, the impact of the biased motion estimates will be lessened when other GS-related regressors are used in the processing. When comparing results from existing aging studies with different processing methodologies, the potential effect of the bias also needs to be considered. Concerning this matter, future work that further examines the potential effect of the bias on the young versus old differences in rsFC measures would be of interest.

Our findings revealed that the effect of the bias on the young versus old rsFC analysis may be caused by different levels of aGS and head motion between the groups, with the old subjects showing a lower level of aGS and a higher level of head motion than the young subjects. The higher levels of motion in the old group have been consistently reported in prior rsfMRI studies (Hausman et al., 2022; Saccà et al., 2021), which reflects declines in executive functioning with aging (Hausman et al., 2022). A direct comparison of aGS between young and old groups does not appear to have been addressed in prior studies. However, previous studies have reported that old subjects show lower rsFC and BOLD variability (measured as the standard deviation of BOLD timeseries) in large-scale brain networks, including the DMN, compared to young subjects (Ferreira & Busatto, 2013; Grady & Garrett, 2018; Kumral et al., 2020; Nomi et al., 2017), supporting the observed lower aGS in the old subjects.

In conclusion, we found that resting-state global brain activity can lead to bias in Ty and Tz motion estimates obtained with a widely used motion correction algorithm. Furthermore, the bias in the motion estimates can lead to reductions in the rsFC measures obtained after motion regression, with an increasing level of bias for runs showing higher global signal amplitude and smaller head motion levels. This GS-related decrease in rsFC values is similar to the reduction seen with global signal regression (T. T. Liu et al., 2017), and therefore concerns about the effects of global signal regression may also apply when regressing out motion parameters estimated from rsfMRI data acquired at typical echo times (e.g., 30 ms). Moreover, our results show that regression with biased motion estimates can reduce group-level rsFC differences between young and old subjects. Similar effects may also be present in other rsfMRI studies in which the groups exhibit different mean levels of global signal amplitude. For rsFC studies in which the processing includes additional GS-related nuisance regressors, such as white matter and cerebrospinal fluid signals and their temporal derivatives, the potential bias due to motion regression is likely to be shared among the various GS-related regressors, especially when partial voluming effects are not adequately minimized (Behzadi et al., 2007; T. T. Liu et al., 2017). As a result when using motion regression together with either global signal regression or approaches such as CompCor (Behzadi et al., 2007), the unique contribution of the bias in the motion estimates is likely to be attenuated. On the other hand, for studies that make use of rsfMRI global activity to estimate vigilance levels (T. T. Liu & Falahpour, 2020), the potential bias in the motion estimates must be taken into account as it can affect the estimation of the metrics. Some studies, such as multi-echo fMRI studies that do not employ motion regression and instead use ICA-based approaches to remove motion-related artifacts (Setton et al., 2022), are unlikely to be affected by the bias. Overall, our results suggest that a greater degree of caution should be used when interpreting rsfMRI metrics obtained when motion regression is used in the processing of the data.

The study that collected the open-source dataset was approved by the Institutional Review Board at Cornell University and the Research Ethics Board at York University. Written informed consent was obtained from all study participants.

The dataset can be downloaded from openNeuro (dataset ds003592). Analysis code and files to generate the figures and results presented in this paper will be made available upon publication through the Open Science Framework DOI 10.17605/OSF.IO/SH79V.

Y.M. and T.T.L. designed the experiments. Y.M. carried out the experiments and prepared the results. Y.M. and T.T.L. wrote and edited the main manuscript text. All authors discussed the results and provided interpretations of the results. All authors reviewed and provided feedback on the manuscript. All authors contributed to the article and approved the submitted version.

The authors have declared no competing interests.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00151.

Avants
,
B. B.
,
Tustison
,
N. J.
,
Song
,
G.
,
Cook
,
P. A.
,
Klein
,
A.
, &
Gee
,
J. C.
(
2011
).
A reproducible evaluation of ANTs similarity metric performance in brain image registration
.
NeuroImage
,
54
,
2033
2044
. https://doi.org/10.1016/j.neuroimage.2010.09.025
Avants
,
B. B.
,
Tustison
,
N. J.
,
Stauffer
,
M.
,
Song
,
G.
,
Wu
,
B.
, &
Gee
,
J. C.
(
2014
).
The Insight ToolKit image registration framework
.
Frontiers in Neuroinformatics
,
8
,
44
. https://doi.org/10.3389/fninf.2014.00044
Behzadi
,
Y.
,
Restom
,
K.
,
Liau
,
J.
, &
Liu
,
T. T.
(
2007
).
A component based noise correction method (Comp-Cor) for BOLD and perfusion based fMRI
.
NeuroImage
,
37
,
90
101
. https://doi.org/10.1016/j.neuroimage.2007.04.042
Belaroussi
,
B.
,
Milles
,
J.
,
Carme
,
S.
,
Zhu
,
Y. M.
, &
Benoit-Cattin
,
H.
(
2006
).
Intensity non-uniformity correction in MRI: Existing methods and their validation
.
Medical Image Analysis
,
10
,
234
246
. https://doi.org/10.1016/j.media.2005.09.004
Bright
,
M. G.
, &
Murphy
,
K.
(
2013
).
Removing motion and physiological artifacts from intrinsic BOLD fluctuations using short echo data
.
NeuroImage
,
64
,
526
537
. https://doi.org/10.1016/j.neuroimage.2012.09.043
Buur
,
P. F.
,
Poser
,
B. A.
, &
Norris
,
D. G.
(
2009
).
A dual echo approach to removing motion artefacts in fMRI time series
.
NMR in Biomedicine
,
22
,
551
560
. https://doi.org/10.1002/nbm.1371
Chang
,
C.
,
Leopold
,
D. A.
,
Schölvinck
,
M. L.
,
Mandelkow
,
H.
,
Picchioni
,
D.
,
Liu
,
X.
,
Ye
,
F. Q.
,
Turchi
,
J. N.
, &
Duyn
,
J. H.
(
2016
).
Tracking brain arousal fluctuations with fMRI
.
Proceedings of the National Academy of Sciences
,
113
,
4518
4523
. https://doi.org/10.1073/pnas.1520613113
Cox
,
R. W.
(
1996
).
AFNI: Software for analysis and visualization of functional magnetic resonance neuroimages
.
Computers and Biomedical Research
,
29
,
162
173
. https://doi.org/10.1006/cbmr.1996.0014
Cox
,
R. W.
, &
Jesmanowicz
,
A.
(
1999
).
Real-time 3D image registration for functional MRI
.
Magnetic Resonance in Medicine
,
42
,
1014
1018
. https://doi.org/10.1002/(sici)1522-2594(199912)42:6<1014::aid-mrm4>3.0.co;2-f
Dijk
,
K. R. A. V.
,
Hedden
,
T.
,
Venkataraman
,
A.
,
Evans
,
K. C.
,
Lazar
,
S. W.
, &
Buckner
,
R. L.
(
2010
).
Intrinsic functional connectivity as a tool for human connectomics: Theory, properties, and optimization
.
Journal of Neurophysiology
,
103
,
297
321
. https://doi.org/10.1152/jn.00783.2009
Dijk
,
K. R. A. V.
,
Sabuncu
,
M. R.
, &
Buckner
,
R. L.
(
2012
).
The influence of head motion on intrinsic functional connectivity MRI
.
NeuroImage
,
59
,
431
438
. https://doi.org/10.1016/j.neuroimage.2011.07.044
Falahpour
,
M.
,
Chang
,
C.
,
Wong
,
C. W.
, &
Liu
,
T. T.
(
2018
).
Template-based prediction of vigilance fluctuations in resting-state fMRI
.
NeuroImage
,
174
,
317
327
. https://doi.org/10.1016/j.neuroimage.2018.03.012
Falahpour
,
M.
,
Wong
,
C. W.
, &
Liu
,
T. T.
(
2016
).
The resting state fMRI global signal is negatively correlated with time-varying EEG vigilance
. In
Proceedings of the 24th Annual Meeting of the ISMRM
(p.
641
). https://cds.ismrm.org/protected/16MProceedings/PDFfiles/0641.html
Ferreira
,
L. K.
, &
Busatto
,
G. F.
(
2013
).
Resting-state functional connectivity in normal brain aging
.
Neuroscience & Biobehavioral Reviews
,
37
,
384
400
. https://doi.org/10.1016/j.neubiorev.2013.01.017
Freire
,
L.
, &
Mangin
,
J. F.
(
2001
).
Motion correction algorithms may create spurious brain activations in the absence of subject motion
.
NeuroImage
,
14
,
709
722
. https://doi.org/10.1006/nimg.2001.0869
Freire
,
L.
,
Roche
,
A.
, &
Mangin
,
J. F.
(
2002
).
What is the best similarity measure for motion correction in fMRI time series
?
IEEE Transactions on Medical Imaging
,
21
,
470
. https://doi.org/10.1109/tmi.2002.1009383
Friston
,
K. J.
,
Ashburner
,
J.
,
Frith
,
C. D.
,
Poline
,
J.
,
Heather
,
J. D.
, &
Frackowiak
,
R. S. J.
(
1995
).
Spatial 560 registration and normalization of images
.
Human Brain Mapping
,
3
,
165
189
. https://doi.org/10.1002/hbm.460030303
Friston
,
K. J.
,
Williams
,
S.
,
Howard
,
R.
,
Frackowiak
,
R. S. J.
, &
Turner
,
R.
(
1996
).
Movement-related effects in fMRI time-series
.
Magnetic Resonance in Medicine
,
35
,
346
355
. https://doi.org/10.1002/mrm.1910350312
Goodale
,
S. E.
,
Ahmed
,
N.
,
Zhao
,
C.
,
Zwart
,
J. A. d.
,
Özbay
,
P. S.
,
Picchioni
,
D.
,
Duyn
,
J.
,
Englot
,
D. J.
,
Morgan
,
V. L.
, &
Chang
,
C.
(
2021
).
fMRI-based detection of alertness predicts behavioral response variability
.
eLife
,
10
,
e62376
. https://doi.org/10.7554/elife.62376
Grady
,
C. L.
, &
Garrett
,
D. D.
(
2018
).
Brain signal variability is modulated as a function of internal and external demand in younger and older adults
.
NeuroImage
,
169
,
510
523
. https://doi.org/10.1016/j.neuroimage.2017.12.031
Gu
,
Y.
,
Han
,
F.
,
Sainburg
,
L. E.
, &
Liu
,
X.
(
2020
).
Transient arousal modulations contribute to resting-state functional connectivity changes associated with head motion parameters
.
Cerebral Cortex
,
30
,
5242
5256
. https://doi.org/10.1093/cercor/bhaa096
Gu
,
Y.
,
Sainburg
,
L. E.
,
Kuang
,
S.
,
Han
,
F.
,
Williams
,
J. W.
,
Liu
,
Y.
,
Zhang
,
N.
,
Zhang
,
X.
,
Leopold
,
D. A.
, &
Liu
,
X.
(
2021
).
Brain activity fluctuations propagate aswaves traversing the cortical hierarchy
.
Cerebral Cortex
,
31
,
5242
5256
. https://doi.org/10.1093/cercor/bhab064
Hausman
,
H. K.
,
Hardcastle
,
C.
,
Kraft
,
J. N.
,
Evangelista
,
N. D.
,
Boutzoukas
,
E. M.
,
O’Shea
,
A.
,
Albizu
,
A.
,
Langer
,
K.
,
Etten
,
E. J. V.
,
Bharadwaj
,
P. K.
,
Song
,
H.
,
Smith
,
S. G.
,
Porges
,
E.
,
Hishaw
,
G. A.
,
Wu
,
S.
,
DeKosky
,
S.
,
Alexander
,
G. E.
,
Marsiske
,
M.
,
Cohen
,
R.
, &
Woods
,
A. J.
(
2022
).
The association between head motion during functional magnetic resonance imaging and executive functioning in older adults
.
Neuroimage: Reports
,
2
,
100085
. https://doi.org/10.1016/j.ynirp.2022.100085
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
. https://doi.org/10.1006/nimg.2002.1132
Jenkinson
,
M.
, &
Smith
,
S.
(
2001
).
A global optimisation method for robust affine registration of brain images
.
Medical Image Analysis
,
5
,
143
156
. https://doi.org/10.1016/s1361-8415(01)00036-6
Kumral
,
D.
,
Sansal
,
F.
,
Cesnaite
,
E.
,
Mahjoory
,
K.
,
Al
,
E.
,
Gaebler
,
M.
,
Nikulin
,
V.
, &
Villringer
,
A.
(
2020
).
BOLD and EEG signal variability at rest differently relate to aging in the human brain
.
NeuroImage
,
207
,
116373
. https://doi.org/10.1016/j.neuroimage.2019.116373
Kundu
,
P.
,
Inati
,
S. J.
,
Evans
,
J. W.
,
Luh
,
W. M.
, &
Bandettini
,
P. A.
(
2012
).
Differentiating BOLD and non-BOLD signals in fMRI time series using multi-echo EPI
.
NeuroImage
,
60
,
1759
1770
. https://doi.org/10.1016/j.neuroimage.2011.12.028
Kundu
,
P.
,
Voon
,
V.
,
Balchandani
,
P.
,
Lombardo
,
M. V.
,
Poser
,
B. A.
, &
Bandettini
,
P. A.
(
2017
).
Multi-echo fMRI: A review of applications in fMRI denoising and analysis of BOLD signals
.
NeuroImage
,
154
,
59
80
. https://doi.org/10.1016/j.neuroimage.2017.03.033
Laurita
,
A. C.
,
DuPre
,
E.
,
Ebner
,
N. C.
,
Turner
,
G. R.
, &
Spreng
,
R. N.
(
2020
).
Default network interactivity during mentalizing about known others is modulated by age and social closeness
.
Social Cognitive and Affective Neuroscience
,
15
,
537
549
. https://doi.org/10.1093/scan/nsaa067
Liu
,
T. T.
, &
Falahpour
,
M.
(
2020
).
Vigilance effects in resting-state fMRI
.
Frontiers in Neuroscience
,
14
,
321
. https://doi.org/10.3389/fnins.2020.00321
Liu
,
T. T.
,
Nalci
,
A.
, &
Falahpour
,
M.
(
2017
).
The global signal in fMRI: Nuisance or information
?
NeuroImage
,
150
,
213
229
. https://doi.org/10.1016/j.neuroimage.2017.02.036
Liu
,
X.
,
Zwart
,
J. A. d.
,
Schölvinck
,
M. L.
,
Chang
,
C.
,
Ye
,
F. Q.
,
Leopold
,
D. A.
, &
Duyn
,
J. H.
(
2018
).
Subcortical evidence for a contribution of arousal to fMRI studies of brain activity
.
Nature Communications
,
9
,
395
. https://doi.org/10.1038/s41467-017-02815-3
Nomi
,
J. S.
,
Bolt
,
T. S.
,
Ezie
,
C. C.
,
Uddin
,
L. Q.
, &
Heller
,
A. S.
(
2017
).
Moment-to-moment BOLD signal variability reflects regional changes in neural flexibility across the lifespan
.
The Journal of Neuroscience
,
37
,
5539
5548
. https://doi.org/10.1523/jneurosci.3408-16.2017
Power
,
J. D.
,
Barnes
,
K. A.
,
Snyder
,
A. Z.
,
Schlaggar
,
B. L.
, &
Petersen
,
S. E.
(
2012
).
Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion
.
NeuroImage
,
59
,
2142
2154
. https://doi.org/10.1016/j.neuroimage.2011.10.018
Power
,
J. D.
,
Mitra
,
A.
,
Laumann
,
T. O.
,
Snyder
,
A. Z.
,
Schlaggar
,
B. L.
, &
Petersen
,
S. E.
(
2014
).
Methods to detect, characterize, and remove motion artifact in resting state fMRI
.
NeuroImage
,
84
,
320
341
. https://doi.org/10.1016/j.neuroimage.2013.08.048
Power
,
J. D.
,
Schlaggar
,
B. L.
, &
Petersen
,
S. E.
(
2015
).
Recent progress and outstanding issues in motion correction in resting state fMRI
.
NeuroImage
,
105
,
536
551
. https://doi.org/10.1016/j.neuroimage.2014.10.044
Saccà
,
V.
,
Sarica
,
A.
,
Quattrone
,
A.
,
Rocca
,
F.
,
Quattrone
,
A.
, &
Novellino
,
F.
(
2021
).
Aging effect on head motion: A machine learning study on resting state fMRI data
.
Journal of Neuroscience Methods
,
352
,
109084
. https://doi.org/10.1016/j.jneumeth.2021.109084
Satterthwaite
,
T. D.
,
Wolf
,
D. H.
,
Loughead
,
J.
,
Ruparel
,
K.
,
Elliott
,
M. A.
,
Hakonarson
,
H.
,
Gur
,
R. C.
, &
Gur
,
R. E.
(
2012
).
Impact of in-scanner head motion on multiple measures of functional connectivity: Relevance for studies of neurodevelopment in youth
.
NeuroImage
,
60
,
623
632
. https://doi.org/10.1016/j.neuroimage.2011.12.063
Schaefer
,
A.
,
Kong
,
R.
,
Gordon
,
E. M.
,
Laumann
,
T. O.
,
Zuo
,
X. N.
,
Holmes
,
A. J.
,
Eickhoff
,
S. B.
, &
Yeo
,
B. T. T.
(
2018
).
Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI
.
Cerebral Cortex
,
28
,
3095
3114
. https://doi.org/10.1093/cercor/bhx179
Setton
,
R.
,
Mwilambwe-Tshilobo
,
L.
,
Girn
,
M.
,
Lockrow
,
A. W.
,
Baracchini
,
G.
,
Hughes
,
C.
,
Lowe
,
A. J.
,
Cassidy
,
B. N.
,
Li
,
J.
,
Luh
,
W. M.
,
Bzdok
,
D.
,
Leahy
,
R. M.
,
Ge
,
T.
,
Margulies
,
D. S.
,
Misic
,
B.
,
Bernhardt
,
B. C.
,
Stevens
,
W. D.
,
Brigard
,
F. D.
,
Kundu
,
P.
, …
Spreng
,
R. N.
(
2022
).
Age differences in the functional architecture of the human brain
.
Cerebral Cortex
,
33
,
114
134
. https://doi.org/10.1093/cercor/bhac056
Speck
,
O.
, &
Hennig
,
J.
(
2001
).
Motion correction of parametric fMRI data from multi-slice single-shot multi-echo acquisitions
.
Magnetic Resonance in Medicine
,
46
,
1023
1027
. https://doi.org/10.1002/mrm.1291
Spreng
,
R. N.
,
Setton
,
R.
,
Alter
,
U.
,
Cassidy
,
B. N.
,
Darboh
,
B.
,
DuPre
,
E.
,
Kantarovich
,
K.
,
Lockrow
,
A. W.
,
Mwilambwe-Tshilobo
,
L.
,
Luh
,
W. M.
,
Kundu
,
P.
, &
Turner
,
G. R.
(
2022
).
Neurocognitive aging data release with behavioral, structural and multi-echo functional MRI measures
.
Scientific Data
,
9
,
119
. https://doi.org/10.1038/s41597-022-01231-7
Wong
,
C. W.
,
DeYoung
,
P. N.
, &
Liu
,
T. T.
(
2016
).
Differences in the resting-state fMRI global signal amplitude between the eyes open and eyes closed states are related to changes in EEG vigilance
.
NeuroImage
,
124
,
24
31
. https://doi.org/10.1016/j.neuroimage.2015.08.053
Wong
,
C. W.
,
Olafsson
,
V.
,
Tal
,
O.
, &
Liu
,
T. T.
(
2013
).
The amplitude of the resting-state fMRI global signal is related to EEG vigilance measures
.
NeuroImage
,
83
,
983
990
. https://doi.org/10.1016/j.neuroimage.2013.07.057
Yan
,
C. G.
,
Cheung
,
B.
,
Kelly
,
C.
,
Colcombe
,
S.
,
Craddock
,
R. C.
,
Martino
,
A. D.
,
Li
,
Q.
,
Zuo
,
X. N.
,
Castellanos
,
F. X.
, &
Milham
,
M. P.
(
2013
).
A comprehensive assessment of regional variation in the impact of head micromovements on functional connectomics
.
NeuroImage
,
76
,
183
201
. https://doi.org/10.1016/j.neuroimage.2013.03.004

In this section, we review the theory and implementation of AFNI 3dvolreg. Let Y=[y1,y2,...,yK]N×K be the functional data, where N is the number of voxels, K is the number of volumes, and ykN×1 is the kth column of Y, that is, the kth volume of the data with k=1,2,...,K. Let yrN×1 be the reference image used in the registration. To align yk to yr, the algorithm first estimates the motion parameters ak6×1 (consisting of 3 rotational and 3 translational motion parameters) by iteratively minimizing the weighted least square cost function:

ak=arg mina||W12(ykT(yr,a))||22
(A. 1)

where W is a diagonal weight matrix, and T:N×1N×1 is the function that performs rigid body transformation. Because we disabled the weights, we have W=I, where I is the identity matrix. After obtaining ak, the algorithm then transforms yk with ak, which can be written as T(yk,ak), to register yk to yr. To solve the minimization problem in Eq. A.1, the algorithm linearizes the rigid body transformation with a first-order Taylor approximation, which can be written as

T(yr,a)yr+j=16ajT(yr,a)aj
(A. 2)

where aj is the jth element of a. Furthermore, the partial derivatives in Eq. A.2 are approximated with finite differences:

djT(yr,a)aiT(yr,ϵsj)T(yr,ϵsj)2ϵ,j=1,2,...,6
(A. 3)

where ϵ is a scalar with small magnitude and sj is a 6×1 vector whose jth element is one while the other elements are zeros. Here, we use djN×1 to represent the partial derivative with respect to (w.r.t.) aj and refer to it as the spatial derivative image w.r.t. the jth motion axis. The linearization of the transformation can then be written as

T(yr,a)yr+Da
(A. 4)

where D=[d1,d2,...,d6]N×6. With the above linearization and W=I, the minimization problem in Eq. A.1 can be written as

ak=arg mina||ykyrDa||22
(A. 5)

In AFNI 3dvolreg, the motion parameters are estimated by iteratively minimizing the cost function in Eq. A.5. Let ak,n denote the motion parameters estimated at the nth iteration,

ak,n=(DTD)1DT(yk,nyr)n=1,2,...,Nk
(A. 6)

where Nk is the total number of iterations for the kth volume, and yk,n can be written as

yk,n={yk,n=1T(yk,n1,ak,n1),n=2,3,...,Nk
(A. 7)

The algorithm stops the iterative process when ak,n is smaller than a fixed threshold (for the dataset used in this study, the translational and rotational thresholds are 0.03 mm and 0.02 degrees, respectively) or Nk exceeds the maximum number of iterations, denoted as Nmax. The algorithm uses this iterative process to deal with the approximation error in the linearization of the rigid transformation described in Eq. A.4. In practice, for most volumes, the iterative process stops at the second iteration, suggesting that the parameters estimated in the second iteration are smaller than the thresholds and Eq. A.4 is a valid approximation. After the iterative estimation process stops, the motion parameters of the kth volume are calculated by summing the estimated parameters over all iterations:

ak=n=1Nkak,n
(A. 8)

Stacking the motion parameters over volumes yields

M=[a1T;a2T;...;aKT]
(A. 9)

Here, the jth column of M represents the time series of the motion parameters of the jth motion axis (denoted as mj). Furthermore, in order to explicitly show the relation between M and the functional data Y, we define yk,n=yr for Nk<nNmax so that ak,n=0 for Nk<nNmax. Therefore, the motion parameters can be written as

M=[n=1Nmaxa1,nT;n=1Nmaxa2,n;;n=1NmaxaK,nT]=n=1Nmax[(y1,nyr)T;(y2,nyr)T;;(yK,nyr)T]D(DTD)1=n=1Nmax(YnT1KyrT)D(DTD)1
(A. 10)

where Yn=[y1,n,y2,n,...,yK,n] and 1K is a K by 1 all-one vector. Here, (YnT1KyrT)D(DTD)1 represents the parameters estimated at the nth iteration.

In practice, the algorithm directly estimates and outputs the transformation parameters used to align yk to yr (i.e., ak) by using D. Therefore, in this work, the transformation parameters generated by AFNI 3dvolreg were multiplied by 1 to represent the motion parameters ak.

To explore the key terms underlying r(Δm,GS), we derived the following empirical approximations of the motion estimates:

M=(YT1KyrT)D(DTD)11stiteration+n=2Nmax(YnT1KyrT)D(DTD)1otheriterationsMA1=(YT1KyrT)D(DTD)1(Approx.1)MA2=YTD(DTD)1(Approx.2)MA3=YTD(diag(diag(DTD)))1(Approx.3)
(B. 1)

where 1KK×1 is an all-one vector and diag() denotes the Matlab function diag(). The first approximation of the motion estimates, denoted as MA1, represents the motion parameters estimated at the first iteration and reflects the empirical observation noted above that in practice the iteration process typically stops after the first iteration. Then, because the reference image is almost orthogonal to the derivative images (i.e., yrTD0), we can make the second approximation, denoted as MA2, which represents the motion parameters estimated without subtracting the reference image from the functional data. The third approximation of the motion estimates, denoted as MA3, represents the motion parameters estimated after zeroing out the off-diagonal terms in DTD, reflecting the fact that the covariance between the derivative images in different axes is relatively small compared to the variance of the derivative images. To evaluate the approximations, the temporal correlations between M and the approximated motion estimates were calculated for each motion axis. As shown in Figure S23, M is highly correlated with MA1 (mean r=0.99), MA2 (mean r=0.99), and MA3 (mean r=0.91)). Figure S24 displays the motion estimates before and after the approximations in the Tz and Ty axes from three example runs.

With the third approximation, the motion estimates of one motion axis, denoted as mA3K×1, can be written as

mA3=YTdd2
(B. 2)

where dN×1 is the spatial derivative image w.r.t. that motion axis.

To interpret r(Δm,GS), we derived a mathematical approximation that reveals the key factors. We begin with the expression

r(gs,Δm)=gsTΔmΔmgs=gsT(me2me1)Δmgs
(C. 1)

where gsK×1, Δm=me2me1K×1 represent the GS and the difference between the e2 and e1 motion estimates, respectively, and meiK×1 represents the motion estimates from the ith echo. Then, using the approximation described in Eq. B.2, we can write

r(Δm,GS)=gST(me2me1)||Δm||||gs||gST(Ye2Tde2||de2||2Ye1Tde1||de1||2)||Δm||||gs||=||gs||||Δm||(gsTYe2Tde2||gs||2||de2||2gsTYe1Tde1||gs||2||de1||2)=||gs||||Δm||(βgs,e2Tde2||de2||2βgs,e1Tde1||de1||2)
(C. 2)

where βgs,ei=YeiTgsgs2 is the GS beta coefficient map calculated from the linear fit of the GS to the unregistered functional data of the ith echo.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC BY 4.0) license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode.

Supplementary data