## Abstract

During orthostatic stress, arterial and cardiopulmonary baroreflexes play a key role in maintaining arterial pressure by regulating heart rate. This study presents a mathematical model that can predict the dynamics of heart rate regulation in response to postural change from sitting to standing. The model uses blood pressure measured in the finger as an input to model heart rate dynamics in response to changes in baroreceptor nerve firing rate, sympathetic and parasympathetic responses, vestibulo-sympathetic reflex, and concentrations of norepinephrine and acetylcholine. We formulate an inverse least squares problem for parameter estimation and successfully demonstrate that our mathematical model can accurately predict heart rate dynamics observed in data obtained from healthy young, healthy elderly, and hypertensive elderly subjects. One of our key findings indicates that, to successfully validate our model against clinical data, it is necessary to include the vestibulo-sympathetic reflex. Furthermore, our model reveals that the transfer between the nerve firing and blood pressure is nonlinear and follows a hysteresis curve. In healthy young people, the hysteresis loop is wide, whereas, in healthy and hypertensive elderly people, the hysteresis loop shifts to higher blood pressure values, and its area is diminished. Finally, for hypertensive elderly people, the hysteresis loop is generally not closed, indicating that, during postural change from sitting to standing, baroreflex modulation does not return to steady state during the first minute of standing.

- mathematical modeling
- heart rate control
- baroreflex function
- sympathetic and parasympathetic responses
- vestibulo-sympathetic reflex

baroreflexes play a significant role in short-term cardiovascular control in adaptation to orthostatic stress. Postural change is a physiological stimulus that activates the baroreflex feedback control system. Standing up is clinically relevant for elderly people who are at risk for orthostatic hypotension due to baroreflex impairment. Blood pressure rapidly decreases upon standing and then increases to or above the baseline within the first minute of standing. Heart rate increases in response to a blood pressure decline and decreases as a response to a blood pressure increase. These responses are often diminished in elderly people and in people with cardiovascular diseases.

Current methodologies to assess baroreflex sensitivity include both time and frequency domain techniques. In the time domain (45), baroreflex sensitivity is calculated by relating the change in heart period to beat-to-beat blood pressure or muscle sympathetic nerve activity (MSNA) in response to physiological (e.g., postural change) and pharmacological stimuli (e.g., injection of vasoactive agents). In the frequency domain, baroreflex sensitivity is determined from spectral analysis of spontaneous fluctuations in heartbeat and blood pressure, which are thought to reflect parasympathetic and sympathetic responses (30, 48).

With aging, baroreflex sensitivity (24) and vagally mediated cardiac variability (9) become attenuated, and sympathetically mediated fluctuations in vascular tone become prominent. Chronic elevation of sympathetic tone and blood pressure in hypertension is associated with further impairment of baroreflex function, reduction of heart rate variability, increased vascular resistance, and the shift of baroreflex operating range toward higher blood pressure values (3, 16). Endothelial dysfunction and circulating vasoconstrictor factors contribute significantly to decreased baroreflex sensitivity and disinhibition of sympathetic activity that occurs with aging and hypertension (7).

From observations discussed above, it becomes clear that baroreflex feedback control plays a crucial role in heart rate regulation, and that development of reliable measures for baroreflex function is needed to assess the integrity of the autonomic nervous system. The main difficulty in assessing autonomic function based on heart rate and blood pressure data is that an inverse problem must be solved, because only the output (heart rate) can be measured while internal variables, i.e., baroreceptor firing rate, sympathetic and parasympathetic tone, or concentrations of neural hormones catecholamine (epinephrine and norepinephrine) and acetylcholine have to be estimated. Most of the previous studies mentioned above use either spectral analysis or simple linear regression to assess autonomic function. While these simple models may provide valuable insight into heart rate regulation, the effects of nonlinear feedback control and the contributions from each of the intermediate steps remain unknown. A few studies use mathematically more advanced methods that are based on nonlinear control theory and nonlinear differential equations (20, 31, 33, 35, 40). However, these models were not based on well-established physiological theory, and none of these models was validated against clinical data.

To study the relation between heart rate and intermediate controls, we developed a nonlinear mathematical model that explicitly represents each component of the baroreflex control system, and we used this model to study effects of aging and hypertension on baroreflex feedback control of heart rate dynamics. This study is focused on modeling dynamics of heart rate regulation; however, the modeling philosophy and technique as presented here can be applied to more complex cardiovascular models that incorporate all mechanisms involved in short-term blood pressure and blood flow velocity regulation (32). That is, we modeled heart rate dynamics as a function of blood pressure, but ignored the feedback that heart rate changes and regulation of peripheral vascular resistance, vascular tone, and cardiac contractility have on blood pressure. Furthermore, our model does not include effects of respiration, and, hence, it does not account for cardiopulmonary reflexes, which are known to play a role during heart rate regulation.

