Verification of a Central Pacemaker in Brain Stem by Phase-Coupling Analysis Between HR Interval- and BOLD-Oscillations in the 0.10–0.15 Hz Frequency Band

The origin of slow intrinsic oscillations in resting states of functional magnetic resonance imaging (fMRI) signals is still a matter of debate. The present study aims to test the hypothesis that slow blood oxygenation level-dependent (BOLD) oscillations with frequency components greater than 0.10 Hz result from a central neural pacemaker located in the brain stem. We predict that a central oscillator modulates cardiac beat-to-beat interval (RRI) fluctuations rapidly, with only a short neural lag around 0.3 s. Spontaneous BOLD fluctuations in the brain stem, however, are considerably delayed due to the hemodynamic response time of about ∼2–3 s. In order to test these predictions, we analyzed the time delay between slow RRI oscillations from thorax and BOLD oscillations in the brain stem by calculating the phase locking value (PLV). Our findings show a significant time delay of 2.2 ± 0.2 s between RRI and BOLD signals in 12 out of 23 (50%) participants in axial slices of the pons/brain stem. Adding the neural lag of 0.3 s to the observed lag of 2.2 s we obtain 2.5 s, which is the time between neural activity increase and BOLD increase, termed neuro-BOLD coupling. Note, this time window for neuro-BOLD coupling in awake humans is surprisingly of similar size as in awake head-fixed adult mice (Mateo et al., 2017).

Introduction

Slow BP oscillations in the frequency range around 0.1 Hz are known as Mayer waves and well documented [see, review by Julien (2006)]. However, it is still unknown how they are generated. Baroreflex mechanisms and intrinsic activity of sinoatrial node cells may play a dominant role (Van Roon et al., 2004; Moen et al., 2019), but also a central (neural) pacemaker may contribute independently to their generation (Preiss and Polosa, 1974; Zhang et al., 1998; Eckberg et al., 2016; Ghali and Ghali, 2020). This has been controversely discussed a couple of years ago (Eckberg and Karemaker, 2009). Support for the central pacemaker hypothesis comes from the work of Perlitz et al. (2004). They observed a neural brain stem rhythm with frequencies around 0.15 Hz (emerging in reticular neurons in lower brain stem) which was phase coupled with respiration and heart rate beat-to-beat interval (RRI) fluctuations in dogs as well as humans (e.g., Lambertz and Langhorst, 1998; Perlitz et al., 2004). The brain stem seems to be a reasonable location for a central pacemaker, because it accommodates the cardiovascular and respiratory centers and is the origin of the ARAS (Moruzzi, 1958), a network arising from the lower brain stem. The ARAS is dynamically organized and modulates distant cortical areas. Several studies indicate that the “0.15-Hz rhythm” of Perlitz et al. (2004) may not be confined to this single frequency, but comprises slow rhythms with frequency components of 0.15 ± 0.03 Hz. Kuusela et al. (2003) analyzed BP and HR signals and observed two distinct frequency components, one around 0.08 Hz and another around 0.12 Hz. Moreover, Rassler et al. (2018) found alternating 0.10-Hz and 0.15-Hz waves in RRI and respiration signals, and Yuen et al. (2019) reported on blood oxygenation level-dependent (BOLD) oscillations that most likely are associated with vasomotor activity centered at 0.15 Hz. Together, these studies document the importance of frequency components >0.1 Hz in the brain.

In the present study, we analyze time delays (using the phase locking value, PLV; Lachaux et al., 1999) between RRI fluctuations and BOLD signals in the frequency range of 0.10–0.15 Hz. Due to the close interaction between brain and heart (see review, Thayer and Lane, 2009), slow changes in neural activity are linked to slow RRI changes. Thus, RRI oscillations can be used to index central nervous system activity. Because neural oscillations cannot be recorded from the brain stem of humans, we instead record BOLD signals from the brain stem and predict the time lag on the basis of the following working hypothesis:

If slow rhythmic neural activity (central pacemaker) exists in brain stem, then two responses are expected: slow oscillating cardiac RR intervals and slow oscillating BOLD signals in brain stem. Specified: If a neural pacemaker generates oscillations in a particular frequency band (e.g., 0.10–0.15 Hz), it should modulate slow cardiac interval (RRI) fluctuations rapidly due to fast vagal innervation (within 0.2–0.3 s; Smith, 1974). Moreover, it should affect BOLD fluctuations in the brain location of interest, which are, however, delayed by a few seconds due to the hemodynamic response. This hemodynamic response is often termed neurovascular coupling (Arthurs and Boniface, 2002), however, this term can be misleading. Therefore, it is suggested to use the term neuro-BOLD coupling (NBC), which defines the time between neural activity (e.g., gamma power) increase and BOLD increase. Mateo et al. (2017) reported a NBC time of 2.6 s in mice. This means, the pattern of the cardiac interval time series (RRI signal) and the pattern of the BOLD signal should coincide wave-by-wave with a time-shift of 2–3 s.

If this hypothesis is correct, a method is available to assess sources of neural oscillations in the brain stem indirectly by analyzing the phase coupling between RRI and BOLD signals in the sense of a surrogate indicator.

Taken together, the two main questions guiding this research were: (i) Does a “central pacemaker” exist in the brain stem and contribute to the generation of RRI oscillations? (ii) Have BOLD signals in a frequency range above 0.10 Hz a neural origin?

Materials and Methods

Participants

A total of 23 participants (12 female, 22 right-handed) between 19 and 34 years (mean ± SD: 24 ± 3.2 years) took part in the study. All were naïve to the purpose of the study, had no former MRI experience, had normal or corrected-to-normal vision and were without any record of neurological or psychiatric disorders, as assessed by self-report. All participants gave informed written consent to the protocol of the study, which was approved by the local Ethics Committee at the University of Graz (number: GZ. 39/75/63 ex 2013/14).

The experimental task consisted of four resting states (R1, R2, R3, R4) and four within-scanner questionnaires (AS1, AS2, AS3, AS4), carried out in two sessions separated by about 50 min. The first session (a) started with the first questionnaire (AS1) and was followed by the first resting state (R1). Thereafter, two motor tasks, unrelated to this study, were performed. The first session ended with the second resting state (R2) and second questionnaire (AS2). Filling out each questionnaire took approximately 5 min and each resting state lasted for about 350 s. The second session (b) was quasi a duplicate of the first one, with two resting states (R3, R4) and two within-scanner questionnaires (AS3, AS4). Individuals were requested to keep their eyes open, stay awake, and avoid movements during the resting states.

State anxiety was assessed with the state-trait anxiety and depression inventory (STADI; Laux et al., 2013), which was presented on a screen within the scanner. The STADI is an instrument constructed to assess both state and trait aspects of anxiety and depression. It is based on the State-Trait Anxiety Inventory (Spielberger et al., 2009), but allows a reasonable separation of anxiety and depression symptoms. Items were answered with a trackball following each resting state.

Recording of Physiological Signals (ECG, Respiration) and RRI Time Courses

Electrocardiogram and respiration were recorded inside the scanner using the PMU (Siemens Physiological Measurement Unit). For the positioning of the ECG electrodes on the thorax, standard channels (Siemens Standard, lead 1) were used. The respiratory data were acquired using a pneumatic cushion connected via an air hose to a pressure sensor on the recording unit. The cushion was attached to the subject by using a respiration belt. This technique allowed a wave-by-wave determination of respiratory rate. Sampling rate was 400 Hz. QRS detection and subsequent computation of RRI time series was performed using fMRI plug-in for EEGLAB (Niazy et al., 2005). Each interval was referred to the previous R-wave (instantaneous interval; De Boer et al., 1985). In order to assure a reasonable quality of the ECG recording during a fMRI high scanning rate, RRI signals were re-processed using the Kubios HRV Premium Package (Kubios Ltd., Finland; version 3.0.2; Tarvainen et al., 2014). The RRI time series were subsequently resampled to the same sampling frequency of the BOLD time series. For the display and identification of local peaks in BOLD and RRI signals a resampling rate of 10 Hz was used.

