Resting state global brain activity induces bias in fMRI motion estimates

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


Introduction
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 (Friston et al., 1996;Power et al., 2012;Satterthwaite et al., 2012;Dijk et al., 2012;Yan et al., 2013;Power et al., 2014Power et al., , 2015)).Prior taskbased fMRI studies demonstrated that task-evoked brain activations can induce bias in motion estimates that is temporally correlated with brain activations (Freire and 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 TEdependent 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 (Buur et al., 2009;Bright and Murphy, 2013;Kundu et al., 2012Kundu et al., , 2017)).For example, Burr et al. effectively used the first echo data acquired at a short TE with minimal BOLD-weighting to model and correct for the motion artifacts in the BOLD-weighted second echo data (Buur et al., 2009).
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 TEindependent 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 and 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 and 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 (Liu and Falahpour, 2020)), there is growing evidence suggesting that the GS is linked to global neural activity (Wong et al., 2013(Wong et al., , 2016;;Falahpour et al., 2018;Liu et al., 2018;Gu et al., 2021).
In this work, we characterized the bias in rsfMRI motion parameters as estimated by AFNI 3dvolreg.
Furthermore, we investigated the consequences of using biased estimates as regressors in resting-state functional connectivity (rsFC) analyses.

Data preprocessing
AFNI was used for data preprocessing (Cox, 1996).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.For each run and each echo, the data were normalized so that the mean signal over all voxels and all volumes was equal to 100.Since the normalization is equivalent to multiplication by a global scaling factor, the spatial and temporal information of the data in each run is preserved.

Calculation of the motion estimates and global signal
In this study, we examined the motion estimates from AFNI 3dvolreg (Cox and 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.For each run, the motion was estimated using both the e1 and e2 data.Then, 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 unregistered e2 data.Before calculating the GS, each voxel's signal was converted to a percent change signal.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.

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 m e1 ∈ R K×1 and m e2 ∈ R 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 m h ∈ R K×1 represents head motion, α ei is the regression weight corresponding to head motion for the ith echo, g s ∈ R K×1 is the GS, TE 1 and TE 2 are the first and second echo times, respectively, β ei is the regression weight of the BOLD-weighted GS bias for the ith echo and ϵ ei ∈ R K×1 is the estimation error for the ith echo.With the above model, subtracting (2) Since prior findings suggest that head motion can be accurately estimated from both e1 and e2 data (Buur et al., 2009;Speck and Hennig, 2001), we assume that 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 correlations 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, i.e. 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.

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 Appendix C, where β g s ,ei ∈ R N ×1 is the GS beta coefficient map for the ith echo, d ei ∈ R N ×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 β g s ,ei = Y ei g s ∥g s ∥ 2 was calculated from the linear fit of the GS to the unregistered functional data Y ei ∈ R N ×K of the ith echo.
We calculated the spatial derivative images w.r.t.motion axes following the calculation 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 y r,ei ∈ R N ×1 , the spatial derivative image of the ith echo w.r.t the jth motion axis was calculated as where T j : R N ×1 → R N ×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, 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 • .The exact values were used in AFNI 3dvolreg when creating the derivative images and were calculated based on the spatial resolution of the functional data.

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 task positive network (TPN).The regions of interest (ROIs) within these networks were defined in (Dijk et al., 2010), with four ROIs in the DMN (the posterior cingulate cortex (PCC), lateral parietal cortex (LatPar), medial prefrontal cortex (mPFC) and Hippocampal formation (HF)) and three ROIS in the TPN (frontal eye field (FEF), intraparietal cortex (IPS) and middle temporal area (MT+)).Seed ROIs were created using a sphere with a radius of 12 mm centered about each seed coordinate.The left and right ROIs were combined to form bilateral ROIs.Prior to averaging signals within each ROI, the e2 data after volume registration were transferred to MNI space and spatially smoothed with a 4mm FWHM Gaussian kernel.
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 showed in the supplementary materials (Fig. S4 and S5) 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 Pearson's correlation coefficients between the ROI average signals.Correlation values were converted to z-scores using the Fisher-z transformation.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.
Furthermore, we examined the effect of the bias on FC within four groups of runs based on levels of head motion and GS.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.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.

Effect of the bias on young vs. old group-level FC analysis
In this study, we used a public dataset consisting of young and old subjects (181 younger and 120 older adults as described in Setton et al. (2022)).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 × 10 7 random permutations to allow us to examine the commonly used p-values 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.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.34 as compared to the low motion run, and reflects motion artifacts in the GS.After motion regression (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.

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.For the other motion axes (Rx, Ry, Rz and Tx), the percent of significant r(∆m, GS) values fluctuates around 5%, ranging from 1.2% to 6.6%, indicating minimal BOLD-weighted GS bias for these axes.

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 multiplying the difference of β T g s d ∥d∥ 2 from e2 and e1 by a scaling factor 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 is a result of β g s and d sharing a relatively similar top-bottom asymmetric spatial pattern.
In the Ty axis, as shown in Fig. S3, d ∥d∥ 2 exhibits negative and positive values at the anterior and posterior edges of the brain, respectively, whereas β g s 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 β g s ⊙d ∥d∥ 2 difference 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 relfect the asymmetries in β g s maps, which show higher values at the superior and posterior parts of the brain as compared to the inferior and anterior parts of the brain.

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.
Note that motion censoring was applied before motion regression.There is no significant (p > 0.01) difference between the group of low motion and low aGS runs and the group of high motion and high aGS runs.The distributions of the mean ∆z values over ROI pairs for the four groups can be found in Fig. S5.
Fig. S6 and S7 demonstrate that the observed effect of the bias is dominated by the motion estimates in the Ty and Tz axes.In Fig. S6, when including only the Ty and Tz motion regressors, the results are similar to those shown in Fig. 4. In contrast, when excluding the Ty and Tz motion regressors (Fig. S7), 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, and there is a stronger effect for runs with higher aGS and lower head motion.

Effect of the bias on young vs. old group-level analysis
As shown in Fig. 5 (a) and (b), young subjects show significantly (p < 1 × 10 −6 ) higher aGS and lower mean FD values as compared to the old subjects.Consequently, as shown in Fig. 5 (c), the young subjects show significantly (p < 1 × 10 −6 ) more negative mean ∆z values over ROI pairs as compared to the old subjects.The mean ∆z values over runs and ROI pairs are -0.69 and -0.17 for the young and old subjects, respectively.Furthermore, Fig. 5 (d) and (e) visualize the mean ∆r and ∆z values over the young and old subjects for all ROI pairs.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.
Moreover, we examined whether the bias affects the group-level analysis between the young and the old subjects.Fig. 6 (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 < 1 × 10 −3 ) higher FC than the old subjects for all ROI pairs, except for the pair between IPS and LatPar.However, after e2 motion regression, as shown in Fig. 6 (b), the significant differences for five pairs of ROIs become insignificant (p > 0.01), including the pairs between (1) PCC and FEF, (2) PCC and IPS, (3) LatPar and FEF, (4) LatPar and MT+ and (5) FEF and IPS.Additionally, Fig. 6 (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.Furthermore, as shown in Fig. S8 and S9, we verified that the observed effect is dominated by regressing out the Ty and Tz motion regressors.Together, these results demonstrate that motion regression with biased motion estimates can lead to underestimation of FC differences between the young and old subjects.

Discussion
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.In contrast, in (Freire and Mangin, 2001;Freire et al., 2002), the existence of the bias in task-based fMRI is primarily illustrated using simulation with experimental data limited to one scan.Moreoever, 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 While the interpretation of the GS is still controversial (reviewed in (Liu and Falahpour, 2020)), there is growing evidence suggesting that the GS is linked to vigilance (also known as arousal level) (Wong et al., 2013(Wong et al., , 2016;;Falahpour et al., 2016Falahpour et al., , 2018;;Liu et al., 2018;Gu et al., 2021).For example, Falahpour et al. (2016) found that the GS is negatively correlated with the EEG vigilance.Furthermore, studies have shown that the vigilance templates, calculated as voxel-wise correlations between the EEG vigilance and the fMRI signal, can be used to estimate the EEG vigilance fluctuations during the fMRI scans (Chang et al., 2016;Falahpour et al., 2018;Goodale et al., 2021).To investigate whether the asymmetric spatial patterns in β g s are associated with the vigilance, Fig. S10 compares β g s from the example run shown in Fig. 3 to the eyes-open vigilance template estimated in (Falahpour et al., 2018) (both spatial maps were transferred to the MNI152 standard space).As shown in Fig. S10 (b), the estimated vigilance template exhibits similar superior-inferior and posterior-anterior asymmetric spatial patterns to the example β g s map, resulting in a strong negative correlation (r = −0.38,p < 1 × 10 −6 ) between β g s and the vigilance template.Hence, the asymmetric spatial patterns in β g s may be attributed to resting-state global brain activity linked to vigilance changes.
We also found that the posterior-anterior asymmetry pattern in β g s may be associated with the spatial intensity inhomogeneity in the fMRI data, a slow and smooth variation related to a combination of radio-frequency coil homogeneity, the pulse sequences used and the imaged object itself (Belaroussi et al., 2006).To illustrate the effect of the intensity inhomogeneity, Fig. S11 shows the reference images and estimated bias field maps from two low motion example runs, one with a high Ty r(∆m, GS) value of 0.87 (the example run shown in Fig. S3) and one with a low Ty r(∆m, GS) value of 0.07 (an example run has not shown before).For the example run with a strong r(∆m, GS) in the Ty axis shown in Fig. S11 (a), its reference image and estimated bias field map exhibit higher values at the posterior part of the brain and lower values at the anterior part of the brain.Whereas, for the run with r(∆m, GS) close to zero (shown in Fig. S11 (b)), its reference image and estimated bias field map show comparable values at the anterior and posterior parts of the brain.
While this work primarily focused on motion estimates from AFNI 3dvolreg, we also performed a supplementary examination of the r(∆m, GS) values calculated with motion estimates from different algorithms, including FSL (Jenkinson and Smith, 2001;Jenkinson et al., 2002), ANTS (Avants et al., 2011(Avants et al., , 2014) ) and SPM (Friston et al., 1995) and cost functions, including the least-squares (LS), normalized correlation (NC) and mutual information (MI) cost functions.(See sections S1.1 and S2.2 in the supplementary materials for details).In conclusion, both the algorithms and cost functions affect the level of the bias in the Ty and Tz axes, with the algorithms having a greater impact than the cost functions.
Specifically, AFNI (LS) and SPM (LS) show a relatively greater level of the bias, FSL with different cost functions (LS, NC and MI) shows a moderate level of the bias, and ANTS (MI) shows the lowest level of the bias.Future work is needed to thoroughly examine how different algorithms and cost functions affect the level of the bias in the motion estimates.
In this study, we demonstrated that regression with biased motion estimates can reduce rsFC differences between the young and old subjects with the DMN and TPN ROIs.Importantly, we showed in Fig. 6 that the significant rsFC differences for five pairs of ROIs became insignificant (p > 0.01) due to regression with the biased motion estimates.Among these five pairs of ROIs, four of them are connections between DMN and TPN.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 aged brain (reviewed in (Ferreira and 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 can use motion estimates from the first echo data (if available) or consider excluding the Ty and Tz motion estimates from the motion regression.Moreover, our findings indicate that the biased motion estimates may lead to discrepancies in rsFC results via regression.Therefore, when comparing results from existing aging studies with different processing methodologies, the potential effect of the bias needs to be considered.Concerning this matter, future work that thoroughly examines the effect of the bias on the young vs. old rsFC analysis would be helpful.
Our findings revealed that the effect of the bias on the young vs. 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.It has been reported that old subjects tend to have greater head motion during rsfMRI scans (Saccà et al., 2021;Hausman et al., 2022), which reflects declined executive functioning with aging (Hausman et al., 2022).On the other hand, a direct comparison of aGS between young and old groups appears to be outside of prior studies.Nevertheless, 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 and Busatto, 2013;Grady and Garrett, 2018;Nomi et al., 2017;Kumral et al., 2020), 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 ROI-ROI rsFC obtained after motion regression, with an increasing level of effect for runs showing higher aGS and lower head motion.Moreover, our results show that the bias can reduce group-level FC differences between young and old subjects.To minimize the effect of the bias, investigators can use motion estimates from the first echo data (if available) or consider excluding the Ty and Tz motion estimates from the motion regression.Future work is needed to investigate the dependence of the bias on the registration algorithm and scan parameters.In addition, future work focused on developing registration algorithms with minimal bias would be helpful.A positive value (red color) indicates that the young subjects show higher connectivity as compared to the old subjects.

Acknowledgements
(c) Scatter plot comparing 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.(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.e1 MotReg: e1 motion regression; e2 MotReg: e2 motion regression.
can be written as where d ∈ R N ×1 is the spatial derivative image w.r.t. that motion axis.

Fig. 1
Fig.1 (a,b) show the motion estimates in the Tz axis, including m e1 (blue), m e2 (green) and ∆m (red) from two example runs with (a) low and (b) high levels of motion.Note that for the high motion run, m e1 and m e2 are plotted at 1/20th scale to facilitate comparison with m e1 and m e2 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 m e1 and m e2 from the high motion run (std(m e1 ) = 0.583, std(m e2 ) = 0.602) are an order of magnitude larger than the standard deviations of m e1 and m e2 from the low motion run (std(m e1 ) = 0.027, std(m e2 ) = 0.034).Fig.1 (c) and (d) show the GS and ∆m from the low motion and high motion runs, respectively.
Fig. 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, 61.1% 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.55 to 0.58 after e1 motion regression.Remarkably, 96 out of 602 total runs show r(∆m, GS) values larger than 0.8 after e1 motion regression.In the Ty axis, 30.6% and 38.2% 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.29 to 0.34 after e1 motion regression.Fig.S1 and S2show examples of the GS and the motion estimates, including m e1 , m e2 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.
∥g s ∥ ∥∆m∥ .The approximation provides us a unique angle to interpret the GS bias by looking at the relation of the β g s and d spatial maps.Fig. 3 visualizes β g s , d ∥d∥ 2 , β g s ⊙d ∥d∥ 2 , and the β g s ⊙d ∥d∥ 2 difference maps in the Tz and Tx axes from an example run, where ⊙ represents element-wise multiplication.In the Tz axis, d ∥d∥ 2 exhibits high positive and negative values at the superior and inferior edges of the brain, respectively, whereas β g s 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 β ∆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 β g s shows high positive values across the cortex, including both the left and right edges.The resulting high positive and negative values observed in the β g s ⊙d ∥d∥ 2 Fig.4shows ∆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.09 and ∆z ranging from -0.63 to -1.41), 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.05 to -0.33).Furthermore, a one-way ANOVA on the mean ∆z values calculated over ROI pairs showed that there is a significant (p < 1 × 10 −6 , F 3,598 = 40.53)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 −6 ) less negative mean ∆z values as compared to the other three groups.

Figure 1 :
Figure 1: Motion estimates in the Tz axis, including m e1 (blue), m e2 (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 (MotReg) from the low and high motion runs, respectively.

Figure 3 :
Figure 3: Visualization of the spatial maps underlying r(∆m, GS) in the Tz and Tx 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 β g s , d ∥d∥ 2 , and β g s ⊙d ∥d∥ 2 maps from e1 and e2, where ⊙ represents element-wise multiplication.Row 4 shows

Figure 4 :
Figure 4: 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 averaged differences in r-values and a lower left triangle showing the the averaged differences in z-scores.

Figure 5 :
Figure 5: (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 × 10 −6 ) 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 averaged differences in r-values and a lower left triangle showing the the averaged differences in z-scores.