Our mathematical model for baroreflex feedback control of heart rate responses to postural change is divided into four submodels connected in series (see Fig. 1). The first submodel is an afferent trigger model, which uses blood pressure measured in the finger as an input to predict the firing rates of baroreflex afferent fibers. For the remainder of this study, we assume that the finger blood pressure can be used in place of carotid blood pressure, which was not measured in this study. This assumption is reasonable, since the finger and the carotid pressure waveforms have similar shape, even though the amplitude of the finger pressure may be higher and the mean finger pressure may be lower than the carotid pressure. Furthermore, we assume that the arterial wall mainly displays elastic behavior. Hence, the change of stretch of the arterial wall is proportional to the change in blood pressure.

The second submodel, which represents the central nervous system, uses baroreceptor afferent nerve activity as an input to predict sympathetic and parasympathetic firing in response to the rate of change of the mean blood pressure. Within one to two cardiac cycles following postural change, the parasympathetic response is decreased, leading to an immediate increase in heart rate. This response is mediated through a decreased release of acetylcholine at the postsynaptic nerve terminal, which leads to cardiac acceleration. Subsequently, the sympathetic firing rate increase leads to an increase in norepinephrine release, which, in turn, increases heart rate. For healthy young people, the sympathetic response is delayed 6–10 s following postural change, whereas for healthy and hypertensive elderly people the sympathetic response may be further delayed. In our mathematical model, we explicitly incorporate the delay in the sympathetic response. This submodel also includes the mechanisms underlying the initial heart rate increase that precedes the blood pressure decline observed during standing. It is not clear what is the main mediator of this early response. A study by Borst et al. (2) attributed this initial heart rate increase to the response of an exercise reflex, which may be mediated by central command and activation of MSNA. A recent study (49) has shown that vestibular input to the otolith organs modulates muscle MSNA and that this activation is related to blood pressure changes. An additional study (21) during head rotation with tilt showed that the MSNA reflex might be one of the earliest mechanisms to sustain blood pressure upon standing. Similar observations have been found during head-down rotation (38) and during climbing (53). A common conclusion drawn in these studies is that the vestibulo-sympathetic reflex contributes to blood pressure maintenance, and, because of its short latency, this reflex may be one of the earliest mechanisms that maintain blood pressure upon standing, and it precedes the baroreflex-mediated response. The studies summarized above indicate that vestibulo-sympathetic reflex has an impact on initial heart rate regulation; however, other factors, in particular central command, may also play a significant role. The main effect of these mechanisms is their contribution to the increase in heart rate that precedes the baroreflex-mediated response. Because of the contributions from the vestibular nerves and from central command, we chose to model this initial heart rate increase by adding an impulse function to the baroreflex-mediated sympathetic response. We could possibly have modeled the vestibular feedback using a more complex model. However, at the present, sufficient biological information is not readily available to develop a physiologically based model. In the remainder of the paper, we will treat this early control as a vestibulo-sympathetic response, but we do note that it is possible that the response has some component of central command and cardiopulmonary reflexes embedded in the response.

The third submodel uses sympathetic and parasympathetic responses as an input to predict concentrations of the neurotransmitters norepinephrine and acetylcholine. The fourth submodel, the effector model, uses concentrations of neurotransmitters as an input to predict heart rate. Finally, we provide quantitative comparisons of our model with physiological heart rate data from healthy young, healthy elderly, and hypertensive elderly people during postural change from sitting to standing.

## METHODS

### Data Collection

Model validation used data from 30 people in three groups: *1*) 10 healthy men and women ages 20–40 yr; *2*) 10 healthy men and women ages 60–80 yr with no known systemic disease, no cardiovascular medications, no history of head or brain injury, and no history of more than one episode of syncope; and *3*) 10 men and women ages 60–80 yr, with a diagnosis of systolic hypertension (with a systolic blood pressure >160 mmHg), who were not previously treated for hypertension. All data collection procedures were reviewed and approved by the institutional review board of Hebrew SeniorLife, and all subjects provided written informed consent.

Each subject was instrumented with a three-lead ECG to obtain heart rate. A photoplethysmographic device on the middle finger of the nondominant hand was used to obtain noninvasive beat-to-beat blood pressure (Finapres device, Ohmeda Monitoring Systems, Englewood, CO). To eliminate effects of gravity, the hand was held at the level of the right atrium and supported by a sling. All physiological signals were digitized at 50 Hz (Windaq, Dataq Instruments) and stored for offline analysis.

After instrumentation, subjects sat in a straight-backed chair with their legs elevated at 90° in front of them. After 1 min of stable signals, the subjects were asked to stand. Standing was defined as the moment both feet touched the floor. Two 5-min sit-to-stand trials were performed for each subject. These data were used for validation of our mathematical model, but have been published earlier (see Ref. 29).

### Modeling Baroreflex Feedback Control of Heart Rate

