We recruited 13 participants at hypercapnic, normocapnic and hypocapnic baselines, with further exclusion of two participants due to the absence of physiological measurements, resulting in 11 participants (25–38 years old; 2 males and 9 females) being included in the study. Participants were recruited through the Baycrest Participants Database, consisting of individuals from the Baycrest and local communities. The study was approved by the research ethics board (REB) of Baycrest, the experiments were performed with the understanding and written informed consent of each participant, according to REB guidelines.

Image acquisition was performed with a Siemens TIM Trio 3 T System (Siemens, Forchheim, Germany), which employed 32-channel phased-array head coil reception and body-coil transmission. BOLD data was acquired using a gradient-echo EPI pulse sequence (TR = 380 ms, TE = 30 ms, FA = 40°, 15 slices, 3.44 × 3.44 × 5 mm³ with 20% slices gap, 950 volumes). T1-weighted MPRAGE anatomical image was acquired (TR = 2,400 ms, TE = 2.43 ms, FOV = 256 mm, TI = 1,000 ms, readout bandwidth = 180 Hz/px, voxel size = 1 × 1 × 1 mm³). A finger oximeter built into the scanner was used to monitor heart rate, and a pressure-sensitive belt was used to monitor respiration during the BOLD scan. Data sets without coverage of the aqueduct (n = 1) and fourth ventricle (n = 4) were excluded from subsequent analyses of aqueduct and fourth ventricle signals.

We administered mixtures of O₂, CO₂ and medical air using the RespirAct^TM breathing circuit (Thornhill Research, Toronto, Canada) for all gas manipulations. With the sequential gas delivery method (Slessarev et al., 2007), the end-tidal partial pressure of O₂ (PETO₂) and CO₂ (PETCO₂) pressures were targeted by computerized and independent means. This setup allowed us to manipulate the basal PETCO₂ level of each participant precisely, without altering PETO₂. We separately targeted a normocapnic baseline (participant’s nature baseline), a hypercapnic baseline, and a hypocapnic baseline, separated by 4 mmHg CO₂ (to avoid metabolic change (Chen and Pike, 2010)). There was no change in capnia during the recording in order to avoid the possibility of CSF directly being affected by vasodilation caused by the transition between different capnias (Liu et al., 2025; Zimmermann et al., 2023). The RR was self-regulated throughout the respiratory challenges. Hypocapnia is achieved with RespirAct^TM primarily through increased breathing depths, and hypercapnia is achieved by increasing CO₂ levels in the air supply. A pseudo-randomized sequence of capnic conditions was used across different participants, with approximately two minutes between each condition. During the study, breath-by-breath CO₂ levels were recorded at a rate of 50 Hz using the RespirAct.

