In order to simulate the aliasing effect of the respiratory motion in the fMRI signal acquired at a certain TR, a cosine wave oscillating at frequencies near the common rat ventilatory rate was created using Matlab (Mathworks, Natick, MA, United States). The generated signal was then down-sampled at frequencies equivalent to potential volume acquisition times (TR). The resulting waves were plotted to illustrate the effect of sampling an intruder ventilatory-driven signal at a certain TR.

The function f(TR) calculates the frequency freqS of a signal S that results from sampling an original signal S0 at a given TR.

with k = 1,2,3 … .

Supplementary Figure S1B shows an example of how a signal (e.g., motion from breathing at 60.1 breaths per minute) can be aliased by using typical TRs (i.e., TR > 0.5 s). When using a sampling frequency equal to the peak frequency S0, the resampled signal (S) is flat (frequency = 0) (see Supplementary Figure S1A middle inferior panel). In the practice, complete cancelation of the signal is often not possible, when dealing with analog devices, but by getting very close to the peak frequency, an extremely slow aliased signal can remain, which would lie out of the rs-fMRI analysis (freqS < < < 0.01 Hz, Supplementary Figure S1B). By knowing the specific rate at which the animals breathe (e.g., by imposing a respiration rate with ventilatory support) it is possible to choose an appropriate TR that minimizes the interference of the aliased signal with the spectrum of interest (e.g., 0.01–0.1 Hz).

In order to test whether aliasing occurs in the fMRI signal upon spontaneous breathing, nine adult (300–400 g) male Sprague Dawley rats were included in the study and their fMRI signal was assessed under anesthesia at three different conditions: (1) spontaneous breathing through a mask; (2) ventilated, not paralyzed; and (3) ventilated and paralyzed. All animal procedures were approved by the Animal Protection Committee of Tübingen (Regierungspräsidium Tübingen).

All animals were initially anesthetized with 5% isoflurane in chamber. In order to track changes in the blood pressure, the femoral artery was catheterized and connected to a blood pressure transducer (Biopac, Goleta, CA, United States). As in further steps the effect of a muscle relaxant was to be assessed, the femoral vein was also catheterized and connected to an infusion pump. The inner ear cavity was filled with paste to avoid air-tissue interfaces near the brain. A closed-loop heating system was used to keep the rectal temperature at 37°C during the entire experiment, ensuring stable physiological conditions. The end-tidal CO2 was also controlled throughout the experiment using a Respironics-Novametric monitor. A chest transducer was placed under the chest of the animal and constantly used to monitor breathing through a Biopac acquisition system, at 5000 samples per second.

During acquisition of the fMRI, 6 animals were kept anesthetized with 1.5–2% isoflurane, 1 animal was anesthetized with alpha-chloralose (25 mg/Kg⋅h), 1 with urethane (1 dose of 1–2 g/Kg) and 1 with medetomidine (0.05–0.5 mg/Kg⋅h), to rule out the contribution of a specific anesthetic to the observed effects. First, the not-yet-intubated animals were transferred to the 12 cm horizontal bore of a 14T MRI scanner and breathed oxygen-enhanced air (∼30% oxygen) through a mask. In the second part of the study, animals were intubated and maintained at a constant ventilatory rate that could be varied between scans. A last study was carried out in animals anesthetized with isoflurane with concurrent intravenous infusion (∼1–2 mg/Kg⋅h) of the muscle paralyzer agent Pancuronium (typical fMRI setup). All animals were euthanized right after the fMRI study based on the termination procedure from the approved protocols.

To detect the BOLD signal in the anesthetized animals, 2D and 3D Echo Planar Imaging (EPI) sequences were used, covering a field of view (FOV) of 2.24 cm × 1.92 cm × 1.92 cm with a matrix size of 56 × 48 × 32, achieving a resolution of 400 μm × 400 μm ×600 μm (all the data shown in the main figures correspond to 3D EPI, unless stated otherwise -Supplementary Figure S7 is an example of a 2D EPI-). In order to determine whether different TRs modified the frequency of the aliased signal, the TR was varied between scans (TE was kept at 9.75 ms). A phantom containing Phosphate Buffered Saline (PBS) was located over the animal’s head to compare the signal detected inside the brain with a signal experiencing similar motion in an inert tissue. A surface coil of diameter 2.2 cm was used to transmit and receive the RF signals. For anatomical registration purposes, RARE images were also acquired using a TR of 4000 ms, TE of 9.025 ms, FOV of 2.24 × 1.92 × 1.92 and matrix size of 128x128x32 (resolution of 175 μm × 150 μm × 600 μm). Table 1 summarizes the conditions and fMRI parameters used in our study.

