## Abstract

The most important goal of this study is to enhance our understanding of the crucial functional relationships that determine the behavior of the systemic circulation and its underlying physiological regulatory mechanisms with minimal modeling. To the present day, much has been said about the indirect hydraulic effects of right atrial pressure (P_{RA}) via cardiac output (CO) on arterial pressure (P_{a}) through the heart and pulmonary circulation or the direct regulatory effects of P_{RA} on P_{a} through the cardiopulmonary baroreflex; however, very little attention has been given to the hydraulic influence that P_{RA} exerts directly through the systemic circulation. The experimental data reported by Guyton et al. in 1957 demonstrated that steady-state P_{RA} and the rate at which blood passes through the systemic circulation are locked in a functional relationship independent of any consequence of altered P_{RA} on cardiac function. With this in mind, we emphasize the analytic algebraic analysis of the systemic circulation composed of arteries, veins, and its underlying physiological regulatory mechanisms of baroreflex and autoregulatory modulation of total peripheral resistance (TPR), where the behavior of the system can be analytically synthesized from an understanding of its minimal elements. As a result of this analysis, we present a novel mathematical method to determine short-term TPR fluctuations, which accounts for the entirety of observed P_{a} fluctuations, and propose a new cardiovascular system identification method to delineate the actual actions of the physiological mechanisms responsible for the dynamic couplings between CO, P_{a}, P_{RA}, and TPR in an individual subject.

- autonomic regulation
- arterial and cardiopulmonary baroreceptors
- local vascular autoregulation
- system identification
- cardiovascular regulatory physiology

the experimental data reported by Guyton et al. (13) demonstrated that steady-state right atrial pressure (P_{RA}) and the rate at which blood passes through the systemic circulation are locked in a functional relationship independent of any consequence of altered P_{RA} on cardiac function. This relationship has been reproduced in independent laboratories and with different methodologies (11, 15, 20). Accordingly, the lowest pressure point downstream the systemic circulation, P_{RA}, can hydraulically affect the highest pressure point upstream, thus arterial pressure (P_{a}), through two distinctly different pathways subsequently referred to as “direct” and “indirect” pathways. The indirect pathway via cardiac output (CO) is composed of the heart and pulmonary circulation and contains a time delay. The direct pathway is instantaneous and independent of the heart and pulmonary circulation, as demonstrated by Guyton et al. in 1957. As a result, the direct pathway depends only on the viscoelastic properties of the systemic circulation composed of the arterial and venous vasculature and its underlying physiological regulatory mechanisms.

The reflex regulation of the circulation by arterial baroreflexes originating in the carotid sinus and aortic arch regions has interested many physiologists for more than a century (18, 23, 35). A fall in pressure at these baroreceptors causes a reflex decrease in parasympathetic discharge and an increase in sympathetic discharge. These reflex changes lead to an increase in heart rate (HR), myocardial contractility, venoconstriction, and total peripheral resistance (TPR), all of which tend to maintain P_{a}. The cardiopulmonary baroreceptors, which lie mostly in the cardiac chambers, can also exert important effects on the circulation (1, 5, 24). An increase in pressure at these receptors causes reflex vasodilation, whereas a decrease causes reflex vasoconstriction. Central control of TPR by the arterial and cardiopulmonary baroreflexes, however, can be readily opposed by distally initiated control mechanisms. Metabolic modulation within active skeletal muscle attenuates sympathetically mediated vasoconstriction, so that oxygen delivery to the active muscles is preserved and local metabolic demands are met (14, 34, 36, 38, 41). Furthermore, blood flow through almost any organ changes very little if the pressure perfusing the arteries of the organ is varied (7, 17, 22). This locally controlled mechanism is termed vascular autoregulation and can contribute significantly to short-term control of TPR.

In this study, we emphasize the analytic algebraic analysis of the systemic circulation composed of arteries, veins, and its underlying physiological control mechanisms of baroreflex and autoregulatory modulation of TPR with minimal modeling. We use Frank's (10) two-element windkessel model of the arterial circulation, which successfully explains diastole as the discharging of a blood volume integrator in which P_{a} falls exponentially and uniformly. It is the most popular model of the arterial system for academic purposes (4, 19) and has been widely utilized to estimate arterial compliance (*C*_{a}) and stroke volume. The two-element windkessel has formed the basis of different methods to estimate *C*_{a}, such as the decay time method (10), the area method (30), the two-area method (37), and the pulse pressure method (40). The most common criticism to Frank's windkessel is that the high-frequency components of the aortic pressure wave are poorly represented. This criticism has led to many windkessel models and many concerns about their applications (6, 43, 44). We believe that more complexity is required only if the use of a simple model does not suffice to describe the properties of the system essential for determining its functioning. Stergiopulos et al. (40) showed that Frank's windkessel's unrealistic prediction of aortic pressure waveforms was due to its poor medium- to high-frequency representation of the aortic input impedance. They concluded, however, that the two-element windkessel describes the low-frequency characteristics of the arterial system correctly and can predict pulse pressure well for the right choice of *C*_{a}. Accordingly, Frank's simple windkessel model suffices for our analytic analysis of the arterial circulation where only the low-frequency characteristics are of interest to us. To include the venous circulation as part of our analysis, we modify Frank's windkessel model as to incorporate Rowell's (33) view of the properties of the venous circulation.

As a result of the aforementioned analysis, we developed a novel mathematical method able to not only track down steady-state changes in TPR very effectively but also short-term TPR fluctuations caused by baroreflex and autoregulatory modulation of TPR. Consequently, short-term P_{a} fluctuations can be successfully explained almost in their entirety when the contributions of ΔTPR to Pa are taken into account. This method, however, does not provide any information regarding the source of TPR fluctuations. With this in mind, we complemented this method with a separate set of quantitative tools, which can be applied to delineate the actual actions of the physiological mechanisms responsible for the observed short-term TPR fluctuations in an individual subject, without altering the underlying regulatory mechanisms. A previous publication by the same group (26), which stresses HR modulation by the arterial baroreflexes and respiration, presents a cardiovascular system identification method intended for the analysis of HR derived from the surface electrocardiogram, noninvasively measured P_{a}, and instantaneous lung volume signals. For that purpose, Perrot and Cohen (29) proposed a model order reduction technique optimized for that particular problem. In this study, we deal with a different problem. We stress TPR modulation by the baroreflexes and autoregulation and use standard model order selection techniques (31). As a result, we propose a novel cardiovascular system identification method intended for the analysis of the independent dynamic closed-loop contributions of CO and P_{RA} to short-term P_{a} fluctuations as well as the independent dynamic closed-loop contributions of P_{a} and P_{RA} to short-term TPR fluctuations, which serves as a natural complement to the cardiovascular system identification method previously proposed by the same group (26), where HR fluctuations play a major role and the regulatory mechanisms responsible for short-term TPR fluctuations could not be accounted for. Furthermore, in a recent publication by the same group, Mukkamala and Cohen (25) estimated short-term TPR fluctuations simply as TPR=P_{a}/CO without considering the effects of arterial elasticity and were unable to directly determine the dynamic transfer relation P_{a}→TPR. By including an analysis on the effects of arterial elasticity, we are now able to successfully formulate a direct system identification approach capable of estimating this transfer relation.

### Glossary

