18 healthy young subjects (age = 26.7 ± 6.5 years, 9F/9M) and 18 healthy older subjects (age = 74.2 ± 7.0 years, 11F/7M) were imaged using a Siemens TIM Trio 3 T scanner (Siemens Healthineers, Erlangen, Germany). All participants were recruited from the Greater Toronto Area, and provided written informed consent as per the policy of our institutional research ethics board. rs-fMRI scans were collected with a 32-channel head coil using simultaneous multi-slice (SMS) echo-planar imaging (EPI) BOLD (TR/TE = 380/30 ms, FA = 40°, 21 5-mm slices, 64×64 matrix, 4x4x5 mm voxels, multiband factor = 3, 1,950 volumes, left–right phase encoding direction, in-plane acceleration of 2). During each scan, cardiac pulsation was recorded using the scanner pulse oximeter, whereas the respiratory signal was recorded using a Biopac^™ system (Biopac Systems Inc. California, United States). A T1-weighted 3D anatomical data set (1 mm isotropic resolution) was also acquired for each participant (MPRAGE, TR = 2,400 ms, TE = 2.43 ms, field-of-view = 256 mm, voxel size = 1 mm isotropic, TI = 1,000 ms, BW = 180 Hz/vox).

The anatomical segmentation was performed using fMRIPrep, as shown in Figure 1. Specifically, the T1 anatomical images were skull-stripped using a Nipype implementation of Advanced Normalization Tools (ANTs) (Avants et al., 2011), following which tissue segmentation was performed using FMRIB Software Library (FSL) and spatial normalization to standard MNI152 space was performed using ANTs. This step resulted in grey matter (GM), WM and CSF segmentations for each data set.

fMRI data preprocessing was performed using fMRIPrep (Version 22.0) (Esteban et al., 2019), and includes motion correction, spatial smoothing (5 mm FWHM), high-pass filtering (>0.01 Hz) and brain extraction. This is a common processing pipeline applied prior to all physiological denoising methods (no correction), though we also obtained framewise displacement (FD) time series for each data set based on Power et al. (2012). Temporal SNR (tSNR) on the preprocessed data (pre-denoising) was also computed for each subject, as the ratio between the mean and standard deviation of each time series. The preprocessed data was channeled separately through five data-driven denoising methods for comparison, as implemented through fMRIPrep. These methods include:

Regressors of each method are estimated using fMRIprep (Esteban et al., 2019), publicly available at fmriprep.org. We did not regress out estimated motion parameters in addition to the noise regressors.

Cardiac pulsation and respiration are the two chief sources of physiological noise in the fMRI signal, especially problematic in rs-fMRI. Cardiac pulsations generate localized brain motion as well as creating inflow effects in and around blood vessels and CSF. On the other hand, thoracic movement during breathing results in magnetic field alteration and phase shift in the image. Moreover, respiration induces fluctuation in the level of arterial CO₂, which has a vasodilating effect and can induce changes in the fMRI BOLD signal. The frequency of physiological signals in an adult human typically falls within the range of 0.2 to 0.3 Hz for respiratory and 0.8 to 1.2 Hz for cardiac signals. In a conventional fMRI scan with a TR of 2 s, these signals alias into the low-frequency range (0.01 to 0.1 Hz), which is known to contain information about fcMRI (Cordes et al., 2001). Physiological denoising in a conventional fMRI dataset may result in a decrease in signal power in the low-frequency range. However, it is impossible to determine whether the removed power is related to the respiratory and cardiac contributions or to neuronal connectivity. The low-frequency range of fMRI signals has signatures of several non-neuronal sources, including systemic low frequency oscillations (Tong et al., 2015, 2019), respiratory variations (Birn et al., 2006), and cardiac rate variations (Chang et al., 2009). Also, the fact that most of the denoising methods exhibit the highest impact on signal spectrum in the low-frequency band may be due to the nature of ROI-based noise regressor generation, whereby high-frequency regressor content may be averaged out through different phase shifts across the ROIs, leading to the retention of low-frequency regressor power.

For context, the methods evaluated in this study are among the most commonly used in the resting-state fMRI literature. Specific to the study of the effect of aging on resting-state connectivity, all of the denoising methods have found broad application, including GSR (Betzel et al., 2014; Chan et al., 2014; Zhang et al., 2014; Geerligs et al., 2015; Siman-Tov et al., 2016; Stumme et al., 2020), WM-CSF (Koch et al., 2010; Betzel et al., 2014; Chan et al., 2014; Song et al., 2014; Zhang et al., 2014; Grady et al., 2016; Siman-Tov et al., 2016; Farras-Permanyer et al., 2019; Varangis et al., 2019; Mancho-Fora et al., 2020; Xie et al., 2020; Zhong and Chen, 2022), CompCor (Onoda and Yamaguchi, 2013; Hausman et al., 2020; Hamada et al., 2021; Patil et al., 2021; Podgórski et al., 2021), and ICA-AROMA (Stumme et al., 2020). Moreover, some recent studies of aging used no physiological denoising (Damoiseaux et al., 2008; Huang et al., 2015; Ao et al., 2022).