The purpose of the analysis was to determine the influence of the respiration in the fMRI time course. Therefore, both signals (i.e., the chest–derived motion trace and the fMRI time-course), were considered for analysis.

The amplitude of the fluctuations in the BOLD signal was considered as an indicator of the level of contamination in the data. In order to detect the brain areas most vulnerable to the respiration-related interference, a power map was created in AFNI (Cox, 1996), by calculating the BOLD signal change at each voxel. The function 3dPeriodogram was used to obtain 1024 frequency components and their corresponding power in each functional scan. 3dTstat was used to average the power of the frequencies in the range of 0.005 to 0.4 Hz in each voxel, resulting in a color-coded power map. In order to compare different animals and to average maps together, all scans were registered to a template using 3dAllineate. First, the anatomical datasets were aligned, and later the registration matrix was applied to the functional scans.

In order to compare the respiratory-driven chest motion and the fMRI time-course, both signals were plotted in temporal and frequency domains. fMRI time courses were extracted by manually drawing a mask of 4 voxels from inside the brain and from a phantom placed over the animal’s head. The AFNI function 3dMaskave was used to compute the average of the 4 voxels and output a 1D file. To extract single voxel time courses, the AFNI function 3dmaskdump was implemented. The time courses were then read in matlab using dlmread. In order to create a signal that simulates the chest motion detected from the fMRI, chest movement traces were down-sampled to 1/TR for each fMRI acquisition and included in the analysis. The fast Fourier transform (fft) function of matlab was used to calculate the power spectral density (PSD) from the original chest movement trace, from the down-sampled chest trace, and from the fMRI time courses. The full width at half maximum (FWHF) of the peak aliased frequency was calculated on the enveloped chest and fMRI spectra using matlab, to assess their similarity in terms of bandwidth and differentiate between mask and ventilation conditions. A two-tailed not-paired t-test was performed to assess the differences between groups. Errors in figures represent standard error of the mean.

To quantify the degree of variability that is incorporated to a group of scans solely based on the aliasing effect (i.e., variability of the fMRI signal upon using different ventilatory rates or TRs), spectral signal similarity was evaluated for each pair of conditions (e.g., time course acquired at TR = 1 vs. TR = 1.2 s, or time course acquired from the animal ventilated at 60 bpm vs. 61 bpm) by using the corr2 function in matlab applied to the PSDs. Correlation matrices were built to visualize the effect across all different conditions. Signal amplitude (i.e., BOLD signal change) was calculated as the standard deviation within the time course. A two-tailed not-paired t-test was performed to assess the difference in averaged signal amplitude (indicative of signal contamination) between fMRI time courses from spontaneously breathing animals vs. paralyzed animals. Errors in figures represent standard error of the mean.

In order to have a measure of the spectral specificity of the respiration-induced interference at different conditions (i.e., how wide the contaminated frequency range is or how much it spreads from a center frequency in an animal breathing through mask or mechanically ventilated), a map showing the fractional amplitude of low frequency fluctuations (fALFF) at specific frequency bands was generated using the function 3dRSFC (Zou et al., 2008), to cover a bandwidth of 0.01 Hz at the peak aliased frequency and the immediate lower and higher spectral components. To further verify the effect of ventilating the animals on the width of the frequency range subjected to respiratory-driven interference, the matlab function bandpower was used to calculate the average power within a 0.01 Hz bandwidth centered at the peak aliased frequency or the immediate superior and inferior neighbors. A two-tailed paired t-test was performed to assess the difference in power between the peak and the neighbor frequencies. Errors in figures represent standard error of the mean.

Table 2 provides the p-values obtained from the different t-tests performed during analysis.

A 3D EPI-fMRI method was used to study the potential interference between the animal respiration and the rs-fMRI signal in anesthetized rats. Figure 1 illustrates the potential disruption of B0 field homogeneity due to the respiration-induced movement of the chest/neck, which results in distorted EPI images with characteristic periodic patterns in ventilated rats without muscle relaxant. It is worth noting that the periodic motion artifacts can be detected in animals under different anesthesia (e.g., medetomidine, urethane, isoflurane, and alpha-chloralose) (Supplementary Figure S2). Since muscle relaxants avoid the movement of the chest/neck, they reduce the respiration-induced B0 offset during fMRI acquisition, therefore, 3D EPI images show little distortion in paralyzed rats under ventilation (Figure 1C).

