## Abstract

The causal, directed interactions between brain regions at rest (brain–brain networks) and between resting-state brain activity and autonomic nervous system (ANS) outflow (brain–heart links) have not been completely elucidated. We collected 7 T resting-state functional magnetic resonance imaging (fMRI) data with simultaneous respiration and heartbeat recordings in nine healthy volunteers to investigate (i) the causal interactions between cortical and subcortical brain regions at rest and (ii) the causal interactions between resting-state brain activity and the ANS as quantified through a probabilistic, point-process-based heartbeat model which generates dynamical estimates for sympathetic and parasympathetic activity as well as sympathovagal balance. Given the high amount of information shared between brain-derived signals, we compared the results of traditional bivariate Granger causality (GC) with a globally conditioned approach which evaluated the additional influence of each brain region on the causal target while factoring out effects concomitantly mediated by other brain regions. The bivariate approach resulted in a large number of possibly spurious causal brain–brain links, while, using the globally conditioned approach, we demonstrated the existence of significant selective causal links between cortical/subcortical brain regions and sympathetic and parasympathetic modulation as well as sympathovagal balance. In particular, we demonstrated a causal role of the amygdala, hypothalamus, brainstem and, among others, medial, middle and superior frontal gyri, superior temporal pole, paracentral lobule and cerebellar regions in modulating the so-called central autonomic network (CAN). In summary, we show that, provided proper conditioning is employed to eliminate spurious causalities, ultra-high-field functional imaging coupled with physiological signal acquisition and GC analysis is able to quantify directed brain–brain and brain–heart interactions reflecting central modulation of ANS outflow.

## 1. Introduction and state of the art

### (a) Imaging the autonomic nervous system

The autonomic nervous system (ANS) owes its name to John Newport Langley [1], a British physiologist widely regarded as the father of autonomic neurophysiology. As later communicated by Langley himself, ‘it was a local autonomy that [he] had in mind. The word “Autonomic” does suggest a much greater degree of independence from the central nervous system that in fact exists’ [2]. Hence, since its discovery and despite its nomenclature, the ANS has never been considered or treated as completely independent from the central nervous system (CNS). More than a century after its discovery and initial description, there is still substantial debate about the nature as well as functioning of the control pathways through which the CNS interacts with the ANS.

While functional magnetic resonance imaging (fMRI) based on the blood-oxygen-level dependent (BOLD) effect [3] allows access to non-invasive, spatially and temporally resolved functional information about brain activity at rest as well as during task execution, direct non-invasive *in vivo* recordings of ANS activity are not accessible. In this context, the analysis of heart rate variability (HRV) has emerged as a valid signal-processing approach to the estimation of ANS outflow, which also allows us to partially disentangle sympathetic and parasympathetic activity. In particular, band-specific spectral powers in the RR interval time series (i.e. the series of intervals between two consecutive ‘R’ peaks in the electrocardiogram (ECG)) are thought to provide a reliable static estimate of the parasympathetic (high frequency (HF) 0.15–0.40 Hz) activity, as well as mixed information about sympathetic and parasympathetic tones (low frequency (LF) 0.04–0.15 Hz) [4]. While both time- and frequency-domain RR analysis have been widely employed to assess the strength of these two efferent ANS pathways, only a few techniques (e.g. the pseudo Wigner–Ville transform [5] or wavelet decompositions [6,7]) have allowed a time-resolved estimation of ANS activity, which is mandatory to allow joint analysis with brain functional imaging. In a pioneering paper, Critchley and co-workers [8] jointly analysed the fMRI BOLD signal and HRV indices calculated from concurrently acquired ECG signals, demonstrating a relationship between the dorsal anterior cingulate cortex and sympathetic modulation. Furthermore, using a sliding window approach to RR series analysis, Chang *et al.* [9] presented the first joint fMRI/HRV analysis on a resting-state dataset, investigating the role of the amygdala and dorsal anterior cingulate cortex (which are involved in the central mediation of vigilance and emotional arousal) in modulating undirected brain–brain resting-state functional connectivity.

In order to improve the reliability of time-resolved ANS activity estimation, we have proposed a novel recursive algorithm to estimate instantaneous heart rate (HR) and HRV through a probabilistic point-process (PP) framework [10–13]. Contrary to classical time–frequency decomposition approaches, this framework does not require interpolation of RR series, is inspired by a physiological integrate-and-fire model and provides an instantaneous probabilistic estimate of the occurrence of the next R peak at arbitrary time-resolution from which a number of linear and higher-order HRV-related estimates can be obtained [3,14–16]. The dynamical estimates of ANS activity provided by this framework have been employed in conjunction with BOLD fMRI to obtain *in vivo*, non-invasive information about the functional CNS correlates of the ANS. In this context, using a dynamic grip task Napadow *et al.* [17] demonstrated that fMRI BOLD activity in several cortical and subcortical brain regions exhibits significant correlation with cardiovagal activity, as assessed by the HF component of HRV. In a similar set-up, Sclocco *et al.* [18] demonstrated that the activity of a number of brain regions, such as the insula, regions belonging to the default mode network and those belonging to visual motion areas, were strongly associated with both sympathetic and parasympathetic activity during nauseogenic visual stimulation.

