Robust efficiency and actuator saturation explain
healthy heart rate control and variability
Na Li
a
, Jerry Cruz
b
, Chenghao Simon Chien
c,d
, Somayeh Sojoudi
e
, Benjamin Recht
f
, David Stone
g
, Marie Csete
h
,
Daniel Bahmiller
b
, and John C. Doyle
b,c,i,1
a
School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138;
b
Department of Computing and Mathematical Science, California
Institute of Technology, Pasadena, CA 91125;
c
Department of Electrical Engineering, California Institute of Technology, Pasadena, CA 91125;
d
Advanced
Algorithm Research Center, Philips Healthcare, Thousand Oaks, CA 91320;
e
Department of Neurology, New York University Comprehensive Epilepsy Center,
New York University School of Medicine, New York, NY 10016;
f
Department of Electrical Engineering and Computer Sciences and Department of Statistics,
University of California, Berkeley, CA 94720;
g
Departments of Anesthesiology and Neurosurgery and the Center for Wireless Health, University of Virginia
School of Medicine, Charlottesville, VA 22908;
h
Huntington Medical Research Institutes, Pasadena, CA 91101; and
i
Department of BioEngineering, California
Institute of Technology, Pasadena, CA 91125
Edited* by Michael S. Gazzaniga, University of California, Santa Barbara, CA, and approved June 27, 2014 (received for review January 30, 2014)
The correlation of healthy states with heart rate variability (HRV)
using time series analyses is well documented. Whereas these
studies note the accepted proximal role of autonomic nervous
system balance in HRV patterns, the responsible deeper physiolog-
ical, clinically relevant mechanisms have not been fully explained.
Using mathematical tools from control theory, we combine mech-
anistic models of basic physiology with experimental exercise
data from healthy human subjects to explain causal relationships
among states of stress vs. health, HR control, and HRV, and more
importantly, the physiologic requirements and constraints under-
lying these relationships. Nonlinear dynamics play an important
explanatory role
––
most fundamentally in the actuator saturations
arising from unavoidable tradeoffs in robust homeostasis and
metabolic efficiency. These results are grounded in domain-specific
mechanisms, tradeoffs, and constraints, but they also illustrate
important, universal properties of complex systems. We show that
the study of complex biological phenomena like HRV requires
a framework which facilitates inclusion of diverse domain specifics
(e.g., due to physiology, evolution, and measurement technology)
in addition to general theories of efficiency, robustness, feedback,
dynamics, and supporting mathematical tools.
system identification
|
optimal control
|
respiratory sinus arrhythmia
B
iological systems display a variety of well-known rhythms in
physiological signals (1
–
6), with particular patterns of vari-
ability associated with a healthy state (2
–
6). Decades of research
demonstrate that heart rate (HR) in healthy humans has high
variability, and loss of this high HR variability (HRV) is corre-
lated with adverse states such as stress, fatigue, physiologic se-
nescence, or disease (6
–
13). The dominant approach to analysis
of HRV has been to focus on statistics and patterns in HR time
series that have been interpreted as fractal, chaotic, scale-free,
critical, etc. (6
–
17). The appeal of time series analysis is un-
derstandable as it puts HRV in the context of a broad and
popular approach to complex systems (5, 18), all while requiring
minimal attention to domain-specific (e.g., physiological) details.
However, despite intense research activity in this area, there is
limited consensus regarding causation or mechanism and mini-
mal clinical application of the observed phenomena (10). This
paper takes a completely different approach, aiming for more
fundamental rigor (19
–
24) and methods that have the potential
for clinical relevance. Here we use and model data from exper-
imental studies of exercising healthy athletes, to add simple
physiological explanations for the largest source of HRV and its
changes during exercise. We also present methods that can be
used to systematically pursue further explanations about HRV
that can generalize to less healthy subjects.
Fig. 1 shows the type of HR data analyzed, collected from
healthy young athletes (
n
=
5). The data display responses to
changes in muscle work rate on a stationary bicycle during mostly
aerobic exercise. Fig. 1
A
shows three separate exercise sessions
with identical workload fluctuations about three different means.
With proper sleep, hydration, nutrition, and prevention from
overheating, trained athletes can maintain the highest workload
in Fig. 1 for hours and the lower and middle levels almost in-
definitely. This ability requires robust efficiency: High workloads
are sustained while robustly maintaining metabolic homeostasis,
a particularly challenging goal in the case of the relatively large,
metabolically demanding, and fragile human brain.
Whereas mean HR in Fig. 1
A
increases monotonically with
workloads, both slow and fast fluctuations (i.e., HRV) in HR are
saturating nonlinear functions of workloads, meaning that both
high- and low-frequency HRV component goes down. Results
from all subjects showed qualitatively similar nonlinearities (
SI
Appendix
). We will argue that this saturating nonlinearity is the
simplest and most fundamental example of change in HRV in
response to stressors (11, 12, 25) [exercise in the experimental
case, but in general also fatigue, dehydration, trauma, infection,
even fear and anxiety (6
–
9, 11, 12, 25)].
Physiologists have correlated HRV and autonomic tone (7, 11,
12, 14), and the (im)balance between sympathetic stimulation
and parasympathetic withdrawal (12, 26
–
28). The alternation in
autonomic control of HR (more sympathetic and less para-
sympathetic tone during exercise) serves as an obvious proximate
cause for how the HRV changes as shown in Fig. 1, but the
ultimate question remains as to why the system is implemented
Significance
Reduction in human heart rate variability (HRV) is recognized
in both clinical and athletic domains as a marker for stress or
disease, but previous mathematical and clinical analyses have
not fully explained the physiological mechanisms of the vari-
ability. Our analysis of HRV using the tools of control mathe-
matics reveals that the occurrence and magnitude of observed
HRV is an inevitable outcome of a controlled system with
known physiological constraints. In addition to a deeper un-
derstanding of physiology, control analysis may lead to the
development of timelier monitors that detect control system
dysfunction, and more informative monitors that can associate
HRV with specific underlying physiological causes.
Author contributions: N.L., J.C., B.R., and J.C.D. designed research; N.L., J.C., C.S.C., B.R.,
D.B., and J.C.D. performed research; N.L., J.C., S.S., B.R., and J.C.D. contributed new
reagents/analytic tools; N.L., J.C., C.S.C., S.S., and J.C.D. analyzed data; and N.L., D.S.,
M.C., and J.C.D. wrote the paper.
The authors declare no conflict of interest.
*This Direct Submission article had a prearranged editor.
1
To whom correspondence should be addressed. Email: doyle@caltech.edu.
This article contains supporting information online at
www.pnas.org/lookup/suppl/doi:10.
1073/pnas.1401883111/-/DCSupplemental
.
E3476
–
E3485
|
PNAS
|
Published online August 4, 2014
www.pnas.org/cgi/doi/10.1073/pnas.1401883111
this way. It could be an evolutionary accident, or could follow
from hard physiologic tradeoff requirements on cardiovascular
control, as work in other systems suggests (1). Here, the expla-
nation of HRV similarly involves hard physiological tradeoffs in
robust efficiency and employs the mathematical tools necessary
to make this explanation rigorous in the context of large mea-
surement and modeling uncertainties.
Physiological Tradeoffs
The central physiological tradeoffs in cardiovascular control (27
–
32), shown schematically in Fig. 2, involve interconnected organ
systems and four types of signals that are very different in both
functional role and time series behavior, but together define the
requirements for robust efficiency of the cardiorespiratory sys-
tem. The main control requirement is to maintain (
i
) small, ac-
ceptable
“
errors
”
in internal variables for brain homeostasis [e.g.,
cerebral blood flow (CBF), arterial O
2
saturation (
SaO
2
)] and
efficient working muscle
O
2
utilization (
Δ
O
2
) using (
ii
) actuators
(heart rate
H
, minute ventilation
_
V
E
, vasodilation and systemic
peripheral resistance
R
s
, and brain autoregulation) in response
to (
iii
) external disturbances (workload
W
), and (
iv
) internal
sensor noise and perturbations (e.g., pressure changes from
different respiratory patterns due to pulsatile ventilation
V
).
In healthy fit subjects, keeping errors in CBF,
SaO
2
,and
Δ
O
2
suitably small while responding to large, fast variations
in
W
disturbance necessitates compensating and coordinated
changes in actuation via responses in
H
,
_
V
E
,
R
s
, and cerebral
autoregulation. Thus, healthy response involves low error but
high control variability whereas loss of health is exactly the op-
posite. We will show that the observed striking changes in HRV
such as those seen in Fig. 1 result from tradeoffs between these
errors combined with various actuator saturations.
A challenge to this approach lies in managing the necessary
but potentially bewildering complexity inherent in the physio-
logical details, mathematical methods, and measurement tech-
nology. To achieve this, we make each small step in the analysis
as simple, accessible, and reproducible as possible from analysis
of experimental data to modeling to physiological and control
theoretic interpretation. In addition, we restrict the physiology
(shown schematically in Fig. 2) (27
–
32) and control theory
(32
–
35) to basic levels and all software is standard and open
source. We also make several passes through the analysis and
modeling with increasing complexity, sophistication, and depth,
to aid intuition while highlighting the need for rigorous, scal-
able methods.
In addition to mechanistic physiological models, we also use
systems identification techniques (referred to as
“
black-box
”
fits
in this paper) (25, 33, 36, 37) as intermediate steps to identify
parsimonious canonical dynamical input
–
output models relating
HR as an output variable to input disturbances such as workload
and ventilation. These techniques establish causal deterministic
links between input and output variables, highlight the aspects of
time series and dynamic relationships that are explored further,
and give some indication of the degree of complexity of their
dynamics. Then, we use physiologically motivated models (re-
ferred to as
“
first-principle
”
models in this paper) (29
–
32) to
study the mechanisms that drive the dynamics. The two approaches
are complementary: Black-box fits highlight essential relation-
ships that may be hard to intuit from data alone and can be
obscured in both complex data sets and mechanistic models,
Fig. 1.
HR responses to simple changes in muscle work rate on a stationary
bicycle: Each experimental subject performed separate stationary cycle
exercises of
∼
10 min for each workload profile, with different means but
nearly identical square wave fluctuations around the mean. A typical result is
shown from subject 1 for three workload profiles with time on the hori-
zontal axis (zoomed in to focus on a 6-min window). (
A
) HR (red) and
workload (blue); linear local piecewise static fits (black) with different
parameters for each exercise. The workload units (most strenuous exercise
on top of graph) are shifted and scaled so that the blue curves are also the
best global linear fit. (
B
) Corresponding dynamics fits, either local piecewise
linear (black) or global linear (blue). Note that, on all time scales, mean HR
increases and variability (HRV) goes down with the increasing workload.
Breathing was spontaneous (not controlled).
Fig. 2.
Schematic for cardiovascular control of aerobic metabolism and
summary of main variables: Blue arrows represent venous beds, red arrows
are arterial beds, and dashed lines represent controls. Four types of signals,
distinct in both functional role and time series behavior, together define the
required elements for robust efficiency. The main control requirement is to
maintain (
i
) small errors in internal variables for brain homeostasis (e.g.,
arterial O
2
saturation
SaO
2
, mean arterial blood pressure
P
as
, and CBF), and
muscle efficiency (oxygen extraction
Δ
O
2
across working muscle) despite (
ii
)
external disturbances (muscle work rate
W
), and (
iii
) internal sensor noise
and perturbations (e.g., pressure changes from different respiratory patterns
due to pulsatile ventilation
V
) using (
iv
) actuators (heart rate
H
, minute
ventilation
V
E
, vasodilatation and peripheral resistance
R
, and local cerebral
autoregulation).
Li et al.
PNAS
|
Published online August 4, 2014
|
E3477
PHYSIOLOGY
PNAS PLUS
whereas first-principles models give physiological interpretations
to these dynamical relationships.
Results
Static Fits.
Table 1 lists the minimum root-mean-square (rms)
error
jj
H_data-H_ fit
jj
(where
k
x
k
=
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
P
N
t
=
1
ð
x
t
Þ
2
=
N
q
for a time se-
ries
x
t
of length
N
) for several static and dynamic fits of in-
creasing complexity for the data in Fig. 1. Not surprisingly,
Table 1 shows that the rms error becomes roughly smaller with
increased fit complexity (in terms of the number of parameters).
Rows 2 and 5 of Table 1 are single global linear fits for all of
the data, whereas the remaining rows have different parameters
for each cell and are thus piecewise linear when applied to all
of the data. The
“
best
”
piecewise linear models balancing error
with complexity are further highlighted in yellow in Table 1.
We will initially focus on static linear fits (first four rows) of
the form
h
(
W
)
=
b
·
W
+
c
, where
b
and
c
are constants that
minimize the rms error
jj
H_data-h
(
W)
jj
, which can be found
easily by linear least squares. Static models have limited ex-
planatory power but are simple starting points in which con-
straints and tradeoffs can be easily identified and understood,
and we use only methods that directly generalize to dynamic
models (shown later) with modest increase in complexity. Row 1
of Table 1 is the trivial
“
zero
”
fit with
b
=
c
=
0
; row 2 is the best
global linear fit with (
b,c
)
=
(0.35,53) which is used to linearly
scale the units of
W
(blue) to best fit the HR data (red) in Fig.
1
A
; row 3 is a piecewise constant fit with
b
=
0
and
c
being the
mean of each data set; row 4 is the best piecewise linear fits
(black dashed lines in Fig. 1
A
) with quite different values (
b,c
)of
(0.44,49), (0.14,82), and (0.04,137) at 0
–
50, 100
–
150, and 250
–
300 W. The piecewise linear model in row 4 has less error than
the global linear fit in row 2. At high workload level, HR in Fig. 1
does not reach steady state on the time scale of the experiments,
the linear static fit is little better than constant fit, and so these
data are not considered further for static fits and models.
Both Table 1 and Fig. 1 imply that HR responds somewhat
nonlinearly to different levels of workload stressors. The solid
black curve in Fig. 3
A
shows idealized (i.e., piecewise linear) and
qualitative but typical values for
h
(
W
) globally that are consistent
with the static piecewise linear fits at the two lower watts levels in
Fig. 1
A
. The change in slope of
H
=
h
(
W
) with increasing
workload is the simplest manifestation of changing HRV and is
now our initial focus. A proximate cause is autonomic nervous
system balance, but we are looking for a deeper
“
why
”
in terms
of whole system constraints and tradeoffs.
Static Models.
As we mentioned earlier, in healthy fit subjects, the
central physiological tradeoffs in cardiovascular control require
keeping errors such as CBF,
SaO
2
, and
Δ
O
2
suitably small in
response to variations in
W
disturbance through changes in
actuations such as
H
. To better understand the tradeoff, we
derive a steady-state model (
P
as
,
Δ
O
2
)
=
f
(
H
,
W
) from standard
physiology that constrains the relationship between (
P
as
,
Δ
O
2
)
and (
H
,
W
) independent of how
H
is controlled (details below).
Here
P
as
is mean systemic arterial blood pressure, which is an
important variable affecting the CBF (28, 38) and
Δ
O
2
is the
drop in oxygen content across working muscle [Notice that the
model already assumes constant
SaO
2
, which is consistent with
data measurement and literature (27).] The mesh plot in Fig. 3
C
is the image on the (
P
as
,
Δ
O
2
) plane of the Fig. 3
B
(
H, W
) mesh
plot under this function
f
(
H,W
) for generic, plausible values of
physiological parameters (
SI Appendix
). Thus, any function
H
=
h
(
W
) can be mapped from the (
H
,
W
) plane using model (
P
as
,
Δ
O
2
)
=
f
(
H
,
W
) to the (
P
as
,
Δ
O
2
) plane to determine its con-
sequences for the most important tradeoffs, which involve
P
as
and
Δ
O
2
. These results are shown with the black lines in Fig. 3
B
,
which give
H
=
h
(
W
) curves consistent with Fig. 3
A
and then
are mapped onto Fig. 3
C
.
Hidden complexity is unavoidable in the model (
P
as
,
Δ
O
2
)
=
f
(
H
,
W
), but we temporarily defer these details to focus on the
general shape of the color-coded curves in Fig. 3
B
and
C
, which
have an intuitively clear explanation highlighted by the dashed
red and purple lines. At constant workload, increased HR would
greatly increase
P
as
while slightly decreasing
Δ
O
2
due to greater
flow rate through the muscle. For constant HR, increased
workload would greatly increase
Δ
O
2
while slightly reducing
P
as
due to greater oxygenation and peripheral vasodilatation. The
cardiovascular control system adjusts HR as a function
H
=
h
(
W
)
of workload to tradeoff increasing
P
as
with increasing
Δ
O
2
, both
of which are undesirable. The modest curvature of the colored
meshes in Fig. 3
C
demonstrates a small nonlinearity in the
function (
P
as
,
Δ
O
2
)
=
f
(
H
,
W
). One source of this nonlinearity is
the nonlinear relationship between cardiac output and HR due
to less diastolic filling time as HR increases. However, the solid
black lines in Fig. 3 manifest a much larger nonlinearity in the
control function
H
=
h
(
W
). We will argue that the essential
sources of this nonlinearity are the tradeoff in robust homeo-
stasis and metabolic efficiency and how it changes at different
HR levels.
The hypothetical linear response at low workload in Fig. 3 can
be explained in terms of purely metabolic tradeoffs. Healthy
athletes can maintain the low workload almost indefinitely even
in adverse (e.g., heat) conditions, a feature of human physiology
thought to be an important adaptation for a successful hunter
(39). Prolonged exercise necessarily requires steeply increased
HR to provide sufficient tissue O
2
(low
Δ
O
2
), to maintain aer-
obic lipid metabolism in muscles and preserve precious carbo-
hydrates for the brain.
The nonlinear response in Fig. 3 (solid lines) reflects addi-
tional tradeoffs that arise at higher workload and HR, when the
resulting high
P
as
becomes dangerous mainly due to actuator
saturation of cerebral autoregulatory control. In healthy humans,
CBF is autoregulated to be quite constant (28, 38) over a rela-
tively wide range of
P
as
(50
<
P
as
<
150 mm Hg), so that no new
tradeoffs at moderate exercise levels are required, because
P
as
is
within this range. A new tradeoff does arise at
P
as
above 150 mm
Hg when cerebral autoregulation saturates, and CBF begins to
rise with the severe possible consequences of edema and/or
hemorrhage. Thus, for the dashed black linear response in Fig. 3
Table 1. rms error for models of different complexity for data in Fig. 1
W
0
0
3
–
0
5
2
W
0
5
1
–
0
0
1
W
0
5
–
0
e
r
u
t
c
u
r
t
s
l
e
d
o
M
s
r
e
t
e
m
a
r
a
p
.
o
N
w
o
R
1
0
Zero:
h
(
W
9
5
.
8
4
1
7
9
.
9
9
0
4
.
0
6
0
=
)
2
2
Global static:
h
(
W
)
= b·W + C
10.1
6.9
10.4
3
3 = 1 × 3
Piecewise constant
h
i
(
W
)
= C
i
14.8
4.7
6.7
4
6 = 2 × 3
Piecewise static
h
i
(
W
)
= b
i
·W + C
i
9.6
3.2
6.6
5
3
Global fi
rst order
Δ
h
(
t
)
= ah
(
t
)
+ bW + C
11.6
3.2
6.5
6
9 = 3 × 3
Piecewise fi
rst order
Δ
h
i
(
t
)
= ah
i
(
t
)
+ b
i
W + C
i
8.9
2.1
0.9
Items in yellow highlight the best models balancing fitting error with model complexity. For the piecewise model,
i
=
1,2,3 stands for
each exercise level, respectively; i.e., 0
–
50, 100
–
150, 250
–
300 W.
E3478
|
www.pnas.org/cgi/doi/10.1073/pnas.1401883111
Li et al.
B
and
C
, the resulting
P
as
would be elevated to potentially
pathologic levels, and some nonlinearity as in the solid black
line is necessary. Moreover, in many subjects there may be
diminishing metabolic benefit of high tissue O
2
(low
Δ
O
2
) at high
workloads because muscle mitochondria saturate. Although many
details of cerebral autoregulat
ion (as well as the mitochondrial
saturation) are poorly understood, the
P
as
at which autoregulation
saturates is well-known in healthy adults, and helps to explain an
important change in HRV with stressors. Ultimately, cardiac output
itself saturates at sufficiently high HR due to compromised diastolic
filling time with subsequent dramatic falls in stroke volume.
Mathematically, all these factors can be quantitatively reflec-
ted in a static optimization model using linear least squares, with
H
=
h
(
W
) chosen to minimize a weighted penalty on increasing
P
as
,
Δ
O
2
, and
H
:
min
q
2
P
P
as
−
P
*
as
2
+
q
2
o
2
Δ
O
2
−
Δ
O
p
2
2
+
q
2
H
ð
H
−
H
p
Þ
2
subject to linearization of the constraint (
P
as
,
Δ
O
2
)
=
f
(
H
,
W
)
at 0 and 100 W. Here
P
p
as
;
Δ
O
p
2
;
H
p
are the steady values for
P
as
;
Δ
O
2
;
H
at 0 and 100 W, respectively. Different values for
ð
q
P
;
q
O
2
;
q
H
Þ
reflect different tradeoffs between
P
as
,
Δ
O
2
,and
H
at different workloads. In particular,
q
P
is higher at high work-
loads and high HR, reflecting the greater impact of
P
as
on CBF
due to saturation of autoregulation, and
q
H
is higher to reflect
the saturation of HR itself, which becomes more acute at higher
workload levels. Straightforward, standard computations easily
reproduce the piecewise linear features in Fig. 3 with higher
penalty on
P
as
and
H
at higher workload levels.
An important feature of this approach is that it allows sys-
tematic exploration of models that are both simple and explan-
atory. We have systematically moved from the data in Fig. 1 to
the fit in Fig. 3
A
, and then from very simple well-understood
physiological mechanisms to how healthy HR should behave and
be controlled, reflected in Fig. 3
B
and
C
. The nonlinear be-
havior of HR is explained by combining explicit constraints in the
form (
P
as
,
Δ
O
2
)
=
f
(
H
,
W
) due to well-understood physiology
with constraints on homeostatic tradeoffs between rising
P
as
and
Δ
O
2
that change as
W
increases. The physiologic tradeoffs
depicted in these models explain why a healthy neuroendocrine
system would necessarily produce changes in HRV with stress,
no matter how the remaining details are implemented. Taken
together this could be called a
“
gray-box
”
model because it
combines hard physiological constraints both in (
P
as
,
Δ
O
2
)
=
f
(
H
,
W
) and homeostatic tradeoffs to derive a resulting
H
=
h
(
W
). If new tradeoffs not considered here are found to be
significant, they can be added directly to the model as additional
constraints, and solutions recomputed. The ability to include
such physiological constraints and tradeoffs is far more essential
to our approach than what is specifically modeled (e.g., that
primarily metabolic tradeoffs at low HR shift priority to limiting
P
as
as cerebral autoregulation saturates at higher HR). This ex-
tensibility of the methodology will be emphasized throughout.
The most obvious limit in using static models is that they omit
important transient dynamics in HR, missing what is arguably the
most striking manifestations of changing HRV seen in Fig. 1.
Fortunately, our method of combining data fitting, first-princi-
ples modeling, and constrained optimization readily extends
beyond static models. The tradeoffs in robust efficiency in
P
as
and
Δ
O
2
that explain changes in HRV at different workloads also
extend directly to the dynamic case as demonstrated later.
Dynamic Fits.
In this section we extract more dynamic information
from the exercise data. The fluctuating perturbations in work-
load (Fig. 1) imposed on a constant background (stress) are
targeted to expose essential dynamics, first captured with
“
black-
box
”
input
–
output dynamic versions of above static fits. Fig. 1
B
shows the simulated output
H
(
t
)
=
HR (in black) of simple local
(piecewise) linear dynamics (with discrete time
t
in seconds)
Δ
H
ð
t
Þ
=
H
ð
t
+
1
Þ
−
H
ð
t
Þ
=
Hh
ð
t
Þ
+
bW
ð
t
Þ
+
c
;
[1]
where the input is
W
(
t
)
=
workload (blue). Constants (
a
,
b
,
c
)
are fit to minimize the rms error between
H
(
t
) and HR data
as before (Table 1). The optimal parameter values (
a
,
b
,
c
)
∼
(
−
0.22, 0.11, 10) at 0 W differ greatly from those at 100 W
(
−
0.06, 0.012, 4.6) and at 250 W (
−
0.003, 0.003,
−
0.27), so a
single model equally fitting all workload levels is necessarily
nonlinear. This conclusion is confirmed by simulating HR
(blue in Fig. 1
B
) with one best global linear fit (
a
,
b
,
c
)
∼
(0.06,0.02,2.93) to all three exercises, which has large errors
at high and low workload levels.
The changes of the large, slow fluctuations in both HR (red)
and its simulation (black) in Fig. 1
B
are consistent with well-
understood cardiovascular physiology, and illustrate how the
physiologic system has evolved to maintain homeostasis despite
stresses from workloads. Our next step in modeling is to mech-
anistically explain as much of the HRV changes in Fig. 1 as
possible using only standard models of aerobic cardiovascular
physiology and control (27
–
31). This step focuses on the changes
in HRV in the fits in Fig. 1
B
(in black) and Eq.
1
, and we defer
Fig. 3.
Static analysis of cardiovascular control of aerobic metabolism as workload increases: Static data from Fig. 1
A
are summarized in
A
and the physi-
ological model explaining the data is in
B
and
C
. The solid black curves in
A
and
B
are idealized (i.e., piecewise linear) and qualitatively typical values for
H
=
h
(
W
) that are globally consistent with static piecewise linear fits (black in Fig. 1
A
) at the two lower workload levels. The dashed line in
A
shows
h
(
W
) from
the global static linear fit (blue in Fig. 1
A
) and in
B
shows a hypothetical but physiologically implausible linear continuation of increasing HR at the low
workload level (solid line). The mesh plot in
C
depicts
P
as
–
Δ
O
2
(mean arterial blood pressure
–
tissue oxygen difference) on the plane of the
H
–
W
mesh plot in
B
using the physiological model (
P
as
,
Δ
O
2
)
=
f
(
H
,
W
) for generic, plausible values of physiological constants. Thus, any function
H
=
h
(
w
) can be mapped from
the
H
,
W
plane
(B)
using model
f
to the (
P
,
Δ
O
2
) plane (
C
) to determine the consequences of
P
as
and
Δ
O
2
. The reduction in slope of
H
=
h
(
W
) with increasing
workload is the simplest manifestation of changing HRV addressed in this study.
Li et al.
PNAS
|
Published online August 4, 2014
|
E3479
PHYSIOLOGY
PNAS PLUS
modeling of the high-frequency variability in Fig. 1 until later
(i.e., the differences between the red data and black simulations
in Fig. 1
B
).
The black-box fits allow us to plausibly conjecture that work-
load disturbances cause most of the variability in Fig. 1
B
(black
curves). Here the rigor of the black-box fits is important, as
highlighted by three features: (
i
) no comparably good fits exist
for the data in Fig. 1 without the input of workload, (
ii
) within
the limits of the sensors used and subject fitness we can other-
wise experimentally manipulate the input independently and
over a wide range to make it truly a
“
causal
”
input, and (
iii
) the
fits accurately predict the HR output response to new experi-
ments (i.e., cross-validation; see
SI Appendix
).
First-Principles Models.
Our first-principles model is based on the
circulatory circuit diagram in Fig. 2, using standard mathematical
descriptions of circulation, and with a focus on modeling purely
aerobic exercise. That is, we only model blood flow, blood
pressure, and O
2
in several compartments, and yet the model
captures the overall physiologic HR response during moderate
exercise in young, fit adults. In standard models of aerobic car-
diovascular control (27
–
31) the neuroendocrine system controls
peripheral vasodilation, minute ventilation, and cardiac output
to maintain blood pressure and oxygen saturation within ac-
ceptable physiological limits.
Several features of these control systems allow substantial
simplification of the model. Minute ventilation
_
V
E
alone can
tightly control arterial oxygenation [
O
2
]
a
, so we assume [
O
2
]
a
is
maintained nearly constant (27). Moreover, peripheral resistance
R
s
is decreased during exercise and the decrease is determined by
local metabolic control. The purpose of decreasing
R
s
in the
arterioles is to increase blood flow and regional delivery of O
2
,
glucose, and other substrates as needed. Because the venous
oxygenation [
O
2
]
v
serves as a good signal for oxygen consump-
tion, we also assume that control of peripheral vascular re-
sistance
R
s
is a function only of venous oxygenation [
O
2
]
v
(31).
Combined with those models for blood circulation and oxygen
consumption, we have the following physiological model:
V
as
=
c
as
·
P
as
V
vs
=
c
vs
·
P
vs
V
ap
=
c
ap
·
P
ap
V
vp
=
c
vp
·
P
vp
V
tot
=
V
as
+
V
vs
+
V
ap
+
V
vp
½
O
2
a
=
0
:
2
Q
l
=
c
l
·
H
·
P
vp
Q
r
=
c
r
·
H
·
P
vs
F
s
=
ð
P
as
−
P
vs
Þ
=
R
s
F
p
=
P
ap
−
P
vp
R
p
M
=
ρ
·
W
+
M
0
R
s
=
A
·
½
O
2
v
+
R
s
0
:
[2]
Here
V
and
P
are (subscripts
a
=
arterial,
v
=
venous,
s
=
sys-
temic,
P
=
pulmonary) blood volume and blood pressure, re-
spectively. All of the
c
variables are constants. The main
elements of the model are (more details in
SI Appendix
): (
i
)
arterial and venous compartments of systemic and pulmonary
circulations are treated as compliant vessels, modeled in the
form
V
=
c
·
P
, with the total blood volume a constant
V
tot
;(
ii
)
cardiac output of the left (
Q
l
) and right (
Q
r
) ventricles; (
iii
) blood
flow for systemic (
F
s
) and pulmonary (
F
p
) circulation; (
iv
) the
metabolic consumption
M
;(
v
)[
O
2
]
a
and
R
s
are modeled accord-
ing to the previous description of the control mechanism. Note
that we need not model these control systems in detail, but
simply extract their most well-known features and use them to
constrain the model.
In steady state the follow additional constraints hold:
Q
r
=
Q
l
=
F
s
=
F
p
M
=
F
s
½
O
2
a
−
½
O
2
v
[3]
The first equation is total blood circulation balance and the
second one is based on the oxygen circulation balance, where
F
s
([
O
2
]
a
−
[
O
2
]
v
) is the net change in the arterial and venous
blood O
2
content. The oxygen drop
Δ
O
2
across the muscle bed is
defined as
Δ
O
2
=
[
O
2
]
a
−
[
O
2
]
v
.
Combining
2
and
3
plus simple
algebra (
SI Appendix
) gives the steady-state model (
P
as
,
Δ
O
2
)
=
f
(
H
,
W
) shown in Fig. 3 that constrains the relationship between
(
P
as
,
Δ
O
2
) and (
H
,
W
).
In general the circulatory system is far from steady state in our
experiments. Modeling the blood volume change for each cir-
culatory compartment and the oxygen change in the tissue keeps
the constraints from Eq.
2
but replaces
3
with the following dy-
namic model:
c
as
_
P
as
=
Q
l
−
F
s
c
vs
_
P
vs
=
F
s
−
Q
r
c
ap
_
P
ap
=
Q
r
−
F
p
v
T
;
O
2
_
O
2
v
=
−
M
+
F
s
·
½
O
2
a
−
½
O
2
v
c
vp
P
vp
=
V
total
−
c
as
P
as
+
c
vs
P
vs
+
c
ap
P
ap
:
[4]
Here
v
T
;
O
2
denotes the effective tissue O
2
volume and we assume
that tissues and venous blood gases are in equilibrium, namely
that tissue oxygenation
½
O
2
T
is the same as venous oxygenation
½
O
2
v
(
SI Appendix
). The previous static analysis (and the purely
static tradeoffs it highlights) directly extends to the dynamic case
with modestly increased complexity. The simplest extension is to
use an optimal linear quadratic state feedback controller (34) for
linearizations of
4
at 0 and 100 W, with controller
H
=
u
ð
·
Þ
chosen to minimize a weighted penalty on integrated elevation
of
P
as
,
Δ
O
2
, and
H
:
min
Z
q
2
P
P
as
−
P
p
as
2
+
q
2
o
2
Δ
O
2
−
Δ
O
p
2
2
+
q
2
H
ð
H
−
H
p
Þ
2
dt
[5]
subject to linearizations of the state dynamic constraint
4
. Fig. 4
compares HR and workload data versus simulations of applying
the linear controllers to the model in
4
for two experiments
Fig. 4.
Optimal control model response using first-principle model to two
different workload (blue) demands, approximately square waves of 0
–
50 W
(
Lower
) and 100
–
150 W (
Upper
): For each data set (using subject 2
’
s data),
a physiological model with optimal controller is simulated with workload as
input (blue) and HR (black) as output, and compared with collected HR data
(red). Simulations of blood pressure (
P
as
, purple) and tissue oxygen satura-
tion ([
O
2
]
T
, green) are consistent with the literature but data were not col-
lected from subjects. Breathing is spontaneous (not controlled).
E3480
|
www.pnas.org/cgi/doi/10.1073/pnas.1401883111
Li et al.
(similar to Fig. 1 but with a different subject) with higher penalty
on
P
as
and
H
at higher workloads as in the static case. (See
SI
Appendix
for more details.) Also shown are simulations of
P
as
and [
O
2
]
T
, which are consistent with the literature but were not
measured. The same methods and results apply to other subjects
’
data and new experiments (e.g., cross-validation; see
SI Appen-
dix
.). Also note that for different subjects, we use the same nom-
inal values for most parameters except the constants
c
r
,
c
l
used
in the cardio-output formulation and the weighting parameters
q
P
;
q
O
2
;
q
H
used in the tradeoff function
5
. The different
c
r
,
c
l
reflect the different sizes of subjects and different
q
P
;
q
O
2
;
q
H
reflect the qualitative differences in the control objectives.
The change in the tradeoff as workload increases is consistent
with what we observed using the static model. At low workload
and low HR, the main tradeoff is metabolic because both
P
as
and HR are at safe and sustainable levels. High HR and thus
high [
O
2
]
T
, (low
Δ
O
2
), maintains aerobic muscle metabolism,
extending the potential duration of exercise while preserving
carbohydrate resources for the brain. At higher workloads, this
strategy would produce unsustainably high and potentially dam-
aging
P
as
and possibly HR, so the optimal controller penalizes
these factors more, at the price of reduced [
O
2
]
T
. HRV (slow time
scale) in Fig. 1 (and
P
as
in Fig. 4) decreases with increasing
workload because of the changing tradeoffs between metabolic
overhead and
P
as
,
Δ
O
2
,and
H
as their means increase. These
explanations in HRV derived from the dynamic aerobic model
are richer and more complete but due to the same tradeoffs as
in the simpler static model.
Importantly, although the mathematics and physiology re-
quired are relatively elementary and the resulting explanation is
intuitively clear and mechanistic, they nonetheless highlight the
rigor and scalability of this approach. The simplicity of the black-
box fits in
1
and Fig. 1 helps establish causal relationships be-
tween variables and suggests physiological mechanisms to model
in more detail, and highlights features in the signals that are not
modeled (i.e., we have not explained the high frequency of the
signals at low watts in Fig. 1, considered in the next sections).
The hard homeostatic tradeoffs and the actuator effects of HR,
ventilation, and vasodilation were included in the physiology
model in
2
–
5
but the neuroendocrine implementation details
were not. Also, the impact of cerebral autoregulatory saturation
was included, but the details of implementation were not.
Nonetheless, this approach allows for clinically actionable explan-
ations that do not depend on poorly understood mechanisms
peripheral to the component being modeled, and provides a
framework for systematically refining such models using a similar
(but presumably vastly more complex) combination of black- and
gray-box models and physiology. Again, if new tradeoffs not
considered here are found to be significant, they can be added
directly as additional constraints are recognized and solutions
recomputed. Further, tradeoffs may well change as organ sys-
tems fail, when these models are extended to disease states.
High-Frequency HRV.
The high-frequency fluctuations in HR that
are particularly large at low mean workloads and low HR cannot
be explained by the simple fits or models above, and thus addi-
tional signals and mechanisms must be included that can be
causally related to this setting. Figs. 5 and 6 shed light on
breathing as a cause of much of the high-frequency HRV. Fig.
5
A
shows HR (red) during natural breathing at rest, with the
Fig. 5.
HR response (red) to ventilation
V
(blue) at rest (0 W): The ventilatory data are raw speed of inhalation and exhalation measured at the mouthpiece.
In each case the units for
V
(blue) are chosen to show the optimal static fit
h
(
V
)
=
b
·
V
+
c
to the collected HR data.
A
and
B
show natural breathing, with
B
zoomed in to focus on a smaller window to help visualize the data and fit.
C
and
D
are similar focused smaller windows from a longer controlled breathing
experiment at resting (0 W) where the subject followed a frequency sweep from fast to slow breathing (see Fig. 6 for the full frequency range).
C
focuses on
breathing frequencies close to natural breathing, whereas
D
focuses on frequencies slower than natural. Both
C
and
D
show simulated dynamic fits (black),
and optimal static fits
h
(
V
)
=
b
·
V
+
C
(blue). The dynamic fits improve on the static fits more for the controlled sweep than for natural breathing (Table 2).
Li et al.
PNAS
|
Published online August 4, 2014
|
E3481
PHYSIOLOGY
PNAS PLUS
optimal static fit
H
=
h
(
V
)
=
b
·
V
+
c
, where
V
(blue) is measured
ventilation flow rates (inhalation and exhalation) at the mouth-
piece. The static linear fit can be used to scale the units of
V
(blue) in Fig. 5
A
to visualize the best fit to HR (red) and its
error, shown in Table 2. That HR and ventilation match so well
in frequency is consistent with the observation that under certain
conditions, inspiration is accompanied by an acceleration of
HR, and expiration by a deceleration, a phenomenon called
respiratory sinus arrhythmia (RSA) (16, 40
–
46). However, be-
cause ventilation and HR are both generated by neuroendocrine
control, this fit (i.e., correlation) by itself does not suggest a
specific mechanistic explanation of the resulting ventilation
–
HR
correlation or HRV.
To sharpen this picture, Fig. 6 shows data from subjects
instructed to control respiratory rate (RR, shown in blue) by
following a computer-generated RR frequency sweep (tidal
volumes not controlled), repeated with a background of 0- and
50-W exercises, respectively. (Fig. 5
B
shows HR and zoomed in
for the 0-W exercise data.) For each exercise taken separately,
HR is fit with static (blue in Fig. 5
B
) and one-state models, as
well as a two-state, five-parameter linear model (shown in black
in Figs. 5
B
and 6, Table 2):
Δ
H
ð
t
Þ
=
a
1
H
ð
t
Þ
+
b
1
V
ð
t
Þ
+
x
ð
t
Þ
Δ
x
ð
t
Þ
=
a
2
x
ð
t
Þ
+
b
2
V
ð
t
Þ
+
c
;
where
V
is ventilatory flow rates,
x
is an internal black-box state,
and the parameters depend on workload. Whereas breathing
cannot be varied as systematically and widely as workload, these
black-box fits provide strong evidence that ventilation is the main
factor causing high-frequency HRV. The underlying physiologi-
cal mechanisms remain unclear, but we now know where to look
next. In Fig. 6, minute ventilation naturally increases at 50 W, yet
HRV goes down, a nonlinear pattern consistent with the trends
in Figs. 1 and 3, but more dramatic. As Table 2 shows, dynamic
fits have little benefit for natural breathing at rest, but modestly
reduce the error for the controlled respiratory sweeps at low- and
high-frequency breathing. In all cases, the frequency of HR oscil-
lations is fit better than the magnitude, suggesting both a dynamic
and nonlinear dependence of HR on ventilatory flow rates.
Although RSA magnitude has been used as a measure of vagal
function, after many years of research the mechanism of RSA,
e.g., whether RSA is due to a central or a baroreflex mechanism,
is still debated (16, 40
–
46). Moreover, the data and dynamic fits
show a small resonant peak in the frequency response at around
0.1 Hz at 0 W, and the significance of the peak is unclear. Of
note, this characteristic peak occurred in the fits for every subject
(although the exact RR at which it occurs varied), and is con-
sistent with observations in the literature (16, 40). (In
SI Appendix
,
we also use both workload and ventilation data as inputs to fit
HR data during the easy workout in Fig. 1.)
Resolution of these mysteries requires additional measure-
ments such as arterial blood pressure (more invasive human
studies), more sophisticated physiological modeling including the
mechanical effects of breathing on arterial blood pressure and
pulmonary stretch reflexes, plus changing tradeoffs in control of
arterial blood pressure at different workloads and HR levels. In
particular, model
2
–
5
above assumed continuous ventilation and
HRs (i.e., no intrabreath or -beat dynamics), so more detailed
modeling of physiological respiratory patterns and their me-
chanical and metabolic effects is needed.
Discussion
Robust Efficiency and Actuator Saturation.
We showed how HR
fluctuations in healthy athletes can be largely explained as non-
linear dynamic, but not chaotic, responses to either external (e.g.,
workload) or internal (e.g., ventilation implemented by pulsatile
breathing) disturbance. We provided mechanistic explanations
and plausible conjectures for essentially all of the HRV in Fig. 1,
and showed that changes in HRV per se, no matter how mea-
Fig. 6.
HR response to sweep ventilation on different workload levels: Two experiments with
(A)
HR (red) and dynamic fit (black) to input of controlled
ventilation frequency sweeps with measured ventilatory flow rate (blue) on a fixed background workload of (
B
)0Wor(
C
) 50 W. Ventilatory flow (with
spontaneous ventilation magnitude) was necessarily larger at 50 W and the subject was unable to breathe slowly enough to complete the entire frequenc
ysweep.
E3482
|
www.pnas.org/cgi/doi/10.1073/pnas.1401883111
Li et al.
sured, are much less important mechanistically than the tradeoffs
that produce them. The tradeoffs we highlight between robust-
ness, homeostasis, and metabolic efficiency are universal and
essential (1, 23) features of complex systems but can remain
hidden and cryptic (47) without an appropriate mathematical
framework (4, 48).
“
Universal
”
features illustrated by this phys-
iological (HR) control system include how efficiently maintain-
ing robust homeostasis (e.g., small errors in CBF,
S
a
O
2
, and
Δ
O
2
)
in the presence of large disturbances requires correspondingly
large actuator (e.g., HR, ventilat
ion, and cerebral autoregulation)
responses to compensate, and how nonlinearities in actuator satu-
rations lead to reductions in actuator activity (e.g., HRV) under
increased load or stress.
HR control and HRV highlight layered control, actuator
constraints, and hard tradeoffs of the type that pervade physi-
ology and are generally fundamental in complex control systems.
In summary, actuators are the mechanisms by which controllers
act on the system to provide efficient performance and robust
homeostasis. In the cardiovascular models so far, the most im-
portant saturation is in cerebral autoregulation, which forces
a nonlinear change in
H
=
h
(
W
) as workload increases to avoid
high
P
as
that leads to intracerebral pathology (edema, hemor-
rhagic stroke).
Further understanding of control complexity and the role of
actuator (e.g., HR) variability and saturation in physiologic
control comes from examination of other technological examples
of complex systems. Familiar examples include automobiles,
particularly new autonomous robotic versions. A car is moved
and controlled via the actuators that produce and deliver power,
braking, and steering that result in accelerations in forward,
backward, and lateral directions. Most complexity in an auton-
omous robot car is in the control system for robust efficiency to
uncertainty in real traffic environments and to intrinsic vari-
ability in components. If scientists were forced to
“
reverse en-
gineer
”
such a car without access to the forward engineering
process, they would likely study
“
knockouts
”
of components to
infer their function, and also push the vehicles to extremes to
find the limits of their robust performance. It would be surprising
if
“
reverse engineering
”
cardiovascular control would be easier
than a robot car, or could be accomplished with less sophisti-
cated tools and without domain-specific details.
Loss of car actuator variability due to
“
stress
”
parallels loss of
HRV, in that it is loss of actuator responsiveness that causes
deterioration of function, and loss of variability is only a symp-
tom of actuator dysfunction. Currently, human drivers cause
most crashes when they reduce their actuator responsiveness
because of multitasking, alcohol consumption, fatigue, or poor
visibility, or when surface conditions make the actuators less
effective. Automatic collision avoidance and antilock and anti-
slip traction control systems mitigate these effects, and augment
human control in emergencies. However, even in fully auto-
mated robotic cars with robust control systems, at extremes of
speed, acceleration, braking, and turning (such as a race scenario
or in icy conditions), actuators would frequently saturate and
lose variability, resulting in less maneuverability, and an increase
in errors and risk of crashes. Thus, actuator saturation causing
changes in variability is a
“
signature
”
of a wide variety of dan-
gerous scenarios, and essential to understanding vehicle limits
and malfunction. However, variability per se is unimportant, and
analyzing the statistics of individual signals (e.g., fuel or air rates,
braking, acceleration, turning rate, etc.) in isolation is relatively
less diagnostic compared with understanding integrated, mech-
anistic dynamic models of signal interactions.
Similar tradeoffs to those resulting in HRV are found throughout
technology and biology. For example, glycolytic oscillations were
one of the most persistent mysteries involving dynamics in cell
biology (1). The proximal role of how autocatalytic and regula-
tory feedbacks make oscillations possible was well understood,
but the unresolved deeper
“
why
”
question was the purpose of the
oscillations, or alternatively, whether they were just frozen acci-
dents of evolution. Oscillations are neither functional nor acci-
dental but are a side effect of provably hard tradeoffs involving
efficiency and robust control (1). The glycolysis circuit must
maintain adequate ATP concentrations that are robust to fluc-
tuations in demand and to enzyme and other metabolite levels. It
must also be metabolically efficient, in the sense of not requiring
excessive enzyme concentrations. Any circuit that aims to bal-
ance these competing requirements has the potential to oscillate,
particularly when enzymes saturate.
Mathematical Framework.
It has been difficult to characterize
multilayered aspects of biological control, but our approach is
aimed at providing tools for biologists and clinicians, combining
established principles of system identification fits and control
theory with basic physiological models. The fits in Figs. 1, 5, and
6 highlight causal input
–
output relationships between variables
and help suggest the relevant physiology. By comparison, even
the most sophisticated statistical analysis of individual HR
signals taken out of physiologic context is mechanistically un-
informative. In contrast, the static and dynamic models mecha-
nistically explain Figs. 3 and 4 and most of the variability in Fig. 1
(i.e., the black curves at the lower two watts). Our explanation in
terms of aerobic metabolism is simple and intuitive as well as
mechanistic, and requires only basic mathematics and physiol-
ogy. The main requirement of the models is some mechanistic
relationship between control actuation and its limits in main-
taining robust homeostasis. Thus, we did not need detailed
understanding of neuroendocrine control implementation or
peripheral autoregulation, but only that they adequately manage
the tradeoffs and saturation effects in muscle, brain, and heart as
described above. Moreover, specific details of the computational
approach, e.g., piecewise linear least squares used in this paper,
are not essential to understanding the underlying system control.
What is important is that the right constraints are properly
reflected in the computation, so that the resulting controller
function is constrained by the right physiological mechanisms
plus appropriately changing penalties
–
constraints on vital phys-
iological variables due to metabolic tradeoffs and limit.
Our approach also importantly highlights where mechanisms
are missing. The model in
2
–
5
and Fig. 4 assume continuous
ventilation and HRs (i.e., no intrabreath or -beat dynamics) and
do not capture any of the higher frequency HRV in Fig. 1. The
fits in Figs. 5 and 6 suggest that the dynamics of physiological
respiratory patterns and their mechanical and metabolic effects
Table 2. rms error for models of different complexity for data in Figs. 5 and 6
Row no.
No. prms
Model structure
Resting natural
Resting sweep
50-W sweep
8
9
7
6
6
5
o
r
e
Z
0
1
9
.
7
1
.
3
1
6
.
1
1
t
n
a
t
s
n
o
C
1
2
9
.
5
7
.
9
4
.
6
c
i
t
a
t
S
2
3
4
3
First-order dynamic
6.0
9.7
5.5
5
6
Second-order dynamic
6.0
7.5
4.0
Items in yellow highlight the best models balancing fitting error with model complexity.
Li et al.
PNAS
|
Published online August 4, 2014
|
E3483
PHYSIOLOGY
PNAS PLUS
are necessary and perhaps sufficient to explain the high-frequency
HRV seen in Fig. 1. There may also be connections between
robust efficiency and oscillations as in ref. 1 to explain the origin
of the peaks in frequency response of the breath-to-HR fits.
An essential feature of this project is that our tools be robust
and scalable to more complex signals and models, and that if new
mechanisms and/or tradeoffs are discovered that are important,
they can be added as additional constraints. Fortunately, we can
leverage enormous recent advances in engineering theory and
practice, although these remain largely unknown in mainstream
science outside of the most advanced parts of systems biology.
The general models and methods, particularly moving from
1
to
2
–
5
(and Fig. 1 to Figs. 3 and 4), used for this relatively straight-
forward study should serve well as a foundational framework for
the evaluation of even more complex physiologic (disease) sit-
uations in which the diagnostic possibilities are broader.
Clinical Correlates: Linking the Behavior of Control Systems and
Pathophysiology.
Clinicians know that changes in actuator sig-
nals (e.g., HR increases) can signify a great variety of potentially
important derangements such as hypovolemia, congestive heart
failure, inflammation, sepsis, etc. (6
–
13). Even without specific
diagnostic content, alerts to clinicians that HRV is changing can
be useful. Such an alert has been incorporated into monitoring of
premature neonates. In this case, monitoring of HR character-
istics is used to report a statistic that reflects the performance of
the actuator mechanisms. The utilization of this statistic by
clinicians was sufficient to change practice (i.e., reevaluation of
the patient with consideration of need for additional therapies)
and achieve a significant improvement in outcomes (10, 13). This
particular example of integration of mathematical analysis into
monitors represents a special situation because the alert also
incorporates the occurrence of
cardiac decelerations, a phe-
nomenon not observed in adults. Furthermore, because sepsis is
a common problem in this setting, the clinicians chose to ad-
minister antimicrobial therapy on the basis of the alert which
produced the outcome benefit. In fact, the efficacy of the mon-
itor in early sepsis detection was subsequently demonstrated in
a post hoc analysis (13). However, in most clinical scenarios,
actuator changes alone are usually so generic that they lack
specific diagnostic value, and extensive analyses of individual
time series have not yielded mechanistic explanations that can
narrow the diagnosis (6
–
12).
In contrast, the general type of models and methods used here
and the application of control theory to physiology present an
enormous opportunity to reexamine this area with powerful
mathematical tools and a systems engineering approach (48).
This is important because system dysfunction is manifested
earlier in the behavior of the control system than in any metric
associated with the system
’
s output. The long-term goal of this
research is earlier diagnosis afforded by monitoring control
elements in addition to individual signal outputs.
Materials and Methods
After Caltech Institutional Review Board approval, five fit athletic subjects
(ages 25
–
35 y, three men and two women) performed a series of experi-
mental exercise regimens, each on two different days. The intensity and
durations were less than routine training for these athletes, and they had
used the laboratory equipment before so were familiar and comfortable
with the environment. In all experiments, the subject pedaled a Life Fitness
stationary recumbent bicycle at near-constant speed with the pedaling re-
sistance controlled by a preprogrammed protocol. In the respiration rate
(RR) sweep experiments, a sinusoid signal was preprogrammed in the com-
puter with frequency from 2 Hz to 0.06 Hz and each subject watched the
signal and controlled RR to follow the frequency of the signal until they
were unable to continue.
In all exercise tests, 1-Hz work rate data were recorded from a Life Fitness
stationary recumbent exercise bicycle interfaced with a computer running
MATLAB via the Communications Specification for Fitness Equipment pro-
tocol. Other exercise testing data were collected using commercially avail-
able noninvasive monitors. (
i
) RR interval HR data were recorded with
a Polar heart rate monitor and converted to 1-Hz HR data using the integral
pulse frequency modulation model (49). (
ii
) In the tests shown in Figs. 5 and
6, 100-Hz ventilation (inhalation and exhalation) flow rate data were
recorded with a Philips NiCO
2
monitor, and down-sampled to 1-Hz ventila-
tion flow rate data. (
iii
) In the test shown in Fig. 7, 1-Hz gas data including
minute ventilation
_
V
E
, oxygen consumption
_
VO
2
and carbon dioxide
generation
_
VCO
2
were collected with a Vacumetrics monitor and TurboFit
software. No other preprocessing of data was performed.
A detailed discussion of mathematical methods is given in
SI Appendix
.
ACKNOWLEDGMENTS.
We thank Pamela B. Pesenti for her gift in establish-
ing the John G. Braun Professorship, which supported this research, and
Philips for providing equipment used in the experiments. The research
progress has been presented and discussed at several meetings, including
the International Conference on Complexity in Acute Illness of the Society
for Complexity in Acute Illness (SCAI). Comments from many SCAI members
greatly influenced this paper. We also thank the athletes who were the
subjects for this study. The theoretical aspects of this work and the
connections with other complex systems challenges were supported in part
by Air Force Office of Scientific Research and National Science Foundation.
Preliminary exploration in this research direction was funded by Pfizer,
National Institutes of Health (R01 GM078992), and the Institute of Collabo-
rative Biotechnologies (ARO W911NF-09-D-0001).
1. Chandra FA, Buzi G, Doyle JC (2011) Glycolytic oscillations and limits on robust effi-
ciency.
Science
333(6039):187
–
192.
2. Pool R (1989) Is it healthy to be chaotic?
Science
243(4891):604
–
607.
3. Buzsáki G (2006)
Rhythms of the Brain
(Oxford Univ Press, New York).
4. Glass L (2001) Synchronization and rhythmic processes in physiology.
Nature
410(6825):
277
–
284.
5. Buchman TG (2002) The community of the self.
Nature
420(6912):246
–
251.
6. Goldberger AL, et al. (2002) Fractal dynamics in physiology: alterations with disease
and aging.
Proc Natl Acad Sci USA
99(Suppl 1):2466
–
2472.
7. Camm AJ, et al.; Task Force of the European Society of Cardiology and the North
American Society of Pacing and Electrophysiology (1996) Heart rate variability: Standards
of measurement, physiological interpretation and clinical use.
Circulation
93(5):1043
–
1065.
8. Ivanov PC, et al. (1999) Multifractality in human heartbeat dynamics.
Nature
399(6735):461
–
465.
9. Poon CS, Merrill CK (1997) Decrease of cardiac chaos in congestive heart failure.
Nature
389(6650):492
–
495.
10. Moorman JR, et al. (2011) Cardiovascular oscillations at the bedside: Early diagnosis of
neonatal sepsis using heart rat
e characteristics monitoring.
Physiol Meas
32(11):1821
–
1832.
11. Tulppo MP, Mäkikallio TH, Takala TE, Seppänen T, Huikuri HV (1996) Quantitative
beat-to-beat analysis of heart rate dynamics during exercise.
Am J Physiol
271(1 Pt 2):
H244
–
H252.
12. Yamamoto Y, Hughson RL, Peterson JC (1991) Autonomic control of heart rate during
exercise studied by heart rate variability spectral analysis.
J Appl Physiol (1985)
71(3):
1136
–
1142.
13. Moorman JR, et al. (2011) Mortality reduction by heart rate characteristic monitoring
in very low birth weight neonates: A randomized trial.
J Pediatr
159(6):900
–
906.e1.
14. Akselrod S, et al. (1981) Power spectrum analysis of heart rate fluctuation: A quan-
titative probe of beat-to-beat cardiovascular control.
Science
213(4504):220
–
222.
15. Mackey MC, Glass L (1977) Oscillation and chaos in physiological control systems.
Science
197(4300):287
–
289.
16. Novak V, et al. (1993) Influence of respiration on heart rate and blood pressure
fluctuations.
J Appl Physiol
74(2):617
–
626.
17. Giera
ł
towski J,
_
Zebrowski JJ, Baranowski R (2012) Multiscale multifractal analysis of
heart rate variability recordings with a large number of occurrences of arrhythmia.
Phys Rev E Stat Nonlin Soft Matter Phys
85(2 Pt 1):021915.
18. Barabási AL (2009) Scale-free networks: A decade and beyond.
Science
325(5939):
412
–
413.
19. Glass L (2009) Introduction to controversial topics in nonlinear science: Is the normal
heart rate chaotic?
Chaos
19(2):028501.
20. Kluttig A, Kuss O, Greiser KH (2010) Ignoring lack of association of heart rate variability
with cardiovascular disease and risk factors: Response to the manuscript
“
The relation-
ship of autonomic imbalance, heart rate variability cardiovascular disease risk factors
”
by
Julian F. Thayer, Shelby S. Yamamoto, Jos F. Brosschot.
Int J Cardiol
145(2):375
–
376.
21. Young NS, Ioannidis JPA, Al-Ubaydli O (2008) Why current publication practices may
distort science.
PLoS Med
5(10):e201.
22. Stumpf MPH, Porter MA (2012) Mathema
tics. Critical truths about power laws.
Science
335(6069):665
–
666.
23. Alderson DL, Doyle JC (2010) Contrasting views of complexity and their implications
for network-centric infrastructures.
IEEE Transactions on Systems, Man and Cyber-
netics, Part A: Systems and Humans
40(4):839
–
852.
24. Ioannidis JPA (2005) Why most published research findings are false.
PLoS Med
2(8):e124.
E3484
|
www.pnas.org/cgi/doi/10.1073/pnas.1401883111
Li et al.
25. Wigertz O (1970) Dynamics of ventilation and heart rate in response to sinusoidal
work load in man.
J Appl Physiol
29(2):208
–
218.
26. Warner HR, Cox A (1962) A mathematical model of heart rate control by sympathetic
and vagus efferent information.
J Appl Physiol
17(2):349
–
355.
27. Brooks GA, Fahey TD, Baldwin K (2004)
Exercise Physiology: Human Bioenergetics
and Its Applications
(McGraw
–
Hill, New York).
28. Rowell LB (1993)
Human Cardiovascular Control
(Oxford Univ Press, New York).
29. Grodins FS, Buell J, Bart AJ (1967) Mathematical analysis and digital simulation of the
respiratory control system.
J Appl Physiol
22(2):260
–
276.
30. Guyton A, Hall J (2006)
Textbook of Medical Physiology
(Elsevier Saunders, Phila-
delphia), 11th Ed.
31. Hoppensteadt FC, Peskin CS (2002)
Modeling and Simulation in Medicine and the Life
Sciences
(Springer, New York).
32. Batzel JJ, Schneditz D (2007)
Cardiovascular and Respiratory Systems: Modeling,
Analysis, and Control
(Society for Industrial and A
pplied Mathematics, Phila-
delphia).
33. Ljung L (1999)
System Identification: Theory for the User
(Prentice Hall, Upper Saddle
River, NJ).
34. Kirk D (2004)
Optimal Control Theory: An Introduction
(Dover, New York).
35. Doya K, Ishii S, Pouget A, Rao RPN, eds. (2007)
Bayesian Brain: Probabilistic Ap-
proaches to Neural Coding
. (MIT Press, Cambridge, MA).
36. Cheng TM, Savkin AV, Celler BG, Su SW, Wang L (2008) Nonlinear modeling and
control of human heart rate response during exercise with various work load in-
tensities.
IEEE Trans Biomed Eng
55(11):2499
–
2508.
37. Mullen TJ, Appel ML, Mukkamala R, Mathias JM, Cohen RJ (1997) System identifica-
tion of closed-loop cardiovascular control: Effects of posture and autonomic block-
ade.
Am J Physiol
272(1 Pt 2):H448
–
H461.
38. Paulson OB, Strandgaard S, Edvinsson L (1990) Cerebral autoregulation.
Cerebrovasc
Brain Metab Rev
2(2):161
–
192.
39. Bramble DM, Lieberman DE (2004) Endurance running and the evolution of Homo.
Nature
432(7015):345
–
352.
40. Mehlsen J, Pagh K, Nielsen JS, Sestoft L, Nielsen SL (1987) Heart rate response to
breathing: Dependency upon breathing pattern.
Clin Physiol
7(2):115
–
124.
41. Hirsch JA, Bishop B (1981) Respiratory sinus arrhythmia in humans: How breathing
pattern modulates heart rate.
Am J Physiol
241(4):H620
–
H629.
42. Eckberg DL (1983) Human sinus arrhythmia as an index of vagal cardiac outflow.
J Appl Physiol
54(4):961
–
966.
43. Schäfer C, Rosenblum MG, Kurths J, Abel HH (1998) Heartbeat synchronized with
ventilation.
Nature
392(6673):239
–
240.
44. Grossman P, Karemaker J, Wieling W (1991) Prediction of tonic parasympathetic
cardiac control using respiratory sinus arrhythmia: The need for respiratory control.
Psychophysiology
28(2):201
–
216.
45. Berntson GG, Cacioppo JT, Quigley KS (1993) Respiratory sinus arrhythmia: Auto-
nomic origins, physiological mechanisms, and psychophysiological implications.
Psychophysiology
30(2):183
–
196.
46. Eckberg DL (2009) Point:counterpoint: Respiratory sinus arrhythmia is due to a central
mechanism vs. respiratory sinus arrhythmia is due to the baroreflex mechanism.
J Appl Physiol (1985)
106(5):1740
–
1742, discussion 1744.
47. Doyle JC, Csete M (2011) Architecture, constraints, and behavior.
Proc Natl Acad Sci
USA
108(Suppl 3):15624
–
15630.
48. Winslow RL, Trayanova N, Geman D, Miller MI (2012) Computational medicine:
Translating models to clinical care.
Sci Trans Med
4(158):158rv111.
49. Han K, Nagel JH, Schneiderman N (1992) A continuous representation of heart rate.
Engineering in Medicine and Biology S
ociety, 1992. 14th Annual International
Conference of the IEEE
2:785
–
786.
Li et al.
PNAS
|
Published online August 4, 2014
|
E3485
PHYSIOLOGY
PNAS PLUS