Abstract
Nasal breathing generates a rhythmic signal which entrains cortical network oscillations in widespread brain regions on a cycle-to-cycle time scale. It is unknown, however, how respiration and neuronal network activity interact on a larger time scale: are breathing frequency and typical neuronal oscillation patterns correlated? Is there any directionality or temporal relationship? To address these questions, we recorded field potentials from the posterior parietal cortex of mice together with respiration during REM sleep. In this state, the parietal cortex exhibits prominent θ and γ oscillations while behavioral activity is minimal, reducing confounding signals. We found that the instantaneous breathing frequency strongly correlates with the instantaneous frequency and amplitude of both θ and γ oscillations. Cross-correlograms and Granger causality revealed specific directionalities for different rhythms: changes in θ activity precede and Granger-cause changes in breathing frequency, suggesting control by the functional state of the brain. On the other hand, the instantaneous breathing frequency Granger causes changes in γ frequency, suggesting that γ is influenced by a peripheral reafference signal. These findings show that changes in breathing frequency temporally relate to changes in different patterns of rhythmic brain activity. We hypothesize that such temporal relations are mediated by a common central drive likely to be located in the brainstem.
Significance Statement
The study of the interactions between respiration and brain activity has been focused on phase entrainment relations, in which cortical networks oscillate phase-locked to breathing cycles. Here we discovered new and much broader interactions which link breathing frequency to different patterns of oscillatory brain activity. Specifically, we show that the instantaneous breathing frequency strongly correlates with the instantaneous frequency and amplitude of θ and γ oscillations, two major network patterns associated with cognitive functions. Interestingly, changes in breathing frequency follow θ, suggesting a central drive, while in contrast, γ activity follows changes in breathing frequency, suggesting the role of a reafferent signal. Our results reveal new mechanisms by which nasal breathing patterns may relate to brain functions.
Introduction
Several recent studies have shown that nasal breathing induces local field potential (LFP) oscillations at the same frequency as respiration in multiple regions of the rodent and human brain (for reviews, see Tort et al., 2018a; Heck et al., 2019). LFP phase-locking to respiratory cycles includes regions not primarily related to olfaction (Ito et al., 2014; Lockmann et al., 2016; Nguyen Chi et al., 2016; Zelano et al., 2016; Biskamp et al., 2017; Herrero et al., 2018; Karalis and Sirota, 2018; Köszeghy et al., 2018; Rojas-Líbano et al., 2018; Tort et al., 2018b), suggesting that respiration-entrained oscillations aid the integration of widespread information (Heck et al., 2017, 2019; Tort et al., 2018a), similar to the proposed function of other slow network rhythms (Isomura et al., 2006; Canolty and Knight, 2010). Given the alleged beneficial effects of some respiratory practices to mood, behavioral performance, and cognitive abilities (Pascoe and Bauer, 2015; Zelano et al., 2016; Melnychuk et al., 2018; Nakamura et al., 2018; Perl et al., 2019; Novaes et al., 2020), the fact that nasal breathing can modulate brain oscillations in nonolfactory regions has sparked large interest, including in popular science media, since it could provide a way through which respiration would affect brain functions (Heck et al., 2017; Varga and Heck, 2017; Tort et al., 2018a).
The research on relations between respiration and neuronal network activity has so far mainly focused on phase entrainment relations, in which a given cortical LFP is shown to oscillate phase-locked to breathing cycles, typically accompanied by a peak in the LFP power spectrum at the breathing frequency (Tort et al., 2018b), or by modulation of action potential probability by the breathing phase (Zhong et al., 2017; Karalis and Sirota, 2018; Köszeghy et al., 2018). Other studies have shown that the breathing phase is also capable of modulating the instantaneous amplitude of faster, γ-frequency oscillations (Ito et al., 2014; Yanovsky et al., 2014; Rojas-Líbano et al., 2018; Cavelli et al., 2019), especially in frontal regions (Biskamp et al., 2017; Zhong et al., 2017). All these findings are, however, confined to the time domain of individual oscillation cycles. Nevertheless, both breathing frequency and the pattern of neuronal network oscillations vary over time, and it remains unclear whether there is any interaction between these rhythmic phenomena at the time scale of such variations. Of note, we have recently reported such an interaction by showing that the strength of θ-γ coupling depends on respiratory frequency (Hammer et al., 2020).
We therefore performed a systematic analysis of relations between the varying patterns of cortical oscillations and breathing. In particular, we searched for correlations and Granger causality between the instantaneous breathing frequency and the instantaneous amplitude and frequency of neocortical network oscillations. To that end, we recorded respiration through a whole-body plethysmograph while simultaneously recording field potentials from the posterior parietal cortex of mice. We focused on periods of REM sleep, which is characterized by prominent activity in the θ and γ bands (Montgomery et al., 2008; Scheffzük et al., 2011). Moreover, this state guarantees a stable behavior: that is, any observed changes in breathing are not because of changes in motor demands, and cognitive-behavioral interference with the environment is minimal. Our results reveal strong correlations between the instantaneous breathing frequency and the instantaneous amplitude and frequency of θ and γ oscillations. We also report a clear temporal relation, with θ oscillations being Granger-causal for changes in breathing frequency which, in turn, is Granger-causal for γ activity.
Materials and Methods
Ethics statement
This study was approved by the Governmental Supervisory Panel on Animal Experiments of Baden-Württemberg (35-9185.81 G84/13 and G-115/14). All experiments were conducted in line with the guidelines of the European Science Foundation (2001) and the U.S. National Institutes of Health’s Guide for the care and use of laboratory animals (1996).
Animal care
We used 22 C57/Bl6 mice (9 females and 13 males; age: 13-19 weeks; weight: 22-37 g). The animals were housed in groups of 4 with access to food and water ad libitum. After surgery for electrode implantation, animals were kept individually until completion of the experiments. Mice were housed on a 12 h light-dark-cycle (lights off at 8:00 A.M., except 2 animals on the opposite circadian phase). Data from these same animals were used in a recent study of ours (Hammer et al., 2020).
Surgery
Mice were implanted with epidural surface electrodes. Buprenorphine (0.1 mg/kg) was administered for pain control before and after surgery (every 8 h as necessary). Surgery was performed under isoflurane anesthesia (4% for induction, 1.5% during surgery) (for further details, see Jessberger et al., 2016; X. Zhang et al., 2016). After skull exposure, 0.5- to 1-mm-diameter holes were drilled above the parietal cortex according to stereotactic coordinates (2 mm posterior bregma, 1.5 mm lateral to the midline). Reference and ground electrodes were screwed into the skull above the cerebellum. Electrodes consisted of stainless-steel watch screws. During surgery, body temperature was monitored and maintained at 37°C-38°C. Animals were allowed 7 d of recovery before experiments.
Electrophysiological recordings and behavioral staging
Mice were placed in a whole-body plethysmograph (EMKA Technologies), which was customized to allow simultaneous recordings of respiration and LFPs. Of note, in this work, we use the terms “respiration” and “breathing” interchangeably. The plethysmograph consisted of a transparent cylinder (78 mm inner diameter, 165 mm height) connected to a reference chamber (see Fig. 1A). Individual recording sessions lasted from 2 to 6 h (mean: 3 h). To avoid slow signal drifts, the plethysmograph-derived respiration signal was high-pass filtered >1 Hz. Monopolar electrophysiological signals were filtered (1-500 Hz), amplified and digitized at 2.5 kHz (RHA2000, Intan Technologies). A three-dimensional accelerometer was custom-mounted on the amplifier board located on the head of the mice to allow movement detection. The three signals of the accelerometer were fed to three channels of the amplifier using the same bandpass filter as for the LFPs (see above), therefore removing the gravity-induced sustained potentials of the accelerometers.
Respiration drives phase-entrained oscillations in the parietal cortex during REM sleep. A, Plethysmograph recording of respiration (Resp) along with the parietal LFP during REM sleep. B, Mean (±SEM) power spectra for the two signals (n = 22 mice; the individual spectra were normalized by the maximal value before computing the mean). In addition to prominent theta (θ, 6-12 Hz) oscillations, the LFP spectrum has a smaller power bump at a similar frequency range as respiration (2-5 Hz), which corresponds to the respiration-entrained LFP rhythm (RR) (see Tort et al., 2018b). C, Mean phase coherence spectrum between the respiration and LFP signals, showing marked coherence at the respiration (and not theta) frequency. D, Mean CCG between respiration and LFP. E, Mean Granger causality for the respiration → LFP (blue) and LFP → respiration (black) directions. Left, Frequency spectrum. Right, Overall (time domain) causality levels of each direction. Here and in subsequent figures, black squares represent the mean (±SEM) and connected circles represent individual animal data. Respiration highly Granger-causes LFP activity at the breathing frequency (i.e., respiration Granger-causes RR). The raw LFP and respiration signals were used in this analysis (as opposed to results shown in the subsequent figures). C, E, Surrogate spectra (mean + 2 SD) for coherence or Granger causality were obtained by pairing respiration and LFP signals from different animals (Tort et al., 2018b); these show that the coherence and causality peaks at the respiration frequency are not because of the power peaks present in the individual time series. ***p < 0.001.
REM sleep was manually staged by a senior researcher in the field (J.B.) and identified as follows: (1) minimal accelerometer activity and (2) continuous θ rhythm in the parietal cortex following a sleep stage with slow waves typical for non-REM sleep (Brankačk et al., 2010).
Data analysis
We used built-in and custom-written routines in MATLAB (The Mathworks). For each animal, all identified REM sleep epochs were first detrended and then concatenated before analysis.
Power and coherence
Respiration and LFP power spectra were computed using the pwelch.m MATLAB function from the Signal Processing Toolbox. Phase coherence spectra between raw respiration and LFP signals were computed using the mscohere.m function (Signal Processing Toolbox). For both functions, we used 4 s windows with 50% overlap and nfft = 216.
Time-frequency power spectrograms were computed using the spectrogram.m function (Signal Processing Toolbox). For the analyses shown in Figures 2, 4, and 5, we used 500 ms windows with no overlap. The numerical frequency resolution was set to 0.05 Hz for slower frequencies (0.5-20 Hz) and 0.2 Hz for faster frequencies (25-200 Hz). The power spectra shown in Figure 2C were obtained from the spectrogram in Figure 2A (mean over two 500 ms windows). For each 500 ms window, the corresponding peak frequency for θ or respiration was taken as the frequency with maximal power. For estimating slow- and fast-γ peak frequencies, the power spectrum was first smoothed (15 Hz moving average) and the fit of the power decay with increasing LFP frequencies (see below) was subsequently subtracted (see Fig. 4A) (Brankačk et al., 2012; Scheffzük et al., 2013). The absolute slow- and fast-γ power analyzed in Figure 7 were taken as the mean power in the corresponding frequency range; the “detrended” slow- and fast-γ power analyzed in Figure 8 was obtained as the mean power at the same frequency ranges but computed after subtracting the power decay fit. For the directionality analyses shown in Figures 3, 6, 7, 8 (see below), the time series of instantaneous frequency and power were estimated from spectrograms computed using 500 ms windows but with 90% overlap (i.e., 50 ms step or 20 Hz sampling).
Tracking respiration during REM sleep reveals a link between breathing frequency and theta activity. A, Example time-frequency power spectrogram for the respiration signal using nonoverlapping 500 ms windows. The power spectrum of each window was normalized by the maximal value. Black line indicates the breathing frequency, as defined by the peak power frequency. B, The 9 s of a respiration spectrogram (top) and the corresponding raw signal (bottom). C, Example of respiration power spectra (top) and raw signals for two 1 s epochs differing in breathing frequency. D, Histogram counts of theta and respiration peak frequency in 500 ms windows (pool over 22 animals). E, Top, Mean respiration power spectra for (overlapping) subsets of 500 ms windows binned by breathing frequency (bin intervals: 1-3, 2-4, 3-5, 4-6, 5-7, 6-8, 7-9, 8-10, 9-11, and >10 Hz). Results obtained for a representative animal. Bottom, Average LFP power spectra for the same respiration-binned time windows, as labeled by the color code. Both theta power and theta frequency increase with faster breathing. Insets, Theta peak frequency and power as a function of breathing frequency. F, Group data of theta frequency and normalized theta power as a function of respiration frequency (mean ± SEM; n = 22 mice). The normalization was performed within animals by dividing power values by the value in the first respiration frequency bin.
Estimating γ frequency and power. A, Left, Average LFP power spectrum in logarithm (dB) scale. Notice power bumps at the slow- and fast-γ sub-bands. Orange represents the power decay fit (see Materials and Methods). Inset, Example raw (black) and filtered (colored, as labeled) LFP traces. Middle, Power values of fast LFP frequencies (40-200 Hz) in linear scale (gray trace). Thick blue curve indicates smoothed power values. Right, Detrended power spectrum, obtained by subtracting the power decay fit from the smoothed power. The peak frequency of slow- and fast-γ was estimated from the detrended power curve. B, Left, Mean respiration power spectra for respiration-binned time windows as in Figure 2E (different animal). Middle, Mean LFP power spectra for the same respiration-binned time windows plotted in logarithmic scale and focused in the fast frequency bands. Inset, Power decay fits. The frequency and absolute power of both γ sub-bands increase with faster breathing. Right, Mean detrended power spectra. Inset, Power decay fits in linear scale.
Changes in theta frequency precede changes in breathing frequency. A, Mean CCG (±SEM) between the time series of the instantaneous frequency of theta and respiration during REM sleep (n = 22 mice). B, Frequency (right) and time domain (left) Granger causality between the instantaneous frequency of theta and respiration, showing that theta frequency changes Granger-cause breathing frequency changes during sleep. Inset, Mean power spectra of the analyzed time series. C, D, Same as above, but for theta power. For these and subsequent directionality analyses, the peak LFP frequency and power and the peak respiration frequency were estimated using overlapping 500 ms windows (see Materials and Methods). *p < 0.05. ***p < 0.001.
γ activity depends on breathing frequency. A, B, Group data of peak frequency (left), absolute (middle), and detrended power (right) as a function of respiration frequency (mean ± SEM; n = 22 mice) for slow (A) and fast (B) γ. Power values were normalized by the first respiration frequency bin.
Changes in slow γ frequency follow changes in breathing frequency. A, Mean CCG (±SEM) between the instantaneous frequency of slow (top) and fast (bottom) γ frequencies and respiration frequency during REM sleep (n = 22 mice). B, Frequency (right) and time domain (left) Granger causality between the instantaneous frequency of γ and respiration. Insets, Mean power spectra of the γ frequency time series. Respiration frequency Granger-causes slow γ frequency while having no directionality relation with fast γ frequency. ***p < 0.001, ns, not significant.
Broadband power changes Granger-cause breathing frequency. A, Mean CCG (±SEM) between the absolute power of slow (top) and fast (bottom) γ and respiration frequency during REM sleep. B, Frequency (right) and time domain (left) Granger causality between absolute γ power and respiration. Insets, Mean power spectra of the γ power time series. Absolute γ power changes Granger-cause respiration frequency. C, D, As above, but for the mean value of the power decay fit, which tracks broadband changes in power. Note high Granger-causal relations between the broadband power changes and respiration frequency. ***p < 0.001.
Detrended slow γ power follows breathing changes. A, Mean CCG (±SEM) between the detrended power of slow (top) and fast (bottom) γ and respiration frequency during REM sleep (n = 22 mice). B, Frequency (left) and time domain (right) Granger causality between detrended γ power and respiration. Insets, Mean power spectra of the detrended γ power time series. Respiration frequency Granger-causes detrended slow γ power and has no directionality relation with detrended fast γ power, thus similar to its effect over γ frequency (compare Fig. 6). *p < 0.05, ns, not significant.
Power decay fit
In order to fit the decay of LFP power with increasing LFP frequencies, we used smoothed power values in three sub-bands: 30-40, 80-90, and 190-200 Hz (i.e., a sub-band lower than slow-γ, a sub-band between slow- and fast-γ, and a sub-band higher than fast-γ). For fitting power values, we tested both 1/f (i.e., Power(f) = constant/fa) and exponential decay fits (i.e., Power(f) = constant*e-a*f), and found lower fitting deviations (mean absolute error) when using the latter. To perform the exponential decay fit, power values of the three sub-bands were transformed using the natural logarithm, and the resulting linear relation Y = aX + b (where Y = ln(Power(f)) and X = f) was fitted using the polyfit.m function in MATLAB, which uses a least squares algorithm. The fitting parameters were then fed into the polyval.m function to generate the linear fit for frequencies between 25 and 200 Hz. Finally, the exponential was taken to convert back ln(Power(f)) to Power(f) values. The power decay fit tracks broadband (i.e., frequency unspecific) power changes and was subtracted from the power spectrum to estimate the γ frequencies and detrended γ power (see above and Fig. 4).
Directionality
For estimating directionality relations, the respiration and LFP raw signals (see Fig. 1) or the time series of instantaneous frequency and power (see Figs. 3, 6, 7, 8) were first z-scored. Directionality was estimated by means of cross-correlograms (CCGs) and Granger causality. CCGs were obtained using the xcorr.m function with a maximal lag of 1 s, using the “coeff” option to normalize CCG values between −1 and 1 (equivalent to the correlation coefficient). Granger causality was computed through the Multivariate Granger Causality MATLAB Toolbox (Barnett and Seth, 2014). In Figure 1B, raw LFP and respiration signals were downsampled to 250 Hz before computing Granger causality. The Variational Autoregressive model order was fixed as 20 for all analyses. After computing Variational Autoregressive model parameters (tsdata_to_var.m) and associated auto-covariances (var_to_autocov), the pairwise conditional spectral Granger causality was obtained through the Multivariate Granger Causality toolbox function autocov_to_spwcgc.m. The overall, “time domain” Granger causality was obtained from integration (average) of the Granger spectrum through the function smvgc_to_mvgc.m.
Statistical analysis
Data are expressed as mean ± SEM unless otherwise stated. The significance of the difference between the mean Granger causality in each direction was assessed using paired t tests. Of note, for each significant Granger causality difference involving the time series of frequency and/or power values, we corroborated that such difference was also statistically significant against time-reversed data (Winkler et al., 2016). One-sample t tests were used to infer the significance of CCG lags against 0 ms. Repeated-measures ANOVA was used to assess for significant relations between LFP frequency or power and breathing frequency; for these analyses, we used 5 nonoverlapping breathing frequency bins: 1-3 Hz (2 Hz center frequency), 3-5 Hz (4 Hz center frequency), 5-7 Hz (6 Hz center frequency), 7-9 Hz (8 Hz center frequency), 9-11 Hz (10 Hz center frequency). Statistical significance was set at α = 0.05.
Results
Respiration phase Granger causes neocortical oscillations at the same frequency as breathing
We recorded parietal cortex activity of mice during REM sleep while mice were placed in a plethysmograph to simultaneously monitor their respiratory activity (Fig. 1A). REM episodes were identified by behavioral immobility and prominent θ activity in the cortical field potential (for details, see Materials and Methods). The power spectrum of the plethysmograph-derived respiration signal exhibited a peak at 3.2 ± 0.35 Hz (mean ± SD, n = 22 mice) with a long tail toward faster frequencies up to values ∼11 Hz. At the same time, the neocortical LFP displayed, as expected, a prominent peak at θ frequency (7.7 ± 0.37 Hz; Fig. 1B). In addition, a smaller power bump at a similar frequency as the respiration power peak could also be observed in the LFP spectrum (Fig. 1B); moreover, the phase coherence spectrum between respiration and LFP peaked at the breathing frequency (Fig. 1C). Therefore, combined, the power and coherence spectra show that, in addition to θ oscillations, the parietal cortex LFP also displays the so-called respiration-entrained rhythm (RR) during REM sleep, consistent with previous reports (Zhong et al., 2017; Tort et al., 2018b). Also consistent with these reports, there was no phase coherence between the parietal field potential and respiration at θ frequency (Fig. 1C). Furthermore, there was also no phase coherence when the LFP of one animal was paired with the respiration signal of another animal (Fig. 1C); this surrogate analysis thus shows that the actual respiration-LFP coherence is not because of the power peaks of the individual signals.
After corroborating RR presence in the parietal cortex, we next characterized directionality relations between respiration and LFP signals. To that end, we computed both CCGs and Granger causality between their raw traces (see Materials and Methods). Consistent with the phase coherence spectrum (Fig. 1C), the CCGs exhibited an oscillatory pattern with a period matching the respiration cycle (321 ± 30 ms [mean ± SD; n = 22 mice]; Fig. 1D). Interestingly, using the respiration signal as reference, the average CCG exhibited a negative peak at 110 ms, indicating that the LFP phase entrainment to respiratory cycles lags the respiration signal (Fig. 1D). It should be noted, however, that because of the rhythmical nature of the CCG in this case, directionality relations cannot be unambiguously determined (Bastos and Schoffelen, 2015). We therefore calculated the Granger causality spectrum between raw LFP and respiration signals, which revealed a prominent peak at the breathing frequency specifically for the respiration → LFP direction (Fig. 1E, left). Consistently, the overall, time domain Granger causality was statistically significantly higher in this direction than in the LFP → respiration direction (t(21) = 4.06, p = 0.00056, paired t test; Fig. 1E, right). As for the coherence spectrum, there was no meaningful Granger causality relation for surrogate pairings of LFP and respiration signals (Fig. 1E, left). Together, these results corroborate previous observations that respiration drives neuronal network oscillations of the same frequency (Nguyen Chi et al., 2016; Karalis and Sirota, 2018).
θ frequency and amplitude depend on breathing frequency
Respiration is dynamically modulated (Fig. 1A,B), similar to other, brain-endogenous rhythms, such as θ or γ oscillations. This poses the question of whether there is any influence of breathing frequency on oscillatory network activity or vice versa. We therefore tracked the breathing frequency during REM sleep by computing time-frequency spectrograms of respiratory activity using 500 ms windows with no overlap. For each window, to aid visualization, the power distribution was normalized by its maximal value, and the instantaneous breathing frequency was defined as the frequency of the maximal power value (which equals 1 because of the normalization; Fig. 2A,B). The breathing frequency typically varied between 2 and 5 Hz during the analyzed REM sleep periods, although faster breathing up to 11 Hz could also be observed (Fig. 2A-C). A similar procedure was performed to track the peak frequency of θ in the parietal cortex LFP signal for the same 500 ms windows. Interestingly, the histogram distribution of respiration and θ peak frequencies for the pool of 500 ms windows (Fig. 2D) resembled the average power spectra (compare Fig. 1B). This shows that the respiration power tail toward higher frequencies (Fig. 1B) corresponds to periods in which the animals indeed breathed at faster rates, as opposed to being caused by power leakage from lower breathing frequencies.
We next grouped time windows with similar instantaneous breathing frequency (Fig. 2A,B) into 10 groups (2 Hz bandwidth with 1 Hz overlap: i.e., 1-3, 2-4, 3-5 Hz, etc.). Figure 2E shows power spectra from data for each instantaneous frequency (top) together with the mean LFP power spectra for the same “respiration-binned” window subsets (bottom) in one animal. This alignment revealed a major increase in both θ frequency and θ power with increasing breathing frequency (Fig. 2E, bottom). This result was very robust at the group level of 22 mice, for which repeated-measures ANOVA showed a highly significant effect of breathing frequency on both θ frequency and power (Fig. 2F; θ frequency: F(4,103) = 113.32, p < 10−36; θ power: F(4,103) = 24.3, p < 10−13). Therefore, in addition to the previously described cycle-by-cycle phase-locking effects (Fig. 1), breathing does also relate to oscillatory brain activity through its instantaneous frequency (Fig. 2F).
Changes in θ activity precede changes in breathing frequency
The results in Figure 2F show that both instantaneous frequency and amplitude of θ exhibit a strong positive relationship with the instantaneous breathing frequency. In order to assess potential temporal relations, we next investigated directionality between these oscillatory features. To that end, we performed a similar approach as above, but this time, to allow for CCG and Granger causality calculations, we used a 90% overlap among the 500 ms windows (i.e., 50 ms step size). Therefore, for the directionality analyses, we used time series of instantaneous frequency and/or power sampled at 20 Hz. Surprisingly, the average CCG computed using the breathing frequency as the time lag reference showed that the changes in θ frequency occurred before the changes in breathing frequency (−125 ± 53 ms [mean ± SD], t(21) = −11.08, p < 10−9, one-sample t test against 0 ms; Fig. 3A). This finding was corroborated by Granger causality analysis in the frequency and time domain (Fig. 3B), in which causality in the direction θ frequency → respiration frequency was statistically significantly higher than in the opposite direction (t(21) = −4.77, p = 0.0001, paired t test). Thus, the instantaneous frequency of θ oscillations is Granger-causal for the instantaneous frequency of respiration. We note that this result appears to contrast the opposite Granger causality between the raw respiration and parietal cortex LFP signals shown in Figure 1E. Here, however, we analyzed directionality relations between the time series of instantaneous frequencies, a derived feature from the raw signals reflecting the variance in leading frequency, rather than the actual oscillation itself. In line with this notion, the high Granger causality at slow frequencies in Figure 3B (left) does not correspond to respiration or θ frequencies but is likely because of the slower time scale on which variations of θ and breathing occur (see Fig. 3B, left, inset plot).
We then asked whether variations in θ power would also be temporally related to respiration. Using the time series of instantaneous breathing frequency as reference, we found that θ power changes tended to lead breathing frequency, as inferred by both the mean CCG lag (−80 ± 10 ms [mean ± SD], t(21) = −3.54, p = 0.0019, one-sample t test against 0 ms; Fig. 3C) and Granger causality (t(21) = −2.12, p = 0.046, paired t test between time domain Granger causality levels in each direction; Fig. 3D). The magnitude of this effect was, however, much weaker than the directionality of variations in frequency (t(21) = 8.11, p < 10−7, paired t test between Granger causality levels). In all, as opposed to the phase entrainment relations, which are driven by respiration (Fig. 1), these results show that the instantaneous breathing frequency does not lead, but actually follows, intrinsic changes in brain activity in the θ range.
Slow-γ frequency increases with breathing frequency
In addition to RR and θ oscillations (Fig. 1B), the parietal cortex was also characterized by activity in two γ sub-bands during REM sleep. These could be clearly inferred by the power bumps at ∼50 to ∼80 Hz and ∼120 to ∼160 Hz in the LFP spectrum, especially when plotted at a logarithmic scale (Fig. 4A, left). In this work, we refer to these frequency ranges as slow- and fast-γ activity, respectively.
We first call attention to a technical point that will influence the interpretation of the results: to properly track the activity of genuine γ oscillations, as opposed to unspecific activity in the γ range, we estimated γ frequency and power for each sub-band after first subtracting the power decay fit from the LFP power spectrum (see Fig. 4A,B, right panels). Therefore, by “genuine” activity, we mean the excess power in the analyzed γ sub-band above the power decay level (also referred to as “detrended power”). This approach controls for broad, unspecific changes of overall power at fast frequencies. In any case, to allow comparison with previous approaches (e.g., Chen et al., 2011; Zheng et al., 2015), we also present results obtained from the absolute (i.e., unsubtracted) power spectra, regardless of the presence or not of prominent γ power bumps.
Figure 4B shows power spectra computed for REM sleep epochs binned by respiration frequency as in Figure 2E, but for another example animal. The middle panel shows the absolute LFP power, and the right panel shows the detrended power. Note in these panels an increase in the peak frequency of both slow- and fast-γ oscillations for epochs associated with faster breathing (dark green). The relation between slow- and fast-γ power with breathing frequency, on the other hand, depended on how γ power was estimated. Considering the absolute power level (Fig. 4B, middle), there was a clear increase during faster breathing for both γ sub-bands. However, this apparent power increase did not seem to be specific for the γ sub-bands since there was a broadband shift toward higher power levels for all fast frequencies. Indeed, this was reflected in the upshift also observed for the power decay fits (see Fig. 4B, insets). When removing this general trend by subtracting the corresponding broadband power fit of each spectrum, the relation between detrended γ power and respiration frequency became less clear (Fig. 4B, right).
A systematic analysis of group data revealed that the increase of slow-γ frequency with breathing frequency was statistically significant (F(4,103) = 8.53, p < 0.0001, repeated-measures ANOVA), while no significant relation was observed for the fast-γ frequency (F(4,103) = 1.26, p = 0.29) (Fig. 5A,B, left). Regarding γ power, the increase in the absolute power of both slow- and fast-γ with breathing frequency was highly statistically significant (slow γ: F(4,103) = 30.77, p < 10−16; fast γ: F(4,103) = 70.69, p < 10−27) (Fig. 5A,B, middle). On the other hand, detrended slow- and fast-γ power exhibited large variability among animals and were not statistically significantly related to breathing frequency (slow γ: F(4,103) = 0.35, p = 0.84; fast γ: F(4,103) = 0.05, p = 0.99) (Fig. 5A,B, right).
Changes in slow-γ frequency and power follow changes in breathing frequency
Similar to θ oscillations, we next analyzed directionality relations between breathing frequency and the instantaneous frequency of both γ sub-bands. For slow-γ oscillations, using respiration as reference, the CCG peaked at 114 ± 15 ms (mean ± SD), meaning that the increases in instantaneous slow-γ frequency lagged the increases in breathing frequency (t(21) = 3.54, p = 0.0019, one-sample t test against 0 ms; Fig. 6A, top). Granger causality analysis revealed a consistent finding, in which causality levels in the respiration → slow-γ direction were statistically significantly higher than in the slow-γ → respiration direction (Fig. 6B, top; t(21) = 8.68, p < 10−7, paired t test). For the fast-γ band, the CCG exhibited a peak at −202 ± 115 ms (t(21) = −8.25, p < 10−7, one-sample t test against 0 ms; Fig. 6A, bottom), suggesting a possible lead by fast-γ; however, the Granger analysis showed no difference in causality directions (Fig. 6B, bottom; t(21) = −0.36, p = 0.72, paired t test). In addition, Granger causality levels for fast-γ and respiration frequency were much lower than the causality levels found between respiration and slow-γ frequency (t(21) = −10.16, p < 10−8, paired t test).
Finally, we investigated directionality relations between γ power and breathing frequency. For the absolute power values (Fig. 7A,B), we found that the breathing frequency was Granger-caused by both slow-γ and fast-γ power (slow γ: t(21) = −6.07, p < 10−5; fast γ: t(21) = −9.82, p < 10−8, paired t tests between directions; Fig. 7B), a result confirmed by inspection of the CCG peak lags (slow γ: −175 ± 72 ms [mean ± SD], t(21) = −11.41, p < 10−9; fast γ: −120 ± 45 ms, t(21) = −12.44, p < 10−10, one-sample t tests against 0 ms; Fig. 7A). However, we did not take this result at face value because some corollaries called our attention. For one, it sounded odd that slow-γ power would drive the breathing frequency while its frequency would follow (compare Fig. 6). For another, it also seemed strange that Granger causality levels were much higher for the fast-γ band, which usually has the lowest signal-to-noise ratio; indeed, fast-γ Granger causality was even higher than θ power Granger causality levels (compare Fig. 3D). We thus suspected that these causality relations could be because of broadband, unspecific variations in the absolute power of fast frequency activity. Consistent with this possibility, in Figure 4 the overall power of the fitting of the power decay with LFP frequency increases as a function of respiration frequency. Therefore, we next investigated for temporal relations between breathing frequency and the mean level of broadband power, as inferred by the mean value of the power decay fit (Fig. 7C). The results revealed that the changes in broadband LFP power levels are highly Granger-causal to breathing frequency (t(21) = −14.45, p < 10−11, paired t tests between directions), and even more so than slow- and fast-γ power (F(2,63) = 54.57, p < 10−13, one-way ANOVA; compare the right panels in Fig. 7B,D). Thus, we conclude that the Granger causality between the absolute γ power and respiration frequency is not specific to the γ sub-bands but because of broadband changes in the power of fast LFP activity.
The picture was much different when we analyzed the detrended power levels (Fig. 8). In this case, the Granger causality relations resembled those found for the instantaneous frequency: namely, detrended slow γ power changes were Granger-caused by changes in breathing frequency (t(21) = 2.80, p = 0.0108, paired t test between directions; Fig. 8B, top), whereas detrended fast γ power was not related to breathing frequency under the Granger causality analysis (t(21) = 0.75, p = 0.46, paired t test between directions; Fig. 8B, bottom). Interestingly, and although the CCG lag was not statistically different from zero (−11 ± 171 ms [mean ± SD], t(21) = −0.31, p = 0.76), the respiration-slow γ CCG was negative (Fig. 8A), suggesting an inverse relationship between genuine slow γ power and respiration. In all, we conclude that changes in respiration frequency precede and Granger-cause changes in slow-γ frequency as well as slow-γ power after correcting for broad, frequency-unspecific power shifts. This finding contrasts with the opposite temporal relation found between respiration and θ, in which the LFP θ activity leads respiration (compare Fig. 3).
Combined, therefore, our result suggests that changes in the frequency of cortical θ oscillations precede and Granger-cause changes in respiration frequency (Fig. 3). In contrast, the real-time respiration-entrained LFP component (Fig. 1D,E) and slow-γ activity (Figs. 6 and 8) depend on, and follow, nasal breathing.
Discussion
The effect of respiration on brain activity has been the subject of much recent work. It is by now clear that respiration adds a breathing-frequency component to neuronal network oscillations (for reviews, see Tort et al., 2018a; Heck et al., 2019). Importantly, this effect occurs in widespread brain regions, far beyond olfactory networks, a finding reproduced here (Fig. 1). The present work reveals a new, much broader interaction between respiration and different network oscillations. Both activities do not only interact through cycle-by-cycle phase entrainments but also through mutual dependencies involving the instantaneous frequency. Specifically, we show that the frequency and amplitude of θ and γ oscillations detected in the neocortex covary with breathing frequency. Moreover, our directionality analyses revealed coherent temporal relations among the electrophysiological and respiratory rhythms, in which changes in θ activity precede changes in breathing frequency that, in turn, precede changes in slow γ. Therefore, the present study unveils novel relationships between respiration and brain activity beyond the well-established phase-locking relations.
A myriad of recent studies convincingly showed that the reafferent sensory signal brought about by nasal breathing drives LFP oscillations phase-locked to respiration (RR). This has been demonstrated in multiple ways, including experimental (e.g., naris closure, olfactory epithelium ablation, bulbectomy, and tracheotomy) (Ito et al., 2014; Yanovsky et al., 2014; Lockmann et al., 2016; Bagur et al., 2018; Karalis and Sirota, 2018; Moberly et al., 2018) and computational approaches (e.g., Granger causality analysis) (Nguyen Chi et al., 2016; Karalis and Sirota, 2018). Here we show different relations between variations in instantaneous frequency of network oscillations and respiration. Indeed, cortical LFP activity in the θ range Granger-causes respiration frequency, rather than vice versa. Nevertheless, we believe that this result is consistent with a causal role of a central drive underlying changes in both network and breathing activity. Our working hypothesis, illustrated in Figure 9, is that REM-controlling neuronal circuits in the brainstem send signals to θ-generator circuits and to nearby respiration-controlling nuclei. Thereby, they determine changes in θ as well as respiration frequency during REM sleep. Therefore, in this hypothesized scenario, the central drive for accelerating θ and respiration would be the same. However, the changes in θ activity would precede changes in respiration frequency, and thus look Granger-causal to them under directionality analyses because of a presumed faster activation of θ-generating circuits, located a few synapses away, compared with the longer delays required for modulating respiration, which include modulation of the underlying networks, conduction delays to the diaphragmatic muscle, and subsequent thoracic volume changes. Such a model is compatible with the raw respiration signal being Granger-causal to RR (Fig. 1) since this is determined by activation of the olfactory bulb by the passage of air through the nostrils (Fig. 9) regardless of changes in respiration and θ frequencies. On a technical aspect, we note that here we refer to oscillations that can be empirically detected in the parietal cortex (e.g., RR or θ) regardless of their true regional origin (i.e., regardless of volume conduction) (see also Tort et al., 2018b).
Schematics of hypothesized system interactions leading to the current results.
Consistent with the model in Figure 9, there is experimental evidence for mutual connections among the brainstem microcircuits controlling REM sleep, arousal, and respiration (Orem, 1980; Luppi et al., 2012; Boutin et al., 2017; Yackle et al., 2017; Benarroch, 2018; Del Negro et al., 2018). Interestingly, there are also connections between breathing centers and the generators of other orofacial rhythms, such as whisking (Moore et al., 2013; Kleinfeld et al., 2014; McElvain et al., 2018), which may become active during REM sleep (Robinson et al., 1977; Tiriac et al., 2012). Therefore, the interplay among brainstem nuclei is likely to control higher arousal states within REM sleep that would be associated with faster θ and respiration frequencies.
A possibility is that these activity states play a role in processing emotional information, one of the main functions attributed to REM sleep (van der Helm et al., 2011; Groch et al., 2013; Hutchison and Rathore, 2015; Tempesta et al., 2018), especially considering that experiencing emotions has been related to changes in breathing patterns (Grossman and Wientjes, 2001; Philippot et al., 2002; Q. Zhang et al., 2017). But regardless of the information content during these states, a recent study has shown that the respiration frequency is the main modulator of cerebral oxygenation (Q. Zhang et al., 2019), pointing to a close link between breathing and metabolic brain demands. Therefore, the higher oxygenation brought about by faster breathing indicates greater energy consumption, suggesting network states of heightened information processing.
Regarding potential behavioral or cognitive consequences of changes in breathing frequency, it is worth noticing that differences in nasal flow rate influence olfactory processing (Buonviso et al., 2006; Courtiol et al., 2011a,b; Esclassan et al., 2012). Based on recent studies, it is well feasible that breathing frequency changes could also influence nonolfactory processes (Zelano et al., 2016; Varga and Heck, 2017; Arshamian et al., 2018; Herrero et al., 2018). Interestingly, using the same dataset analyzed here, we have recently shown that θ-γ coupling, a hallmark network activity previously associated with cognitive processes (Tort et al., 2009; Canolty and Knight, 2010; Lisman and Jensen, 2013), depends on the breathing frequency through an inverted V-shaped relation (Hammer et al., 2020). This novel result suggests that indeed changes in breathing may influence nonolfactory information processing. As a technical remark, we note that the phenomenon of θ-γ coupling could not be subjected to the same directionality analyses as investigated here for the individual oscillations. This is because assessing cross-frequency coupling requires using longer epoch lengths (i.e., sampling several θ cycles to average out noise) (Tort et al., 2010; Scheffer-Teixeira and Tort, 2018) and cannot be performed in the same sub-second time windows used in the present study to assess the instantaneous frequency. Another technical aspect is that, differently from Hammer et al. (2020), here we were unable to analyze phasic REM sleep separately since this sub-state is characterized by short bouts of noncontiguous activity (thus not suitable for Granger analysis). Of note, since phasic REM sleep only comprised 2.48% of the analyzed data, the results were largely the same when analyzing tonic REM sleep separately, except for the θ power influence over respiration frequency, which became nonsignificant (data not show). In addition to influencing θ-γ coupling, another recent study indicated that the breathing frequency may also affect LFP phase entrainment by respiration (Girin et al., 2020; see also Jessberger et al., 2016), suggesting a possible interplay between frequency and phase-driven effects.
Strikingly, while the instantaneous breathing frequency followed changes in θ (Fig. 3), the same directionality analyses revealed that respiration was Granger-causal to slow-γ activity (Figs. 6 and 8). This result suggests that neocortical slow-γ is likely to depend on the olfactory bulb activity, which is majorly driven by sensory stimulation of the olfactory epithelium (Buonviso et al., 2006; Courtiol et al., 2011a,b; Esclassan et al., 2012; Short et al., 2016). It is worth noticing that this drive is also provided by the mechanical stimulation through nasal air movement, even in the absence of odors (Grosmaitre et al., 2007; Connelly et al., 2015). The olfactory bulb is well known to express prominent γ activity (Ravel et al., 2003; Martin et al., 2004; Beshel et al., 2007; Cenier et al., 2008, 2009; Kay et al., 2009), and this region has been shown to produce more than one type of γ with different sensitivities to respiration (Kay, 2003; Zhong et al., 2017; Zhuang et al., 2019). Interestingly, original accounts on γ activity related this rhythm to olfaction (Adrian, 1950; Bressler and Freeman, 1980; Rojas-Líbano and Kay, 2008), even when recorded outside primary olfactory regions (Vanderwolf, 1992, 2001). While the present study measured γ over the parietal neocortex, it would be interesting to study its relation to olfactory bulb γ. Indeed, a recent report indicates that olfactory bulb γ is causal to γ in the medial PFC (Karalis and Sirota, 2018).
As opposed to slow γ, genuine (detrended) fast-γ power was not as temporally related to the respiration frequency (Figs. 6 and 8). It is worth highlighting that the parietal recordings exhibited a clear excess power in the fast-γ range (a “power bump,” see Fig. 4A). Moreover, similar to previous characterizations of neocortical activity during REM sleep (Scheffzük et al., 2011, 2013; Brankačk et al., 2012), θ-fast γ coupling was prominent in the present dataset (Hammer et al., 2020). In other words, although the power levels at the fast frequencies of the LFP spectrum may reflect unspecific, aperiodic activity, including EMG contamination, the current data did exhibit a true oscillatory activity in the fast-γ range (i.e., a genuine fast LFP rhythm). But despite the unequivocal presence, fast-γ oscillations were not related to respiration, suggesting they could underlie internal computations not affected by respiratory reafference. Nevertheless, although the individual oscillatory features were not affected, we have recently shown that the higher-order phenomenon of cross-frequency coupling between θ and fast-γ depends on breathing frequency (Hammer et al., 2020), which underscores the complexity of respiration-brain interactions.
Previous studies in awake animals have shown that LFP activity in the θ and γ ranges depend on running speed. The relationship between θ and speed is in fact an old one; since the first recordings, both θ power and θ frequency have been consistently shown to positively correlate with locomotion speed (McFarland et al., 1975; Sławińska and Kasicki, 1998; Hinman et al., 2011). More recently, speed dependencies have also been reported for the γ sub-bands (Chen et al., 2011; Ahmed and Mehta, 2012; Kemere et al., 2013; Zheng et al., 2015; Lopes-Dos-Santos et al., 2018). For instance, Chen et al. (2011), working with mice, reported a gradual amplitude increase with speed for both slow- and fast-γ oscillations recorded in CA1. The same group later demonstrated changes in the instantaneous γ frequency with speed in the rat CA1 (Ahmed and Mehta, 2012). Interestingly, this study also reported that slow γ power negatively correlates with speed, a result subsequently corroborated by independent laboratories (Kemere et al., 2013; Zheng et al., 2015; but see Lopes-Dos-Santos et al., 2018). While we recorded from immobile mice during REM sleep, which is a particular brain state characterized by a unique pattern of cortical activity and excitation/inhibition balance (Niethard et al., 2016), the speed dependencies of LFP oscillations described in awake animals are reminiscent of the θ and γ relations with the breathing frequency observed here, especially considering that faster locomotion speeds are associated with faster breathing (Q. Zhang et al., 2019). Future research that simultaneously records respiration and locomotion activity in suitable spatial arenas should help to disentangle the contribution of each factor to modulating LFP activity. In particular, from the present results, we already know that changes in LFP activity with breathing may occur without locomotion changes. Likewise, it would be interesting to characterize the influence of speed after controlling for changes in breathing patterns.
Some technical aspects should be observed, though, such as the way power values are estimated, which may substantially vary among studies (i.e., absolute vs relative power, correction or not for broadband power shifts, use of independent components or empirical, nonlinear decompositions, etc.) (for a recent discussion on this subject, see Donoghue et al., 2020). As illustrated here, the precise method for estimating γ power may critically influence the results (compare Figs. 7 and 8): while the absolute γ power was Granger-causal to respiration frequency, such an effect was not evident when correcting for broadband power changes. Indeed, the latter could much better account for the changes in respiration frequency than the power levels in the γ sub-bands. Of note, the cause of broadband changes in the power of fast LFP frequencies has yet to be determined; we speculated that EMG contamination could play a role but found no conclusive evidence in a group of 6 extra animals with EMG recordings (data not shown).
Another technical aspect is the recorded region (neocortex vs hippocampus). While θ is considered a global rhythm, γ represents more local activity and may thus exhibit different dynamics among regions. Indeed, Zheng et al. (2015) showed differences in the speed dependency of γ activity between the entorhinal cortex and the hippocampus. To further complicate matters, Lopes-dos-Santos et al. (2018) called attention to potential differences in speed modulations of hippocampal γ activity between rats and mice. Finally, another technical aspect to be kept in mind is the very definition of γ activity. For instance, the pattern that we call fast-γ in the present work is particularly prominent in the neocortex during REM sleep (Scheffzük et al., 2011, 2013; Brankačk et al., 2012); and, despite the similar nomenclature, it is likely not to correspond to the “fast-γ” activity observed in the hippocampus of awake animals (Belluscio et al., 2012; Schomburg et al., 2014). In this sense, while we speculate that the same respiration frequency relations observed here should hold true for θ detected in other brain regions (given its long-distance coherence in mice) (Tort et al., 2018b), we are less sure about the generality of our results regarding γ oscillations.
In addition to the respiratory system, many other physiological systems display ultradian rhythmicity, including organs of the cardiovascular and digestive systems (Gray et al., 2009; Babo-Rebelo et al., 2016; Richter et al., 2017). Interestingly, recent studies have also started to unveil relationships between oscillatory activity produced by these systems and the brain. For instance, Tallon-Baudry and collaborators have shown that ascending signals from the heart and gastrointestinal tract modulate brain dynamics and affect perception, emotion, and cognition (Babo-Rebelo et al., 2016; Richter et al., 2017; Rebollo et al., 2018; Azzalini et al., 2019). In particular, the phase of the gastric cycle has been shown to modulate the human alpha rhythm (Richter et al., 2017). Therefore, it is increasingly evident that brain activity is coupled to autonomously controlled regions of the body, possibly reflecting interactions between central pattern generators and the sympathetic/parasympathetic nervous system (Guyenet, 2014). Nevertheless, the specific type of coupling (e.g., phase-amplitude coupling, n:m phase-locking, frequency-frequency relations) and oscillatory frequencies involved (e.g., θ, α) seem to differ among species and organ systems.
In conclusion, the present results show that the breathing frequency temporally relates to θ and γ network oscillations and that the relationship between respiration and neocortical brain activity goes beyond phase entrainment relations. These findings motivate the investigation of breathing frequency relations in other brain regions, in different behavioral states and cognitive loads, and at different spatial scales, including the cellular level. Moreover, given the link between breathing and γ shown here, it will be interesting to characterize further how specific sub-bands of γ oscillations in different brain regions relate to olfactory bulb γ and respiration.
Footnotes
This work was supported by Deutsche Forschungsgemeinschaft (SFB 1134/A01; Dr 326/10-1), Bundesministerium für Bildung und Forschung (German-Brazil Cooperation Grant 01DN12098), the Brazilian National Council for Scientific and Technological Development, the Brazilian Coordination for the Improvement of Higher Education Personnel, and the Alexander von Humboldt Foundation. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
The authors declare no competing financial interests.
- Correspondence should be addressed to Adriano B.L. Tort at tort{at}neuro.ufrn.br or Andreas Draguhn at andreas.draguhn{at}physiologie.uni-heidelberg.de