This body of literature has helped to characterize a so-called central autonomic network (CAN) (see for instance [19]) whose existence was first postulated in several animal models [20,21]. In this context, in a large meta-analysis of recent neuroimaging studies, Thayer *et al.* [22] showed that cerebral blood flow as measured with fMRI or positron emission tomography (PET) is associated with HRV and identified a number of brain areas whose association with HRV is modulated by emotional or cognitive/motor tasks. In another recent meta-analysis, Beissner *et al.* [23] studied the CAN and its constituents using an activation-likelihood framework. They showed that the CAN involves some core constituents such as the amygdala, right anterior and left posterior insula and midcingulate cortices in a task-specific manner. In addition, the CAN also includes several other brain regions as components of the ANS regulatory network, such as the insula, ventromedial prefrontal cortex, hippocampal formation, mediodorsal thalamus, hypothalamus, posterior cingulate cortex, lateral temporal cortices and bilateral dorsal anterior insula.

To date, all studies investigating the functional links between the CNS and ANS have employed concepts of undirected associations (e.g. correlations) and did not employ ultra-high-field technology. The directed, causal influences within the brain regions comprising the CAN as well as between the CAN and HRV-related metrics of ANS activity have not yet been investigated. Given that CAN regions interact tightly among themselves, the study of bivariate, undirected correlations may not be sufficient (or may introduce spurious results) when aiming to disentangle the central pathways of ANS modulation.

### (b) Causality in the brain

Arguably the most influential contribution to the definition of causality was proposed by Clive W. J. Granger [24]. In brief, a dynamical system is considered to be Granger-causing (G-causing) if observations from the past of that dynamical system allow better predictions of the future of another dynamical system, compared with predictions based on the past of the latter alone. This general definition allowed detection of time-domain causality between time series with a certain degree of confidence [25]. Following the seminal work of Geweke [26], who proposed a measure of directed feedback between time series, the idea of Granger causality (GC) has been mostly associated with its implementation in the context of autoregressive (AR) models [27] and its expansion to multivariate systems through the use of vector AR (VAR) models [28]. Additionally, when GC is implemented in terms of AR and VAR models, it is also possible to model the influence of a variable over another one by including instantaneous (zero temporal lag) influences [15,29,30].

Along with the classical time-domain approach to GC, a number of spectral methods for estimating causality have been proposed in the recent literature. In one of the earliest developments of spectral GC, directed coherence (DC) was used to characterize brain causality between two electroencephalogram (EEG) signals [31]. The concept was further generalized to the multivariate case by Baccalá *et al.* [32], through the idea of multivariate DC, and in the seminal paper by Kaminski & Blinowska [33], where the concept of directed transfer function (DTF) was introduced. Many other spectral methods for GC estimation have been developed, including the partial directed coherence (PDC) [34], generalized PDC (gPDC) [35], information PDC (iPDC) and information DTF (iDTF) [36]. Baccalá *et al.* [37] presented a unified theoretical framework that encompasses all aforementioned spectral domain techniques. Interestingly, the theoretical development of spectral methods for studying causality has been fuelled by its potential application to the analysis of EEG signals, due to the robust body of biomedical literature relating band-specific EEG powers with distinct physiological functions and states [38].

### (c) Causality in functional magnetic resonance imaging

Goebel *et al.* [39] and Harrison *et al.* [40] pioneered the use of VAR models for the analysis of fMRI data, and GC techniques are now routinely employed in fMRI analysis [41,42] (see also the review by Seth *et al.* [43]). While fMRI offers exceptional local specificity in quantifying dynamical brain activity (i.e. a spatial resolution or 3 mm or less, resulting in 10^{4}–10^{5} signals for every fMRI acquisition), compared with EEG data fMRI signals are typically collected with low sampling frequency (usually 1 Hz or less), short data lengths (of the order of minutes or tens of minutes) and low signal-to-noise ratio (SNR) [44], possibly limiting the applicability of VAR-based GC methods. A valid workaround is to obtain locally averaged BOLD time series according to predefined brain atlases (see [45] for a review). Another very recent approach to overcome the dimensionality problem is to employ principal component analysis (PCA) or independent components analysis (ICA) prior to fitting the VAR models. The feasibility of the ICA-based reduction in the context of time-varying BOLD activity has been demonstrated in [46] and has been used to introduce the concept of large-scale GC (lsGC). Also, in the context of EEG analysis ICA allowed the estimation of large-scale, time-varying PDC [47,48]. A proof-of-concept application of GC after PCA in fMRI has been performed in [49], and in [50] lsGC was used to characterize network modules.

