## Abstract

It is well known that the blood oxygen level-dependent (BOLD) signal measured by functional magnetic resonance imaging (fMRI) is influenced—in addition to neuronal activity—by fluctuations in physiological signals, including arterial CO_{2}, respiration and heart rate/heart rate variability (HR/HRV). Even spontaneous fluctuations of the aforementioned physiological signals have been shown to influence the BOLD fMRI signal in a regionally specific manner. Related to this, estimates of functional connectivity between different brain regions, performed when the subject is at rest, may be confounded by the effects of physiological signal fluctuations. Moreover, resting functional connectivity has been shown to vary with respect to time (dynamic functional connectivity), with the sources of this variation not fully elucidated. In this context, we examine the relation between dynamic functional connectivity patterns and the time-varying properties of simultaneously recorded physiological signals (end-tidal CO_{2} and HR/HRV) using resting-state fMRI measurements from 12 healthy subjects. The results reveal a modulatory effect of the aforementioned physiological signals on the dynamic resting functional connectivity patterns for a number of resting-state networks (default mode network, somatosensory, visual). By using discrete wavelet decomposition, we also show that these modulation effects are more pronounced in specific frequency bands.

## 1. Introduction

An increasing number of functional magnetic resonance imaging (fMRI) studies are examining not only the task-related increases in the fMRI signal within specific brain areas, but also the correlation of time series fluctuations between different brain regions. These studies are driven by the observation that low-frequency (less than 0.1 Hz) fluctuations in the fMRI signal intensity time series, which can take place either as a result of task-induced signal modulations or in the absence of an external stimulus or explicit task (resting state), are frequently correlated between functionally related areas. The general hypothesis is that these correlated fluctuations reveal synchronized variations in the neuronal activity of a network of regions. fMRI could then provide a window into the interaction, or connection, among brain areas. The exploration of these networks by analysing coherent signal fluctuations has therefore become known as functional or effective connectivity analysis [1,2], whereby the latter term implies that causal interactions are sought [2].

In particular, the recent discovery that spontaneous brain activity is not random noise, but is specifically organized in resting-state networks (RSNs) has generated considerable interest, as RSNs provide a way to probe brain function without needing explicit tasks to drive brain activity [3]. Specifically, consistent correlations among low-frequency fluctuations between the fMRI time series corresponding to different areas at rest were initially revealed by Biswal *et al.* [4]. Subsequent studies identified numerous consistent and distinct RSNs, including motor, auditory, visual and attention networks [5,6], as well as the presence of RSNs using additional modalities (electroencephalography—EEG and magnetoencephalography—MEG) [7,8]. Correlated resting-state fluctuations have been also identified between regions of the default mode network (DMN)—a set of brain areas that regularly deactivate through a wide range of cognitive tasks [9]. This supports the view that the brain at rest comprises a number of periodically active and synchronized networks.

However, the subject test–retest reliability of RSN connectivity is lower than what would be expected if patterns were solely anatomically determined [10]. This is partly owing to the reconfiguration of functional networks, but could also reflect time-variability and the effect of non-neuronal sources. In this context, RSN non-stationarities are being increasingly considered, mainly using sliding-window techniques combined with seed voxel/region analysis and independent component analysis (ICA) [8,11–14]. Recent studies have focused on dynamic fMRI-based RSN functional connectivity [11,13,15].

Furthermore, the blood oxygen level-dependent (BOLD) signal measured by fMRI is only indirectly related to the underlying neuronal activity, as determined by neurovascular coupling mechanisms and quantified by the hemodynamic response function (HRF) [16]. The effect of non-neuronal physiological variability (e.g. heart rate, arterial CO_{2}, respiration) on fMRI measurements and connectivity has been well established [17,18]. The cerebral vasculature is exquisitely sensitive to arterial CO_{2} fluctuations; it has been shown that spontaneous arterial CO_{2} fluctuations affect both cerebral blood flow in the middle cerebral artery [19] and the fMRI BOLD signal regionally [17,20]. Additionally, slow changes in respiration depth and rate were found to significantly affect resting-state functional connectivity of the DMN [21] and HRV, as well as its low- and high-frequency power components, have been shown to be correlated to regional BOLD signal variations [18]. RSN studies are particularly susceptible to physiological noise, as they are based on a single experimental trial [22,23].