- ABR
- Arterial baroreflex
*C*_{a}- Arterial capacitance
*C*_{RA}- Right atrial capacitance
*C*_{v}- Venous capacitance
- CBR
- Cardiopulmonary baroreflex
- CO
- Cardiac output
*f*_{c}- Cutoff frequency
*f*_{s}- Sampling frequency
*G*{*H*_{CO→Pa}}- Static gain of
*H*_{CO→Pa} *G*{*H*_{Pa}_{→TPR}}- Static gain of
*H*_{Pa}_{→TPR} *G*{*H*_{PRA}_{→Pa}}- Static gain of
*H*_{PRA}_{→Pa} *G*{*H*_{PRA}_{→TPR}}- Static gain of
*H*_{PRA}_{→TPR} *G*{*K*_{CO→Pa}}- Static gain of
*K*_{CO→Pa} *G*{*K*_{CO→TPR}}- Static gain of
*K*_{CO→TPR} *G*{*K*_{Pa}_{→TPR}}- Static gain of
*K*_{Pa}_{→TPR} *G*{*K*_{PRA}_{→Pa}}- Static gain of
*K*_{PRA}_{→Pa} *G*{*K*_{PRA}_{→TPR}}- Static gain of
*K*_{PRA}_{→TPR} *G*{*K*_{TPR→Pa}}- Static gain of
*K*_{TPR→Pa} *H*_{CO→Pa}- Closed-loop transfer function of CO→P
_{a} *H*_{Pa}_{→TPR}- Closed-loop transfer function of P
_{a}→TPR *H*_{PRA}_{→Pa}- Closed-loop transfer function of P
_{RA}→P_{a} *H*_{PRA}_{→TPR}- Closed-loop transfer function of P
_{RA}→TPR - HR
- Heart rate
*k*- Time sample
*K*_{CO→Pa}- Open-loop transfer function of CO→P
_{a} *K*_{CO→TPR}- Transfer function of vascular autoregulation CO→TPR
*K*_{Pa}_{→TPR}- Transfer function of arterial baroreflex P
_{a}→TPR *K*_{PRA}_{→Pa}- Open-loop transfer function of P
_{RA}→P_{a} *K*_{PRA}_{→TPR}- Transfer function of cardiopulmonary baroreflex P
_{RA}→TPR *K*_{TPR→Pa}- Open-loop transfer function of TPR→P
_{a} - P
_{a} - Arterial pressure
- P
_{c} - Capillary pressure
- P
_{RA} - Right atrial pressure
*R*_{a}- Arterial resistance
*R*_{v}- Venous resistance
*s*- Analog complex frequency
*t*- Continuous time
*T*_{ar}- Time delay of VAR
*T*_{br}- Time delay of ABR and CBR
- TPR
- Total peripheral resistance
- VAR
- Vascular autoregulation
*z*- Digital complex frequency
- τ
_{1} - Characteristic
*time constant 1*of*K*_{CO→Pa},*K*_{PRA}_{→Pa}, and*K*_{TPR→Pa} - τ
_{2} - Characteristic
*time constant 2*of*K*_{CO→Pa},*K*_{PRA}_{→Pa}, and*K*_{TPR→Pa} - τ
_{br} - Time constant of ABR, CBR, and VAR
- τ
_{eff} - Effective time constant of systemic circulation

## MATHEMATICAL ANALYSIS OF THE SYSTEMIC CIRCULATION

The systemic circulation composed of arterial and venous vasculature can be represented by a windkessel model as illustrated in Fig. 1, where we characterize the blood storage capacity of the large and small arterial vessels and terminal arterial branches (arterioles) as a capacitance element, *C*_{a}. Thepressure gradient from the aorta across the large and small arterial vessels and the arterioles to the junction between arterial and venous circulation is caused by the arterial resistance, *R*_{a}. The capillaries as well as the terminal venous branches (venules) and small veins responsible for venous blood storage constitute the venous capacitance, *C*_{v}, whereas the large conduit veins responsible for the volume of blood available for ventricular filling are characterized as a capacitance element, *C*_{RA}. Finally, we include an additional resistance to venous return, *R*_{v}, responsible for the pressure gradient from the capillary pressure, P_{c}, across the venules and small veins to the large conduit veins. It is this pressure gradient that determines the passive effects of blood flow on venous volume and the volume of blood available for ventricular filling. As a result, our modified windkessel model describes *1*) the drainage of blood into the periphery (P_{RA}) through TPR (TPR=*R*_{a} + *R*_{v}), *2*) the elasticity of the arteries (*C*_{a}) as responsible for the continued egress of blood from the peripheral terminal arteries during ventricular diastole, and *3*) the pressure gradient across the venous system (P_{c} − P_{RA}) as responsible for the volume of blood available for ventricular filling. Accordingly, the aortic pressure waveform P_{a}(*t*) can be represented as the viscous blood flow through *R*_{a} times *R*_{a} plus the capillary pressure waveform P_{c}(*t*) yielding (1) where CO(*t*) represents the aortic flow. Similarly, P_{c}(*t*) can be represented as the blood flow through *R*_{v} times *R*_{v} plus the P_{RA} waveform P_{RA}(*t*) yielding (2) At first, assume for simplicity that *R*_{a} and *R*_{v} are constant; therefore, Laplace transformation of both sides of *Eq. 1* from time-domain to its frequency-domain representation yields (3) where *s* denotes the complex frequency. The Laplace transform is an operator, which replaces time-domain operations such as differentiation with multiplication by *s* and time delays (*T*) with multiplication by *e*^{−Ts} (28); as a result, an algebraic equation arises rather than a differential equation. Thus solving for P_{a}(*s*) in *Eq. 3* results in the algebraic expression for P_{a}(*s*) as (4) Likewise, Laplace transformation of *Eq. 2* from time-domain into its frequency-domain representation yields (5) and solving for P_{c}(*s*) in *Eq. 5* results in (6) Substitution of *Eq. 6* into *Eq. 4* and subsequently solving for P_{a}(*s*) finally leads to (7) where *K*_{CO→Pa}(*s*) represents the direct hydraulic effects of CO→P_{a} and *K*_{PRA}_{→Pa}(*s*) represents the direct hydraulic effects of P_{RA}→P_{a}. Figure 2 illustrates the block diagram representation of *Eq. 7* with the impulse-response function representations of the linear and time-invariant (LTI) systems *K*_{CO→Pa}(*s*) and *K*_{PRA}_{→Pa}(*s*) when, e.g., *R*_{a}=1.25 mmHg·s·ml^{−1}, *C*_{a}=1 ml/mmHg, *R*_{v}=*R*_{a}/5, and *C*_{v}=16·*C*_{a} (4, 12).

### Baroreflex Control of TPR

If *R*_{a} and *R*_{v} are not constant due to baroreflex modulation of TPR (TPR=*R*_{a} + *R*_{v}), one has to appropriately expand the block diagram representation shown in Fig. 2 as to account for the additional effects of TPR changes (ΔTPR) on P_{a} as illustrated in Fig. 3*A*, where the transfer relation *K*_{Pa}_{→TPR}(*s*) represents the delayed autonomically mediated effects of P_{a}→TPR (arterial baroreflex) and the transfer relation *K*_{PRA}_{→TPR}(*s*) represents the delayed autonomically mediated effects of P_{RA}→TPR (cardiopulmonary baroreflex). Accordingly, arterial and cardiopulmonary baroreflex modulation of TPR can be represented as (8) The transfer relations *K*_{CO→Pa}(*s*), *K*_{PRA→Pa}(*s*), and *K*_{TPR→Pa}(*s*) represent the immediate direct hydraulic effects of CO, P_{RA}, and TPR on P_{a}, respectively. Thus P_{a}(*s*) can be represented as (9) As a result, it is possible to determine the closed-loop effects of CO and P_{RA} on P_{a} representing both hydraulic as well as autonomic interactions simply by substitution of *Eq. 8* into *Eq. 9* and subsequently solving for P_{a}(*s*) as (10) where *H*_{CO→Pa}(*s*) and *H*_{PRA}_{→Pa}(*s*) represent the independent dynamic closed-loop effects of CO→P_{a} and P_{RA}→P_{a} as illustrated in Fig. 3*B*, which is equivalent to Fig. 3*A*. The solid curves in Fig. 3*B* illustrate the dynamic properties, composed of hydraulic and autonomic interactions, of the closed-loop transfer relations *H*_{CO→Pa}(*s*) and *H*_{PRA}_{→Pa}(*s*) depicted in their time-domain form of impulse-response function representations. In contrast, the dashed curves illustrate the dynamic properties of *K*_{CO→Pa}(*s*) and *K*_{PRA}_{→Pa}(*s*) in the absence of any autonomic interactions.

