Assessment of Baroreflex Control of Heart Rate During General Anesthesia Using a Point Process Method

doi:10.1109/ICASSP.2009.4959588

Proc IEEE Int Conf Acoust Speech Signal Process. Author manuscript; available in PMC 2010 May 26.

Published in final edited form as:

Proc IEEE Int Conf Acoust Speech Signal Process. 2009 May 26; 2009: 333–336.

doi: 10.1109/ICASSP.2009.4959588

PMCID: PMC2867254

NIHMSID: NIHMS136233

Assessment of Baroreflex Control of Heart Rate During General Anesthesia Using a Point Process Method

Z Chen, PL. Purdon, ET. Pierce, G Harrell, EN. Brown, and R Barbieri

Neuroscience Statistics Research Lab, Massachusetts General Hospital, Harvard Medical School / Harvard-MIT Division of Health Science and Technology, Boston, MA, USA

Address for correspondence Dr. Zhe Chen, 43 Vassar Street, Rm 46-6057 Massachusetts Institute of Technology Cambridge, MA 02139, USA ; Email: zhechen/at/mit.edu

See other articles in PMC that cite the published article.

Abstract

Evaluation of baroreflex control of heart rate (HR) has important implications in clinical practice of anesthesia and postoperative care. In this paper, we present a point process method to assess the dynamic baroreflex gain using a closed-loop model of the cardiovascular system. Specifically, the inverse Gaussian probability distribution is used to model the heartbeat interval, whereas the instantaneous mean is identified by a linear or bilinear bivariate regression on the previous R-R intervals and blood pressure (BP) measures. The instantaneous baroreflex gain is estimated in the feedback loop with a point process filter, while the RR→BP feedforward frequency response is estimated by a Kalman filter. In addition, the instantaneous cross-spectrum and cross-bispectrum (as well as their ratio) can also be estimated. All statistical indices provide a valuable quantitative assessment of the interaction between heartbeat dynamics and hemodynamics during general anesthesia.

Keywords: point processes, adaptive filters, Volterra series, bilinear systems, baroreflex control

I. Introduction

A change of heart rate (HR), in response to an inverse change in arterial blood pressure (BP), is a measure of the gain or sensitivity of the baroreceptor-cardiac reflex (baroreflex), which is essential in characterizing cardiovascular control and explaining both heartbeat dynamics and hemodynamics [9]. Evaluation of baroreflex control of heart rate during or after anesthesia has important implications for clinical practice (e.g. [5], [15]) and has attracted research attention to determine how anesthetic drugs alter cardiovascular control, and to develop quantitative measures that can be used to monitor patients under general anesthesia [12].

In previous work [3,4,6,7], we have successfully applied probabilistic point process models for estimating instantaneous indices of HR, HR variability (HRV), as well as respiratory sinus arrhythmia. The point process framework enabled us to estimate these physiological indices in a dynamic fashion with milliseconds timescale. Since the cardiovascular system has a closed-loop interaction between many variables including R-R interval and BP, research efforts have been devoted to estimating the baroreflex gain with a closed-loop system identification approach [1], [2]. In this paper, we extend the point process method to model the heartbeat interval, allowing for a dynamical assessment of the baroreflex gain in the feedback BP→RR loop. At the same time, a Kalman filter is used to track the parameters for estimating the RR→BP frequency response in the feedforward cardiovascular loop. To characterize a potential nonlinear interaction, we also model the heartbeat interval mean with a bilinear system, which allows us to estimate the dynamic cross-bispectrum between the R-R and BP measurements. We apply our point process model to real physiological recordings of two subjects during general anesthesia and conduct statistical assessment of the results.

II. Heartbeat Interval Point Process Model

Given a set of R-wave events equation M1

detected from the electrocardiogram (ECG), let RR_j=u_j−u_j−1>0 denote the j^th R-R interval. By treating the R-waves as discrete events, we may develop a probabilistic point process model in the continuous-time domain. Assuming history dependence, the waiting time t−u_t (as a continuous random variable) until the next R-wave event can be modeled by an inverse Gaussian model [3], [4], [6]:

where u_j denotes the previous R-wave event occurred before time t, θ>0 denotes the shape parameter, and μ_t [equivalent]

μ_RR(t) denotes the instantaneous R-R mean. Note that when the mean μ_t is much greater than the variance, the inverse Gaussian can be well approximated by a Gaussian model with a variance equal to equation M3

. In point process theory, the inter-event probability p (t) is related to the conditional intensity function (CIF) λ (t) by a one-to-one transformation: equation M4

. The estimated CIF can be used to evaluate the goodness-of-fit of the probabilistic heartbeat model.

A. Instantaneous Indices of HR and HRV

Heart rate is defined as the reciprocal of the R-R intervals. For time t measured in seconds, the new variable r=c(t−u_t)⁻¹ (where c = 60 s/min) can be defined in beats per minute (bpm). By the change-of-variables formula, the HR probability p(r)=p(c(t−u_t)⁻¹) is given by equation M5

, and the mean and the standard deviation of HR r can be derived [3]:

(1)

where

and

. Essentially, the instantaneous indices of HR and HRV are characterized by the mean μ_HR and standard deviation σ_HR, respectively.

B. Modeling of Instantaneous Heartbeat Interval’s Mean

In general, let us consider a causal, continuous-time nonlinear mapping F between an output variable y(t) and two input variables x(t) and u(t). Expanding the Wiener-Volterra series of function F (up to the second order) with respect to inputs x(t) and u(t) yields [13]:

where

, and a(·),b(·),h₁(·,·) and h₃(·,·) are Volterra kernels with appropriate orders. In our case, y(t) will be replaced by μ_RR(t), x(t) will be replaced by previous R-R intervals, u(t) will be replaced by BP, and the continuous-time integral will be approximated by a finite and discrete approximation.

Case 1

Dropping off all of second-order terms in the Volterra series expansion, we obtain a bivariate discrete-time linear system:

(2)

where the first two terms represent a linear autoregressive (AR) model of the past R-R intervals, a0 compensates the nonzero mean effect of the R-R measurements, and BP_t−j denotes the previous j^th BP value prior to time t. The BP in (2) can be either the systolic or the diastolic value.

Case 2

Dropping off the last two quadratic terms in the Volterra series expansion, we obtain

(3)

which yields a bivariate bilinear system (here the BP measurements are assumed to have a zero mean value) [16].

C. Adaptive Point Process Filtering

Let

denote the vector that contains all unknown parameters in the probabilistic model, we can recursively estimate them via adaptive point process filtering [4]:

where P and W denote the parameter and noise covariance matrices, respectively; Δ=5ms denotes the time bin size. Diagonal noise covariance matrix W that determines the level of parameter fluctuation at the timescale of Δ can be initialized either empirically from a random-walk theory or from a maximum likelihood estimate. Symbols equation M15

and

denotes the first- and second-order partial derivatives of the CIF w.r.t. ξ at time t = kΔ, respectively. The indicator variable n_k=1 if a heart beat occurs in time ((k−1)Δ, kΔ] and 0 otherwise.

D. Closed-loop Cardiovascular Control

a) Modeling Baroreflex Gain (BP→RR Loop)

First, we aim to assess the BP→RR feedback loop, which is directly related to the HR baroreflex. Unlike traditional methods (e.g. [8], in which the baroreflex gain was estimated based on the coherence between HR and BP), the baroreflex is estimated through a closed-loop bivariate AR model [2]. Moreover, our point process model is adaptive and sidesteps the local stationarity assumption, therefore it is capable of capturing the non-stationary nature of the physiological signals due to the drastic cardiovascular control compensatory changes. Specifically, in light of (2) we can compute the frequency response for the baroreflex (BP→RR loop)

(4)

where f₁and f₂ denote the rate for the R-R and BP-BP intervals, respectively; here we assume f₁≈ f₂ [equivalent]

