The impact of vasomotion on analysis of rodent fMRI data

Introduction: Small animal fMRI is an essential part of translational research in the cognitive neurosciences. Due to small dimensions and animal physiology preclinical fMRI is prone to artifacts that may lead to misinterpretation of the data. To reach unbiased translational conclusions, it is, therefore, crucial to identify potential sources of experimental noise and to develop correction methods for contributions that cannot be avoided such as physiological noise. Aim of this study was to assess origin and prevalence of hemodynamic oscillations (HDO) in preclinical fMRI in rat, as well as their impact on data analysis.

Methods: Following the development of algorithms for HDO detection and suppression, HDO prevalence in fMRI measurements was investigated for different anesthetic regimens, comprising isoflurane and medetomidine, and for both gradient echo and spin echo fMRI sequences. In addition to assessing the effect of vasodilation on HDO, it was studied if HDO have a direct neuronal correlate using local field potential (LFP) recordings. Finally, the impact of HDO on analysis of fMRI data was assessed, studying both the impact on calculation of activation maps as well as the impact on brain network analysis. Overall, 303 fMRI measurements and 32 LFP recordings were performed in 71 rats.

Results: In total, 62% of the fMRI measurements showed HDO with a frequency of (0.20 ± 0.02) Hz. This frequent occurrence indicated that HDO cannot be generally neglected in fMRI experiments. Using the developed algorithms, HDO were detected with a specificity of 95%, and removed efficiently from the signal time courses. HDO occurred brain-wide under vasoconstrictive conditions in both small and large blood vessels. Vasodilation immediately interrupted HDO, which, however, returned within 1 h under vasoconstrictive conditions. No direct neuronal correlate of HDO was observed in LFP recordings. HDO significantly impacted analysis of fMRI data, leading to altered cluster sizes and F-values for activated voxels, as well as altered brain networks, when comparing data with and without HDO.

Discussion: We therefore conclude that HDO are caused by vasomotion under certain anesthetic conditions and should be corrected during fMRI data analysis to avoid bias.

1. Introduction

Functional magnetic resonance imaging (fMRI) allows for non-invasive mapping of neuronal activity as well as investigation of functional connectivity. However, caution must be taken with interpretation of fMRI, as it is an indirect representation of neuronal activity. fMRI is based on blood oxygenation level dependent (BOLD) contrast (Ogawa et al., 1990; Buxton, 2013), as the hemodynamic response to neuronal activity causes changes in the MR signal. Therefore, in addition to technical noise and motion, physiological noise including heartbeat, respiration, arterial carbon dioxide (CO₂) concentration and vasomotion can distort the functional signal (Murphy et al., 2013). This is crucial for functional network analysis, since any non-neuronal related process that affects at least one brain region impacts analysis of functional connectivity. Furthermore, non-neuronal signals may affect calculation of BOLD activation maps, since models used for fMRI data analysis with the general linear model (GLM) usually ignore physiological noise. Especially in small animal fMRI, physiological noise must be studied carefully, since measurements are often acquired under anesthesia, which alter physiological processes (Gao et al., 2017; Paasonen et al., 2018; van Alst et al., 2019) compared to the awake condition. Small animal fMRI is valuable for translational research, since it offers numerous options such as studying specific disease models (Tenney et al., 2003; Tristao Pereira et al., 2021; Wachsmuth et al., 2021), pharmacological (Haensel et al., 2015) and surgical interventions (Chen et al., 2019) as well as combining MRI with electrophysiological (Pan et al., 2010; Kosten et al., 2022) and optical recordings (Schulz et al., 2012; Albers et al., 2018; Lake et al., 2020; Ioanas et al., 2022). To avoid incorrect translational conclusions, it is crucial to investigate physiological noise and to develop correction methods for commonly performed preclinical experiments.

Although fMRI can be conducted on awake animals, small animal fMRI studies are often performed under anesthesia to reduce motion and prevent stress for the animal. A variety of anesthetic drugs is used for fMRI studies, including α-chloralose, isoflurane (ISO), ketamine-xylazine, propofol, urethane or medetomidine/dexmedetomidine (MED/DMED). Frequently, modifications of an anesthetic protocol originally proposed by Weber et al. (2006) are used, whereby animal preparation is performed under an inhalation anesthetic followed by fMRI measurements under MED sedation. This protocol offers a number of advantages, since application is non-invasive, can be maintained over extended durations and leads to temporally stable neuronal activity and functional connectivity in rats (Weber et al., 2006; Sirmpilatze et al., 2019). Using a similar anesthesia protocol (animal preparation under ISO, fMRI under DMED), Magnuson et al. (2010, 2014) detected signal oscillations with a frequency of approximately 0.2 Hz in CBV and in BOLD weighted fMRI measurements. The authors hypothesized that these oscillations resulted either from infraslow neuronal activity or from modulations in vascular tone. Close examination of these oscillations is important, because they may represent a physiological noise source that distorts fMRI data analysis. We refer to these oscillations in the following as hemodynamic oscillations (HDO), since both neuronal activity and vascular tone model the hemodynamic system.

This study aimed to investigate origin and prevalence of HDO in fMRI as well as their impact on data analysis (Figure 1). For this purpose, we performed fMRI with two different MR contrasts under ISO and MED anesthesia. We developed algorithms for detection and suppression of HDO, and investigated the effect of vasodilation on HDO. Further, we studied whether HDO are physiological noise or whether they have a direct neuronal correlate, using local field potential (LFP) recordings. Finally, we investigated the impact of HDO on fMRI data analysis using the GLM (Friston et al., 1995) as well as brain network analysis using network based statistics (NBS) (Zalesky et al., 2010).

FIGURE 1

Figure 1. Study design. fMRI was performed on Fischer rats. Algorithms for detection and suppression of HDO were developed. Afterward, HDO prevalence was investigated using the detection algorithm. In detail, it was examined in which brain regions and after which time HDO occurred after switching from ISO anesthesia to MED sedation. Further, the effect of vasodilation was studied and it was examined whether HDO occur in both SE-EPI and GE-EPI scans. Subsequently, the impact of HDO to GLM-based analysis and NBS analysis of fMRI data was examined using the developed algorithms. Additionally, LFP signals were recorded with implanted electrodes for examination whether HDO have a direct neuronal correlate.

2. Materials and methods

All experiments were conducted according to the German Animal Welfare Act and were approved by the State Agency for Nature, Environment and Consumer Protection of North Rhine-Westphalia, Germany (LANUV approval IDs: 87-51.04.2010.A274, 84-02.04.2015.A427, 84-02.04.2016.A135, and 81-02.04.2018.A426). Experiments were performed with 71 adult Fischer rats (>3 month). A total of 70 animals were female and had a weight of (184 ± 14) g, one animal was male and weighed 261 g. Rats were housed in groups of two to five animals under a regular light/dark schedule (12/12 h) with food and water ad libitum.

2.1. Animal handling and data overview

Respiratory gas was 1 L/min O₂ for anesthesia induction and during LFP experiments. For MRI, a gas mixture of 0.25 L/min O₂ and 0.75 L/min air was used. Animal preparation and LFP electrode implantations were performed under ISO anesthesia (5% induction, 2–3% maintenance). Long-term ISO experiments were performed under 1.2% ISO anesthesia. Other experiments were executed after switching to MED sedation initiated by a subcutaneous bolus injection of 0.04 mg/kg followed by a continuous infusion of 0.05 mg/kg/h. After MED bolus, ISO was discontinued within 10–15 min and experiments were started at least 40 min after bolus injection, because after this waiting period a stable physiology can be expected (Weber et al., 2006; Sommers et al., 2009; Amirmohseni et al., 2016; Sirmpilatze et al., 2019; Kint et al., 2020). A total of 30 min prior to stereotaxic surgery for electrode implantation, the animals received analgesia [Metacam (1 mg/kg s.c.) or Metamizol (100 mg/kg s.c.)]. A total of 56 animals were ventilated (MRI-1 Ventilator, CWE, Inc., Ardmore, Pennsylvania, United States). A total of 24 of them received a muscle relaxant [Pancuronium (2 mg/kg bolus followed by continuous injection of 1.5 mg/kg/h) or Atracurium (5 mg/kg bolus followed by continuous injection of 5 mg/kg/h)]. For ventilated animals, the end-expiratory CO₂ was continuously monitored using a CO₂ analyzer (Micro CapStar End-Tidal CO₂ Analyzer, CWE, Inc., Ardmore, Pennsylvania, United States). Respiration rate of ventilated animals was set to 53, 57, or 60 breaths per minute (bpm). The respiratory minute volume was between 95 and 127 mL/min. Respiration rate of spontaneously breathing animals was (56 ± 5) bpm and (68 ± 11) bpm for LFP and fMRI, respectively. Animals were placed on a temperature-controlled bed and fixated with bite and ear bars. Rectal temperature was kept at (36.7 ± 0.4) °C. Functional echo-planar imaging (EPI) and LFP data from seven experiments were investigated: (1) For up to 4 h, 5 or 10 min-long gradient echo-EPI (GE-EPI) measurements were performed repeatedly under 1.2% ISO anesthesia. (2) For up to 8 h short GE-EPI scans were acquired repeatedly under MED sedation. (3) Spin echo-EPI (SE-EPI) measurements were recorded under MED sedation. GE-EPI scans were recorded during execution of (4) CO₂ or (5) ISO challenges. (6) GE-EPI measurements were recorded during electrical paw stimulation. (7) LFP recordings were conducted. A total of 10% of the data were included retrospectively from previous studies (Albers et al., 2019; van Alst et al., 2019; Lambers et al., 2022). Therefore, temporal resolution of MR measurements varied. For characterization of oscillation frequency, optimization and evaluation of the oscillation detection algorithm, a training and test dataset was generated: For both training and test data, 240 time courses each were randomly selected from 2,254 time courses of high temporal resolution GE-EPI data. As a previous study (Lambers et al., 2020) implied that hemodynamics were not gender-specific, measurements of both sexes were combined. Grouping information is provided in Supplementary Table 1.