Resting State fMRI and BOLD Signals

Resting-state functional images were acquired with a 3T scanner (Magneton Skyra, Siemens, Erlangen, Germany) using a multiband GE-EPI sequence (Moeller et al., 2010) with a simultaneous six-band acquisition with TE/TR = 34/871 ms, 52° flip angle, 2 × 2 × 2 mmł voxel size, 66 contiguous axial slices (11 × 6), acquisition matrix of 90 × 104 and a FOV of 180 × 208 mm², 400 volumes.

Pre-processing resorted to the DPARSF toolbox (Yan and Zang, 2010) and included removal of the first 10 volumes, slice timing correction adapted for multiband acquisitions (Woletz et al., 2014), rigid-body motion correction, normalization to Montreal Neurological Institute (MNI) space, resampling to 2-mm isotropic voxels, spatial smoothing with a 4-mm FWHM Gaussian kernel, linear detrending and extraction of time series from various regions of interest (ROIs) using the AAL atlas (Tzourio-Mazoyer et al., 2002).

Averaging of Breathing, RRI and BOLD Waves in Resting States

Averaging can be used to enhance the signal-to-noise ratio of physiological signals but it requires a trigger. Because stimuli-based triggers are not available in resting state data, the periodic positive peaks (maximal beat-to-beat intervals) of the RRI signal were used to define triggers (Pfurtscheller et al., 2017, 2019). This procedure includes the following steps: First, inspection of the RRI time course with slow spontaneous fluctuations and marking of a few positive peaks preferably occurring in intervals larger than 6 s and shorter than 12 s. Second, using these peaks as trigger to calculate averages (mean ± SE) of breathing, RRI and BOLD waves over a time frame between 6 s before and 6 s after the trigger (examples see Figures 1, 4).

FIGURE 1

Figure 1. RRI peak-triggered slow BOLD waves (mean ± SE) from 116 ROIs during rest. Respiratory artifacts of comparably large amplitude are prevalent in ROI 36 (right cingulum, supplied by anterior cerebral artery), ROIs 95, 96 (pons/brain stem, supplied by basilar artery) and ROIs 109–116 (vermis cerebelli, also supplied by basilar artery). The majority of ROIs exhibits artifact-free signals. Surprisingly, even a few ROIs in pons/brain stem (e.g., ROI 91 and 93) are of satisfying signal quality.

Search for Artifact-Free ROIs

For the anatomical leveling, the AAL atlas (Tzourio-Mazoyer et al., 2002) was used, and 116 ROIs were defined. The main benefit of the ROI over the voxels-level approach is the increased signal-to-noise ratio. Each ROI signal is computed by averaging the BOLD signal over hundreds of voxels. ROIs 1–90 correspond to cortical and subcortical structures, ROIs 91 – 108 to regions in the cerebellum including brain stem, and supplementary ROIs 109–116 to the vermis cerebelli. Due to the proximity of cerebellum and brain stem in axial slice and due to the contiguity of a major artery, it is reasonable to assume that the signals collected from AAL ROIs labeled “cerebellum” include the upper brain stem with the pons. Therefore, in this paper the denotation pons/brain stem is used for recordings from ROIs 91–108 instead of “cerebellum.”

ROIs with even numbers correspond to the right, and those with odd numbers to the left hemisphere. A major challenge is the vascular pulsing related to the main cerebral arteries, especially in the pons/brain stem. Therefore, a search for respiratory artifacts in all 116 ROIs was conducted through calculating averaged BOLD waves based on the slow waves in RRI signals. Images of 116 averaged BOLD waves (mean ± SE) from cortex, subcortex and pons/brain stem from one representative subject (9R1) with a spontaneous respiratory rate of ∼ 22 breaths/minute are shown in Figure 1. This Figure demonstrates that ROIs without respiratory artifacts are present in a few ROIs in pons/brain stem (e.g., ROIs 91 and 93) similar to ROIs in cortical areas (e.g., precentral gyrus, ROI 1) although the basilar artery is tight-fitting in the pons.

Based on this type of artifact detection, ROIs 95 and 96 were discarded, and ROIs 107 and 108 were excluded due to the small number of voxels (<200). From the ROIs representing cerebellum and brain stem (ROIs 91–108) 14 ROIs from pons/brain stem were considered for further analyses. These ROIs represent a distance of 16 mm between rostral (ROIs 91, 92) and caudal (ROIs 107,108) parts of pontine structure in upper brain stem (see sagittal and axial slices, Figure 2).

FIGURE 2

Figure 2. T1-images from medial, sagittal and axial slices from rostral pons ROIs 91/92 (left, Talairach space 0, 16, −28) and caudal pons ROIs 107/108 (right, Talairach space 0, 16, −44). The difference between ROIs 91/92 and 107/108 is 16 mm.

Calculation of Phase-Locking Between Oscillating RRI and BOLD Signals

Wavelet transform coherence was applied to the BOLD and RRI time series using the “Cross-wavelet and Wavelet Coherence” toolbox (Grinsted et al., 2004). We focus on the phase component, which allows to compute the PLV throughout the acquisition interval, except for small sections at the beginning and end where results are known to be unreliable (Torrence and Compo, 1998). PLV is a normalized measure of how much the phase difference between two signals changes in a user-chosen time window, regardless of the actual phase difference value; the reader is referred to Lachaux et al. (1999) and Pfurtscheller et al. (2017) for more details. The PLV was calculated in the band 0.10–0.15 Hz between RRI time course and the BOLD time series from the remaining 14 ROIs and ROI 1 (left precentral gyrus) for control. Two PLV-based parameters were computed: time delay (TD), which is the phase delay converted to time, and significant length of phase coupling (sigbin), which is the percentage of time samples within the time series that survive a 0.05 significance threshold (i.e., the null hypothesis of no phase-locking rejected with 95% certainty; Pfurtscheller et al., 2017). A positive time delay (pTD) indicates a lead and a negative delay (nTD) a lag of RRI oscillations compared to BOLD oscillations in the brain stem.

Neural processes (e.g., activity of a central pacemaker) are accompanied by changes in blood oxygenation and cerebral blood flow, both detectable through BOLD contrast measurements (Obrig et al., 2000; Buxton et al., 2004). BOLD is a complex function of cerebral metabolism, cerebral blood flow and cerebral blood velocity and sensitive to different types of artifacts as, for example, respiratory motion (Murphy et al., 2013).

The BOLD signals from pons/brainstem are expected to reflect the contribution of local vascular network and one major blood vessel, the basilar artery, which is close to the medulla and the pons. Similarly, important nuclei involved in cardiac control such as, for example, solitary tract nucleus and vagal motor neurons in the nucleus ambiguus and dorsal vagal motor nucleus are localized in the medulla and the pons (Verberne and Owens, 1998; Thayer and Lane, 2009). This means that two important sources for rhythmic BOLD signals are concentrated in the pons/brain stem: cerebral blood flow and neural activity changes.

Statistical Analysis

T-tests for independent measures were conducted in order to compare pTDs and nTDs between recording sites. Furthermore, the distribution of TDs was analyzed using Kolmogorov-Smirnov tests and Chi² tests. Finally, we compared the highest (HA) and lowest state anxiety (LA) across resting states for each individual using independent t-tests. Of note, TDs for each individual anxiety score were treated as independent samples. Strictly, these measures constitute dependent measures for each individual, however, aggregating reliable TDs (%sigbin > 10) across locations within individuals would result in a large number of missing cases, which severely attenuates statistical power. It should be noted though that our data analytic approach inflates degrees of freedom, thus constituting a rather liberal statistical approach. We deem this acceptable for a rather exploratory study design, but are aware that the findings warrant replication with larger sample sizes, thus allowing more conservative statistical approaches.

Results

Distribution of Time Delays Between Slow RRI and BOLD Oscillations in the Brain Stem