The FreeSurfer reconstruction was performed on the T1 anatomical data for all participants using FreeSurfer 6.0 (available at: https://surfer.nmr.mgh.harvard.edu). This reconstruction provided tissue segmentation of gray matter, white matter structures, as well as ventricles, which can then be used to delineate regions of interest.

Signal metrics were compared across capnic conditions, and all associated results are shown as Cohen’s D values, as summarized quantitatively in Supplementary Figure S3 as well as verbally in Supplementary Table S1. When comparing hypocapnia to normocapnia, the fMRI signal in vascular and CSF ROIs behaved in very similar ways, and the same is true when comparing hypercapnia to normocapnia. Relative to normocapnia, both hypocapnia and hypercapnia are associated with increased fMRI signal power and shifts in BOLD signal frequency in a band-dependent manner. The behaviour of CSF velocity, only defined for frequencies up to 0.1 Hz, also mimics that of the fMRI signals in the vascular and CSF ROIs. Moreover, GGMS signal metrics behaved very similarly to those from vascular ROIs. For visualization purposes, the power spectra shown below was computed using Welch’s method (Hamming window, 41 sample overlaps), but the statistical results are based on the Fourier spectra.

We see obvious differences across the infraslow frequency bands and non-infraslow frequency bands (Figure 4). Interestingly, basal CO₂ was not a strong driver of the differences across capnic conditions for all frequency bands. Below are the details of the effect for each band:

Band IS1: Positive associations between frequency and RR in the artery, vein, lateral ventricle, third ventricle, and GGMS.

It is also worth noting that extra-cranial sources in upper frequencies (Bands 2 and 3) do not significantly modulate neurofluid dynamics in our modeling framework. Moreover, we did not observe a direct modulatory effect of PETCO2 on neurofluid dynamics.

In this work, we use the RRF and CRF to characterize the time-dependent coordination between the ANS and CSF/vascular signals. Based on the fact that RR, HR and HRV are the most predominant mediators of the relationship between capnic condition and fMRI signal changes in neurofluid ROIs (as shown in Figure 4), we estimated the RRF linking respiratory volume per unit time to the fMRI data, as well as the CRF linking HRV to the fMRI data. As shown in Figure 5, there were distinguishable differences in the shapes of the RRF across three capnias. Further statistical analysis revealed that hypercapnia showed a higher first peak FWHM than normocapnia for all ROIs (excepting aqueduct; p_Artery = 0.03, p_vein = 0.02, p_{LateralVentricle} = 0.03, p_{ThirdVentricle} = 0.008, p_{FourthVentricle} = 0.03; D_Artery = 1.00, D_vein = 1.10, D_{LateralVentricle} = 1.00, D_{ThirdVentricle} = 1.10, D_{FourthVentricle} = 1.03) and GGMS (p_GGMS = 0.03; D_GGMS = 1.00), but not for CSF flow. The hypercapnia also showed higher first peak lags than hypocapnia for lateral and third ventricles (p_{LateralVentricle} = 0.04, p_{ThirdVentricle} = 0.02; D_{LateralVentricle} = 1.01, D_{ThirdVentricle} = 1.02).

As shown in Figure 6, there were distinguishable differences in the shapes of the CRF across three capnias. Further statistical analysis revealed that the first peak intensity for the CRF for hypocapnia was significantly higher than normocapnia for CSF flow (p_CSFFlow = 0.02; D_CSFFlow = 0.90). Moreover, normocapnia showed higher second peak intensity than hypocapnia for CSF flow (p_CSFFlow = 0.02; D_CSFFlow = 1.07).

As the CSF flow dynamic in the resting state (after the capnic effect reaches a steady state) originates primarily from vascular activity, especially arteries, given their ability to dilate and contract, ATF was estimated in order to better understand how ANS-induced vascular tone differences impact CSF dynamics. The ATF quantitatively describes how arterial fluctuations are converted to CSF flow fluctuations, a process which is directly influenced by tissue elasticity and vascular tone, as described earlier. As outlined in the biophysical model (Zhong et al., 2024), the increase in arterial signal is mainly due to the contraction of the arteries (less R₂’ decay caused by arterial-derived local-field inhomogeneity). As shown in Figure 7, ATFs from different CSF ROIs and CSF flow behave very similarly at different capnias. Further statistical analysis revealed there is no significant difference in the ATF parameters across different capnias, although qualitative differences exist in the plots for CSF ROIs.

Given that the ATF did not differ significantly across capnias, whereas the RRF and CRF did, and given that we mainly uncovered differences in signal frequency across capnias (Figures 8a–d) that are mediated by HRV(t) and RVT(t) frequencies (Figure 4), we wanted to clarify if these frequency differences in Figures 8a–d are reflected in the RRF and CRF frequencies. As shown in Figures 8e–l, frequency differences across capnias were identified in the RRF estimates associated with all signals, but not in the CRF estimates. This suggests that despite HRV(t) frequency (and not RVT(t) frequency) being the major ANS-related mediator of neurofluid signal frequency differences across capnias, the frequency governing the propagation of HRV variations to neurofluids variations remains unaffected by capnia, whereas the frequency governing the propagation of RVT variations to neurofluids variations is capnia-dependent. The frequency differences for CSF ROIs can be found in Supplementary Figure S4.

To further understand the contribution of respiratory and cardiac variability to neurofluid dynamics, we calculated mutual information measures relating HRV and RVT waveforms to neurofluid waveforms. Compared to HRV(t), RVT(t) showed higher (p_Hypocapnia = 0.03, p_Normocapnia<0.001, p_Hypercapnia = 0.02; D_Hypocapnia = 1.1, D_Normocapnia = 2.0, D_Hypercapnia = 1.1) amounts of mutual information with CSF flow (Figure 9). Nevertheless, the mutual information analysis did not reveal significant differences across capnic conditions. There was, however, a trend of increasing mutual information from hypocapnia to hypercapnia. Details for all ROIs can be found in Supplementary Figure S2.

Aside from the ANS-mediated modulation of CSF dynamics, baseline CO₂ may modulate CSF dynamics through other potential pathways. Baseline CO₂ may alter the heart and respiratory rate directly through stimulating chemoreceptors in the brainstem, independently of the ANS (Kronenberg and Drage, 1973). At a group level, while not significant, HR is higher at both hyper- and hypocapnia, while RR is highest at hypercapnia and lowest at hypocapnia. While the group comparisons did not support significance of these differences, HR and RR still differed across conditions on an individual level (Supplementary Figure S1). Such shifts in HR and RR can result in seemingly counterintuitive observations in terms of fMRI frequency. For example, a shift of the HR-related fMRI signal peak to a higher frequency can result in the aliasing of the HR peak into Band 1. However, based on our HR measurements, this potential for aliasing only applies to 3 of the participants (Supplementary Figure S1). Accordingly, based on the limited group effects observed in our experiment, baseline CO₂ modulation of neurofluid via HR and RR does not appear comparable to the effect of ANS-mediated modulation.

While CO₂ is widely used to alter vascular tone and hemodynamics (Ainslie and Duffin, 2009; Cohen et al., 2002), which are in turn tightly linked to CSF dynamics, the mechanisms through which basal CO₂ influences CSF dynamics are less well documented. Of particular relevance to the present study, CO₂ can modulate vascular tone, and hence hemodynamic fluctuations, independently of the ANS. However, it has been unclear how much basal CO₂ influences CSF dynamics solely through its biomechanical effect on CVR.

Vascular tone is conventionally defined as the baseline level of constriction in a blood vessel relative to its maximally dilated state. However, the ability of fMRI signals to fluctuate depends not only on dilatory capacity but also on constrictive capacity. Thus, in the context of BOLD based neurofluid fluctuation measurements, a more neutral vascular tone (i.e. during normocapnia) is biomechanically associated with the highest level of signal fluctuations, as was shown in our previous work (Halani et al., 2015). This is consistent with previous work by Biswal et al., who uncovered a suppression of 0.08 Hz band rs-fMRI signal fluctuations during 5% hypercapnia, which was attributed to the reduced sensitivity of hemodynamic modulation during a vasodilated state (Biswal et al., 1997). These observations can in part explain our own observations of significant positive correlations between rs-fMRI signal fluctuation amplitude and CVR, with CVR being diminished at hyper- and hypocapnia compared to normocapnia (Golestani et al., 2016). Moreover, previous work by Cohen et al. (2002) found that the BOLD peak height was inversely correlated with basal CO₂, while the BOLD response width was directly correlated with basal CO₂, implicating a fMRI frequency modulation by CO₂ as well.

As hemodynamic fluctuation amplitude reflects CVR, the biomechanical hypothesis would predict that the higher CVR at normocapnia be associated with higher CSF fluctuations than at hyper- or hypocapnia, and that the hypocapnic baseline be associated with faster CSF fluctuations (as shown in Figure 1). The results of our study, however, were very different. Instead, we found CSF as well as vascular fluctuation amplitude and frequency to both be higher at hyper- and hypocapnia (Figure 8). Our findings are also in accordance with previous findings in which it was found that ANS rather than (or in addition to) local CO₂ plays a crucial role in modulating cerebral arterial blood flow (Özbay et al., 2019). Moreover, given the ATF is unchanged across capnias (Figure 7), implying that the ability of arterial pulsations to drive CSF pulsations is independent of vascular tone and CVR, the biomechanical mechanism is unlikely to dominate the differences in the dynamics of CSF across capnias.

An alternate or complementary mechanism for the effect of CO₂ on the rs-fMRI signal is through the effect of CO₂ on neural activity. Elevated CO₂ has been associated with suppressed steady-state amplitude of the band-limited EEG Hilbert envelope amplitude in multiple EEG and MEG bands (Driver et al., 2016; Xu et al., 2011). By the same token, reduced CO₂ should be associated with enhanced EEG amplitude. In this regard, previous work suggests that neuronal activity variations can influence hemodynamics, which in turn influence CSF fluctuations as measured using BOLD rs-fMRI. Fultz et al. (2019) found that the global grey-matter BOLD signal was coupled with the CSF fMRI signal from the fourth ventricle, which was in turn coupled with low-frequency delta activity. Thus, if hypercapnia and hypocapnia act to suppress and enhance synchronized neuronal activity fluctuations, respectively, they can also enhance and suppress the CSF signal fluctuations, respectively.

In this experiment, contrary to this expectation, we did not see a reduction in CSF fluctuations at higher CO₂ levels (also illustrated in Figure 1, and shown in (as seen in Supplementary Figure S3; Supplementary Table S1)). Although we also investigated the capnia-dependence of the GGMS, which has been associated with global neuronal activity (Schölvinck et al., 2010; Huber et al., 2024), the GGMS has also been linked with systemic physiological, hemodynamic fluctuations (Power et al., 2017; Tong et al., 2013) as well as ANS activity (Özbay et al., 2019), particularly within the low-frequency range (Band 1).

Current literature has been focused on the amplitude (i.e., power) of neurofluid fluctuations (Yang et al., 2024; Picchioni et al., 2022), so it remains unclear what the frequency of these fluctuations mean. As mentioned earlier, due to the vasogenic drivers, the infraslow bands (IS1 and IS2) are most strongly associated with ANS control (Figure 4). What is the role of slow CSF flow fluctuations? While the fast, cardiac-driven pulsations are considered the primary engine for the convective bulk flow that powers the glymphatic system, the slower oscillations in CSF, such as those related to ANS modulation, might have a modulatory effect on this convective flow (Kedarasetti et al., 2022). Moreover, slow-wave flow can maximize the exchange between the CSF and interstitial fluid that enables waste removal.

A recurring finding in this work is the influence of capnic condition on neurofluid fluctuation frequency in this slow band, as shown in Figure 8. In this low-frequency range, modulating CO₂ levels appears to modulate the frequency response of the respiratory control of neurofluid fluctuations. Since our RR was maintained across all capnic conditions (Table 1), and since RR itself does not necessarily correlate with low-frequency variations in RR (which is measured by RVT), a more interesting question is “what underlies the link between frequencies of HRV(t), RVT(t) and neurofluid signal frequency?”, which remains to be clarified in future research.

One of the major limitations of this study may be the limited number of participants. However, although our dataset provides unique insight into the physiological activity during different capnic conditions, only 13 participants were included, which limits the generalizability of the results from linear mixed effect models. In addition, participants were excluded for the aqueduct (n = 1) and fourth ventricle (n = 4) ROIs, which necessitates more cautious interpretation of results from these two regions. Further, we used a relatively low spatial resolution in order to maximize the temporal sampling rate. While we believe that a high sampling rate is essential for an analysis of this nature, we realize the potential for partial-volume effects in our analysis. Additionally, due to the limited field of view, ROIS for the aqueduct and fourth ventricle analysis may not be available in all participants (also see Methods section). Finally, despite our efforts to maximize the sampling rate, our temporal resolution of 380 ms may not fully capture all cardiac harmonics.

In our study, we chose macrovascular and CSF ROIs that are distant from sites of local neural activity. In spite of this, the impact of the capnic condition on global neural activity cannot be entirely ignored (Xu et al., 2011). There is a possibility that variations in global neural activity can be transmitted into macrovasculature, which supplies and drains brain blood. In our linear-mixed effect model, these factors have not been included. In addition, several other factors, such as vasomotion and vascular tone, are difficult to monitor directly and have not been incorporated into the linear mixed effect model. It is also noteworthy that the macrovascular ROIs used in this study are based on regions with high fMRI-signal power, which may differ from real macrovascular structures based on recent biophysical models (Zhong et al., 2024).

Peripheral physiological recordings provide insight into cardiac and respiratory activity, and we used them as a surrogate for the activity of the autonomic nervous system. Numerous studies have demonstrated the link between ANS activity and cardiac and respiratory activity; however, the modulation of cardiorespiratory processes depends on far more than ANS activity alone (Gordan et al., 2015). Yet, as we included only healthy young participants and used short capnic periods, the ANS may be the most likely pathway for cardiac and respiratory variation in the present study. Nevertheless, the interplay of different physiological systems on neurofluids dynamics across different capnias should not be overlooked and should be investigated in future studies.