AJP - Heart Ad Instruments
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Am J Physiol Heart Circ Physiol 294: H293-H301, 2008. First published November 2, 2007; doi:10.1152/ajpheart.00852.2007
0363-6135/08 $8.00
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow All Versions of this Article:
294/1/H293    most recent
00852.2007v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Chen, X.
Right arrow Articles by Mukkamala, R.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Chen, X.
Right arrow Articles by Mukkamala, R.

Estimation of the total peripheral resistance baroreflex impulse response from spontaneous hemodynamic variability

Xiaoxiao Chen,1 Jong-Kyung Kim,2 Javier A. Sala-Mercado,2 Robert L. Hammond,2,3 Rafat I. Elahi,1 Tadeusz J. Scislo,2 Gokul Swamy,1 Donal S. O'Leary,2 and Ramakrishna Mukkamala1

1Department of Electrical and Computer Engineering, Michigan State University, East Lansing; and 2Departments of Physiology and 3Surgery, Wayne State University School of Medicine, Detroit, Michigan

Submitted 21 July 2007 ; accepted in final form 29 October 2007


   ABSTRACT
TOP
 ABSTRACT
MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX A
 APPENDIX B
 GRANTS
 REFERENCES
 
We previously developed a mathematical analysis technique for estimating the static gain values of the arterial total peripheral resistance (TPR) baroreflex (GA) and the cardiopulmonary TPR baroreflex (GC) from small, spontaneous beat-to-beat fluctuations in arterial blood pressure, cardiac output, and stroke volume. Here, we extended the mathematical analysis so as to also estimate the entire arterial TPR baroreflex impulse response [hA(t)] as well as the lumped arterial compliance (AC). The extended technique may therefore provide a linear dynamic characterization of TPR baroreflex systems during normal physiological conditions from potentially noninvasive measurements. We theoretically evaluated the technique with respect to realistic spontaneous hemodynamic variability generated by a cardiovascular simulator with known system properties. Our results showed that the technique reliably estimated hA(t) [error = 30.2 ± 2.6% for the square root of energy (EA), 19.7 ± 1.6% for absolute peak amplitude (PA), 37.3 ± 2.5% for GA, and 33.1 ± 4.9% for the overall time constant] and AC (error = 17.6 ± 4.2%) under various simulator parameter values and reliably tracked changes in GC. We also experimentally evaluated the technique with respect to spontaneous hemodynamic variability measured from seven conscious dogs before and after chronic arterial baroreceptor denervation. Our results showed that the technique correctly predicted the abolishment of hA(t) [EA = 1.0 ± 0.2 to 0.3 ± 0.1, PA = 0.3 ± 0.1 to 0.1 ± 0.0 s–1, and GA = –2.1 ± 0.6 to 0.3 ± 0.2 (P < 0.05)] and the enhancement of GC [–0.7 ± 0.44 to –1.8 ± 0.2 (P < 0.05)] following the chronic intervention. Moreover, the technique yielded estimates whose values were consistent with those reported with more invasive and/or experimentally difficult methods.

autonomic nervous system; hemodynamics; modeling; system identification; transfer function

THE TOTAL PERIPHERAL RESISTANCE (TPR) baroreflex is a well-appreciated feedback control mechanism for rapid arterial blood pressure (ABP) regulation. Consistent with negative feedback, the arterial TPR baroreflex is widely known to respond to a step increase (decrease) in ABP by decreasing (increasing) steady-state TPR. Similarly, the cardiopulmonary TPR baroreflex has been shown to respond to a step increase (decrease) in central venous pressure (CVP) by decreasing (increasing) steady-state TPR (24). In this way, the cardiopulmonary TPR baroreflex is able to anticipate and thereby quickly buffer the imminent change in ABP that will arise from the preload alteration.

Much of our knowledge of the TPR baroreflex systems is owed to open-loop studies in which the steady-state or static gain properties of one baroreflex was experimentally determined by selectively exciting it and measuring the steady-state TPR response via Ohm's law [i.e., Formula divided by mean cardiac output (Formula)] (see, e.g., Refs. 19 and 27). Indeed, we found only a handful of studies aiming to individually quantify each of the TPR baroreflex systems while they are simultaneously active, presumably because of the need for complex experimental methods and data analysis. In the first of these studies, conducted nearly two decades ago, Raymundo et al. (24) determined individual static gain values of both TPR baroreflex systems through multiple regression analysis of steady-state ABP, CVP, and TPR measurements obtained during a set of step perturbations to the ventricular pacing rate and total blood volume. Although this method is fundamentally sound, it is complex and has yet to be repeated. Furthermore, we found few studies related to quantification of the dynamic properties of the TPR baroreflex (e.g., the time course TPR takes to reach its steady-state value in response to step excitation), perhaps due to the inability to measure sufficiently rapid TPR changes in which Ohm's law may no longer be valid as well as the enhanced demands on experimentation and/or data analysis. In earlier studies, Rosenbaum et al. and other investigators (5, 25, 26) experimentally maintained flow or pressure to a particular regional circulation so as to measure peripheral resistance changes through pressure or flow at that region and then quantified the dynamic properties of one baroreflex by applying step excitation, sinusoidal stimulations at various frequencies, or randomized perturbation in conjunction with system identification analysis of the measured variations. However, this approach is obviously difficult and not as relevant to overall ABP regulation as those studies that have observed TPR. In very recent studies, Aljuri et al. (1) and Hughson et al. (9) estimated TPR changes from beat-to-beat measurements of ABP, CVP, and CO and then sought to identify, for the very first time, the individual linear dynamic properties of both simultaneously active TPR baroreflex systems. In particular, Aljuri et al. essentially generalized the approach of Raymundo et al. by applying dual-input system identification analysis to fluctuations in ABP, CVP, and estimated TPR obtained during randomized pacing and venous occlusion, whereas Hughson et al. applied a similar analysis to these fluctuations as obtained during randomized lower body negative pressure. However, these methods are also difficult, and, as we have previously shown in Ref. 17, estimation of TPR fluctuations could introduce artifacts in the identified TPR baroreflex. Thus, a more convenient technique is needed to fully elucidate the integrated, dynamic functioning of the arterial and cardiopulmonary TPR baroreflex systems.