The baroreflex feedback control maintains a normal blood pressure by regulating heart rate and vascular resistance in response to transient changes in stretch of arterial baroreceptors (17). During postural change from sitting to standing, ∼500 ml (19) of blood are pooled in the legs as the result of gravitational force. As a consequence, there is a transient decline in aortic pressure and in cardiac output. To compensate for these changes, the frequency of afferent impulses in the aortic and carotid sinus nerves is reduced, leading to parasympathetic withdrawal and sympathetic activation. Throughout this paper, the nerve activity will be referred to as the baroreceptor firing rate, or simply the firing rate. Sympathetic activation leads to an increased release of neural hormonal catecholamines (norepinephrine) and by release of epinephrine from the adrenal gland, which contribute to restoration of blood pressure by increasing heart rate, cardiac contractility, vasoconstrictor, and venoconstrictor tone. In addition, descending impulses from the cerebral cortex and central autonomic network further modulate adaptation to upright posture by reducing parasympathetic outflow. This is a rapid process capable of cardiac acceleration within a heartbeat.

In this section, we describe a comprehensive mathematical model for the baroreflex regulation of heart rate. The basic structure of the model (see Fig. 1) consists of the following components: the arterial blood pressure (the input), the afferent baroreceptor nerve fibers going to the central nervous system, the central nervous system, the efferent sympathetic and parasympathetic nerve fibers going to the heart, and the effects of neurotransmitters (norepinephrine and acetylcholine) on heart rate.

#### Afferent baroreceptor activity.

This part of our model describes the dynamics of the baroreceptors afferent signaling to a change in arterial blood pressure. The baroreceptors are sensitive to stretch of the carotid arterial wall, which is caused by changes in arterial pressure (18, 28). Due to the complex composition of the arterial wall, experiments that studied baroreceptor response to a change in arterial pressure show several nonlinear characteristics (5, 43). The frequency of firing rate increases with increased arterial pressure. There is a wide variation in threshold for different receptors. Sufficiently fast decreases in pressure cause a complete cessation of firing rate, which, after a few seconds, begins to discharge again at a frequency characteristic of the new pressure level. Finally, the firing rate response curve for elderly people and, in particular, for people with hypertension translates to the right along the pressure axis, compared with curves for healthy young people (8). Despite the above-mentioned knowledge of the baroreceptors' nonlinear characteristics, the mechanoelectrical transduction that takes place in the baroreceptors themselves is not well known (5, 6, 23). Whether it is the instantaneous pressure, an average pressure, or a combination of both that triggers the baroreceptors' afferent signaling is not clear from measurements (47). Earlier, there has been some effort to model the nonlinear characteristics of the baroreceptors using a different assumption on various forms of input. For example, Landgren (25) and Robinson and Sleight (41) suggested simple functional descriptions of the response to a step change in pressure only. Other studies (11, 35, 36, 41, 42, 47, 50, 51) suggested various models based on ordinary differential equations. Our model, including initial choices for parameter values, follows our earlier work described in Refs. 34–36. More specifically, we model the nonlinear response in firing rate *n* as a function of the rate of change of mean arterial pressure using the following system of ordinary differential equations (see Table 1) (1) where *i* is short S, intermediate I, and long L, and *n* = *n*_{S} + *n*_{I} + *n*_{L} + *N*; *M* is the maximum firing rate; *p̄* is pressure; *t* is time; and τ_{i} is the time constant associated with *n*_{i}. In the above equation, the change in nerve firing rate is directly proportional to the rate of change of the mean arterial pressure d*p̄*/d*t*. To account for the wide variation in thresholds for the different receptors, three characteristic time scales were included: short [≈1 (s)], intermediate [≈5 (s)], and long [≈250 (s)], represented by the time constants τ_{S}, τ_{I}, and τ_{L}, respectively. The arterial wall consists of viscoelastic tissue embedded with nerve fibers. The relaxation and stretch of the arterial wall involves different mechanisms (14); thus in *Eq. 1* it is assumed that the arterial walls constrict more easily than they dilate. Hence the response in firing rate occurs faster during pressure decrease than during pressure increase. The change of sign in the d*p̄*/d*t* term in *Eq. 1*, which corresponds to the increase and decrease in the mean pressure, accounts for this hysteresis effect. The variables *n*_{i}, *i* = S, I, L, represent deviations from the baseline firing rate *N*. Finally, the first term on the right-hand side, *k*_{i} (d*p̄*/d*t*)*n*(*M* − *n*)/(*M*/2)^{2}, implies that the effect of the rate of change in the mean pressure in firing rate lies between the physiological threshold values *n* = 0 and *n* = *M*, respectively. The maximum firing rate *M* is chosen to be 120 in all simulation studies. Values for the parameters τ_{S}, τ_{I}, τ_{L}, *k*_{S}, *k*_{I}, *K*_{L}, and *M* are based on estimates for animals reported in the literature (4, 12, 13, 25, 41).

Because the instantaneous pressure oscillates from beat to beat, our motivation for using the mean arterial pressure is to obtain stable numerical simulation results. Using ideas from Ref. 32, mean pressure values are computed as weighted averages, where the present is weighted higher than the past, according to the following expression (2) where α (1/s) is the history weighing parameter, which accounts for the number of cardiac cycles included in the calculation of mean pressure, i.e., a small value of α gives rise to a slow exponential decay, and thus more cardiac cycles are included in the calculation of mean pressure, while a large value of α gives rise to a fast exponential decay, and thus only a few cardiac cycles are included in the mean pressure calculation. To incorporate the mean arterial blood pressure into the rest of our mathematical model, which is described by ordinary differential equations, we differentiate the integral *Eq. 2* to obtain a differential equation for the mean pressure (3)