No previous studies have examined the effect of physiological factors on *dynamic* (time-varying) resting-state functional connectivity to our knowledge. Therefore, the aim of this study is to determine whether dynamic resting functional connectivity is modulated by the time-varying properties of simultaneously recorded physiological signals by using end-tidal CO_{2} (PETCO_{2}) and HR measurements obtained during scanning. In order to do so, we assess dynamic functional connectivity and its corresponding network degree for the DMN, visual and somatosensory RSNs and quantify its relation to time-varying physiological signal power. The results reveal a modulatory effect of the latter on the obtained dynamic functional connectivity patterns for both CO_{2} and HR. By using the discrete wavelet decomposition, we also show that this effect is more pronounced in specific frequency bands. Finally, we show that the observed modulation effects were considerably clearer when processing was done in the individual functional space, compared to when it was done in standard (Montreal Neurological Institute, MNI) space.

## 2. Material and methods

### (a) Experimental data

Twelve healthy volunteers (seven male, aged 29.2±4.6 years) underwent resting-state fMRI at the Cardiff University Brain Imaging Centre. Imaging was performed using a 3 T scanner (General Electric). Head movement was minimized with a bite bar. For each subject, a T1-weighted FSGPR structural scan (256×256×172 voxels of 1×1×1 mm) was acquired and used to assist in placing individual subjects' data into a common stereotactic space. During the resting-state scan, the subjects were instructed to keep their eyes closed and remain awake. The duration of the experiment was 630 s, corresponding to 210 time points (gradient echo EPI sequence, TR=3 s, TE=35 ms, FOV/slice=20.5 cm/3.2 mm, flip=90°, 53 slices with 91×109×91 voxels of 3.2×3.2×3.2 mm). PETCO_{2} was recorded using a nasal cannula attached to a gas analyser (AEI Technologies).

The HR and HRV signals were extracted from photoplethysmography waveforms obtained during the scan by the built-in scanner pulse oximeter. Specifically, before proceeding to compute these signals, we checked the raw signal quality using a signal quality index (SQI) which labels every 10 s of data as ‘reliable’ or ‘unreliable’, based on a set of physiologically relevant checks, followed by QRS template matching [24]. For the subsequent computation and analysis, we used only time series which were labelled as ‘reliable’ for more than 80% of the time. Subsequently, QRS detection was performed using the Hamilton & Tompkins [25] algorithm. The HRV signal at each beat was defined as the periods between consecutive R-peaks, whereas the HR signal was defined as the inverse of these periods multiplied by 60—in units of beats per minute. The final HR and HRV signals were obtained by applying spline interpolation at a frequency of 4 Hz and subsequently filtering with a median filter in order to remove spurious spike artefacts. Note that the results obtained by the HR and HRV signals were found to be almost identical; therefore, we present only results obtained using HR.

### (b) Functional magnetic resonance imaging data pre-processing

We implemented a processing pipeline that was developed using FSL [26] and in-house-developed Matlab scripts. Data pre-processing steps included rigid body motion correction and skull removal for functional images using BET. For each of the subjects, motion parameters were calculated to reflect head motion in six directions with respect to the mean image. The resting-state functional images were normalized into MNI standard space using linear registration implemented using FSL-FLIRT [26].

### (c) Resting-state functional connectivity analysis

#### (i) Mask-based analysis