To this end, in a previous study (17), we introduced a technique for estimating static gain values of the arterial TPR baroreflex (GA) and cardiopulmonary TPR baroreflex (GC) systems by mathematical analysis of beat-to-beat fluctuations in ABP, CO, and stroke volume (SV). Rather than estimating TPR fluctuations, the technique accounts for these variations by physiological modeling of identified systems relating the measured fluctuations. In that study (17), we theoretically validated the technique with respect to realistic hemodynamic variability generated by a cardiovascular simulator with exactly known system properties. In a followup study (15), we experimentally validated the technique by applying it to spontaneous hemodynamic variability measured from seven conscious dogs before and after chronic arterial baroreceptor denervation and showed that the estimated TPR baroreflex static gain values changed in the expected manner following the chronic intervention. In the present study, we extended the mathematical analysis so as to estimate the TPR baroreflex impulse response (i.e., the time derivative of the step response) as well as the lumped arterial compliance (AC). In contrast to previous methods, the extended technique may therefore provide a linear dynamic characterization of TPR baroreflex systems during normal physiological conditions using potentially noninvasive transducers (e.g., Doppler ultrasound and applanation tonometry). In addition, we theoretically and experimentally evaluated this technique with respect to spontaneous hemodynamic variability obtained from the aforementioned cardiovascular simulator and seven conscious dogs.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 MATERIALS AND METHODS
RESULTS
 DISCUSSION
 APPENDIX A
 APPENDIX B
 GRANTS
 REFERENCES
 
Extended Mathematical Analysis Technique

Our initial mathematical analysis technique for estimating GA and GC from beat-to-beat fluctuations in ABP, CO, and SV has been described in detail, including justification of its underlying assumptions, in Refs. 15 and 17. We briefly review this technique below and then describe an extension of the mathematical analysis so as to also estimate the TPR baroreflex impulse response and AC.

The initial technique is based on the compelling study of Raymundo et al. (24). These investigators defined the arterial TPR baroreflex as the system that couples fluctuations in ABP to TPR and the cardiopulmonary TPR baroreflex as the system that couples fluctuations in central venous transmural pressure (CVTP) to TPR. (Note that fluctuations in CVTP are nearly equal to fluctuations in right atrial transmural pressure, which may represent the more precise input to the cardiopulmonary baroreflex.) They demonstrated that these TPR baroreflex systems may be regarded as linear and time invariant (LTI) over a wider hemodynamic range than that attained by the resting variability considered herein. To estimate GA and GC as just defined from only beat-to-beat measurements of ABP, CO, and SV, the initial technique first identifies impulse responses coupling the beat-to-beat measurements (step 1) and then computes static gain values based on physiological models of the identified impulse responses (step 2). Importantly, prior to these two steps are executed, each beat-to-beat measurement (X) has its mean value (X) removed ({Delta}X = X – X) and is then divided by X so as to arrive at normalized, beat-to-beat fluctuations ({Delta}X/X).

More specifically, in step 1 of the technique, the impulse responses coupling {Delta}CO/Formula to {Delta} ABP/Formula (CO->ABP) and {Delta} SV/Formula to {Delta} ABP/Formula (SV->ABP) are identified, as indicated in the block diagram of Fig. 1. In the process, the noise term NABP in the block diagram, which reflects the residual variability in {Delta}ABP/Formula not accounted for by the two inputs, is also estimated. This step is mathematically achieved with the following dual-input, autoregressive exogenous input (ARX) model:

Formula 1(1)
where n and square brackets indicate discrete time. The three sets of unknown coefficients {ai, bi, ci} completely specify the CO->ABP and SV->ABP impulse responses, whereas the unmeasured residual error (WABP) together with the set of coefficients {ai} fully define NABP (12). The unknown model order parameters (m, p, and q) limit the number of coefficients retained in the ARX model. The coefficient sets are estimated from 5- to 10-min intervals of {Delta}ABP/Formula 1, {Delta}CO/Formula 1, and {Delta}SV/Formula 1 sampled to 0.5 Hz by least-squares minimization of WABP in conjunction with a model order reduction algorithm (21). (Typically, this algorithm selects each model order parameter to be between 3 and 5.) Note that an implicit assumption here is that the two inputs are uncorrelated enough for reliable identification. This assumption will be violated when heart rate (HR) variability is absent or tightly correlated to CVTP fluctuations [and thus SV fluctuations (15)] via, for example, a potent Bainbridge reflex.


Figure 1
View larger version (7K):
[in this window]
[in a new window]

  		Fig. 1. Step 1 of the mathematical analysis technique. Spontaneous, beat-to-beat measurements in cardiac output (CO), stroke volume (SV), and arterial blood pressure (ABP) are analyzed to estimate the impulse responses of the two systems, CO->ABP and SV->ABP, as well as the noise term NABP, which represents the residual ABP variability not accounted for by the two systems. For each of the three beat-to-beat measurements X, the notation {Delta}X/Formula 4 is precisely defined as (X – Formula 4)/Formula 4, where Formula 4 denotes its mean value. hCO->ABP(t), CO->ABP impulse response.

 