#### Central nervous system.

This submodel describes effects of the afferent baroreceptor nerve activity on efferent parasympathetic and sympathetic responses. This model is empirical, but it is based on well-known experimental facts. The parasympathetic response *T*_{par} is known to follow the “direct law”, i.e. (4) whereas the sympathetic response *T*_{sym} follows the reciprocal law (34, 35). In addition, experimental studies have revealed that there is a time delay of between 6 and 10 s for the peak response to appear in the sympathetic nervous system and almost instantaneous in the parasympathetic nervous system (50, 52). This is accounted for by introducing a time delay τ_{d}. Furthermore, the parasympathetic response has an inhibitory effect on the sympathetic response (27). Finally, experiments in humans indicated that, during postural change, the vestibulo-sympathetic reflex acts to defend against a possible sustained hypotensive episode before a drop in arterial pressure is sensed by the baroreceptors (21, 38, 39, 49, 53). It should be emphasized that the central command can contribute to this response. To account for these effects, we have modeled the sympathetic response as (5) where β represents the dampening factor of the parasympathetic response, and *u*(*t*) is the impulse function accounting for the regulation of the sympathetic nerve activity by the vestibulo-sympathetic system. This impulse response is modeled as a parabolic function given by (6) The parameters *t*_{start} and *t*_{stop} are the start and stop time of the impulse, and *u*_{0} is the amplitude of the response. This empirically based model for the vestibulo-sympathetic reflex is designed so that its maximal effect on the sympathetic response occurs halfway between the onset and the end of the vestibulo-sympathetic reflex. It is important to note that we were not able to accurately predict measured heart rate without inclusion of this impulse response function.

#### Neurotransmitters norepinephrine and acetylcholine.

Presynaptic and postsynaptic regulation of cardiac response is modulated by several neurotransmitters that have inhibitory, excitatory, or both effects on cardiac function. Neurotransmitter release is not limited to centrally mediated neural traffic but may be triggered in response toneurotransmitters and paracrine substances from blood or nearby tissue. In this model, we have limited the parasympathetic neurotransmitter response to acetylcholine and of the sympathetic response to norepinephrine and their effects on heart rate. We have not assessed the interactions between them, nor the effects of other substances that may affect the heart rate. We model concentrations of norepinephrine C_{nor} and acetylcholine C_{ach} as linear differential equations with forcing terms that depend on the sympathetic (T_{sym}) and parasympathetic (T_{par}) responses. Equations for concentrations are kept on nondimensional form to reduce the total number of parameters to be predicted, (7) where τ_{nor} and τ_{ach} are time constants.

#### Heart rate.

As a response to the decrease in blood pressure, heart rate is increased by the sympathetic system through an increased release of norepinephrine and by parasympathetic withdrawal through a decreased release of acetylcholine. In this study, we modeled heart rate in response to these chemical concentrations using an integrate-and-fire model of the form (8) where *M*_{S} and *M*_{P} are scaling factors for sympathetic and parasympathetic response, respectively, and H_{0} = 100 (beats/min) is the intrinsic heart rate when denervated, i.e., without any sympathetic or parasympathetic stimulation (1, 10, 15, 17). In this model, a heartbeat occurs every time the time-dependent function φ passes the value 1. When this event occurs, φ is reset to zero. If the heart beats at consecutive times *t*_{i} and *t*_{i+1}, then the heart rate (HR) is given by HR = 1/(*t*_{i+1} − *t*_{i}) (see Fig. 2).

### Parameter Estimation

The mathematical models for heart rate regulation described in *Eqs. 1*–*8* contain a total of 16 unknown parameters denoted by *q⃗* = (*k*_{S}, *k*_{I}, *k*_{L}, τ_{S}, τ_{I}, τ_{L}, *N*, α, β, *u*_{0}, *t*_{start}, *t*_{stop}, τ_{nor}, τ_{ach}, *M*_{S}, *M*_{P})^{T}. To estimate these parameters from the experimental data, we formulate an inverse least squares problem minimizing the difference between the measured and computed values of heart rate. That is, we seek to find *q⃗* to minimize the cost functional (*J*) (9) where *N*_{d} is the number of measurements, HR^{d}(*t*_{i}) are the measured values of heart rate obtained at time *t*_{i}, and HR^{m}(*t*_{i}) are heart rate values obtained from solving the mathematical models (*Eqs. 1*–*8*) at the same times where the data are recorded.

The differential equations are solved using MATLAB's (The MathWorks, Natick, MA) solver ode15s, which is a variable-order, variable-step solver for stiff ordinary differential equations. It is based on the numerical differentiation formulas with the option of using the backward differentiation formula, also known as Gear's method. As it is a variable-step solver, the numerical solutions are, in general, computed at times that are, in general, not the same as those of the experimental data. Hence, to obtain computed values of HR at the same times where experimental data are recorded, numerical interpolation is required to evaluate the computed value at the same temporal points where the experimental data were measured. Finally, because our mathematical models are divided into four parts in series, they are solved sequentially with the output of one part being used as the input for the next submodel, as shown in Fig. 1. However, it should be emphasized that, during the parameter estimation process, the model was considered as one integrated component. That is, all unknown parameters in the model were identified simultaneously.

