Resting-state global brain activity induces bias in fMRI motion estimates

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 S₀ 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.

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 mm³; 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 mm³; 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.

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.

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).

To characterize BOLD-weighted GS bias in the motion estimates, we considered a simple signal model for the e1 and e2 motion estimates. Let $me1∈ℝK×1$ and $me2∈ℝK×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,

where $mh∈ℝK×1$ represents head motion, $αei$ is the regression weight corresponding to head motion for the $i$th echo, $gs∈ℝK×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 $i$th echo, and $ϵei∈ℝK×1$ is the estimation error for the $i$th echo. With the above model, subtracting $me1$ from $me2$ yields

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

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 $Δm$and 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.

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,

where $βgs,ei∈ℝN×1$ is the GS beta coefficient map for the $i$th echo, $dei∈ℝN×1$ is the spatial derivative image with respect to (w.r.t.) one motion axis of the $i$th echo, $N$ is the number of voxels, and $∥⋅∥$ denotes the $L$2-norm. For each run, the GS beta coefficient map $βgs,ei=Yeigs∥gs∥2$ was calculated from the linear fit of the GS to the unregistered functional data $Yei∈ℝN×K$ of the $i$th 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 $i$th echo as $yr,ei∈ℝN×1$⁠, the spatial derivative image of the $i$th echo w.r.t the $j$th motion axis was calculated as

where $Tj:ℝN×1→ℝN×1$ is the function that transforms the reference volume along the $j$th 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.

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.

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×10−3$ and $1×10−6$⁠. 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.

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/20$th 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$⁠).

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.

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.

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.

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 $βgsTd∥d∥2$ 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$⁠, $d∥d∥2$⁠, $βgs​⊙d∥d∥2$⁠, and the $βgs​⊙d∥d∥2$ 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 $d∥d∥2$ maps indicate voxels with relatively higher spatial derivatives; larger magnitudes in the element-wise multiplication $βgs​⊙d∥d∥2$ 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,e2⊙de2∥de2∥2$ are greater than the corresponding values in the element-wise map for the first echo $βgs,e1⊙​de1∥de1∥2$⁠, where the differences are driven for the most part by the greater $βgs$ values observed for the second echo.

In the Tz axis, $d∥d∥2$ 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 $βgs​⊙d∥d∥2$ are positive with greater amplitudes for e2, resulting in high positive values in the $βgs​⊙d∥d∥2$ difference map. Consequently, summing up the $βgs​⊙d∥d∥2$ difference map leads to a high $r(Δm,GS)$ value. In contrast, for the Tx axis, $d∥d∥2$ 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 $βgs​⊙d∥d∥2$ 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, $d∥d∥2$ 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 $βgs​⊙d∥d∥2$ 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 $βgs​⊙d∥d∥2$ 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 $βgs​⊙d∥d∥2$ 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$⁠.

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.

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).

Furthermore, a one-way ANOVA on the mean $Δz$ values calculated over ROI pairs showed that there is a significant (⁠$p<1×10−6$⁠, $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×10−4$⁠) more negative mean $Δz$ values as compared to the other three groups, whereas the high motion and low aGS runs exhibit significantly (⁠$p<1×10−3$⁠) 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×10−6$⁠) 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.

As shown in Figure 6 (a) and (b), young subjects show significantly (⁠$p<1×10−6$⁠) 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×10−6$⁠) 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.

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.

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.

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 $yk​∈ℝN×1$ is the $k$th column of $Y$⁠, that is, the $k$th volume of the data with $k=1,2,...,K$⁠. Let $yr∈ℝN×1$ be the reference image used in the registration. To align $yk$ to $yr$⁠, the algorithm first estimates the motion parameters $ak​∈ℝ6×1$ (consisting of 3 rotational and 3 translational motion parameters) by iteratively minimizing the weighted least square cost function:

where $W$ is a diagonal weight matrix, and $T:ℝN×1→ℝN×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

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

where $ϵ$ is a scalar with small magnitude and $sj$ is a $6×1$ vector whose $j$th element is one while the other elements are zeros. Here, we use $dj∈ℝN×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 $j$th motion axis. The linearization of the transformation can then be written as

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

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 $n$th iteration,

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

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 $k$th volume are calculated by summing the estimated parameters over all iterations:

Stacking the motion parameters over volumes yields

Here, the $j$th column of $M$ represents the time series of the motion parameters of the $j$th 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​<n≤Nmax$ so that $ak,n=0$ for $Nk<n≤Nmax$⁠. Therefore, the motion parameters can be written as

where $Yn=[y1,n,y2,n,...,yK,n]$ and $1K$ is a $K$ by $1$ all-one vector. Here, $(YnT−1KyrT)D(DTD)−1$ represents the parameters estimated at the $n$th 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:

where $1K∈ℝK×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., $yrTD≈0$⁠), 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 $mA3∈ℝK×1,$ can be written as

where $d∈ℝN×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

where $gs∈ℝK×1$⁠, $Δm=me2−me1∈ℝK×1$ represent the GS and the difference between the e2 and e1 motion estimates, respectively, and $mei∈ℝK×1$ represents the motion estimates from the $i$th echo. Then, using the approximation described in Eq. B.2, we can write

where $βgs​,ei=YeiTgs∥gs∥2$ is the GS beta coefficient map calculated from the linear fit of the GS to the unregistered functional data of the $i$th echo.