In step 2 of the technique, GA and GC are computed by representing the identified CO->ABP and SV->ABP impulse responses with the physiologic models shown in Fig. 2. In these models, the systemic arterial tree (SAT) is defined to be the system that couples {Delta}CO/Formula 1 and {Delta}TPR/Formula 1 to {Delta}ABP/Formula 1, while the inverse heart-lung unit is defined to be the system that couples {Delta}SV/Formula 1 to {Delta}CVTP/Formula 1. A key point is that each of these two systems, which are likewise assumed to be LTI, has a static gain value of unity due to its input and output reflecting normalized fluctuations. (For example, it is easy to show this property for the SAT through Ohm's law.) As a result, GA and GC may be specifically computed from the static gain values of the identified CO->ABP and SV->ABP impulse responses based on the physiological models shown in Fig. 2. This step is achieved with the following formulas involving the ARX coefficient sets estimated from step 1:

Formula 2(2)
See APPENDIX A for a derivation of these two formulas.


Figure 2
View larger version (18K):
[in this window]
[in a new window]

  		Fig. 2. Step 2 of the mathematical analysis technique. The static gain values of the arterial and cardiopulmonary total peripheral resistance (TPR) baroreflex systems (GA and GC, respectively) are computed from the static gain values of the CO->ABP and SV->ABP impulse responses estimated in step 1 according to the depicted physiological models, wherein the systemic arterial tree (SAT) and inverse heart-lung unit have unity static gain values. A: physiological model used to derive the GA formula; B: physiological model used to derive the GC formula. CVTP, central venous transmural pressure; hSAT(t), SAT impulse response; hA(t), arterial TPR baroreflex impulse response.

 
To extend the technique so as to estimate the arterial TPR baroreflex impulse response [hA(t)] (step 3), as shown in Fig. 2A, the SAT impulse response (hSAT) must likewise be known. However, while the static gain value of the SAT is known to be unity, its system dynamics are generally unknown and must first be estimated. To do so, the following two assumptions are made based on known physiology. First, as justified in Refs. 16 and 18, the SAT is assumed to be well represented by a lumped parameter model governed by a single time constant ({tau}SAT) equal to the product of TPR and AC for the slow, beat-to-beat fluctuations considered herein. In other words, the sampled hSAT (hSAT[n]) is assumed to take on the following form:

Formula 3(3)
where t and round brackets indicate continuous time; u is the unit step function; and Fs is the sampling frequency, likewise set to 0.5 Hz. Note that the integrand here represents the continuous-time hSAT (with unity static gain) that arises from the above assumption, whereas the integration preserves unity static gain after sampling (by virtue of the sum of hSAT[n] being 1 for all physical {tau}SAT) as well as attenuates any artifacts in the sampling process. Second, as justified from Refs. 2–4 and 25, TPR baroreflex dynamics are assumed to be delayed and slower than SAT dynamics. Thus, as shown in Fig. 3, the initial (<5 s) time evolution of the CO->ABP impulse response (which represents the {Delta}ABP/Formula 3 response to a transient unit increase in {Delta}CO/Formula 3) may be attributed solely to the SAT, as the more sluggish arterial TPR baroreflex has yet to be activated. So, in step 3 of the extended technique, {tau}SAT in Eq. 3 is first solved for by fitting hSAT[0] to the corresponding sample of the identified 0.5-Hz CO->ABP impulse response (hCO->ABP[n]). In this way, hSAT[n] is fully defined, and AC may be determined via {tau}SAT/(Formula 3/Formula 3). Then, hA[n] is calculated from hSAT[n] and hCO->ABP[n] through the following formula arising from the physiological model shown in Fig. 2A:

Formula 4(4)
where {otimes} is the convolution operation governing the input-output relationship of LTI systems and the indicated inverse is defined as h[n] {otimes} h[n]–1 = {delta}[n] (where {delta} represents the unit impulse function). While hA[n] here may be determined through standard computation of the involved inverse or exactly solved via an ARX formulation, any error or noise in the estimates of hSAT[n] and hCO->ABP[n] would be greatly amplified by either of these naïve inversion calculations. Thus, hA[n] in Eq. 4 is computed through a noise robust method that we developed in which the impulse response is compactly represented with damped sinusoidal basis functions whose parameters are estimated using regularized least squares with the constraint that its static gain value is equal to GA from Eq. 2. Finally, estimated hA[n] is converted to physical units via hA(t) = FshA[tFs] (20). See APPENDIX B for a complete mathematical specification of the extension.


Figure 3
View larger version (9K):
[in this window]
[in a new window]

  		Fig. 3. Step 3 (extension) of the mathematical analysis technique. hSAT(t) in the physiologic models shown in Fig. 2 (dashed line) is estimated from the initial portion of hCO->ABP(t) estimated in step 1 (solid line) in which the TPR baroreflex has yet to be activated. The lumped arterial compliance is then specified by hSAT(t) and measured Formula 4, and hA(t) is estimated from hSAT(t) and hCO->ABP(t) according to the physiologic model shown in Fig. 2A.

 
Note that the inverse heart-lung unit impulse response is neither known nor estimated. The extended technique therefore does not provide a direct estimate of the cardiopulmonary TPR baroreflex impulse response. However, since both TPR baroreflex systems are governed by the {alpha}-sympathetic nervous system (7), it may be reasonable to assume that their dynamic properties are quite similar. To the extent that this assumption holds, an estimate of the cardiopulmonary TPR baroreflex impulse response may then be obtained by scaling estimated hA(t) to have a static gain value equal to GC from Eq. 2.

Technique Evaluation

Theoretical evaluation. We theoretically evaluated the extended mathematical analysis technique based on a cardiovascular simulator, which, in contrast to an experimental system, offers precisely known TPR baroreflex system properties and AC values for direct comparison of the estimates. We specifically employed a cardiovascular simulator that we previously developed, demonstrated to generate realistic short-term, beat-to-beat variability, and utilized to theoretically validate our initial technique. The simulator is described in detail in Refs. 14 and 17. Briefly, the major components of the simulator are a circulatory system, a short-term regulatory system, and resting perturbations. The circulatory system consists of contracting left and right ventricles, systemic and pulmonary arterial trees, and systemic and pulmonary veins. Each of these six compartments is specifically modeled as a first-order system accounting for both viscous and compliant effects. The regulatory system comprises arterial and cardiopulmonary baroreflex control of HR, TPR, systemic venous unstressed volume, and maximum ventricular elastance as well as a direct neural coupling between respiration and HR. Each of the baroreflex effector systems is specifically modeled as a static nonlinearity to account for saturation phenomenon followed by LTI dynamics. The resting perturbations include spontaneous respiratory activity, which impinges on the circulatory system through a first-order model of the airways and lung as well as the direct neural coupling, low-frequency disturbances to TPR and systemic venous unstressed volume, and a 1/(f, frequency) disturbance to HR. These three disturbances are randomly generated. Thus, each simulator run or realization will actually differ under the same set of simulator parameter values. For the present study, we replaced the first-order model of the SAT with a well-known third-order model accounting for inertial effects in addition to viscous and compliant effects (6, 13). (This modification, for example, has the effect of introducing a bump in the diastolic interval of the simulated ABP waveform.) In this way, the SAT of the simulator was more complex than what is assumed by the technique.

For this theoretical evaluation, our specific aim was to determine if the technique could reliably estimate hA(t), GC, and AC when applied to hemodynamic variability simulated under different TPR baroreflex static gain values (implemented via system amplitude scaling), TPR baroreflex overall time constant values (implemented via system time scaling), and AC values. To achieve this aim, we first established the actual reference quantities against which the estimates could be compared for each investigated set of simulator parameter values. More specifically, we determined the actual reference hA(t) by separating the arterial TPR baroreflex model from the simulator, applying an impulse input to the isolated model, and measuring the output response (see Ref. 17 for details); the actual reference GC value by first establishing the actual reference cardiopulmonary TPR baroreflex impulse response in an analogous manner and then computing its sum; and the actual reference AC value by summing the two individual compliances in the SAT model. Then, we applied the technique to 25 different realizations of ABP, CO, and SV simulated under each investigated set of simulator parameter values to determine the mean and SE of the estimates. Next, we succinctly characterized the estimated and actual reference hA(t) in terms of the following parameters: the square root of the energy of hA(t) (EA), absolute peak amplitude of hA(t) (PA), GA [sum of hA(t) values], and overall time constant of hA(t) determined with a robust rectangular-based method (11) ({tau}A). Finally, we quantitatively compared the estimates of these four impulse response parameters as well as GC and AC with their actual reference quantities in terms of the root mean square (RMS) of the error (E) normalized (N) by the actual reference quantity (i.e., RMSNE).

Experimental evaluation. We also experimentally evaluated the extended technique based on canine hemodynamic data, which were originally collected by us to address different specific aims and previously utilized to experimentally validate our initial technique. The materials and methods used in the collection of these data are described in detail in Ref. 10. As we have done in Ref. 15, we briefly describe the most relevant aspects of the experimental procedures below. All procedures were reviewed and approved by the Wayne State University Animal Investigation Committee. We studied seven conscious dogs (20–25 kg) of either gender. For each dog, we installed chronic instrumentation through a series of aseptic surgeries for measurement of continuous ABP (via an aortic catheter, Abbott Laboratories), CO (via an aortic flow probe, Transonic Systems), and HR as well as other hemodynamic variables. After dogs had recovered from the instrumentation surgeries, we recorded the hemodynamic variables on a beat-to-beat basis for ~10 min while the dog was standing quietly. Then, we achieved complete baroreceptor denervation by transection of the aortic depressor and carotid sinus nerves. We verified the completeness of the arterial baroreceptor denervation by observing essentially no HR response to an intravenous bolus infusion of phenylephrine, which increased Formula 4 by ~40 mmHg. Finally, ~2 wk after the completion of the arterial baroreceptor denervation surgeries, we likewise recorded beat-to-beat hemodynamic variables.

For this experimental evaluation, our main specific aim was to determine if the technique could detect the alterations in TPR baroreflex functioning that are known to occur following chronic intervention. To address this aim, we first applied the technique to the beat-to-beat ABP, CO, and SV (= CO/HR) recorded from each dog before and after the chronic intervention. Then, we again parameterized each estimated hA(t) with the EA, PA, GA, and {tau}A. Finally, we performed paired t-tests to determine if the group averages of the impulse response parameters as well as estimated GC and AC and measured Formula 4 were significantly altered by chronic arterial baroreceptor denervation. We considered a P value of <0.05 to be statistically significant.


   RESULTS
TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
DISCUSSION
 APPENDIX A
 APPENDIX B
 GRANTS
 REFERENCES
 
Theoretical Evaluation

Figures 4 and 5, respectively, show the estimated hA(t) (means ± SE) along with the actual reference hA(t) for three different simulator TPR baroreflex static gain values (but with the same overall time constant) and three different simulator TPR baroreflex overall time constants (but with the same static gain value). Visually, the estimated hA(t) closely agreed with the reference hA(t) for each of the investigated simulator parameter values. Quantitatively, this agreement corresponded to overall EA, PA, GA, and {tau}A RMSNEs of 30.2 ± 2.6%, 19.7 ± 1.6%, 37.3 ± 2.5%, and 33.1 ± 4.9%, respectively. Table 1 shows the estimated GC and AC (means ± SE) and their corresponding actual reference values for three different simulator TPR baroreflex static gain values and AC values, respectively. The overall GC RMSNE was 48.4 ± 12.0%, whereas the overall AC RMSNE was17.6 ± 4.2%. Although the estimated GC was the least accurate, it was able to reliably detect changes in the simulator GC value. In other words, the GC RMSNE had a significant bias component.


Figure 4
View larger version (18K):
[in this window]
[in a new window]

  		Fig. 4. Theoretical evaluation results of the extended mathematical analysis technique in which the estimated hA(t) [mean (dashed line) ± SE (dashed-dotted line)] is shown along with the corresponding actual reference hA(t) (solid line) for 3 different simulator TPR baroreflex gain values. The quantitative values at the top are the root mean squared normalized error (RMSNE) of the following 4 parameters of the estimated hA(t): square root of energy (EA), absolute peak amplitude (PA), static gain value or sum (GA), and overall time constant ({tau}A).

 

Figure 5
View larger version (18K):
[in this window]
[in a new window]

  		Fig. 5. Theoretical evaluation results of the extended mathematical analysis technique in which the hA(t) [mean (dashed line) ± SE (dashed-dotted line)] is shown along with the corresponding actual reference hA(t) (solid line) for 3 different simulator TPR baroreflex overall time constant values.

 

View this table:
[in this window]
[in a new window]

  		Table 1. Theoretical evaluation results of the extended mathematical analysis technique in which two of the estimated parameters from the cardiovascular simulator are shown along with the actual reference values and corresponding quantitative errors

 
Experimental Evaluation

Figure 6 shows the group average estimated hA(t) (means ± SE) from the seven conscious dogs before and after chronic arterial baroreceptor denervation. Before the baroreceptor denervation (i.e., the control condition), the group average estimated hA(t) indicated that, in response to a hypothetical impulse increase in ABP at time 0, TPR would immediately decrease and then return toward baseline via a small amplitude oscillation, with {tau}A equal to 7.6 s. Following the baroreceptor denervation, the group average estimated hA(t) was blunted essentially to 0 and therefore indicated that TPR would remain virtually unaltered in response to an ABP change. Table 2 shows the group averages (means ± SE) of three parameters of the estimated hA(t) (EA, PA, and GA), estimated GC and AC, and measured Formula 4 from the dogs before and after chronic arterial baroreceptor denervation. [Note that {tau}A is not included in Table 2 because it not well defined when hA(t) is nearly zero.] Group average estimated EA, PA, and GA were all significantly blunted following baroreceptor denervation, thereby statistically buttressing the visual results shown in Fig. 6. In addition, the magnitude of the group average estimated GC was significantly enhanced after the chronic intervention. [Note that the GA and GC results here were reported earlier in the context of experimentally evaluating our initial technique (15).] On the other hand, the group average estimated AC was reduced but not to a level of statistical significance after the baroreceptor denervation, with a value of 2.69 ± 0.44 ml/mmHg over both conditions. The group average measured Formula 4 was likewise not significantly altered.


Figure 6
View larger version (14K):
[in this window]
[in a new window]

  		Fig. 6. Experimental evaluation results of the extended mathematical analysis technique in which the group average estimated hA(t) [mean (solid line) ± SE (dashed line)] from 7 conscious dogs is shown before and after chronic arterial baroreceptor denervation.

 

View this table:
[in this window]
[in a new window]

  		Table 2. Experimental evaluation results of the extended mathematical analysis technique in which group averages of the estimated parameters from seven conscious dogs are shown before and after chronic arterial baroreceptor denervation along with the level of statistical significance

 

   DISCUSSION
TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
APPENDIX A
 APPENDIX B
 GRANTS
 REFERENCES
 
In two previous studies (15, 17), we introduced and demonstrated the validity of a technique based on system identification analysis and physiological modeling for estimating GA and GC from spontaneous beat-to-beat fluctuations in ABP, CO, and SV (see Figs. 1 and 2). In the present study, we extended this technique through additional physiological modeling and analysis so as to estimate hA(t) as well as AC (see Fig. 3). In addition, with the assumption that the two {alpha}-sympathetically mediated TPR baroreflex systems have the same dynamic properties, the cardiopulmonary TPR baroreflex impulse response may then be estimated by scaling hA(t) to have a static gain value equal to GC. However, we have left investigation of the veracity of this assumption for the future.

The extended technique introduced herein is closely related to a few mathematical methods proposed very recently in the cardiovascular physiology literature. As discussed above, Aljuri et al. (1) and Hughson et al. (9) also proposed to estimate the linear dynamic properties of both simultaneously active TPR baroreflex systems. Potential advantages of the extended technique are that it does not require invasive measurements (i.e., CVP), any external perturbations, or estimation of beat-to-beat fluctuations in TPR, which, as mentioned above, could introduce artifacts in the identified TPR baroreflex impulse response. On the other hand, with a CVP measurement, direct estimation of the cardiopulmonary TPR baroreflex impulse response is, in principle, feasible. In addition, Quick et al. (23) proposed to estimate AC by likewise exploiting slow, beat-to-beat fluctuations in ABP and CO via the calculated impedance of the SAT (23). The extended technique here is different in that it attempts to be more precise by also accounting for baroreflex dynamics.

To evaluate the extended technique, ideally, the estimated hA(t) (with a static gain value of GA), GC, and AC from an experimental system would be directly compared with high-fidelity, independent reference quantities. However, as we have argued in Ref. 14, it would be very difficult to establish such reference quantities. For example, to experimentally determine reference hA(t), beat-to-beat TPR would have to be measured in response to an impulsive change in ABP while CVTP was held constant. On the other hand, a direct comparison could be readily conducted using a cardiovascular simulator with precisely known system properties. However, this theoretical evaluation would only be as meaningful as the extent to which the simulator coincided with an experimental system. Alternatively, an experimental evaluation in terms of detecting changes in the estimated quantities to a known intervention could be more easily achieved. Although this sort of experimental evaluation would not establish the extent to which the estimated quantities are correct, it would at least demonstrate the utility of the technique with respect to real data. Thus, in the present study, we performed theoretical and experimental evaluations of the technique to reap the advantages of both types of evaluations.

We first theoretically evaluated the extended technique with respect to a cardiovascular simulator that we previously developed and demonstrated to generate realistic spontaneous beat-to-beat variability. Our results showed that, under various simulator parameter values, the technique was able to reliably estimate hA(t) and AC (see Figs. 4 and 5 and Table 1), even though the simulator included more complex dynamics than those assumed by the technique (e.g., baroreflex saturation and inertial effects in the SAT). However, the technique did underestimate the magnitude of GC (see Table 1) due to the fact that SV fluctuations did not entirely account for CVTP fluctuations in the simulator, as assumed in the physiological model shown in Fig. 2B (see Ref. 15 for a further discussion of this assumption). Importantly, the underestimation was consistent, and the technique was therefore able to reliably track changes in the simulator GC value. These theoretical evaluation results help validate various aspects of the mathematical extension introduced herein. In particular, the results indicate that 1) the specific sampling method in Eq. 3 is justified; 2) the basis function representation for hA(t) in Eq. B3 is flexible enough under variations to the TPR baroreflex static gain value and overall time constant; and 3) the bias introduced by the use of regularization (see APPENDIX B), which is needed to reduce the error variance due to noise amplification resulting from the inversion in Eq. 4, is modest.