The autonomically mediated transfer relations *K*_{Pa}_{→TPR}(*s*) and *K*_{PRA}_{→TPR}(*s*) can be characterized by the first order LTI systems of the form (11) and (12) where τ_{br} denotes the time constant of the sympathetic nervous system acting through the α-adrenergic pathway (32); the time delay *T*_{br} corresponds to the time it takes for the sympathetic nervous system to receive the information, process it, and start a reaction to it (8, 42); and the coefficients *G*{*K*_{Pa}_{→TPR}} and *G*{*K*_{PRA}_{→TPR}} represent the static gains of the arterial and cardiopulmonary baroreflexes, respectively. The hydraulically mediated transfer relations *K*_{CO→Pa}(*s*), *K*_{PRA}_{\#8594;Pa}(*s*), and *K*_{TPR→Pa}(*s*) can be characterized by the following second-order LTI systems (13) and (14) and (15) respectively, where τ_{1} and τ_{2} represent the time constants characterizing the combined viscoelastic properties of the arterial and venous circulation. As a result, it is possible to show by substitution of *Eqs. 11*–*15* into *Eq. 10* that the closed-loop transfer relations *H*_{CO→Pa}(*s*) and *H*_{PRA}→P_{a}(*s*) must be characterized by the third-order closed-loop systems of the form (16) and (17) respectively, where the closed-loop characteristic frequencies are determined by the dynamic properties of the feedback system that characterizes the arterial baroreflex, *K*_{Pa}_{→TPR}(*s*). For a detailed discussion on this issue, please see the appendix.

#### Steady-state analysis.

Steady-state analysis of the modified windkessel model in Fig. 1 yields (18) where the overbars denote mean values. Therefore, the steady-state changes in P_{a} can be described by (19) where Δ denotes steady-state changes in that signal around its mean value. Normalization of *Eq. 19* by and multiplication by 100% to describe percent changes finally results in (20) where *G*{*K*_{CO→Pa}}, *G*{*K*_{PRA}_{→Pa}}, and *G*{*K*_{TPR→Pa}} represent the static gains of the open-loop transfer relations *K*_{CO→Pa}(*s*), *K*_{PRA}_{→Pa}(*s*), and *K*_{TPR→Pa}(*s*), respectively, and the Δ denotes percent changes in that signal around its mean value. ΔP_{RA}, however, denotes changes in mmHg around its mean value, because, unlike the higher pressure in P_{a}, a 1-mmHg error in absolute P_{RA} resulting from deviations between the precise location of the right atrium and the pressure transducer such as in the companion paper (2) can introduce a substantial degree of error if given in percent. Fluctuations in P_{RA} given in mmHg remain unaffected and offer a solution to this problem. A simple mathematical expression can be finally attained for the steady-state relation between autonomically and hydraulically mediated interactions by evaluation of *H*_{CO→Pa}(*s*) in *Eq. 16* at *s*=0 as (21) where *G*{*H*_{CO→Pa}} represents the static gain of *H*_{CO→Pa}(*s*). Thus it is possible to determine *G*{*K*_{Pa}_{→TPR}} by solving *Eq. 21* for *G*{*K*_{Pa}_{→TPR}} as (22) Likewise, evaluation of *H*_{PRA}_{→Pa}(*s*) in *Eq. 17* at *s*=0 results in (23) where *G*{*H*_{PRA}_{→Pa}} represents the static gain of *H*_{PRA}_{→Pa}(*s*). Accordingly, the static gain of *H*_{PRA}_{→Pa}(*s*) is said to depend on the separate static effects of two physically isolated sensory regions (arterial and cardiopulmonary baroreceptors) on a common effector region (TPR): the arterial and the cardiopulmonary baroreflexes. Consequently, it is possible to determine *G*{*K*_{PRA}_{→Pa}} by substitution of *G*{*K*_{Pa}_{→TPR}} from *Eq. 22* into *Eq. 23* and subsequently solving for *G*{*K*_{PRA}_{→Pa}} as (24) Essentially, *Eqs. 21* and *23* demonstrate that *G*{*H*_{CO→Pa}} is independent of *G*{*K*_{PRA}_{→Pa}} (cardiopulmonary baroreflex), whereas *G*{*H*_{PRA}_{→Pa}} depends on both *G*{*K*_{Pa}_{→TPR}} and *G*{*K*_{PRA}_{→TPR}} (arterial and cardiopulmonary baroreflexes).

#### Physical interpretation of the mathematical theory.

Consider the following hypothetical experiment where the right atrium is decoupled from the left ventricle and CO is maintained constant by externally controlled volume infusion into the left ventricle while the resulting venous return drains freely through the right atrium into an infinite blood reservoir. If the pressure downstream the systemic circulation (P_{RA}) changes due to, for example, occlusion of the right atrium, the pressure upstream the systemic circulation (P_{a}) has to adjust itself to overcome those pressure changes while “CO remains constant,” as illustrated by the upward arrows in Fig. 1. The dynamic nature of this adjustment depends on the viscoelastic properties of the system that couples P_{RA} and P_{a}, which is not only composed of the arterial vasculature (*R*_{a} and *C*_{a} with *R*_{a}*C*_{a}=τ_{a}) but also the venous vasculature (*R*_{v} and *C*_{v} with *R*_{v}*C*_{v}=τ_{v} ≠ τ_{a}). As a result, the dynamic open-loop relationship between P_{RA} and P_{a}, represented by *K*_{PRA}_{→Pa}(*s*), may be quantitatively described in terms of the windkessel elements *R*_{a}, *C*_{a}, *R*_{v}, and *C*_{v}, as illustrated algebraically in *Eq. 7* or, i.e., in terms of the “combined viscoelastic properties” of the arteries and veins. The open-loop characteristic frequencies of *K*_{PRA}_{→Pa}(*s*) may be determined from the roots of their denominator polynomials as *s*_{1,2}=1/τ_{1,2}, where both τ_{1} and τ_{2} depend on all windkessel parameters; thus τ_{1}=*f*(*R*_{a}, *C*_{a}, *R*_{v}, *C*_{v}) ≠ τ_{a} and τ_{2}=*f*(*R*_{a}, *C*_{a}, *R*_{v}, *C*_{v}) ≠ τ_{v}. Because TPR=*R*_{a} + *R*_{v} is not constant, but rather modulated by the centrally mediated baroreflexes *K*_{Pa}_{→TPR}(*s*) and *K*_{PRA}_{→TPR}(*s*), these reflex interactions must therefore, in addition to τ_{1} and τ_{2}, inherently influence the dynamic closed-loop transfer relationship between P_{RA} and P_{a}, represented by *H*_{PRA}_{→Pa}(*s*), as demonstrated in *Eq. 10* [note the presence of *K*_{Pa}_{→TPR}(*s*) in the denominator and *K*_{PRA}_{→TPR}(*s*) in the numerator] and as illustrated graphically in Fig. 3*B*, where the gray area in *H*_{PRA}_{→Pa}(*s*) denotes the change in P_{a} due to arterial and cardiopulmonary baroreflex regulation of TPR given an impulse change in P_{RA}. As a result, *H*_{PRA}_{→Pa}(*s*) can be said to also depend, in addition to the combined viscoelastic properties of arteries and veins, on the separate dynamic effects of two physically isolated sensory regions (arterial and cardiopulmonary) on a common effector region (TPR); hence, the arterial and the cardiopulmonary baroreflexes. Any hydraulic influence of P_{RA} on P_{a} via the heart and pulmonary circulation is implicit in CO and thus already incorporated into the dynamic closed-loop transfer relationship between CO and P_{a}, represented by *H*_{CO→Pa}(*s*). Consequently, the closed-loop transfer relations *H*_{CO→Pa}(*s*) and *H*_{PRA}_{→Pa}(*s*) are meant to represent the independent dynamic closed-loop contributions of CO and P_{RA} on Pa fluctuations. The shaded area in *H*_{CO→Pa}(*s*) in Fig. 3*B* denotes the change in P_{a} due to arterial baroreflex regulation of TPR given an impulse change in CO.

Figure 4*A* displays *H*_{CO→Pa}(*s*) and *H*_{PRA}_{→Pa}(*s*) in their time-domain form of step-response function representations, whereas Fig. 4*B* illustrates the dynamic properties of the arterial and cardiopulmonary baroreflexes in their time-domain form of step-response function representations. Δ denotes the percent fluctuations in that signal around its mean value denoted by the overbar with the exception of P_{RA}, where Δ denotes fluctuations in mmHg. In essence, Fig. 4*A* illustrates how *H*_{CO→Pa}(*s*) depends solely on the arterial baroreflex *K*_{Pa}_{→TPR}(*s*), whereas *H*_{PRA}_{→Pa}(*s*) depends on both baroreflexes, *K*_{Pa}_{→TPR}(*s*) and *K*_{PRA}_{→TPR}(*s*). To further elaborate on the previous statement, three cases will be individually addressed in this paragraph. First, the case where both baroreflexes are absent shall be considered, as illustrated by the dashed curves in Fig. 4*A*. Therefore, if CO were to increase by a step change of 1% while P_{RA} remained unchanged, then P_{a} would rapidly increase *G*{*K*_{CO→Pa}}% in <10 s, which is the approximate time necessary to reach a steady state in P_{a} given a step change in CO, i.e., about 10 times the dominant time constant of the open-loop system *K*_{CO→Pa}(*s*) mainly determined by the viscoelastic properties of the arterial circulation. Likewise, if P_{RA} were to increase by a step change of 1 mmHg while CO remained constant, P_{a} would increase *G*{*K*_{PRA}_{→Pa}}% in <20 s, which is the approximate time necessary to reach a steady state in P_{a} given a step change in P_{RA}, i.e., about 10 times the dominant time constant of the open-loop system *K*_{PRA}_{→Pa}(*s*) determined by the combined viscoelastic properties of the arterial and venous circulation. The case where both baroreflexes are present shall be considered next, as illustrated by the solid curves in Fig. 4*A*. Thus if CO were to increase by a step change of 1% while P_{RA} remained unchanged, P_{a}'s rapid increase would be truncated and reversed by the nonimmediate negative feedback effects of the arterial baroreflex on TPR (see time delay in Fig. 4*B*). In this example, the arterial baroreflex kicks in at *t*=2 s forcing P_{a} to settle down in about 30 s into a steady-state value smaller than *G*{*K*_{CO→Pa}} denoted as *G*{*H*_{CO→Pa}}. The static gain *G*{*H*_{CO→Pa}} is determined by *Eq. 21*, which states the mathematical relation between *G*{*H*_{CO→Pa}} and the static gain of the arterial baroreflex *G*{*K*_{Pa}_{→TPR}}. Likewise, if P_{RA} were to increase by a step change of 1 mmHg while CO remained constant, P_{a}'s increase would be truncated and reversed by the nonimmediate negative feedback regulation of the arterial baroreflex in addition to the nonimmediate negative cardiopulmonary baroreflex regulation of TPR (see time delay in Fig. 4*B*). In this example, both baroreflexes kick in at *t*=2 s and force P_{a} to settle down in about 30 s into a steady-state value remarkably lower than *G*{*K*_{PRA}_{→Pa}} denoted as *G*{*H*_{PRA}_{→Pa}}. The static gain *G*{*H*_{PRA}_{→Pa}} is determined by *Eq. 23*, which states the mathematical relation between *G*{*H*_{PRA}_{→Pa}} and the static gains of the arterial and cardiopulmonary baroreflexes *G*{*K*_{Pa}_{→TPR}} and *G*{*K*_{PRA}_{→TPR}}. Finally, the special case where the arterial baroreflex is present but the cardiopulmonary baroreflex is absent shall be considered, as illustrated by the light gray dashed-dotted curve in Fig. 4*A* [present only in *H*_{PRA}_{→Pa}(*s*) because *H*_{CO→Pa}(*s*) is independent of the cardiopulmonary baroreflex]. Accordingly, *H*_{PRA}_{→Pa}(*s*) resembles a scaled version of *H*_{CO→Pa}(*s*). With this particular case, it becomes apparent that only the presence of the cardiopulmonary baroreflex can result in *G*{*H*_{PRA}_{→Pa}} < 0. The mathematical relation given in *Eq. 23* states that *G*{*H*_{PRA}_{→Pa}} depends on the separate static effects of two physically isolated sensory regions (arterial and cardiopulmonary baroreceptors) on a common effector region (TPR). Specifically, *Eq. 23* limits the influence of arterial baroreflex modulation of TPR to only yield positive values in *G*{*H*_{PRA}_{→Pa}} smaller than *G*{*K*_{PRA}_{→Pa}}, because *G*{*K*_{Pa}_{→TPR}} is always negative. As a result, a negative value in *G*{*H*_{PRA}_{→Pa}} can only be explained by cardiopulmonary baroreflex modulation of TPR because *G*{*K*_{PRA}_{→TPR}} is always negative and expressed in the numerator.

### Effects of Local Vascular Autoregulation

Vascular autoregulation is characterized by local vasoconstriction after an increase of the vessel's transmural pressure and local vasodilation after a decrease. As a result, autoregulation implies a constant blood flow in the face of altered perfusion pressure, and thus, if present, autoregulation contributes to short-term control of TPR. Local vascular autoregulation is probably determined by both mechanical (P) and metabolic (Q) changes. Because the local pressures P_{i} (which are not in the model) do not add up to a systemic pressure P, and while all local Q_{i} (which are not in the model either) do add up to a systemic flow Q made up of the sum of all local Q_{i} as Q=∑ Q_{i}, one can think of autoregulation as being determined by both local metabolic (Q_{i}) and local mechanical (P_{i} which results in Q_{i} ∼ P_{i} anyway) changes, where CO=Q can be viewed as a plausible measure of autoregulation. Consequently, the presence of autoregulatory modulation of TPR in addition to autonomic baroreflex regulation makes it necessary to expand the block diagram representation of the systemic circulation depicted in Fig. 3*A* as to account for the short-term dynamic autoregulatory effects of CO on TPR as illustrated in Fig. 5*A*. Accordingly, *Eq. 10* becomes (25) where *K*_{CO→TPR}(*s*) represents the dynamic autoregulatory effects of CO on TPR. Consequently, *Eq. 25* describes how *H*_{CO→Pa}(*s*) is affected by *K*_{CO→TPR}(*s*) while *H*_{PRA}_{→Pa}(*s*) remains independent of *K*_{CO→TPR}(*s*), which can be characterized, e.g., by the first-order LTI system of the form (26) where *G*{*K*_{CO→TPR}} denotes the static gain of *K*_{CO→TPR}(*s*) and *T*_{ar} corresponds to the time it takes for the local vasculature to sense a change in blood flow and subsequently react to it (22, 38). The time constant of this system can be assumed to be very similar to τ_{br}, the time constant of the baroreflexes (22, 32). Likewise, from the block diagram shown in Fig. 5*A*, TPR can then be determined as (27) Thus solving *Eq. 27* for CO(*s*) and subsequently substituting into *Eq. 25* yields (28) after appropriately solving for TPR(*s*). Consequently, *Eq. 28* describes how the baroreflexes *K*_{Pa}_{→TPR}(*s*) and *K*_{PRA}_{→TPR}(*s*) are influenced by autoregulation *K*_{CO→TPR}(*s*). As a result, the dynamic closed-loop transfer relation *H*_{Pa}_{→TPR}(*s*) describes the interaction between the centrally mediated arterial baroreflex and local vascular autoregulation, whereas the dynamic closed-loop transfer relation *H*_{PRA}_{→TPR}(*s*) describes the interaction between the centrally mediated cardiopulmonary baroreflex and local vascular autoregulation. As a result, the independent dynamic closed-loop contributions of P_{a} and P_{RA} on TPR can differ considerably from the simple first-order open-loop transfer relations *K*_{Pa}_{→TPR}(*s*) and *K*_{PRA}_{→TPR}(*s*) depending on the degree of local vascular autoregulation present.

#### Steady-state analysis.

The static properties of the closed-loop transfer relation *H*_{CO→Pa}(*s*) in *Eq. 25* represented as (29) together with the static properties of *H*_{Pa}_{→TPR}(*s*) and *H*_{PRA}_{→TPR}(*s*) in *Eq. 28*, respectively, represented as (30) and (31) and the static gain *G*{*H*_{PRA}_{→Pa}} from *Eq. 23* provide simple mathematical expressions for the steady-state relations between regulatory and hydraulically mediated interactions. For instance, it can be easily verified that in the presence of autoregulation, as graphically illustrated in Fig. 6*B*, the amplitudes of G{*H*_{Pa}_{→TPR}} and *G*{*H*_{PRA}_{→TPR}} must always be smaller than the amplitudes of the individual baroreflexes (*G*{*K*_{Pa}_{→TPR}} and *G*{*K*_{PRA}_{→TPR}}) due to the fact that *G*{*K*_{CO→TPR}} and *G*{*H*_{CO→Pa}} must always have positive numbers and the baroreflexes negative values. Finally, solving *Eq. 30* for *G*{*K*_{Pa}_{→TPR}} and consequently substituting the resulting *G*{*K*_{Pa}_{→TPR}} into *Eq. 29* yields (32) after appropriately solving *Eq. 29* for *G*{*H*_{Pa}_{→TPR}}. Likewise, the solution of *Eq. 31* for *G*{*K*_{PRA}_{→TPR}} and consequent substitution of the resulting *G*{*K*_{PRA}_{→TPR}} and *G*{*K*_{Pa}_{→TPR}} into *Eq. 23* results in (33) after appropriately solving *Eq. 23* for *G*{*H*_{PRA}_{→TPR}}.

#### Physical interpretation of the mathematical theory.

Figure 5*B* graphically illustrates the dynamic closed-loop transfer relations *H*_{CO→Pa}(*s*) and *H*_{PRA}_{→Pa}(*s*) in their time-domain form of impulse-response function representations in the presence of substantial local vascular autoregulation in addition to autonomic baroreflex modulation of TPR. The dashed-dotted curves represent the case where TPR is constant; thus no TPR regulation takes place. The dashed curve represents the case with substantial autoregulatory modulation of TPR but no baroreflexes [present only in *H*_{CO→Pa}(*s*) because *H*_{PRA}_{→Pa}(*s*) is independent of autoregulation]. The solid curves represent the case with baroreflex and autoregulatory modulation of TPR. The dynamic properties of *H*_{CO→Pa}(*s*) reveal in all three cases an immediate initial increase in P_{a} given a positive impulse change in CO demonstrating the direct positive hydraulic effects of CO on P_{a}. In this example, local vascular autoregulation kicks in at ∼1 s to slow down P_{a}'s natural time decay in *H*_{CO→Pa}(*s*) as a result of the longer time constant of autoregulation. Up to that point in time (*t*=1 s), P_{a}'s natural time decay is mainly determined by the viscoelastic properties of the arteries. At ∼2 s, the arterial baroreflex kicks in to speed up P_{a}'s natural time decay and further force P_{a} into negative territory as a result of its negative feedback regulation. Similarly, the dynamic properties of *H*_{PRA}_{→Pa}(*s*) reveal an initial increase in P_{a} given a positive impulse change in P_{RA} demonstrating the direct positive hydraulic effects of P_{RA} on P_{a}. P_{a}'s natural time decay is determined by the combined viscoelastic properties of arteries and veins up to the time when the arterial and the cardiopulmonary baroreflexes kick in to dramatically speed up P_{a}'s natural time decay and quickly lead P_{a} into negative values as a result of negative feedback regulation of the arterial baroreflex in addition to cardiopulmonary baroreflex regulation of TPR.

Figure 6 graphically illustrates with solid curves the dynamic closed-loop transfer relations *H*_{Pa}_{→TPR}(*s*) and *H*_{PRA}_{→TPR}(*s*) in the presence of substantial local vascular autoregulation in addition to autonomic baroreflex modulation of TPR. The dashed curves represent the case where autoregulation is practically absent; thus *H*_{Pa}_{→TPR}(*s*)=*K*_{Pa}_{→TPR}(*s*) and *H*_{PRA}_{→TPR}(*s*)=*K*_{PRA}_{→TPR}(*s*). Figure 6*A* depicts impulse-response function representations, whereas Fig. 6*B* displays step responses. The effects of autoregulation become especially apparent in *H*_{Pa}_{→TPR}(*s*), where the initial increase in TPR at 1 s clearly demonstrates the nonimmediate positive regulatory effects of P_{a} on TPR due to local vascular autoregulation. If the pressure perfusing the arteries of almost any organ is varied, blood flow through the organ changes very little. Therefore, an increase in P_{a} results in a likewise increase in TPR opposing the imminent increase in organ blood flow. In this example, local vascular autoregulation kicks in at ∼1 s and remains unchallenged until ∼2 s when the arterial baroreflex kicks in to oppose the increase in P_{a}, thereby reversing TPR and forcing it to remain negative until the end of its natural time decay. The apparent conflict between vascular autoregulation and baroreflexes results in *1*) an effective shortening of the characteristic time constant compared with τ_{br} as illustrated in Fig. 6*A*, and *2*) smaller static gains *G*{*H*_{Pa}_{→TPR}} and *G*{*H*_{PRA}_{→TPR}} compared with *G*{*K*_{Pa}_{→TPR}} and *G*{*K*_{PRA}_{→TPR}}, respectively, as illustrated in Fig. 6*B*. Furthermore, depending on the value of *G*{*H*_{Pa}_{→TPR}}, it is possible to clearly determine the dominance of one reflex over the other as follows: if *G*{*H*_{Pa}_{→TPR}} > 0, autoregulation dominates over the arterial baroreflex. If *G*{*H*_{Pa}_{→TPR}} < 0, the arterial baroreflex dominates over autoregulation. If *G*{*H*_{Pa}_{→TPR}} ≈ 0, the arterial baroreflex and autoregulation cancel each other out, as graphically illustrated in Fig. 6*B* by the solid curve representing *H*_{Pa}_{→TPR} with *G*{*K*_{CO→TPR}}=−*G*{*K*_{Pa}_{→TPR}}.