A potentially limiting factor in applications of GC in fMRI is the very nature of the BOLD signal, which stems from a convolution product between the neuronal activation and an unknown haemodynamic response function (HRF) reflecting the local homeostatic adjustment of blood flow at the tissue level [51]. Characterization of the relationship between neuronal activation and the BOLD signal is far from complete (see [52] and [53] for reviews), and poses practical difficulties due to inter-subject variability [54], is tissue-dependent both in human [55] and in animal models [56], and may also be time-dependent [57] as well as task-dependent [58] (although inter-subject variability seems to dominate all other sources of fluctuations [59]). In this context, some authors have suggested addressing the HRF-related confound through a blind data-driven deconvolution method under the assumption that neural activation is sparse and of a binary nature [60]. This method has been employed in resting-state fMRI to characterize transfer entropy between PCA components (note that transfer entropy is equivalent to GC under the assumption of Gaussian noise [61]) in computing GC through a lagged version of structural equation modelling [62]. While the assumptions behind blind deconvolution may appear arbitrary, in general they have proved useful in repeated tasks paradigms and when the neural activation can be described by a purely threshold-trigger mechanism [63]. In an effort to investigate the true influence of the HRF-related confounds, recent papers employed direct intracranial measurements of neural activity to explicitly investigate the effect of HRF convolution on GC estimates *in vivo* [64], at the population level [46] and with synthetic data simulations [65] using analytically solvable VAR models as well as VAR models with non-trivial topologies and realistic haemodynamics. They showed that GC is ‘extremely resilient’ against HRF convolution even when local variations in HRF parameters strongly degrade information about time precedence between neural signals. Instead, GC estimates are sensitive to signal down-sampling in the time domain and to low SNR. With respect to the problem of time resolution (which may be the most important limiting factor in GC applications of fMRI) Hemmelmann *et al.* [66] employed both simulated (through the balloon model [67]) and real BOLD signals to demonstrate that combining the complementary information provided by time-domain GC and frequency-domain PDC allows correct detection of GC in spite of HRF-related confounds. It should be noted that, when using causality techniques to analyse the directed interaction between brain-derived BOLD time series and external signals, an alternative approach is the convolution of the external signal with the best assumed HRF.

### (d) Conditioning of causality

While the problem of causality detection in the brain is clearly multivariate, in line with most classical, correlation-based connectivity studies the use of bivariate frameworks dominates the recent literature on directed intra-brain interactions. However, when quantifying brain–brain interactions, the high inter-dependence of brain-derived signals can lead to the appearance of false positives (i.e. spurious, unrealistic causal links). In this context, the importance of ‘conditioning’ GC between two brain regions (i.e. estimating their ‘true’ directed interaction while accounting for indirect influences mediated by additional brain regions) has been proposed by Guo *et al.* [68], and these also demonstrate the analogy between conditioned GC and the ideas behind partial correlation. Also Marinazzo *et al.* [69] defined a partial conditioning approach based on *a priori*, mutual information-based identification of the subset of signals which share the most information with the target signal under study, and demonstrated its application to the brain. While in [69] a prior downselection within conditioning variables was deemed necessary due to the limited length of typical fMRI time series, in GC analyses of the brain the widest possible set of conditioning variables (i.e. brain regions) should be included within a well-defined multivariate conditioning framework [28,41–43,65]. This is equally relevant when studying the causal links between multiple brain regions and external signals, like time-varying estimates of ANS activity.

### (e) Aims of the study

In this paper, we provide the first causality-based investigation of directed brain–heart interactions in the resting state and associate our results with estimates of directed, conditioned brain–brain resting-state networks. In order to mitigate fMRI confounds related to sensitivity, temporal resolution and physiological noise, we employed state-of-the-art, ultra-high-field fMRI combined with simultaneous multi-slice (SMS) echo-planar imaging (EPI) techniques and physiological signals acquisition. We employed time-varying estimates of heartbeat dynamics and ANS outflow obtained through our state-of-the-art probabilistic PP framework and investigated how proper conditioning of GC estimates, which accounts for the intra-subject interdependence between local, fMRI-derived brain activity across the whole brain, affected inference about directed brain–brain and directed brain–heart interactions.

## 2. Theory