To specify the aliasing effect that occurs when the respiratory cycle is sampled at a given TR (see simulated data in Supplementary Figure S1), the respiratory signal (i.e., the chest movement) was recorded simultaneously with the rs-fMRI signal in anesthetized rats under three conditions: (1) not-ventilated, spontaneously breathing through mask, (2) ventilated, without muscle relaxant (i.e., not paralyzed), and (3) ventilated, with muscle relaxant (i.e., paralyzed). In Figure 2, the TR of the 3D EPI-fMRI was 1.2 s, i.e., 0.83 Hz. As shown in the PSD of the respiratory signal, the freely breathing rat (i.e., not ventilated -breathing through mask-) exhibits a broad frequency peak centered at ∼1 Hz (“chest”, Figure 2A). The FWHM of this broad peak in non-ventilated rats is ∼0.06 Hz ± 0.004 Hz (Figure 2E). The aliasing effect of this respiratory signal at a sampling rate of 0.83 Hz causes a slow oscillatory signal at the 0.17 Hz ± 0.03 Hz, which is observed in both, brain voxels and phantom voxels positioned at the top of the rat head (Figure 2A). A strong correlation between the peak frequency observed from the aliased (i.e., down-sampled) chest motion and the fMRI time course of both, brain and phantom voxels, suggests that the artifacts originate from the respiratory-induced B0 offset and not from brain functional dynamic signals. Similarly, the fMRI signal from rats ventilated at 1 Hz shows clear aliasing effect at 0.17 Hz, but the bandwidth of the aliased signal, assessed as the FWHM of the aliased peak, is much sharper than the observed in not-ventilated animals (Figure 2E, p-value = 7^∗10^∧-6), in both brain and phantom voxels (Figure 2B). The aliased signal is significantly reduced in paralyzed ventilated rats (Figures 2C,F).

To specify the spatial distribution of the respiration-related motion artifacts on the 3D EPI images, the mean voxel-wise power map (0.005–0.4 Hz) was calculated for each group (Figure 2D), showing that dorsal cortical areas and ventral brain regions are more sensitive to the motion artifacts. However, no generalized spatial patterns of motion artifacts can be achieved; instead, the respiratory-related motion and, hence, the interferential signal, varies across animals (Supplementary Figure S3). The significantly lower power estimate at the peak frequency in the paralyzed group (p-value = 2.9^∗10^∧-4, compared to the non-paralyzed ventilated condition) suggests a reduced aliasing effect in animals without spontaneous breathing (Figure 2D, bottom map, and Figure 2F). Figure 2F and Supplementary Figure S4 demonstrate the differences in power between the aliased center frequency and the collateral frequency bands at three conditions (mask, ventilated and ventilated paralyzed). Figure 2F shows the averaged power within a bandwidth of 0.01 Hz calculated at the peak and the neighbor frequencies across different animals, showing a significantly sharper effect (i.e., less contribution of the neighbor frequencies) in ventilated animals, compared to freely breathing rats (p-value = 9.68^∗10^∧-11). Supplementary Figure S4 provides the fALFF maps specific for the peak and neighbor frequency bands in a representative animal. The power estimates in the adjacent frequency bands are much lower when the animal is ventilated, in comparison to the free-breathing condition. This indicates that ventilation can prevent the spread of respiratory-motion aliasing effects along a broader frequency bandwidth (Figures 2E,F and Supplementary Figure S4). This is particularly important when the conventional TR at 1s is used, which causes the aliasing oscillatory effect at < 0.1 Hz to directly interfere with the low-frequency rs-fMRI signal fluctuation (Supplementary Figure S5). This result demonstrates that ventilation can be used to confine the aliasing oscillatory effect to a sharp frequency range.

The features of the aliased oscillation were evaluated in non-paralyzed animals at multiple voxels across different parts of the brain and a phantom located over the animal’s head. Different voxels exhibited different signal interference, in both, mask-breathing and non-paralyzed ventilated animals (Figures 3B,C, 4A,B), which demonstrates a highly varied effect of the respiration on the rat fMRI signal. In order to reduce the potential interference of the aliased respiration on the LFF (< 0.1 Hz), different TRs were used (0.5 s to 2 s), demonstrating that the aliasing oscillation can be shifted along the spectrum in non-ventilated animals (e.g., to be out of the < 0.1 Hz low-frequency range) (Figure 3D and Supplementary Figures S1B, S5). Similarly, to avoid the interference to the < 0.1 Hz rs-fMRI signal fluctuation, the ventilation rate could be altered in mechanically ventilated animals with fixed TR (1 s, Figure 4C and Supplementary Figure S1A), as an alternative strategy for anesthetized animal rs-fMRI studies without employing muscle relaxant agents.