2.2. Electrical paw stimulation and vasodilation challenges

For electrical paw stimulation, two electrodes were inserted into one forepaw and 1-ms pulses were applied at a frequency of 9 Hz with an amplitude of 1 or 1.5 mA according to block paradigms. Further details are given in the respective figure legends. Approximately 50% of GE-EPI scans used for algorithm optimization and investigation of oscillation prevalence were recorded during electrical paw stimulation, the other 50% were recorded without stimulation (resting state). However, both resting state and stimulation fMRI data were pooled for analysis. To examine the effect of HDO on analysis of fMRI data, stimulation was applied with the following paradigm: 5-s stimulation, 25-s rest, 20 repetitions. For CO₂ challenge experiments, the mixture of inspiratory gas was alternately set to standard (0.25 L/min O₂, 0.75 L/min air, and 0.00 L/min CO₂) and 5% CO₂ (0.20 L/min O₂, 0.75 L/min air, and 0.05 L/min CO₂) according to the following paradigm: 3-min standard, 3-min 5% CO₂, 4-min standard. All animals undergoing CO₂ challenge were ventilated to enable recording of CO₂ levels in expiratory air. Starting point of the CO₂ challenge was defined at the moment when the expired CO₂ level exceeded baseline by more than 50%. End point of the CO₂ challenge was defined as the moment when the expired CO₂ level fell below 50% above baseline. Only datasets were evaluated in which the CO₂ challenge was at least 2.5 min long. Additionally, ISO challenge experiments were conducted using the following paradigm: 5-min no ISO, 10-min ISO (0.8, 1.5, or 2%), 10-min no ISO.

2.3. Functional magnetic resonance imaging experiments

MR measurements were acquired using a 9.4 T Bruker Biospec 94/20 small animal scanner equipped with a 720 mT/m gradient system (Bruker BioSpin GmbH, Ettlingen, Germany) and different combinations of transmit coils and receive-only coils. First, an anatomical image was acquired to identify the slice position for fMRI measurements: 2D Turbo spin echo sequence (RARE), TR/TE_Eff = 2,000/50 ms, RARE factor 8, Matrix 256 × 256, field of view (FOV) 28 mm × 26 mm, slice thickness 1.2 mm, 9–12 contiguous slices. Subsequently, B₀ homogenization was performed using the MAPSHIM Bruker routine. Functional GE-EPI and SE-EPI measurements were performed using a single-shot EPI sequence (same FOV as anatomy, Matrix: 80 × 80, 1.2 mm slice thickness; GE-EPI: TE 18 ms, TR 100/125/1,000 ms, flip angle 18/21/60 or 65°, 1/2/9 or 12 slices; SE-EPI: TE 35.9 ms, TR 250 ms, flip angle 90°, 3 slices). One fMRI measurement lasted between 5 and 25 min (Supplementary Table 1).

2.4. Functional magnetic resonance imaging data processing

This section describes the processing of all fMRI data with exception of the data analyzed using the GLM or NBS (see Section “2.5. Functional magnetic resonance imaging data analysis using the GLM or NBS”). Data processing was done using MATLAB (Release 2021b, The MathWorks, Inc., Natick, Massachusetts, United States). The first 5 s of each measurement were discarded to avoid pre-steady-state artifacts. A template of 14 regions based on a rat brain atlas (Paxinos, 1997) was registered to the data (Supplementary Figure 1). If the measurements contained only one slice of the template, seven regions were registered. These regions were selected since they were not expected to contain strong susceptibility artifacts, which would affect comparison of SE-EPI and GE-EPI data. For vasodilation challenge experiments, only the right forelimb region of the primary somatosensory cortex (S1Fl) was examined. Signals from voxels located in one region were summed up, downsampled to a temporal resolution of 1 s and normalized to their mean. For determination of the oscillation frequency, no downsampling was performed, as a higher temporal resolution was required for the fitting procedure.

2.5. Functional magnetic resonance imaging data analysis using the GLM or NBS

Data analysis using both, GLM or NBS was performed on data recorded for 10 min during electrical paw stimulation. For preprocessing, the first five iterations were removed, data were realigned using SPM12¹ and brains were manually masked using MRIcroGL² or MagnAN (BioCom, Uttenreuth). Subsequently, analysis was performed using GLM or NBS. GLM-based analysis was done with SPM 12 by using the 3rd order canonical basis set: convolution of the stimulation paradigm and a rat HRF (Lambers et al., 2020) and its derivatives were used as regressors. For each dataset, two BOLD activation maps were calculated using an F-test (p < 0.05, family wise error correction): before and after voxel-wise application of the HDO suppression algorithm. Effect of HDO on GLM-based analysis was only examined on datasets with more than 12 and a maximum of 100 activated voxels in the activated S1Fl (before suppression algorithm application). Data were divided into non-oscillatory and oscillatory by analyzing time courses of the activated S1Fl with the HDO detection algorithm. The evaluation of the algorithm was reviewed using manual classifications. Two datasets were excluded from analysis due to large differences between manual and automated classification.

Additionally, NBS analysis was conducted with MagnAn. In most cases, stimulation was applied to the left paw. Scans were mirrored, if the right paw was stimulated. Data were smoothed (gaussian kernel, 3 × 3 pixel, full width half max 0.6 mm) and lowpass filtered at 0.3 Hz (Fourier filter). Subsequently, regression of the global mean signal was performed and a template of 28 brain regions based on the Paxinos and Watson rat brain atlas (Paxinos, 1997) was registered to the data (Supplementary Table 2). Non-registerable datasets (e.g., measurements that did not cover both hemispheres) were excluded from analysis. If an animal was measured twice, registration was done on the first measured scan. Data were divided into non-oscillatory and oscillatory by analyzing time courses of all 28 brain regions with the HDO detection algorithm. For each dataset, Pearson correlation coefficients were calculated between pairs of brain regions. Resulting 28 × 28 correlation matrices represented undirected, weighted functional connections between brain regions. Only positive correlation coefficients were considered. Per dataset, correlation matrices were calculated twice (before and after voxel-wise application of the HDO suppression algorithm). Subsequently, mean correlation coefficients were calculated group wise. For this purpose, Pearson’s r-values were converted to Fisher’s z-values to provide normal distribution, averaged, and converted back to Pearson’s r-values. Resulting connectivity matrices (Supplementary Figure 2) clearly showed interhemispheric connections, indicating that the analysis detected meaningful functional connections. Subsequently, for the 140 strongest connections that appeared in the averaged matrices, NBS was used to find differences in functional connectivity: a paired t-test was performed between correlation matrices of HDO-suppressed und unsuppressed data and the largest components of significantly modulated connections (p < 0.05) were retained. Subsequently, permutation testing was applied with 1,000 iterations for correcting the t-test for multiple testing. Finally, mean differences of correlation coefficients (MDCC) were calculated separately for oscillating and non-oscillating data. For this purpose, differences of correlation matrices (suppressed--unsuppressed) were calculated for each significantly different connection and averaged over datasets. For visualizing results of network analysis, the software Cytoscape was used.³