Estimating GC from the *j*th variable to the *i*th variable (*j*→*i*) amounts to testing the null hypothesis that knowledge of past values of *j* improves the prediction of (or restricts the uncertainty about) the future of *i*. In AR-model-based bivariate GC, we consider two models. First, the so-called restricted (order *P*)-AR model of the *i*th variable, which includes only the past values of *i* itself. Second, the so-called unrestricted AR model, which includes both variables *i* and *j*. The unrestricted model has order *P* with respect to the *i*th variable, and order *Q* with respect to the *j*th variable,
2.1
and
2.2
where *a*′_{p} and *a*_{k,p} are the AR model parameters, and *ε*_{t} and *ε*′_{t} are uncorrelated white noise processes. The restricted model represents the null hypothesis (*Q*=0) that should be rejected if *j* significantly G-causes *i*. The unrestricted model is the alternative hypothesis that *j* significantly G-causes *i*. Prior to model identification, optimal AR orders *P* and *Q* have to be inferred. As far as we are aware, all applications of GC to brain signals assume identical AR orders *P*=*Q*. However, given that *P* can only be chosen based on merit, quantified through some sort of information criterion (e.g. final prediction error of the Akaike [70], the Bayesian information criterion (BIC) [71] or the Pagano and Hartley *t*-test [72]) of the restricted model (and should not be varied in the unrestricted model), *Q* should be chosen independently while adopting the same figure of merit which has been adopted for the restricted model. Information criteria, rules of thumb and extensive searches are all strategies that have been employed in the past. One simple approach entails conducting a search of all possible *Q*>0 (up to a reasonable limit) to be employed as the alternative hypothesis, and Geweke [25] suggests that the number of lags *Q* of the independent variable should be kept smaller than the number of lags *P* of the dependent variable.

We now extend the bivariate GC model (equations (2.1) and (2.2)) to the so-called ‘conditioned’ GC approach. If the causing variable *x*_{j} is highly correlated with a set of variables (*x*_{1}…*x*_{L}) in order to correctly disentangle (*j*→*i*) in both the restricted and unrestricted model the influence of variables (*x*_{1}…*x*_{L}) should be included (i.e. the causality estimate is ‘conditioned’ to exclude the influence of the latter variables); the restricted and unrestricted conditioned models now read
2.3
and
2.4
As before, given *P* and *R* as the optimal regression orders for the restricted model, the order *Q* is to be chosen independently.

In order to accept or reject the null hypothesis, in either approach (bivariate or conditioned), the residual sum of squares (RSS) is used to compute the so-called *f*_{ratio} which will follow an *F*-distribution with *Q* and (*N*_{obs}−*P*−*LR*−*Q*) degrees of freedom [73,74],
2.5
where *N*_{obs} is the number is observations of *x*_{i}, and RSS_{r} and *RSS*_{ur} refer to the restricted and unrestricted models, respectively. We can then reject the null hypothesis and assume that a significant GC relation (*j*→*i*) exists if the observed *f*_{ratio} corresponds to a *p*-value defined as
2.6
which is less than 0.05 [74]. This means that the reduction in the RSS resulting from the inclusion of *j* in the prediction of *i* is significantly higher than would be expected from including a non-G-causing variable in the prediction of *i*.

## 3. Material and methods

### (a) Ultra-high-field magnetic resonance imaging with physiological signal acquisition

Nine healthy volunteers (age 28±3 years) underwent 7 T MRI with simultaneous physiological signal acquisition under institutional review board approval at the Athinoula A. Martinos Center for Biomedical Imaging (MGH) on a Siemens 7 T whole-body scanner (Siemens Healthcare, Erlangen, Germany). We used a single-shot 2D EPI readout for 1.8 mm isotropic T2*-weighted axial images, with matrix size/generalized autocalibrating partially parallel acquisitions (GRAPPA) factor/nominal echo-spacing=136×136/2/0.57 ms and additional parameters echo time/flip angle/repetitions=26 ms/40°/300. The use of an SMS factor equal to 3 enabled whole-brain coverage (number of slices=75) and a repetition time (TR) of 1.5 s which minimized aliasing in ANS frequency bands (see below). With TR=1.5 s, the sampling frequency of the brain signal was 0.667 Hz. Physiological noise correction (high-pass filtering at 0.01 Hz and removal of second-order RETROspective Image CORrection [75] regressors) was applied using custom MATLAB scripts, while slice timing, motion correction and co-registration to Montreal Neurological Institute (MNI) space were applied to fMRI data using the Oxford Centre for Functional MRI of the Brain (FMRIB) Software Library (FSL). Successively, the average BOLD signal was extracted in 116 regions of interest (ROIs) using the automated anatomical labeling atlas [76] and a 117th region relative to the brainstem was extracted using the Harvard–Oxford subcortical atlas (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases). Cardiac pulsation was recorded by a piezoelectric finger pulse sensor and respiration by a piezoelectric respiratory bellow positioned around the subject's chest (1 kHz sampling).

### (b) Point-process-based model for heart rate variability