f. With the estimated time-varying AR coefficients {a_i(k)}and{b_j(k)} at time t = kΔ, we may evaluate the dynamic frequency response of (4) at different ranges (LF, 0.04-0.15 Hz; HF, 0.15-0.5 Hz). The baroreflex gain, characterized by|H₁₂(f)|, represents the effect of BP on heartbeat, mediated by the neural autonomic reflex. Given the baroreflex gain, we can estimate the cross-spectrum (between BP and RR) as equation M18

. When the coefficients {a_i(t)}and{b_j(t)}are iteratively updated, the point process filter produces an assessment of instantaneous (parametric) baroreflex gain as well as cross-spectrum at a very fine temporal resolution without using the window technique.

b) Modeling RR→BP Feedforward Loop

Simultaneous to baroreflex assessment, we aim to model the RR→BP feedforward loop, which enables us to evaluate the impact of heartbeat on the hemodynamics. Similar to (2), BP is also modeled by a bivariate linear AR model:

(5)

where μ_RR(k − i) represents the estimated instantaneous R-R mean value at the time when BP-events occur. The coefficients equation M20

and

will be dynamically tracked by a Kalman filter or recursive least-squares (RLS) filter. Unlike the point process filter, the update occurs only at the time of BP-events. Similarly, we can also estimate the frequency response of the RR→BP cardiovascular loop:

(6)

where f denotes the sampling rate (beat/sample) for BP-BP intervals. Likewise, we can estimate the dynamic gain and phase of H₂₁(f) at each single BP-event. Similarly, we can estimate the cross-spectrum: equation M23

, where the time-varying autospectrum of R-R series is given by a parametric form:

(7)

Furthermore, the instantaneous normalized cross-spectrum (i.e., coherence) can be further computed as equation M25

c) Estimating Dynamic Cross-Bispectrum

If μ_t is identified by a bilinear system as in (3), we can further estimate the higher-order statistics. For simplicity of derivation, we further assume that the R-R and BP measurements are all Gaussian (such that their own third-order cumulant statistics are zeros), and the cross-bispectrum between BP (input u) and R-R (output y) is given by:

(8)

where

denotes the Fourier transform of the 2nd-order kernel coefficients{h_kl}, and equation M28

denote the power spectra of the BP and R-R series, respectively.

Proof

We only sketch the basic steps here, the details of derivation are omitted due to space limit. For clarity of proof, we assume that two inputs u(t) and x(t) have zero means. Similar to [11], we first decompose the output y(t) into three (two linear and one bilinear) terms and derive that equation M29

Second, we compute the cross third-order cumulant between u(t) and y(t) (viz. cross bicovariance): equation M30

. Third, we compute the two-dimensional Fourier transform of C_uuy(τ₁,τ₂), which finally yields (8).

Let h(t) denote a vector that contains all of 2nd-order coefficients{h_kl(t)}; in light of (8), we may compute an instantaneous index that quantifies the fractional contribution between the cross-spectrum and the cross-bispectrum:

where |·|denotes either the norm of a vector or the modulus of a complex variable. The “ ≈ ” is due to the approximation of a Gaussian assumption used in (8). A small value of ρ implies a presence of significant (nonzero) values in {h_kl} (i.e. nonlinearity), whereas a perfect linear Gaussian model would imply ρ =1.

III. Experimental Data and Results

Two healthy volunteer subjects, ages 20 and 32, gave written consent to participate in this study approved by the Massachusetts General Hospital (MGH). Any subject whose medical evaluation did not allow him or her to be classified as American Society of Anesthesiologists Physical Status I was excluded from the study. Intravenous and arterial lines were placed in each subject. Propofol was infused intravenously using a previously validated computer-controlled delivery system running STANPUMP [14] connected to a Harvard 22 syringe pump (Harvard Apparatus, Holliston, MA). Five effect-site target concentrations (0-4 mcg/ml) were each maintained for 15 minutes respectively. In subject 2, an additional effect-site target concentration of 5 mcg/ml was administered. Capnography, pulse oximetry, ECG, and arterial BP (P1) were recorded and (sampling rate 1 kHz) and monitored continuously by an anesthesiologist throughout the study. Bag-mask ventilation with 30% oxygen was administered as needed in the event of propofol-induced apnea. Because propofol is a potent peripheral vasodilator, phenylephrine was administered intravenously to maintain mean arterial BP within 20% of the baseline value [12].