2.6. Hemodynamic oscillations frequency fitting

To determine the oscillation frequency for all 240 time courses of the test dataset, a Fast Fourier Transformation (FFT) was calculated. A total of 60 of the resulting spectra showed a peak between 0.08 and 0.8 Hz. Corresponding time courses were classified as oscillating and the oscillation frequency was determined. For this purpose, a function consisting of an exponential and a Gaussian term (exponential-Gaussian function, EGF) was fitted to the spectra across the interval from 0.01 to 1 Hz.

$EGF\left(f\right)=a+b\cdot e^{-c\cdot f}+d\cdot e^{-\frac{{(f-f_o)}^2}{2{\mathrm{\sigma }}^2}}.$

The parameter f_o represented the frequency of the HDO while σ was the width of the peak in the spectrum caused by the oscillations. Frequencies and peak widths determined using the test data were averaged and used for optimization of the detection algorithm and the suppression algorithm. Parameters from two EGFs were excluded from averaging since one of the fitted parameters was unphysiological (negative peak amplitude or peak position).

2.7. Hemodynamic oscillations detection and suppression algorithms

In order to systematically investigate the prevalence of HDO and their effect on the analysis of fMRI data, we developed two algorithms using MATLAB. Both require time courses extracted from fMRI data as input and are online available.⁴ The first algorithm detects oscillations and delivers the output whether the inserted time course contains HDO. The second algorithm suppresses HDOs from the inserted time course.

The HDO suppression algorithm performed double filtering of the time courses with a finite impulse response (FIR) bandstop filter. Using the averaged oscillation frequency f_o and peak width σ the passband of the FIR filter was set to [f_o−3σ,f_o + 3σ] while the interval of the stopband was set to [f_o−σ,f_o + σ]. The filter attenuation was adaptively adjusted. For this purpose, the FFT of the time course was calculated, the interval [f_o−4σ,f_o + 4σ] of the spectrum was extracted and downsampled to a spectral resolution of half the peak width σ. The ratio of maximum and baseline of the downsampled spectrum was used for filter attenuation.

The HDO detection algorithm examined the FFT of the inserted time courses as well as the autocorrelation functions, since both are different for oscillating data compared to non-oscillating time courses: The FFT of oscillating data showed a peak while the autocorrelation function showed a periodic pattern which were absent for non-oscillating data. The detection algorithm first calculated the FFT of the extracted MR time courses. The interval [f_o−4σ,f_o + 4σ] was extracted from the resulting spectrum using the averaged oscillation frequency f_o and peak width σ. Baseline drifts were removed from the extracted interval using the function “detrend” and the detrended spectrum was downsampled to a spectral resolution of half the peak width σ. Additionally, the autocorrelation function of the time course was calculated for 15 time lags between 1 and 15 s using the function “autocorr.” Baseline drifts of the autocorrelation function were corrected using the function “detrend.” Subsequently, the detrended autocorrelation was averaged by calculating the mean of points that were 5 s apart from each other. This was done because the averaged autocorrelations of oscillatory and non-oscillatory time courses differed: Values of averaged autocorrelation functions of non-oscillatory time courses were around zero. Values for oscillatory time courses were non-zero, since their autocorrelation showed a periodic pattern with a period of 5 s. Both the averaged autocorrelation function, and the averaged spectrum, were transferred to an artificial neural network (ANN), which classified the data as non-oscillatory or oscillatory. Training of the ANN was done using 226 manually classified time courses of the training dataset (240 time courses), separating oscillating from non-oscillating time courses. Subsequently, the performance of the ANN was evaluated using 231 manually classified time courses of the test dataset (240 separate time courses, Supplementary Figure 3). All MR time courses investigated in this study were screened for HDO using the trained ANN. The developed HDO detection algorithm can be applied to time courses extracted from MR data recorded with a TR of 1 s or less.

2.8. Local field potential recordings

For LFP recordings, two custom made silver electrodes (250 μm dimeter) were implanted in the brain and glued to the skull using dental cement. One electrode was positioned in the S1Fl. The other electrode was used as reference and implanted in the cerebellum. LFP signals were recorded as difference between electric currents detected by the electrode in the S1Fl and the reference electrode using a differential amplifier (DPA-2FX, NPI Electronics, Tamm, Germany). Signals were recorded at 2 kHz using a multifunction data acquisition device (PCIe-6363, National Instruments, Austin, Texas, United States) and a custom-written LabView script (National Instruments). LFP recordings were performed for 3 min during stimulation which evoked neuronal activation in S1Fl. Only LFP recordings with visible response to stimulation were evaluated. LFP data were normalized to their baseline and the FFT was calculated using MATLAB. Time courses and spectra were manually analyzed with respect to a neuronal correlate of HDO, as no HDO detection algorithm was available for LFP data.

2.9. Data representation

Results in the text are represented as mean ± standard deviation. Data clusters are displayed as overlay of bee swarm plots and boxplots. In boxplots, central marks represent the median, boxes include the 25th and 75th percentiles, whiskers extend to the last data points within the 1.5× interquartile and outliers are indicated as crosses. Spectra shown within one figure are scaled equally.

3. Results

We investigated HDO, periodic oscillations with a period of approximately 5 s which were detected in signal time courses of fMRI measurements of anesthetized rats. These oscillations may represent a physiological source of noise for functional imaging. In total, 62% of all 303 fMRI measurements analyzed in this study showed HDO (Figure 2). For detailed investigation, we developed algorithms for detecting and suppressing HDO. Using these algorithms, we investigated HDO prevalence and assessed their impact on analysis of fMRI data.

FIGURE 2

Figure 2. Appearance and occurrence of HDO. Periodic oscillations in fMRI time courses were investigated in 303 measurements. In total, 62% of the measurements showed HDO which had a period of approximately 5 s. Exemplary time courses without and with HDO for resting state data are shown in the (upper row) and data recorded during sensory stimulation in the (middle row). Stimulation periods are indicated by black bars. Enlarged segments of one paradigm repetition (bottom row) showed that HDO (maxima indicated by black arrows) were triggered by the BOLD response (red arrows).

3.1. Hemodynamic oscillations were detected at 0.2 Hz

First, we investigated the frequency domain of MR data using FFT and continuous wavelet transformation (CWT, Supplementary Figure 4). Spectra of datasets with HDO showed increased spectral density around 0.2 Hz. Additionally, respiration often caused a peak at approximately 1 Hz, yet this peak was not related to the occurrence of HDO. No other conspicuous peaks occurred in the spectra, suggesting that HDO were not directly driven by cardiac, respiratory, myogenic, neurogenic or endothelial oscillators. We characterized the HDO frequency using a fitting procedure with an EGF (Figure 3A). An HDO frequency f_o of (0.20 ± 0.02) Hz with an average peak width σ of (0.025 ± 0.015) Hz was obtained (Figure 3B). Both HDO detection and suppression algorithms were adjusted using the determined frequency and peak width. Subsequently, the HDO suppression algorithm eliminated HDO without substantially affecting other signal components as exemplary shown in Figures 3C, D. Finally, the neural network of the detection algorithm was trained using the manually classified training data. After training, the algorithm performance was evaluated using the manually classified test data. The detection algorithm had a specificity of 95% (Supplementary Figure 3).

FIGURE 3

Figure 3. Hemodynamic oscillations frequency and assessment of suppression algorithm performance. (A) An exponential-gaussian function (EGF, blue) was fitted to FFT (gray) of an exemplary GE-EPI time course. Note that respiration caused a peak at 0.9 Hz. (B) HDO frequency f_o and peak width σ, resulting from fitting 58 spectra. Representative (C) non-oscillatory and (D) oscillatory datasets: time courses (left), spectra (middle) and autocorrelation functions (right) before (gray) and after (blue) application of the HDO suppression algorithm. In contrast to the non-oscillatory dataset, the spectrum of the oscillating dataset showed a peak and the autocorrelation function showed a periodic pattern. The suppression algorithm eliminated HDO without substantially affecting non-oscillating signals.

3.2. Vasodilation interrupted hemodynamic oscillations