Owing to the high sampling rate of our fMRI data, we were able to dissect the differential effect of different denoising methods by examining the resultant fcMRI signal power spectrum. All methods resulted in spectral-power reduction in all three frequency bands of the GM signal, removing as high as 60–70% of the signal in the respiratory and cardiac bands. The amount of the removed power spectrum is vastly different among the methods, with GS and AROMA removing the most power followed by aCompCor. ICA-AROMA also removed more signal power in the physiological frequency bands compared to all other denoising methods. These findings are consistent with the finding of signal quality improvement by GSR and ICA-FIX by a previous study (Geerligs et al., 2017). Moreover, consistent with this finding, aCompCor was previously found to be less effective than ICA-FIX at reducing motion artifacts (Scheel et al., 2022).

In our findings, aCompCor and tCompCor removed significantly less power than GS and AROMA in all three frequency bands in both the GM and WM of both age groups, with only WM-CSF regression showing less power change (Tables 2A,B). This finding for aCompCor is in contrast to the findings of Agrawal et al. (2020) who found aCompCor to reduce low-frequency power while adding high-frequency noise power. Conversely, Muschelli et al., found aCompCor to remove more head-motion related signal contribution than alternatives such as WM-CSF regression (Muschelli et al., 2014). However, these latter findings were restricted to young adults, which suggests these discrepancies in findings may be related to the significance of age x denoising-method interactions reported in Table 1. Indeed, as reported in Supplementary Table S2, both mean and max FD were higher in the older adults. Nonetheless, no physiological frequency or head motion differences across age groups reached statistical significance.

Also, as mentioned earlier, Dipasquale et al. (2017) found ICA-AROMA to retained more low-frequency spectral power while reducing more high-frequency noise power than WM-CSF, which is in contrast to our findings. However, we wish to highlight that in order to isolate the effect of physiological-denoising, we did not regress out motion alongside the physiological noise regressors in this study, in contrast to previous studies (Dipasquale et al., 2017; Geerligs et al., 2017; Kassinopoulos and Mitsis, 2022). Thus, we cannot directly compare our results to those of Dipasquale et al., which included WM-CSF alongside motion regression. The interaction between motion and noise is complex (Jones et al., 2008), and the order of each type of correction as well as the number of motion regressors can sway the results (Ciric et al., 2017; Parkes et al., 2018), and thus these questions are beyond the scope of this work. Moreover, additional differences in findings could be attributed to the major differences in the TR of the BOLD acquisition. In the Dipasquale study, TR was 2.75 s, while in our study, it was 0.38 s. This difference in BOLD signal sampling rate could affect the performance of regression-based approaches such as WM-CSF, which do not discriminate between different temporal and spatial constituents in the WM and CSF ROIs. The evaluation of the effects of sampling rate and head motion on different denoising methods is a promising future direction.

Compared to the GM, the results in the WM are more variable across subjects, and WM denoising is much less examined in the literature. Notably, among all the methods, only AROMA appeared to reduce cardiac power (Figure 1D). All methods reduced power in the respiratory and low-frequency frequency bands (Figures 1E,F), with aCompCor and AROMA leading the other methods in this regard. Nevertheless, the power-reduction variability in the WM is much larger than that in GM, especially for GSR, WM-CSF and aCompCor. GSR is also very effective at suppressing the contribution of physiological frequencies in the GM, but not very effective in the WM. Unsurprisingly, both GSR and WM-CSF regression introduced increases to cardiac power, and are thus inappropriate for WM fMRI denoising. This observation could be attributed to the manner in which the noise regressors are derived for these methods, in which the WM signal is indiscriminately designated as noise. Interestingly, WM-CSF did not significantly reduce WM BOLD signal power in the cardiac and respiratory power, although this method can be viewed as involving regressing the WM signal out of itself. Although aCompCor also derives regressors from WM, its performance is superior to that of WM-CSF, as it only designated the leading principal components of the WM signal (and CSF) signal as noise regressors (Figure 1E).

It has been established that head-motion effects interact with physiological denoising (Jones et al., 2008) and age-related fcMRI differences (Mowinckel et al., 2012; Geerligs et al., 2017). Lesser known, but critical, is the additional influence of fMRI sampling rate, especially with the advent of fast fMRI (Chen et al., 2019; Power et al., 2019). We found a sub-significant age-group difference in head motion (FD), with older adults associated with higher FD (Supplementary Table S2).