For biological reasons, three of the unknown parameters, *N*, *M*_{S}, and *M*_{P}, are further constrained by upper and lower bounds. The resting firing rate *N* cannot be larger than the maximum firing rate *M*. Furthermore, we assume that it cannot be smaller than *M*/2. *M*_{S} and *M*_{P} are both assumed to be bounded by 0 and 1 (0 ≤ *M*_{S}, *M*_{P} ≤ 1). These bounds on the parameters imply that the optimization in *Eq. 9* now becomes a constrained minimization problem. To avoid solving a constrained minimization problem, which is computationally more expensive and difficult than an unconstrained problem, we parameterize the constrained parameters *N*, *M*_{S}, and *M*_{P} as follows (10) where 0 ≤ η, ξ_{S}, ξ_{P} < ∞. That is, we optimized the parameterized constants η, ξ_{S}, and ξ_{P} (instead of *N*, *M*_{S}, and *M*_{P}). For biological reasons, all unknown parameters are assumed to be positive. This positive constraint in the parameters can be enforced mathematically by replacing the parameters with their squared values in the mathematical models.

Optimization methods for unconstrained minimization problems fall into two classes: gradient-based methods and sampling methods. It is well known that gradient-based methods perform very well when the optimization landscape is relatively smooth (22). However, for many applications, including heart rate regulation, the nonlinear interaction of biological mechanisms can give rise to nonsmoothness and nonconvexity in the landscape, which can defeat most gradient-based methods. In this paper, we employ the Nelder Mead method developed by Kelley (22), which does not require gradient information, but rather samples the objective function (10) on a stencil or pattern to determine the progress of the iteration and whether or not to change the size, but not the shape of the stencil. In essence, given *N*_{p} number of unknown parameters to be estimated, the Nelder Mead method will form a simplex with *N*_{p} + 1 vertices. It then attempts to minimize the objective function by replacing the simplex point that has the highest function value (i.e., the worst point), with a new vertex point with a lower function value. This process is repeated until the simplex region becomes sufficiently small. For a detailed treatment of the Nelder Mead algorithm, see Ref. 22.

#### Initial conditions and initial parameter values.

