Harvard-Massachusetts Institute of Technology Division of Health
Sciences and Technology, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139
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
 |
INTRODUCTION |
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 [Ra(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.

View larger version (21K):
[in this window]
[in a new window]
|
Fig. 1.
Mean (solid traces) and standard deviation (dotted
traces) of cardiovascular system identification estimates obtained
from a standing human subject breathing randomly. Figure adapted from
raw data described in Ref. 34.
|
|
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,
NHR and NABP, which are determined from
analysis of the measured signals. NHR 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. NABP represents
fluctuations in ABP not caused by PHR or fluctuations in ILV. Such ABP
fluctuations may result, for example, from fluctuations in
Ra(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 NHR) 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 NHR
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 A 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 (
, 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 (Q0). 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.

View larger version (12K):
[in this window]
[in a new window]
|
Fig. 2.
Electrical circuit analog of the human pulsatile heart and
circulation model. Here, current is analogous to blood flow rate
( , ml/s), whereas voltage is analogous to pressure (P, mmHg).
The model includes 6 compartments representing 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 linear or nonlinear resistance (R) and a linear or
nonlinear and time (t)-varying compliance (C) with an
associated unstressed volume. Note that each nonlinear element is
denoted with a box. The reference pressure is intrathoracic (th)
pressure for the ventricular and pulmonary compartments. The 4 ideal
diodes represent the ventricular inflow and outflow valves.
|
|
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 Pa(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 Pa(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
Ra(t). The ANS specifically adjusts
these parameters based on the history of
Pa(t)
P
according to the following nonlinear, dynamical map
|
(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 Gp and
Gs reflect the respective weighting values of
the effector mechanisms.

View larger version (7K):
[in this window]
[in a new window]
|
Fig. 3.
Block diagram of the feedback system depicting the arterial
baroreflex arc. The feedback system is aimed at tracking a setpoint
(sp) pressure through the following sequence of events. The
baroreceptors sense Pa(t). The
autonomic nervous system (ANS) compares the deviation between the
sensed pressure and P with 0 and then responds by
adjusting 4 parameters of the pulsatile heart and circulation to keep
the resulting Pa(t) near P .
The four adjustable parameters are instantaneous heart rate
[F(t)], ventricular contractility
[C (t)], systemic venous unstressed
volume [Q (t)], and systemic arterial
resistance [Ra(t)]. The ANS
controls these parameters through a static saturation mapping, followed
by convolution with the impulse responses in Fig. 4 (see Eq. 1).
|
|

View larger version (11K):
[in this window]
[in a new window]
|
Fig. 4.
Unit-area effector mechanisms representing the fast,
parasympathetic limb p(t) (A) and slower,
sympathetic limb s(t) (B). These effector
mechanisms characterize the dynamical properties of the block labeled
ANS in Fig. 3. The signals p(t) and s(t) are
derived from experimental data in Ref. 5 (figure adapted from Ref.
13).
|
|
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)
Pth(t)] of
the pulsatile heart and circulation model (see Fig. 2), where
P
is the effective right artrial transmural
pressure, Pra is right atrial pressure, and Pth
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
[Qlu(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
Qlu(t). Because ILV strongly resembles a
sinusoid during metronomic breathing, Qlu(t) is
mathematically represented as follows
|
(2)
|
where Qfr is the functional residual volume of the
lungs, Qt is the tidal volume, and
Tr 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 Qlu(t)
|
(3)
|
where Qds is the dead space in the airways (air), the
argument j denotes the jth respiratory cycle,
whereas the subscript j denotes the commencement of the
jth respiratory cycle, and
|
(4)
|
with x(j) being the outcome of the
jth 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 Qlu(t) to
Pth(t) as well as to
Palv(t), which influences
pa(t) (R. Mukkamala, D. A. Sherman,
R. G. Mark and R. J. Cohen, unpublished observations).

View larger version (7K):
[in this window]
[in a new window]
|
Fig. 5.
Model of ventilatory mechanics in terms of its electrical
circuit analog, in which charge is analogous to air volume (Q) and
voltage is analogous to air pressure (P). The resistor (R)
represents the airway (air) between the atmosphere (atm) and the lungs
(lu), whereas the capacitor (C) represents the lung
compartment. Pth(t) and
Palv(t), intrathoracic and alveolar pressures,
respectively.
|
|
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
Ra(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
Ra(t)
[nRa(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) [nF(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,
nF(t) as well as
nRa(t) are considered to
be unmeasurable quantities in the forward model.
Parameter values.
The parameter values characterizing the models of
Qlu(t) and ventilation are provided in Table
2. These values are taken from (21,
41).
Both
2 and
2, which characterize
nRa(t) and
nF(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 Pa(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.
View this table:
[in this window]
[in a new window]
|
Table 3.
Power in three low-frequency bands of HR and ABP as determined from the
forward model and experimental human data during fixed-rate breathing
at 0.25 Hz
|
|
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 Ra(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
NABP, 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 B. 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
and covariance
matrix
x (which reflects the uncertainty in
the estimate) and its corresponding gold standard be vector
x0. The RNMSE is then defined as follows
|
(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
x is assumed to be a matrix of zeros.
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.

View larger version (19K):
[in this window]
[in a new window]
|
Fig. 6.
Power spectra of instantaneous lung volume (ILV;
normalized to unity power), HR, and arterial blood pressure (ABP) at
frequencies below the mean HR from the forward model and a single,
standing human during fixed-rate breathing at 0.25 Hz. The model
spectra represent the average over 10 different realizations of forward
model-generated data, in which alv = 70 ml/s.
|
|
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 Ra(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
Ra(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).

View larger version (16K):
[in this window]
[in a new window]
|
Fig. 7.
Power spectra of left ventricular flow rate (LVFR) at
frequencies below the mean HR from the forward model and a single human
tilted upright 30° with respect to the supine posture during
fixed-rate breathing at 0.25 Hz. The model spectrum represents the
average over 10 different realizations of forward model-generated data,
in which alv = 70 ml/s.
|
|
 |
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. Figure
8 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.

View larger version (21K):
[in this window]
[in a new window]
|
Fig. 8.
Mean (solid traces) and standard deviation (dotted
traces) of cardiovascular system identification estimates along
with their respective gold standards (dashed traces), which
characterize the forward model.
|
|
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 NHR 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 NHR power spectrum, which is no
longer valid in the presence of a cardiopulmonary HR baroreflex, must
be redefined (see APPENDIX C for the specific procedure).
Figure 9, A-C,
illustrates the RNMSE results of the HR baroreflex, ILV
HR, and
NHR 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.

View larger version (15K):
[in this window]
[in a new window]
|
Fig. 9.
Root-normalized-mean-squared error (RNMSE) results (mean ± SD) for HR baroreflex (A), ILV HR (B),
NHR, (C) and ILV HR (D) with
respect to the redefined gold standard (see APPENDIX C) as
a function of the ratio of the static gains of the cardiopulmonary to
arterial HR baroreflexes.
|
|
From Fig. 9A, 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 9B 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 (see
APPENDIX C for the specific procedure). Figure
9D 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 NHR estimates against more narrow
linear ABP ranges.
Figure 10 shows the RNMSE results of
the HR baroreflex, ILV
HR, and NHR 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.

View larger version (13K):
[in this window]
[in a new window]
|
Fig. 10.
RNMSE results (means ± SD) for HR baroreflex (A),
ILV HR (B), and NHR (C) as a
function of the linear ABP range (as defined in the text) of the
arterial baroreflex.
|
|
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 NHR 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.

View larger version (20K):
[in this window]
[in a new window]
|
Fig. 11.
Sensitivity results (solid trace; means ± SD) in terms of
estimated versus actual percentage of change in the static gain of
ILV HR (A), static gain of HR baroreflex (B),
absolute peak amplitude of ILV HR (C), and absolute peak
amplitude of HR baroreflex (D). The percentage of change is
with respect to the nominal values and the dashed trace is the identity
line.
|
|
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 NHR 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 NHR 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 NHR, unlike ILV
HR, is not
normalized for inputs, e.g., cerebral activity.

View larger version (19K):
[in this window]
[in a new window]
|
Fig. 12.
Sensitivity results (solid trace; mean ± SD) in terms of
estimated versus actual percentage of change in the total power
(A), low-frequency (LF) power (0.0-0.15 Hz)
(B), and high-frequency (HF) power (C)
(0.15-0.4 Hz) of the NHR power spectrum. The
percentage of change is again with respect to the nominal values, and
the dashed trace again denotes the identity line.
|
|
 |
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