We then experimentally evaluated the extended technique with respect to spontaneous hemodynamic variability measured from seven conscious dogs during a control condition and following chronic arterial baroreceptor denervation. Our results generally showed that the technique was able to correctly predict the changes in TPR baroreflex functioning that are known to occur following the chronic intervention as well as yield estimates whose values are consistent with those previously reported with more invasive and/or experimentally difficult methods.

More specifically, the group average estimated hA(t) from the control condition indicated negative feedback dynamics (see Fig. 6), which is congruent with known baroreflex physiology. Moreover, the overall time constant of this group average impulse response of 7.6 s is consistent with the value of 9 s reported in the study by Rosenbaum et al. referred to above. However, the lack of a delay in the impulse response is not supported by current evidence (2, 3). Poor agreement between estimated hSAT(t) and hCO->ABP(t) at t = 2 s (in which the TPR baroreflex should yet to be activated; see Fig. 3) due to imperfect estimation of the latter impulse response (i.e., step 1 of the technique) certainly contributed to the lack of delay. However, note that this estimation error was partially attenuated by averaging the estimates over all the dogs. Another contributing factor may have been imperfect estimation in step 3 of the technique, as the time delay was also missed during the theoretical evaluation (see Figs. 4 and 5) wherein the agreement between hSAT(t) and hCO->ABP(t) at t = 2 s was excellent.