CO₂ and ISO challenges were used to investigate whether short-time vasodilation interrupts HDO. Time courses were extracted from S1Fl, and three intervals (before, during and after challenge) were examined (Figure 4). In 15 vasodilation challenge scans (CO₂: n = 6, ISO: n = 9), time intervals preceding the challenge showed HDO. In all time courses, the MR signal increased during the challenge and subsequently returned to baseline, indicating that the challenge caused vasodilation. No HDO were detected during challenge in any of these measurements. We therefore concluded that vasodilation interrupted HDO. In 50% of CO₂ challenge measurements, oscillations returned immediately after end of challenge. In 90% of ISO challenge scans, oscillations were detected 5 min after challenge had ended. Furthermore, it was investigated whether HDO occur during long-term ISO anesthesia, which is known to cause permanent vasodilation. Ten animals were measured for up to 4 h under ISO anesthesia without MED. From these, in a total of 56 GE-EPI scans 14 different brain regions were examined: Only 14 of the resulting 784 region-specific time courses were classified as oscillatory by the detection algorithm. This corresponded to 1.8%, which is considerably below the false positive rate of the detection algorithm. Thus, we concluded that the occurrence of HDO under long-term ISO anesthesia, i.e., permanent vasodilation, was negligible. In an analysis, described in Supplementary Figure 5, we found that short-term vasodilation induced by stimulation of neuronal activation also interrupted HDO briefly, resulting in HDO being stimulus-locked. This effect is consistent with two studies by Sirotin and Das (2009) and Sirotin et al. (2012), which reported stimulation-locked oscillations in hemodynamic signals.

FIGURE 4

Figure 4. Vasodilation interrupts HDO. Vasodilation was induced by administration of (A) CO₂ and (B) ISO. First row: exemplary time courses recorded during CO₂ and ISO challenge. Three segments were analyzed using the HDO detection algorithm (gray background): before, during and after the respective challenge (black bars). Time course segments classified as oscillating were colored blue. Second and third row show spectra and autocorrelation functions, respectively, of the corresponding sections above.

3.3. Brain-wide hemodynamic oscillations started within 1 h after end of ISO anesthesia

Next, we investigated after which duration HDO started when switching from ISO anesthesia to MED sedation. We further assessed whether HDO prevalence changed with prolonged MED sedation. For this purpose, the HDO prevalence in 135 GE-EPI scans from 23 animals over a scan time period of maximum 8 h was examined. In 85% of the measurements, at least one time course of the investigated brain regions was oscillatory. Data acquisition started approximately 1 h after switching anesthesia. At this time, ten of 14 animals showed HDO. HDO prevalence increased over time (Figure 5A): Around 20% of early recorded (1–3 h after switching anesthesia) time courses oscillated. For late measurements (4 h or more), the portion of oscillating time courses increased to approximately 50%. HDO prevalence maps showed that HDO occurred in all examined regions (Figure 5B). After assessing HDO prevalence using GE-EPI measurements, HDO occurrence in SE-EPI data was investigated. HDO were detected in both GE-EPI and SE-EPI scans (Figure 6). To assess HDO prevalence in SE-EPI measurements, 27 SE-EPI scans from seven animals were examined. SE-EPI scans of five animals showed HDO. In total, at least one oscillating time course was found in 63% of the SE-EPI measurements. In SE-EPI data, HDO were observed only in cortical regions (Supplementary Figure 6).

FIGURE 5

Figure 5. Hemodynamic oscillations prevalence increases over time. HDO prevalence was assessed in 135 GE-EPI scans of 23 animals measured after switching from ISO to MED anesthesia. (A) Time line illustrating relative prevalence of HDO. For each animal, all measurements acquired within 1 h were combined and the percentage of oscillating signal time courses was plotted. Not every animal was measured at all time points. Measurements recorded after more than 6 h were displayed together. (B) HDO prevalence maps separated into early (recorded 1–3 h after MED start) and late (recorded 4 h or later) measurements. Color codes relative number of scans with HDO in each of the 14 brain regions.

FIGURE 6

Figure 6. Hemodynamic oscillations are detected with both GE-EPI and SE-EPI. Exemplary oscillating (A) SE-EPI and (B) GE-EPI scans. Time courses (left), spectra (middle) and autocorrelation functions (right). Both datasets were recorded in the same animal consecutively without delay. Despite lower temporal signal to noise ratio (tSNR) of SE-EPI scans, oscillations were similarly detected in both SE-EPI and GE-EPI scans.

3.4. Hemodynamic oscillations significantly impact analysis of fMRI data

Before we evaluated the impact of HDO on fMRI analysis, we assessed whether HDO are directly driven by neuronal activity. To this end, LFP recordings were performed in S1Fl, and scrutinized for electric activity patterns with frequencies around 0.2 Hz. A total of 32 LFP recordings from six animals were examined. The absence of a prominent peak around 0.2 Hz in the spectra indicated that no oscillations were present (Figure 7). The according absence of a direct neuronal correlate of HDO suggested that HDO may represent a physiological noise source for BOLD fMRI measurements. To evaluate the influence of this physiological noise on the analysis of fMRI data, BOLD datasets, recorded during electrical paw stimulation, were examined using two methods. First, statistical analysis was performed for calculation of BOLD activation maps using the GLM. Secondly, brain network analysis was performed using NBS. Both analyses were conducted before and after voxel-wise HDO removal using the suppression algorithm.

FIGURE 7

Figure 7. Hemodynamic oscillations do not have a direct neuronal correlate. (A) Exemplary LFP recording, acquired upon sensory stimulation with 10-s long stimulation periods (black bars). (B) Individual spectra (light blue) and their average (dark blue). The absence of a prominent peak around 0.2 Hz indicated that no oscillations were present.

Statistical analysis was performed on 14 fMRI datasets (14 animals) using the GLM with the 3rd order canonical basis set. Seven datasets each were classified as oscillatory or non-oscillatory. BOLD activation maps showed larger clusters of activated voxels and larger maximum F-values after applying the suppression algorithm, as compared to the original data (Figure 8). Cluster sizes and maximum F-values obtained for HDO suppressed data were normalized to values extracted from original data. Resulting relative cluster sizes were 1.13 ± 0.08 and 1.27 ± 0.25 for non-oscillating and oscillating data, respectively. Similarly, for F-values, the suppression algorithm had a larger impact on oscillating data, with relative F-values of 1.21 ± 0.06 and 1.63 ± 0.48 for non-oscillating and oscillating data, respectively. Statistical comparisons of the normalized values (U-test, Bonferroni corrected) showed that the difference between oscillating and non-oscillating data was significant for normalized F-values (p = 0.008). We therefore concluded that HDO significantly affected the GLM analysis of BOLD activation.

FIGURE 8

Figure 8. Application of the suppression algorithm improves GLM-based analysis of fMRI data. GLM-based analysis was performed on 14 fMRI datasets before (gray) and after HDO suppression (blue). (A) Periods (5 min) of BOLD time courses and (B) activation maps of one exemplary oscillating dataset. Duration of electrical paw stimulation were highlighted by black bars. (C) Cluster sizes and (D) maximum F-values resulting from GLM-based analysis of BOLD activation are shown separately for data with and without HDO (n = 7, each). Data points for identical measurements are connected by lines. Bold lines indicate the dataset shown in panels (A,B). U-tests were used to test for significant differences between relative cluster sizes and relative F-values of data with and without HDO.

Additionally, brain network analysis was performed on 18 GE-EPI fMRI datasets recorded in 16 animals. After data pre-processing (including lowpass filtering at 0.3 Hz and registration of a template with 28 brain regions) analysis was divided into five steps. First, nine datasets each were classified as non-oscillatory or oscillatory. In non-oscillating data sets, no oscillations were detected in any region. For oscillating data, oscillations were found in at least five brain regions. Secondly, we examined which of the brain regions showed frequent oscillations (oscillatory nodes). In total, 12 oscillatory nodes were detected, located in association cortex, motor cortex, sensory cortex, and thalamus (Supplementary Table 2). Third, for each dataset two correlation matrices were calculated (before and after the application of the HDO suppression algorithm, Supplementary Figure 2). Fourth, for both groups (oscillating and non-oscillating), correlation matrices of the HDO-unsuppressed data were compared with those of the suppressed data using NBS. Significant differences were found for both oscillating and non-oscillating data (p < 0.001 and p = 0.007, respectively, Figures 9A–C). In oscillating data, 44 correlation coefficients were significantly different. For 39 of these coefficients, at least one of the connected regions was an oscillatory node (Figures 9A, B), suggesting that differences were directly related to HDO suppression. Finally, differences in correlation coefficients were examined. Absolute MDCC of non-oscillating data were smaller when compared to oscillating data (0.03 ± 0.01 and 0.05 ± 0.02, respectively, Figure 9D). This difference was highly significant (U-test, p < 0.000001). The strong impact of the suppression algorithm on oscillatory data indicated that HDO strongly affected the network analysis of fMRI data.

FIGURE 9

