## Abstract

We propose two identification algorithms for quantitating the total peripheral resistance (TPR) baroreflex, an important contributor to short-term arterial blood pressure (ABP) regulation. Each algorithm analyzes beat-to-beat fluctuations in ABP and cardiac output, which may both be obtained noninvasively in humans. For a theoretical evaluation, we applied both algorithms to a realistic cardiovascular model. The results contrasted with only one of the algorithms proving to be reliable. This algorithm was able to track changes in the static gains of both the arterial and cardiopulmonary TPR baroreflex. We then applied both algorithms to a preliminary set of human data and obtained contrasting results much like those obtained from the cardiovascular model, thereby making the theoretical evaluation results more meaningful. This study suggests that, with experimental testing, the reliable identification algorithm may provide a powerful, noninvasive means for quantitating the TPR baroreflex. This study also provides an example of the role that models can play in the development and initial evaluation of algorithms aimed at quantitating important physiological mechanisms.

- baroreceptors
- autonomic nervous system
- hemodynamics
- system identification
- mathematical modeling

the arterial and cardiopulmonary baroreflex systems attempt to maintain blood pressures on a time scale of seconds to minutes by dynamically controlling heart rate (HR), ventricular contractility (VC), total peripheral resistance (TPR), and systemic venous unstressed volume (SVUV) via autonomic efferent pathways (7). Among these four distinct baroreflex control pathways, the HR baroreflex is by far the most extensively studied and understood (12, 13) due to the relative simplicity of measuring HR. Although the circulatory baroreflex pathways are less understood, they may be more significant to arterial blood pressure (ABP) regulation than the cardiac baroreflex pathways. For example, according to Guyton and Hall (7), venous return is nearly saturated at nonelevated right atrial (RA) pressures, and thus ABP could not be substantially increased by enhancing cardiac function. The TPR baroreflex, in particular, may be the most important short-term contributor to ABP regulation because TPR affects ABP directly via Ohm's law and indirectly via venous return. Because of the importance of the TPR baroreflex, several invasive techniques have been previously developed for its study in animal preparations (e.g., Refs. 20 and 23). However, to our knowledge, there is currently no noninvasive technique available that could be applied to humans. Such a technique could advance the basic understanding of blood pressure regulation during normal and pathological conditions (e.g., peripheral autonomic neuropathies due to diabetes mellitus or Parkinson's disease) and under various environmental factors (e.g., microgravity). A noninvasive technique for quantitating the TPR baroreflex could also potentially be employed to guide therapy in patients with symptomatic orthostatic hypotension.

To this end, we have developed a general approach for quantitatively characterizing autonomic cardiovascular control mechanisms by mathematically analyzing the beat-to-beat fluctuations in multiple cardiorespiratory signals over intervals of ∼6 min. The mathematical analysis is based on the methods of system identification (11), which aim to determine the dynamic transfer properties that couple the fluctuations between the measured signals rather than simply quantitating the measured fluctuations (e.g., power spectral analysis). To obtain a complete characterization over the physiological range of frequencies, the approach includes a broadband excitation protocol in which subjects breathe according to a sequence of randomly spaced auditory tones (1). Importantly, this excitation protocol does not significantly disrupt normal system operation (1). In contrast, conventional techniques for studying autonomic control mechanisms perturb the cardiovascular system in a more nonphysiological manner (e.g., lower body negative pressure and neck chamber suction) to elicit a compensatory response (12,13). We have found that when the subjects are following this random-interval breathing protocol and are otherwise at rest that the measured fluctuations are sufficiently small and stationary such that the autonomic coupling mechanisms may be represented by linear, time-invariant (LTI) transfer functions (5). We have previously applied this general approach to identify the impulse responses (intuitive time-domain representations of transfer functions) characterizing the HR baroreflex and other important physiological mechanisms from noninvasively measured fluctuations in HR, ABP, and instantaneous lung volume (ILV) (16-18).

To apply this general approach for the quantitation of the TPR baroreflex, we formulated the block diagram shown in Fig.1, which is based on the work of Raymundo et al. (20) and specifies the particular coupling mechanisms we seek to identify. This block diagram includes the feedback pathways of both the arterial and cardiopulmonary baroreflex arcs, which respectively couple ABP fluctuations to TPR fluctuations (arterial TPR baroreflex) and RA transmural pressure (TP) fluctuations to TPR fluctuations (cardiopulmonary TPR baroreflex). The block diagram also includes the perturbing noise source N_{TPR}, which reflects the residual variability in TPR not accounted for by the two baroreflex coupling mechanisms. Such variability may be due to, for example, the autoregulation of local vascular beds. Ideally, one would want to obtain beat-to-beat measurements in ABP, RATP, and TPR to identify the impulse responses charactering the arterial TPR baroreflex and cardiopulmonary TPR baroreflex as well as the power spectrum of N_{TPR}. However, in practice, techniques for directly measuring beat-to-beat fluctuations in TPR are not available. Furthermore, although it is possible to measure RATP in humans via a central venous catheter and an esophageal balloon [for intrathoracic pressure (ITP)] (14), these measurement techniques are invasive and typically do not reliably account for RATP on a beat-to-beat basis (e.g., swallowing artifacts).

