## Abstract

We present a theoretical evaluation of a cardiovascular system identification method that we previously developed for the analysis of beat-to-beat fluctuations in noninvasively measured heart rate, arterial blood pressure, and instantaneous lung volume. The method provides a dynamical characterization of the important autonomic and mechanical mechanisms responsible for coupling the fluctuations (inverse modeling). To carry out the evaluation, we developed a computational model of the cardiovascular system capable of generating realistic beat-to-beat variability (forward modeling). We applied the method to data generated from the forward model and compared the resulting estimated dynamics with the actual dynamics of the forward model, which were either precisely known or easily determined. We found that the estimated dynamics corresponded to the actual dynamics and that this correspondence was robust to forward model uncertainty. We also demonstrated the sensitivity of the method in detecting small changes in parameters characterizing autonomic function in the forward model. These results provide confidence in the performance of the cardiovascular system identification method when applied to experimental data.

- beat-to-beat model
- mathematical modeling
- cardiovascular modeling
- autonomic nervous system
- hemodynamics

when one considers modeling the cardiovascular system, one usually envisions constructing a model based on physical principles that is capable of generating realistic data. This type of modeling approach, which we refer to as forward modeling, is useful for enhancing one's understanding of cardiovascular physiology. One may also consider an inverse modeling approach, in which models are built from measured data. This type of modeling approach is referred to as system identification when system dynamics or memory is being considered. System identification may potentially provide a powerful means for the intelligent patient monitoring of cardiovascular function. Rather than simply recording hemodynamic signals, the signals are mathematically analyzed so as to provide a dynamical characterization of the physiological mechanisms responsible for coupling them.

To this end, we (26, 27, 29) have previously developed a cardiovascular system identification method for the analysis of short-term, beat-to-beat fluctuations in heart rate signals derived from the surface electrocardiogram (ECG) and noninvasively measured arterial blood pressure (ABP) and instantaneous lung volume (ILV) signals. The method provides a dynamical characterization of the important autonomic and mechanical mechanisms that couple the measured fluctuations specifically in terms of impulse responses. The method also provides power spectra of perturbing noise sources, which represent the residual variability in each of the signals not attributed to the variability generated through the impulse responses.

We have evaluated the cardiovascular system identification method with experimental data collected during conditions of pharmacological autonomic blockade, postural changes, and diabetic autonomic neuropathy (27, 29). We found that these three conditions altered the impulse response and power spectra of perturbing noise source estimates in a manner consistent with known physiological mechanisms. This suggests, but does not directly demonstrate, the validity of the estimated system dynamics.

To evaluate directly the system dynamics estimated by the method, it is necessary to be able to establish the impulse responses and power spectra of perturbing noise sources in a manner independent of system identification. These impulse responses and power spectra of perturbing noise sources may then be regarded as the gold standards against which the corresponding estimates from system identification may be compared. Ideally, one would validate the method based on experimental data. However, the establishment of gold standard impulse responses and power spectra of perturbing noise sources would require extreme experimental conditions that would be virtually impossible to implement in practice. Consider, for example, the establishment of an impulse response, which would require, according to mathematical definition, applying an arbitrarily narrow, unit-area input to the appropriate point of the cardiovascular system and measuring the output response of interest while all other perturbations to the output are held constant. Moreover, even if this type of experiment could be implemented, it is unreasonable to expect that the cardiovascular state and system operating point would be precisely maintained, which renders the meaning of such a comparison to be questionable.

By contrast, the method could be readily evaluated on a theoretical basis with a forward model of the cardiovascular system, because the gold standard impulse responses and power spectra of perturbing noise sources would either be precisely known or easily determined for the relevant cardiovascular state and system operating point. That is, the gold standards could be readily established according to their mathematical definitions. A forward model would also provide a powerful means to analyze the sensitivity of the cardiovascular system identification method. That is, we would be able to determine precisely how much the dynamical properties of the forward model would have to be altered before we would see a corresponding change in the estimates. The major limitation of this type of evaluation is that the results would be only as meaningful as the extent to which the forward model coincides with the actual cardiovascular system. This limitation can be attenuated by determining the robustness of the estimates over a set of models that reflect the uncertainty in the relevant properties of the forward model. We note that the general concept of analyzing inverse modeling algorithms based on forward models of the cardiovascular system is not novel. For example, investigators (9, 16) have utilized complex forward models of the systemic circulation to evaluate the estimation of lumped parameters representing systemic arterial resistance [*R*
_{a}(*t*)] and systemic arterial compliance.

The principal goal of this study is to present a theoretical validation of our previously developed cardiovascular system identification method. To realize this aim, this study includes the following five components: *1*) a brief review of the cardiovascular system identification method, which has been previously described in detail (26, 27, 29); *2*) a description of a forward model of the human cardiovascular system; *3*) the procedure for evaluating the system identification method against the forward model, which includes the establishment of gold standards based on their mathematical definitions; *4*) the validation of the forward model in terms of power spectra of beat-to-beat variability; and *5*) the evaluation of the impulse response and power spectra of perturbing noise source estimates derived by applying system identification to beat-to-beat variability generated from the forward model against the corresponding gold standards in terms of accuracy, robustness, and sensitivity. Note that the final component amounts to assessing the equivalence of the impulse responses and power spectra of perturbing noise sources determined from executing the forward model under two conditions: *1*) resting conditions by applying system identification to analyze beat-to-beat variability, and*2*) extreme conditions established by the literal meaning of these quantities. Because the extreme conditions cannot be implemented in practice, this study seeks to determine whether the cardiovascular system identification method is able to estimate reliably the physiologically important impulse responses and power spectra of perturbing noise sources during experimental conditions that could be readily implementable.

## CARDIOVASCULAR SYSTEM IDENTIFICATION REVIEW

Figure 1 illustrates the model upon which the cardiovascular system identification method is based. The model includes five physiological coupling mechanisms relating ECG-derived heart rate signals, ABP, and ILV: circulatory mechanics, heart rate (HR) baroreflex, ILV→HR, ILV→ABP, and sinoatrial (SA) node.