Figure 9. Hemodynamic oscillations elimination in fMRI data improves NBS network analysis. Brain network analysis was performed before and after application of the HDO suppression algorithm and compared using NBS. Significantly different links were displayed in circular layout, separately for (A) non-oscillating and (B) oscillating data (n = 9, each). Line thickness is proportional to absolute mean differences of correlation coefficients (MDCC). Nodes where HDO were frequently detected were highlighted (bold frame). (C) Brain regions were arranged and color coded by affiliation to six anatomical groups. (D) Absolute MDCC of oscillating data were significantly larger than for non-oscillating data (U-test).

4. Discussion

In this study, we explored the prevalence of HDO in rat fMRI measurements and their impact on analysis of fMRI data. In total 62% of all 303 fMRI measurements analyzed in this study showed HDO. This number is not a precise, general estimation of probability for the occurrence of HDO in preclinical fMRI experiments, because the measurements were conducted under widely established, specific experimental conditions, suspected to induce or inhibit HDO. However, their frequent occurrence clearly showed that HDO are relevant for functional MRI and need to be examined in more detail. For this purpose, we developed algorithms for HDO detection and suppression and evaluated fMRI measurements recorded with GE-EPI or SE-EPI sequences as well as LFP recordings.

4.1. Vasomotion leads to hemodynamic oscillations

We detected HDO in anesthetized rats. Several studies reported similar oscillatory signals in measurements sensitive to blood flow or arterial vessel diameter in various species such as mice (Drew et al., 2010; Mateo et al., 2017; Winder et al., 2017; Fan et al., 2020), rats (Fujii et al., 1990; Mayhew et al., 1996; Kleinfeld et al., 1998), rabbits (Hundley et al., 1988), cats (Rivadulla et al., 2011) and humans (Mitra et al., 1997; Obrig et al., 2000; Rayshubskiy et al., 2014; Liu et al., 2022), and identified vasomotion as their origin. Vasomotion refers to rhythmic oscillations of blood vessel diameter, which causes oscillations of blood flow (Aalkjaer et al., 2011). However, the physiological sources of vasomotion are still not fully understood. We also considered vasomotion as origin of the detected oscillations. Analysis of hemodynamic recordings often involves dividing signals into five frequency bands, attributed to cardiac, respiratory, myogenic, neurogenic and endothelial oscillators, respectively (Li et al., 2006; Aleksandrin et al., 2018). We divided our signals into corresponding frequency bands and found that these oscillators do not directly drive vasomotion-induced HDO. However, we cannot completely exclude that these oscillators interact with vasomotion and thus influence fMRI measurements.

4.2. Anesthesia exerts diverse effects on vasomotion

Most preclinical fMRI experiments use anesthesia to minimize motion and imaging artifacts. Here a well accepted anesthesia protocol was used (ISO followed by MED/DMED), for which Magnuson et al. (2010, 2014) reported the occurrence of oscillations. We intended to investigate the prevalence of these oscillations (referred as HDO) in fMRI data and identified oscillations with a frequency around 0.2 Hz in MED sedated animals. In contrast, studies applying other anesthesia or measuring awake subjects often reported a vasomotion frequency around 0.1 Hz (Fujii et al., 1990; Mayhew et al., 1996; Mitra et al., 1997; Kleinfeld et al., 1998; Obrig et al., 2000; Rivadulla et al., 2011; Rayshubskiy et al., 2014; Mateo et al., 2017; Fan et al., 2020; Liu et al., 2022). The higher oscillation frequency under MED sedation can be explained by the fact that MED is a vasoconstrictor (Sinclair, 2003; Fukuda et al., 2013) and vasoconstriction causes an increase in vasomotion frequency (Colantuoni et al., 1984; Fujii et al., 1990). In contrast to MED, ISO causes vasodilation (Flynn et al., 1992; Farber et al., 1997; Iida et al., 1998). We showed that both CO₂ and ISO administration led to dilation of blood vessels and interrupted HDO. This is consistent with other studies showing that vasodilation decreases frequency of vasomotion, and that pronounced vasodilation abolishes vasomotion (Colantuoni et al., 1984; Fujii et al., 1990; Oude Vrielink et al., 1990; Obrig et al., 2000). We accomplished the CO₂ challenge by replacing a portion of the oxygen in the breathing gas by CO₂. Although the difference in oxygen supply was low (5%), we cannot explicitly assign the observed vasodilation to the resulting hypercapnia alone, a slight hypoxia may have contributed to the reduced vascular tone as well. Yet, vasodilation caused by hypercapnia can suspend vasomotion (Hudetz et al., 1998). This is relevant when using MED sedation as MED can induce hypercapnia (Brynildsen et al., 2017). However, vasomotion was not inhibited under MED sedation applied in the current study. In contrast to our results, a study in cats detected vasomotion during long-term ISO anesthesia (Rivadulla et al., 2011). However, a lower dose of ISO was used in that study, suggesting a dose-dependent effect of ISO. Accordingly, vasomotion is expected to occur not only during measurements with pure MED sedation, but also with a combination of MED and low ISO doses, a popular anesthetic regimen which provides temporally stable brain states in mice (Pradier et al., 2021), resembles the awake condition in rat brains more closely than MED or ISO alone (Paasonen et al., 2018) and enhances functional connectivity specificity as shown by Grandjean et al. (2022). Another important point is that ISO depresses sympathetic nervous activity (Skovsted and Sapthavichaikul, 1977) and that the latter was found to be a prerequisite for vasomotion (Colantuoni et al., 2001; Nilsson and Aalkjaer, 2003). Accordingly, vasomotion can only restart after end of ISO administration when both sympathetic nervous system activity has returned and vasodilation has regressed. In our experiments, the vasomotion suppressing effect of ISO had no long-term effect on consecutive measurements under MED. In contrast to Magnuson et al. (2014), we detected vasomotion already 1 h after switch from ISO to MED. Earlier times were not investigated as there was no stable physiological baseline directly after switching from ISO to MED. The early appearance of vasomotion implies that limiting experimental duration alone is unlikely to prevent the occurrence of oscillations.

4.3. Vasomotion significantly impacts analysis of fMRI data

Functional magnetic resonance imaging data are interpreted under the assumption that detected signal changes reflect neuronal activity due to neurovascular coupling. However, MR signal fluctuations of non-neuronal origin can distort data analysis. No direct correlation between our LFP and fMRI data was to be expected as the measurements were performed separately and LFP and hemodynamics show time-dependent variations. Nevertheless, the absence of a prominent signal component of the LFP signals at 0.2 Hz strongly suggests that the vasomotion-induced HDO have no direct neuronal correlate. This is consistent with other studies (Rivadulla et al., 2011; Winder et al., 2017). Mateo et al. (2017) reported a correlation between an envelope of γ-band (30–80 Hz) and vasomotion. In contrast to our study, these authors performed simultaneous hemodynamic and LFP recordings in awake mice. They interpreted their findings as evidence for a potential linkage between changes in brain oxygenation by modulation of vasomotion and ultra-slow variability in neuronal activity. Mechanistic insight supporting this notion is not available yet. Other studies found that vasomotion can also occur independently of neuronal activity. Winder et al. (2017) detected vasomotion after pharmacological blocking of neuronal activity, and further, spontaneous vasomotion was observed in isolated, pressurized vessels (Osol and Halpern, 1988). All together we concluded that vasomotion can introduce physiological noise in fMRI measurements and must be considered when evaluating these measurements.

While correction procedures for other physiological noise sources such as for example breathing and heart rate are available (Murphy et al., 2013), physiological noise from vasomotion is mostly neglected. Models used for GLM analysis usually ignore vasomotion and network analysis is based on the evaluation of temporal signal fluctuations irrespective of their origin. To investigate the influence of vasomotion on the analysis of fMRI data, we performed fMRI data analysis before and after HDO correction with a new suppression algorithm. Both GLM-based analysis and NBS network analysis showed that application of this algorithm had an impact on all data (measurements with and without vasomotion). However, in both analysis methods, HDO suppression had a significantly higher impact on data with vasomotion compared to vasomotion-free data. The relatively smaller effect on analysis of data without vasomotion indicated that the algorithm removed random noise. The much larger influence on data with HDO showed that vasomotion strongly affected GLM-based and NBS analysis of fMRI data.

4.4. Hemodynamic oscillations suppression complements preprocessing