In this paper, we investigate the possibility of quantitating the TPR baroreflex mechanisms shown in Fig. 1 with beat-to-beat measurements of TPR and RATP unavailable for analysis. We specifically propose two system identification algorithms that require only beat-to-beat measurements of ABP and left ventricular flow rate [cardiac output (CO)]. Both of these measurements may be obtained noninvasively in humans with, for example, Finapres (10) and Doppler ultrasound (6) techniques, respectively. As will be shown, the tradeoff in considering only these measurements is that the impulse responses to be identified reflect the dynamic properties of other distinct physiological mechanisms in addition to TPR baroreflex mechanisms. Importantly, however, the static gains (integral of the impulse responses) characterizing each TPR baroreflex mechanism may be computed from the identified impulse responses through the formulation of physiological models. These static gains specifically indicate the steady-state TPR change that would occur if each TPR baroreflex mechanism were stimulated by a unity step increase in its sensed pressure. It may also be possible to recover additional quantitative information characterizing the TPR baroreflex mechanisms from the identified impulse responses with further assumptions.

We evaluated the performance of the two system identification algorithms with respect to a realistic test bed of beat-to-beat variability generated from a previously developed computational model of the human cardiovascular system (16). In contrast to an animal model, an independent measure of the impulse responses to be identified can be easily obtained in the computational model by applying, according to mathematical definition, an arbitrarily narrow, unit-area input to the appropriate point in the model and then measuring the output response of interest while all other perturbations to the output are held constant. These impulse responses may then be regarded as the gold standards against which the corresponding impulse responses, which are estimated from the model beat-to-beat variability, may be compared. The computational model also provides a powerful means to analyze the sensitivity of the system identification algorithms. That is, it is possible to determine how much the dynamic properties of the computational model would have to be altered before we would see a corresponding change in the estimates. In addition to this theoretical evaluation, we also applied both system identification algorithms to a preliminary set of noninvasively measured signals obtained from 10 healthy human subjects in the supine posture who were then tilted to 30° with respect to the horizontal position.

## IDENTIFICATION ALGORITHMS

The two system identification algorithms for quantitating the TPR baroreflex that we propose here are conceptually different in terms of their approach for accounting for TPR fluctuations from the fluctuations in ABP and CO, the signals that are assumed to be available for analysis. The approach of the first algorithm is the most obvious and involves estimating TPR fluctuations directly from the available signals. We refer to this algorithm as direct identification. The approach of the second algorithm is based on the concept that the dynamic relationship between fluctuations in ABP, CO, and RATP indirectly reflect the fluctuations in TPR that are caused by the TPR baroreflex. We refer to this algorithm as indirect identification.

The two system identification algorithms are similar in that they each assume that the left ventricular stroke volume (SV) signal, which may be computed from the available CO signal, is an acceptable surrogate for the unavailable RATP signal. Some experimental evidence that supports this assumption are as follows: *1*) steady-state SV is determined exclusively by steady-state RATP provided that average ABP does not exceed ∼180 mmHg (9); *2*) pulmonary ABP is essentially constant due to recruitment and distension (3, 9); and *3*) ventricular contractility changes little during stable, resting conditions (8, 21). For small fluctuations, SV variability is specifically assumed to be completely accounted for by the fluctuations in RATP through a LTI relationship. The impulse response that precisely couples RATP fluctuations to SV fluctuations characterizes the input-output relationship of the heart-lung unit (9). The static gain of the heart-lung unit is equal to the right ventricular diastolic compliance (9). By further assuming that this impulse response is invertible, RATP fluctuations may be completely determined from the convolution of SV fluctuations with an impulse response, which may be thought of as characterizing the inverse, dynamic properties (output-input relationship) of the heart-lung unit (inverse heart-lung unit). The static gain of the inverse heart-lung unit is thus equal to the reciprocal of the right ventricular diastolic compliance. Importantly, however, when SV and RATP fluctuations are normalized by their respective mean values (as will be the case, henceforth, for all considered signals), the static gain of the inverse heart-lung unit is always equal to one and is no longer dependent on the diastolic properties of the right ventricle. (See Ref. 16 for a more formal, mathematical justification of the assumptions presented here.)

### Direct Identification

The block diagram shown in Fig. 2illustrates the coupling mechanisms that the direct identification algorithm seeks to identify. This block diagram is derived from the ideal block diagram shown in Fig. 1 by substituting actual TPR with estimated TPR (
) and RATP with SV. As a consequence of the latter substitution, SV→TPR, which couples SV fluctuations to TPR fluctuations (assuming they are perfectly estimated), is sought to be identified rather than the desired cardiopulmonary TPR baroreflex. However, SV→TPR encompasses the dynamic properties of the cardiopulmonary TPR baroreflex as well as the inverse heart-lung unit through a cascade combination, as depicted in the physiological model shown in Fig.3. Moreover, the static gain of SV→TPR is equal to that of the cardiopulmonary TPR baroreflex, because the fluctuations in the considered signals are normalized by their respective mean values (as discussed above). Because of the signal normalization, the static gain of SV→TPR as well as the arterial TPR baroreflex are unitless quantities. So, for example, with the unitless static gain of SV→TPR, one could determine the steady-state percent change in TPR with respect to its mean value that would occur if cardiopulmonary TPR baroreflex were stimulated by an *X*% step change in RATP with respect to its mean value by multiplying*X* with the unitless static gain value.