The control group average estimated hA(t) appears similar to the one determined by Aljuri et al. in conscious sheep in terms of reflecting negative feedback dynamics, overall time constant, and general shape, but different in detail with a smoother initial downstroke and a somewhat larger overshoot. On the other hand, the impulse response reported here is grossly different from the one obtained by Hughson et al. in humans, which reflected positive feedback dynamics that were attributed to a myogenic response. Differences between our control group average estimated hA(t) and the ones previously reported could be entirely due to variations in experimental preparations. As discussed above, it is also possible that imperfect estimation of TPR fluctuations could have rendered their impulse response estimates to be different.

Following chronic arterial baroreceptor denervation, the group average estimated hA(t) along with its parameters (EA, PA, and GA) were largely abolished, as expected (see Fig. 6 and Table 2). [Note that a GA value of 0 does not necessitate that hA(t) be 0 for all t but rather that only its sum over t be 0.] In addition, the magnitude of the group average estimated GC was increased after the chronic intervention, perhaps as a result of a central compensatory mechanism (see Table 2). These changes are entirely consistent with the pattern of alterations in GA and GC reported in the aforementioned study by Raymundo et al. after the same intervention in conscious dogs.

On the other hand, the group average estimated AC did not change following chronic arterial baroreceptor denervation (see Table 2). This lack of change is consistent with the belief that the baroreflex only controls HR, TPR, ventricular contractility, and systemic venous unstressed volume (7). To give further credence to this result, the group average estimated AC over both conditions of 2.69 ± 0.44 ml/mmHg is on par with the previously reported value of ~2 ml/mmHg determined with the flow stop method (4). Finally, the group average measured Formula 4 did not change after the baroreceptor denervation (see Table 2), as baroreflex functioning does not affect mean hemodynamic values (7).

