## Abstract

In order to examine the effect of changes in heart rate (HR) upon cerebral perfusion and autoregulation, we include the HR signal recorded from 18 control subjects as a third input in a two-input model of cerebral haemodynamics that has been used previously to quantify the dynamic effects of changes in arterial blood pressure and end-tidal CO_{2} upon cerebral blood flow velocity (CBFV) measured at the middle cerebral arteries via transcranial Doppler ultrasound. It is shown that the inclusion of HR as a third input reduces the output prediction error in a statistically significant manner, which implies that there is a functional connection between HR changes and CBFV. The inclusion of nonlinearities in the model causes further statistically significant reduction of the output prediction error. To achieve this task, we employ the concept of principal dynamic modes (PDMs) that yields dynamic nonlinear models of multi-input systems using relatively short data records. The obtained PDMs suggest model-driven quantitative hypotheses for the role of sympathetic and parasympathetic activity (corresponding to distinct PDMs) in the underlying physiological mechanisms by virtue of their relative contributions to the model output. These relative PDM contributions are subject-specific and, therefore, may be used to assess personalized characteristics for diagnostic purposes.

## 1. Introduction

Following the advent of transcranial Doppler ultrasound (TCD) for measuring cerebral blood flow velocity (CBFV) with high temporal resolution [1], the study of cerebral haemodynamics has often benefited from data-based input–output predictive models that describe quantitatively how measured physiological variables of interest interact dynamically with each other. Over the past 25 years, a particular focus has been the study of cerebral autoregulation, which is defined as the ability of the cerebrovascular bed to maintain a constant cerebral blood flow in response to pressure changes [2]. Specifically, dynamic autoregulation has been studied by examining the relation between step (e.g. thigh cuff deflations) or spontaneous arterial blood pressure (ABP) changes and the corresponding CBFV variations [3–12]. Furthermore, it is well known that the cerebrovascular bed is exquisitely sensitive to arterial CO_{2} changes; spontaneous fluctuations of end-tidal CO_{2} (ETCO_{2}) have been shown to be correlated to both global blood flow using TCD [13,14] as well as regional blood flow using functional magnetic resonance imaging [15]. Therefore, an approach that has gained considerable scientific traction in recent years is to quantify how spontaneous changes in ABP and ETCO_{2}, viewed as two concurrent inputs, cause changes in the output variable of CBFV [13,16–20]. For instance, this approach has been recently used to assess the effect of orthostatic stress [17], as well as possible impairments of cerebral autoregulation or vasomotor reactivity in patients with Alzheimer disease [19] and amnestic mild cognitive impairment [20].

A key concept used in this data-based modelling approach is the set of principal dynamic modes (PDMs) that are the main characteristics of each given dynamic nonlinear system and make the representation, as well as estimation, of nonlinear Volterra-type models feasible and reliable [21]. The PDMs, furthermore, reveal the key functional features of the subject system and, consequently, facilitate model interpretation and offer unique insights into the functional mechanisms of the system. The latter benefit may prove valuable in deciphering the possible physiological causes of a given disease and, therefore, assist clinical diagnosis as well as suggest potential treatments. Last, but not least, this capability can be useful in assessing the effects of treatments or pharmaceuticals in a cohort of patients relative to appropriate controls. This multitude of potential benefits has given recent impetus to this approach, although obstacles remain in terms of data availability (owing to practical constraints) and peer acceptance owing to its relative (perceived) mathematical complexity. This paper seeks to extend the aforementioned model to include changes in heart rate (HR) as a third input and explore the concurrent dynamic effects of HR changes upon CBFV using a three-input model of cerebral haemodynamics (inputs: ABP, ETCO_{2}, HR; output: CBFV).

The rationale for this three-input/one-output modelling objective is twofold: (i) the desire to assess quantitatively/statistically the existence of dynamic effects of HR changes upon CBFV and, if present, to elucidate the particular characteristics of these dynamic effects with regard to plausible physiological hypotheses by examining the form of the obtained PDMs; (ii) the extent to which HR and the related dynamic changes in cardiac output may have direct influence on cerebral blood flow independent of changes in ABP and ETCO_{2}. It is evident that the motivation for pursuing these objectives is the development of model-based diagnostic procedures for diseases with a significant cerebrovascular component, as well as the potential discovery of improved treatments that rely on a better understanding of the underlying physiological mechanisms.