A time delay (TD) was considered reliable, when the associated significant length of phase coupling (sigbin) was equal to, or larger than 10% (percentage of samples within the time series; arbitrary threshold). TD values and sigbins from all subjects were summarized in two matrices with 23 rows (subjects) and 14 columns (ROIs) for each resting state. For each individual resting state and all states together, histograms of TDs were plotted, and differences were calculated between both brain sides (Table 1). Significant side differences were found for the first resting state (p = 0.039) and all four resting states together (p = 0.002) with more pTDs in the left side. The AS showed the highest level in the first resting state (R1) with AS = 20 ± 5 and a decline thereafter (Table 1). Interestingly, this high AS in R1 coincides with the left-sided predominance of pTDs in R1. The histograms across all resting states are presented in Figures 3A,B and show a clear dominance of pTDs.

TABLE 1

Table 1. Resting state (R), state anxiety (AS) and time delays in (s) for the left and right brain side for each resting state and all states together (1–4).

FIGURE 3

Figure 3. Histograms of TDs averaged over all 4 resting states, plotted separately for the left (A), and right brain side (B). Data from 7 ROIs and sigbin ≥10% are accumulated. ROIs with odd numbers represent the left, ROIs with even numbers the right brain side. Histograms representing pTDs (between 0.1 and 3 s) are indicated by dark gray color.

The histograms in Figures 3A,B are based on ∼ 800 TDs from 23 subjects, 4 resting states and 14 ROIs. They are clearly dominated by pTDs with largest accumulation (peaks) close to +2 s in both sides of the cerebellum/brain stem. The distributions of nTDs and pTDs were significantly different for both the left (Kolmogorov–Smirnov Z = 8.18, p < 0.001) and right side (Kolmogorov–Smirnov Z = 8.83, p < 0.001). Furthermore, a Chi² test revealed that the numbers of pTDs and nTDs differed by side (Chi² = 5.65, p = 0.02), suggesting that pTDs in the left side were overrepresented and nTDs underrepresented. Most importantly, our findings – as depicted in Figure 3 – show that none of the observed TDs was larger than +/−3 s. Remarkable is not only the predominance of pTDs on the left side but also the general dominance of pTDs in the pons/brain stem.

In summarizing, the analysis of TD histograms revealed the following findings: (i) No TDs were found larger than +/−3 s. (ii) A clear predominance of pTDs was found in the left pons/brain stem. (iii) Significant phase-coupling between RRI and BOLD oscillations was largest and most dominant across subjects in the first resting state and declined thereafter.

The Search for pTDs in the Brain Stem With Maximal Delays Between RRI and BOLD Oscillations

The accumulation of pTDs close to 2 s in the left pons/brain stem (Figure 3A) and its restriction to less than 3 s is a novel observation and raises questions not only about the incidence of pTDs across subjects, but also about the maximal possible time delay. In order to search for such a SSMTD, a threshold has to be defined. Based on the histogram across all resting states from the left side (Figure 3A) with a sharp drop of the columns to the right of the peak close to 2 s (dark gray), an arbitrary threshold of TD = 1.9 s was defined. For the search process the following three conditions had to be fulfilled: TD ≥ 1.9 s, sigbin ≥20% and at least two positive findings within the 14 ROIs.

The results of the search process in the four resting states are summarized in Table 2. Seven subjects complying with the above conditions were found in R1, three subjects in R2, two subjects in R3 and three subjects in R4. Altogether, 12 different subjects were identified with significant phase-coupling between RRI and BOLD oscillations in the pons/brain stem and a SSMTD = +2.2 ± 0.2 s. These findings show a clear predominance for R1. They primarily comprise the ROIs 93, 99 and 103 and depict slow RRI oscillations leading those in BOLD signals by about 2.2 s. Remarkably, this SSMTD was very consistent across these subjects despite of widely different breathing periods. In these selected subjects, the sigbin values varied between 23 and 82%. With the exception of two subjects (14R4, 3R2) all ROIs with dominant pTDs were left-sided. Noteworthy, the small inter-subject variability of the RRI vs. BOLD time delay (pTD = 2.2 ± 0.2 s) indicates an interesting “landmark.”

TABLE 2

Table 2. Age (years), state anxiety (AS), breath duration (Tresp [seconds], given as means ± SD), time delay (TD [seconds]) and significant phase coupling (%sig) for ROIs with a subject-specific maximal pTD (SSMTD) in pons/brain stem, ROI 1 (precentral gyrus) and ROI 93 in pons/brain stem.

For a comparison, pTD calculation was also performed for ROI 93, a region in the pons/brain stem with large sigbins, and ROI 1, representing the left precentral gyrus (Table 2). In both cases, the mean pTD was smaller and its inter-subject variability (SD) higher as compared to SSMTD.

Slow RRI Oscillations Lead BOLD Oscillations

RRI and BOLD signals have quite different origins. One is recorded from the thorax (RRI) and composed of successive cardiac beat-to-beat intervals measured in seconds. The other signal (BOLD) is recorded from the pons/brain stem and is composed of different neural and non-neural components (Obrig et al., 2000; Buxton et al., 2004; Murphy et al., 2013). It may be helpful to illustrate this novel approach by depicting the signals (also including respiration) in one characteristic subject (11R1), thus highlighting the dynamic and coincidence of various physiological signals (Figure 4). Subject 11R1 exhibits a 1:1 coupling between RRI and respiration signals. The graph shows a dominant coincidence of slow waves in RRI, respiration and BOLD signals from ROI 100 (right side). Furthermore, a delay of 2–3 s in BOLD waves from pons/brain stem (ROI 103, left side) and precentral gyrus (ROI 1) in comparison to RRI waves is shown. This delay of averaged BOLD waves (Figure 4, right panel) corresponds quite well with the pTD values in Table 2 (ROI 103: pTD = 2.5 s; ROI1: pTD = 2.4 s). The sigbin values (82 and 54%) underline the relative long-lasting stability of the phase coupling in this subject, although the dominant period of breathing waves (about 80% of all recorded breaths) varied between 7 and 10 s.

FIGURE 4

Figure 4. Examples of ongoing BOLD, RRI and respiration signals (1:1 coupling) and corresponding averaged waves (±SE) for one subject (11R1) with dominant 0.10-Hz and 0.15-Hz oscillations. The vertical dashed lines indicate maxima (peaks) of waves with a period of ∼7 s of RRI oscillations used as trigger for averaging. Peaks of the averaged waves are indicated. Note, the peak difference between RRI and BOLD (ROI 103) waves of about 2.6 s, the coincidence of RRI, respiration and BOLD (ROI 100) waves and the cessation of respiratory sinus arrhythmia (nonRSA). Modified from Pfurtscheller et al. (2019).

The examples shown in Figure 4 demonstrate that two types of spontaneous slow BOLD oscillations in the pons/brain stem can be observed, one of presumably neural origin and another of non-neural origin. The former is delayed (relative to the RRI signal) by 2–3 s in the case of 1:1 coupling (one breath per RRI cycle), the latter is in-phase with respiration (see also, Pfurtscheller et al., 2019). In subjects with a normal resting respiration rate (e.g., of 0.3 Hz and about three breaths per RRI cycle of 9 s), a similar pattern was found with a delay of 2–3 s between slow RRI and BOLD oscillations in pons/brain stem and precentral gyrus.

pTDs and State Anxiety

An interesting question is the functional meaning of increased pacemaker activity. To clarify this, the AS of the 23 subjects was investigated in relation to the pTD values dominating in both sides. Across resting states, the highest and lowest anxiety score was identified for each subject, thus resulting in two anxiety categories (high anxiety and low anxiety). In most cases, the highest AS score was found in the first resting state (Chapman et al., 2010; Pfurtscheller et al., 2018). However, there were 6 subjects with elevated AS in one of the subsequent resting states. The pTDs from both anxiety categories are summarized in Table 3 and visualized as histograms in Figure 5. Anxiety categories significantly differed in AS (p < 0.001), and there was a significant difference between the pTD values on the left side (p < 0.001). High anxiety was accompanied by a higher pTD as compared to low anxiety. The difference of pTDs between left and right brain side was also significant for high anxiety only (p = 0.006), suggesting a left side dominance. Additional analyses were conducted to compare differences in the distribution of TDs between high and low anxiety categories using Kolmogorov-Smirnov tests for two samples. Results revealed that pTDs from low and high anxiety category, respectively, seem to stem from different populations (Z = 2.40, p < 0.001). These findings may suggest that a pacemaker activity (large pTDs) may be associated with the processing of anxiety.