In conclusion, we have introduced a novel technique to mathematically estimate linear dynamic properties of the arterial and cardiopulmonary TPR baroreflex systems as well as AC from spontaneous hemodynamic variability that could potentially be measured noninvasively. We have demonstrated the validity of the technique using both realistic simulated data and experimental measurements. In the future, we plan to continue evaluating the technique in animals as well as begin testing in humans with interventions of known effects. Ultimately, we hope that the practical technique is employed to advance the basic understanding of the normal, integrated, and dynamic functioning of the TPR baroreflex systems in humans and animals in health and disease.


   APPENDIX A
TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX A
APPENDIX B
 GRANTS
 REFERENCES
 
We derive below the formulas for GA and GC in Eq. 2. These formulas generally arise by equating the static gain values of the CO->ABP and SV->ABP impulse responses expressed in terms of the dual-input ARX model of Eq. 1 to the corresponding static gain values expressed with the physiologic models shown in Fig. 2 through the application of basic linear systems theory (11, 20, 22).

To derive the GA formula, the z-transform is first applied to the dual-input ARX model to compute the system function of the CO->ABP impulse response [HCO->ABP(z)] as follows:

Formula A1(A1)
Then, the static gain value of the CO->ABP impulse response (GCO->ABP) is expressed in terms of ARX coefficients by evaluating HCO->ABP(z) at z = 1 as follows:

Formula A2(A2)
Note that the rightmost equality here merely indicates that the static gain value is defined to be the steady-state output response divided by the steady-state input. Next, linear systems theory is applied to the physiological model shown in Fig. 2A to give the following steady-state relationships:

Formula A3(A3)

Formula A4(A4)
where unity static gain has been invoked for the SAT in Eq. A3. Then, Eq. A3 is substituted into Eq. A4 to give the following expression:

Formula A5(A5)
Finally, Eq. A2 is substituted into Eq. A5 so as to arrive at the formula for GA in Eq. 2.

To derive the GC formula, the static gain value of the SV->ABP impulse response (GSV->ABP) is first similarly expressed in terms of ARX coefficients as follows:

Formula A6(A6)
Then, linear systems theory is applied to the physiological model shown in Fig. 2B to give the following steady-state relationships:

Formula A7(A7)

Formula A8(A8)

Formula A9(A9)
where unity static gain has again been invoked for SAT in Eq. A7 as well as for the inverse heart-lung unit in Eq. A9. Next, Eqs. A7 and A9 are substituted into Eq. A8 to give the following expression:

Formula A10(A10)
Finally, the GA formula in Eqs. 2 and A6 is substituted into Eq. A10 so as to arrive at the formula for GC in Eq. 2.


   APPENDIX B
TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX A
 APPENDIX B
GRANTS
 REFERENCES
 
We outline below a complete mathematical description of step 3 of the extended technique. See Refs. 8, 22, and 28 for related background material.