Given our high sampling rate, we will elaborate on the unique effect of fast sampling on the treatment of head motion and the measurement of denoising outcome. In the acquisition used in this study, we used an SMS acceleration factor of 3 and 21 slices at TR = 0.38 ms, translating into a maximum samplable frequency of 1.3 Hz. This frequency is sufficient to capture the fundamental cardiac frequency, preventing its aliasing into lower-frequency bands. The SMS acquisition scheme also results in the shuffling of slice acquisition ordering across 7 slice groups. If residual head motion chiefly manifests as through-plane movement between neighboring voxels, then all slices in each slice group would shift simultaneously, expanding the influence of motion across the brain volume. Additionally, Power et al. (2019) reported on lower-frequency motion being subdivided by fast sampling (TR of 720 ms) and thus appearing as lower-amplitude higher-frequency FD that overlaps in frequency with respiratory noise. This effect of fast sampling is further complicated by slice groups being ordered in an interleaved fashion. That is, the slice group ordering was [1, 8, 15], [3, 10, 17], [5, 12, 19], … [2, 9, 16], [4, 11, 18], and so on, such that physiological and motion effects that should be similar across neighboring slices are sub-sampled at up to TR/2 apart. Each slice group would thus sub-sample the head-motion time course with a different phase shift, resulting in potentially elevated within-ROI signal heterogeneity for methods such as GSR and WM-CSF regression.

Relating this to different sampling rates, note that the use of a TR of 720 ms in Power et al.’s work sets the maximum samplable BOLD signal frequency at 0.7 Hz. At this sampling rate, the cardiac frequency (typically ~1 Hz) would alias into the 0.3 Hz frequency band, which coincides with the typical respiratory frequency. Thus, manifestation of fast-sampled motion, and consequently the effect of different denoising methods on the fast-sampling related motion-subsampling effects, which has previously been reported (Burgess et al., 2016; Ciric et al., 2017; Dipasquale et al., 2017), also depends on the actual sampling frequency and SMS factor.

The low-frequency band’s interpretation remains complex, however. On the one hand, it is the established spectral location for the majority of the neurovascular contribution to brain functional connectivity (Cordes et al., 2001), and retention of fcMRI strength has been used as a metric for evaluating retention of meaningful BOLD signal (Dipasquale et al., 2017). Nonetheless, variability in the respiratory and cardiac pattern can also induce BOLD signal variations in the low-frequency frequency band (Wise et al., 2004; Birn et al., 2006; Chang et al., 2013). Recent work, however, demonstrated that respiratory-volume variability may have neuronal associations (Yuan et al., 2013; Shams et al., 2021), rendering its interpretation uncertain.

Our results showed that GSR and AROMA removed the most low-frequency BOLD signal power. Although low-frequency power did not significantly modulate GSR-based functional connectivity (Pearson’s correlation) or FCC as a measure of network segregation (Kassinopoulos and Mitsis, 2022), GSR and AROMA resulted in two of the lowest FCC values in comparison to the other methods (see Table 3B). Our downsampled data (Supplementary Figure S3 and Supplementary Table S4) followed very similar trends as the original data, suggesting that our findings are generalizable to different sampling rates, and that the effects of band-limited noise, once shifted to different frequency bands due to aliasing, may retain the same impact on functional connectivity metrics.

Moreover, AROMA-based Pearson’s correlation was significantly mediated by both cardiac and low-frequency power (Supplementary Table S5A). In contrast, WM-CSF and aCompCor were associated with the highest FCC values as demonstrated in (Figure 5B; Table 3B). Still, based on the power spectral analysis alone, it is impossible to understand whether the removed low-frequency spectrum can alter brain connectivity as measured by rs-fMRI. Furthermore, it is possible that a higher fcMRI measurement does not result from higher low-frequency power or lower noise power (Agrawal et al., 2020). Lastly, the reduced observable age-related differences in fcMRI measures associated with GSR and AROMA was not significantly mediated by spectral-power changes in any of the three frequency bands (Supplementary Table S5C), despite strong correlations between low-frequency power change and modularity in both GSR and AROMA (0.46 and − 0.74, respectively). Interestingly, in aCompCor alone, was spectral power change significantly associated with modularity index (positively for respiratory power change, and negatively for both cardiac and low-frequency power). These correlations are reflected also in the findings for tCompCor and AROMA, suggesting that these methods may share the property whereby higher age-related differences in connectivity may be linked to lower reductions in low-frequency and cardiac power and greater reduction in respiratory power. This is an interesting scenario for future investigations.

We also evaluated the denoising methods based on their effects on the observable age-related rs-fcMRI differences. Although aCompCor and tCompCor did not result in the greatest retention of low-frequency BOLD signal power (Figure 2; Table 2), they generated maps with the highest modularity indices (a measure of sensitivity to age-related differences), as demonstrated in Figure 5A and Table 3A. Moreover, we observe a significant age x denoising-method interaction when examining the fractional spectral power changes in the low-frequency and respiratory bands (Table 1).