TABLE 3

Table 3. Anxiety category (AC), state anxiety (AS), positive time delays (pTD) in (s) for the left and right brain side for “high anxiety” (HA) and “low anxiety” (LA).

FIGURE 5

Figure 5. The upper row depicts histograms of TDs from the “high anxiety” category (23 subjects, 7 ROIs out of ROIs 91–106) for the left (A) and right side (B). The lower row shows TD histograms from the “low anxiety” category for left (C) and right side (D). ROIs with odd numbers represent the left, ROIs with even numbers the right brain side. Histograms representing pTDs (between 0.1 and 3 s) are indicated by dark gray color.

Discussion

Phase locking analysis between RRI and BOLD signals in the pons/brain stem revealed several interesting findings that support our pacemaker hypothesis. The findings that the time windows of the time delays (TDs) between RRI and BOLD never exceeded 3 s and that their maximal time delay (“landmark”) showed a quite precise mean ± SD of pTD = 2.2 s ± 0.2 s in 50% of the sample with a left-sided dominance in the first resting state are particularly noteworthy.

BOLD Signals in Brain Stem

Our interest is directed to the upper brain stem with the pons and the basilar artery as one of the major brain arteries directly ahead. The pons has a rostro-caudal length of ∼2.5 cm and a diameter of ∼ 3 cm in humans (Koehler et al., 1985), includes specific nuclei (e.g., the parabrachial nuclei and locus coeruleus, two structures pivotal for control of cardiovascular activity; Davern, 2014). Furthermore, nuclei of the reticular formation network originating in the lower brain stem seem particularly important. The AAL atlas (Tzourio-Mazoyer et al., 2002) allows to study BOLD signals in 14 ROIs (7 on each side) composed of many voxels in the upper brain stem with a spacing of 2 mm, of which ROIs 91and 92 represent the rostral part and ROIs 107 and 108 the caudal part of the pontine structure (see Figure 2). Studying the specific nuclei in pontine structure needs disentangling signals from noise carefully and cleaning-up of BOLD signals (Brooks et al., 2013).

Unspecific reticular neurons might be involved in the generation of slow rhythms in the pons/brain stem. Lambertz and Langhorst (1998) studied hundreds of reticular neurons in the lower brain stem (medulla oblongata) of adult dogs and observed spontaneous discharges at almost 0.15 Hz in more than 50%. Reticular neurons have also been found in the upper brain stem (pons). Projections from pontine reticular formation to the thalamus have been localized by diffuse tensor imaging (Yeo et al., 2013). All these reticular neurons are part of a diffuse ascending network in lower and upper brain stem and play an important role not only in maintaining behavioral arousal and consciousness (Moruzzi, 1958) but also in control of cardiovascular activity. Pontine structures are also activated during slow breathing (Critchley et al., 2015) and it is concluded that pontine neural activity is essential for physiological coupling of centrally generated respiratory and cardiovascular activities (Eckberg and Karemaker, 2009). The presence of Mayer waves in brain stem sympathetic-related and non-sympathetic-related cells was discussed by Ghali and Ghali (2020).

Prerequisites for the verification of a slow rhythmic neural source in the pons/brain stem are the availability of an associated slow hemodynamic oscillation (neural BOLD signal) delayed by the NBC and the existence of a type of reference signal. Such a reference signal can be either (in rare occasions such as during slow breathing) a non-neural BOLD signal time-locked to slow rhythmic motion of the basilar artery (vascular artifact) or the time-course of cardiac beat-to-beat intervals (RRI) linked to the neural activation and associated with vagal function. Examples of respiratory artifacts are accentuated in Figure 1 close to the cingulum (ROI 36) and in the pons/brain stem (ROIs 95 and 96). The cingulum is surrounded by large blood vessels of the anterior cerebral artery and shows high perfusion (Birn et al., 2006). An almost identical artifact is generated by the basilar artery in the pons/brain stem. A BOLD signal can be used as reference, however, only in the case of spontaneous slow breathing and cessation of RSA (also termed as “switch-off” of RSA, nonRSA; Pfurtscheller et al., 2019). Figure 4 depicts an example with coherent patterns of breathing and BOLD signals in the pons/brain stem (ROI 100). This respiratory artifact in the pons/brain stem can be used to monitor the cyclic start of vasomotion of the basilar artery. Notably, at the same moment of time a delayed neural BOLD signal can be observed in a short distance (e.g., ROI 103 in Figure 4). This example in Figure 4 highlights the existence of two different types of slow BOLD signals in close proximity in the upper brain stem in the case of slow spontaneous internally controlled breathing.

When rhythmic neural activation and control of cardiac function are closely linked, hemodynamic oscillations provoked by rhythmic neural activity changes are also coupled with fluctuations of the cardiac beat-to-beat intervals. Measuring of phase-coupling between RRI and BOLD signals in the 0.10–0.15 Hz band is independent of the respiratory rate and is not based on the selection of epochs with slow RRI fluctuations. Therefore, it can be applied to all BOLD data. Nevertheless, it is remarkable that such similar pTDs with a small variance (pTD = 2.2 s ± 0.2 s; N = 12) were observed in different pons/brain stem ROIs separated by few millimeters (see Table 2). Note, data of 5 subjects (part of our cohort) revealed a time delay between two BOLD signals in pons/brain stem of 2.6 ± 0.4 s (Pfurtscheller et al., 2019), meaning that two completely different approaches revealed the about same result.

The availability of various BOLD signals in a close proximity of a few millimeters in pons/brain stem suggests that the neural oscillatory source is not focused and not stable in a specific localization but has its origin very likely in a dispersive neural network (Ramirez et al., 2004) in the brain stem. Neural network dynamics are not fixed but depend on the state of the network or on the action of neuromodulators (Peña-Ortega, 2012), which is in accordance with our observations. The SSMTD of 2.2 ± 0.2 s documents that the time delay between RRI and almost undisturbed neural BOLD oscillations in the pons/brain stem is relatively rare. In the majority of cases, there exists a superposition of neural and non-neural waves, resulting in histogram peaks of pTDs between 1 and 2 s (Figure 3). This may be, at least in part, an explanation for the great inconsistency of time delays between RRI and BOLD signals.

Although a maximal pTD was documented in pons/brain stem, a time delay of similar size was observed in distant ROIs as for example in precentral gyrus. This is not unexpected considering that the ARAS is localized in the brain stem. The observation of similar intrinsic slow waves in ROIs representing pons/brain stem and precentral gyrus suggests a similar neural source. However, further studies are necessary to search for the origin of this neural source. Activation of cortical areas cannot completely be excluded during fMRI-induced anxiety and associated emotional breathing. Behavioral influences during uncomfortable supine position with the head in a limited noisy space may activate primary sensory areas in cortex prior to respiratory neurons in brain stem (Homma and Masaoka, 2008).

Hemodynamic Delay and pTD

It is important to distinguish between the hemodynamic response function (HRF) used for BOLD corrections in fMRI connectivity studies and the hemodynamic delay characteristic for NBC during rest. In the former case the participant is requested to focus attention to a stimulus of different modality. Therefore, the stimulus-evoked response is affected by stimulus processing, depends on stimulus modality and peaks at 4 – 5 s (West et al., 2019). In contrast, NBC indicates the relationship between slow spontaneous neural and hemodynamic (BOLD) oscillations. Hence, it could be expected that the NBC is shorter than HRF. Interestingly, short delays in humans were found by Fox et al. (2007) and, moreover, in mice (Mateo et al., 2017). Fox et al. (2007) found spontaneous BOLD signals of significantly different magnitude 2 to 3 s after soft vs. hard button pressing in humans. In an experimental trial with awake head-fixed adult mice (Mateo et al., 2017), the entrainment of arteriole vasomotor fluctuations at ∼ 0.1 Hz by neural activity was accompanied by BOLD oscillations about 2.6 s later.