First, hSAT[n] is estimated from the identified hCO->ABP[n] based on physiological knowledge (see Extended Mathematical Analysis Technique). In particular, hSAT[n] in Eq. 3 is determined by numerically searching for the value of {tau}SAT over a physiological range (0.1–10 s) that equates hSAT[0] to hCO->ABP[0]. Then, hSAT[n] in Eq. 3 is computed from n = 0 to n = N – 1, where N represents its effective length (i.e., the number of samples or time duration for which the impulse response is essentially nonzero).

Second, hA[n] is computed from the newly estimated hSAT[n] and identified hCO->ABP[n] according to the physiological model shown in Fig. 2A. More specifically, hCO->ABP[n], which is strictly of infinite length due to the ARX representation in Eq. 1, is truncated from n = 0 to n = M – 1, where M likewise represents its effective length. In addition, hA[n] is assumed to be effectively of finite length and nonzero from n = 1 to n = L. [The value of L is set to 50 so that hA(t) can extend up to 100 s.] Now, invoking the definitions of the inverse and convolution, Eq. 4 may be rewritten as follows:

Formula 1(B1)
where g[n] = hSAT[n] {otimes} hCO->ABP[n] and e[n] has been introduced to account for any estimation and/or modeling error. This equation represents a linear set of equations with N + M + L – 2 equations corresponding to the discrete times n isin [1, N + M + L – 2] for which the equation is nonzero and L unknowns corresponding to the samples of hA[n]. Since the static gain value of hA[n] has already been estimated in step 2 to be GA in Eq. 2, the following additional equation arises:

Formula 2(B2)
Taken together, Eqs. B1 and B2 represent a linear set of equations with N + M + L – 1 equations and L unknowns. This set of equations may be conveniently solved based on least-squares minimization of e[n] via the closed-form linear least-squares solution with Tikonov regularization to attenuate any noise amplification in the involved inverse calculation. However, since L is relatively large, the resulting estimates of the unknowns may still appear to be noisy.

To reduce the number of unknowns so as to improve the estimation accuracy, hA[n] is further assumed to be compactly represented with damped sinusoidal basis functions as follows:

Formula 3(B3)
where {{lambda}k} is a set of unknown damping parameters, {dk, fk} is a set of unknown coefficient parameters, {{omega}k} is a set of unknown frequency parameters, and r is an unknown number of basis functions. This assumption is similar to representing hA[n] with an ARX equation, which was previously done by Aljuri et al. (1) and Hughson et al. (9) and generally permits a significant reduction in the number of parameters to be estimated. Substituting Eq. B3 into Eqs. B1 and B2 yields the following equations:

Formula 4(B4)
These equations represent a nonlinear set of equations with N + M + L – 1 equations and 4r unknowns corresponding to the basis function parameters. This set of equations may also be estimated by regularized least-squares minimization of e[n]. However, this optimization problem is considerably more difficult, because the equations are nonlinear, particularly in the unknown parameter sets {{lambda}k, {omega}k}. To simplify this problem, each damping factor in the parameter set {{lambda}k} is assumed to be equal to {lambda} = exp(–3/L) so that hA[n] approximately decays to 0, while the frequencies in the parameter set {{omega}k} are allowed to take on only discrete values according to the Fourier series {i.e., 2{pi}l/L for l = 0, 1,..., ceil[(L – 1)/2], where ceil(x) is the smallest integer ≥ x}. Then, for each set of frequency parameters {{omega}k} considered (see below), the corresponding coefficient parameter sets are estimated through the linear least-squares solution with Tikonov regularization. More specifically, Eq. B4 may be expressed in matrix form by stacking each individual equation, corresponding to each n, one on top of the other as follows:


Formula 5

(B5)
where

Formula 4
with j = 1, 2 and f[x] = cos[x] for j = 1 and sin[x] for j = 2. Then, after the last row of matrix A and vector b is scaled by a factor of 10 so as to ensure that the estimated hA[n] will have a static gain value of nearly GA, the regularized linear least-squares estimate of vector x comprising the coefficient parameters is obtained via x = [(ATA + {alpha}I)–1]ATb, where T means matrix transpose, {alpha} represents the regularization factor and I is the 2r x 2 r identity matrix. (The value of {alpha} was empirically set to 2, which essentially eliminated spurious noise in hA[n]. However, hA[n] was not very sensitive to this specific value.) Among all such estimates computed for each considered set of frequency parameters {{omega}k}, the one that results in the least-squares value of the vector Ax – b is selected so as to provide the optimal estimate of the parameter set {dk, fk, {omega}k}. Finally, the number of basis functions r is determined iteratively by starting with a single basis function representation and then adding one basis function at a time until the least-squares value of the vector Ax – b no longer decreases by >5%. For further simplicity, in the kth iteration, only the frequency parameter of the newly added basis function is added with the frequency parameters of the previous (k – 1) basis functions set to the estimates obtained from the (k – 1)th iteration. In the kth iteration, the number of sets of frequency parameters considered is therefore {ceil[(L – 1)/2] – k + 2}.


   GRANTS
TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX A
 APPENDIX B
 GRANTS
REFERENCES
 
This work was supported by National Institutes of Health Grants EB-004444 and HL-55473.


   FOOTNOTES
 

Address for reprint requests and other correspondence: R. Mukkamala, Dept. of Electrical and Computer Engineering, Michigan State Univ., 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.


   REFERENCES