The fMRI signal is an indirect measure of neuronal activation and therefore, any variability introduced by external sources is undesired and should be minimized. In order to assess the degree of variability between data acquired under different ventilatory rates or different TRs (i.e., different aliasing effects), and to determine whether the use of a paralyzer agent could reduce the differences across scans, fMRI time courses were acquired at different sampling and ventilatory conditions with and without the muscle relaxant Pancuronium (Figure 5). In ventilated rats anesthetized with constant 2% isoflurane and without muscle relaxant, the fMRI time courses acquired at different TRs exhibited a large between-scan spectral variability, compared to groups of scans performed during infusion of the paralyzer agent (Figure 5B). Similarly, the spectral correlation between scans acquired at different ventilatory rates was lower in spontaneously breathing animals (Figure 5F). The magnitude of the aliased-motion artifact, estimated as the standard deviation within the fMRI time course, was calculated across animals, revealing a significantly lower contamination of the fMRI signal in the paralyzed condition (Figure 5E, p-value = 1.13^∗10^∧-5). These experiments demonstrate that the respiration-induced interference on the fMRI signal and the subsequent inter-scan variability can be ameliorated by providing a muscle relaxant during the fMRI study. In agreement with this, the seed-based correlation analysis of rs-fMRI acquired from ventilated rats without muscle relaxant (no denoising procedure) showed highly varied correlation patterns throughout the brain, in contrast to the more specific correlation pattern observed from rats ventilated and paralyzed with a muscle relaxant, which identifies RSNs more reliably (Supplementary Figure S6), as previously reported (Williams et al., 2010; Lu et al., 2012; Hsu et al., 2016; Bajic et al., 2017).

Human cardiovascular and respiratory signals in the young adult fluctuate at 0.66–1.66 Hz and 0.16–0.33 Hz, respectively (Fleming et al., 2011). Human respiration can, therefore, be sampled if relatively short TRs are used (< 1.5 s), preventing aliasing. In contrast, the cardiovascular oscillations constitute a clear source of aliasing in human fMRI (TR needs to be ∼0.3 s to sample the human heart rate appropriately) (Hu et al., 1995; Kiviniemi et al., 2005; Feinberg et al., 2010; Hugger et al., 2011). Despite the similarity between the neurological signals of human and rodents (e.g., EEG spectrum), the heart and respiratory rates are approximately 5 times faster in the rat (5–8 Hz and 0.9–1.5 Hz, respectively), with their respiratory cycle matching the human cardiac cycle (Zhao et al., 2008; Williams et al., 2010). The spectral features of the rat cardio-respiratory signals imply the need of correction methods to avoid aliasing in the 0.01–0.1 Hz fMRI signal. Here, we have shown how the respiratory interference can be detected and shifted away from the rs-fMRI analysis by changing the ventilatory rate or the rate of fMRI sampling (i.e., TR). This is possible due to the fact that the aliased frequency depends on the frequency of the interference (i.e., the respiratory rate) and on the sampling rate (i.e., the TR), which is a consequence of the Nyquist Sampling Theorem, stating that if the rate of sampling is less than twice the signal frequency, aliasing occurs, with f_alias = | f_resp-n^∗f_sample|, n = 1, 2, …, f_sample = 1/TR (Holst, 1996; Kiviniemi et al., 2005). Future studies using fast imaging schemes (e.g., line scanning with 50 ms TR Yu et al., 2014) could be used to characterize the cardiovascular contributions to the fMRI signal in rodent studies.