The reported time delay of pTD = 2.2 ± 0.2 s indicates that RRI oscillations precede those in BOLD signals by 2–3 s or, in other words, the longest RRI is followed by a BOLD maximum about 2–3 s later. Inasmuch stimulation of the vagus nerve requires 0.2–0.3 s to produce HR slowing (Smith, 1974), a hemodynamic delay between neural and BOLD oscillations of almost 3 s is therefore realistic.

Although mice and humans differ considerably in body size, vasomotion (a natural collective oscillation of contractile tone in smooth muscles of arterioles) and neural fluctuations (beta/alpha power or envelope over gamma power) seem to be similar and occur at frequencies around 0.1 Hz both in mice (Drew et al., 2011; Mateo et al., 2017) and humans (Obrig et al., 2000; Pfurtscheller et al., 2012). Therefore, it is reasonable to assume that hemodynamic delay is similar, too.

For the investigation of neural-BOLD coupling, the in-vivo animal models have been widely used and essentially increased physiological knowledge (Shtoyerman et al., 2000). The background is that the topological matrices of specific animal brain networks are comparable with those of humans as shown by studies with various designs in rats (Liang et al., 2011; Ma et al., 2018), mice (Scholtens et al., 2018) and Macaque monkeys (Vincent et al., 2007). One of the pioneers in the study of neurovascular coupling by recording of ECoG and regional cerebral blood flow in anesthetized rats was Golanov et al. (1994). The time delay between the onsets of an ECoG burst and slow cerebrovascular wave varied between 1 and 3 s.

Time Delays Between Slow RRI and BOLD Oscillations Never Exceed ±3 s

The calculation of phase-coupling between physiological signals (RR-interval, pulse-interval, respiration) and fMRI signals (either voxel-based or ROI-based) is an established method (Lachaux et al., 2012; Pfurtscheller et al., 2017; Shokri-Kojori et al., 2018). Investigations of time delays in the 0.07–0.13 Hz band between RRI signals and ROIs in cortical and subcortical areas revealed delays between +3 s and −3 s (Pfurtscheller et al., 2017). Phase lag studies for low frequency fluctuations (<0.1 Hz) between pulse interval and respiratory signals, and voxels in cortical and subcortical areas exhibited mean phase lags between +1.21 ± 1.71 rad (mean ± SD) and −1.55 ± 1.07 rad (this corresponds to +1.9 s and −2.7 s for frequencies of 0.1 Hz) in 203 participants. These studies show that time delays between fMRI and physiological signals include a limited range of a few seconds.

Histograms of TDs are an appropriate way to visualize different features as for example the occurrence of dominant oscillations and specific limits. The histograms in Figure 3 are a catchy document of the prevalence of positive TDs with a peak between 1 and 2s and of the clear limitation within the range of ±3 s in pons/brain stem. This is notable because histograms of TDs from cortical and subcortical ROIs represent two clear peaks, one close to + 2s and another close to – 1s (not shown). While two peaks indicate the existence two dominant BOLD oscillations of neural and vascular origin, one peak of positive TDs refers to a dominance of neural BOLD oscillations interconnected with intrinsic rhythmic neural activity.

In order to visualize phase relationships between slow BOLD, RRI and neural oscillations, the phasor representation (Figure 6A) may be helpful (Zheng et al., 2010; Pfurtscheller et al., 2017). One cycle of 360° represents one 10-s period of a 0.1-Hz oscillation, and 90° phase-shift corresponds to TD of 2.5 s. Figure 6A shows the phasor representation of the relationship between neural, RRI and BOLD oscillations in 12 subjects, and Figure 6B depicts patterns of RRI signals from thorax and BOLD oscillations from pons/brain stem during rest. The strong coincidence between both oscillatory patterns with the delayed BOLD oscillations by 2–3 s points to a neural pacemaker in brain stem.

FIGURE 6

Figure 6. Phasor representation of neural, RRI and BOLD oscillations from 12 subjects with a time delay of 2.2 ± 0.2 s between RRI and BOLD oscillations (A). Neural oscillatory activity leads RRI oscillations by ∼0.3 s. B: Patterns of RRI recorded from thorax and BOLD oscillations in pons/brain stem from subject 11R1. Note, RRI oscillations recorded on thorax precede those in BOLD signals in left pons/brain stem (ROI 103) by ∼2.2 s (dashed lines).

Preponderance of pTDs in Left Pons/Brain Stem

The comparison of TD histograms across all resting states and subjects (Figure 3) revealed a significant side-difference (p < 0.001) with a dominance of pTDs on the left side. This asymmetry is largest in the first resting state and in accordance with the highest anxiety level (Table 1). The left-sided predominance of pTDs appears amplified in the “high anxiety” category and is accompanied by an elevated HR variability (HRV), and hence might provide resources for the processing of anxiety (Thayer and Lane, 2009). In contrast, a dominant activity of the sinoatrial node is associated with a reduced HRV and a diminished response to external autonomous stimulation (Moen et al., 2019).

The dominance of neural BOLD oscillations >0.1 Hz in the left side of the pons/brain stem is particularly noteworthy. Note, all subjects, except two, were right-handed. The vagus nerve plays an important role in fast cardiac control, and the parasympathetic (vagal) cardiac activity seems to be controlled by the left hemisphere (Craig, 2005). Amygdala studies also confirm the left hemisphere dominance (Baas et al., 2004).

Why Did Only 50% of the Subjects Show a Maximal pTD of 2.2 ± 0.2 s?

This question warrants some further thoughts. Note, Table 3 and Figure 5 provide evidence that a pacemaker activity (large pTDs) may be associated with the processing of anxiety during unpleasant settings, such as fMRI scanning. However, individuals may adopt different behavioral strategies to cope with such situations. For example, a breathing rhythm around 0.10 Hz has been reported to increase HRV and to facilitate emotional well-being (Lehrer, 2013; Mather and Thayer, 2018). Such a slow breathing pattern associated with a “switch-off” of RSA (Rassler et al., 2018) has also been observed in a minority of fMRI subjects. Therefore, it could be speculated that low levels of anxiety do not require additional HRV increase and thus no pacemaker activity. From this point of view, we may speculate that a portion of 50% of subjects hint toward a neural pacemaker operating as a kind of “emergency system” in threat-evoking situations, when the organism is in need of additional resources.

Limitatons and Future Prospect

This study has several limitations. A methodological limitation refers to artifacts in the ECG recording during scanning. When the scanning rate is high (more than one scan per second), artifacts might compromise the ECG recording (Kugel et al., 2003). Consequently, the detection of beat-to-beat-intervals is limited and hence, the calculation of phase-locking between BOLD and RRI signals becomes questionable. Therefore, BP beat-to-beat interval signals are also used instead of cardiac beat-to-beat intervals (Shokri-Kojori et al., 2018). However, it might be questioned whether BP and HR time series are equivalent. Note, the temporal resolution of the BOLD signal is limited to about 1 s (Buxton, 2013).

Another limitation of this protocol is that we did not measure end-tidal partial pressure of CO₂ and can’t exclude changes of vascular responses caused by episodes of hyperventilation.

In contrast to respiration, BP and BOLD signals, which represent real analog changes in the time domain, the cardiac interval time course (RRI signal) represents beat-to-beat intervals in seconds in relation to heart beats. The two extreme cases are that the actual RR interval is dedicated either to the start of the interval or to the end of the interval. The former case is known as instantaneous RRI signal and the latter case as delayed RRI signal (De Boer et al., 1985). Therefore, a time delay of about 1 s exists between both measurements. This means, in the case of a mean pTD = 2.2 s between RRI and BOLD signals, the actual time can be longer by fractions of a second.