TOP
 ABSTRACT
 MATERIALS AND METHODS
 RESULTS
 DISCUSSION
 APPENDIX A
 APPENDIX B
 GRANTS
 REFERENCES
 

  1. Aljuri N, Marini R, Cohen RJ. Test of dynamic closed-loop baroreflex and autoregulatory control of total peripheral resistance in intact and conscious sheep. Am J Physiol Heart Circ Physiol 287: H2274–H2286, 2004.[Abstract/Free Full Text]
  2. Berger RD, Saul JP, Cohen RJ. Transfer function analysis of autonomic regulation. I. Canine atrial rate response. Am J Physiol Heart Circ Physiol 256: H142–H152, 1989.[Abstract/Free Full Text]
  3. Borst C, Karemaker JM. Time delays in the human baroreceptor reflex. J Auton Nerv Syst 9: 399–409, 1983.[CrossRef][Web of Science][Medline]
  4. Brunner MJ, Bishop GG, Shigemi K. Arterial compliance and its control by the baroreflex in hypertensive dogs. Am J Physiol Heart Circ Physiol 265: H616–H620, 1993.[Abstract/Free Full Text]
  5. Dampney RA, Taylor MG, McLachian EM. Reflex effects of stimulation of carotid sinus and aortic baroreceptors on hindlimb vascular resistance in dogs. Circ Res 29: 119–127, 1971.[Abstract/Free Full Text]
  6. Guarini M, Urzúa J, Cipriano A, González W. Estimation of cardiac function from computer analysis of the arterial pressure waveform. IEEE Trans Biomed Eng 45: 1420–1428, 1998.[CrossRef][Web of Science][Medline]
  7. Guyton AC, Hall JE. Textbook of Medical Physiology. Philadelphia, PA: Saunders, 1996.
  8. Hansen PC. Deconvolution and regularization with Toeplitz matrices. Numerical Algorithms 29: 323–378, 2002.[CrossRef][Web of Science]
  9. Hughson RL, O'Leary DD, Shoemaker JK, Lin DC, Topor ZL, Edwards MR, Tulppo MP. Searching for the vascular component of the arterial baroreflex. Cardiovasc Eng 4: 155–162, 2004.[CrossRef]
  10. Kim JK, Sala-Mercado JA, Rodriguez J, Scislo TJ, O'Leary DS. Arterial baroreflex alters strength and mechanisms of muscle metaboreflex during dynamic exercise. Am J Physiol Heart Circ Physiol 288: H1374–H1380, 2005.[Abstract/Free Full Text]
  11. Lathi BP. Signal Processing and Linear Systems. Carmichael, CA: Berkeley Cambridge, 1998.
  12. Ljung L. System Identification: Theory for the User. Englewood Cliffs, NJ: Prentice Hall, 1987.
  13. Manning TS, Shykoff BE, Izzo JL. Validity and reliability of diastolic pulse contour analysis (windkessel model) in humans. Hypertension 39: 963–968, 2002.[Abstract/Free Full Text]
  14. Mukkamala R, Cohen RJ. A forward model-based validation of cardiovascular system identification. Am J Physiol Heart Circ Physiol 281: H2714–H2730, 2001.[Abstract/Free Full Text]
  15. Mukkamala R, Kim J, Li Y, Sala-Mercado J, Hammond RL, Scislo TJ, O'Leary DS. Estimation of arterial and cardiopulmonary total peripheral resistance baroreflex gain values: validation by chronic arterial baroreceptor denervation. Am J Physiol Heart Circ Physiol 290: H1830–H1836, 2006.[Abstract/Free Full Text]
  16. Mukkamala R, Reisner AT, Hojman HM, Mark RG, Cohen RJ. Continuous cardiac output monitoring by peripheral blood pressure waveform analysis. IEEE Trans Biomed Eng 53: 459- 467, 2006.[CrossRef][Web of Science][Medline]
  17. Mukkamala R, Toska K, Cohen RJ. Noninvasive identification of the total peripheral resistance baroreflex. Am J Physiol Heart Circ Physiol 284: H947–H959, 2003.[Abstract/Free Full Text]
  18. Noordergraaf A. Circulatory System Dynamics. New York: Academic, 1978.
  19. Olivier NB, Stephenson RB. Characterization of baroreflex impairment in conscious dogs with pacing-induced heart failure. Am J Physiol Regul Integr Comp Physiol 265: R1132–R1140, 1993.[Abstract/Free Full Text]
  20. Oppenheim AV, Schafer RW, Buck JR. Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice Hall, 1999.
  21. Perrott MH, Cohen RJ. An efficient approach to ARMA modeling of biological systems with multiple inputs and delays. IEEE Trans Biomed Eng 43: 1–14, 1996.[Web of Science][Medline]
  22. Proakis JG, Manolakis DK. Digital Signal Processing: Principles, Algorithms and Application. Upper Saddle River, NJ: Prentice Hall, 1996.
  23. Quick CM, Berger DS, Hettrick DA, Noordergraaf A. True arterial system compliance estimated from apparent arterial compliance. Ann Biomed Eng 28: 291–301, 2000.[CrossRef][Web of Science][Medline]
  24. Raymundo H, Scher AM, O'Leary DS, Sampson PD. Cardiovascular control by arterial and cardiopulmonary baroreceptors in awake dogs with atrioventricular block. Am J Physiol Heart Circ Physiol 257: H2048–H2058, 1989.[Abstract/Free Full Text]
  25. Rosenbaum M, Race D. Frequency-response characteristics of vascular resistance vessels. Am J Physiol 215: 1397–1402, 1968.[Free Full Text]
  26. Scher AM, Young AC. Servoanalysis of carotid sinus reflex effects on peripheral resistance. Circ Res 12: 152–162, 1963.[Free Full Text]
  27. Schmidt RM, Kumada M, Sagawa K. Cardiac output and total peripheral resistance in carotid sinus reflex. Am J Physiol 221: 480–487, 1971.[Free Full Text]
  28. Strang G. Introduction to Linear Algebra. Cambridge, MA: Wellesley-Cambridge, 1998.



This article has been cited by other articles:


Home page
 Phil Trans R Soc AHome page
A. Porta, F. Aletti, F. Vallais, and G. Baselli
Multimodal signal processing for the analysis of cardiovascular variability
Phil Trans R Soc A, January 28, 2009; 367(1887): 391 - 409.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow All Versions of this Article:
294/1/H293    most recent
00852.2007v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Chen, X.
Right arrow Articles by Mukkamala, R.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Chen, X.
Right arrow Articles by Mukkamala, R.

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Visit Other APS Journals Online
Copyright © 2008 by the American Physiological Society.