We implemented two different versions of seed-based resting-state analysis. According to the first, we used anatomical masks defined in the MNI space (FSL atlas) in order to define regions of interest (ROIs). The following grey matter ROIs, which define the DMN, were defined: medial prefrontal cortex (mPFC), posterior cingulate cortex (PCC), anterior cingulate cortex (ACC), thalamus and precuneus. The average BOLD fMRI time series in each region were extracted and correlated with each other in order to obtain a functional connectivity matrix in sliding windows of 150 s duration (50 time lags), overlapped by 30 s (10 time lags). In order to select the window length, we obtained correlation matrices for different window lengths. Specifically, the cross-correlation between two time-series is estimated from finite data samples using the following sum
where *x*, *y* are the two time series (in our case, fMRI time series corresponding to different regions), *N* is data length and *m* is cross-correlation lag. For stationary time series, this estimate tends to the true underlying cross-correlation sequence as *N* tends to infinity (consistent estimate), whereas it becomes worse for smaller *N*, particularly for large values of *m*. A window length of 50 s was eventually selected, as it achieved a good balance between being able to track time-varying functional connectivity while providing reliable estimates (i.e. closer to their values for longer time windows) of the cross-correlation function.

The power spectral density of the PETCO_{2} and HR signals was calculated in the same sliding windows by using the Welch method [27]. To quantify correlation between different ROIs, we calculated both the absolute maximum and the absolute average cross-correlation values between 0 and ±5 time lags. In addition to performing these calculations in standard MNI space, we repeated the procedure in the functional space of each subject in order to examine the effects of registration (of the BOLD time series to the MNI space) on the results. To do so, we registered the ROI masks originally defined in MNI space to the functional space of each subject (FLIRT). After obtaining correlation values between all areas in the examined RSNs, the overall time-varying connectivity of the network was quantified by computing the graded network degree (i.e. we did not convert brain networks to binary graphs) [28] for all 150 s windows. Specifically, the graded degree is defined as
where *N* is the number of brain regions in the network and *W*_{ij} are the absolute maximum (or average) correlation values (within±5 lags) between the network nodes *i* and *j*.

#### (ii) Seed-based analysis

We performed standard seed-based analysis by first selecting a seed voxel, correlating all voxels of the brain to this voxel and applying a threshold in order to define the spatial pattern of the corresponding RSN. After examining different threshold values, we present results for a threshold value of 0.65, as they yielded spatial RSN patterns relatively similar to the mask-based analysis, whereby similarity is defined as spatial overlap between the regions defined by mask-based and seed-based analysis (i.e. number of overlapping voxels over number of voxels of the smaller mask which in this case was yielded by the seed-based analysis). Note that the precise threshold value was not found to affect the overall results, as discussed below. For the DMN, the seed voxel was selected within the PCC, whereas for the visual and somatosensory networks, the seed voxels were selected in the occipital cortex and postcentral gyrus, respectively. The time-varying RSN-graded degree for all networks was computed as described above.

#### (iii) Wavelet transformation

Wavelet transforms are multi-resolution decompositions that can be used to analyse signals and images offering good time- and scale-resolution [29,30]. In the case of one-dimensional signals, such as the ones we are dealing with, the discrete wavelet transform is essentially a decomposition into a sum of basis functions. These basis functions are shifted and dilated versions of a (bandpass) wavelet function and shifted versions of a (low-pass) scaling function. This decomposition, which is achieved via a fast filter bank algorithm and includes a digital filtering and dyadic subsampling step at each level, results in a set of wavelet and scaling coefficients. The signal can then be reconstructed via the inverse scheme [29]. By decomposing the fMRI and physiological signals via the discrete wavelet decomposition and reconstructing them using the coefficients of a single decomposition level at each time, we are basically ‘breaking down’ the signal into its different frequency sub-bands which then makes it possible to determine whether the effect of the physiological signals on the fMRI signal is more pronounced within a specific frequency sub-band. An advantage of the approach compared with other frequency-based analysis approaches is that the temporal information of the signal is preserved. In addition, application of the wavelet filter bank allows us to search for modulatory effects at different frequency sub-bands, whereas still taking the entire signal into consideration. Using the techniques described and using the wavelet toolbox in Matlab, we applied a five-level wavelet decomposition to the PETCO_{2} and HR signals, obtained the corresponding frequency decompositions and calculated correlations with the corresponding time-varying RSN degree at all the resulting frequency sub-bands. Because the frequency sub-bands resulting from a wavelet decomposition are the result of digital filtering and dyadic subsampling, they depend on the sampling rate of the signal analysed. They were, therefore, different for PETCO_{2} and HR because they were sampled at 1/3 and 4 Hz, respectively, as shown in figure 3. Lastly, it is worth mentioning that the choice of basis function did not make a significant difference in our analysis. The results presented were obtained by using the Coiflets wavelet basis [30] (figure 1).