The Mayer waves are composed of two components originating in the baroreflex loop and in an independent central oscillator, respectively (Julien, 2006). Both components belong to the low frequency band (0.04–0.15 Hz), but there is evidence that components below and above 0.10 Hz behave clearly differently (Schwerdtfeger et al., 2020). This was an important reason for studying the band 0.10–0.15 Hz. This particular frequency band is often disregarded due to the recommendation that BOLD recordings should preferably be analyzed after low-pass filtering at 0.10 or 0.09 Hz. This low-pass filtering suppresses all frequency components above 0.10 Hz and thus, reduces artifacts (Snyder and Raichle, 2012), but also precludes the assessment of rhythmic brain stem activity.

Our findings suggest that a central pacemaker is localized in pontine or related structures. Whether such a pacemaker activates cortical structures via the ARAS (Moruzzi, 1958) remains unclear. In order to investigate such a possibility, the Granger causality principle could be used to determine causal coupling between different BOLD and RRI oscillations (Blinowska et al., 2004). The observation of time delays of similar size between RRI oscillations and distant slow BOLD signals in precentral gyrus and pons/brain stem (see Table 2) points toward a similar neural source but this has to be studied in more detail.

Remarkably, in some subjects with increased state anxiety, inspiration was associated with a HR deceleration indicating a “switch-off” of RSA (nonRSA; Rassler et al., 2018). Whether this cessation of RSA is related to the activity of a central pacemaker or to the intrinsic activity of sinoatrial node cells (Moen et al., 2019) deserves further research.

Conclusion

(1) It could be shown that not only oscillations originating in the baroreflex loop contribute to the Mayer waves, but also central pacemaker oscillations. Such a pacemaker localized very likely in the pons/brain stem drives RRI oscillations and is associated with delayed BOLD oscillations due to the NBC time.

(2) RRI intervals are divided in LF (0.04 – 0.15 Hz) and HF (0.15 – 0.4 Hz) bands (e.g., Schwerdtfeger et al., 2020) with a boundary at 0.15 Hz. Note, that just this frequency, which characterizes the “0.15 Hz rhythm” of Perlitz et al. (2004) in brain stem was frequently observed in RRI and breathing signals (Rassler et al., 2018). Intrinsic frequency clusters, centered at 0.08 Hz and 0.15 Hz in BOLD signals, were reported by Yuen et al. (2019). Important is to note that physiological rhythms (and their pacemakers) operate over frequencies that are not controlled in rigid manner. While oscillations of BP and RRI operating in the LF band with origin in the baroreflex loop (Julien, 2006) are well established, the “0.15 Hz rhythm” in brain stem needs special attention. In our fMRI subjects, we observed two main slow frequencies in the RRI and BOLD oscillations, one with an average frequency of about 0.10 Hz, the other one with a frequency around 0.15 Hz. These oscillations represent the activity of the hypothesized pacemaker in the brain stem. In contrast, oscillations in the frequency band below 0.10 Hz were almost completely absent. Therefore, we defined the limits of our frequency band of interest by the averages of the two dominant slow rhythms. However, we plan to reassess the frequency limit at 0.15 Hz in the near future. Moreover, these two main slow frequencies document that the recommended lowpass filtering of BOLD signals at 0.10 Hz (Snyder and Raichle, 2012) may lead to a loss of information about slow oscillations that are of neural origin.

(3) The time delay between RRI and BOLD oscillations at ∼ 0.10 Hz with pTD = 2.2 ± 0.2 s found in 12 out of 23 subjects is particularly interesting. This finding strongly supports not only the existence of a central neural oscillator, but also documents that the neurovascular coupling time of 2–3 s in humans matches the delay between gamma power and BOLD signals of ^∼2.6 s in awake mice (Mateo et al., 2017). Thus, our study provides the first evidence that the hemodynamic delay between intrinsic neural and BOLD oscillations is of similar size in different species and ranges between 2 and 3 s in awake humans.

(4) The results suggest that a method is available to assess sources of neural oscillations in the pons/brain stem (central pacemaker) indirectly by analyzing phase coupling between cardiac beat-to-beat-intervals recorded from thorax and BOLD signals as a surrogate indicator.

Data Availability Statement

All datasets presented in this study are included in the article/supplementary material.

Ethics Statement

The studies involving human participants were reviewed and approved by Ethics Committee at the University of Graz. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

GP conceptualized and drafted the original manuscript. BR and AA contributed to the methodology, data processing, statistics, writing, and visualization. AS, BR, GS, and WK reviewed and edited the manuscript. All authors contributed to the article and approved the submitted version.

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.

Acknowledgments

We like to thank Thomas Zussner, David Fink, and Karl Koschutnig for support in data acquisition and Clemens Brunner for support in data processing. All University of Graz. We acknowledge support from the German Research Foundation (DPG) and Leipzig University within the program of Open Access Publishing.

Abbreviations

(f)MRI, (functional) magnetic resonance imaging; AAL, automated anatomy labeling; ARAS, ascending reticular activation system; AS, state anxiety; BOLD, blood-oxygenation-level-dependent; BP, blood pressure; ECG, electrocardiogram; EcoG, electrocorticogram; EEG, electroencephalogram; HR, heart rate; HRV, heart rate variability; NBC, neuro-BOLD coupling; nTD, negative time delay; PLV, phase-locking value; pTD, positive time delay; ROI, region of interest; RRI, beat-to-beat interval; RSA, respiratory sinus arrhythmia; SSMTD, subject-specific maximal pTD.

References

Arthurs, O. J., and Boniface, S. (2002). How well do we understand the neural origins of the fMRI bold signal? Trends Neurosci. 25, 27–31. doi: 10.1016/s0166-2236(00)01995-0

CrossRef Full Text | Google Scholar

Baas, D., Aleman, A., and Kahn, R. S. (2004). Lateralization of amygdala activation: a systematic review of functional neuroimaging studies. Brain Res. Rev. 45, 96–103. doi: 10.1016/j.brainresrev.2004.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Birn, R. M., Diamond, J. B., Smith, M. A., and Bandettini, P. A. (2006). Separating respiratory-variation-related fluctuations from neuronal-activity-related fluctuations in fMRI. Neuroimage 31, 1536–1548. doi: 10.1016/j.neuroimage.2006.02.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Blinowska, K. J., Kuś, R., and Kaminski, M. (2004). Granger causality and information flow in multivariate processes. Phys. Rev. E 70:050902. doi: 10.1103/PhysRevE.70.050902

PubMed Abstract | CrossRef Full Text | Google Scholar

Brooks, J. C. W., Faull, O. K., Pattingson, K. T. S., and Jenkinson, M. (2013). Physiological noise in brainstem fMRI. Front. Hum. Neurosci. 7:623. doi: 10.3389/fnhum.2013.00623

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

CrossRef Full Text | Google Scholar

Buxton, R. B., Uludag, K., Dubowitz, D. J., and Liu, T. T. (2004). Modeling the hemodynamic response to brain activation. Neuroimage 23, 220–233. doi: 10.1016/j.neuroimage.2004.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Chapman, H. A., Bernier, D., and Rusak, B. (2010). MRI-related anxiety levels change within and between repeated scanning sessions. Psych. Res. Neuroimaging 182, 160–164. doi: 10.1016/j.pscychresns.2010.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Craig, A. B. (2005). Forebrain emotional asymmetry: a neuroanatomical basis? Trends Cog. Sci. 9, 566–571. doi: 10.1016/j.tics.2005.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Critchley, H. D., Nicotra, A., Chiesa, P. A., Nahai, Y., Gray, A., Minati, L., et al. (2015). Slow breathing and hypoxic challenge: cardiorespiratory consequences and their central neural substrate. PLoS One 10:e0127082. doi: 10.1371/journal.pone.0127082

PubMed Abstract | CrossRef Full Text | Google Scholar