The block diagram shown in Fig. 2 is mathematically represented by an autoregressive moving average (ARMA) difference equation—a highly flexible subclass of LTI systems—of the following form
Equation 1
where *t* is discrete time; the terms *m*,*n*, and *p* limit the number of terms in the equation (model order); and W_{TPR} is the unmeasured, residual error. The three sets of parameters {*a _{i}
*,

*b*,

_{i}*c*} completely specify the impulse responses characterizing the arterial TPR baroreflex and SV→TPR, whereas the residual error together with the set of parameters {

_{i}*a*} fully define the power spectrum of N

_{i}_{TPR}(11). The values of the parameters are determined from 6-min segments of zero-mean ABP, SV, and normalized by their mean values and sampled at 0.5 Hz by minimizing the variance of the residual error in conjunction with an ARMA parameter reduction algorithm (19). These three signals are determined from continuous records of CO and ABP measurements as follows.

Among the potential techniques for determining from ABP and CO signals, we choose a straightforward method that is relatively robust to the high-frequency wave reflections in noninvasively measured, peripheral ABP signals. This method specifically determines the fluctuations in by computing the ratio of average ABP to average CO over intervals in which TPR changes little and the net flow through the large, compliant arteries is small. It should be possible to choose such intervals, because the dominant time constant of the systemic arteries [∼2 s (22)], which is established by the product of average TPR and total arterial compliance, is smaller than the time constant governing changes in TPR [∼10 s (2)], which is established by the relatively slow α-sympathetic nervous system. We specifically estimated a value of TPR for each cardiac cycle by computing the ratio of average ABP to average CO over the interval that includes the five previous and five subsequent cardiac cycles. We then formed a stepwise continuous process whose value corresponds to the estimated TPR value of the current cardiac cycle for the time period of that cycle. We finally sampled the stepwise continuous process to 0.5 Hz with an antialiasing filter whose impulse response is a unit-area boxcar of 4-s duration. The ABP and SV signals that are utilized for identification were similarly processed. The signals were determined for each cardiac cycle by respectively averaging and integrating the continuous measurements of ABP and CO over the five previous and five subsequent cardiac cycles (and dividing the integrated CO by 11). Stepwise continuous processes were then analogously formed and sampled to 0.5 Hz.

### Indirect Identification

The block diagram shown in Fig. 4illustrates the coupling mechanisms that the indirect identification algorithm seeks to identify. As will be shown below via physiological models, these coupling mechanisms reflect the dynamic properties of the coupling mechanisms in the ideal block diagram shown in Fig. 1 that we are ultimately aiming to quantitate.

CO→ABP, which couples CO fluctuations to ABP fluctuations, encompasses the dynamic properties of the arterial TPR baroreflex as well as the systemic arterial tree according to the physiological model shown in Fig. 5
*A*. (This model and that shown in Fig. 5
*B* assume that the considered ABP variability is small and primarily generated by only the depicted physiological mechanisms.) As its name suggests, the systemic arterial tree characterizes the mechanical properties of the systemic arteries and specifically couples CO fluctuations to ABP fluctuations as well as TPR fluctuations to ABP fluctuations. The physiological model shown in Fig. 5
*A* indicates that, for example, an increase in CO would initially cause ABP to increase via the systemic arterial tree. This would, in turn, excite the arterial TPR baroreflex/systemic arterial tree arc to decrease TPR so as to maintain ABP. Importantly, the static gain of the arterial TPR baroreflex may be exactly computed from the static gain of CO→ABP, because the static gain of the systemic arterial tree is identical to one due to the normalization of the analyzed signals with their respective mean values (see *Eqs. 2
* and *
3
*). Again, due to this normalization, the static gain of CO→ABP, and thus the arterial TPR baroreflex, will be unitless and may be interpreted analogously to SV→TPR (see above). According to the physiological model shown in Fig. 5
*A*, the static gain of the arterial TPR baroreflex will be computed as a negative value, indicating a negative feedback mechanism provided that the static gain of CO→ABP is <1. Moreover, assuming that the arterial TPR baroreflex does not reduce TPR by a greater percentage (with respect to mean values) than the increase in CO, which would indicate overcompensation, the static gain of CO→ABP will also be >0. Finally, it may also be possible to recover additional dynamic properties characterizing the arterial TPR baroreflex from CO→ABP by making additional assumptions, which are outlined in discussion, *Indirect Identification*.