Circulatory mechanics represents the mechanical feedforward effects of a pulsatile heart rate (PHR) signal [defined to be a train of unit-area impulses occurring at the times of ventricular contraction (R wave)] on the pulsatile ABP waveform. HR baroreflex represents the autonomically mediated feedback effects of fluctuations in ABP on fluctuations in an HR tachogram, which is derived from the R-R intervals (3). ILV→HR represents the autonomically mediated coupling between respiration and HR, which is responsible for governing respiratory sinus arrhythmia. ILV→ABP represents the mechanical effects of respiration on ABP due to the alterations in venous return and the filling of intrathoracic vessels and heart chambers associated with the changes in intrathoracic pressure. SA node maps HR (the net autonomic input to the sinoatrial node) to PHR (the onset times of each ventricular contraction) through an “integrate and fire” device referred to as the integral pulse frequency modulation (IPFM) model (3). Unlike the other four physiological coupling mechanisms, SA node is predefined and not identified from the measured signals.

The model also incorporates two perturbing noise sources, N_{HR} and N_{ABP}, which are determined from analysis of the measured signals. N_{HR} represents the fluctuations in HR not caused by fluctuations in ABP or ILV. Such fluctuations may result, for example, from autonomically mediated perturbations driven by cerebral activity. N_{ABP} represents fluctuations in ABP not caused by PHR or fluctuations in ILV. Such ABP fluctuations may result, for example, from fluctuations in*R*
_{a}(*t*) as tissue beds adjust local vascular resistance to match local blood flow to demand.

To obtain a complete characterization of each of the coupling mechanisms over their physiological range of frequencies, we employ a broadband excitation protocol, in which subjects breathe according to a random sequence of auditory tones (4). We have found that when the subjects are breathing according to this protocol and are otherwise at rest that the fluctuations are sufficiently small and stationary such that the coupling mechanisms may be characterized by linear time-invariant (LTI) transfer functions around a given system operating point (11). We specifically represent each of the coupling mechanisms with autoregressive moving average difference equations, a highly flexible subclass of LTI models, whose parameters may be conveniently identified with the analytic methods of linear least-squares estimation. Further details of the method are provided in (26, 29).

Figure 1 includes an example of the cardiovascular system identification results computed from a standing, healthy subject breathing randomly, in which ABP and ILV are noninvasively measured with a continuous blood pressure monitor (model 2300 Finapres; Ohmeda; Englewood, CO) and a Respitrace System two-belt chest-abdomen inductance plethysmograph (Ambulatory Monitoring; Ardsley, NY) (34). The identified transfer functions are presented in their time domain form of impulse responses. For example, the HR baroreflex impulse response estimate indicates that HR would immediately decrease and then return to baseline if a unit-area impulse of ABP were applied at *time 0* while the other inputs to HR (ILV and N_{HR}) were set to 0. The estimated perturbing noise sources in the figure are depicted in their frequency domain form of power spectra. For example, the power spectrum of N_{HR}indicates that ABP and ILV do not fully account for HR variability within ∼0.05 Hz.

## FORWARD MODEL DESCRIPTION