#### (iv) Correlation coefficients

In each case, we quantified the modulating effect of time-varying physiological signal power for both PETCO_{2} and HR on the time-varying RSN degree by calculating the Spearman rank correlation coefficient, which is suitable in our case due to the fact that our results are paired with respect to time (network degree versus time-varying physiological signal power). We performed the analysis on all 12 subjects (11 for HR analysis, as one subject did not satisfy the SQI mentioned above) and then calculated the mean and standard deviation of the resulting correlation coefficients over all subjects. This was done both for the mask-based analysis (functional and MNI space), seed-based analysis as well as for different wavelet sub-bands. Finally, statistical significance for the computed correlation coefficients was assessed by calculating *p*-values for all subjects, accounting for multiple comparisons by using Hommel's correction method [31].

## 3. Results

### (a) Physiological signals and functional magnetic resonance imaging time-series

First, we note that according to the methodology outlined above for assessing signal quality (SQI), we discarded the HR data from one subject, owing to poor quality; therefore, we present results for 11 of 12 subjects for the HR signal. Representative global BOLD time series, along with the corresponding PETCO_{2} and HR time series are shown in figure 2. Note that in functional space the correlations between spontaneous PETCO_{2} fluctuations and the BOLD signal (figure 2*a*, top panel) are clearer compared with the MNI space (figure 2*a*, bottom panel), owing to the blurring effect of registering the BOLD time series to the MNI space.

### (b) Time-varying connectivity patterns and signal properties

The time-varying DMN resting functional connectivity matrices for one representative subject are shown in figure 3*a* for successive windows of 50 s length, where considerable variability can be observed. The results shown in the figure were obtained using the mask-based analysis in the individual functional space. As we can see from the corresponding colour maps, the correlation values between DMN areas are overall high. In figure 3*b*, we show the power spectral density of the PETCO_{2} and HR signals in the same successive 50 s windows, where it can be seen that there exist time variations in the spectral content of the signals.

### (c) Default-mode network

#### (i) Mask-based analysis

The DMN is shown in standard MNI space in figure 4*a* and in a representative functional space in figure 4*b*.

Figure 5 illustrates the time-varying network degree as well as the time-varying PETCO_{2} and HR signal power for a representative subject in the case of mask-based analysis in the individual functional and MNI spaces. In the case shown here, the degree was computed for the case in which correlations between areas were quantified using the average absolute cross-correlation value between 0 and ±5 time lags (similar results were obtained when the absolute maximum cross-correlation value was used). The results were similar for smaller/longer window sizes but a window length of 50 s, as mentioned in Materials and methods, achieved a good balance between time resolution and reliable cross-correlation function estimation. We can observe considerable differences between the results in the functional and MNI spaces. Specifically, in figure 5*a*,*c*, it can be seen that the time-varying DMN degree follows the time-varying HR and PETCO_{2} power, respectively. On the other hand, in figure 5*b*,*d*, the same relation is not as clear, which is reflected on the obtained correlation coefficients (table 1).

In figure 5, we also present results obtained when using the wavelet transform, for the frequency sub-bands within which the correlation between time-varying network degree and band-limited signal power was found to be highest. In both functional and MNI spaces, the A4 level (0–0.25 Hz) yielded the highest correlations for the HR signal, whereas for the PETCO_{2} signal, the A2 level (0–0.08 Hz) yielded the highest correlations in the functional space, owing to that resting PETCO_{2} fluctuations exhibit most of their power below 0.05 Hz [32]. In the MNI space, level A4 (0–0.02 Hz) yielded the stronger correlations for the PETCO_{2} signal, possibly due to the fact that registration of the BOLD time series to the MNI space essentially results in low pass filtering. Overall, the correlations between time-varying degree and signal power were found to be stronger for the time-varying band-limited power compared with total signal power (see also tables 1 and 2).