SV→ABP, which couples SV fluctuations to ABP fluctuations, encompasses the dynamic properties of the arterial TPR baroreflex and cardiopulmonary TPR baroreflex as well as the inverse heart-lung unit and systemic arterial tree according to the physiological model shown in Fig. 5
*B*. This physiological model suggests that, for example, an increase in SV would indicate that an increase in RATP had occurred through the inverse heart-lung unit. This RATP increase would excite the cardiopulmonary TPR baroreflex to decrease TPR, which would then stimulate the arterial TPR baroreflex/systemic arterial tree arc to increase TPR and maintain ABP. The physiological model here may appear to be counterintuitive, because the increase in SV does not cause an increase in CO and thus ABP through the systemic arterial tree. The reason for this is that SV→ABP is mathematically defined to characterize the effects of SV fluctuations on ABP fluctuations while all other considered inputs to ABP fluctuations (which includes CO fluctuations here) are held perfectly constant. This implies that the increase in SV must be accompanied by a commensurate decrease in HR. Importantly, a unitless static gain of the cardiopulmonary TPR baroreflex (which may be interpreted analogously to SV→TPR; see above) may be exactly computed from the unitless static gains of both SV→ABP and CO→ABP (see *Eqs. 2
* and *
4
*). According to the physiological models shown in Fig. 5, the static gain of the cardiopulmonary TPR baroreflex will be computed as a negative value (indicating a negative feedback mechanism) if the static gain of SV→ABP is negative and the static gain of CO→ABP is positive (see above).

In addition to the two coupling mechanisms, the block diagram shown in Fig. 4 also includes the perturbing noise source N_{ABP}. This quantity is unmeasured and represents the residual variability in ABP not attributed to CO and SV fluctuations. Such fluctuations may be due to, for example, the autoregulation of local vascular beds. The block diagram shown in Fig. 4 is also mathematically represented by an ARMA difference equation, which is given as follows
Equation 2
where the terms *q*, *r*, and*s* limit the model order and W_{ABP} is the unmeasured, residual error. The three sets of parameters {*d _{i}
*,

*e*,

_{i}*f*} and the residual error term, which completely specify the impulse responses characterizing CO→ABP and SV→ABP and the power spectrum of N

_{i}_{ABP}, are determined from the measured signals similarly to direct identification (seeidentification algorithms,

*Direct Identification*). According to the physiological models shown in Fig. 5, the static gains of the arterial TPR baroreflex and cardiopulmonary TPR baroreflex may be respectively computed from the three sets of parameters as follows Equation 3and Equation 4

## METHODS OF EVALUATION

We applied the two proposed system identification algorithms to*1*) a realistic test bed of data generated from a cardiovascular model (theoretical evaluation) and *2*) noninvasively measured data obtained from healthy human subjects in the supine posture who were then tilted to 30° with respect to the horizontal position (preliminary experimental analysis).

### Cardiovascular Model

We based the theoretical evaluation on a slightly modified version of a computational model of the human cardiovascular system that we recently developed and demonstrated to generate realistic short-term, beat-to-beat variability (16). We previously employed the model for the similar purpose of theoretically evaluating a system identification technique for quantitatively characterizing the HR baroreflex and other important physiological mechanisms (16).

#### Description.

The block diagram shown in Fig. 6illustrates the major components of the original cardiovascular model. For the present study, a key property of this model is the signal-to-noise ratio (SNR) of TPR fluctuations–the ratio of the SD of the baroreflex-induced TPR fluctuations to the SD of the unmeasurable TPR disturbance (N_{TPR} in Fig. 6). Because the actual SNR of TPR fluctuations is unknown, we arbitrarily set it to a value of ∼5 by adjusting the original value assigned to the SD of N_{TPR}. (The SNR of TPR fluctuations could have been set to as low as ∼2 without significantly altering the results of our study.) Note that this chosen SNR value may be approximately validated if the impulse responses identified from the model-generated data resemble the corresponding impulse responses identified from the preliminary set of experimental human data. To maintain realistic beat-to-beat hemodynamic variability, we also introduced a white disturbance to SVUV that was bandlimited to 0.1 Hz with a SD of 3.125 ml (see Ref. 16for a detailed discussion of this topic).

Because continuous CO is somewhat difficult to measure in practice, we also accounted for measurement noise in the present study. (Note that N_{TPR} is not measurement noise but rather a physiological disturbance to TPR that is also unmeasurable.) We specifically added zero-mean Gaussian white noise to the model-generated CO, ABP, and cardiac cycle length beat sequences with SD equal to 25%, 12.5%, and 12.5%, respectively, of the SDs of the corresponding beat sequence. Importantly, although these SDs were chosen arbitrarily, the conclusions of our study are not affected by their precise values.

#### Gold standard results.

To evaluate the performance of the two proposed system identification algorithms applied to data generated from the cardiovascular model, it is necessary to establish, in a manner independent of system identification, the impulse responses and power spectra of perturbing noise sources characterizing the physiological mechanisms shown in Figs. 1, 2, and 4. These impulse responses and power spectra of perturbing noise sources may then be regarded as the gold standard results against which the corresponding identification results obtained from the model beat-to-beat variability may be compared.

