# Robust efficiency and actuator saturation explain healthy heart rate control and variability

^{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

See allHide authors and affiliations

Edited* by Michael S. Gazzaniga, University of California, Santa Barbara, CA, and approved June 27, 2014 (received for review January 30, 2014)

## 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 variability. Our analysis of HRV using the tools of control mathematics 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 understanding 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.

## Abstract

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 physiological, clinically relevant mechanisms have not been fully explained. Using mathematical tools from control theory, we combine mechanistic 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 underlying 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.

Biological systems display a variety of well-known rhythms in physiological signals (1⇓⇓⇓⇓–6), with particular patterns of variability 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 correlated with adverse states such as stress, fatigue, physiologic senescence, 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 understandable 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 minimal 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 experimental 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 indefinitely. 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 parasympathetic 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 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 explanation 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 measurement 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 system. The main control requirement is to maintain (*i*) small, acceptable “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 *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*, *R*_{s}, and cerebral autoregulation. Thus, healthy response involves low error but high control variability whereas loss of health is exactly the opposite. 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 physiological details, mathematical methods, and measurement technology. 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, scalable 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 (referred 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 relationships that may be hard to intuit from data alone and can be obscured in both complex data sets and mechanistic models, whereas first-principles models give physiological interpretations to these dynamical relationships.

## Results

### Static Fits.

Table 1 lists the minimum root-mean-square (rms) error ||*H_data-H_fit*|| (where *x*_{t} of length *N*) for several static and dynamic fits of increasing 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 ||*H_data-h*(*W)*||, which can be found easily by linear least squares. Static models have limited explanatory power but are simple starting points in which constraints 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 consequences 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 homeostasis 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 aerobic lipid metabolism in muscles and preserve precious carbohydrates for the brain.

The nonlinear response in Fig. 3 (solid lines) reflects additional 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 relatively 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 *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 autoregulation (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 reflected 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*:*P*_{as}, ∆*O*_{2}) = *f*(*H*, *W*) at 0 and 100 W. Here *P*_{as}, *∆O*_{2}, and *H* at different workloads. In particular, *P*_{as} on CBF due to saturation of autoregulation, and *P*_{as} and *H* at higher workload levels.

An important feature of this approach is that it allows systematic exploration of models that are both simple and explanatory. 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 behavior 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 extensibility 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-principles 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 workload (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)*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 mechanistically 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 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 workload 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 otherwise 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 experiments (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 cardiovascular control (27⇓⇓⇓–31) the neuroendocrine system controls peripheral vasodilation, minute ventilation, and cardiac output to maintain blood pressure and oxygen saturation within acceptable physiological limits.

Several features of these control systems allow substantial simplification of the model. Minute ventilation *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 consumption, we also assume that control of peripheral vascular resistance *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* and *P* are (subscripts *a* = arterial, *v* = venous, *s* = systemic, *P* = pulmonary) blood volume and blood pressure, respectively. 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 according 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:*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 circulatory compartment and the oxygen change in the tissue keeps the constraints from Eq. **2** but replaces **3** with the following dynamic model:_{2} volume and we assume that tissues and venous blood gases are in equilibrium, namely that tissue oxygenation *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 *P*_{as}, *∆O*_{2}, and *H*:**4**. Fig. 4 compares HR and workload data versus simulations of applying the linear controllers to the model in **4** for two experiments (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 Appendix*.). Also note that for different subjects, we use the same nominal values for most parameters except the constants *c*_{r},*c*_{l} used in the cardio-output formulation and the weighting parameters **5**. The different *c*_{r},*c*_{l} reflect the different sizes of subjects and different

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 damaging *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 required 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 between 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 explanations 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 systems 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 additional 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 optimal static fit *H* = *h*(*V*) = *b·V* + *c*, where *V* (blue) is measured ventilation flow rates (inhalation and exhalation) at the mouthpiece. 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, because 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):*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 physiological 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 oscillations 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 consistent 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 measurements 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 mechanical and metabolic effects is needed.

## Discussion

### Robust Efficiency and Actuator Saturation.

We showed how HR fluctuations in healthy athletes can be largely explained as nonlinear 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 measured, are much less important mechanistically than the tradeoffs that produce them. The tradeoffs we highlight between robustness, 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 physiological (HR) control system include how efficiently maintaining 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, ventilation, and cerebral autoregulation) responses to compensate, and how nonlinearities in actuator saturations 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 physiology 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 important 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, hemorrhagic 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 autonomous robot car is in the control system for robust efficiency to uncertainty in real traffic environments and to intrinsic variability in components. If scientists were forced to “reverse engineer” 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 sophisticated 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 symptom 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 antislip traction control systems mitigate these effects, and augment human control in emergencies. However, even in fully automated 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 dangerous 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, mechanistic 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 regulatory 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 accidents of evolution. Oscillations are neither functional nor accidental 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 fluctuations 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 balance 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 uninformative. In contrast, the static and dynamic models mechanistically 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 physiology. The main requirement of the models is some mechanistic relationship between control actuation and its limits in maintaining 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 physiological 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 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 straightforward study should serve well as a foundational framework for the evaluation of even more complex physiologic (disease) situations in which the diagnostic possibilities are broader.

### Clinical Correlates: Linking the Behavior of Control Systems and Pathophysiology.

Clinicians know that changes in actuator signals (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 characteristics 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 phenomenon not observed in adults. Furthermore, because sepsis is a common problem in this setting, the clinicians chose to administer antimicrobial therapy on the basis of the alert which produced the outcome benefit. In fact, the efficacy of the monitor 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 experimental 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 resistance controlled by a preprogrammed protocol. In the respiration rate (RR) sweep experiments, a sinusoid signal was preprogrammed in the computer 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 protocol. Other exercise testing data were collected using commercially available 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 ventilation flow rate data. (*iii*) In the test shown in Fig. 7, 1-Hz gas data including minute ventilation

A detailed discussion of mathematical methods is given in *SI Appendix*.

## Acknowledgments

We thank Pamela B. Pesenti for her gift in establishing 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 Collaborative Biotechnologies (ARO W911NF-09-D-0001).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: doyle{at}caltech.edu.

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.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1401883111/-/DCSupplemental.

## References

- ↵.
- Chandra FA,
- Buzi G,
- Doyle JC

- ↵.
- Pool R

- ↵.
- Buzsáki G

- ↵
- ↵
- ↵.
- Goldberger AL,
- et al.

- ↵.
- Camm AJ,
- et al.,
- Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology

- ↵
- ↵
- ↵
- ↵.
- Tulppo MP,
- Mäkikallio TH,
- Takala TE,
- Seppänen T,
- Huikuri HV

- ↵.
- Yamamoto Y,
- Hughson RL,
- Peterson JC

- ↵
- ↵.
- Akselrod S,
- et al.

- ↵.
- Mackey MC,
- Glass L

- ↵.
- Novak V,
- et al.

- ↵
- ↵.
- Barabási AL

- ↵
- ↵.
- Kluttig A,
- Kuss O,
- Greiser KH

- ↵
- ↵.
- Stumpf MPH,
- Porter MA

- ↵
- ↵
- ↵.
- Wigertz O

- ↵.
- Warner HR,
- Cox A

- ↵.
- Brooks GA,
- Fahey TD,
- Baldwin K

- ↵.
- Rowell LB

- ↵.
- Grodins FS,
- Buell J,
- Bart AJ

- ↵.
- Guyton A,
- Hall J

- ↵.
- Hoppensteadt FC,
- Peskin CS

- ↵.
- Batzel JJ,
- Schneditz D

- ↵.
- Ljung L

- ↵.
- Kirk D

- ↵Doya K, Ishii S, Pouget A, Rao RPN, eds. (2007).
*Bayesian Brain: Probabilistic Approaches to Neural Coding*. (MIT Press, Cambridge, MA). - ↵
- ↵.
- Mullen TJ,
- Appel ML,
- Mukkamala R,
- Mathias JM,
- Cohen RJ

- ↵
- ↵
- ↵
- ↵.
- Hirsch JA,
- Bishop B

- ↵.
- Eckberg DL

- ↵
- ↵
- ↵
- ↵.
- Eckberg DL

- ↵.
- Doyle JC,
- Csete M

- ↵.
- Winslow RL,
- Trayanova N,
- Geman D,
- Miller MI

- ↵.
- Han K,
- Nagel JH,
- Schneiderman N

*Engineering in Medicine and Biology Society, 1992. 14th Annual International Conference of the IEEE*2:785–786.

## Citation Manager Formats

## Article Classifications

- Biological Sciences
- Physiology