## 2. Methods

### (a) Data collection and preparation

Time-series data were collected in 18 healthy subjects (nine men and nine women, age 66.8±7.4 years), who participated voluntarily in this study at the Institute for Exercise and Environmental Medicine of the University of Texas Southwestern Medical Center, Dallas, TX, USA, and signed the informed consent form that was approved by the institutional review board. The data were collected in a quiet, environmentally controlled laboratory under supine resting conditions. After 20 min of supine rest, 5 min of recordings were made at an initial sampling rate of 1 kHz. ABP was measured continuously with finger photo-plethysmography (Finapres), ETCO_{2} tension was obtained via a nasal cannula using capnography (Criticare Systems), and HR was monitored by electrocardiography. CBFV was measured in the middle cerebral arteries using a 2 MHz transcranial Doppler (TCD) probe (Multiflow, DWL) placed over the temporal window and fixed at a constant angle with a custom-made holder. All measurements are non-invasive, safe and comfortable for the subjects. The highly sampled data were subsequently reduced to beat-to-beat measurements using the procedures that have been outlined in our previous publications [19,20]. Figure 1 shows illustrative time-series data over 5 min for one of the subjects.

### (b) Data analysis and modelling

The beat-to-beat data of each control subject were analysed using the method of PDMs to obtain predictive models of the dynamic effects of changes in the three inputs of ABP, ETCO_{2} and HR upon the output CBFV. We briefly outline below the PDM-based modelling approach. For the many mathematical and technical details of PDM and Volterra-type modelling, the reader is referred to the monograph [21] and to our publications presenting its application to cerebral haemodynamics [19,20].

The PDM-based modelling methodology has its rigorous foundations in the general input–output nonlinear modelling approach initially proposed by Norbert Wiener for Gaussian white-noise inputs in 1958 [22] that was based on the general mathematical foundation of the Volterra functional expansion [23]
where *y*(*t*) and *x*(*t*) denote the system output and input, respectively, and {*k*_{r}} is the *r*th-order Volterra kernel that describes the *r*th-order nonlinearities of the system. The modelling task is the estimation of the Volterra kernels from input–output data. The Volterra modelling approach is applicable to all finite-memory dynamic nonlinear systems, which covers almost all physiological systems (except for chaotic systems or non-dissipating oscillators). The Volterra model has been extended to the case of multiple inputs [21], as in this study. The PDM-based modelling methodology was developed over the past 30 years as a succession of practical adaptations of the original Volterra–Wiener theory to the constraints imposed by actual physiological applications [21]. The driving goal was the reliable estimation of dynamic nonlinear models from relatively short records of spontaneous (or evoked) physiological time-series data. The PDM modelling task commences with the estimation of Volterra models using the Laguerre expansion of kernels [24]. This yields the modified Volterra model for single input and single output [21]
and
where *n* denotes the discrete-time index (*t*=*nT*), *T* is the sampling interval, *Q* is the order of system nonlinearity, *M* is the system memory, *L* is the number of employed Laguerre basis functions {*b*_{j}} and {*c*_{r}} denotes the *r*th-order kernel expansion coefficients. The residual * ε*(

*n*) is the model prediction error, which we seek to minimize in the normalized mean-square error (NMSE) sense, defined as The Laguerre expansion coefficients of the modified Volterra model can be estimated through least-squares regression and yield estimates of the Volterra kernels. The form of the modified Volterra model can be extended to include multiple inputs—e.g. three inputs in this study. Upon estimation of the Volterra kernels for each input, a rectangular matrix is composed of the estimated kernels of each input over the entire cohort for which global PDMs are sought. Singular value decomposition (SVD) of this rectangular cohort kernel matrix reveals the significant singular vectors (corresponding to the significant singular values according to a threshold criterion) that are the waveforms that can represent all the kernels in a manner that balances accuracy with parsimony. These singular vectors are the global PDMs.