Our approach for establishing the gold standard impulse responses is based on their literal mathematical meaning. In particular, the impulse response is defined to represent the output response of a system to an arbitrarily narrow, unit-area input. (Note that the impulse response also completely characterizes the input-output relationship of the coupling mechanisms here because of our LTI assumption.) Hence, in establishing each of the gold standard impulse responses, we applied an impulse input to the desired point in the cardiovascular model and measured the output response of interest while holding all other perturbations to the output constant. We specifically employ this approach to establish the gold standard impulse responses characterizing the arterial TPR baroreflex, cardiopulmonary TPR baroreflex, systemic arterial tree, and heart-lung unit. We then applied the LTI system theory to the physiological models shown in Figs. 3 and 5 to compute the gold standard impulse responses characterizing SV→TPR, CO→ABP, and SV→ABP. The gold standard power spectra of perturbing noise sources were likewise determined according to their literal mathematical meanings. (See Ref.16 for a detailed description of the specific types of techniques employed for establishing the gold standard results.)

#### Monte Carlo simulations.

Although the impulse responses, as determined by the two system identification algorithms, include a measure of uncertainty in terms of a covariance matrix, this uncertainty measure is only an estimate itself (19). To account more accurately for estimation error variance, we compared the means with 95% confidence intervals for the impulse responses as well as for the power spectra of perturbing noise sources as identified from 20 different realizations of model beat-to-beat variability with the corresponding gold standard results.

### Experimental Data

For the preliminary experimental analysis, we collected a set of data from 10 healthy human subjects [5 men and 5 women, age: 25.2 ± 3.7 yr (means ± SD)]. The experimental protocol was approved by the Massachusetts Institute of Technology Committee on the Use of Humans as Experimental Subjects. Written informed consent was obtained from each subject. All experiments were performed in the Clinical Research Center at the Massachusetts Institute of Technology.

One lead of the surface ECG, ABP, and CO were measured continuously and noninvasively and were interfaced on-line to a microcomputer running a dedicated data collection and analysis program (BVA, Andiamo; Oslo, Norway). The ECG signal was measured with a Hewlett-Packard EKG Monitor 78203A (Andover, MA), and the ABP signal was measured at the finger with a 2300 Finapres Continuous Blood Pressure Monitor (Ohmeda; Englewood, CO). The CO signal was measured according to a previously described Doppler ultrasound technique (6). Briefly, aortic velocity was measured with a bidirectional ultrasound Doppler velocimeter (CFM-750 Vingmed Sound; Horten, Norway), which was operated in pulsed mode at 2 MHz with the hand-held transducer placed on the suprasternal notch. The aortic diameter of the rigid aortic ring was determined in a separate session by parasternal sector-scanner imaging (CFM-750 Vingmed Sound). By assuming that the aortic valve orifice is circular, its area may be computed from the measured diameter. Instantaneous CO can then by calculated from the product of the measured instantaneous maximum velocity and the area of the orifice. This calculation is based on the additional assumptions that the velocity profile in the aortic valve orifice is rectangular and that this velocity is conserved as the central maximum velocity of a jet 3–4 cm upward in the aortic root. The digitized signals (quantized at 12 bits and sampled at 100 Hz) were collected for two ∼6-min time periods in which each subject was following the random-interval breathing protocol (see Introduction). Each subject was in the supine posture during the first 6-min time period and tilted 30° with respect to the supine position during the second 6-min time period. A minimum of 5 min was allowed for hemodynamic equilibration after the change in posture.

## RESULTS

### Direct Identification

Figure 7 illustrates the results of applying the direct identification algorithm to the cardiovascular model. The two direct identification coupling mechanisms, arterial TPR baroreflex and SV→TPR, are each characterized here in terms of a step response, which is mathematically defined to be the running integral of an impulse response. Note that the asymptotic value of the step response (value of the step response at ∼50 s here) represents the static gain of the coupling mechanism. Also, recall fromidentification algorithms, *Direct Identification*that the static gain of the gold standard SV→TPR equals the static gain of the gold standard cardiopulmonary TPR baroreflex. The perturbing noise source N_{TPR} in Fig. 7 is shown in terms of its power spectrum. Figure 7 shows large deviations between the direct identification estimates and their respective gold standards. For example, the estimated dynamics are much faster than the gold standard dynamics, and the estimated arterial TPR baroreflex step response indicates a positive feedback mechanism. That is, this step response erroneously suggests that TPR would increase if a step increase in ABP were applied to the cardiovascular model.

Figure 8 shows the analogous results obtained by applying the direct identification algorithm to the preliminary set of experimental data obtained from the 10 healthy human subjects in the supine posture. Note that the step responses and power spectrum here are similar in structure to the corresponding estimates shown in Fig. 7, which were obtained from the cardiovascular model. By structural similarity, we are referring to the fast, initial transient, sign of initial transient, sign of asymptotic value, and time to reach asymptotic value for the step responses and spectral power concentrated at low frequencies (<0.05 Hz) for the perturbing noise sources. The same structural similarity was also seen when the direct identification algorithm was applied to the experimental data obtained from the human subjects while they were tilted 30° with respect to the supine posture.

### Indirect Identification