The forward model of the human cardiovascular system to which we apply the cardiovascular system identification method includes three major components: a pulsatile heart and circulation, a short-term regulatory system, and resting physiological perturbations. The pulsatile heart and circulation is a nonlinear, lumped parameter model (R. Mukkamala, D. A. Sherman, R. G. Mark, and R. J. Cohen, unpublished observations). The short-term regulatory system consists of three subcomponents: an arterial baroreflex, a cardiopulmonary baroreflex, and a direct neural coupling mechanism between respiration and HR. The resting physiological perturbations also include three subcomponents: respiratory activity, the autoregulation of local vascular beds, and 1/*f* HR fluctuations. We describe all of these components and subcomponents and include parameter values below. See appendix
for a description of forward model implementation.

### Pulsatile Heart and Circulation

The nonlinear, lumped parameter model of the pulsatile heart and circulation is illustrated in Fig. 2 in terms of its electrical circuit analog, in which charge is analogous to blood volume (Q, ml), current, to blood flow rate (q˙, ml/s), and voltage, to pressure (P, mmHg). The model consists of six compartments, which represent the left and right ventricles (l and r), systemic arteries and veins (a and v), and pulmonary arteries and veins (pa and pv). Each compartment consists of a conduit for viscous blood flow, which is characterized by either a linear or nonlinear resistance (*R*) and a volume storage element, which is characterized by either a linear or nonlinear compliance (C) with an associated unstressed volume (Q^{0}). The reference (ref) pressure is atmospheric pressure (or ground) for the systemic compartments and intrathoracic (th) pressure for the ventricular and pulmonary compartments. The compliances of the ventricles vary periodically over time (*t*) according to the model of Suga and colleagues (38, 39) and are responsible for driving the flow of blood. The four ideal diodes represent the ventricular inflow and outflow valves and ensure unidirectional blood flow. We have demonstrated (R. Mukkamala, D. A. Sherman, R. G. Mark, and R. J. Cohen, unpublished observations) that this model behaves reasonably in terms of pulsatile waveforms and limiting static terminal pressure-flow rate properties of the systemic circulation and heart-lung unit. However, by itself, the model cannot generate the short-term, low-frequency hemodynamic variability that the cardiovascular system elicits during resting conditions. To account for this variability, it is necessary to incorporate a short-term regulatory system and resting physiological perturbations with the model here.

### Short-Term Regulatory System

#### Arterial baroreflex.

The short-term regulatory system includes a baroreflex system, which is based on a previously developed model (13). This model, which is fully described here, conceptualizes the arterial (A) baroreflex arc as the feedback system illustrated in the block diagram of Fig. 3. The feedback system is aimed at tracking a setpoint (sp) pressure through the following sequence of events. The baroreceptors sense P_{a}(*t*) and then relay this pressure via autonomic afferent fibers to the autonomic centers in the brain. Here, the autonomic nervous system (ANS) compares the deviation between the sensed pressure and P
with zero and then sends commands via autonomic efferent fibers to adjust four parameters of the pulsatile heart and circulation such that the ensuing P_{a}(*t*) is near P
. The four adjustable parameters are instantaneous HR [*F*(*t*)], ventricular contractility [in terms of C
(*t*), the differential compliances that the ventricles assume at their unstressed volumes during end systole] (R. Mukkamala, unpublished communications), systemic venous unstressed volume [Q
(*t*)], and*R*
_{a}(*t*). The ANS specifically adjusts these parameters based on the history of P_{a}(*t*) − P
according to the following nonlinear, dynamical map
Equation 1
where *X* may represent any of the four controllable parameters, the atan function (parameterized by the constant*S*) imposes arterial baroreflex saturation, p(*t*) and s(*t*) are unit-area effector mechanisms, which, respectively, represent the fast, parasympathetic limb of the ANS and the slower, sympathetic limb (see Fig.4), and *G*
_{p} and*G*
_{s} reflect the respective weighting values of the effector mechanisms.

To map *F*(*t*) to the times of onset of ventricular contraction, an IPFM model of the sinoatrial node is also incorporated. This model integrates *F*(*t*) (in units of beats/s) over time until the integral is equal to one. At this point, systole is initiated through the ventricular variable compliance model (with the systolic duration determined from the preceding cardiac cycle length), the integral is set to zero, and the integration is repeated.

#### Cardiopulmonary baroreflex.

The baroreflex system also includes a cardiopulmonary (*CP*) baroreflex arc, which is conceptualized with a feedback diagram analogous to Fig. 3. However, the sensed pressure here is defined to be the effective right atrial transmural pressure [P
(*t*) = P_{“ra”}(*t*)−P_{th}(*t*)] of the pulsatile heart and circulation model (see Fig. 2), where P
is the effective right artrial transmural pressure, P_{ra} is right atrial pressure, and P_{th}is thoractic pressure.

#### Direct neural coupling mechanism between respiration and HR.

The short-term regulatory system model also includes a direct neural coupling mechanism between respiration and HR. This coupling mechanism is implemented in terms of an LTI impulse response, which maps fluctuations in ILV to fluctuations in *F*(*t*). In accordance with previous studies (29, 40) with experimental data, the impulse response is defined here by a linear combination of s(*t*) and p(*t*), each of which are advanced in time by 1.5 s. The model ILV signal [Q_{lu}(*t*)] that is necessary for the implementation of this mechanism is described below.

#### Parameter values.

The numerical values of the parameters characterizing the baroreflex system are taken from (13 and T. Heldt, E. B. Shim, R. D. Kamm, and R. G. Mark, unpublished observations) (see below for exceptions) and are provided in Table 1. Note that the cardiopulmonary baroreflex control of*F*(*t*) and*c*
(*t*) are both neglected in the model due to either controversial data or a lack of available data. Although the α- and β-sympathetic sublimbs are both characterized by s(*t*), the table specifies the particular sublimb corresponding to each s(*t*). This permits the modeling of, for example, the effects of propranolol, which selectively blocks the β-sympathetic sublimb.

The weighting values for the direct neural coupling mechanism are taken from (29) and are as follows (in ml): −0.0002 beats/s for s(*t*) and 0.00012 beats/s for p(*t*). These values indicate that upon inspiration, there is a withdrawal of parasympathetic activity, followed by a withdrawal of β-sympathetic activity.

### Resting Physiological Perturbations

#### Respiratory activity.

The resting physiological perturbations model includes both metronome (met) and random-interval respiratory activity represented in terms of Q_{lu}(*t*). Because ILV strongly resembles a sinusoid during metronomic breathing, Q_{lu}(*t*) is mathematically represented as follows
Equation 2where Q_{fr} is the functional residual volume of the lungs, Q_{t} is the tidal volume, and*T _{r}
* is the respiratory period. The model here may also be considered to represent ILV during normal, spontaneous breathing, which, at least to first approximation, appears to be sinusoidal as well.

In the random-interval breathing protocol, the period of each respiratory cycle is determined from the outcome of an independent probability experiment governed by a modified exponential probability density function, in which the respiratory period ranges from 1–15 s with a mean of 5 s (4). The subjects are allowed to control their tidal volume so that normal ventilation is not disrupted. Each respiratory cycle during random-interval breathing is therefore modeled to be one period of a sinusoid with the period determined from the experimental outcome of the probability experiment and the tidal volume chosen such that the alveolar (alv) ventilation rate is identical to that from metronome breathing. This amounts to the following mathematical representation for random-interval Q_{lu}(*t*)
Equation 3
where Q_{ds} is the dead space in the airways (air), the argument *j* denotes the *j*th respiratory cycle, whereas the subscript *j* denotes the commencement of the*j*th respiratory cycle, and
Equation 4with *x*(*j*) being the outcome of the*j*th independent probability experiment defined by a uniform distribution ranging from 0–1.

To account for the mechanical effects of respiratory activity on intrathoracic pressure, we incorporate the simple model of ventilation illustrated in Fig. 5 in terms of its electrical circuit analog. The electrical components may be interpreted similarly to those in Fig. 2 by considering air here rather than blood. Hence, the resistor may be thought of as a conduit for airflow between the atmosphere and the lungs whereas the capacitor may be interpreted as an air volume container representing the lung compartment, which is parameterized by an unstressed volume in addition to a linear compliance. The governing differential equations of this model provide a precise mapping from Q_{lu}(*t*) to P_{th}(*t*) as well as to P_{alv}(*t*), which influencesq˙_{pa}(*t*) (R. Mukkamala, D. A. Sherman, R. G. Mark and R. J. Cohen, unpublished observations).

#### Autoregulation of local vascular beds.

Although respiratory-induced hemodynamic fluctuations are fairly well understood, the mechanisms responsible for hemodynamic variability at frequencies below ∼0.1 Hz have not been adequately elucidated. It has been previously demonstrated that ABP fluctuations are not determined by HR fluctuations at these lower frequencies, and it has therefore been postulated that fluctuations in*R*
_{a}(*t*) caused by the autoregulation of local vascular beds are responsible for these ABP fluctuations, which, in turn, induces HR fluctuations through the baroreflex (1,24). However, little is known about the system dynamics characterizing autoregulatory processes except that these processes are relatively slow with a characteristic time constant from ∼5 to 20 s (6). We therefore include into the model of resting physiological perturbations a zeroth-order model of autoregulatory processes. This model represents these processes as a stochastic, exogenous disturbance to*R*
_{a}(*t*) [n_{Ra}(*t*)] defined by a Gaussian white noise process with zero mean and ς^{2}variance that is bandlimited to 0.1 Hz.

#### 1/f HR fluctuations.

Figure 1 demonstrates that ABP and ILV are not the only factors responsible for perturbing HR at frequencies within ∼0.05 Hz. However, the other perturbing factors, which may include cerebral activity impinging on the ANS, are essentially unknown. It is well known that HR fluctuations demonstrate fractal behavior over at least four decades of frequency, in which the highest frequency decade is within the frequency band of interest here (7, 37). We therefore include in the resting physiological perturbations model an exogenous, stochastic 1/*f* disturbance to*F*(*t*) [n_{F}(*t*)] The disturbance is created by passing a zero-mean, Gaussian white noise process with variance λ^{2} through a cascade combination of two LTI filters. One filter is of unit-area and is designed to have a 1/*f* magnitude squared frequency response over four decades of frequency starting at 10^{−4} Hz (12, 20,36). The other filter is defined by a linear combination of s(*t*) (β-sympathetic) and p(*t*) (with arbitrary weighting values of 12 beats/s each), because HR fluctuations are almost exclusively mediated by the ANS (2). Importantly, n_{F}(*t*) as well as n_{Ra}(*t*) are considered to be unmeasurable quantities in the forward model.

#### Parameter values.

The parameter values characterizing the models of Q_{lu}(*t*) and ventilation are provided in Table2. These values are taken from (21,41).

Both ς^{2} and λ^{2}, which characterize n_{Ra}(*t*) and n_{F}(*t*), respectively, are considered to be free parameters whose values are chosen such that the model low-frequency hemodynamic variability matched that determined from a set of experimental data previously obtained from 12 healthy standing humans breathing at a fixed rate of 0.25 Hz (10). The data set consists of simultaneous, noninvasive measurements of ECG and ABP (Finapres), as well as uncalibrated ILV derived from the ECG for each subject. The specific procedure for choosing the parameter values based on this data set is as follows. We computed the mean ± SD of the power in three nonoverlapping frequency bands for HR and ABP over the group of 12 subjects (see Table 3) by using autoregressive spectral analysis (22). The free parameters ς^{2} and λ^{2} were then tuned during fixed-rate breathing excitation at 0.25 Hz, such that the power in these frequency bands for *F*(*t*) and P_{a}(*t*) were near their respective human values. To satisfy sufficiently this matching procedure, we also had to increase the weighting values of the parasympathetic and β-sympathetic effector mechanisms by 33% from their original values given in (13). However, this increase did not move these values outside of their physiological range, as established by experimental data (14). On the basis of this procedure, we set ς = 0.16 mmHg · s · ml^{−1} and λ = 0.0035 beats/s.

On the basis of this procedure, the hemodynamic variability generated by the forward model may be thought of as representing an “average” subject. We did not tune the model to represent an individual subject because there is tremendous intersubject variability in experimental hemodynamic spectra. If we had chosen to represent an individual subject, the results of this study would be biased toward the chosen subject.

## FORWARD MODEL-BASED EVALUATION PROCEDURE

### Gold Standard Cardiovascular System Identification Results

To evaluate the performance of the cardiovascular system identification method applied to data generated from the forward model, it is necessary to establish, in a manner independent of system identification, the impulse responses and power spectra of perturbing noise sources in Fig. 1, which characterize the forward model. These impulse responses and power spectra of perturbing noise sources may then be regarded as the gold standard cardiovascular system identification results against which the estimates may be compared.

We establish the gold standards according to the literal mathematical meaning of the impulse responses and the power spectra of the perturbing noise sources. In particular, the impulse response is defined to represent the output response of a system to an arbitrarily narrow, unit-area input. The impulse response also completely characterizes the input-output relationship of the physiological coupling mechanisms here because of our LTI assumption. Hence, in establishing each of the gold standard impulse responses, we apply an impulse input to the desired point in the forward model and measure the output response of interest while holding all other perturbations to the output constant. For example, to establish the gold standard ILV→ABP impulse response, we measure the ABP response to an impulse input of ILV while holding HR constant and setting the exogenous disturbance to *R*
_{a}(*t*) to zero.

By definition, the perturbing noise source represents the residual variability in a measured output not due to the fluctuations in the other measured inputs. So, in establishing each of the gold standard power spectra of perturbing noise sources, we first measure the output of interest in the forward model while holding the two measured input fluctuations constant and then compute the power spectrum of the output. For example, to establish the gold standard power spectrum of N_{ABP}, we first measure ABP in the forward model while holding HR and ILV constant and then compute the power spectrum of the resulting ABP.

The specific procedures for establishing each of the gold standard impulse responses and power spectra of perturbing noise sources in Fig.1 are described in appendix . This description includes the signal processing necessary for the estimated and gold standard impulse responses to be sampled identically.

### Root-Normalized-Mean-Squared Error

To provide a compact means for evaluating the similarity of each of the estimates with respect to its corresponding gold standard, we utilize a statistic, which we refer to as the root-normalized-mean-squared error (RNMSE). Let the estimate be vector **x** with mean vector **x̂** and covariance matrix **Λ _{x}
** (which reflects the uncertainty in the estimate) and its corresponding gold standard be vector

**x**

_{0}. The RNMSE is then defined as follows Equation 5 where

*E*( · ) is the expectation operator and trace( · ) is the operator, which sums the diagonal elements of its matrix argument. This statistic may be thought of as representing the percentage of error of the estimate with respect to its gold standard. In the case of the power spectral estimates, in which no measure of uncertainty is available, the matrix

**Λ**is assumed to be a matrix of zeros.

_{x}### Monte Carlo Simulations

Although the impulse response estimates include a measure of uncertainty in terms of a covariance matrix, this uncertainty measure is only an estimate itself (31). To account more accurately for estimation error variance, we report the mean ± SD for the impulse response estimates as well as for the power spectral estimates and RNMSE statistics as determined from 20 different realizations of forward model generated data.

## FORWARD MODEL VALIDATION

The forward model-based evaluation of the cardiovascular system identification method is only meaningful provided that the low-frequency power spectral content of the signals analyzed by the method is realistic. Because the forward model is tuned to represent what one may think of as an “average” subject, we initially considered demonstrating that the model spectra resemble the group average spectra of the 12 healthy humans breathing at a fixed rate of 0.25 Hz (see *Parameter Values* in *Resting Physiological Perturbations*). However, averaging tended to smear out the spectral peaks, which were usually not centered at the same frequencies for each subject. We therefore demonstrate that the model spectra resemble what one may find experimentally through a comparison with the spectra from one of the healthy subjects. This does not seem unreasonable if we keep in mind that the power in three nonoverlapping frequency bands of the model spectra quantitatively match that computed from the group average of the 12 subjects (see Table 3).

Figure 6 illustrates that the power spectra for ILV, HR, and ABP at frequencies below the mean HR obtained from the forward model during fixed-rate breathing excitation at 0.25 Hz indeed resemble spectra obtained from one of the 12 subjects. Differences at the respiratory frequency in Fig. 6 and in Table 3 may be attributed to the fact that the subject here did not precisely follow the fixed-rate breathing protocol as well as discrepancies in tidal volume, which was not measured. Note that although the forward model was tuned to represent the “average” subject, it could have been tuned to match more precisely the subject here.

Figure 6 also demonstrates that the model spectra exhibit a spectral peak at ∼0.07 Hz, which is near the center frequency of the “posture” peak in humans (32). Because the exogenous disturbances within 0.1 Hz are broadband (see *Autoregulation of Local Vascular Beds* and *1/f HR Fluctuations*), this peak implies a system resonance in the forward model. This result supports the simple computer simulation of de Boer (14), which demonstrated that the “posture” peak could be due to a system resonance. de Boer specifically implicated the system resonance to the arterial*R*
_{a}(*t*) baroreflex arc. However, based on a few simple experiments, in which the weighting values of the open-loop effector mechanisms were varied, we have found that the resonance peak is substantially diminished in the absence of the arterial baroreflex feedback pathways responsible for controlling Q
(*t*) as well as*R*
_{a}(*t*). That is, the arterial baroreflex controlling α-sympathetic activity in general is involved in eliciting the resonance peak.

It is also necessary to demonstrate that the cross spectra between each of these signals are realistic. This essentially amounts to demonstrating that the transfer functions, which characterize the couplings between the signals are consistent with those determined from experimental data. However, this has been virtually taken care of by the very construction of the model, which was based largely on physiological findings published in the literature.

Although it is only necessary for the purposes of this study that the low-frequency spectral content of ILV, HR, and ABP be reasonable, we assume that the low-frequency content of the remaining model signals are accounted for as well by virtue of reasonably representing these three signals. Figure 7 gives some credence to this assumption by illustrating that the model spectrum of left ventricular flow rate (LVFR) at frequencies below the mean HR resembles the corresponding spectrum measured with a Doppler ultrasound technique (17) from an individual subject breathing at a fixed rate of 0.25 Hz and tilted upright (30° with respect to the supine posture).

## FORWARD MODEL-BASED EVALUATION

### Accuracy Analysis

We applied the cardiovascular system identification method to the beat-to-beat variability generated from the forward model during random-interval breathing excitation. Figure8 illustrates that there is a close agreement between the resulting cardiovascular system identification estimates and the corresponding gold standards. These results indicate the validity of the cardiovascular system identification method to the extent that the forward model coincides with the actual cardiovascular system.

### Robustness Analysis

Perhaps the major source of uncertainty in the relevant properties of the forward model involves the system dynamics responsible for controlling HR. We therefore consider here the robustness of the HR baroreflex, ILV→HR, and N_{HR} estimates over a set of forward models, which reflect some of our uncertainty in these system dynamics.

#### Cardiopulmonary HR baroreflex.

We first consider how these estimates would be altered if a cardiopulmonary HR baroreflex were included in the forward model. Inclusion of this reflex may be achieved simply by adjusting its weighting values to nonzero values (see Table 1). To proceed with this analysis, the gold standard N_{HR} power spectrum, which is no longer valid in the presence of a cardiopulmonary HR baroreflex, must be redefined (see appendix for the specific procedure).

Figure 9, *A*–*C*, illustrates the RNMSE results of the HR baroreflex, ILV→HR, and N_{HR} estimates as a function of the ratio of the static gains of the cardiopulmonary to arterial HR baroreflexes. The range of this ratio is determined from published experimental data (15,35). Note that a negative ratio indicates a Bainbridge type of cardiopulmonary baroreflex, whereas a positive ratio indicates a type of cardiopulmonary baroreflex, which contributes to ABP regulation. The results indicate that the ILV→HR impulse response estimate is most significantly and adversely influenced by the presence of the cardiopulmonary HR baroreflex. This result implies that low-frequency effective right atrial transmural pressure fluctuations are determined largely by ILV fluctuations.

From Fig. 9
*A*, we see that the HR baroreflex impulse response estimate is not affected by a positive ratio of the static gains of the cardiopulmonary to arterial HR baroreflexes. However, for a negative ratio, the estimate is adversely affected but not as significantly as the ILV→HR impulse response estimate (consider RNMSE standard deviation). We hypothesize that this discrepancy between positive and negative ratios is due to the size of low-frequency ABP fluctuations. That is, in contrast to the positive ratio, which results in tighter ABP regulation, the negative ratio leads to sufficiently large low-frequency ABP fluctuations such that they contribute somewhat to the generation of low-frequency effective right atrial transmural pressure fluctuations in the model.

Figure 9
*B* illustrates that the ILV→HR estimate becomes progressively worse with respect to its corresponding gold standard with increasing absolute static gain of the cardiopulmonary HR baroreflex. However, this does not necessarily imply that the estimate is meaningless for large absolute static gain values. In fact, the estimate is quite meaningful regardless of the static gain value of the cardiopulmonary HR baroreflex provided that we redefine the gold standard such that it encompasses the dynamics of both the direct neural coupling mechanism and the cardiopulmonary HR baroreflex (seeappendix
for the specific procedure). Figure9
*D* shows that the RNMSE results for the ILV→HR estimate with respect to the new gold standard is essentially independent of the ratio of the static gain values of the cardiopulmonary to arterial HR baroreflexes. The results here suggest that when considering the ILV→HR impulse response identified from experimental data, it is probably more accurate to interpret the estimate analogously to the gold standard as redefined here.

#### Arterial baroreflex saturation.

The *S* parameter in *Eq. 1
* characterizes the degree of arterial baroreflex saturation in the forward model. This parameter provides an upper bound on the deviation between sensed ABP and its setpoint pressure so that the manipulated variables cannot be controlled to arbitrarily high or low values. The maximum deviation permitted by this parameter in the forward model is ∼28 mmHg. This pressure is approximately equal to the range of ABP, over which the atan mapping in *Eq. 1
* differs from true linearity by ∼10% (that is, the linear ABP range). However, experimental data from one study (25) indicate that the linear ABP range may actually be only ∼10 mmHg. We therefore consider the robustness of the HR baroreflex, ILV→HR, and N_{HR} estimates against more narrow linear ABP ranges.

Figure 10 shows the RNMSE results of the HR baroreflex, ILV→HR, and N_{HR} estimates as a function of the linear ABP range of the arterial baroreflex. These results indicate that only the HR baroreflex estimate is significantly and adversely affected by increasing degrees of saturation. Furthermore, the deviation of the RNMSE from its nominal value only becomes significant when the extent of saturation is no longer substantiated by the experimental data in (25). Hence, provided that arterial baroreflex saturation is the only significant nonlinearity, then the assumption that the fluctuations in cardiovascular system identification data are small enough such that the couplings between the fluctuations are related linearly seems quite tenable.

### Sensitivity Analysis

The ultimate potential of the cardiovascular system identification method is to provide a clinician with an intelligent, sensitive tool for monitoring a patient's cardiovascular state over time so as to guide therapy. An analysis of the resolving power of the method is germane to the realization of this potential. We therefore consider here an analysis of the sensitivity of the method particularly in terms of detecting changes in autonomic activity, as reflected by the weighting values for p(*t*) and β-sympathetic s(*t*) in the forward model, through the HR baroreflex, ILV→HR, and N_{HR} estimates.

Figure 11, *A* and*B*, shows the sensitivity results to changes in autonomic function in terms of the estimated versus actual percentage of change in the static gains (sum of the weighting values) of the autonomically mediated impulse responses, where the percentage of change here is with respect to the nominal static gain values. These results indicate that a change of 25% in the static gain of the gold standard HR baroreflex impulse response is required to detect reliably a change in the static gain of the corresponding estimate. On the other hand, a 25% change in the static gain of the gold standard ILV→HR impulse response is not sufficient to detect reliably a change in the static gain of the corresponding estimate.

On the basis of our experience with experimental human data (27,29), we have found that the absolute peak amplitude of the impulse response estimates seems to be quite reliable in detecting changes in autonomic function. Figure 11, *C* and*D*, shows the sensitivity results in terms of the estimated versus actual percentage of change in the absolute peak amplitude of the autonomically mediated impulse responses. Because we varied the static gains of the impulse responses simply via scaling, the actual percentage of change in the static gain of the impulse responses is equal to the actual percentage of change in the absolute peak amplitude. The results indicate substantial improvement in sensitivity with respect to the ILV→HR estimate and some improvement with respect to the HR baroreflex estimate. That is, as small as a 10% change in the absolute peak amplitude of the gold standard ILV→HR impulse response is sufficient to detect reliably a change in the absolute peak amplitude of the corresponding estimate.

The substantial improvement in the sensitivity of the ILV→HR impulse response may be explained by considering the ratio of the spectra of ILV to the spectra of the unmeasured HR disturbance (signal-to-noise ratio as a function of frequency) in the model. This ratio is smallest at frequencies near 0 Hz because of the 1/*f* character of the unmeasured disturbance as well as the reduced ILV spectral content near 0 Hz due to the probability density defining the random-interval respiratory pattern (see *Respiratory Activity*). The static gain of the ILV→HR impulse response estimate, which reflects the 0 Hz estimate, is therefore not reliable. However, the ratio is larger at higher frequencies, and the absolute peak amplitude, which reflects the wider bandwidth parasympathetic filter, is therefore quite reliable. On the other hand, the spectral content of ABP is sufficiently significant with respect to that of the unmeasured 1/*f* HR fluctuations at frequencies near 0 Hz such that the static gain of the HR baroreflex impulse response is fairly reliable.

Figure 12 illustrates the sensitivity results in terms of the estimated and actual percentage of change in total, low-frequency (0–0.15 Hz), and high-frequency (0.15–0.4 Hz) power of the N_{HR} spectra. These results emphasize the superior reliability of the high-frequency components of the estimates due to the 1/*f* character of the unmeasured HR disturbance. Note that the high-frequency power of N_{HR} is just as sensitive a measure of parasympathetic function as the absolute peak amplitude of ILV→HR. However, we believe that this absolute peak amplitude may be a better measure of parasympathetic function, because the perturbing noise source N_{HR}, unlike ILV→HR, is not normalized for inputs, e.g., cerebral activity.

## DISCUSSION

### Previous Forward Models

A few forward models of the cardiovascular system have been previously developed for the purposes of eliciting and analyzing low-frequency hemodynamic variability. de Boer (14) presented a beat-to-beat model, which incorporated windkessel and Starling properties with arterial baroreflex control of HR and systemic arterial resistance. The model was perturbed with sinusoidal respiration and stochastic disturbances to HR and pulse pressure to generate resting hemodynamic fluctuations. Madwed et al. (23) developed a model to study the generation of Mayer waves. Their model included windkessel properties, constant stroke volume, arterial baroreflex control of HR and systemic arterial resistance, and sinusoidal respiration. Cavalcanti (8) presented a model that incorporated windkessel properties with nonlinear arterial baroreflex control of HR and assumed stroke volume and systemic arterial resistance to be constant. Despite the absence of exogenous perturbations, the model was demonstrated to elicit HR fluctuations due to the nonlinear delayed structure of closed-loop control.

None of these simple models are nearly as comprehensive as our forward model in terms of accounting for pulsatility, the short-term regulatory system, and resting physiological perturbations. Our forward model thus provides a more complete means for studying the mechanisms responsible for eliciting beat-to-beat variability compared with these oversimplified models. Consider, for example, our more thorough analysis of the physiological mechanisms responsible for the system resonance phenomenon with respect to the study of de Boer (14) (see forward model validation).

### Forward Model Limitations

The mechanisms responsible for eliciting low-frequency hemodynamic variability, particularly at frequencies below ∼0.1 Hz, have not been fully elucidated. Nevertheless, we constructed a forward model of the human cardiovascular system that is capable of generating reasonable low-frequency hemodynamic variability and is largely consistent with previous experimental findings. However, the forward model presented here is not the unique model for simultaneously realizing this behavior and being consistent with experimental findings. That is, we could alter the properties of the forward model and still generate realistic low-frequency hemodynamic variability. For example, the resonance peak at ∼0.07 Hz could also be elicited with a stochastic, exogenous disturbance to Q
(*t*), which may reflect, for example, fast-acting hormonal loops (19). In the forward model, this disturbance would induce ABP fluctuations, and consequently HR fluctuations, through stroke volume fluctuations. This is consistent with experimental findings that have hypothesized systemic arterial resistance fluctuations, but have not precluded stroke volume fluctuations, in the genesis of low-frequency ABP and HR fluctuations (1, 24). The forward model as a general beat-to-beat variability model should therefore be considered with caution. Importantly, however, the forward model as a means for evaluating our cardiovascular system identification method is more reliable by virtue of the realistic autospectra of ILV, HR, and ABP at frequencies below the mean HR (Fig. 6) and the realistic cross spectra between these signals (Fig. 8).

### Robustness to Other Forward Model Uncertainties

Although we evaluated how the presence of a cardiopulmonary HR baroreflex and varying degrees of arterial baroreflex saturation would affect the autonomically mediated estimates, we did not address the effects of nonstationarities or other types of nonlinearities on these estimates. However, the stationarity property of the forward model seems quite tenable given that the data are considered over short time periods during stable experimental conditions. The stationarity assumption is further supported by preliminary studies that we have conducted, in which the impulse response estimates do not change much from adjacent time periods of data collection. Although hysteresis and nonlinear interaction between the arterial and cardiopulmonary baroreflexes have been previously reported (25, 30), such complex behaviors are usually elicited during extreme experimental conditions. Furthermore, nonlinear models have been shown to provide only a modest improvement with respect to linear models in accounting for HR fluctuations measured during the relatively stable experimental conditions of system identification data collection (11). Hence, the inclusion of only the ubiquitously reported arterial baroreflex saturation as a regulatory system nonlinearity in the forward model seems quite reasonable. However, we acknowledge that if significant system nonstationarities and/or nonlinearities excluding baroreflex saturation were present, then the autonomically mediated estimates could be substantially altered with respect to their gold standards.

We also did not include an analysis of the effects of the size of the 1/*f* HR fluctuations on the estimates. This disturbance may be thought of as the unmeasured noise in the identification of the autonomically mediated impulse responses, HR baroreflex and ILV→HR. The relative contribution of this noise term with respect to HR fluctuations due only to ABP and ILV fluctuations reflects the signal-to-noise ratio of this identification problem. This relative contribution, specifically in terms of the ratio of the standard deviation of those HR fluctuations not attributed to linear couplings with ABP and ILV to the standard deviation of all HR fluctuations, has been reported to be 0.31 ± 0.14 based on experimental data obtained from humans breathing randomly (11). However, the ratio of the standard deviation of the 1/*f* HR fluctuations to the standard deviation of HR fluctuations in the forward model is 0.49 ± 0.07. In fact, we have found that none of the autonomically mediated estimates in Fig. 8 are significantly altered for ratios up to ∼0.6.

We also did not analyze the robustness of the mechanically mediated estimates (circulatory mechanics, ILV→ABP, and N_{ABP}) to uncertainties in the pulsatile heart and circulation model (i.e., nonlinearities in systemic arterial compliance). However, we believe that this is not too critical because linear lumped parameter models are quite reasonable when small, low-frequency hemodynamic fluctuations are being analyzed.

### Sensitivity Analysis Results Supported By Experimental Data Studies

From experimental human data (27, 29), we have found the absolute peak amplitude of ILV→HR to be the most sensitive measure of autonomic function, followed by the power in the N_{HR} spectra (especially the high-frequency power), the absolute peak amplitude of HR baroreflex, and the static gain of HR baroreflex. In contrast, we have never found the static gain of ILV→HR to be predictive of changes in autonomic function. The fact that the more precise sensitivity analysis here (see *Sensitivity Analysis*) retrospectively predicts the resolving power of the estimates determined from experimental human data helps confirm the validity of the low-frequency hemodynamic fluctuations generated by the forward model during random-interval breathing excitation.

In addition to the static gain and absolute peak amplitude of the impulse response estimates, we have considered another parameter for detecting autonomic changes in experimental human data, which we refer to as the characteristic time parameter and define as follows
Equation 6where *h*(*t*) is the estimated impulse response (27, 29). Changes in this parameter essentially reflect shifts in balance between parasympathetic and β-sympathetic function. For example, simply scaling the impulse response, as we have done (see*Sensitivity Analysis*), does not change the characteristic time because there is no shift in balance in the activity of the two autonomic limbs. On the basis of our sensitivity analysis results, we expect the ILV→HR impulse response estimate to be unable to measure true shifts in balance through the characteristic time parameter, because only the parasympathetic limb can be estimated accurately. However, we do expect the characteristic time parameter of the HR baroreflex impulse response estimate to be a somewhat sensitive measure of true shifts in balance between parasympathetic and β-sympathetic function. The validity of these expectations is supported by our analysis with experimental human data (27, 29), in which only the characteristic time parameter of the HR baroreflex impulse response estimate was found to be sensitive to changes in autonomic function.

In conclusion, based on the results of the forward model-based evaluation of the cardiovascular system identification method, we make the following inferences about the performance of the method with respect to experimental data. First, each of the cardiovascular system identification estimates is likely to reflect the system dynamics of actual physiological mechanisms. Second, the ILV→HR impulse response estimate encompasses both direct neural coupling and cardiopulmonary HR baroreflex mechanisms. Third, arterial baroreflex saturation and the relative contributions of the HR fluctuations independent of ILV and ABP are unlikely to affect significantly the autonomically mediated estimates. Fourth, the absolute peak amplitude of the ILV→HR impulse response estimate is a very sensitive measure of parasympathetic function. Finally, the HR baroreflex impulse response estimate is a reasonably sensitive measure of both parasympathetic function (through absolute peak amplitude and static gain) and β-sympathetic function (through static gain) and consequently shifts in balance between the two autonomic limbs (through characteristic time).

The theoretical evaluation therefore provides confidence in the performance of the method with respect to experimental data and demonstrates the power of the cardiovascular system identification method. That is, the method is able to estimate reliably physiologically important impulse responses and perturbing noise sources by analyzing beat-to-beat hemodynamic variability obtained during near normal conditions rather than extreme experimental conditions, which would be required by the literal mathematical meaning of these quantities. This study supports system identification as a powerful approach for the intelligent patient monitoring of cardiovascular function. Moreover, the forward model, on which the theoretical validation is based, provides a convenient test bed of data, which may facilitate the development of new methods that could be incorporated with the cardiovascular system identification method so as to provide a more detailed picture of cardiovascular state (26).

## Acknowledgments

The authors thank Dr. Janice Meck, National Aeronautics and Space Administration (NASA) Johnson Space Center, for providing the experimental data utilized in this work, and Roger G. Mark for useful comments.

## Appendix

### Forward Model Implementation

The model of the pulsatile heart and circulation may be mathematically represented by a set of coupled, first-order differential equations, which are obtained with Kirchoff's Current Law and the constitutive relationship of each of the elements (R. Mukkamala, D. A. Sherman, R. G. Mark, and R. J. Cohen, unpublished observations). Given an initial set of pressures, the ensuing pressures, volumes, and flow rates are determined by numerically integrating this set of equations with a fifth-order Runge-Kutta method with adaptive step size (33). The average time step is ∼0.008 s; however, the models of the short-term regulatory system and resting physiological perturbations are bandlimited to frequencies within the mean HR. For the purposes of reducing unnecessary computations, these models are implemented at a uniform sampling period of 0.0625 s. To implement simultaneously the entire forward model, it is therefore necessary to transform nonuniformly sampled signals to uniformly sampled signals and vice versa. The former transformation is achieved by averaging over the most recent 0.25 s every 0.0625 s, whereas the latter transformation is accomplished with linear interpolation.

## Appendix

### Specific Procedure For Establishing Gold Standards

The gold standard circulatory mechanics impulse response, which represents the ABP waveform that results from a single ventricular contraction, is established through the superposition principle. That is, the forward model is first executed for *n*, and then*n* − 1, ventricular contractions while n_{Ra}(*t*) is set to 0 and Q_{lu}(*t*) and *F*(*t*) are held constant. The two resulting nonuniformly sampled P_{a}(*t*) signals are resampled to 90 Hz, and their difference is defined to be the gold standard.

To establish the gold standard HR baroreflex impulse response, the feedback pathway of the arterial HR baroreflex arc, which is defined by*Eq. 1
* and Table 1, is removed from closed-loop operation. The input applied to this isolated feedback pathway resembles an impulse and is mathematically expressed as follows
Equation B1where *ς*
_{ABP}, the area of the “impulse,” is set to the standard deviation of the ABP fluctuations analyzed by the cardiovascular system identification method. The resulting *F*(*t*) is then resampled to 1.5 Hz and normalized by *ς*
_{ABP} to arrive at the gold standard.

The gold standard ILV→HR impulse response is predefined by virtue of its forward model implementation (see *Direct Neural Coupling Mechanism Between Respiration and HR*). It is only necessary to resample this impulse response to 1.5 Hz to establish the gold standard.

The gold standard ILV→ABP impulse response is determined by applying an impulse of Q_{lu}(*t*) to the forward model while n_{Ra}(*t*) is set to 0 and*F*(*t*) is held constant. The precise input applied here is analogous to *Eq. EB1
* with the area of the “impulse” set to the standard deviation of the ILV fluctuations analyzed by the cardiovascular system identification method (ς_{ILV}). The resulting P_{a}(*t*) is resampled to 1.5 Hz with its mean removed and normalized by ς_{ILV} to arrive at the gold standard.

The gold standard N_{HR} perturbing noise source is precisely given by n_{F}(*t*). The power spectrum of this signal may therefore be given by the product of λ^{2} and the magnitude squared frequency response of the 1/*f* filter and the filter defined by the linear combination of s(*t*) and p(*t*) (see *1/f HR Fluctuations*).

The gold standard N_{HR} perturbing noise source is determined by executing the forward model while *F*(*t*) and Q_{lu}(*t*) are held constant and only n_{Ra}(*t*) is active. The resulting P_{a}(*t*) is resampled to 1.5 Hz with its mean removed, and then the power spectrum of this signal is computed by autoregressive spectral analysis. This procedure is repeated for 10 different realizations of forward model data, with the average power spectrum defined to be the gold standard.

## Appendix

### Specific Procedure For Redefining Gold Standards With a Cardiopulmonary HR Baroreflex

The gold standard N_{HR} perturbing noise source in the presence of a cardiopulmonary HR baroreflex is established by first executing the forward model, in which P_{a}(*t*) is fixed to P
by replacing the systemic arterial capacitor in Fig. 2 with a constant voltage source, while Q_{lu}(*t*) is held constant. The resulting*F*(*t*) is resampled to 1.5 Hz with its mean removed, and the power spectrum of this signal is computed analogously to the gold standard N_{ABP} power spectrum (seeappendix
) to establish the gold standard.

The gold standard ILV→HR impulse response in the presence of a cardiopulmonary HR baroreflex is established by applying an impulse of Q_{lu}(*t*) (see *Eq. EB1
*) to the forward model with a voltage source in lieu of the systemic arterial capacitor as described above, while *n _{F}
*(

*t*) and

*n*

_{Ra}(

*t*) are set to zero. The resulting

*F*(

*t*) is resampled to 1.5 Hz with its mean removed and normalized by ς

_{ILV}to arrive at the newly defined gold standard.

## Footnotes

This work was supported by NASA Grant NAG5-4989, a National Space Biomedical Research Institute grant, and National Institutes of Health Grant P41-RR-13622.

Address for reprint requests and other correspondence: R. Mukkamala, MIT, E25-505, 45 Carleton St., Cambridge, MA 02139 (E-mail: rmukkama{at}mit.edu).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2001 the American Physiological Society