The cardiac pulsation signal was used to detect cardiac peaks from which a peak-to-peak interval series (henceforth called the RR series) was derived and modelled as a PP [11], hence circumventing the need for RR interpolation. Briefly, assuming history dependence, the probability distribution of the waiting time *t*−*u*_{j} until the next R-wave event *u*_{i−1} follows an inverse Gaussian model:
3.1
with the instantaneous mean RR interval defined as:
3.2
where is a left continuous function denoting the index of the previous R-wave event that occurred before time *t*, *H*_{t}=(*u*_{j},RR_{j},RR_{j−1},…) is the history of events, *Δ*RR_{i}=(RR_{N(t)−i}−RR_{N(t)−i−1}), *ξ*(*t*)=[*ξ*_{0}(*t*),*γ*_{0}(*t*),…*γ*_{1}(*i*,*t*)…] is the parameter vector, and *ε*(*t*) are independent, identically distributed Gaussian random variables. Parameter estimation was carried out using a local maximum-likelihood method within a sliding window of duration *W* and instantaneous, time-varying indices of HF power [77] (0.15–0.4 Hz, thought to reflect mainly parasympathetic activity), LF power (0.04–0.15 Hz, thought to reflect both sympathetic and parasympathetic activity) and sympathovagal balance (Bal=LF/HF) were estimated from the first-order regression terms. In this paper the parameters of the PP were evaluated at a frequency of 1000 Hz, and then downsampled to match the frequency of the BOLD signal. In this paper, *W* was varied in 20 regular intervals within a previously validated plausible range *W*=[30 s,120 s], after which for each regressor (HF, LF, Bal) the median value across parameter sets was computed at each time point to minimize bias in model parameter selection. The resulting ANS regressors were thresholded under the 95th percentile in order to enhance sensitivity to the full dynamics of the HRV time series and to ensure robustness to outliers. Finally, time-varying indices of HF, LF and Bal were convolved with a standard double gamma HRF (*α*_{1}=6, *α*_{2}=16, *β*_{1}=1, *β*_{2}=2, *c*=0.17) [78]. Example signals generated by this procedure are shown in figure 1.

### (c) Brain–brain Granger causality

We computed both bivariate and globally (i.e. to all other brain regions) conditioned within-brain (i.e. brain–brain) causality between locally averaged brain signals, using equations (2.1) and (2.2) and equations (2.3) and (2.4), respectively, and .

For the bivariate GC model, the AR orders were *P*_{b}=5 and *P*_{c}=5 for both the causing variable and the caused variable, corresponding to the minimum of the BIC [71] after independent optimizations as described above. For the globally conditioned GC (GCGC) model
and
the model orders for both the causing variable and the caused variable were *P*_{b}=5 and *P*_{c}=5 and *L*=117 regions. The optimal model order (corresponding to the minimum BIC value) for all conditioning variables was *R*=1—it is worth noting that, using 300 data points per subject, the maximum *R* we could have employed would have been *R*=2. The *f*_{ratio} followed an *F*-distribution with *P*_{c} and (*N*_{obs}−*P*_{b}−*P*_{c}+117) degrees of freedom.

### (d) Brain–heart Granger causality

As above, we computed both bivariate and globally (i.e. to all other brain regions) conditioned brain–heart causality between locally averaged brain signals and the PP-derived dynamical band-specific powers (henceforth called regressors). To estimate bivariate GC between a brain region *x*_{i} and a heart regressor *H*_{sig} (where ‘sig’ labels HF, LF and Bal) an AR model of order *P*_{b} is adopted for the brain signal, and an AR model of order *P*_{h} is adopted for the heart regressor. The AR order for the heart regressor was *P*_{h}=6, corresponding to the minimum of the BIC the restricted model. The AR order for the brain signal was *P*_{b}=5, corresponding to the minimum of the BIC for the unrestricted model given the order of the BIC for the restricted model. When evaluating the GCGC (*x*_{j}→*H*_{sig}) between a brain region *x*_{j} and a heart regressor *H*_{sig} both the restricted and unrestricted model (equations (2.3) and (2.4)) included the influence of all other brain variables, i.e. the restricted and unrestricted models are
3.3
and
3.4
where *C*_{k,p}, *C*′_{k,p}, *D*_{k,p}, *D*′_{k,p}, are VAR model parameters and *η*_{t} and *η*′_{t} are uncorrelated white noise processes. In the GCGC the AR model for the heart regressor was *P*_{h}=6. The AR order *P*_{m} for conditioning brain variables was set to *P*_{m}=1, corresponding to the minimum value of the BIC with respect to *H*_{sig,ur}(*t*). The AR model order *P*_{q} for causing variable *x*_{j} was set to *P*_{g}=3 to minimize the BIC with respect to *H*_{sig,ur}(*t*).

### (e) Group analysis