In conclusion, autoregulatory modulation of TPR slows down P_{a}'s natural time decay as a result of its longer time constant, whereas arterial baroreflex modulation of TPR speeds up P_{a}'s natural time decay as a result of its negative feedback regulation. Finally, the ongoing conflicting effects between the locally controlled vascular autoregulation and the centrally regulated autonomic baroreflexes manifest in the dynamic closed-loop transfer relations *H*_{Pa}_{→TPR}(*s*) and *H*_{PRA}_{→TPR}(*s*), which no longer represent sole control of TPR by the baroreflexes but also by local vascular autoregulation.

### Closed-Loop Computational Model of Heart and Circulation

To apply the methods for TPR determination and cardiovascular system identification presented next to the analysis of experimentally acquired animal data such as in the companion paper (2), it is necessary to first evaluate the proposed methods on a theoretical basis where the TPR fluctuations and the transfer relations of interest are known a priori. To carry out an evaluation with CO, P_{RA} and P_{a} signals assembled in closed loop, we implement a simple closed-loop computational model of the heart and circulation where the feedforward component characterizes cardiac function and the feedback component constitutes the peripheral vascular system properties represented by the lumped parameter model of the systemic circulation illustrated in Fig. 1 and its underlying regulatory mechanisms of baroreflex and autoregulatory modulation of TPR. Accordingly, the output variable, CO, is linked to its own input variables, P_{a} and P_{RA}, by the negative feedback interaction of the vascular subdivision. As a result, other things being equal, an increase in P_{a} causes decreased CO in the cardiac subdivision and a decrease in CO causes a decrease in P_{a} in the vascular subdivision. Similarly, other things being equal, an increase in P_{RA} causes increased CO in the cardiac subdivision and an increase in CO causes decreased P_{RA} in the vascular subdivision. In view of our minimal model approach, we implement a rather simple pacemaker-driven impulse model of cardiac ejection, where the implemented Starling mechanism is not offered as a new model for heart dynamics but simply to close the loop between the heart and circulation and thereby provide computationally generated data sets of CO, P_{RA}, P_{a}, and TPR signals assembled in closed loop, which could be subsequently used for system identification purposes and as the gold standard for TPR determination.