In this work, we demonstrated the aliasing effect of the respiratory-driven artifacts in the rat rs-fMRI acquired as 3D EPI, which is a highly efficient sequence to acquire whole-brain fMRI signal. Nevertheless, a significant portion of the animal and human fMRI studies are performed with multi-slice (2D) EPI (Keilholz et al., 2004; Bruno et al., 2010; Magnuson et al., 2010; Williams et al., 2010; Lu et al., 2012; D’Arceuil et al., 2013; Chen et al., 2015b), which is less sensitive to physiological noise than its 3D counterpart (Lutti et al., 2013). While 2D EPI slices are acquired consecutively within ∼30 ms, with each slice sensing a different dynamic B0 field during a small portion of the respiratory cycle, 3D EPI schemes sample the whole brain slab in an additional phase encoding step, which incorporates the B0 field distortion through the whole respiratory cycle. This difference will certainly make the 3D EPI method more sensitive to the B0 field distortions due to respiration. To account for the potential differences between both imaging schemes, multi-slice 2D EPI was acquired in animals to demonstrate the presence of a similar aliased respiratory signal in both rs-fMRI paradigms (Supplementary Figure S7). Although the abovementioned aliased signals could be detected in both imaging schemes, the aliased signal interference remains but is reduced in the 2D EPI method in comparison to the 3D EPI [a more exhaustive comparison between the sensitivity of both methodologies is beyond the purpose of the current study, for a discussion on the topic see: (Goerke et al., 2005)]. While physiological noise can be regressed out retrospectively from 2D EPI datasets by slice-selective correction methods (Glover et al., 2000; Hutton et al., 2011), other strategies, e.g., ICA-based artifact removal (Griffanti et al., 2014), may be applied to 3D EPIs to diminish the noise generated by the aliased cardiorespiratory interference (Caballero-Gaudes and Reynolds, 2017).

The B0 offset induced by respiration-related motion causes dynamic B0 spatial distortion through the respiratory phasic cycles. Figures 2D, 3C, 4B show how the oscillatory signals due to motion-related artifacts dominate the fMRI signal fluctuation of the cerebellum, cortex and hypothalamic area across a frequency bandwidth (0.005–0.4Hz). Despite a major susceptibility of the dorsal, occipital and ventral regions, the distribution of the artifactual signals can vary between animals with different respiratory motion (Supplementary Figure S3), which shows the challenges of the generalized correction method. Paralysis of the ventilated animal resulted in a significantly reduced interference, as passive ventilation does not involve recruitment of accessory muscles in the chest or neck of the animal, which are the major contributors to the B0 inhomogeneities (Raj et al., 2001; Van de Moortele et al., 2002). However, other physiological variables not coupled to motion, such as the arterial CO2 signal (Wise et al., 2004), can change with each respiratory cycle and its direct vascular interference can be aliased into the rs-fMRI signal, even in paralyzed animals. The aliasing of these signals would cause a strong correlation between large vessels (Dagli et al., 1999). Besides the breathing-related physiological variables, an enveloped blood pressure signal of frequency ∼0.01–0.04 Hz (Meyer waves, dependent on sympathetic activity) may also appear and contaminate the BOLD signal near large vessels (Diehl et al., 1995; Pfurtscheller et al., 2017). Importantly these signals would not be a consequence of aliasing (e.g., generated due to insufficient sampling) but may be caused by real cardiac-induced slow fluctuations (Bhattacharyya and Lowe, 2004). The ventral area (i.e., near the circle of Willis) highlighted in the power map of Figures 3, 4 may also potentially indicate a contribution of an aliased cardiorespiratory signal in the rs-fMRI data.

Free-breathing subjects may exhibit a broad respiratory frequency bandwidth (FWHM≅0.06 Hz) along the acquisition time (Figures 2A,E and Supplementary Figures S5, S7, black and red traces of the mask condition). This produces an equally broad bandwidth of frequencies aliased into the fMRI signal (Figures 2A,E and Supplementary Figures S5, S7, blue traces of the mask condition). In 2D EPI, different slices within the volume may be acquired at different phases of the varied respiratory cycles along the scan time, which leads to mild aliasing periodic patterns that could be diminished with slice-selective denoising methods (Glover et al., 2000; Hutton et al., 2011). In 3D EPI, the B0-offset throughout the varied respiratory cycles may produce artifacts along the second phase encoding direction that is harder to be corrected using the existing regression methods. Ventilation of animals at constant rates during scans forces the respiratory cycle to remain unvaried, therefore avoiding phase difference between cycles and sharpening the spectrum of the respiratory-derived motion within functional scans (Figure 2E and Supplementary Figure 4b). By providing ventilation at specific rates, the process of identification of the aliased frequency is optimized, allowing to exclude the given spectral components with band-stop filters. Furthermore, paralyzer agents could be used to further dampen the respiratory interference.