Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Neuroimage. Author manuscript; available in PMC 2009 August 1.
Published in final edited form as:
PMCID: PMC2603289

Brain Correlates of Autonomic Modulation: Combining Heart Rate Variability with fMRI


The central autonomic network (CAN) has been described in animal models but has been difficult to elucidate in humans. Potential confounds include physiological noise artifacts affecting brainstem neuroimaging data, and difficulty in deriving non-invasive continuous assessments of autonomic modulation. We have developed and implemented a new method which relates cardiac-gated fMRI timeseries with continuous-time heart rate variability (HRV) to estimate central autonomic processing. As many autonomic structures of interest are in brain regions strongly affected by cardiogenic pulsatility, we chose to cardiac-gate our fMRI acquisition to increase sensitivity. Cardiac-gating introduces T1-variability, which was corrected by transforming fMRI data to a fixed TR using a previously published method (Guimaraes et al. 1998). The electrocardiogram was analyzed with a novel point process adaptive-filter algorithm for computation of the high-frequency (HF) index, reflecting the time-varying dynamics of efferent cardiovagal modulation. Central command of cardiovagal outflow was inferred by using the HF timeseries resampled at as a regressor to the fMRI data. A grip task was used to perturb the autonomic nervous system. Our combined HRV-fMRI approach demonstrated HF correlation with fMRI activity in the hypothalamus, cerebellum, parabrachial nucleus/locus ceruleus, periaqueductal gray, amygdala, hippocampus, thalamus, and dorsomedial/dorsolateral prefrontal, posterior insular, and middle temporal cortices. While some regions consistent with central cardiovagal control in animal models gave corroborative evidence for our methodology, other mostly higher cortical or limbic-related brain regions may be unique to humans. Our approach should be optimized and applied to study the human brain correlates of autonomic modulation for various stimuli in both physiological and pathological states.

Keywords: HF, HRV, cardiovagal, exercise, handgrip


The autonomic nervous system (ANS) is comprised by the sympathetic, parasympathetic and enteric nervous systems (Langley 1921). While the physiology of the ANS has been well researched in animal models, the central processing and control of the ANS in the human have been more difficult to map and understand. This is because physiological artifacts degrade neuroimaging data, especially in brainstem structures. Furthermore, a non-invasive and continuous estimate of efferent sympathetic or parasympathetic activity has been elusive.