## DETERMINATION OF TPR

The differential equation represented in *Eq. 1* describes *1*) the drainage of blood into the periphery (P_{RA}) through the outflow resistance (TPR) and *2*) the elasticity of the arteries (*C*_{a}) as responsible for the continued egress of blood from the peripheral vasculature during ventricular diastole. Thus to adequately determine not only steady-state but also dynamic changes in TPR, it is necessary to first obtain an estimate for *C*_{a}. With this in mind, it is possible to statistically fit the windkessel model described in Fig. 1 directly to measured CO and P_{a} signals sampled at frequency *f*_{s}, where CO is chosen as the input and P_{a} is taken as the response. For that purpose, it is first necessary to adequately transform the differential equation that characterizes the model dynamics from continuous time to discrete time. For simplicity purposes, fluctuations in P_{c}, which are very small compared with fluctuations in P_{a}, shall be considered negligible. Accordingly, *Eq. 2* can be rewritten as (34) Subsequent substitution of *Eq. 34* into *Eq. 1* yields (35) where τ_{eff} represents the effective time constant characterizing the dynamic behavior of the system. Adequate transformation of *Eq. 35* from continuous time to discrete time (see *Conversion From Continuous-Time to Discrete-Time Systems* in the appendix) results in the linear constant-coefficient difference equation of the form (36) where *k* represents the time sample and *E*(*k*) represents the residual error, and the estimation parameters *X*_{1} and *X*_{2} are chosen so that the least squared error between fitted and measured P_{a} is minimized. As a result, *C*_{a} may then be determined from the estimated *X̂*_{1} and *X̂*_{2} as (37) The discrete-time difference equation represented in *Eq. 36* sets out the mathematical relationship among the sampled signals CO, P_{a}, and P_{RA}, the model parameters TPR and *C*_{a}, and the sampling frequency *f*_{s}. Hence, a dynamic TPR(*k*) that fluctuates between *k* samples but that changes at a slower time rate than the sampled cardiovascular signals can be obtained by solving *Eq. 36* for TPR as (38) where *C*_{a} is the previously estimated arterial compliance from *Eq. 37*.