One could argue that lowpass filtering strategies have already been implemented to minimize the impact of physiological noise in fMRI data. Complex preprocessing pipelines, such as the independent component analysis-based method FIX (Salimi-Khorshidi et al., 2014; Zerbi et al., 2015), presumably diminish the effects of vasomotion. For brain network analysis often a 0.1-Hz-lowpass filter is used for preprocessing. Indeed, its application removes HDO as effectively as the developed suppression algorithm (data not shown). However, Pan et al. (2013) showed a high correlation between LFP and fMRI ranging from 0.1 to 0.2 Hz in DMED sedated rats. Following this line of evidence, a study by Grandjean et al. (2014) set the filter to 0.3 Hz for MED/DMED anesthesia. For anesthetic regimens other than MED, we expect a lower oscillation frequency (see Section “4.2. Anesthesia exerts diverse effects on vasomotion”). Therefore, a 0.1 Hz lowpass filter may not be sufficient to remove HDO and vasomotion correction may still be needed. Accordingly, we recommend application of lowpass filter together with an HDO suppression algorithm.

4.5. Hemodynamic oscillations occur in both small and large blood vessels

The functional GE-EPI signal is most sensitive to large venous blood vessels, while the SE-EPI signal is most sensitive to small vessels (Buxton, 2013; Baez-Yanez et al., 2017). Since we detected oscillations in both GE-EPI and SE-EPI measurements, our results suggest that HDO occur in both small and large venous blood vessels. Since all regions of the brain are densely vascularized, oscillations are likely to occur everywhere in the brain. Therefore, it appears unlikely that it will be possible to avoid recordings of oscillations by excluding individual brain regions from the measurements. This is not only relevant for fMRI, but for all techniques that are sensitive to hemodynamic changes. In addition to optical recordings of blood volume or flow and laser doppler flowmetry, this also applies to optical recordings of fluorescent dyes, as hemodynamic changes cause artifacts in these measurements (Kozberg et al., 2016; Ma et al., 2016; Lambers et al., 2022; Zhang et al., 2022).

5. Conclusion

Here, we have investigated prevalence of HDO and their impact on analysis of fMRI data under anesthetic regimens that are widely used for rodent fMRI experiments. HDO were not detected under vasodilatory conditions such as for example ISO anesthesia. Under the commonly applied MED anesthesia, vasomotion-induced HDO occurred already during short experimental duration and in most regions of the brain. Since HDO significantly influenced analysis of fMRI data, HDO should be corrected during fMRI data analysis. We provide HDO detection and suppression algorithms that can be readily applied to fMRI data.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The animal study was reviewed and approved by the Landesamt für Natur, Umwelt und Verbraucherschutz of Nordrhein-Westfalen, Germany.

Author contributions

HL and CF designed the study and wrote the manuscript. HL and LW conducted the experiments. HL and CL performed the analysis. CL developed the HDO algorithms. All authors edited and approved the manuscript.

Funding

This work was funded by the German Research Foundation DFG (project 406818964) and the Interdisciplinary Center for Clinical Research Münster (PIX).

Acknowledgments

We thank Bruno Pradier for valuable comments and Florian Breuer for performing brain masking.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2023.1064000/full#supplementary-material

Abbreviations

EGF, exponential-Gaussian function; HDO, hemodynamic oscillations; MDCC, mean differences of correlation coefficients.

Footnotes

1. ^ https://www.fil.ion.ucl.ac.uk/spm/
2. ^ https://www.nitrc.org/frs/?group_id=889
3. ^ https://cytoscape.org/
4. ^ https://github.com/TheFaberLab/HDO

References

Aalkjaer, C., Boedtkjer, D., and Matchkov, V. (2011). Vasomotion - what is currently thought? Acta Physiol. 202, 253–269. doi: 10.1111/j.1748-1716.2011.02320.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Albers, F., Wachsmuth, L., Schache, D., Lambers, H., and Faber, C. (2019). Functional MRI readouts from BOLD and diffusion measurements differentially respond to optogenetic activation and tissue heating. Front. Neurosci. 13:1104. doi: 10.3389/fnins.2019.01104

PubMed Abstract | CrossRef Full Text | Google Scholar

Albers, F., Wachsmuth, L., van Alst, T. M., and Faber, C. (2018). Multimodal functional neuroimaging by simultaneous BOLD fMRI and fiber-optic calcium recordings and optogenetic control. Mol. Imaging Biol. 20, 171–182. doi: 10.1007/s11307-017-1130-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Aleksandrin, V. V., Ivanov, A. V., Virus, E. D., Bulgakova, P. O., and Kubatiev, A. A. (2018). Application of wavelet analysis to detect dysfunction in cerebral blood flow autoregulation during experimental hyperhomocysteinaemia. Lasers Med. Sci. 33, 1327–1333. doi: 10.1007/s10103-018-2485-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Amirmohseni, S., Segelcke, D., Reichl, S., Wachsmuth, L., Gorlich, D., Faber, C., et al. (2016). Characterization of incisional and inflammatory pain in rats using functional tools of MRI. Neuroimage 127, 110–122. doi: 10.1016/j.neuroimage.2015.11.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Baez-Yanez, M. G., Ehses, P., Mirkes, C., Tsai, P. S., Kleinfeld, D., and Scheffler, K. (2017). The impact of vessel size, orientation and intravascular contribution on the neurovascular fingerprint of BOLD bSSFP fMRI. Neuroimage 163, 13–23. doi: 10.1016/j.neuroimage.2017.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Brynildsen, J. K., Hsu, L. M., Ross, T. J., Stein, E. A., Yang, Y., and Lu, H. (2017). Physiological characterization of a robust survival rodent fMRI method. Magn. Reson. Imaging 35, 54–60. doi: 10.1016/j.mri.2016.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Buxton, R. B. (2013). The physics of functional magnetic resonance imaging (fMRI). Rep. Prog. Phys. 76:096601. doi: 10.1088/0034-4885/76/9/096601

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Pais-Roldan, P., Chen, X., Frosz, M. H., and Yu, X. (2019). MRI-guided robotic arm drives optogenetic fMRI with concurrent Ca(2+) recording. Nat. Commun. 10:2536. doi: 10.1038/s41467-019-10450-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Colantuoni, A., Bertuglia, S., and Intaglietta, M. (1984). The effects of alpha- or beta-adrenergic receptor agonists and antagonists and calcium entry blockers on the spontaneous vasomotion. Microvasc. Res. 28, 143–158. doi: 10.1016/0026-2862(84)90014-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Colantuoni, A., Bertuglia, S., and Marchiafava, P. L. (2001). Phentolamine suppresses the increase in arteriolar vasomotion frequency due to systemic hypoxia in hamster skeletal muscle microcirculation. Autonomic Neurosci. 90, 148–151. doi: 10.1016/S1566-0702(01)00281-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Drew, P. J., Shih, A. Y., Driscoll, J. D., Knutsen, P. M., Blinder, P., Davalos, D., et al. (2010). Chronic optical access through a polished and reinforced thinned skull. Nat. Methods 7, 981–984. doi: 10.1038/nmeth.1530

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, J. L., Rivera, J. A., Sun, W., Peterson, J., Haeberle, H., Rubin, S., et al. (2020). High-speed volumetric two-photon fluorescence imaging of neurovascular dynamics. Nat. Commun. 11:6020. doi: 10.1038/s41467-020-19851-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Farber, N. E., Harkin, C. P., Niedfeldt, J., Hudetz, A. G., Kampine, J. P., and Schmeling, W. T. (1997). Region-specific and agent-specific dilation of intracerebral microvessels by volatile anesthetics in rat brain slices. Anesthesiology 5, 1191–1198. doi: 10.1097/00000542-199711000-00024

PubMed Abstract | CrossRef Full Text | Google Scholar

Flynn, N. M., Buljubasic, N., Bosnjak, Z. J., and Kampine, J. P. (1992). Isoflurane produces endothelium-independent relaxation in canine middle cerebral arteries. Anesthesiology 76, 334–41. doi: 10.1097/00000542-199203000-00021

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J., Holmes, A. P., Worsley, K. J., Poline, J.-P., Frith, C. D., and Frackowiak, R. S. J. (1995). Statistical parametric maps in functional imaging: A general linear approach. Hum. Brain Mapp. 2, 189–210. doi: 10.1002/hbm.460020402

CrossRef Full Text | Google Scholar

Fujii, K., Heistad, D. D., and Faraci, F. M. (1990). Vasomotion of basilar arteries in vivo. Am. J. Physiol. 258, H1829–H1834. doi: 10.1152/ajpheart.1990.258.6.H1829

PubMed Abstract | CrossRef Full Text | Google Scholar