The peripheral afferent and efferent arms of the ANS do not constitute a simple feedback system, as feed-forward control by the brain can significantly affect autonomic outflow (Loewy et al. 1990). Afferent-efferent integration and central command of efferent activity is managed by the brain via the purported Central Autonomic Network (CAN) (Loewy et al. 1990; Benarroch 1993; Saper 2002; Janig 2006). Afferent (visceral) sensations are carried via cranial, sacral (parasympathetic), and thoracolumbar (sympathetic) nervous pathways to the nucleus tractus solitarius (NTS) in the medulla. This nucleus then transfers information to the parabrachial nucleus (PBN) in the pons, which then relays signals to the thalamus, hypothalamus, amygdala, and insula (sometimes considered the visceral primary sensorimotor cortex, see (Cechetto et al. 1990), and (Craig 2002). Several of these afferent processing regions (e.g. NTS and insula) can send information to efferent (premotor) ANS nuclei located in the medulla. These efferent regions include the rostral ventrolateral medulla (rVLM, sympathetic), dorsal motor nucleus of the vagus (DMNX, parasympathetic) and the nucleus ambiguus (NAmb, parasympathetic). Importantly, efferent regions can also be modulated by excitatory and inhibitory inputs (e.g. central command) from a set of regions including the rostral ventromedial medulla (rVMM), midbrain periaqueductal gray (PAG), hypothalamus, amygdala, and dorsomedial prefrontal and anterior cingulate cortices (Loewy et al. 1990; Benarroch 1993; Saper 2002; Janig 2006). While these regions have been well described in animal models, it is not clear if analogous and/or different connections exist in the human.

In specific regard to the heart, beat-to-beat dynamics are mediated via efferent outflow from the medulla. Bradycardic parasympathetic (cardiovagal) outflow to the sinoatrial (SA) node comes mainly from preganglionic neurons in the NAmb, with some input from the DMNX (McAllen et al. 1978; Janig 2006). Tachycardic response is mediated via sympathetic outflow from the rVLM (Dampney et al. 2002). Interestingly, these efferent outflows can be estimated non-invasively by analysis of heart rate variability (HRV) (Luczak et al. 1973; Sayers 1973; Akselrod et al. 1981; Malik et al. 1996). This method typically involves deriving the frequency spectrum of a continuous estimate of heart rate (HR). While some controversy in interpretation remains, the spectral peak in a low frequency band (LF, 0.01–0.15Hz) is thought to be influenced by both cardiovagal and sympathetic activity, while the peak in a high frequency band (HF, 0.15–0.40Hz) is influenced solely by cardiovagal activity (Malik et al. 1996). Cardiovagal and sympathetic activity are not necessarily reciprocal or antagonistic (Berntson et al. 1991; Paton et al. 2005) and cannot be inferred from heart rate changes alone. For example, tachycardia may be produced by either increased sympathetic activity, decreased parasympathetic activity, both, or even peripheral chemoreflexes.

Non-invasive neuroimaging can be applied to study the central ANS in human subjects. Several groups have attempted to relate estimates of efferent autonomic activity using HF power with concurrent brain response using Positron Emission Tomography (PET) (Lane et al. 2001; Gianaros et al. 2004) or functional Magnetic Resonance Imaging (fMRI) (Critchley et al. 2003; Napadow et al. 2005) during behavioral tasks that are thought to modulate the ANS. Correlating fMRI task response with out-of-scanner baseline HF (O'Connor et al. 2007) or HF response to the in-scanner task (Matthews et al. 2004) has also been performed. Unfortunately, neuroimaging in many central ANS structures (e.g. brainstem, hypothalamus) has been hampered by physiological motion artifact from cardiorespiratory sources. In addition, relating HRV to neuroimaging data is difficult given that HRV analyses using fast Fourier transform (FFT) approaches must assume stationarity over a finite temporal window (minutes), which is at least an order of magnitude greater than the temporal resolution of neuroimaging methods such as fMRI (several seconds). Furthermore, FFT methods are not causal in nature, i.e. the output does not depend solely on present and past inputs. Hence, data driven approaches to incorporating validated HRV regressors into fMRI general linear modeling have been elusive.

In our study, we have collected concurrent fMRI and ECG data during a dynamic exercise task chosen to perturb the ANS. We used cardiac-gated fMRI and improved group coregistration to enhance the fidelity of results in brainstem regions. We have also employed a novel adaptive recursive algorithm based on point process methods to compute instantaneous and causal estimates of HR and HRV in order to infer autonomic activity. Combining these two methodological advances allowed us to estimate the central command network for exercise-induced cardiovagal activity.


This study was performed with seven (7) healthy, right-handed (Edinburgh Inventory (Oldfield 1971)) subjects (4 male, age: 21–33years), and was approved by the Massachusetts General Hospital Subcommittee on Human Studies. Functional MRI and ECG data were collected simultaneously in order to assess brain correlates of autonomic modulation, perturbed by a dynamic grip task.

Dynamic Grip Task Protocol and Instrumentation

Subjects were asked to lie supine on the scanner bed holding an MRI-compatible grip dynamometer in their non-dominant (left) hand. This grip device, similar to other MRI-compatible grip devices (Dettmers et al. 1996), has been used in past studies in our Center to assess brain response to grip tasks (Cramer et al. 2002). The subjects were asked to perform a single maximum voluntary contraction (MVC) in order to set the target force level for the following fMRI scan. The voltage output from the force transducer was acquired at 200Hz by our custom data acquisition system (Labview 7.0, DAQCard 6024E, National Instruments, Austin, Texas). During fMRI scanning, subjects were instructed to look at a projected mirror reflection of a bar-scale where real-time grip force was visualized as percentage of their MVC and to follow instructions presented on the screen. In total, there were 15 1-minute ON blocks interspersed between 16 1-minute REST blocks. The first REST block was 1.5 minutes in order to allow for stabilization of the point-process HRV estimation algorithm (see below). During REST blocks, subjects saw the word “rest” next to the bar-scale. During ON blocks, subjects saw the word “squeeze” appear next to an arrow pointing to the 60% level for 2 seconds followed by the word “rest” for 1 second. This 2-second grip / 1-second relax paradigm was repeated for the duration of the 1-minute ON block. A 1-minute dynamic grip task was chosen to avoid the potential confound of excessive muscle metaboreflex-induced changes in HR, MAP and HRV (Victor et al. 1989; Laughlin 1999), while a 60% MVC effort was chosen for a robust modulation of both parasympathetic and sympathetic outflow (Kluess et al. 2000). It should be noted that relative to dynamic grip exercise, isometric grip exercise produces sustained increases in intramuscular pressure and increased vascular resistance leading to increased mean arterial pressure (MAP), as well as muscle metaboreflex mediated increases in muscle sympathetic nerve activity (MSNA). As we were most interested in localizing central command mediated autonomic activity, we chose to use a relatively short-term dynamic grip stimulus to mitigate the effects of non-central ANS modulation.

MRI and Cardiorespiratory Physiology Data acquisition

Subjects were outfitted with a nasal cannula to measure end tidal CO2, and four MRI-compatible chest ECG electrodes. Electrocardiogram (ECG) and respiratory data were collected at 200Hz (InVivo Magnitude CV, Invivo Research Inc., Orlando, Florida). These cardiorespiratory physiology output data were recorded on the same laptop acquisition system as the grip device data (Labview 7.0, DAQCard 6024E, National Instruments, Austin, Texas).

Functional and structural scans were acquired using a 3.0 Tesla Siemens Trio MRI System with an 8-channel head array coil equipped for echo planar imaging. Structural images were collected prior to functional imaging using a T1-weighted MPRAGE sequence (TR/TE = 2.73/3.19ms, flip angle = 7°, FOV = 256×256mm; slice thickness = 1.33 mm) in order to facilitate visualization.

BOLD functional imaging was performed using a gradient echo T2*-weighted pulse sequence (TE = 30msec, matrix = 64 × 64, FOV = 200mm, flip angle = 90°). Volume acquisition was gated to the subjects’ every other ECG R-wave; thus, TR varied with each acquisition (approximately 2sec). Cluster acquisition of 17 coronal slices was completed within 926ms (about a single RR-interval). Each slice was 3mm thick with 0.6mm gap (voxel size 3.13 × 3.13 × 3.6mm) and was slightly tilted to be parallel to the brainstem longitudinal axis (Figure 1). This was done in order to minimize through-slice brainstem motion. Motion correction was incorporated into the fMRI scanning sequence with Prospective Acquisition Correction (PACE, Siemens Medical Systems, Erlangen, Germany) – a critical addition due to increased opportunity for head motion during our long duration scan run.

Figure 1
Functional MRI slice localization and schematic of the combined HRV-fMRI methodology used to derive the brain correlates of cardiovagal activity.

Cardiac-gated acquisitions were performed to minimize physiological motion artifact due to pressure wave pulsatility in arteries within and around the brain. Brain parenchyma motion is mostly cranio-caudal and is characterized by rapid displacement in systole, followed by a slow diastolic recovery (Poncelet et al. 1992). Poncelet et al. used cardiac-gated echoplanar MR imaging to estimate brainstem velocities to be 2mm/sec. A different study demonstrated initial caudal and subsequent cranial brainstem displacement of 2–3mm (Maier et al. 1994). This cardiogenic motion is greater with caudal distance away from the cerebrum. In addition, cardiogenic arterial pulsatility also introduces motion artifact and loss of sensitivity in several brain regions (such as the pons, tectum, and portions of the ACC) that are close to major arteries (Dagli et al. 1999). These regions may be plagued by diminished SNR and spurious fMRI results if cardiac pulsatility is not accounted for. As many autonomic structures of interest are in brain regions strongly affected by cardiogenic pulsatility, we chose to cardiac-gate our fMRI acquisition in order reduce blurring and indeterminate voxel localization, thereby increasing fMRI sensitivity (Guimaraes et al. 1998).

HRV Analysis

Several time-frequency approaches have been used to study non-stationary heart rate spectral dynamics, each of which has their advantages and disadvantages. For instance, the pseudo Wigner-Ville transform can respond well to fast transients but is dogged by artifactual cross terms which can impact interpretation (Mainardi et al. 2002). Wavelet approaches can set temporal resolution to match different frequency components of interest, but are still limited by the existence of a temporal window which, similar to Wigner-Ville, violates causality (i.e. response depends on future data values as well as past) (Mainardi et al. 2002; Toledo et al. 2003). An interesting approach which attempted to incorporate continuous HRV estimates with fMRI data used band-pass filter methods (Critchley et al. 2003); however, simply filtering the R-R time series over different frequency bands does not produce an easily interpretable metric, and has not been validated to correspond to either sympathetic or parasympathetic activity. Our previously validated approach uses point process methods to obtain continuous estimates of the heartrate power spectrum (see below). It has excellent temporal resolution and preserves causality.

ECG data (acquired at 200Hz) were processed with the WFDB (WaveForm DataBase) Software Package (Goldberger et al. 2000) and MATLAB 6.5 (The Mathworks, Inc. Natick, MA). Data were automatically annotated with careful manual correction for QRS peak detection in order to form an R-R interval time series (less than 1% of R-wave annotations needed to be corrected). A novel adaptive recursive algorithm was applied to the R-R series to compute instantaneous estimates of heart rate and heart rate variability from electrocardiogram recordings of R-wave events. This approach is based on the point process methods for neural spike train data analysis (Brown et al. 2001; Barbieri et al. 2004; Eden et al. 2004) already used to develop both local likelihood (Barbieri et al. 2005) and adaptive (Barbieri et al. 2006) heart rate estimation algorithms. The stochastic structure in the R-R intervals is modeled as an inverse Gaussian renewal process. The inverse Gaussian probability density is derived directly from an elementary, physiologically-based integrate-and-fire model (Barbieri et al. 2005; Barbieri et al. 2006). The model also represents the dependence of the RR interval length on the recent history of parasympathetic and sympathetic inputs to the SA node by modeling the mean as a linear function of the last p RR intervals. This set of p coefficients allows for decomposition of the spectral power (HRV) into classic spectral components: Very Low frequency (VLF, 0 to 0.04 Hz), Low Frequency (LF, 0.04 to 0.15 Hz), and High Frequency (HF, 0.15 to 0.5 Hz) (Malik et al. 1996). The point process recursive algorithm is able to estimate the dynamics of the model parameters, and consequently the time-varying behavior of each spectral index, at any time resolution. This statistical model for deriving the HRV timeseries has been cross-validated with standard time-frequency domain approaches for HRV analysis (Barbieri et al. 2005). The dynamic response for our method was found to be an improvement on the more conventional RLS algorithm (Barbieri et al. 2006). Figure 1 shows a diagram of the algorithm and how it is incorporated into a GLM framework with the fMRI data.

We considered the instantaneous heart rate index computed from the point process model as our measure of heart rate (HR). Of the different instantaneous HRV indices computed by the point process technique, we selected the instantaneous high frequency power (HF) as our representative HRV metric. Heart rate HF power has garnered the greatest consensus for reflecting respiratory sinus arrhythmia (RSA) modulation by parasympathetic activity (Malik et al. 1996; Lane et al. 2001; Thayer et al. 2002). Our choice was further validated by spectral analysis of the respiratory signal, which confirmed general stationarity of the respiratory variance in the HF range, in particular with no significant difference between ON and OFF epochs (Students’ paired t-test; p=0.3162). Furthermore, the relatively short ON blocks (1-minute) required algorithm settings for fast adaptive estimation which increased uncertainty of the LF power estimates, and consequently of the LF/HF index as well. The original HR and HF series were estimated using model order p = 8 at a sampling of Δ = 5ms, low-pass filtered at 1 Hz to prevent aliasing, and were resampled at the fMRI timepoints.

In order to test whether the dynamic grip task modulated heart rate and HF power, we performed two separate repeated measures two-way analyses of variance (ANOVA) using an algorithm developed for MATLAB (Trujillo-Ortiz et al. 2004). Our analysis contained two independent factors, TASKSTATE (ON or OFF) and BLOCK (1 to 15). The partial p-value for the TASKSTATE factor was used to test if either HR or HF power were significantly modulated by the grip task. The mean percent signal change comparing HR and HF during the ON block to its preceding OFF block for each subject was calculated. A Student’s t-test was performed to test if the mean change for any individual (across blocks) was significantly different from zero.

Functional MRI Analysis

Structural MRI data were co-registered to the MNI152 brain template using the Automated Brainstem Co-registration (ABC) method (Napadow et al. 2006) in order to maximize sensitivity in important brainstem and subcortical brain regions. Functional MRI data were first motion corrected with AFNI (NIH) through an iterated, linearized, weighted least-squares method with Fourier interpolation (Cox 1996). Data runs were excluded if gross translational motion exceeded 3mm on any axis. Any residual motion was then removed by performing probabilistic independent component analysis (ICA) and removing components related to motion artifact (e.g. positive / negative fMRI response on opposing edges of the brain, IC timeseries spikes consistent with motion correction timeseries spikes) (Beckmann et al. 2004).

We performed cardiac-gated fMRI in order to reduce cardiogenic motion artifact. However, cardiac gating produces data with variable TR, and hence variable T1-weighting. If the temporal variability of T1-weighting is significant it will overwhelm any BOLD fMRI T2*-weighting. Data were corrected for T1 variability (due to variable TR) using the method detailed by Guimaraes et al. (Guimaraes et al. 1998). Their study clearly demonstrated improved SNR (i.e. mitigated T1-variability) in brainstem regions (the inferior colliculus for auditory stimuli in their case) compared to non-gated fMRI. Briefly, T1 was estimated by minimizing the variance in fitting an exponential T1 decay curve to fMRI signal strength within each voxel. Data were then corrected to a fixed TR value for each time point. This procedure uses the fMRI data itself to perform the correction, thereby avoiding spatial registration errors which would be incurred if T1 were estimated in a separate MRI scan. Minimal spatial smoothing (FWHM=3mm) was done on the fMRI data as some brainstem nuclei are on the order of image voxel size. The data were also high-pass filtered in the temporal domain (fhigh = 0.0083 Hz) to remove baseline signal drifts and cardiorespiratory artifacts aliased this lower frequency range. Corfield et al. found that high-pass filtering with this cutoff frequency yielded similar results to more intensive artifact mitigation procedures such as regressing out cardiac and respiratory phase timings in a general linear model (Corfield et al. 1999).

Statistical parametric mapping was completed via a generalized linear model (GLM) by using FMRI Expert Analysis Tool (FEAT, FSL). Time-series statistical analysis was carried out using FILM (FMRIB's Improved Linear Model) with local autocorrelation correction. The hemodynamic response function utilized in the GLM analysis was a prescribed gamma function (standard deviation = 3s, mean lag = 6s). The block paradigm was defined to be “1” during the 1-minute ON blocks and “0” during the 1-minute OFF blocks.

HRV – f MRI Combined Analysis

In order to derive the brain correlates of autonomic modulation via dynamic exercise, we incorporated an HRV-derived regressor into the fMRI GLM (Figure 1). Two regressors were used in the GLM: 1) time-resolved HF modulation re-sampled at the fMRI TR, and 2) the ON-OFF block design orthogonalized with respect to the HF regressor. The block design regressor was included in the GLM to account for brain response to the motor (exercise) task. Orthogonalization of the block design regressor relative to the HF regressor was performed to insure that all of the variance in the fMRI time series associated with HF modulation would be attributed to that regressor, while orthogonalization was not performed on the HF regressor to maintain the physiological interpretation of this regressor (i.e. cardiovagal modulation). While HF is not, strictly speaking, an index of parasympathetic or cardiovagal tone (Malik et al. 1993), under physiological conditions and with healthy adult subjects we do expect HF to represent an estimate of parasympathetic output to the heart.

Parameter estimates in the brainstem using the HF regressor were transformed to the ABC-adjusted standard space and imported to a group-level fixed-effects analysis with FSL. Note that as we used a fixed-effects model, our results may apply only to the individuals comprising our cohort. The MNI space was used for group results outside the brainstem. Resultant statistical parametric maps were corrected for multiple comparisons at a FDR of 0.05 (Genovese et al. 2002). Data were then modestly clustered with a minimum volume of 2 image voxels to allow sensitivity for activation within smaller brainstem and subcortical nuclei. Additionally, individual subject maps were threshold at p<0.005 and any activation/deactivation was noted for brain regions reported by the group map result.

Anatomical Localization

Anatomical localization for forebrain neuroimaging data was aided by a stereotaxic (Talairach coordinate system) atlas (Mai et al. 2004). Given that the Talairach or MNI coordinate system is not applicable for the brainstem, we adopted a previously validated approach (DaSilva et al. 2002) under supervision of one of the investigators (NM) of this study. Briefly, localization was performed using a landmark-based topographical methodology wherein the human brainstem has been subdivided into 28 distinct regions of interest or parcellation units (PUs). These PUs included the rostral dorsal medulla (site of DMNX, NTS, and NAmb) and rostral dorsal pons (site of PBN, and locus ceruleus, LC), which were of principal interest in the present study. Activity within the PUs was assessed for consistency with known brainstem nuclei (Paxinos et al. 1995).


We found that our dynamic exercise handgrip task was associated with increased heart rate for every individual tested. The block-averaged response was calculated from the average percent change for each ON block relative to the preceding OFF block. The average percent increase in HR ranged from 6.2% (subject 2) to 37.6% (subject 3). While all subjects demonstrated significant increase in HR, not all subjects demonstrated consistent change in HF power. HF power was consistently decreased in three subjects and more varied in others. The average percent change in HF power ranged from −48.1% (subject 7) to 54.5% (subject 3).

Group level analysis demonstrated that the grip task produced activation in both contralateral and ipsilateral primary sensorimotor cortex. Activation was also observed in typical motor system regions including the supplementary motor area, bilateral thalamus, and ipsilateral cerebellum. These results were consistent with other neuroimaging studies evaluating brain response to a grip task (Dai et al. 2001; Cramer et al. 2002).

The continuous estimate for HF was incorporated as a regressor into the GLM analysis of fMRI data. The group results demonstrated that brain regions with fMRI time courses correlating with HF power fluctuation included brainstem, mesencephalon, and diencephalon, as well as other subcortical and cortical brain regions (Figure 2 and Figure 3, Table 1). Specifically, a positive correlation was found in the left hypothalamus, left amygdala, right anterior hippocampus, and right DMPFC and dorsolateral prefrontal cortex (DLPFC). Conversely, a negative correlation was found in the right PBN/LC, left PAG, right posterior hippocampus, right mediodorsal thalamus, left caudate and right septal nuclei, left posterior insula, and right medial temporal gyrus.

Figure 2
Parasympathetic (HF) modulation correlated with fMRI signal in the hypothalamus. (A) Group map derived by using the HF regressor in the GLM demonstrated significant correlation with fMRI signal in the left hypothalamus. (B) HF power time series and fMRI ...
Figure 3
Brain correlates of parasympathetic (HF) modulation. Group map using the HF regressor in a GLM demonstrated significant correlation of fMRI signal in several brainstem, subcortical, and cortical regions. (PBN = parabrachial nucleus; LC = locus ceruleus; ...
Table 1
Neural Correlates of HF Variability


Our study has introduced a combined HRV-fMRI method which derives the central neural correlates of a continuous and causal estimate of efferent cardiovagal activity. The brain regions implicated are likely components of the network modulating central command of parasympathetic outflow to the heart. Our data analysis demonstrated that fMRI activity in several important CAN brain regions including the parabrachial nucleus/locus ceruleus, cerebellum, periaqueductal gray, hypothalamus, amygdala, and insular and dorsomedial prefrontal cortices demonstrated significant correlation with cardiovagal activity assessed by HF power. While these regions have been implicated by invasive animal studies, our non-invasive approach also implicated other brain regions which may be unique to the human: dorsolateral prefrontal cortex, mediodorsal thalamus, hippocampus, caudate, septal n. and middle temporal gyrus.

Our results should be interpreted from the context of the known relationship between HF power and cardiovagal physiology. Under physiological conditions and without a significant change in respiratory signal power within the HF power band (confirmed in our data), HF power of the HR signal is an expedient marker of parasympathetic modulation (Malik et al. 1996). The SA node responds very quickly to vagal (~150msec latency, steady state at 1–2sec), in contrast to sympathetic (~1–2sec latency, steady state at 30–60sec), influence (Spear et al. 1979; Levy 1997). Hence, continuous HF power is an excellent metric to use as a regressor within the fMRI GLM framework, given the temporal resolution of fMRI.

Some of the brain regions where fMRI signal was correlated with HF were components of the CAN as implicated by more invasive animal studies, serving as corroborative evidence for our methodology. These regions included several brainstem and cerebellar regions. Sensitivity in these regions was enhanced by cardiac-gated fMRI and brainstem-specific coregistration techniques. The PBN / LC in the right pons were found to be negatively correlated to the HF regressor. The PBN is an important viscerosensory nucleus, relaying information from the NTS to a wide array of target regions including the PAG, hypothalamus, amygdala, and other forebrain nuclei (Gauriau et al. 2002). The LC is the source of noradrenergic input in the brain and is associated with arousal and the hypothalamic-pituitary-axis (HPA) stress response, relaying afferent signaling to the hypothalamus and amygdala (Benarroch 1993; Ziegler et al. 1999). Thus, a negative correlation to HF in the LC would be consistent with decreased stress response during increased cardiovagal modulation. The right cerebellum was also negatively correlated with HF. The cerebellum has been implicated in autonomic modulation arising from motor (e.g. grip task) and vestibular activity (Balaban 1999) via direct connections with the hypothalamus (Dietrichs et al. 1994). The PAG, which was also negatively correlated with HF, has been found by intracranial recording to respond to exercise with spectral power increase across multiple frequency bands (Green et al. 2007). It is directly connected with both sensory (NTS) and premotor (DMNX, NAmb) cardiovagal nuclei (Farkas et al. 1997), as well as sympathetic premotor nuclei.

FMRI signal in the hypothalamus was found to be positively correlated with HF. This brain structure is considered the principal controller of autonomic and neuroendocrine homeostasis in the body and is the organizational center for many behavioral autonomic response patterns, including exercise-induced ANS outflow (Loewy et al. 1990). FMRI signal in the amygdala was also found to be positively correlated with HF. Behavioral-associated (e.g. emotion) input from the amygdala to the hypothalamus is transmitted via the ventral amygdalofugal pathway (Parent 1996) and translates affect into sympathetic and/or parasympathetic response (Thayer et al. 2006). In addition, fMRI signal in several cortical areas was also either positively or negatively correlated with HF power. In monkeys, autonomic modulation by the DMPFC may result from interconnections with both the PAG (An et al. 1998) and hypothalamus (Ongur et al. 1998). The posterior insula is a known ANS sensorimotor region (Saper 2002). This region is interconnected with the hypothalamus and its stimulation in the rat can induce either a pressor or depressor response, depending on the particular site of stimulation (Yasui et al. 1991).

Due to phylogenetic divergence, higher cortical regions implicated by our method in humans are the most likely to be different from CAN regions evaluated in lower mammalian models. The DLPFC, positively correlated with HF, has not been strongly implicated in studies using animal models and may be unique to the more developed human neocortex. Furthermore, fMRI signal in the mediodorsal thalamus, hippocampus, caudate, septal n. and middle temporal gyrus correlated negatively with HF power. These regions have also not received extensive attention from animal studies of the CAN, but are strongly associated with limbic and/or higher cortical processing and may play a greater role for humans compared to lower mammalian models. For instance, the DLPFC is an important brain region governing human cognition (see Duncan et al. 2000 for review) and connects mono- and multi-synaptically with paralimbic CAN regions (Petrides 2005). Hence, the DLPFC may influence outflow commensurate to perceived grip task effort in a similar manner as was found for the ACC in a SPECT study which used hypnotism to suggest a perceived grip task (Williamson et al. 2002).

Several other studies have correlated or covaried cardiac metrics with neuroimaging data. Gianaros et al. covaried HF with H215O PET during graded memory tasks (Gianaros et al. 2004) and found that regional cerebral blood flow (rCBF) covaried positively with HF in the right VMPFC, left insula, and left amygdala-hippocampal complex, while covarying negatively in the right cerebellum. Lane et al. similarly covaried HF with PET neuroimaging during an emotion-laden task (Lane et al. 2001). Their results showed positive covariation in the left DMPFC, and left anterior insula/orbitofrontal cortex. A study by Critchley et al. used a bandpass filtered RR-series as a regressor to an fMRI GLM during both cognitive and isometric grip tasks (Critchley et al. 2003). The results found positive correlation in the cerebellum and several cortical regions including the supplementary motor area, cingulate and DMPFC. Our methodology benefited from this last approach, but incorporated important methodological advances such as improved brainstem sensitivity for fMRI and, perhaps most importantly, an accepted measure of continuous cardiovagal modulation (Malik et al. 1996), as opposed to forming regressors by bandpass filtering the RR time-series in the LF or HF band – an approach which has not been cross-validated or correlated with cardiovagal activity. These methodological differences may have resulted in unique HF-correlated activity in the hypothalamus, PAG by our method, while the cingulate, implicated by Critchley et al. and Matthews et al. (Matthews et al. 2004) and anatomically connected with the insula and amygdala (implicated by our approach), could not be assessed by our methodology due to partial brain coverage. However, when the methodology reported by Critchley et al. was applied to our own data, no brain regions passed an identical post-processing threshold cutoff used for our methodology. Other studies have attempted to correlate HR or MAP with fMRI (Kimmerly et al. 2005) or PET (Critchley et al. 2000). However, measures such as mean heart rate (or mean RR interval) do not proportionally quantify either sympathetic or parasympathetic modulation (Malik et al. 1996), and are thus less likely to relate to activity in discrete brain regions with known sympathetic or parasympathetic connections.

The most direct interpretation of brain regions derived by our method is that they are involved with central command of cardiovagal modulation – either premotor or via multi-synaptic higher cognitive control. However, other indirect possibilities exist and should be discussed. As our method is correlational, if afferent vagal signaling or efferent sympathetic signaling is correlated with efferent cardiovagal modulation, then the brain regions implicated by our method may also reflect viscerosensory processing. However, sympathetic outflow is not strictly antagonistic to vagal outflow (Janig 2006), while the afferent/efferent cardiovagal feedback loop is even less consistent and remains to be more fully described (Fallen 2005). Thus, the more likely interpretation of our results is that the regions implicated are indeed related to central command of what HF is a measure of – efferent parasympathetic modulation to the heart.

Several limitations in our study design should be discussed. Firstly, as we used cardiac-gated fMRI, our field of view could not cover the entire brain. As our hypotheses were focused on brainstem and subcortical autonomic regions, we chose to exclude the most anterior and posterior portions of the brain. We preferred to compromise brain coverage for improved fidelity in the brain regions we did image. Secondly, we made an assumption that brain regions correlating with HF power were somehow involved with central cardiovagal activity. We feel this assumption was justified given that humoral factors (renin-angiotensin, adrenomedullary catecholamines, muscle metaboreflex etc.) mainly influence lower frequencies, while the influence of respiratory-induced mechanical stretching of the SA node should be minimal in our subject population (Berntson et al. 1991). In addition, artifactual “activity” was found in the lateral ventricles. While our methods were focused on mitigating sources of physiological artifact in the brainstem, other sources, such as CSF pulsatility, may play a greater role in regions next to the lateral ventricles, diminishing our confidence for results in these regions (septal nucleus, mdThalamus). Finally, our derived brain correlates may be specific to cardiovagal modulation by exercise and may not be generalizable to other autonomic challenges (e.g. cognitive tasks) nor to autonomic outflow to other effector systems (e.g. sudomotor, pupillary, etc.).

In conclusion, we have described a combined HRV-fMRI approach to derive the neural correlates of exercise-induced cardiovagal outflow in humans. Our methodology incorporated a continuous estimate of HF power that better matched the temporal resolution of fMRI leading to a more accurate assessment of central command for cardiovagal outflow. Furthermore, cardiac gating and improved brainstem coregistration were used to improve the fidelity of measurements in brainstem and hypothalamic regions. This approach should be optimized and applied to study the human brain correlates of ANS modulation for various stimuli in physiological and pathological states.


We would like to thank the NIH for funding support (VN: K01-AT002166, P01-AT002048; RB: R01-HL084502; EB: R01-DA015644), as well the NCRR (P41RR14075), and the Mental Illness and Neuroscience Discovery (MIND) Institute. We would also like to thank Dr. Judith Schaechter for providing the grip device used in this study.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


  • Akselrod S, Gordon D, Ubel FA, Shannon DC, Berger AC, Cohen RJ. Power spectrum analysis of heart rate fluctuation: a quantitative probe of beat-to-beat cardiovascular control. Science. 1981;213(4504):220–222. [PubMed]
  • An X, Bandler R, Ongur D, Price JL. Prefrontal cortical projections to longitudinal columns in the midbrain periaqueductal gray in macaque monkeys. J Comp Neurol. 1998;401(4):455–479. [PubMed]
  • Balaban CD. Vestibular autonomic regulation (including motion sickness and the mechanism of vomiting) Curr Opin Neurol. 1999;12(1):29–33. [PubMed]
  • Barbieri R, Brown EN. Analysis of heartbeat dynamics by point process adaptive filtering. IEEE Trans Biomed Eng. 2006;53(1):4–12. [PubMed]
  • Barbieri R, Frank LM, Nguyen DP, Quirk MC, Solo V, Wilson MA, Brown EN. Dynamic analyses of information encoding in neural ensembles. Neural Comput. 2004;16(2):277–307. [PubMed]
  • Barbieri R, Matten EC, Alabi AA, Brown EN. A point-process model of human heartbeat intervals: new definitions of heart rate and heart rate variability. Am J Physiol Heart Circ Physiol. 2005;288(1):H424–H435. [PubMed]
  • Beckmann CF, Smith SM. Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans Med Imaging. 2004;23(2):137–152. [PubMed]
  • Benarroch EE. The central autonomic network: functional organization, dysfunction, and perspective. Mayo Clin Proc. 1993;68(10):988–1001. [PubMed]
  • Berntson GG, Cacioppo JT, Quigley KS. Autonomic determinism: the modes of autonomic control, the doctrine of autonomic space, and the laws of autonomic constraint. Psychol Rev. 1991;98(4):459–487. [PubMed]
  • Brown EN, Nguyen DP, Frank LM, Wilson MA, Solo V. An analysis of neural receptive field plasticity by point process adaptive filtering. Proc Natl Acad Sci U S A. 2001;98(21):12261–12266. [PubMed]
  • Cechetto DF, Chen SJ. Subcortical sites mediating sympathetic responses from insular cortex in rats. Am J Physiol. 1990;258(1 Pt 2):R245–R255. [PubMed]
  • Corfield DR, Murphy K, Josephs O, Fink GR, Frackowiak RS, Guz A, Adams L, Turner R. Cortical and subcortical control of tongue movement in humans: a functional neuroimaging study using fMRI. J Appl Physiol. 1999;86(5):1468–1477. [PubMed]
  • Cox RW. AFNI: software for analysis and visualization of functional magnetic resonance neuroimages. Comput Biomed Res. 1996;29(3):162–173. [PubMed]
  • Craig AD. How do you feel? Interoception: the sense of the physiological condition of the body. Nat Rev Neurosci. 2002;3(8):655–666. [PubMed]
  • Cramer SC, Weisskoff RM, Schaechter JD, Nelles G, Foley M, Finklestein SP, Rosen BR. Motor cortex activation is related to force of squeezing. Hum Brain Mapp. 2002;16(4):197–205. [PubMed]
  • Critchley HD, Corfield DR, Chandler MP, Mathias CJ, Dolan RJ. Cerebral correlates of autonomic cardiovascular arousal: a functional neuroimaging investigation in humans. J Physiol. 2000;523(Pt 1):259–270. [PubMed]
  • Critchley HD, Mathias CJ, et al. Human cingulate cortex and autonomic control: converging neuroimaging and clinical evidence. Brain. 2003;126(Pt 10):2139–2152. [PubMed]
  • Dagli MS, Ingeholm JE, Haxby JV. Localization of cardiac-induced signal change in fMRI. Neuroimage. 1999;9(4):407–415. [PubMed]
  • Dai TH, Liu JZ, Sahgal V, Brown RW, Yue GH. Relationship between muscle output and functional MRI-measured brain activation. Exp Brain Res. 2001;140(3):290–300. [PubMed]
  • Dampney RA, Coleman MJ, et al. Central mechanisms underlying short- and long-term regulation of the cardiovascular system. Clin Exp Pharmacol Physiol. 2002;29(4):261–268. [PubMed]
  • DaSilva AF, Becerra L, Makris N, Strassman AM, Gonzalez RG, Geatrakis N, Borsook D. Somatotopic activation in the human trigeminal pain pathway. J Neurosci. 2002;22(18):8183–8192. [PubMed]
  • Dettmers C, Connelly A, Stephan KM, Turner R, Friston KJ, Frackowiak RS, Gadian DG. Quantitative comparison of functional magnetic resonance imaging with positron emission tomography using a force-related paradigm. Neuroimage. 1996;4(3 Pt 1):201–209. [PubMed]
  • Dietrichs E, Haines DE, Roste GK, Roste LS. Hypothalamocerebellar and cerebellohypothalamic projections--circuits for regulating nonsomatic cerebellar activity? Histol Histopathol. 1994;9(3):603–614. [PubMed]
  • Duncan J, Owen AM. Common regions of the human frontal lobe recruited by diverse cognitive demands. Trends Neurosci. 2000;23(10):475–483. [PubMed]
  • Eden UT, Frank LM, Barbieri R, Solo V, Brown EN. Dynamic analysis of neural encoding by point process adaptive filtering. Neural Comput. 2004;16(5):971–998. [PubMed]
  • Fallen EL. Vagal afferent stimulation as a cardioprotective strategy? Introducing the concept. Ann Noninvasive Electrocardiol. 2005;10(4):441–446. [PubMed]
  • Farkas E, Jansen AS, Loewy AD. Periaqueductal gray matter projection to vagal preganglionic neurons and the nucleus tractus solitarius. Brain Res. 1997;764(1–2):257–261. [PubMed]
  • Gauriau C, Bernard JF. Pain pathways and parabrachial circuits in the rat. Exp Physiol. 2002;87(2):251–258. [PubMed]
  • Genovese CR, Lazar NA, Nichols T. Thresholding of statistical maps in functional neuroimaging using the false discovery rate. Neuroimage. 2002;15(4):870–878. [PubMed]
  • Gianaros PJ, Van Der Veen FM, Jennings JR. Regional cerebral blood flow correlates with heart period and high-frequency heart period variability during working-memory tasks: Implications for the cortical and subcortical regulation of cardiac autonomic activity. Psychophysiology. 2004;41(4):521–530. [PubMed]
  • Goldberger AL, Amaral LA, et al. PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals. Circulation. 2000;101(23):E215–E220. [PubMed]
  • Green AL, Wang S, et al. Identifying cardiorespiratory neurocircuitry involved in central command during exercise in humans. J Physiol. 2007;578(Pt 2):605–612. [PubMed]
  • Guimaraes AR, Melcher JR, et al. Imaging subcortical auditory activity in humans. Hum Brain Mapp. 1998;6(1):33–41. [PMC free article] [PubMed]
  • Janig W. The Integrative Action of the Autonomic Nervous System. Cambridge: Cambridge University Press; 2006.
  • Kimmerly DS, O'Leary DD, Menon RS, Gati JS, Shoemaker JK. Cortical regions associated with autonomic cardiovascular regulation during lower body negative pressure in humans. J Physiol. 2005;569(Pt 1):331–345. [PubMed]
  • Kluess HA, Wood RH, Welsch MA. Vagal modulation of the heart and central hemodynamics during handgrip exercise. Am J Physiol Heart Circ Physiol. 2000;278(5):H1648–H1652. [PubMed]
  • Lane R, Reiman E, Ahern G, Thayer J. Activity in medial prefrontal cortex correlates with vagal component of heart rate variability during emotion. Brain and Cognition. 2001;47:97–100.
  • Langley J. The Autonomic Nervous System. Part I. Cambridge: W.Heffer; 1921.
  • Laughlin MH. Cardiovascular response to exercise. Am J Physiol. 1999;277(6 Pt 2):S244–S259. [PubMed]
  • Levy MN. Neural control of cardiac function. Baillieres Clin Neurol. 1997;6(2):227–244. [PubMed]
  • Loewy A, Spyer K. Central Regulation of Autonomic Functions. New York: Oxford University Press; 1990.
  • Luczak H, Laurig W. An analysis of heart rate variability. Ergonomics. 1973;16(1):85–97. [PubMed]
  • Mai J, Assheuer J, Paxinos G. Atlas of the Human Brain. San Diego: Elsevier Academic Press; 2004.
  • Maier SE, Hardy CJ, Jolesz FA. Brain and cerebrospinal fluid motion: real-time quantification with M-mode MR imaging. Radiology. 1994;193(2):477–483. [PubMed]
  • Mainardi LT, Bianchi AM, Cerutti S. Time-frequency and time-varying analysis for assessing the dynamic responses of cardiovascular control. Crit Rev Biomed Eng. 2002;30(1–3):175–217. [PubMed]
  • Malik M. A. Camm and and members of the Task Force of the European Society of Cardiology. Heart rate variability: standards of measurement, physiological interpretation and clinical use. Circulation. 1996;93(5):1043–1065. [PubMed]
  • Malik M, Camm AJ. Components of heart rate variability--what they really mean and what we really measure. Am J Cardiol. 1993;72(11):821–822. [PubMed]
  • Matthews SC, Paulus MP, Simmons AN, Nelesen RA, Dimsdale JE. Functional subdivisions within anterior cingulate cortex and their relationship to autonomic nervous system function. Neuroimage. 2004;22(3):1151–1156. [PubMed]
  • McAllen RM, Spyer KM. Two types of vagal preganglionic motoneurones projecting to the heart and lungs. J Physiol. 1978;282:353–364. [PubMed]
  • Napadow V, Dhond R, Kennedy D, Hui KK, Makris N. Automated brainstem co-registration (ABC) for MRI. Neuroimage. 2006;32(3):1113–1119. [PubMed]
  • Napadow V, Dhond R, Purdon P, Kettner N, Makris N, Kwong K, Hui K. Correlating Acupuncture fMRI in the Human Brainstem with Heart Rate Variability. 27th Annual International IEEE EMBS Conference; Shanghai, China. 2005. [PubMed]
  • O'Connor MF, Gundel H, McRae K, Lane RD. Baseline vagal tone predicts BOLD response during elicitation of grief. Neuropsychopharmacology. 2007;32(10):2184–2189. [PubMed]
  • Oldfield RC. The assessment and analysis of handedness: the Edinburgh inventory. Neuropsychologia. 1971;9(1):97–113. [PubMed]
  • Ongur D, An X, Price JL. Prefrontal cortical projections to the hypothalamus in macaque monkeys. J Comp Neurol. 1998;401(4):480–505. [PubMed]
  • Parent A. Carpenter's Human Neuroanatomy. Baltimore: Williams & Wilkins; 1996.
  • Paton JF, Boscan P, Pickering AE, Nalivaiko E. The yin and yang of cardiac autonomic control: vago-sympathetic interactions revisited. Brain Res Brain Res Rev. 2005;49(3):555–565. [PubMed]
  • Paxinos G, Huang X-F. Atlas of the Human Brainstem. San Diego: Academic Press; 1995.
  • Petrides M. Lateral prefrontal cortex: architectonic and functional organization. Philos Trans R Soc Lond B Biol Sci. 2005;360(1456):781–795. [PMC free article] [PubMed]
  • Poncelet BP, Wedeen VJ, Weisskoff RM, Cohen MS. Brain parenchyma motion: measurement with cine echo-planar MR imaging. Radiology. 1992;185(3):645–651. [PubMed]
  • Saper CB. The central autonomic nervous system: conscious visceral perception and autonomic pattern generation. Annu Rev Neurosci. 2002;25:433–469. [PubMed]
  • Sayers BM. Analysis of heart rate variability. Ergonomics. 1973;16(1):17–32. [PubMed]
  • Spear JF, Kronhaus KD, Moore EN, Kline RP. The effect of brief vagal stimulation on the isolated rabbit sinus node. Circ Res. 1979;44(1):75–88. [PubMed]
  • Thayer JF, Siegle GJ. Neurovisceral integration in cardiac and emotional regulation. IEEE Eng Med Biol Mag. 2002;21(4):24–29. [PubMed]
  • Thayer JF, Sternberg E. Beyond heart rate variability: vagal regulation of allostatic systems. Ann N Y Acad Sci. 2006;1088:361–372. [PubMed]
  • Toledo E, Gurevitz O, Hod H, Eldar M, Akselrod S. Wavelet analysis of instantaneous heart rate: a study of autonomic control during thrombolysis. Am J Physiol Regul Integr Comp Physiol. 2003;284(4):R1079–R1091. [PubMed]
  • Trujillo-Ortiz A, Hernandez-Walls R, Trujillo-Perez RA. RMAOV2:Two-way repeated measures ANOVA. 2004.
  • Victor RG, Seals DR. Reflex stimulation of sympathetic outflow during rhythmic exercise in humans. Am J Physiol. 1989;257(6 Pt 2):H2017–H2024. [PubMed]
  • Williamson JW, McColl R, Mathews D, Mitchell JH, Raven PB, Morgan WP. Brain activation by central command during actual and imagined handgrip under hypnosis. J Appl Physiol. 2002;92(3):1317–1324. [PubMed]
  • Yasui Y, Breder CD, Saper CB, Cechetto DF. Autonomic responses and efferent pathways from the insular cortex in the rat. J Comp Neurol. 1991;303(3):355–374. [PubMed]
  • Ziegler DR, Cass WA, Herman JP. Excitatory influence of the locus coeruleus in hypothalamic-pituitary-adrenocortical axis responses to stress. J Neuroendocrinol. 1999;11(5):361–369. [PubMed]