### Important Considerations

One of the main goals of this study is to introduce a robust mathematical method able to determine from measured CO, P_{RA}, and P_{a} signals not only steady-state changes in TPR but also short-term TPR fluctuations caused by baroreflex and autoregulatory modulation of TPR. With this in mind, it is necessary to demonstrate to what extent sampling effects and omission of arterial elasticity can affect the estimation of the desired short-term TPR fluctuations when the differential equation (*Eq. 1*) corresponding to a dynamic model of the circulation based on sound and relatively simple physical principles is further simplified (*Eq. 35*) and adequately transformed from continuous time to discrete time (*Eq. 36*). Solving *Eq. 36* for TPR leads to *Eq. 38*, which defines the present value of TPR as a function of *f*_{s}, present and past values of the sampled CO, P_{RA}, and P_{a} signals, and the previously estimated *C*_{a} from *Eq. 37*. As a result, the relationship represented in *Eq. 38* is defined over two equidistant time samples and describes TPR fluctuations as a function of *f*_{s} and *C*_{a}. Figure 7*A* illustrates the effectiveness of *Eq. 38* to accurately determine short-term TPR fluctuations from P_{a}, P_{RA}, and CO fluctuations computationally generated in closed loop. The true TPR depicted in red are the actual short-term TPR fluctuations generated by the closed-loop computational model of the heart and circulation, whereas the calculated TPR depicted in blue are the short-term TPR fluctuations determined via *Eq. 38* using the computational model CO, P_{RA}, and Pa signals generated in closed loop and resampled to 0.5 Hz and the estimated *C*_{a} from *Eq. 37*. This method not only provides for an accurate representation of ΔTPR but also for the correct values in and *C*_{a}. Despite this method's effectiveness (actual and calculated TPR in Fig. 7*A* are almost indistinguishable), this mathematical approach has to be taken with caution as to make sure that the measured signals are resampled to sufficiently low frequencies after appropriate antialiasing filtering (e.g., 10th-order FIR filter with cutoff frequency *f*_{c} ≤ *f*_{s}/2). The sampling frequency *f*_{s} must be very carefully chosen because *f*_{c} > 1/(2πτ_{eff}) yields rather inaccurate results. Figure 7*B* illustrates the case of incorrect resampling with a comparison in time (*top*) and frequency (*bottom*) between actual (red curves) and calculated TPR determined via *Eq. 38* (blue curves) using the computational model CO, P_{RA}, and P_{a} signals resampled to 1 Hz (see text below for description of gray curves). The energy observed at *f*_{resp} and the spectral leakage from the HR frequency component are present only in the calculated TPR and not in the actual TPR fluctuations, which have a bandwidth [*f* < *f*_{br}=1/(2πτ_{br}) ≈ 0.016Hz] limited by the much lower characteristic frequency of the baroreflexes and autoregulation, which modulate the actual TPR fluctuations and which, in this particular example, are solely responsible for their being. These artifactual fluctuations in calculated TPR arise because *Eq. 38* determines TPR from signals resampled to 1 Hz, thereby allowing for the frequency peak at the respiratory frequency *f*_{resp}=0.4 Hz, which is quite significant in CO, P_{RA}, and P_{a} signals, and the spectral leakage from the very high energy content at the HR frequency to be transmitted into the calculation of TPR. Nevertheless, the representation of TPR at the lower frequencies [*f* < 1/(2πτ_{eff}) ≈ 0.1 Hz] remains extremely accurate and well separated from the artifactual high-frequency content, as illustrated by the almost indistinguishable difference in power spectral densities between the red and blue curves at frequencies lower than 0.1 Hz. Accordingly, resampling to sufficiently low frequencies after appropriate antialiasing filtering provides for an accurate representation of the actual TPR fluctuations by the calculated TPR determined via *Eq. 38*.

In contrast to the presented method for TPR determination via *Eq. 38*, the calculation of TPR simply as the ratio of P_{a} and CO, as illustrated by the gray curves in Fig. 7B, yields grossly inaccurate results. A comparison in frequency (Fig. 7*B*, *bottom*) clearly illustrates how the residual error spreads out along all frequencies, thereby corrupting the low frequencies as well. The yellow area represents to some extent the error caused by the disregarding of *C*_{a} in the calculation of TPR; however, it is for the most part due to the disregarding of low-frequency P_{RA} fluctuations, whereas the error represented by the gray area is attributable almost in its entirety to the disregarding of *C*_{a}. Furthermore, filtering of CO and P_{a} to very low frequencies (*f*_{c}=0.02 Hz) still inevitably results in an inaccurate time representation of TPR as illustrated in Fig. 8, *A* and *B*, by the gray curves. In contrast, TPR fluctuations calculated via *Eq. 38* depicted in blue remain, as expected, almost indistinguishable from the actual TPR fluctuations represented in red. The step changes depicted in Fig. 8*A* demonstrate how the disregarding of P_{RA} fluctuations in the calculation of TPR produces a substantial degree of steady-state error because induced steady-state changes in P_{a} caused by pacing, for example, necessarily lead to opposite steady-state changes in P_{RA}, which contribute significantly to ΔTPR. The error due to *C*_{a}, mostly visible at the beginning of a sudden change, can be rightly described as an unavoidable dynamic error and remains quite evident despite aggressive low-pass filtering of the CO and P_{a} signals. Calculation of TPR as P_{a}/CO − P_{RA} corrects for the steady-state error but does not eliminate the dynamic error clearly illustrated in Fig. 8*B* by the random TPR fluctuations depicted in gray. As a result, reliable short-term TPR fluctuations cannot be accurately assessed from measured CO, P_{RA}, and P_{a} signals if the elasticity of the arteries is omitted in their estimation. It is not in the least surprising that the actual TPR fluctuations differ so much when components of the differential equation (*Eq. 1*) used to generate the actual fluctuations are omitted in the calculation of TPR=P_{a}/CO; however, it is essential to recognize which elements play a significant role. For example, disregarding of *C*_{v} does not seem to play any role in the calculation of TPR, whereas the disregarding of *C*_{a} clearly does. Nevertheless, it is still a common practice to estimate TPR fluctuations as TPR=P_{a}/CO and hence anticipated when researchers are not able to successfully formulate a system identification approach for the quantification of the baroreflexes using TPR fluctuations simply as TPR=P_{a}/CO.