The correlation coefficient values between time-varying DMN degree and the total/band-limited power of the PETCO_{2} and HR signals are given in table 1 for all subjects (mean±s.d.). These suggest the presence of temporal correlations between network degree and signal power. Correlations were found to be much higher in the functional space compared with the MNI space. The correlation coefficients for the band-limited power are given for the wavelet sub-bands that yielded the highest correlation values (A4 (0–0.25 Hz) for HR in both spaces, A2 (0–0.08 Hz) and A4 (0–0.02 Hz) for PETCO_{2} in the functional and MNI spaces, respectively). In both cases, using the band-limited time-varying signal power yielded higher mean correlations. For instance, the correlation coefficients for PETCO_{2} increased from around 0.6 (total signal power) to almost 0.8 (A2 level) in functional space. The values obtained in functional space are also considerably higher, suggesting that registration to the MNI space blurs the modulatory effect of PETCO_{2} on time-varying network degree.

Table 2 illustrates the obtained *p*-values for the Spearman correlation coefficients for all subjects in each space/sub-band. In agreement to the above observations, statistical significance (*p*<0.05) was observed for most subjects in the individual functional space when total signal power was used and for almost all subjects (10 out of 12 for PETCO_{2}, 11 out of 11 for HR) when the corresponding optimal sub-band was used. On the other hand, statistical significance was not achieved for most subjects in the MNI space, even when the optimal frequency sub-band was used.

### (d) Seed-based analysis

In the present and next sections, we present results obtained from the individual functional space of each subject, because it was found that registration to the MNI space and the resulting blurring of the corresponding BOLD time series considerably affects the results, as shown above. The spatial overlap between mask-based and seed-based analysis for the DMN areas was equal to 94.82% for the ACC, 44.65% for the PCC, 76.33% for the thalamus, 83.15% for the precuneus and 92.85% for the medial prefrontal cortex ROIs. The Spearman correlation coefficients between time-varying DMN degree and band-limited PETCO_{2}/HR power are given in table 3 for all subjects (mean±standard deviation). These suggest the presence of temporal correlations between network degree and the PETCO_{2}/HR power, as in the case of mask-based analysis (somewhat lower values—not shown—were obtained for total power, as above). As before, correlations were found to be much higher in the wavelet sub-band A2 for PETCO_{2} (0–0.08 Hz) and in the wavelet sub-band A4 for HR (0–0.25 Hz). Table 4 illustrates the *p*-values for the obtained Spearman correlation coefficients for all subjects in each space/sub-band. Overall, the results from the seed-based analysis are similar to the mask-based analysis, with statistical significance attained in most subjects (nine out of 12 for PETCO_{2}, nine out of 11 for HR; figure 6).

### (e) Visual and somatosensory networks

We applied the procedure described above (seed-based analysis only) on two additional RSNs, namely the visual and somatosensory networks (figure 7). The results were overall similar to those obtained for the DMN and revealed a modulatory effect for both PETCO_{2} and HR on the time-varying degree for both networks. The wavelet-based analysis yielded higher correlations than those achieved for the total signal power for the same frequency sub-bands as above (tables 5 and 6). Statistically significant correlations between the time-varying visual and somatosensory RSN degree and band-limited HR and PETCO_{2} power were obtained for all subjects in the case of the visual cortex, whereas, interestingly, for the somatosensory cortex significance was obtained for all subjects in the case of HR but for only two out of 12 subjects for PETCO_{2}—even though results were marginally significant for all the remaining subjects (table 6).

## 4. Discussion