In this study, four discrete Laguerre functions were used for each input with the Laguerre alpha parameters 0.5, 0.8 and 0.85 for ABP, HR and ETCO_{2}, respectively, which yield the minimum average prediction error over all subjects following a search procedure over a grid of alpha values with increments 0.05. The appropriate number of Laguerre functions was selected on the basis of the Bayesian information criterion (BIC), which takes into consideration the total number of free parameters in the respective Laguerre-based Volterra model. The latter is (*L*+*Q*)!/*L*!*Q*!, where *L* denotes the number of Laguerre functions and *Q* is the nonlinear order of the Volterra model. Like all other criteria derived from statistical hypothesis testing, this criterion is a sufficient (but not necessary) condition for rejecting the null hypothesis of non-significant reduction in residual variance.

Having obtained the kernel estimates for the cohort subjects, we apply the aforementioned SVD procedure on the rectangular cohort kernel matrix (simply composed of all kernel estimates) to obtain the PDMs. Three PDMs were selected in this case, using the threshold criterion of the respective singular values being larger than 10% of the maximum singular value. The obtained PDMs form a filter-bank for each input that receives the respective input signal and generates (via convolution) signals, which are subsequently transformed by the associated nonlinear functions (ANFs) that represent the nonlinear characteristics of the system for each PDM. The coefficients of the ANFs (cubic polynomials in this application) are estimated via multi-linear regression of the PDM outputs and their powers upon the output signal—because the outputs of the ANFs are summed to form the model output prediction as (also depicted schematically in figure 2)
where {*u*_{m},*j*}, {*u*_{h}(*t*)} and {*u*_{m}(*t*)} are the PDM outputs (i.e. convolutions of the *m*th input with the *j*th PDM) for the three inputs, and *f*_{m,j} is the ANF associated with the PDM output *u*_{m,j}. Thus, the PDM-based model separates the dynamics (PDMs) from the nonlinearities (ANFs). Because the separability of the system nonlinearities cannot be always assumed, we must generally include in the PDM-based model cross-terms in the form of pair-products of ANF outputs that have a correlation coefficient with the output signal that exceeds the 99% significance level in a statistical hypothesis test using the *w-statistic* (based on the Fisher transformation of the correlation coefficient estimate). In this study, only four subjects out of 18 exhibited a single statistically significant pair-product between ANF outputs (see Results section). Because statistically significant cross-terms were not found in most of the subjects, we have not included cross-terms in the structure of the PDM-based model used for all subjects (figure 2).

As indicated in figure 2, the output of the PDM-based model is formed by additive signal components that are generated by cascaded operations of convolutions of each input signal with its respective PDMs and subsequent static nonlinear (cubic in this case) transformations by the respective ANFs. The obtained PDMs constitute the common functional basis for the representation of the model dynamics for all subjects (often termed global PDMs), resulting in compact dynamic nonlinear model representations. The ability of the global PDMs to represent the dynamics of the entire ensemble of subjects is validated by the predictive accuracy of the global PDM-based model for each subject.

Three global PDMs and cubic ANFs were found to be adequate for each input of this model, according to the BIC. The use of three global PDMs for each input and cubic ANFs keeps the total number of free parameters for the three-input PDM-based model low (i.e. 27 plus an offset constant). We emphasize that the global PDMs are common for all subjects, but the estimated ANFs associated with each global PDM are *subject-specific* and can be used to characterize uniquely the cerebral haemodynamic characteristics of each subject. This is the basis for the potential utility of the PDM-based modelling approach for clinical diagnosis.

## 3. Results

We begin with the statistical assessment of the effect of including the third HR input and/or the transition from a linear to a nonlinear model. We pursue this by using the BIC, which balances the reduction in NMSE of the output prediction by the three-input and/or nonlinear model with the number of free parameters contained in the respective models. Under the assumption that the model prediction errors are independent and identically distributed according to a normal distribution, the BIC is defined as

where is the prediction error variance, *N* is the number of the time-series data and *P* is the number of free parameters in the model. We note that the number of free parameters is 7 for the two-input linear models, 19 for the two-input (cubic) nonlinear model, 10 for the three-input linear model and 28 for the three-input (cubic) nonlinear model. Table 1 summarizes the mean and standard deviation (s.d.) values of the BIC difference between various model pairs over all 18 subjects, as well as the corresponding *p*-values for the paired *t*-test of BIC reduction. As indicated in table 1, the BIC reduction is statistically significant for the three-input linear or nonlinear model relative to the respective linear or nonlinear two-input model (*p*=0.00525 for linear and *p*=0.00012 for nonlinear), as well as for the transition from linear to nonlinear model with two inputs or three inputs (*p*=0.00016 for two inputs and *p*=0.00002 for three inputs). Therefore, the BIC corroborates the potential utility of using three-input nonlinear models, instead of two-input linear or nonlinear models.

Following the procedure outlined in the Methods, we obtained the global PDMs from the cohort of 18 control subjects that are shown in figure 3 for the ABP input figure 4 for the ETCO_{2} input and figure 5 for the HR input (in the time and frequency domains). We observe that the global PDMs exhibit distinctive spectral characteristics in terms of resonant peaks, which may attain importance in the interpretation of the PDM-based model and the functional properties of this system (see Discussion). Specifically, we observe the following.

For the ABP input (figure 3)

— the first PDM exhibits a high-pass characteristic (above 0.1 Hz) akin to the well-known Windkessel model of passive fluid-mechanical admittance of the cerebral vasculature;

— the second PDM exhibits a resonant peak around 0.2 Hz, which is probably related to the respiratory (mechanical) effect on the pulmonary vasculature; and

— the third PDM exhibits a resonant peak around 0.12 Hz, which is posited to be related to the sympathetic activity, based on the accepted view of the origin of HR spectral peaks.

For the ETCO_{2} input (figure 4)

— the first PDM exhibits a low-pass characteristic (below 0.03 Hz), combined with a resonant peak around 0.08 Hz;

— the second PDM also exhibits a low-pass characteristic (below 0.05 Hz), combined with a resonant peak around 0.15 Hz; and

— the third PDM exhibits a resonant peak around 0.045 Hz.

For the HR input (figure 5)

— the first PDM exhibits a low-pass characteristic (below 0.03 Hz) similar to the first PDM of ETCO

_{2}, combined with a resonant peak around 0.12 Hz;— the second PDM exhibits a low-pass characteristic (below 0.05 Hz) similar to the second PDM of ETCO

_{2}, combined with a resonant peak around 0.18 Hz probably related to modulation by parasympathetic activity; and— the third PDM exhibits a resonant peak around 0.045 Hz similar to the third PDM of ETCO

_{2}.

The PDMs are obtained by performing singular value decomposition on the matrix containing the first- and second-order Volterra kernels from all subjects. In order to illustrate their variability, we used the bootstrap method to obtain 90% confidence intervals [25]. The results are shown in figure 6, where it can be seen that the PDMs corresponding to ABP exhibit the least variability, whereas the PDMs corresponding to HR exhibit the most.

Unlike the PDMs, the (cubic) ANFs are estimated via least-squares regression and, therefore, the three coefficients (linear, quadratic and cubic) are statistical estimates with some probability density function. The average (cubic) ANFs defined by the average coefficient estimates for each of the PDM branches of the three-input model are shown in the model diagram of figure 2. It is evident that ETCO_{2} and HR have the strongest nonlinearities. The ANFs for the HR input exhibit symmetry about the origin (operating point); however, only one ANF for the ETCO_{2} input exhibits such symmetry (see Discussion). Our previous studies have shown that most of the nonlinear characteristics of the cerebral haemodynamics are exhibited in the low-frequency range—particularly below 0.05 Hz [13,17,19,20]—with regard to the ETCO_{2} input. Therefore, the nonlinearities are expected to be associated primarily with the PDMs that exhibit resonant peaks in that frequency range. This is, in fact, observed in the ANFs shown in the PDM-based model of figure 2. However, when we assess the statistical significance of the obtained estimates of the ANF coefficients using *p*-values (under the assumption of system stationarity), we observe significant nonlinearity only in the quadratic coefficient of the first PDM of the ETCO_{2} input (*p*=0.039), as indicated in table 2 where the average (s.d.) values of these coefficient estimates and the corresponding *p*-values are shown for the three PDMs of each input. It is also seen in table 2 that several linear coefficient estimates do not rise to the level of statistical significance. In the Discussion section, we consider the possible requirement of longer data records and/or larger sets of subjects with respect to the statistical significance of nonlinear estimates, as well as the possible effects of system non-stationarities on these results.

Although the PDM-based model separates the dynamics (PDMs) from the nonlinearities (ANFs), it cannot be generally assumed that the system nonlinearities are ‘separable’ in terms of the specified PDM–ANF pathways—i.e. cross-interactions may generally exist between ANF outputs of the model. For this reason, we must generally include in the PDM-based model cross-terms in the form of pair-products of ANF outputs, which have a statistically significant correlation coefficient with the output signal, based on hypothesis testing of the *w-statistic* which is the Fisher transformation of the correlation coefficient estimate. In this study, only two subjects out of the 18 exhibited a single statistically significant pair-product between ANF outputs of two different inputs (PDM2 of HR with PDM1 of ETCO_{2} in subject 9 and PDM1 of ABP with PDM3 of ETCO_{2} in subject 14). Two other subjects also exhibited a single statistically significant pair-product between two ANF outputs of the ETCO_{2} input (PDM1 with PDM2 of ETCO_{2} in subject 4 and PDM1 with PDM3 of ETCO_{2} in subject 7). Because such statistically significant cross-terms were not consistently found in most of the subjects, we have not included cross-terms in the structure of the PDM-based model used for all subjects (figure 2).

To illustrate the relative contributions of the PDM–ANF branches of each input to the model prediction of the output, we show in figure 7 the separate contribution of each input (i.e. the sum of the respective three ANF outputs for each input) in the time and frequency domains, along with the total model prediction and the actual output signal for one subject. It is evident that all three inputs make significant contributions, although with different spectral characteristics. The ETCO_{2} contribution is more pronounced in lower frequencies and does not show the effects of respiratory sinus arrhythmia evident in the contributions of the other two inputs. The ABP contribution is the strongest (average root-mean-square (RMS) value of 1.86), with the ETCO_{2} contribution second (average RMS value of 1.41) and the HR contribution third (average RMS value of 0.92).

## 4. Discussion and conclusions

The PDM modelling study of cerebral haemodynamics in 18 control subjects (nine men and nine women) has shown that the inclusion of HR as a third input (along with ABP and ETCO_{2}) in the model that predicts the CBFV output causes a statistically significant reduction in the resulting BIC value for all subjects. We note that the BIC balances the changes in prediction NMSE and the number of free parameters in the respective model. This implies that there exist contributions of HR changes to CBFV changes that are independent from the contributions of ABP and ETCO_{2} changes, which are included as concurrent inputs of the model. Furthermore, this finding demonstrates that the HR dynamic effects on CBFV changes can be captured and quantified, along with the concurrent effects of ABP and ETCO_{2}, by a three-input PDM-based model that can be estimated from time-series data of each subject. It is noted that the inclusion of HR as a third input does *not* alter the PDMs of the other two inputs obtained by use of the two-input model. This is additional evidence of the independent influence of HR changes.

An important question that arises from critical contemplation of the proposed approach is whether the obtained PDMs (being singular vectors derived from SVD of the kernel matrix, and not statistical estimates) reflect actual physiological entities or they are simply by-products of the computational method. This valid concern can be addressed through statistical analysis of the ANF coefficient estimates that are associated functionally with the PDMs, because these coefficient estimates (obtained via least-squares regression) directly define the contribution of each PDM to the model prediction. When we performed the traditional analysis of statistical significance (*p*-values) on the obtained ANF coefficient estimates, only a subset of these coefficients were shown to be statistically significant (table 2) despite the fact that nonlinear models were found to satisfy the BIC as discussed above. This may imply that a larger sample of subjects and/or longer data records may be required to demonstrate the statistical significance of some of these coefficient estimates (because the statistical hypothesis testing is a sufficient but not necessary condition). It may also indicate that the presence of possible system non-stationarities, which has been investigated recently [26,27], may affect the results of traditional hypothesis testing (*p*-values).

Examination of the relative contributions of the three inputs (in terms of their average RMS values over all subjects) to the prediction of the output reveals that ABP makes the largest contribution (RMS mean = 1.78) with the ETCO_{2} second (RMS mean = 1.37) and the HR third (RMS mean = 0.93). The importance and independence of the contribution of HR changes to the dynamic changes in CBFV is illustrated in figure 5 for one subject, along with the relative contributions of ABP and ETCO_{2} changes. It is evident that the ETCO_{2} contribution is mainly in lower frequencies (less than 0.05 Hz), whereas the other two inputs make contributions over a similar frequency range (as illustrated in the respective spectra shown in figure 5) and exhibit clearly the effects of respiratory sinus arrhythmia—although the latter effects are about six times stronger in the ABP spectrum than in the HR spectrum. These findings demonstrate the potential of this modelling approach to explain the dynamic changes in CBFV under the combined influence of concurrent changes of all three inputs and, thereby, achieve a reliable representation of the haemodynamics of each subject that can be used potentially for clinical diagnosis of cerebrovascular disease.

To the best of our knowledge, one study has considered the effects of three different inputs on CBFV. Specifically, a recent study by Panerai *et al*. [18] examined the combined effects of ABP and ETCO_{2} changes on CBFV changes in response to a motor task (elbow flexion), aiming to disentangle physiological effects from neurovascular coupling (effect of motor task on CBFV). Therefore, the effects of HR were not considered and the effect of the motor task was incorporated as a third input in their multiple-input autoregressive model using a perfect step function, i.e. an idealized signal that is rather different in nature than experimental measurements of physiological fluctuations, which exhibit broadband characteristics. With respect to the effects of ABP and ETCO_{2} on CBFV variability, their results overall agree with this study, in that ABP explains mainly fast CBFV variations and ETCO_{2} explains the slow CBFV variations (e.g. figs 6 and 3 in [18]). ETCO_{2} was found to have a slightly larger contribution to CBFV variability than ABP—albeit more variable across subjects and subsequently less statistically significant (table 1 in [18]). However, this could be due to the fact that there appears to be a strong correlation between the third input of their model (elbow flexion) and ABP, and also to that the mean value of ABP and CBFV increased during elbow flexion. Overall, our results regarding the effects of ETCO_{2} agree with previous two-input studies using Volterra models [13] and multiple coherence [28,29].

A note should be made regarding the use of NMSE as a criterion for model performance and a means for model validation. Although most investigators would accept the ability of a model to reduce the residual variance as evidence of improved model performance, it must be noted that a reduction in residual variance may be accompanied by an increase in estimation variance, which is a detriment to model performance. The balance in this delicate trade-off often depends on the statistical characteristics of the data-contaminating noise, which are generally unknown; one way to achieve a good trade-off is to use a statistical criterion such as the BIC.

The PDMs of this model show promise in revealing the important dynamic characteristics of the underlying physiological mechanisms and allow the formulation of specific quantitative hypotheses about the role of the sympathetic and parasympathetic activity (corresponding to specific PDMs) in cerebral perfusion and autoregulation. Consistent with the widely held view regarding the temporal characteristics of autonomic control of the cardiovascular system, it is posited that the HR PDMs with a resonant spectral peak around 0.2 Hz may be associated with parasympathetic activity, whereas the PDMs with a resonant spectral peak around 0.12 Hz may be associated with sympathetic activity. Hamner & Tan [8] investigated the relative contributions of sympathetic, cholinergic and myogenic mechanisms to cerebral autoregulation using separate blockades of alpha-adrenergic receptors, muscarinic receptors and smooth-muscle calcium channels, respectively. Their data were collected during the application of periodic low-body negative pressure (LBNP) at a frequency of 0.03 Hz under the premise that this forcing will ‘actively engage the mechanisms of cerebral autoregulation’. The delineation of the relative contributions of the three separate blockades was sought through the statistical method of ANCOVA decomposition, combined with projection pursuit regression, and was limited to the study of the steady-state autoregulatory curve at the LBNP frequency of 0.3 Hz—unlike our study that examines the resting-state dynamics over all frequencies. Their main finding was that the primary influence on the slope of the autoregulatory curve (at 0.03 Hz) within the range of active autoregulation is exerted by sympathetic activity. This finding is consistent with our PDM model interpretation that the third PDM of the ABP-to-CBFV relationship may be associated with sympathetic influences on cerebral autoregulation in low frequencies within the region of active autoregulation (resting state), because this third PDM has larger values in the low frequencies (less than 0.05 Hz).

The relative importance of the physiological mechanisms described by the PDMs is quantified by the values of the ANFs. This interpretation of the PDM-based model, if it becomes confirmed in future studies, may allow the development of model-based diagnostic tools for diseases with a significant cerebrovascular component (such as neurodegenerative diseases, stroke, diabetes, etc.) and enable the quantitative assessment of interventions or pharmacological treatments of such diseases.

The obtained ANFs in the PDM-based model are *subject-specific* (whereas the PDMs are common to the entire cohort) and indicate that there exist significant nonlinearities in the dynamic effects of the ETCO_{2} and HR inputs upon the output CBFV. It is seen in figure 2 that some ANFs of the ETCO_{2} and HR inputs are not symmetric about the origin, indicating the asymmetric effect of ETCO_{2} (hypercapnia versus hypocapnia) upon CBFV. No significant nonlinearities were observed for the dynamic effects of the ABP input within the narrow range of ABP variations in the resting-state data analysed in this study. It is surmised that these data remain within the plateau region of the celebrated pressure–flow homeostatic nonlinearity (i.e. much larger variations of ABP are needed in the data before the pressure–flow homeostatic nonlinearity can be captured by the model).

In figure 8, we illustrate the model-predicted CBFV response to concurrent pulse changes of HR and ABP, while keeping ETCO_{2} at baseline (figure 8*a*), and to concurrent pulse changes of ETCO_{2} and ABP, while keeping HR at baseline (figure 8*b*). The pulse strength is the half-average RMS value of the respective time-series data. The autoregulatory response to ABP changes is evident by the rapid return of CBFV to the steady state after the onset of the ABP input pulse, and it is symmetric when ETCO_{2} and HR are at baseline. However, the peak of the autoregulatory response is asymmetric when ETCO_{2} is elevated or HR is not at baseline. We also observe the anticipated increase in CBFV for increased ETCO_{2} (vasomotor reactivity) or increased HR, as well as the asymmetric decrease in CBFV for decreased ETCO_{2} or decreased HR (owing to their respective nonlinearities). The dynamic response of CBFV is faster to HR changes than to ETCO_{2} changes.

The presented data-based modelling approach holds the promise of revealing new and valuable insights into the dynamic characteristics of cerebral haemodynamics and autoregulation in future extended studies. Its ability to yield quantitative subject-specific measures of cerebrovascular characteristics engenders the prospect of improved diagnostic procedures for diseases with a significant cerebrovascular component. To the best of our knowledge, it is the first study to demonstrate that, in addition to the well-established effects of ABP spontaneous fluctuations on CBFV, there are direct cardiac effects on the latter, as shown by the contribution of HR. Therefore, it is a potentially valuable tool in advancing our understanding of heart–brain interactions.

## Author' contributions

V.Z.M. designed and directed the data analysis/modelling work and wrote the first draft of the paper. G.D.M. contributed significantly to the design of the analysis/modelling study and made critical contributions to the manuscript, including the preparation and submission of the final draft of the paper. D.C.S. performed most of the data analysis/processing and took care of the proper presentation of the results, including all figures and tables. R.Z. provided the experimental data and offered valuable physiological/scientific guidance in the design of the analysis/modelling study and the interpretation of the results.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported in part by the Biomedical Simulations Resource at the University of Southern California under NIH/NIBIB grant no. P41-EB001978 and NIA grant no. R01AG033106-01 to the UT-SWMC.

## Acknowledgements

We thankfully acknowledge Dr T. Tarumi for his contribution to the collection of the data. This work was supported in part by the Biomedical Simulations Resource at the University of Southern California under NIH/NIBIB grant P41-EB001978 and NIA R01AG033106-01 grant to the UT-SWMC. No competing interests or conflict of interest exist.

## 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.