The impact of the presented mathematical method for TPR determination via *Eq. 38* can be best appreciated in Fig. 9, where the basic windkessel model of the circulation characterized by *Eq. 35* is depicted together with experimentally acquired signals of CO, P_{RA}, and P_{a} depicted in red. CO was measured during a preliminary animal experiment via an ultrasonic flow probe placed around the aortic root of a sheep externally paced at ∼130 beats/min while P_{a} and P_{RA} were measured via strain-gauge-based pressure transducers connected to catheters placed in the descending aorta and right atrium. Figure 9*A* shows a 20-s comparison between measured P_{a} and the predicted P_{a} depicted in blue. In this example, the predicted P_{a} assumes a constant TPR obtained from the mean value of *Eq. 38* and disregards the contributions of ΔTPR to P_{a} fluctuations. As a result, the predicted P_{a} cannot account for the very low-frequency fluctuations in measured P_{a} caused by changes in the actual TPR yielding a correlation coefficient of *r*=0.89. Nevertheless, direct visual inspection of measured and predicted P_{a} during a smaller time window (∼2 s), where TPR can be assumed to remain constant, suffices to conclude that the predicted P_{a} describes very well the low-frequency characteristics of the experimentally measured P_{a} observed during systole and diastole over a few cycles, which encompass many *k* data samples of the measured signals sampled at 100 Hz. In contrast, Fig. 9*B* displays the block diagram representation of the windkessel model depicted in Fig. 9*A* together with a 15-min comparison between measured P_{a} resampled to 0.5 Hz and predicted P_{a} when ΔTPR is determined via *Eq. 38* and the contributions of ΔTPR to P_{a} are taken into account. The observed steady-state changes in the measured signals were induced via manipulations of pacing rate and venous return as to generate significant steady-state changes in TPR. Consequently, it is possible to clearly illustrate the difference between the blue curve depicting the case where P_{a} is predicted assuming TPR remains constant (*r* ≈ 0.7) and the black curve exemplifying the case where the contributions of ΔTPR to P_{a} fluctuations are taken into account (*r* ≈ 1.0). As a result, we can rightly assume that the estimated model elements and *C*_{a} offer an adequate characterization for the dynamic open-loop transfer relations *K*_{CO→Pa}, *K*_{PRA}_{→Pa}, and *K*_{TPR→Pa} from *Eq. 9* (depicted here in their time-domain form of step-response function representations), which respectively represent the immediate hydraulic effects of CO, P_{RA}, and TPR on P_{a} fluctuations. In conclusion, the differential equation (*Eq. 35*) corresponding to the very basic dynamic two-element ( and *C*_{a}) windkessel model of the systemic circulation illustrated in Fig. 9*A* is good enough to explain the aortic pressure waveforms observed during systole and diastole and suffices to account for the entirety of the observed low-frequency (*f* < 0.25 Hz because *f*_{s}=0.5 Hz) fluctuations in P_{a} when *Eq. 38* is employed to determine ΔTPR.

## CARDIOVASCULAR SYSTEM IDENTIFICATION

Up until now, we have focused on the analytic algebraic analysis of the dynamics of the systemic vasculature composed of arteries, veins, and its underlying physiological control mechanisms in a manner independent of system identification. This forward model approach yields analytically defined impulse-response function representations constructed from basic physical laws and other well-established relationships, where the model parameters can be basically viewed as vehicles that reflect physical considerations in the system. By contrast, system identification is an inverse modeling technique (9, 21, 39) intended to be applied to the analysis of experimentally acquired data such as in the companion paper (2), where the model parameters and model orders are not known a priori. System identification can be divided into three major tasks or components: *1*) a data record, *2*) a choice of a model structure, and *3*) determination of the best model in the set, guided by the data. The presented mathematical analysis incorporates physical insight into a model structure (in the appendix, we readdress the mathematical analysis with a focus on measured data composed of sample data rather than a continuum of times, whether experimentally acquired or computationally generated), whereas the present section is concerned with the determination of the best model in the set, guided by the data.

### Hemodynamic System Identification

Dynamic closed-loop effects of CO and P_{RA} on P_{a} can be modeled via fluctuations in CO and P_{RA} as illustrated in Fig. 5*B* via the autoregressive exogenous model represented by the linear constant-coefficient difference equation of the form (39) where Δ denotes the percent fluctuations in that signal around its mean value indicated by the overbars and *e*_{Pa} is the residual error. The parameter values of the *a*_{i}, *b*_{1i}, and *b*_{2i} coefficients and the model orders *n*, *m*_{1}, and *m*_{2} are determined by the linear least-squares minimization of the residual error in conjunction with Rissanens's minimum description length (MDL) principle (31), a model order selection criterion that evaluates a given model's performance compared with other models. Once these parameters are determined, the transfer relations *H*_{CO→Pa} and *H*_{PRA}_{→Pa} are fully defined. In addition, the estimated *G*{*H*_{CO→Pa}} and *G*{*H*_{PRA}_{→Pa}} can then be used to indirectly determine *G*{*H*_{Pa}_{→TPR}} from *Eq. 32* and *G*{*H*_{PRA}_{→TPR}} from *Eq. 33* as (40) and (41) respectively, where *G*{*K*_{TPR→Pa}}=*G*{*K*_{CO→Pa}}=TPR·CO/P_{a} and *G*{*K*_{PRA}_{→Pa}}=100%/P_{a}.

### Regulatory System Identification

Dynamic autonomic closed-loop control of TPR by the arterial and cardiopulmonary baroreceptors in the presence of local vascular autoregulation can be modeled via fluctuations in P_{a} and P_{RA} as illustrated in Fig. 6*A*. The modulation of TPR can be described by the autoregressive exogenous model represented by the linear constant-coefficient difference equation of the form (42) where *e*_{TPR} is the residual error. The parameter values of the *c*_{i}, *d*_{1i}, and *d*_{2i} coefficients, the model orders *p*, *q*_{1}, and *q*_{2}, and the system delays *nk*_{1} and *nk*_{2} are determined by the linear least-squares minimization of the residual error in conjunction with the MDL principle. Once these parameters are determined, the transfer relations *H*_{Pa}_{→TPR} and *H*_{PRA}_{→TPR} are fully defined.

### Monte Carlo Simulations

System identification requires the input signals to be poorly correlated and sufficiently broadband so that all the modes of the system to be identified are excited and reliable (21, 39). Consequently, the computational model of the heart and circulation was simultaneously excited by two separate sources as to simulate an orthogonal input design in which HR and venous return were independently varied with frequency band limited to 0.1 Hz about their mean values in a nearly uncorrelated fashion while CO, P_{a}, P_{RA}, and TPR were measured. To evaluate the effectiveness of hemodynamic (HSI) and regulatory system identification (RSI) to quantitatively characterize the dynamic closed-loop transfer relations *H*_{CO→Pa}, *H*_{PRA}_{→Pa}, *H*_{Pa}_{→TPR}, and *H*_{PRA}_{→TPR}, the uncertainty of the estimates was determined by calculating the mean and standard error for the impulse-response estimates from 100 different realizations of computer-generated data sampled at *f*_{s}=0.5 Hz. Figure 10*A* displays a graphical comparison between the analytically derived solutions for the closed-loop transfer relations *H*_{CO→Pa} and *H*_{PRA}_{→Pa} depicted in circles in their step-response function representation and mean HSI results depicted in squares. For illustrative purposes, the analytically derived solutions for the open-loop transfer relations *K*_{CO→Pa} and *K*_{PRA}_{→Pa} depicted in triangles are also shown. Figure 10*B* displays a graphical comparison between the analytically derived solutions for *H*_{Pa}_{→TPR} and *H*_{PRA}_{→TPR} depicted in circles in their step-response function representation and mean RSI results depicted in squares. Table 1 displays a numerical comparison between the actual and estimated static gain values obtained via HSI and RSI.

### Summary and Conclusions

In this study, we emphasized the analytic algebraic analysis of the dynamics of the peripheral vasculature composed of arteries, veins, and its underlying physiological control mechanisms of baroreflex and autoregulatory modulation of TPR. We acknowledge the existence of far more complex and complicated models of the circulation. Our goal, however, was to enhance our understanding of the crucial functional relationships that determine the behavior of the systemic circulation and its underlying physiological regulatory mechanisms with minimal modeling. Ultimately, we hope that the presented analytic analysis simultaneously simplified and deepened our understanding of the cardiovascular system. Furthermore, we readdressed the mathematical analysis with a focus on measured data composed of sample data rather than a continuum of times, whether experimentally acquired or computationally generated, and highlight the significance of this practical issue. As a result of this analysis, we developed a novel mathematical method to determine short-term TPR fluctuations, which accounts for the entirety of observed P_{a} fluctuations. Finally, we directed our attention to the development and evaluation of quantitative tools, which could be subsequently utilized to delineate the actual actions of the physiological mechanisms responsible for the couplings among CO, P_{a}, P_{RA}, and TPR without altering the underlying regulatory mechanisms of baroreflex and autoregulatory modulation of TPR. As a result, we proposed a novel cardiovascular system identification method, introduced as two separate model structures referred to as HSI and RSI, able to quantitatively characterize the independent dynamic closed-loop contributions of CO and P_{RA} to P_{a} fluctuations as well as the independent dynamic closed-loop contributions of P_{a} and P_{RA} to short-term TPR fluctuations.