Davern, P. J. (2014). A role of the lateral parabrachial nucleus in cardiovascular function and fluid homeostasis. Front. Physiol. 5:436. doi: 10.3389/fphys.2014.00436

PubMed Abstract | CrossRef Full Text | Google Scholar

De Boer, R. W., Karemaker, J. M., and Strackee, J. (1985). Relationships between short-term blood-pressure fluctuations and heart-rate variability in resting subjects II: a simple model. Med. Biol. Eng. Comput. 23, 359–364. doi: 10.1007/bf02441590

PubMed Abstract | CrossRef Full Text | Google Scholar

Drew, P. J., Shih, A. Y., and Kleinfeld, D. (2011). Fluctuating and sensory-induced vasodynamics in rodent cortex extend arteriole capacity. Proc. Nat. Acad. Sci. U.S.A. 108, 8473–8478. doi: 10.1073/pnas.1100428108

PubMed Abstract | CrossRef Full Text | Google Scholar

Eckberg, D. L., Cooke, W. H., Diedrich, A., Biaggioni, I., Buckey, J. C. Jr., Pawelczyk, J. A., et al. (2016). Respiratory modulation of human autonomic function on Earth. J. Physiol. 594, 5611–5627. doi: 10.1113/JP271654

PubMed Abstract | CrossRef Full Text | Google Scholar

Eckberg, D. L., and Karemaker, J. M. (2009). Respiratory sinus arrhythmia is due to a central mechanism vs. respiratory sinus arrhythmia is due to the baroreflex mechanism. J. Appl. Physiol. 106, 1740–1744. doi: 10.1152/japplphysiol.91107.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Fox, M. D., Snyder, A. Z., Vincent, J. L., and Raichle, M. E. (2007). Intrinsic fluctuations within cortical systems account for intertrial variability in human behavior. Neuron 56, 171–184. doi: 10.1016/j.neuron.2007.08.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghali, M. G. Z., and Ghali, G. Z. (2020). Mechanisms contributing to the generation of mayer waves. Front. Neurosci. 14:395. doi: 10.3389/fnins.2020.00395

PubMed Abstract | CrossRef Full Text | Google Scholar

Golanov, E. V., Yamamoto, S., and Reis, D. (1994). Spontaneous waves of cerebral blood flow associated with pattern of electrocortical activity. Am. J. Physiol. 266, R204–R214. doi: 10.1152/ajpregu.1994.266.1.R204

PubMed Abstract | CrossRef Full Text | Google Scholar

Grinsted, A., Moore, J. C., and Jevrejeva, S. (2004). Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Process Geophys. 11, 561–566. doi: 10.5194/npg-11-561-2004

CrossRef Full Text | Google Scholar

Homma, I., and Masaoka, Y. (2008). Breathing rhythms and emotions. Exp. Physiol. 93, 1011–1021. doi: 10.1113/expphysiol.2008.042424

PubMed Abstract | CrossRef Full Text | Google Scholar

Julien, C. (2006). The enigma of Mayer waves: facts and models. Cardiovasc. Res. 70, 12–21. doi: 10.1016/j.cardiores.2005.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Koehler, P. R., Haughton, V. M., Daniels, D. L., Williams, A. L., Yetkin, Z., Charles, H. C., et al. (1985). MR measurements of normal and pathological brainstem dimeters. AJNR 6, 425–427.

Google Scholar

Kugel, H., Bremer, C., Peschel, M., Fischbach, R., Lenzen, H., Tombach, B., et al. (2003). Hazardous situation in the MR bore: induction in ECG leads causes fire. Eur. Radiol. 13, 690–694. doi: 10.1007/s00330-003-1841-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuusela, T. A., Kaila, T. J., and Kähönen, M. (2003). Fine structure of the low-frequency spectra of heart rate and blood pressure. BMC Physiol. 3:11. doi: 10.1186/1472-6793-3-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Lachaux, J.-P., Axmacher, N., Mormann, F., Halgren, E., and Crone, N. (2012). High-frequency neural activity and human cognition: past, present and possible future of intracranial EEG research. Prog. Neurobiol. 98, 279–301. doi: 10.1016/j.pneurobio.2012.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Lachaux, J., Rodriguez, E., Martinerie, J., and Varela, F. J. (1999). Measuring phase synchrony in brain signals. Hum. Brain Mapp. 208, 194–208. doi: 10.1002/(sici)1097-0193(1999)8:4<194::aid-hbm4>3.0.co;2-c

CrossRef Full Text | Google Scholar

Lambertz, M., and Langhorst, P. (1998). Simultaneous changes of rhythmic organization in brainstem neurons, respiration, cardiovascular system and EEG between 0.05 Hz and 0.5 Hz. J. Auton. Nerv. Syst. 68, 58–77. doi: 10.1016/s0165-1838(97)00126-4

CrossRef Full Text | Google Scholar

Laux, L., Hock, M., Bergner-Koether, R., Hodapp, V., Renner, K. H., and Merzbacher, G. (2013). Das State-Trait-Angst-Depressions-Inventar [The State-Trait Anxiety-Depression Inventory]. Goettingen: Hogrefe.

Google Scholar

Lehrer, P. (2013). How does heart rate variability biofeedback work? Resonance, the baroreflex, and other mechanisms. Biofeedback 41, 26–31. doi: 10.5298/1081-5937-41.1.02

CrossRef Full Text | Google Scholar

Liang, Z., King, J., and Zhang, N. (2011). Uncovering intrinsic connectional architecture of functional networks in awake rat brain. J. Neurosci. 31, 3776–3783. doi: 10.1523/jneurosci.4557-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Z., Ma, Y., and Zhang, N. (2018). Development of brain-wide connectivity architecture in awake rats. Neuroimage 176, 380–389. doi: 10.1016/j.neuroimage.2018.05.009

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, 1–13. doi: 10.1016/j.neuron.2017.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Mather, M., and Thayer, J. F. (2018). How heart rate variability affects emotion regulation brain networks. Curr. Opin. Behav. Sci. 19, 98–104. doi: 10.1016/j.cobeha.2017.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Moeller, S., Yacoub, E., Olman, C. A., Auerbach, E., Strupp, J., Harel, N., et al. (2010). Multiband multislice GE-EPI at 7 Tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI. Mag. Reson. Med. 63, 1144–1153. doi: 10.1002/mrm.22361

PubMed Abstract | CrossRef Full Text | Google Scholar

Moen, J. M., Matt, M. G., Ramirez, C., Tarasov, K. V., Chakir, K., Tarasova, Y. S., et al. (2019). Overexpression of a neuronal type adenylyl cyclase (Type 8) in sinoatrial node markedly impacts heart rate and rhythm. Front. Neurosci. 13:615. doi: 10.3389/fnins.2019.00615

PubMed Abstract | CrossRef Full Text | Google Scholar

Moruzzi, G. (1958). The functional significance of the reticular formation. Arch. Ital. Biol. 96, 17–28.

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

Niazy, R. K., Beckmann, C. F., Iannetti, G. D., Brady, J. M., and Smith, S. M. (2005). Removal of fMRI environment artifacts from EEG data using optimal basis sets. Neuroimage 28, 720–737. doi: 10.1016/j.neuroimage.2005.06.067

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Peña-Ortega, F. (2012). Pacemaker Neurons and Neuronal Networks in Health and Disease, Advances in Clinical Neurophysiology. Available online at: https://www.intechopen.com/books/advances-in-clinical-neurophysiology/pacemaker-neurons-and-neuronal-networks-in-health-and-disease (accessed August 19, 2020).

Google Scholar