This paper examined the modulation of dynamic, resting-state connectivity by physiological signal fluctuations. The results reveal a clear effect of time-varying signal power for both HR and PETCO_{2} on the time-varying network degree for three well-described RSNs: the default-mode, visual and somatosensory RSNs, revealing brain–heart interactions in the context of fMRI-based RSN connectivity studies. This effect was found to be more pronounced for the fluctuations in the physiological spectral content in specific frequency sub-bands (time-varying band-limited power), as revealed by wavelet analysis. In addition, it was found that the observed modulations were not apparent when the analysis was performed in the MNI space, using anatomical masks to define the RSNs of interest. Despite the well-established effect of physiological signals on fMRI connectivity and particular RSN connectivity [17,18,21–23], to the best of our knowledge, this is the first study that demonstrates that the time-varying properties of physiological signals may affect dynamic resting-state functional connectivity. This has important implications, as it suggests that even moderate modulations in the power of these signals can considerably influence RSN analyses and that a significant source of dynamic variations in resting-state connectivity is physiological in nature. Given that resting-state (spontaneous) fluctuations of physiological signals such as HR and PETCO_{2} are of small magnitude, they are not expected to significantly affect neuronal activity *per se*; therefore, the observed modulatory effects are likely physiological in origin.

Fluctuations in fMRI-based functional brain networks have been observed in time scales from tens of seconds to several minutes [11,12,15]. For instance, recent fMRI studies with high temporal resolution [15] have demonstrated the presence of time-varying patterns that are related to large-scale topological properties of the brain. It has been speculated that these fluctuations may achieve more efficient information transfer and energy expenditure. We selected the sliding window length (50 time lags) as a good trade-off between time resolution and obtaining a good cross-correlation estimate. Interestingly, several studies related to dynamic functional connectivity have used sliding windows of around 50 s (some with lower TR values than 3 s) [11,33], whereas in [15], 60 s windows were used, albeit with a subsecond TR (60 s=83 time lags). In addition, in [34], it was suggested that average correlation values within and between RSNs stabilize at approximately 240 s. Therefore, the selected window length is reasonable given these observations and is similar to the window lengths examined in other studies [35].

In order to obtain an overall measure of dynamic functional RSN connectivity, we computed the time-varying average degree of the brain network defined by the corresponding RSN, as this graph-theoretic measure is a straightforward way to quantify overall connectivity in a brain network and to monitor its variations over time. Alternative graph measures such as clustering coefficient or efficiency are not relevant for our purpose, as we are interested in how ‘connected’ a specific RSN is within a given time window. Note that we computed the *graded* degree [28], i.e. we did not convert the networks to binary networks, as often done in practice. This is due to that we were interested in the precise degree of modulation of dynamic connectivity by physiological fluctuations, without being concerned if and which network connections are deemed ‘significant’ (non-zero), which would further complicate matters and importantly would also depend on the method of network binarization (i.e. simple thresholding, surrogate data [36], etc.). Furthermore, we selected to quantify functional and not effective connectivity, as the use of the latter has been questioned in the case of fMRI-based connectivity (particularly RSN connectivity), owing to the low-pass filtering introduced by the HRF, which limits the time resolution and the ability to infer causal effects [37]. However, it is worth noting that methods for estimating the HRF from resting-state data have been recently proposed [38], which could in turn provide more reliable effective connectivity estimates from fMRI RSN data [37,39]. This is an important point that deserves further attention.

We examined two different methods to define RSNs (mask-based analysis and seed-based analysis) as the network node definition procedure has been shown to influence both connectivity analyses [37] and statistical comparisons between different conditions [40] in fMRI studies. Here, it was not found to influence the results to a large degree. On the other hand, registration of the BOLD time series to the MNI atlas was found to blur the modulatory effects of physiological fluctuations to a large extent. This suggests that sensitivity in tracking dynamic functional connectivity changes is lost when working in the MNI space and that possibly the overall connectivity patterns may not be as accurate as in the individual functional space. Overall, our results were found to be similar among all RSNs examined, with HR yielding overall more pronounced modulations of time-varying degree than PETCO_{2} (tables 2, 4 and 6). Interestingly, it was found that the somatosensory network yielded marginally significant results for PETCO_{2} for most subjects in comparison with other networks (table 6), which could be due to vascular anatomy differences between these networks.