In the present experiment, systolic BP (SBP) value was used for baroreflex evaluation. For the linear model (2), the bivariate orders p and q were fitted from 2 to 8 and the optimal order was chosen according to the Akaike information criterion (AIC). For the bilinear model (3), the order r=2 was chosen empirically to avoid demanding computation burden, and the initial h_ij was estimated by fitting the residual error via least-squares. Upon estimating the CIF, the goodness-of-fit of the probabilistic heartbeat model is evaluated with the Kolmogorov-Smirnov (KS) test [3]. For all of data fitted here, our model achieves fairly satisfactory goodness-of-fit: among a total of 11 epochs (5 in subject 1 plus 6 in subject 2), the linear model is able to reach 95% confidence bounds in KS test for 8 epochs.

For each subject, we first examine the mean statistics of baroreflex gain (LF and HF) during each epoch (15 min each, upon reaching the steady state). Specifically, during the steady state we observed a clear reduction of baroreflex gain at HF from control baseline to the start of anesthesia (level 1), and it continued to decrease as the level of the drug concentration increased. This observation is also consistent with some published results in the literature [5], [15]. For HR and HRV, upon the induction of general anesthesia, we observed an opposite trend between two subjects (for subject 1, HR increased and HRV slightly decreased). Results from two subjects across all 11 epochs are summarized in Fig. 1.

Figure1

Summarized estimated mean statistics for two subjects

Furthermore, to evaluate the tracking performance of point process filter, we also examine its performance during transient dynamics. We use subject 1 to illustrate this point in Fig. 2. As seen, the baroreflex responses were triggered by injections of drug around 1890 s, which is accompanied by a drop in the baroreflex gain (about 55%), and the point process filter quickly captures the change. Meanwhile, the instantaneous cross-spectrum (or coherence) between SBP and RR shows that these two series are strongly correlated at the HF range, first staying around 0.3 Hz at the conscious baseline and then shifting around 0.25 Hz at level-1 of drug concentration (now shown here).

Figure 2