The nerve firing rates, *n*_{i}(0) = 0, *i* = S, I, L, are assumed to be zero initially so that the total firing rate *n*(0) = *N*, which is the baseline firing rate at rest. The initial condition for the mean arterial pressure *p̄*(0) is calculated as a mean of the experimental data over the first 10 s. During steady state, there is no change in concentrations of norepinephrine and acetylcholine, dC_{nor}/d*t* = dC_{ach}/d*t* = 0. Consequently, the differential equations in *Eq. 7* give C_{nor}(0) = T_{sym}(0) = (1 − *N*/*M*)/(1 + β*N*/*M*), and C_{ach} = T_{par} = *N*/*M*. Simulations start at the beginning of the cardiac cycle, thus φ(0) = 0.

Since the Nelder Mead algorithm is an iterative method, initial guesses for the 16 parameters must be specified. Whenever possible, we used the literature as a guidance for initial parameter values, From Ref. 35 we get *k*_{S} = *k*_{L} = 2, *k*_{I} = 1.5 (Hz/mmHg), and τ_{S} = 0.5, τ_{I} = 5, τ_{L} = 250 (s). Since the reaction rates and the history-weighting parameter are not generally well known, we used τ_{ach} = τ_{nor} = 0.05 (s) and α = 1 (1/s). Sympathetic response is delayed ∼6–10 s, and, hence, we used an initial value of τ_{d} = 7 (s). An initial value of β = 6 was chosen for the dampening factor in *Eq. 6*. To calculate the initial values for the scaling factors *M*_{S} and *M*_{P} in the heart rate model (8), we made the following assumptions. We assumed that, at *t* = 0, dφ/d*t* = 45 (beats/min), which is the intrinsic heart rate, H_{0} = 100 (beats/min), that C_{nor} = 0, and that C_{ach} = 1. Substituting these expressions into the heart rate *Eq. 8* gives (11) Similarly, when dφ/d*t* = H_{max} (the maximal firing rate), C_{nor} = 1 and C_{ach} = 0. This implies that (12) The maximum heart rate H_{max} depends on age and is computed as H_{max} = 217 − (0.85) × age. As discussed above, it should be emphasized that we optimized the parameters ξ_{S} and ξ_{P} (instead of *M*_{S} and *M*_{P}). Initial values for ξ_{S} and ξ_{P} are computed from initial values for *M*_{S} and *M*_{P} using *Eq. 10*. To calculate the initial value for the baseline firing rate *N*, which is subject dependent, we calculate the mean value of the first five heart rate data points for each individual and denoted this mean value by dφ/d*t*. Substituting this mean value into the heart rate differential *Eq. 8* we obtain (13) where the maximum firing rate *M* = 120 (Hz). Using *Eq. 13*, it is possible to solve for *N* to obtain an initial value for the baseline sympathetic firing rate, which, in turn, is used to obtain an initial value for the optimized parameter η used in *Eq. 10*. Finally, for the impulse function, we used *u*_{0} = 1, *t*_{start} = 58, and *t*_{stop} = 63 (s) as initial values.

## RESULTS

We have validated the model against 60 data sets, i.e., based on the initial estimates for all parameters discussed above, we used the Nelder Mead nonlinear optimization to minimize the discrepancy between computed and measured values for heart rate. This minimization was performed for each data set; as a result, for each group of subjects, we had 20 values for each parameter. Mean values and standard deviations are given in Table 2. From the heart rate data alone, it is not easy to separate the biological mechanisms controlling heart rate regulation, such as parasympathetic withdrawal, sympathetic activation, and sympathetic regulation by the vestibulo-sympathetic system. However, from the results of our submodels to be discussed in this section, these mechanisms are readily identified. Figure 2 shows typical blood pressures (input to the model) for a healthy young (*A*), healthy elderly (*B*), and hypertensive elderly (*C*) subject. Graphs show measured pulsatile (solid lines) and mean (gray lines) blood pressure as a function of time. The mean blood pressure is computed as described in *Eq. 2*. The subject stands up at *t* = 60 (s), which is marked by a dotted vertical line. Shortly after standing, blood is pooled in the legs, and, as a result, blood pressure (systolic, diastolic, and mean value) drops. Blood pressure regulation restores blood pressure ∼20 (s) after standing. It is noted that the blood pressure is significantly higher for the hypertensive subject (Fig. 2*C*) than for the healthy subjects (*A* and *B*).

Figure 3 shows the baroreceptor firing rate modeled as described in *Eq. 1*. Note that the change in firing rate dynamics is greater in the young subject (*A*) than in the elderly subjects (*B* and *C*), and that the total firing rate is reduced for the hypertensive subject (compare *A* and *B* with *C*). Furthermore, it should be noted that the short-term response is almost negligible for the elderly subjects (*B* and *C*) compared with the young subject.

Figure 4 shows the baroreceptor nerve firing rate computed from *Eq. 1* as a function of mean blood pressure. For the young subject (*A*), these graphs clearly show hysteresis effects, where the baroreceptors firing rate takes the lower path when the pressure is decreasing and the upper path when the pressure is increasing. For the healthy (*B*) and hypertensive (*C*) elderly subjects, the hysteresis curves shift toward higher pressures (in particular for the hypertensive subject), and the areas enclosed by the hysteresis curves are diminished, particularly for the healthy elderly people. The areas of the hysteresis curves are given in Table 3. These areas are computed by fitting a sixth-order polynomial to the hysteresis curves and then by integrating between the resulting polynomial functions. Figure 4 also shows the hysteresis curves, the fitted polynomials, and a slope. The slope is computed by placing a linear curve between the top and bottom of the hysteresis curve. The top point marks where pressure is recovered after standing, and the bottom point is placed at the minimum value for both pressure and area.

We used ANOVA analysis to compare the areas of the hysteresis curves among the three groups (see Table 3). The hysteresis area was largest for the young subjects [healthy young vs. healthy elderly (*P* < 0.000001), healthy young vs. hypertensive elderly (*P* < 0.0002), and healthy elderly vs. hypertensive elderly (*P* < 0.006)]. Baroreflex sensitivity is often assessed as linear fit between heart rate and blood pressure (44). This is similar to the slopes shown in Fig. 4, which were different between healthy young and healthy elderly (*P* < 0.006), and between healthy young and hypertensive elderly (*P* < 0.00001), while the difference was borderline between the groups of healthy elderly and hypertensive elderly people (*P* = 0.057). In addition, for some hypertensive elderly subjects, the hysteresis curve is open (Fig. 4*C*), indicating that baroreflex modulation does not return to the steady-state value during the first minute of standing.

The graphs of the parasympathetic and sympathetic responses are plotted as function of time (see Fig. 5) and as a function of the baroreceptor nerve firing rate (see Fig. 6). Since parasympathetic and sympathetic responses are affected by the baroreceptor nerve activity, their dynamics and characteristics are similar to the baroreceptor firing rate (see Fig. 3). Note in Fig. 5*A*, the sympathetic responses are bimodal. The first peak represents the regulation of sympathetic nerve activity by the vestibulo-sympathetic system, and the second peak represents the baroreflex-mediated sympathetic outflow. The vestibulo-sympathetic responses are indeed preceding the decline in parasympathetic response due to the drop in the arterial blood pressure (see also Fig. 7, which shows the heart rate is increased before the drop in arterial blood pressure). In our model, we accounted for the vestibulo-sympathetic reflex activation by adding the impulse function in *Eq. 6* to the sympathetic relation in *Eq. 5*. The onset and duration of the impulse function are treated as unknown parameters. The fact that our model predicts that the impulse response (and thus the vestibulo-sympathetic reflex) is activated before the onset of pressure decrease is really remarkable and, in fact, agrees with experimental studies (21, 38, 49, 53). It is noted that the responses of these mechanisms are significantly reduced for elderly subjects. The damping effect of the sympathetic response by the parasympathetic responses is also visible (see Fig. 5*B*), as the sympathetic response increases more slowly than the parasympathetic increases. The linear relationships of the parasympathetic and sympathetic responses to the baroreceptor firing rates are clearly illustrated in Fig. 6. As nerve firing decreases, the sympathetic response increases and the parasympathetic response decreases linearly.

A sample of heart rate model predictions compared with experimental data for the three groups of subjects is shown in Fig. 8. The graphs clearly show that we obtain excellent agreement between model predictions and measured data. Our mathematical model shows an increase in heart rate from the vestibulo-sympathetic reflex (before the blood pressure starts to decrease), a continuous heart rate increase in response to the parasympathetic withdrawal, and a steep increase in heart rate following the delayed sympathetic response. This is followed by a steep return to a higher resting heart rate when blood pressure is successfully raised and heart rate regulation is completed. The cost function in *Eq. 9* is a measure of how well our model predicts the data (see Table 2). The largest discrepancy between the model and the data (the largest values of *J =* 15.2 ± 6.8) was found for the healthy young subjects, while, for the elderly subjects, the model predicted the data very well (*J =* 3.36 ± 1.9 for healthy elderly subjects and *J =* 4.52 ± 2.76 for hypertensive elderly subjects). The larger cost (error) for the young subjects can be explained from the fact that heart rate variability is larger for young subjects than for elderly subjects (see Fig. 8, data for young subjects display more oscillations than data for the elderly subjects).

The means of the estimated parameter values together with their standard deviations are given in Table 2, these mean values are close to the initial values calculated from physiological observations, as described in the previous section. For each group of parameters, we used ANOVA analysis to compute the *P* values to compare the parameters among the three groups of subjects. Parameters that were statistically significant among the groups include the following: τ_{d}, which increased with age and even more with hypertension; *k*_{I}, which was significantly larger for the hypertensive subjects or for the healthy young subjects than for the elderly subjects (both healthy and hypertensive); *k*_{L}, which was significantly larger for hypertensive subjects than for the healthy young subjects; τ_{ach}, which was significantly smaller for hypertensive subjects than for the healthy young and elderly subjects; β, which was significantly larger for the elderly subjects (both healthy and hypertensive) than for the young subjects; and *M*_{P}, which is significantly different between healthy (young and elderly) and hypertensive subjects.

Finally, we have tested the capability of our model to predict the effects of age and hypertension. The effects of age were simulated using input pressure for a young subject combined with parameters predicted for a healthy elderly subject, and the effects of hypertension were simulated using input pressure for a healthy elderly subject combined with parameters for a hypertensive elderly subject. Only parameters that differed >50% between the two subjects (i.e., between the healthy young and the healthy elderly subject; and between the healthy elderly and the hypertensive elderly subject) were permuted. Results of these parameter mutations showed (see Fig. 9) that it is possible to use our model to predict effects of aging or hypertension. This figure shows sympathetic (baroreflex and vestibular mediated) and parasympathetic responses (*left* side), hysteresis curves (firing rate vs. blood pressure, *center*), and heart rate (*right* side). In all figures, results with permuted parameters are depicted with gray dashed lines and computations with original pressure inputs [i.e., for the healthy (Fig. 9*A*) and hypertensive (*B*) elderly subjects] are depicted with gray solid lines. The figure shows that simulations with permuted parameters (gray dashed lines) closely matched results obtained with original parameters (gray solid lines). The main difference between results with permuted parameters and those from original simulations are the firing rate/blood pressure responses. Because the blood pressure for the healthy young (Fig. 9*A*) and the healthy elderly (*B*) subject was used as inputs, the baroreflex firing rate curves were shifted to the left. However, besides this shift, they had the same shape as originally computed. The limitation of this parameter permutation study is that, if a young person is made elderly or if healthy elderly subject is made hypertensive, then aging and hypertension would also affect the input pressure. But our current model does not account for the feedback that heart rate and neural response have on blood pressure.

## DISCUSSION, SUMMARY, AND SIGNIFICANCE

In this study, we developed a mathematical model that can accurately predict heart rate dynamics during postural change from sitting to standing (see Table 2). The model was successfully validated against 60 data sets separated into three groups: healthy young people, healthy elderly people, and hypertensive people. Arterial baroreflexes play a key role in maintaining arterial pressure by regulating heart rate in the upright position. Heart rate dynamics were modeled in response to baroreceptor firing rate, sympathetic and parasympathetic responses, vestibulo-sympathetic reflex, and concentrations of norepinephrine and acetylcholine, triggered by blood pressure decline upon standing up. To predict heart rate dynamics, we first estimated afferent baroreflex firing rate. One very important result observed by graphing afferent baroreflex firing rate vs. mean blood pressure is the difference of the hysteresis curves. The area spanned by these curves is significantly different among all three groups of subjects, with young people having the largest area, hypertensive elderly people having intermediate areas, and healthy elderly people having very small areas. It seems natural to conclude that a large area indicates high baroreflex sensitivity and also that sensitivities to increases or decreases of blood pressure may differ. Therefore, these area measures may provide a better assessment of baroreflex sensitivity than average slopes of the blood pressure/firing rate relation because they take hysteresis into account. Our model also shows that, in some hypertensive subjects, the hysteresis curve is not closed, indicating that baroreflex modulation does not return to baseline during the first minute of standing. Another conclusion is that the dynamics of the afferent firing rate change decreases with age and that the resting firing rate decreases significantly for the hypertensive subjects. Again, these observations indicate reduced baroreflex sensitivity with aging and even further reduced baroreflex function for hypertensive subjects. These findings are important, as they may reflect combined effects of aging and hypertension on baroreceptor firing.

Most previous models of heart rate dynamics focus on modeling the neural response. Our model has shown that, if we do not explicitly include both the sympathetic time delay and account for vestibulo-sympathetic responses, we could not predict dynamics of heart rate responses to upright posture. In other words, heart rate responses to postural changes cannot be explained by baroreflex regulation alone. To our knowledge, this has not been shown in previous studies.

The fact that we explicitly account for latency of sympathetic vasoconstriction allowed us to observe the impact of aging on this delay. While the latency was not significantly changed between healthy and hypertensive elderly people, we observed a large difference between the young subjects and the elderly subjects. This increase in delay could account for inhibition of the sympathetic response with aging (46).

The vestibular and baroreflex-mediated sympathetic reflexes, which give rise to bimodal distribution, as well as the parasympathetic baroreflex response are separated explicitly in time, as shown in Fig. 5. The maximum or minimum responses are given in Table 4. For all three responses (vestibulo-sympathetic response, parasympathetic response, sympathetic response), the dynamics are diminished with age (no difference is detected between healthy and hypertensive elderly). However, comparing baroreflex sympathetic and vestibulo-sympathetic responses within each group shows that, for healthy young subjects (*P* = 0.487), the two responses cannot be distinguished, whereas for the elderly subjects, the baroreflex sympathetic response is diminished to a larger degree than the vestibulo-sympathetic response (*P* = 0.05) for both healthy and hypertensive elderly subjects. This latter observation is in agreement with the study by Ray and Mohanan (37).

Another conclusion we can draw from our model is that the parasympathetic inhibition of the sympathetic response is increased with age. This follows from the observation that the baseline sympathetic response is increased (although this is not statistically significant) combined with an increase of the dampening factor β (this is statistically significant, see Table 2), whereas the parasympathetic response is diminished.

The predictive simulations shown in Fig. 9 could be used to speculate that if responses of certain drugs are known to affect given combination of parameters then this model could potentially be used to simulate the response to various treatments.

#### Limitations.

It should be noted that only a few select parameters are significantly different between the three groups of people, and it remains to be studied why it is only these parameters that are affected and not all of them. One thing to note is the large standard deviation that is observed within each group. This is due to the fact that we see large variation in both blood pressure and heart rate measurements within each group of subjects. This study did not factor in the effect of gender, which is known to significantly affect heart rate dynamics; however, we did not have information about gender for the analyzed data sets. Furthermore, the pulmonary bed and cardiopulmonary interactions were not included in the model, and, therefore, the contributions from the cardiopulmonary receptors unloading to the initial and steady heart rate responses were not accounted for.

## GRANTS

This work was supported by a US Austria-Denmark cooperative research grant entitled “Modeling and Control of the Cardiovascular-Respiratory System,” Grant 0437037 from the National Science Foundation. Work performed at Beth Israel Deaconess Medical Center General Clinical Research Center was supported by the National Institutes of Health (NIH) under Grants M01-RR-01302, R01-NS-045745–01A2, and P60-AG-08812. Data collection and analysis was supported by a Joseph Paresky Men's Associates grant from Hebrew SeniorLife, a Research Nursing Home Grant AG-004390, and an Alzheimers Disease Research Center Grant AG-05134 from the National Institute on Aging. H. T. Tran was supported in part by the NIH under grant RO1- GM-067299–03. Part of this research was carried out in the REU summer program, which was held at North Carolina State University, May 29-August 4, 2005. Support for the REU program was provided by NSA (H98230–05-1–0284).

## Acknowledgments

Participants in the Research Experiences for Undergraduates (REU) program who contributed to this study are as follows: Robert Benim (Department of Mathematics, The University of Portland), Sarah Joyner (Department of Mathematics and Computer Science, Meredith College), Eamonn Tweedy (Department of Mathematics, North Carolina State University), and Chris Vogl (Department of Mathematical Sciences, Illinois Wesleyan University). Finally, Anthony Dixon (Department of Mathematics, North Carolina State University) has contributed to various aspects of this research including the *P* value calculations.

## Footnotes

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

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

- Copyright © 2006 the American Physiological Society