Figure 9 illustrates the results of applying the indirect identification algorithm to the cardiovascular model. The two indirect identification coupling mechanisms, CO→ABP and SV→ABP, are also characterized in terms of step responses, and the perturbing noise source N_{ABP} is shown in terms of its power spectrum. A comparison between Figs. 7 and 9 indicates that the indirect identification step responses are much more reliably estimated than the direct identification step responses. Note, in particular, the correspondence between the asymptotic values of the estimated and gold standard step responses in Fig. 9. Also, note the excellent overall agreement between the estimated and gold standard CO→ABP step responses.

Because of the reliable estimation of the asymptotic step-response values shown in Fig. 9, we tested the sensitivity of the indirect identification algorithm in detecting actual changes to the TPR baroreflex static gain values of the cardiovascular model. Figure10 shows the sensitivity results in which the estimated static gains (mean compared with 95% confidence intervals) of the arterial TPR baroreflex and cardiopulmonary TPR baroreflex (as computed from the asymptotic values of the estimated CO→ABP and SV→ABP step responses; see identification algorithms, *Indirect Identification*) are plotted against their respective gold standard static gain values. These results indicate that the indirect identification algorithm can reliably estimate the static gain of the arterial TPR baroreflex, but it somewhat underestimates the static gain of the cardiopulmonary TPR baroreflex. Importantly, however, the sensitivity results show that the indirect identification algorithm is able to track changes in the static gains of both the arterial TPR baroreflex and cardiopulmonary TPR baroreflex, which can be as small as 15–20%.

Figure 11 illustrates the results (analogous to Fig. 9) obtained by applying the indirect identification algorithm to the preliminary set of experimental data in which each subject was in the supine posture. Note that these results also have structural similarities to the corresponding estimates in Fig. 9, which were obtained from the cardiovascular model. The structural similarities that we are referring to here are the fast, initial upstroke followed by damping, positive asymptotic value, and time to reach the asymptotic value for the CO→ABP step response and fast, initial change, negative asymptotic value, and time to reach the asymptotic value for the SV→ABP step response. A virtually identical structural similarity was also seen when the indirect identification algorithm was applied to the experimental data obtained from the human subjects while they were tilted 30° with respect to the supine posture. This suggests that the value of the SNR of TPR fluctuations chosen for the cardiovascular model is not unreasonable (see*Description*).

The average static gain values for the supine arterial TPR baroreflex and cardiopulmonary TPR baroreflex, as computed from the asymptotic values of the step-response estimates shown in Fig. 11, were −0.54 ± 0.44 and −0.75 ± 0.51, respectively. The duration of time for the CO→ABP step response shown in Fig. 11 to reach its asymptotic value was ∼30 s. The average static gain values for the 30° tilt arterial TPR baroreflex and cardiopulmonary TPR baroreflex were −2.28 ± 0.88 and −1.10 ± 0.26, respectively, whereas the duration of time for the 30° tilt CO→ABP step response to reach its asymptotic value was also ∼30 s. On the basis of a paired*t*-test, we did not find statistically significant differences between the supine and 30° tilt gain values for the arterial TPR baroreflex (*P* = 0.16) and cardiopulmonary TPR baroreflex (*P* = 0.41).

For comparison, in the cardiovascular model (which is based on independent experimental findings), the gold standard static gain values for the arterial TPR baroreflex and cardiopulmonary TPR baroreflex were −1.09 and −0.56, respectively, and the time it takes for the gold standard CO→ABP step response to reach steady state was ∼40 s. For an additional comparison, in the study of Raymundo et al. (20) in conscious dogs, the average static gain values for the arterial TPR baroreflex and cardiopulmonary TPR baroreflex (normalized accordingly) were −0.69 and −0.16, respectively.

## DISCUSSION

### Direct Identification

The direct identification algorithm we proposed in this paper is probably the most straightforward approach for attempting to quantitate the TPR baroreflex from beat-to-beat measurements of CO and ABP. However, when we applied this algorithm to the cardiovascular model, we found that reliable TPR baroreflex identification was not achieved (see Fig. 7). The intent of the direct identification algorithm was to analyze the relative fluctuations in ABP, SV, and
with respect to their mean values to characterize the coupling mechanisms in Figs. 2 and 3. We chose to estimate TPR fluctuations with a straightforward method that is robust to the high-frequency wave reflections in noninvasively measured ABP waveforms (see identification algorithms,
*Direct Identification*). However, it turns out that this method artifactually imposes a nonphysiological relationship between the direct identification inputs (ABP and SV) and the direct identification output (
) as follows
Equation 5where each fractional quantity here reflects relative fluctuations with respect to mean values. That is, this relationship erroneously suggests that the arterial TPR baroreflex and SV→TPR step responses are unit step functions scaled by 1 and −1, respectively.

The direct identification problem may therefore be thought of as having a physiological solution and a nonphysiological solution, and one may conclude that this problem is not invertible. However, this is not the case here due to the inclusion of measurement noise, which tends to favor the nonphysiological solution, as Fig. 7 suggests. In fact, as the size of the measurement noise is increased, the estimated asymptotic step-response values of arterial TPR baroreflex and SV→TPR tend toward 1 and −1, respectively. This tendency is less marked for SV→TPR, because HR fluctuations are not considered as an input by the direct identification algorithm. However, as HR variability is reduced, the tendency of the estimated asymptotic step-response value of SV→TPR toward −1 is enhanced.