With regard to the effect of threshold value in the seed-based analysis, we repeated our analysis for different threshold values. For values between around the reported value (between 0.55 and 0.75), the resulting areas from seed-based analysis were relatively cohesive, and the results were very similar to those reported above. For larger threshold values, the seed-based areas become more disjointed; however, they can be constrained anatomically to yield the same number of DMN areas (ACC, PCC, etc.), and the physiological modulations were found to be overall similar. For smaller threshold values, the seed-based areas become more extensive (and possibly overlapping), but the effect of physiological modulations is still present. Finally, extreme threshold values yield either a few voxels or very extensive areas and are not of much practical interest.

In order to examine the effects of standard physiological correction methods on the results, we repeated the analysis after implementing two different schemes of physiological correction (RETROICOR [41]) to regress out the effects of (i) HR and respiration and (ii) HR, respiration and PETCO_{2}. The results are summarized in table 7, where we show results (Spearman correlation coefficients) for the broadband signals (i.e. without using wavelets).

When HR and respiration were regressed out, the modulatory effects of both signals were found to be reduced but not completely removed. When all signals are regressed out, the drop in modulation is more substantial for CO_{2}. Therefore, modulating physiological effects do not disappear even after physiological correction—particularly for HR, for which the correlation coefficient remains fairly high (0.53). Moreover, it is worth noting that, whereas in most fMRI studies, HR and respiration are regressed out using RETROICOR or other similar techniques, PETCO_{2} data are typically not recorded and consequently, the effects of CO_{2} are not regressed out. The results suggest that collecting PETCO_{2} data in resting-state studies is important and also that better physiological denoising algorithms for such data are needed. Furthermore, we selected not to perform white matter (WM) or cerebrospinal fluid (CSF) regression, which is a relatively common practice in similar studies, as it is rather difficult to obtain accurate individualized masks for WM and/or CSF and consequently, some fraction of the grey matter signal may be regressed from the data, masking out the pure effects of physiological noise, in which we are interested [22].

The frequency sub-bands that yielded the stronger correlations between time-varying degree and band-limited power are related to each signal's characteristic and suggest that fluctuations in these bands have a clearer modulatory effect. For PETCO_{2} modulations in very low-frequency power (0–0.08 Hz in functional space and 0–0.02 in MNI space) were found to have the clearer effect, which is consistent with the previously described spectral characteristics of PETCO_{2} fluctuations [19,32]. For HR, a wider frequency band was identified (0–0.25 Hz) reflecting the fact that HR has a richer spectral content than PETCO_{2}. Additionally, the identified frequency band contains the low-frequency (LF—around 0.1 Hz) spectral peak and possibly also the high-frequency (HF—around 0.2 Hz) spectral peak of the HR signal. While, often, these two spectral peaks are assumed to correspond to cardiac sympathetic and cardiac parasympathetic neural activity, respectively, the relative contributions of each mechanism still remain the subject of controversy and investigation. What is generally accepted is that both LF and HF components of the HR signal (or HRV—note that fluctuations in these two signals have the same spectral characteristics) are affected by the complex interactions between both parasympathetic and sympathetic nerve fibres as well as mechanical, and other factors on the pacemaker cells usually located in the sinoatrial node [42]. The effects of these two components could be disentangled by obtaining estimates of instantaneous HF and LF HR power by using, e.g. point-process models [43].

## Author' contributions

F.N. performed data analysis, prepared all figures and wrote the paper. C.O. performed data analysis and wrote the paper. P.P. performed data analysis. K.M. collected the data and provided comments on the paper. R.G.W. collected the data and provided comments on the paper. G.D.M. conceived and designed the study, wrote the paper and approved the final version.

## Competing interests

We declare we have no competing interests.

## Funding

We received no funding for this study.

## Footnotes

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

- Accepted February 11, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.