For each subject, both bivariate GC and GCGC between each of the 117 brain regions and every other region were computed. Additionally, bivariate GC and GCGC between each of the 117 brain regions and each heart regressor (HF, LF and Bal) were computed. For each directed connection between any two brain regions, for each connection between each brain region and every ANS regressor, and for both bivariate GC and GCGC the appropriate *f*_{ratio} was used to evaluate the *p*-value of the connection. Successively, to provide a group-wise count for each connection the number of subjects with *p*<0.01 and *p*<0.05 was counted. A measure of relative connection strength was evaluated as the logarithm of *f*_{ratio},
3.5
and used for visualization purposes. All processing was performed through custom-build code in Mathematica 10.

## 4. Results

### (a) Within-brain networks

Figures 2 and 3 show the bivariate GC strengths (see equation (3.5)) of all connections in all nine subjects separately as well as the within-subject coefficient of variations, and the group-wise count of significant bivariate GC connections, respectively. In general, computing within-brain directed resting-state networks using the bivariate (i.e. non-conditioned) approach yields extremely dense networks where each brain region generally appears to cause or be caused by a high number of other regions (see further for a group-wise anatomical depiction of this phenomenon). Additionally, most of the connectivity matrices show an extremely high number of causal connections from and to the cerebellum which span the whole brain. Computing within-brain directed resting-state networks using the globally conditioned approach yields a much more sparse network of causal connections which, however, exhibit an approximately equivalent pattern of inter-subject variability, possibly indicating a model-independent, physiological effect. Figure 4 shows the GCGC strengths (see equation (3.5)) of all connections in all nine subjects separately, and can be compared directly with figure 2. Figure 5 depicts, for each causal connection, the number of subjects in which that particular connection was found to be statistically significant (using a threshold of either *p*<0.01 or *p*<0.05), and can be compared directly with figure 3. Notably, the prominent bidirectional cerebellar involvement in influencing the rest of the brain is significantly decreased when using the global conditioning approach. A more direct comparison between the bivariate and the globally conditioned approach to computing within-brain resting-state networks is shown in figure 6. Here, all connections which were seen to be statistically significant in at least six subjects are depicted in the anatomical space and weighted by their average strength.

### (b) Brain–heart networks

The results of studying directed brain–heart interactions using GCGC are summarized in table 1, while the results of the same analysis using the simple bivariate approach are presented in the electronic supplementary material, table S2. For each of the ANS regressors (HF, LF, Bal), the brain regions where the number of subjects in which GCGC is statistically significant (*p*<0.05) is shown.

## 5. Discussion

To our knowledge this is the first investigation of causal brain correlates of ANS outflow. In this study, we analysed GC-based brain–brain interactions through fMRI BOLD signals, and GC-based brain–heart interactions through combined evaluation of BOLD fMRI signals and ANS-related signals. The latter were obtained by use of a PP approach, which models the waiting time between two consecutive heart beats as an inverse Gaussian distribution and represents the parameters of the distribution through an AR system. One major advantage of the PP model is that there is no need to align the heart rate signal with the sampling rate of the BOLD signal, hence avoiding any type of interpolation which would be necessary when employing any other type of time–frequency decomposition. To this end, it has been shown that resampling after classical interpolation significantly degrades the accuracy of model fitting [11]. In addition, it should be noted that the instantaneous PP model for the heartbeat incorporates physiological knowledge about the processes which underlie and modulate the heartbeat. Specifically, the model caters for the threshold behaviour of the heart's endogenous pacemaker, the history of autonomic inputs to the sinoatrial node as well as their time-dependency.

Brain signals were locally averaged according to a ROI-based approach, while heart-related regressors were computed as described in [11,17,18] and resulted in three signals: HF—an estimate of parasympathetic activity; LF—an estimate of both sympathetic and parasympathetic activity; Bal=LF/HF—an estimate of sympathovagal balance. Brain–brain GC analyses were conducted in two different ways: a ‘classical’ bivariate approach, building monovariate and bivariate AR models for the restricted and unrestricted models, respectively, and a GCGC approach, consisting of a VAR system encompassing the whole state space. Each considered link (either between any two brain regions or between any brain region and an ANS regressor) was considered causal if it was found significant (according to its *p*-value) in a sufficient number of subjects. While correlation-based (i.e. non-causal) fMRI studies at 2 T [8] and 3 T [8,17,18] and longer (approx. 3 s) TR have provided important insight into the human brain correlates of ANS modulation, this study has tackled for the first time the problem of causal link estimation between brain activity and ANS outflow. In the brain–brain part of our analysis, within-brain bivariate GC analysis results in a great number of significant directed edges (figures 3 and 6*a*), with challenging physiological interpretation. Instead, using within-brain GCGC analysis results in a reduced number of significant directed edges (see figures 5 and 6*b*) whose physiological interpretation could be tackled based on current knowledge about correlation-based resting-state networks. The resulting graph spans most brain regions while retaining a certain degree of simplicity; given the high redundancy of brain signals even in the resting state, it is plausible that the high percentage of connections which are removed through the use of the globally conditioned approach correspond to spurious/indirect causalities which should not be interpreted as direct within-brain causal influences.