## APPENDIX

Adding feedback around an existing system yields a new system with different pole-zero locations than the original system and thus a different dynamic behavior. Hence, the poles of the closed-loop systems *H*_{CO→Pa}(*s*) and *H*_{PRA}_{→Pa}(*s*) introduced in *Eq. 10* depend on the dynamic effect of the feedback system that characterizes the arterial baroreflex, *K*_{Pa}_{→TPR}(*s*), which moves the locations of the poles and zeros of the original open-loop systems *K*_{CO→Pa}(*s*) and *K*_{PRA}_{→Pa}(*s*). Consequently, *Eqs. 16* and *17* demonstrate that the static gain *G*{*K*_{Pa}_{→TPR}}, time delay *T*_{br}, and characteristic frequency 1/τ_{br} of the arterial baroreflex determine the poles of *H*_{CO→Pa}(*s*) and *H*_{PRA}_{→Pa}(*s*). For example, if the amount of feedback in *Eq. 16* were negligible, i.e., no arterial baroreflex modulation of TPR were present in *H*_{CO→Pa}(*s*), the poles would simply be the characteristic frequencies of the original system, 1/τ_{1} and 1/τ_{2}; however, if feedback cannot be neglected, the poles are not simply 1/τ_{1}, 1/τ_{2}, and 1/τ_{br} because the exact pole-zero locations depend on *G*{*K*_{Pa}_{→TPR}} and *T*_{br} as well. To further illustrate the legality of this practical matter, it is necessary to convert the continuous-time system functions *H*_{CO→Pa}(*s*) and *H*_{PRA}_{→Pa}(*s*) to their equivalent discrete-time system functions *H*_{CO→Pa}(*z*) and *H*_{PRA}_{→Pa}(*z*), respectively.

### Conversion from Continuous-Time to Discrete-Time Systems

To manipulate signals composed of a sequence of samples rather than a continuum of times, it is necessary to convert the continuous-time system function *F*(*s*) to its equivalent discrete-time system function *F*(*z*). This can be accomplished by the bilinear transform, which is an algebraic transformation between the complex frequency variables *s* and *z* that maps the entire imaginary axis in the analog *s*-plane onto one circumference of the unit circle in the discrete-time *z*-plane (27). This transformation, which also arises from applying the trapezoidal integration rule to the differential equation corresponding to *F*(*s*) (16), corresponds to replacing *s* by (A1) where *f*_{s} denotes the sampling frequency. As a result, a stable continuous-time system *F*_{s}(*s*) will always map to a stable discrete-time system *F*_{s}(*z*) and the frequency response of *F*_{s}(*s*) will be exactly replicated in the frequency response of *F*_{s}(*z*).

Bilinear transformation of the continuous-time closed-loop system *H*_{CO→Pa}(*s*) in *Eq. 16* results in the equivalent discrete-time frequency-domain representation (A2) where the system order *n* (and consequently the number of poles) depends on the intrinsic value of the feedback parameter *T*_{br} and the choice of *f*_{s} because *n*=3 + *f*_{s}*T*_{br}. Thus by rearranging *Eq. A2* as increasing factors of *z*^{−1} and deliberately selecting the sampling frequency so that *f*_{s}*T*_{br}=1, we arrive at (A3) where the feedback parameter *G*{*K*_{Pa}_{→TPR}} embedded in β_{1} determines the exact pole-zero locations of the closed-loop system *H*_{CO→Pa}(*z*). Likewise, bilinear transformation of the continuous-time closed-loop system *H*_{PRA}_{→Pa}(*s*) in *Eq. 17* results in the equivalent discrete-time frequency-domain representation (A4) Thus by rearranging *Eq. A4* for increasing factors of *z*^{−1} and deliberately selecting the sampling frequency so that *f*_{s}*T*_{br}=1, we arrive at (A5) As a result, the equivalent discrete-time frequency-domain representation of *Eq. 10* becomes (A6) Multiplication of *Eq. A6* by the denominator polynomial factors where multiplication in discrete-time frequency with *z*^{−k} corresponds to a discrete-time shift by *k* samples (27) and subsequently solving for P_{a} finally transforms *Eq. A6* into the linear fourth-order constant-coefficient difference equation of the form (A7)

Bilinear transformation of the continuous-time transfer relations *K*_{Pa}_{→TPR}(*s*) and *K*_{PRA}_{→TPR}(*s*) characterized by the first-order LTIsystems in *Eqs. 11* and *12* yields (A8) and (A9) respectively. Consequently, by deliberately selecting the sampling frequency so that *f*_{s}*T*_{br}=1, we arrive at the equivalent discrete-time frequency-domain representation of *Eq. 8* expressed as (A10) where multiplication in discrete-time frequency with *z*^{−k} corresponds to a discrete-time shift by *k* samples, and hence *Eq. A10* finally transforms into the linear constant-coefficient difference equation of the form (A11)

### Effects of Local Vascular Autoregulation

Bilinear transformation of the continuous-time transfer relation *K*_{CO→TPR}(*s*) characterized by the first-order LTI system in *Eq. 26* yields (A12) Similarly, bilinear transformation of the continuous-time transfer relation *H*_{CO→Pa}(*s*) in *Eq. 25* characterized by the third-order closed-loop system (A13) results in (A14) where γ=β_{1}*G*{*K*_{CO→TPR}}/*G*{*K*_{Pa}_{→TPR}} and the system order *n* depends upon the intrinsic value of the feedback parameter *T*_{br} and the choice of *f*_{s} because *n*=3 + *f*_{s}*T*_{br}. Thus by rearranging *Eq. A14* as increasing factors of *z*^{−1} and deliberately selecting for simplicity purposes the sampling frequency *f*_{s}=1/*T*_{br}=1/*T*_{ar}, we arrive at (A15) As a result, the equivalent discrete-time frequency-domain representation of *Eq. 25* becomes (A16) where multiplication in discrete-time frequency with *z*^{−k} corresponds to a discrete-time shift by *k* samples, and consequently *Eq. A16* finally transforms into the linear fourth-order constant-coefficient difference equation of the form (A17) Substitution of *H*_{PRA}_{→Pa}(*z*), *K*_{Pa}_{→TPR}(*z*), *K*_{PRA}_{→TPR}(*z*), *K*_{CO→TPR}(*z*), and *H*_{CO→Pa}(*z*) from *Eqs. A5*, *A8*, *A9*, *A12*, and *A14*, respectively, into the equivalent discrete-time frequency-domain representation of *Eq. 28* expressed as (A18) results in the discrete-time characterization of the independent dynamic closed-loop contributions of P_{a} and P_{RA} on TPR as (A19) and (A20) respectively, where the system order *p* depends on the intrinsic value of the time delay *T*_{ar} and the choice of *f*_{s} because *p*=3 + *f*_{s}*T*_{ar}. Thus by rearranging *Eqs. A19* and *A20* as increasing factors of z^{−1} and deliberately selecting for simplicity purposes *f*_{s}=1/*T*_{br}=1/*T*_{ar}, we arrive at the equivalent discrete-time frequency-domain representation of *Eq. 28* as (A21) where multiplication in discrete-time frequency with z^{−k} corresponds to a discrete-time shift by *k* samples, and thus *Eq. A21* finally transforms into the linear fourth-order constant-coefficient difference equation of the form (A22) *Equations A7*, *A11*, *A17*, and *A22* correspond to a standard autoregressive exogeneous parametric model structure (21). The model is autoregressive because the output looks back on past values of itself and exogenous because it possesses external inputs.

## GRANTS

This work was sponsored by the United States National Aeronautics and Space Administration through a grant from the National Space Biomedical Research Institute and a grant from the Center for the Integration of Medicine and Innovative Technology.

## 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 © 2004 by the American Physiological Society