Fukuda, M., Vazquez, A. L., Zong, X., and Kim, S. G. (2013). Effects of the alpha(2)-adrenergic receptor agonist dexmedetomidine on neural, vascular and BOLD fMRI responses in the somatosensory cortex. Eur. J. Neurosci. 37, 80–95. doi: 10.1111/ejn.12024

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y. R., Ma, Y., Zhang, Q., Winder, A. T., Liang, Z., Antinori, L., et al. (2017). Time to wake up: Studying neurovascular coupling and brain-wide circuit function in the un-anesthetized animal. Neuroimage 153, 382–398. doi: 10.1016/j.neuroimage.2016.11.069

PubMed Abstract | CrossRef Full Text | Google Scholar

Grandjean, J., Desrosiers-Gregoire, G., Anckaerts, C., Angeles-Valdez, D., Ayad, F., Barrière, D. A., et al. (2022). StandardRat: A multi-center consensus protocol to enhance functional connectivity specificity in the rat brain. bioRxiv [preprint] doi: 10.1101/2022.04.27.489658

CrossRef Full Text | Google Scholar

Grandjean, J., Schroeter, A., Batata, I., and Rudin, M. (2014). Optimization of anesthesia protocol for resting-state fMRI in mice based on differential effects of anesthetics on functional connectivity patterns. Neuroimage 102, 838–847. doi: 10.1016/j.neuroimage.2014.08.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Haensel, J. X., Spain, A., and Martin, C. (2015). A systematic review of physiological methods in rodent pharmacological MRI studies. Psychopharmacology 232, 489–499. doi: 10.1007/s00213-014-3855-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Hudetz, A. G., Biswal, B. B., Shen, H., Lauer, K. K., and Kampine, J. P. (1998). Spontaneous fluctuations in cerebral oxygen supply. Adv. Exp. Med. Biol. 454, 551–559. doi: 10.1007/978-1-4615-4863-8_66

PubMed Abstract | CrossRef Full Text | Google Scholar

Hundley, W. G., Renaldo, G. J., Levasseur, J. E., and Kontos, H. A. (1988). Vasomotion in cerebral microcirculation of awake rabbits. Am. J. Physiol. 254, H67–H71. doi: 10.1152/ajpheart.1988.254.1.H67

PubMed Abstract | CrossRef Full Text | Google Scholar

Iida, H., Ohata, H., Iida, M., Watanabe, Y., and Dohi, S. (1998). Isoflurane and sevoflurane induce vasodilation of cerebral vessels via ATP-sensitive K+ channel activation. Anesthesiology 89, 954–960. doi: 10.1097/00000542-199810000-00020

PubMed Abstract | CrossRef Full Text | Google Scholar

Ioanas, H. I., Schlegel, F., Skachokova, Z., Schroeter, A., Husak, T., and Rudin, M. (2022). Hybrid fiber optic-fMRI for multimodal cell-specific recording and manipulation of neural activity in rodents. Neurophotonics 9:032206. doi: 10.1117/1.NPh.9.3.032206

PubMed Abstract | CrossRef Full Text | Google Scholar

Kint, L. T., Seewoo, B. J., Hyndman, T. H., Clarke, M. W., Edwards, S. H., Rodger, J., et al. (2020). The Pharmacokinetics of Medetomidine Administered Subcutaneously during Isoflurane Anaesthesia in Sprague-Dawley Rats. Animals 10:1050. doi: 10.3390/ani10061050

PubMed Abstract | CrossRef Full Text | Google Scholar

Kleinfeld, D., Mitra, P. P., Helmchen, F., and Denk, W. (1998). Fluctuations and stimulus-induced changes in blood flow observed in individual capillaries in layers 2 through 4 of rat neocortex. Proc. Natl. Acad. Sci. U S A. 95, 15741–15746. doi: 10.1073/pnas.95.26.15741

PubMed Abstract | CrossRef Full Text | Google Scholar

Kosten, L., Emmi, S. A., Missault, S., and Keliris, G. A. (2022). Combining magnetic resonance imaging with readout and/or perturbation of neural activity in animal models: Advantages and pitfalls. Front. Neurosci. 16:938665. doi: 10.3389/fnins.2022.938665

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozberg, M. G., Ma, Y., Shaik, M. A., Kim, S. H., and Hillman, E. M. (2016). Rapid postnatal expansion of neural networks occurs in an environment of altered neurovascular and neurometabolic coupling. J. Neurosci. 36, 6704–6717. doi: 10.1523/JNEUROSCI.2363-15.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Lake, E. M. R., Ge, X., Shen, X., Herman, P., Hyder, F., Cardin, J. A., et al. (2020). Simultaneous cortex-wide fluorescence Ca(2+) imaging and whole-brain fMRI. Nat. Methods 17, 1262–1271. doi: 10.1038/s41592-020-00984-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambers, H., Segeroth, M., Albers, F., Wachsmuth, L., van Alst, T. M., and Faber, C. (2020). A cortical rat hemodynamic response function for improved detection of BOLD activation under common experimental conditions. Neuroimage 208:116446. doi: 10.1016/j.neuroimage.2019.116446

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambers, H., Wachsmuth, L., Thomas, D., Boumezbeur, F., Hoesker, V., Pradier, B., et al. (2022). Fiber-based lactate recordings with fluorescence resonance energy transfer sensors by applying an magnetic resonance-informed correction of hemodynamic artifacts. Neurophotonics 9:032212. doi: 10.1117/1.NPh.9.3.032212

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Tam, E. W., Kwan, M. P., Mak, A. F., Lo, S. C., and Leung, M. C. (2006). Effects of prolonged surface pressure on the skin blood flowmotions in anaesthetized rats–an assessment by spectral analysis of laser doppler flowmetry signals. Phys. Med. Biol. 51, 2681–2694. doi: 10.1088/0031-9155/51/10/020

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Zhao, L., and Liu, Y. (2022). The origin of vasomotion and stochastic resonance in vasomotion. Front. Bioeng. Biotechnol. 10:819716. doi: 10.3389/fbioe.2022.819716

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Y., Shaik, M. A., Kim, S. H., Kozberg, M. G., Thibodeaux, D. N., Zhao, H. T., et al. (2016). Wide-field optical mapping of neural activity and brain haemodynamics: considerations and novel approaches. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 371:1705. doi: 10.1098/rstb.2015.0360

PubMed Abstract | CrossRef Full Text | Google Scholar

Magnuson, M., Majeed, W., and Keilholz, S. D. (2010). Functional connectivity in blood oxygenation level-dependent and cerebral blood volume-weighted resting state functional magnetic resonance imaging in the rat brain. J. Magn. Reson. Imaging 32, 584–592. doi: 10.1002/jmri.22295

PubMed Abstract | CrossRef Full Text | Google Scholar

Magnuson, M. E., Thompson, G. J., Pan, W. J., and Keilholz, S. D. (2014). Time-dependent effects of isoflurane and dexmedetomidine on functional connectivity, spectral characteristics, and spatial distribution of spontaneous BOLD fluctuations. NMR Biomed. 27, 291–303. doi: 10.1002/nbm.3062

PubMed Abstract | CrossRef Full Text | Google Scholar

Mateo, C., Knutsen, P. M., Tsai, P. S., Shih, A. Y., and Kleinfeld, D. (2017). Entrainment of Arteriole Vasomotor Fluctuations by Neural Activity Is a Basis of Blood-Oxygenation-Level-Dependent “Resting-State” Connectivity. Neuron 96, 936–948 e933. doi: 10.1016/j.neuron.2017.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayhew, J. E., Askew, S., Zheng, Y., Porrill, J., Westby, G. W., Redgrave, P., et al. (1996). Cerebral vasomotion: a 0.1-Hz oscillation in reflected light imaging of neural activity. 4(3 Pt 1), 186–193. doi: 10.1006/nimg.1996.0069

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitra, P. P., Ogawa, S., Hu, X., and Uğurbil, K. (1997). The nature of spatiotemporal changes in cerebral hemodynamics as manifested in functional magnetic resonance imaging. Magn. Reson. Med. 37, 511–518. doi: 10.1002/mrm.1910370407

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy, K., Birn, R. M., and Bandettini, P. A. (2013). Resting-state fMRI confounds and cleanup. Neuroimage 80, 349–359. doi: 10.1016/j.neuroimage.2013.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Nilsson, H., and Aalkjaer, C. (2003). Vasomotion: mechanisms and physiological importance. Mole. Intervent. 3, 79–51. doi: 10.1124/mi.3.2.79

PubMed Abstract | CrossRef Full Text | Google Scholar