In this study, using superior spatial and temporal resolution enabled by 7 T imaging (which we used in order to improve the reliability of our GC estimates), we demonstrated the existence of significant causal links between cortical as well as subcortical brain regions and ANS outflow for all three heart-related signals. While our findings cannot be directly compared with correlation-based studies, the results of the GCGC approach to studying brain–heart networks appear partially consistent with previous studies. In particular, the amygdala is seen to play a causal role (significant causal connection in at least four subjects) in the modulation of HF power (parasympathetic activity). This finding might seem counterintuitive, since the amygdala might be expected to be associated with increased heart rate, given its involvement in regulating the ‘fear’ response [80]; however, it should be noted that an increased heart rate might be induced by a reduction in the parasympathetic tone and, if the amygdala plays a role in this mechanism, this would justify the results of the amygdala having causal effects in the modulation of HF power (parasympathetic activity). This argument is consistent with a number of studies conducted in the LeDoux laboratory [81] which proposed that the amygdala may serve as a rapid (HF-related) emotional appraisal acting as a detector of potential threats and mediator of adaptive ‘fear’ responses. There is also consistent and robust evidence from animal and human research that the amygdala plays a key role in generating a series of adaptive ‘fear’ responses (e.g. fight, flight, freeze responses) which are all critically mediated via the ANS [22,23]. Our data may also partially fit with animal models of predator defence that have shown that intermediate levels of threat (e.g. when the most adaptive response would be to freeze) may cause heart rate deceleration or fear-induced bradycardia. Studies in rodents have demonstrated that this freezing behaviour and the associated increased parasympatethic tone depend on the amgydala to periacqueductal grey (PAG) projections [82]. A recent study in healthy volunteers [83] has also shown that the same neural circuits (i.e. amgydala–PAG connections) may be critically involved in humans as well to mediate fear bradycardia.

Thayer *et al.* [22] showed how areas belonging to the prefrontal cortex (the medial and ventromedial prefrontal cortex in particular) play an inhibitory role on the amygdala, and this might explain the causal link of frontal brain areas to both HF and BAL regressors; given that our results suggest that the activity of such brain areas significantly affects sympathovagal balance, it is possible that the superior frontal gyrus (left and right) and the medial frontal cortex (left and right) play a key role in such modulation, thus indicating that the modulation of sympathovagal balance is possibly intertwined with self-awareness functions [84]. Additionally, both sympathetic and parasympathetic activities are modulated by activity in cerebellar regions, in the vermis, in the hippocampus, in parahippocampal regions and crucially in the brainstem. The brainstem contains several grey matter nuclei involved in the control of autonomic functions, including the periaqueductal grey, parabrachial complex, solitary nucleus and the dorsal motor nucleus of vagus [85]. The causal relations between the cerebellum and parasympathetic outflow are in accordance with previous studies [23]. In general, our analysis does not show a large number of brain regions directly involved with sympathetic outflow. This is compatible with the idea that, as also stated in [18], the brain regions correlated with sympathetic versus parasympathetic response differ substantially, providing grounds to posit the existence of two distinct CANs with complementary roles. Also, while finding regional concordance in around four subjects (with a maximum of six in one brain region) out of nine may point towards low statistical power, it should be noted that this does not affect the validity of the results, which were detected as significant.

In this paper, we chose to use GC in the time domain as opposed to choosing one of the several frequency-domain approaches to GC, such as multivariate DC [32], DTF [33], PDC [34], gPDC [35], iPDC and iDTF [36]. Given that we used AR models containing both brain-related BOLD fMRI signals [36] (where the full bandwidth of the preprocessed signal was considered) and PP-derived ANS regressors, which represent instantaneous band-specific powers, a frequency-domain GC approach would probably have had limited advantage. Furthermore, it would have suffered from the differences in frequency domains between the two (brain and heart related) sets of signals. Additionally, the use of a frequency-domain method would have required interpolation of the unevenly sampled RR interval series, possibly introducing bias in the resulting estimates. Moreover, while estimates of parasympathetic and sympathetic activity could be readily extracted in the frequency domain, a definition of sympathovagal balance (which is considered a physiologically reliable indicator) would have been problematic in the frequency domain. Also, all VAR-based GC methods rely on *a priori* selection of the model lag length (i.e. the amount of past information to include in the model for the present observation), which, in most VAR applications, coincides with model order selection or, more specifically, the choice of AR order. While this choice can have a significant impact on the estimation of GC, and a number of criteria exist for this choice [86] (also see [27] for an extensive discussion), no general consensus exists on how to estimate AR orders in causality studies. In this context it should be noted that, while the development of algorithms for VAR parameter estimation has coaxed most published studies towards employing identical orders for the caused and for the causing variable, there is no theoretical reason for this choice, which may in fact contradict the basic concept underlying the idea of GC (i.e. that, given the restricted model, one should find the best possible unrestricted model based on some figure of merit). In this study, we, therefore, moved away from this common practice by performing independent, sequential order selection for the caused, causing and conditioning variables.