Perlitz, V., Lambertz, M., Cotuk, B., Grebe, R., Vandenhouten, R., Flatten, G., et al. (2004). Cardiovascular rhythms in the 0.15-Hz band: common origin of identical phenomena in man and dog in the reticular formation of the brain stem? Eur. J. Physiol. 448, 579–591. doi: 10.1007/s00424-004-1291-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., Daly, I., Bauernfeind, G., and Müller-Putz, G. R. (2012). Coupling between intrinsic prefrontal HbO2 and central EEG beta power oscillations in the resting brain. PLoS One 7:e43640. doi: 10.1371/journal.pone.0043640

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., Rassler, B., Schwerdtfeger, A., Klimesch, W., Andrade, A., Schwarz, G., et al. (2019). “Switch-off” of respiratory sinus arrhythmia may be associated with the activation of an oscillatory source (pacemaker) in the brain stem. Front. Physiol. 10:939. doi: 10.3389/fphys.2019.00939

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., Schwerdtfeger, A., Seither-Preisler, A., Brunner, C., Aigner, C. S., Brito, J., et al. (2017). Brain-heart communication: evidence for “central pacemaker” oscillations with a dominant frequency at 0.1 Hz in the cingulum. Clin. Neurophysiol. 128, 183–193. doi: 10.1016/j.clinph.2016.10.097

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., Schwerdtfeger, A., Seither-Preisler, A., Brunner, C., Aigner, C. S., Calisto, J., et al. (2018). Synchronization of intrinsic 0.1-Hz blood-oxygen-level-dependent oscillations in amygdala and prefrontal cortex in subjects with increased state anxiety. Eur. J. Neurosci. 47, 417–426. doi: 10.1111/ejn.13845

PubMed Abstract | CrossRef Full Text | Google Scholar

Preiss, G., and Polosa, C. (1974). Patterns of sympathetic neuron activity associated with Mayer waves. Am. J. Physiol. 226, 724–730. doi: 10.1152/ajplegacy.1974.226.3.724

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramirez, J. M., Tryba, A. K., and Pena, F. (2004). Pacemaker neurons and neural networks: an integrative view. Curr. Opin. Neurobiol. 14, 665–674. doi: 10.1016/j.conb.2004.10.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Rassler, B., Schwerdtfeger, A., Aigner, C. S., and Pfurtscheller, G. (2018). Switch-off” of respiratory sinus arrhythmia can occur in a minority of subjects during functional magnetic resonance imaging (fMRI). Front. Physiol. 9:1688. doi: 10.3389/fphys.2018.01688

PubMed Abstract | CrossRef Full Text | Google Scholar

Scholtens, L. H., Feldman-Barrett, L., and Van den Heuvel, M. P. (2018). Cross-species evidence of interplay between neural connectivity at the micro- and macroscale of connectome organization in human, mouse, and rat brain. Brain Connect. 8, 595–603. doi: 10.1089/brain.2018.0622

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwerdtfeger, A. R., Schwarz, G., Pfurtscheller, K., Thayer, J. F., Jarczok, M. N., and Pfurtscheller, G. (2020). Heart rate variability (HRV): from brain death to resonance breathing at 6 breaths/min. Clin. Neurophysiol. 131, 676–693. doi: 10.1016/j.clinph.2019.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Shokri-Kojori, E., Tomasi, D., and Volkow, N. D. (2018). An autonomic network: synchrony between slow rhythms of pulse and brain resting state is associated with personality and emotions. Cereb. Cortex 28, 3356–3371. doi: 10.1093/cercor/bhy144

PubMed Abstract | CrossRef Full Text | Google Scholar

Shtoyerman, E., Arieli, A., Slovin, H., Vanzetta, I., and Grinvald, A. (2000). Long-term optical imaging and spectroscopy reveal mechanisms underlying the intrinsic signal and stability of cortical maps in V1 of behaving monkeys. J. Neurosci. 20, 8111–8121. doi: 10.1523/JNEUROSCI.20-21-08111.2000

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, O. A. (1974). Reflex and central mechanisms involved in the control of the heart and circulation. Ann. Rev. Physiol. 36, 93–123. doi: 10.1146/annurev.ph.36.030174.000521

PubMed Abstract | CrossRef Full Text | Google Scholar

Snyder, A. Z., and Raichle, M. E. (2012). A brief history of the resting state: the Washington University perspective. Neuroimage 62, 902–910. doi: 10.1016/j.neuroimage.2012.01.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Spielberger, C. D., Gorssuch, R. L., Lushene, P. R., Vagg, P. R., and Jacobs, G. (2009). Manual for the State-Trait Anxiety Inventory. Palo Alto, CA: Consulting Psychologists Press Inc.

Google Scholar

Tarvainen, M. P., Niskanen, J.-P., Lipponen, J. A., Rantaho, P. O., and Karjalainen, P. A. (2014). Kubios HRV – heart rate variability analysis software. Comput. Methods Programs Biomed. 113, 210–220. doi: 10.1016/j.cmpb.2013.07.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Thayer, J. F., and Lane, R. D. (2009). Claude Bernard and the heart-brain connection: further elaboration of a model of neurovisceral integration. Neurosci. Biobehav. Rev. 33, 81–88. doi: 10.1016/j.neubiorev.2008.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Torrence, C., and Compo, G. P. (1998). A practical guide to wavelet analysis. Bull. Am. Meteorol. Soc. 79, 61–78. doi: 10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2

CrossRef Full Text | Google Scholar

Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., et al. (2002). Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage 15, 273–289. doi: 10.1006/nimg.2001.0978

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Roon, A. M., Mulder, L. J. M., Althaus, M., and Mulder, G. (2004). Introducing a baroreflex model for studying cardiovascular effects of mental workload. Psychophysiology 41, 961–981. doi: 10.1111/j.1469-8986.2004.00251.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Verberne, A. J. M., and Owens, N. C. (1998). Cortical modulation of the cardiovascular system. Progr. Neurobiol. 54, 149–168. doi: 10.1016/s0301-0082(97)00056-7

CrossRef Full Text | Google Scholar

Vincent, J. L., Patel, G. H., Fox, M. D., Snyder, A. Z., Baker, J. T., Van Essen, D. C., et al. (2007). Intrinsic functional architecture in the anaesthetized monkey brain. Nature 447, 83–86. doi: 10.1038/nature05758

PubMed Abstract | CrossRef Full Text | Google Scholar

West, K. L., Zuppichini, M. D., Turner, M. P., Sivakolundu, D. K., Zhao, Y., Abdelkarim, D., et al. (2019). BOLD hemodynamic response function changes significantly with healthy aging. Neuroimage 188, 198–207. doi: 10.1016/j.neuroimage.2018.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Woletz, M., Hoffmann, A., Ganger, S., Paul, K., Seiger, R., Pfabigan, D. M., et al. (2014). “Slice-timing correction for multi-band images in SPM,” in Proceedings of the 20th Meeting of the Organization for Human Brain Mapping, Hamburg.

Google Scholar

Yan, C. G., and Zang, Y. F. (2010). DPARSF: a MATLAB toolbox for ‘pipeline’ data analysis of resting-state fMRI. Front. Syst. Neurosci. 4:13. doi: 10.3389/fnsys.2010.00013

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeo, S. S., Chang, P. H., and Jang, S. H. (2013). The ascending reticular activating system from pontine reticular formation to the thalamus in the human brain. Front. Hum. Neurosci. 7:416. doi: 10.3389/fnhum.2013.00416

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuen, N. H., Osachoff, N., and Chen, J. J. (2019). Intrinsic frequencies of the resting state fMRI signal: the frequency dependence of functional connectivity and the effect of mode mixing. Front. Neurosci. 13:900. doi: 10.3389/fnins.2019.00900

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, R., Zuckerman, J. H., Giller, C. A., and Levine, B. D. (1998). Transfer function analysis of dynamic cerebral auroregulation in humans. Am J Physiol. 274, H233–H241. doi: 10.1152/ajpheart.1998.274.1.H233

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, F., Sassaroli, A., and Fantini, S. (2010). Phasor representation of oxy- and deoxyhemoglobin concentrations: what is the meaning of out-of-phase oscillations as measured by near-infrared spectroscopy? J. Biomed. Opt. 15:040512. doi: 10.1117/1.3483466

CrossRef Full Text | Google Scholar