Obrig, H., Neufang, M., Wenzel, R., Kohl, M., Steinbrink, J., Einhaupl, K., et al. (2000). Spontaneous low frequency oscillations of cerebral hemodynamics and metabolism in human adults. Neuroimage 12, 623–639. doi: 10.1006/nimg.2000.0657

PubMed Abstract | CrossRef Full Text | Google Scholar

Ogawa, S., Lee, T. M., Nayak, A. S., and Glynn, P. (1990). Oxygenation-sensitive contrast in magnetic resonance image of rodent brain at high magnetic fields. Magn. Reson. Med. 14, 68–78. doi: 10.1002/mrm.1910140108

PubMed Abstract | CrossRef Full Text | Google Scholar

Osol, G., and Halpern, W. (1988). Spontaneous vasomotion in pressurized cerebral arteries from genetically hypertensive rats. Am. J. Physiol. 254, H28–H33. doi: 10.1152/ajpheart.1988.254.1.H28

PubMed Abstract | CrossRef Full Text | Google Scholar

Oude Vrielink, H. H., Slaaf, D. W., Tangelder, G. J., Weijmer-Van Velzen, S., and Reneman, R. S. (1990). Analysis of vasomotion waveform changes during pressure reduction and adenosine application. Am. J. Physiol. 258, H29–H37. doi: 10.1152/ajpheart.1990.258.1.H29

PubMed Abstract | CrossRef Full Text | Google Scholar

Paasonen, J., Stenroos, P., Salo, R. A., Kiviniemi, V., and Grohn, O. (2018). Functional connectivity under six anesthesia protocols and the awake condition in rat brain. Neuroimage 172, 9–20. doi: 10.1016/j.neuroimage.2018.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, W. J., Thompson, G., Magnuson, M., Majeed, W., Jaeger, D., and Keilholz, S. (2010). Simultaneous FMRI and electrophysiology in the rodent brain. J. Vis. Exp. 42:1902. doi: 10.3791/1901

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, W. J., Thompson, G. J., Magnuson, M. E., Jaeger, D., and Keilholz, S. (2013). Infraslow LFP correlates to resting-state fMRI BOLD signals. Neuroimage 74, 288–297. doi: 10.1016/j.neuroimage.2013.02.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Paxinos, G. W. C. (1997). The rat brain in stereotactic coordinates. United States, MA: Academic Press.

Google Scholar

Pradier, B., Wachsmuth, L., Nagelmann, N., Segelcke, D., Kreitz, S., Hess, A., et al. (2021). Combined resting state-fMRI and calcium recordings show stable brain states for task-induced fMRI in mice under combined ISO/MED anesthesia. Neuroimage 245:118626. doi: 10.1016/j.neuroimage.2021.118626

PubMed Abstract | CrossRef Full Text | Google Scholar

Rayshubskiy, A., Wojtasiewicz, T. J., Mikell, C. B., Bouchard, M. B., Timerman, D., Youngerman, B. E., et al. (2014). Direct, intraoperative observation of ^~0.1 Hz hemodynamic oscillations in awake human cortex: implications for fMRI. Neuroimage 87, 323–331. doi: 10.1016/j.neuroimage.2013.10.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Rivadulla, C., de Labra, C., Grieve, K. L., and Cudeiro, J. (2011). Vasomotion and neurovascular coupling in the visual thalamus in vivo. PLoS One 6:e28746. doi: 10.1371/journal.pone.0028746

PubMed Abstract | CrossRef Full Text | Google Scholar

Salimi-Khorshidi, G., Douaud, G., Beckmann, C. F., Glasser, M. F., Griffanti, L., and Smith, S. M. (2014). Automatic denoising of functional MRI data: combining independent component analysis and hierarchical fusion of classifiers. Neuroimage 90, 449–468. doi: 10.1016/j.neuroimage.2013.11.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulz, K., Sydekum, E., Krueppel, R., Engelbrecht, C. J., Schlegel, F., Schroter, A., et al. (2012). Simultaneous BOLD fMRI and fiber-optic calcium recording in rat neocortex. Nat. Methods 9, 597–602. doi: 10.1038/nmeth.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Sinclair, M. D. (2003). A review of the physiological effects of alpha2-agonists related to the clinical use of medetomidine in small animal practice. Canadian Vet. J. 44, 885–897.

PubMed Abstract | Google Scholar

Sirmpilatze, N., Baudewig, J., and Boretius, S. (2019). Temporal stability of fMRI in medetomidine-anesthetized rats. Sci. Rep. 9:16673. doi: 10.1038/s41598-019-53144-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Sirotin, Y. B., Cardoso, M., Lima, B., and Das, A. (2012). Spatial homogeneity and task-synchrony of the trial-related hemodynamic signal. Neuroimage 59, 2783–2797. doi: 10.1016/j.neuroimage.2011.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Sirotin, Y. B., and Das, A. (2009). Anticipatory haemodynamic signals in sensory cortex not predicted by local neuronal activity. Nature 457, 475–479. doi: 10.1038/nature07664

PubMed Abstract | CrossRef Full Text | Google Scholar

Skovsted, P., and Sapthavichaikul, S. (1977). The effects of isoflurane on arterial pressure, pulse rate, autonomic nervous activity, and barostatic reflexes. Can. Anaesth. Soc. J. 24, 304–314. doi: 10.1007/BF03005103

PubMed Abstract | CrossRef Full Text | Google Scholar

Sommers, M. G., van Egmond, J., Booij, L. H., and Heerschap, A. (2009). Isoflurane anesthesia is a valuable alternative for alpha-chloralose anesthesia in the forepaw stimulation model in rats. NMR Biomed. 22, 414–418. doi: 10.1002/nbm.1351

PubMed Abstract | CrossRef Full Text | Google Scholar

Tenney, J. R., Duong, T. Q., King, J. A., Ludwig, R., and Ferris, C. F. (2003). Corticothalamic modulation during absence seizures in rats: a functional MRI assessment. Epilepsia 9, 1133–1140. doi: 10.1046/j.1528-1157.2003.61002.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tristao Pereira, C., Diao, Y., Yin, T., da Silva, A. R., Lanz, B., Pierzchala, K., et al. (2021). Synchronous nonmonotonic changes in functional connectivity and white matter integrity in a rat model of sporadic Alzheimer’s disease. Neuroimage 225:117498. doi: 10.1016/j.neuroimage.2020.117498

PubMed Abstract | CrossRef Full Text | Google Scholar

van Alst, T. M., Wachsmuth, L., Datunashvili, M., Albers, F., Just, N., Budde, T., et al. (2019). Anesthesia differentially modulates neuronal and vascular contributions to the BOLD signal. Neuroimage 195, 89–103. doi: 10.1016/j.neuroimage.2019.03.057

PubMed Abstract | CrossRef Full Text | Google Scholar

Wachsmuth, L., Datunashvili, M., Kemper, K., Albers, F., Lambers, H., Luttjohann, A., et al. (2021). Retrosplenial Cortex Contributes to Network Changes during Seizures in the GAERS Absence Epilepsy Rat Model. Cereb Cortex Commun. 2:tgab023. doi: 10.1093/texcom/tgab023

PubMed Abstract | CrossRef Full Text | Google Scholar

Weber, R., Ramos-Cabrer, P., Wiedermann, D., van Camp, N., and Hoehn, M. (2006). A fully noninvasive and robust experimental protocol for longitudinal fMRI studies in the rat. Neuroimage 29, 1303–1310. doi: 10.1016/j.neuroimage.2005.08.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Winder, A. T., Echagarruga, C., Zhang, Q., and Drew, P. J. (2017). Weak correlations between hemodynamic signals and ongoing neural activity during the resting state. Nat. Neurosci. 20, 1761–1769. doi: 10.1038/s41593-017-0007-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zalesky, A., Fornito, A., and Bullmore, E. T. (2010). Network-based statistic: identifying differences in brain networks. Neuroimage 53, 1197–1207. doi: 10.1016/j.neuroimage.2010.06.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Zerbi, V., Grandjean, J., Rudin, M., and Wenderoth, N. (2015). Mapping the mouse brain with rs-fMRI: An optimized pipeline for functional network identification. Neuroimage 123, 11–21. doi: 10.1016/j.neuroimage.2015.07.090

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W. T., Chao, T. H., Yang, Y., Wang, T. W., Lee, S. H., Oyarzabal, E. A., et al. (2022). Spectral fiber photometry derives hemoglobin concentration changes for accurate measurement of fluorescent sensor activity. Cell Rep. Methods 2:100243. doi: 10.1016/j.crmeth.2022.100243

PubMed Abstract | CrossRef Full Text | Google Scholar