In aging studies of fcMRI, rs-fcMRI reduction in the DMN is the most commonly reported finding [for a review please refer to Jockwitz and Caspers (2021)], although the extent and location of the connectivity reduction differ among the studies, consistent with our finding (Figure 2). For instance, Stumme et al. (2020) used AROMA denoising and reported reduced connectivity in the visual network and no changes in the DMN and somatomotor networks, whereas other studies that used tCompCor or WM-CSF denoising have reported no change (Zhang et al., 2014) or reduced (Onoda and Yamaguchi, 2013; Chan et al., 2014; Zhang et al., 2014) connectivity in the visual network and reduced connectivity in the DMN (Onoda and Yamaguchi, 2013; Betzel et al., 2014; Chan et al., 2014; Zhang et al., 2014; Grady et al., 2016), as well as no connectivity change in the somatomotor network (Onoda and Yamaguchi, 2013; Betzel et al., 2014; Chan et al., 2014). Here, we demonstrate that one possible source of the inconsistencies across studies is the choice of method for physiological denoising. Despite the fact that the spectral power removed by different denoising methods were at times not significantly different, these differences resulted in observable differences in the resultant observable age-related fcMRI differences. Specifically, as shown in Figure 5A and Table 3A, the GSR method also visibly reduced the observable age difference in fcMRI, while CompCor and WM-CSF minimally impacted the observable age difference in connectivity. Nonetheless, the ability to detect age differences does not equate that of detecting neuronally-driven age differences.

In the systematic comparison of WM-CSF and aCompCor reported by Geerligs et al. (2017), WM-CSF regression was associated with higher connectivity strength and reliability as well as reliability of age effects than aCompCor. Moreover, Geerligs et al. showed that the WM-CSF regression resulted in a negative age-fcMRI association, while aCompCor resulted in a positive age-fcMRI association. We did not test for reliability of our fcMRI metrics, preventing a direct comparison of results in that regard. However, we did not observe a significant difference in modularity index between WM-CSF and aCompCor.

This study presents several limitations. Firstly, we focused only on the denoising methods implemented in fmriprep. More advanced methods have been developed recently (Agrawal et al., 2020; Shin et al., 2022; Bancelin et al., 2023) wechose to compare only the methods implemented in fmriprep, as they may be more commonly used. The same methodology can be used to compare other physiological denoising techniques. Moreover, we only evaluated the methods in their ability to remove physiological signals. The methods can perform differently in removing other types of noise. For example, ICA-AROMA has been shown to be more successful in removing head motion compared to aCompCor (Pruim et al., 2015a). The removed head motion may affect high-frequency as well as low-frequency bands, depending on the head motion type. Low-frequency head motion has been associated with rhythmic motion related to breathing (Power et al., 2019), and may be addressed by physiological regression, while random abrupt motion can be a more broadband disruption that is not captured in physiological noise time courses. In fact, these high-frequency motion signatures have also been associated with useful between-group effects (Yuan et al., 2016).

Secondly, in the absence of a ground-truth, we simply showcase a comparison of different denoising methods for revealing the age-related fcMRI differences. Indeed, head motion (both bulk and localized) differs between age groups, with older adults known to exhibit greater degrees of movement (Van Dijk et al., 2012; Pardoe et al., 2016; Savalia et al., 2017; Madan, 2018; Saccà et al., 2021), although in this study head motion is not significantly different between our young and old subjects (see Supplementary Table S2). Furthermore, we did not include head motion as regressor in our pipeline, fully recognizing the likely complex interactions between motion regression and denoising. Thus, a bigger separation between ages does not imply greater ability to detect functional differences, but may rather serve as a marker of age-related difference, which has applications in its own right. Moreover, we did not formally test BOLD fractional power change in any of the three frequency bands as a mediator of the fcMRI differences associated with the different denoising methods.

Thirdly, we did not explicitly remove head motion as we did not record head motion independently of the image acquisition. The effect of bulk and regional head motion can be profound to our ability to discern neuronally specific connectivity patterns, and will be the focus of our future work.

Lastly, hybrid approaches that combine the test methods are becoming increasingly recognized (Li et al., 2019; Van Schuerbeek et al., 2022; Weiler et al., 2022). While a thorough investigation of these potentially powerful combinations is beyond the scope of this work, we trust that in order to effectively combine different methods, it is even more important to understand the effect of each denoising method in isolation. Alternatively, a multiverse approach (Dafflon et al., 2020; Chen et al., 2023) can be used to systematically explore a diverse range of denoising methods and explore the robustness and reproducibility of findings corresponding to different methodological choices.