It is interesting to consider the case in which measurement noise is not present, which is realizable with the cardiovascular model but not a realistic scenario in practice. In this case, we found that reliable estimation of the static gains of the arterial TPR baroreflex and SV→TPR is achieved provided that the model order in *Eq. 1
*is sufficiently small. This implies, as expected, that the estimated TPR fluctuations are accurately determined at very low frequencies. However, we observed that the direct identification problem becomes ill conditioned with a modest increase in the model order. This observation may be explained as follows. *Equation 5
* is not strictly an adequate second solution here due to the presence of HR fluctuations, which are not analyzed by the direct identification algorithm. However, when the model order is even modestly large, there are enough past values of ABP and SV to account for these HR fluctuations, thereby rendering *Eq. 5
* to be an adequate second solution.

When we applied the direct identification algorithm to the preliminary set of experimental human data, we found the estimates to be similar in structure to the corresponding estimates obtained from the cardiovascular model (see results, *Direct Identification*). Without the cardiovascular model-based analysis, one may have immediately interpreted the estimated coupling mechanism between experimental ABP fluctuations and experimental T̂PR fluctuations to represent the positive feedback autoregulation of local vascular beds rather than the negative feedback arterial TPR baroreflex. However, because of the cardiovascular model-based analysis, we believe that the direct identification step-response estimates obtained from the human data are artifactual as a result of the particular method employed for estimating TPR fluctuations. Moreover, because the asymptotic values of these step-response estimates are closer to 1 and −1 (for both the supine and 30° tilt data) than the corresponding estimates obtained from the cardiovascular model, we conclude that the preliminary set of experimental data is corrupted by more measurement noise and/or has less HR variability than the cardiovascular model-generated data.

On the other hand, if TPR fluctuations could be estimated with a different technique, then reliable estimation by the direct identification approach may be possible. Such a technique must be able to overcome the high-frequency wave reflections that corrupt noninvasively measured peripheral ABP waveforms. However, most currently available techniques are not applicable to peripheral ABP waveforms, because they are based on lumped Windkessel theory (e.g., Ref. 4).

### Indirect Identification

The indirect identification algorithm that we proposed in this paper for quantitating the TPR baroreflex assumes prior physiological models (see Fig. 5) but is not encumbered by the challenging task of estimating TPR fluctuations. When we applied this algorithm to the cardiovascular model, we found that the static gain of the arterial TPR baroreflex was reliably estimated. One could utilize this estimated, unitless static gain value to determine, for example, the steady-state, percent TPR change with respect to its mean value that would occur if the arterial TPR baroreflex were stimulated by an *X*% step increase in ABP with respect to its mean value. Although the unitless static gain of the cardiopulmonary TPR baroreflex was consistently underestimated, we found that the indirect identification algorithm reliably detected changes in its value.

Moreover, the indirect identification algorithm was able to estimate accurately the entire CO→ABP step response. If the dynamic properties of the systemic arterial tree are known, then one could compute the step response characterizing the arterial TPR baroreflex from the reliably estimated CO→ABP step response according to the physiological model shown in Fig. 5
*A*. The dynamic properties of the systemic arterial tree may possibly be determined by another technique (see, e.g., Ref. 4) or estimated from the initial portion of the CO→ABP impulse response (assuming the arterial TPR baroreflex is significantly more sluggish than the systemic arterial tree; seeidentification algorithms, *Direct Identification*). Note that with this assumption, one may conclude that the systemic arterial tree time constant of the human subjects in the supine posture is larger than that of the human cardiovascular model (compare Figs. 9 and 11).

However, the entire SV→ABP step response and, as a consequence, the power spectrum of N_{ABP} were not reliably identified, with the step response being inaccurate specifically in terms of transient dynamics. We hypothesize that this is due to the no-delay property of the closed-loop relationship between SV fluctuations and ABP fluctuations. That is, SV fluctuations influence ABP fluctuations through the compliance properties of large arteries, whereas ABP fluctuations simultaneously influence SV fluctuations through afterload effects. Wellstead and Edmunds (24) previously proved that reliable identification is not possible when the data are obtained in closed-loop with no delay. However, importantly, this simultaneous interaction between ABP and SV fluctuations is an immediate, high-frequency effect. Consequently, the static gain of the SV→ABP estimate is relatively unaffected. We also acknowledge that some of the discrepancy between the SV→ABP estimate and its gold standard (shown in Fig. 9) may be attributed to the fact that SV fluctuations are not perfectly determined by RATP fluctuations. This may explain, at least to some extent, why the asymptotic value of the estimated step response is somewhat underestimated. We finally note that the indirect identification estimates were only slightly sensitive to the SD of the measurement noise.