A snapshot of estimated dynamic statistical indices ranging from the baseline to level-1 (subject 1). The dashed line marks the start of anesthesia at level-1, and the blank region represents the transient period after the baseline (where SBP measures (more ...)

Next, we also investigate the role of nonlinearity that is played by the bilinear model before and during general anesthesia. Specifically, we compare the mean ratio statistic between the conscious baseline and the level-1 drug concentration, and the result is listed in Table I. For both subjects, the ρ value (LF and HF) is significantly greater in the conscious baseline condition (P<0.01, Mann-Whitney test), which suggests that the bilinear interaction between BP and RR became more active during general anesthesia, where the parasympathetic activity is suppressed or attenuated [10] (this phenomenon is also consistent with our observations in another experimental protocol [7]). Meanwhile, the reduction of mean coherence (see Table I) during anesthesia also suggests that the BP-RR relation might have either nonlinear components, or two signals are (relatively) less linearly correlated.

TABLE I

Comparison of mean HR, SBP, coherence, and ρ-value

IV. Conclusion

We have developed a point process method for assessing the baroreflex control of heart rate during general anesthesia using clinical recordings. The proposed point process method enables us to estimate instantaneous HR, HRV, cross-spectrum, and cross-bispectrum, all of which may serve as useful indicators in ambulatory monitoring for clinical practice. The empirical results have demonstrated that the baroreflex responses were reset during general anesthesia to allow faster HR at lower BP than during consciousness and that the quantitative baroreflex gain (esp. at HF) decreased dramatically after administration of anesthesia. The change in HR/BP set point can be attributed to propofol’s systemic vasodilatory effect, whereas the reduction in baroreflex gain is most likely the result of disruption of areas within the central nervous system responsible for cardiac control. These preliminary results encourage future effort to collect more data and evaluate the proposed method in a group study for testing the significance of the result.

Acknowledgement

This work was supported by NIH Grants R01-HL084502, K25-NS05758, DP1-OD003646,.and R01-DA015644.

The authors are with the Neuroscience Statistics Research Laboratory, Massachusetts General Hospital, Harvard Medical School, Boston, MA 02114, USA. E. N. Brown is also with the Harvard-MIT Division of Health Science and Technology and the Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA.

The authors thank L. Citi, K. Habeeb, R. Merhar, A. Salazar and C. Tavares for assistance in collecting and preprocessing the data used in our experiments.

References

[1] Barbieri R, Waldmann RA, Di Virgilio V, et al. Continuous quantification of baroreflex and respiratory control of heart rate by use of bivarate autoregressive techniques. Annals of Noninvasive Electrocardiology. 1996;vol. 3:264–277.

[2] Barbieri R, Parati G, Saul JP. Closed- versus open-loop assessment of heart rate baroreflex. IEEE Eng. Med. Biol. 2001;vol. 20(no. 2):33–42. [PubMed]

[3] Barbieri R, Matten EC, Alabi AA, Brown EN. A pointprocess model of human heartbeat intervals: new definitions of heart rate and heart rate variability. Am J. Physiol. Heart Cicr. Physiol. 2005;vol. 288:424–435. [PubMed]

[4] Barbieri R, Brown EN. Analysis of heart beat dynamics by point process adaptive filtering. IEEE Trans. Biomed. Eng. 2006;vol. 53(no. 1):4–12. [PubMed]

[5] Carter JA, Clarke TNS, Prys-Roberts C, Spelina KR. Restoration of baroreflex control of heart rate during recovery from anaesthesia. British J. Anaesth. 1986;vol. 58:415–421. [PubMed]

[6] Chen Z, Brown EN, Barbieri R. Proc. ICASSP’08. Las Vegas, USA: A study of probabilistic models for characterizing human heart beat dynamics in autonomic blockade control; pp. 481–484.

[7] Chen Z, Brown EN, Barbieri R. Proc. IEEE EMBC’08. Vancouver, Canada: Characterizing nonlinear heartbeat dynamics within a point process framework; pp. 2781–2784.

[8] Clayton RH, Bowman AJ, Murray A. Measurement of baroreflex gain from heart rate and blood pressure spectra. Physiol. Meas. 1995;vol. 16:131–139. [PubMed]

[9] Eckberg DL. Arterial baroreflexes and cardiovascular modeling. Cardiovasc. Eng. 2008;vol. 8:5–13. [PubMed]

[10] Feld J, Hoffman W, Paisansathan C, et al. Autonomic activity during dexmedetomidine or fentanyl infusion with desflurane anesthesia. J. Clinical Anesthesia. 2003;vol. 19(no. 1):30–36. [PubMed]

[11] Nikias C, Petropulu AP. Higher Order Spectra Analysis: A Non-Linear Signal Processing Framework. Prentice Hall; 1993.

[12] Purdon PL, Pierce ET, Bonmassar G, Walsh J, Harrell G, et al. Simultaneous electroencephalography and functional magnetic resonance imaging of general anesthesia. Annals New York Acad. Sci. in press. [PMC free article] [PubMed]

[13] Schetzen M. The Volterra and Wiener Theories of Nonlinear Systems. Wiley; 1980.

[14] Shafer A, Doze VA, Shafer SL, et al. Pharmacokinetics and pharmacodynamics of propofol infusions during general anesthesia. Anesthesiology. 1988;vol. 69:348–356. [PubMed]

[15] Tanaka M, Nagaski G, Nishikawa T. Moderate hypothermia depresses arterial baroreflex control of heart rate during, and delays its recovery after general anesthesia in humans. Anesthesiology. 2001;vol.95:51–55. [PubMed]

[16] Tsoulkas V, Koukoulas P, Kalouptsidis N. Identification of input output bilinear systems using cumulants. IEEE Trans. Signal Processing. 2001;vol. 49(no. 11):2753–2761.