## 6. Limitations and conclusion

Recent evidence points towards the viability of GC analysis in fMRI, provided the data are acquired under a high SNR and high temporal-resolution regime. In this exploratory contribution we chose a conservative, atlas-based approach to spatially averaging brain activity. While it is well known that the mutual influence of brain regions (i.e. ‘connectivity’) plays a crucial role in physiological functioning of the brain, rendering a conditioning approach necessary, a putative mutual influence between HF, LF and Bal would have an unclear physiological basis. We therefore chose to exclude the reverse causality direction from this study. Also, while we cannot exclude that our PP-based approach may introduce some extent of autocorrelation in our time-varying spectral estimates, the latter estimates are band-specific, while the possible autocorrelation would most probably be in frequencies lower than the lower bound in LF. In this sense, the PP model is not akin to more classical strategies used to obtain time-varying estimates (e.g. time–frequency decompositions) or dynamical inference approaches [87,88] nor is it similar to transfer-function-based methods [89] or directly related to mechanistic model-based approaches to analyse ANS activity [90] or the cardiovascular system, in general [91].

Our model uses a relatively large number of parameters (around 120) with respect to the number of time points available, possibly running the risk of overparametrization. In this paper, the difference in the degrees of freedom between the restricted and the unrestricted model (five in this study) is low with respect to the difference between the number of points and the degrees of freedom of the unrestricted model (which in our case is 300–127=173). Therefore, indices calculated by comparing the two VAR residuals exclusively (as opposed to the parameters) could be less affected by overparametrization. Still, the relationship between differences in model orders and the quality of GC estimates remains to be investigated in dedicated future studies involving synthetic data simulations. In this context, in future investigations possible countermeasures could also include assuming a sparse transition matrix and fitting the VAR model using a least-squares approach with a least absolute shrinkage and selection operator (LASSO)-type penalty [92–96]. Also, in this paper we have chosen not to include instantaneous causality effects in our model. While most strategies which have been proposed to include an instantaneous effect in GC involve some degree of *a priori* knowledge about the direction of such an effect, it has been shown that the presence of instantaneous effects, whose presence cannot be excluded in the present investigation, can have detrimental effects on GC estimates [97]. Future work should aim to use those methods which tackle the issue of instantaneous GC in an assumption-free framework [29,98,99].

While including a dimensionality reduction step such as ICA or PCA [47–50] may add robustness to our globally conditioned brain–heart networks, the dimensionality reduction step (which typically estimated hundreds of components) is itself susceptible to sample size. Given the underlying nonlinear nature of interactions which occur in physiological networks, nonlinear GC techniques or kernel-based GC approaches [100] could be implemented within a conditioning approach. While these techniques may be able to capture more subtle aspects of directed brain–brain or brain–heart interactions, their implementation would pose much more severe constraints on data length and possibly temporal resolution. Also, GC techniques defined in the frequency domain (PDC [34], gPDC [35], iPDC and iDTF [36]) may offer complementary information when combined with time-domain GC in situations of low SNR or poor time resolution [66]. It should be noted that the GCGC approach is not to be limited to the study of healthy subjects or of the resting state. Future developments therefore include the study of task-dependent changes of brain–brain and brain–heart causal connectivity as well as their changes in ANS-related stimuli (e.g. pain [101]) or disease, where the study of causal brain–brain and brain–heart links could prove a sensitive instrument both in neurodegenerative disorder and in other disorders involving ANS impairment such as Parkinson's disease [102–104], epilepsy [7,105], multiple system atrophy or certain forms of mild cognitive impairment [106,107].

While the limited number of subjects employed in this study necessarily warrants some caution in casting a strong physiological hypothesis based on our globally conditioned networks, in this feasibility study we have shown that 7 T functional imaging coupled with GCGC estimates is able to quantify directed brain–heart interactions, which can be interpreted in terms of central modulation of ANS outflow while correctly disentangling the high redundancy between locally aggregated brain signals.

## Author' contributions

A.D. designed the study, performed the data analysis and wrote the manuscript. M.B. acquired, preprocessed the data and drafted the manuscript. L.P. drafted the manuscript and interpreted the results. L.L.W. acquired the data. M.G. and R.B. helped to draft the manuscript. N.T. conceived and designed the study, acquired the data, drafted the manuscript and preprocessed the data. All authors read and approved the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the Italian Ministry of University and Research (MIUR), grant no. RBFR08VABD (awarded to N.T.) and NIH NIBIB P41-RR014075.

## Footnotes

One contribution of 16 to a theme issue ‘Uncovering brain–heart information through advanced signal and image processing’.

- Accepted February 5, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.