When we applied the indirect identification algorithm to the preliminary set of experimental human data, we also found structural similarities between the experimental estimates and the corresponding estimates obtained from the cardiovascular model (seeresults, *Indirect Identification*). The only major structural difference between these step responses is the sign of the initial change in SV→ABP. However, in both experimental and simulated data, the initial change is fast. We believe that the fast, initial change in the experimental SV→ABP step response is also an estimation artifact due to the no-delay, closed-loop relationship between SV and ABP. We conclude that the cardiovascular model does not behave well in terms of simulating the initial change in ABP due to a change in SV; however, because the initial change is fast for both the experimentally derived and model-estimated SV→ABP step responses, directional accuracy may not be so important for the present study, which focuses on static gains. Also, the difference between the experimentally derived and model-estimated N_{ABP} power spectra is mainly due to this difference between the experimentally derived and model-estimated SV→ABP step response.

Aside from this one difference, the experimentally derived step responses and the step responses identified from the cardiovascular model as computed by both direct and indirect identification algorithms generally appeared to be similar in structure. Quantitative deviations between the experimentally derived and model-estimated step responses may be due to differing parameter values characterizing the structure. That is, it is possible to tune the parameter values of the cardiovascular model to better match the experimental step responses. For example, an increase in the arterial compliance of the model would lower the peak value of the initial upstroke of the CO→ABP step response, which would better match its experimental counterpart. The general structural similarity between the experimentally derived and model-estimated results is important, because it suggests that the cardiovascular model is, to a large extent, reasonable in terms of being able to simulate the hemodynamic behaviors relevant to the study. Because the theoretical evaluation employed here is based on the cardiovascular model, the results obtained from this evaluation are, in turn, more meaningful.

Although we found no statistically significant differences between the supine and 30° tilt static gain values of the arterial TPR baroreflex and cardiopulmonary TPR baroreflex, the trend was in the expected direction with the magnitude of each TPR baroreflex static gain increasing with the tilt (see results: *Indirect Identification*). The lack of statistical significance may be due to the small perturbation imposed. These encouraging results underscore the need for thorough future experimental validation of the indirect identification algorithm, which we are indeed planning to do.

Finally, it is important to note that the indirect identification algorithm requires CO fluctuations and SV fluctuations to be not perfectly linearly correlated. This implies that the indirect identification algorithm will be unreliable when HR variability is negligible. To overcome this limitation, one may be tempted to consider surrogates other than SV fluctuations to represent RATP fluctuations. Of the limited practical possibilities, one potential surrogate is ILV fluctuations, which may be strongly related to RATP fluctuations through ITP variations. We therefore applied the indirect identification algorithm to the cardiovascular model with ILV fluctuations as a surrogate for RATP fluctuations. However, we found that reliable static gain estimation was not achieved due to the low ILV spectral power, near 0 Hz (1, 16).

### Conclusions

In this paper, we proposed a direct identification algorithm and an indirect identification algorithm for quantitating the TPR baroreflex, which is possibly the most important short-term contributor to blood pressure regulation. The two algorithms require only beat-to-beat measurements of CO and ABP, which may both be obtained noninvasively in humans. The direct identification approach is the most obvious and involves estimating TPR fluctuations from the measured signals. The indirect identification approach does not require the estimation of TPR fluctuations but is instead based on assumed physiological models of the relationship between the measured signals. For a theoretical evaluation, we applied the two identification algorithms to beat-to-beat variability generated by a realistic, human cardiovascular model whose dynamic properties are precisely known or easily determined. The results we obtained from the two algorithms were contrasting, with the direct identification algorithm erroneously suggesting mechanisms that were not implemented in the cardiovascular model and the indirect identification algorithm reliably indicating only the implemented physiological mechanisms. Importantly, the indirect identification algorithm was able to track actual changes in the static gains of both the arterial and cardiopulmonary TPR baroreflex mechanisms. We then applied each of these algorithms to a preliminary set of experimental data and obtained contrasting results much like those obtained from the cardiovascular model. The general structural similarity between the model and experimental results makes the results of the theoretical evaluation more meaningful. Moreover, if it were not for the model-based analysis, it would have been difficult to reconcile the physiologically contrasting experimental results. We also found evidence (albeit not statistically significant) that the indirect identification algorithm may be able to detect average changes in the static gains of both the arterial and cardiopulmonary TPR baroreflex mechanisms caused by small perturbations of 30° tilt in posture. This study suggests that, with experimental testing, the indirect identification algorithm may potentially provide the first, noninvasive means for studying the TPR baroreflex in humans during health, disease (e.g., peripheral autonomic neuropathies), and under various environmental factors (e.g., microgravity). This study also provides an example of the role that models can play in the development and initial evaluation of algorithms aimed at quantitating important physiological mechanisms.

## Acknowledgments

This work was sponsored by United States National Aeronautics and Space Administration Grant NAG5-4989 and by grants from the National Space Biomedical Research Institute.

## Footnotes

Address for reprint requests and other correspondence: R. Mukkamala, Michigan State Univ., Dept. of Electrical and Computer Engineering, 2120 Engineering Bldg., East Lansing, MI 48824 (E-mail: rama{at}egr.msu.edu).

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.First published November 14, 2002;10.1152/ajpheart.00532.2002

- Copyright © 2003 the American Physiological Society