Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2867254

Formats

Article sections

- Abstract
- I. Introduction
- II. Heartbeat Interval Point Process Model
- III. Experimental Data and Results
- IV. Conclusion
- References

Authors

Related links

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.4959588PMCID: PMC2867254

NIHMSID: NIHMS136233

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: ude.tim@nehcehz

See other articles in PMC that cite the published article.

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.

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.

Given a set of R-wave events ${\left\{{u}_{j}\right\}}_{j=1}^{J}$ detected from the electrocardiogram (ECG), let *RR _{j}*=

$$p\left(t\right)={\left(\frac{\theta}{2\pi {t}^{3}}\right)}^{{\scriptstyle \frac{1}{2}}}\mathrm{exp}\left(-\frac{\theta {(t-{u}_{t}-{\mu}_{t})}^{2}}{2{\mu}_{t}^{2}(t-{u}_{t})}\right)(t>{u}_{t}),$$

where *u _{j}* denotes the previous R-wave event occurred before time

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}*)

$${\mu}_{\mathit{HR}}={\stackrel{~}{\mu}}^{-1}+{\stackrel{~}{\theta}}^{-1},{\sigma}_{\mathit{HR}}=\sqrt{\left(2\stackrel{~}{\mu}+\stackrel{~}{\theta}\right)/\stackrel{~}{\mu}{\stackrel{~}{\theta}}^{2}}$$

(1)

where $\stackrel{~}{\mu}={c}^{-1}{\mu}_{\mathit{RR}}$ and $\stackrel{~}{\theta}={c}^{-1}\theta $. Essentially, the instantaneous indices of HR and HRV are characterized by the mean *μ _{HR}* and standard deviation

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]:

$$\begin{array}{c}\hfill y\left(t\right)=F\left(x\right(t),u(t\left)\right)={\int}_{0}^{t}a\left(\tau \right)x(t-\tau )\mathrm{d}\tau +{\int}_{0}^{t}b\left(\tau \right)u(t-\tau )\mathrm{d}\tau +{\int}_{0}^{t}{\int}_{0}^{t}{h}_{1}({\tau}_{1},{\tau}_{2})x(t-{\tau}_{1})u(t-{\tau}_{2})\mathrm{d}{\tau}_{1}\mathrm{d}{\tau}_{2}\hfill \\ \hfill +{\int}_{0}^{t}{\int}_{0}^{t}{h}_{2}({\tau}_{1},{\tau}_{2})x(t-{\tau}_{1})x(t-{\tau}_{2})\mathrm{d}{\tau}_{1}\mathrm{d}{\tau}_{2}+{\int}_{0}^{t}{\int}_{0}^{t}{h}_{3}({\tau}_{1},{\tau}_{2})u(t-{\tau}_{1})u(t-{\tau}_{2})\mathrm{d}{\tau}_{1}\mathrm{d}{\tau}_{2}\hfill \end{array}$$

where $F():{\mathbb{R}}^{2}\mathbb{R}$, and *a*(·),*b*(·),*h*_{1}(·,·) and *h*_{3}(·,·) 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.

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

$${\mu}_{t}={a}_{0}\left(t\right)+\sum _{i=1}^{p}{a}_{i}\left(t\right){\mathit{RR}}_{t-i}+\sum _{j=1}^{p}{b}_{j}\left(t\right){\mathit{BP}}_{t-j}$$

(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

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

$${\mu}_{t}={a}_{0}\left(t\right)+\sum _{i=1}^{p}{a}_{i}\left(t\right){\mathit{RR}}_{t-i}+\sum _{j=1}^{p}{b}_{j}\left(t\right){\mathit{BP}}_{t-j}+\sum _{i=1}^{r}\sum _{j=1}^{r}{h}_{\mathit{ij}}\left(t\right)({\mathit{RR}}_{t-i}-<\mathit{RR}>){\mathit{BP}}_{t-j}$$

(3)

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

Let $\xi ={\left[{\left\{{a}_{i}\right\}}_{i=0}^{p},{\left\{{b}_{j}\right\}}_{j=1}^{p},\left\{{h}_{ij}\right\},\theta \right]}^{T}$ denote the vector that contains all unknown parameters in the probabilistic model, we can recursively estimate them via adaptive point process filtering [4]:

$$\begin{array}{cc}\hfill & {\xi}_{kk-1}\hfill & {P}_{kk-1}\hfill & {\xi}_{kk}\hfill \hfill \hfill \end{array}$$

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 $\Delta {\lambda}_{k}={\scriptstyle \frac{{\lambda}_{k}{\xi}_{k}}{}}$ and ${\Delta}^{2}{\lambda}_{k}={\scriptstyle \frac{{2}^{{\lambda}_{k}}{\xi}_{k}{\xi}_{k}^{T}}{}}$ 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 ((

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)

$${H}_{12}\left(f\right)=\frac{{\phantom{\sum _{j=1}^{q}{b}_{j}\left(k\right){z}^{-j}z={e}^{j2\pi {f}_{2}}}}_{{\phantom{1-\sum _{i=1}^{p}{a}_{i}\left(k\right){z}^{-i}z={e}^{j2\pi {f}_{1}}}}_{},}}{}$$

(4)

where *f*_{1}and *f*_{2} denote the rate for the R-R and BP-BP intervals, respectively; here we assume *f*_{1}≈ *f*_{2} *f*. With the estimated time-varying AR coefficients {*a _{i}*(

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:

$${\mathit{BP}}_{k}={c}_{0}+\sum _{i=1}^{p}{c}_{i}{\mathit{BP}}_{k-i}+\sum _{i=1}^{p}{d}_{i}{\mu}_{\mathrm{RR}}(k-i),$$

(5)

where *μ*_{RR}(*k − i*) represents the estimated instantaneous R-R mean value at the time when BP-events occur. The coefficients ${\left\{{c}_{i}\right\}}_{i=0}^{p}$ and ${\left\{{d}_{i}\right\}}_{i=1}^{p}$ 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:

$${H}_{21}\left(f\right)=\frac{{\phantom{\sum _{i=1}^{p}{d}_{i}\left(k\right){z}^{-i}z={e}^{j2\pi f}}}_{{\phantom{1-\sum _{i=1}^{p}{c}_{i}\left(k\right){z}^{-i}z={e}^{j2\pi f}}}_{},}}{}$$

(6)

where *f* denotes the sampling rate (beat/sample) for BP-BP intervals. Likewise, we can estimate the dynamic gain and phase of *H*_{21}(*f*) at each single BP-event. Similarly, we can estimate the cross-spectrum: ${\mathcal{C}}_{uy}\left(f\right)={H}_{21}\left(f\right){\mathcal{Q}}_{\mathrm{RR}}\left(f\right)$, where the time-varying autospectrum of R-R series is given by a parametric form:

$${\mathcal{Q}}_{\mathrm{RR}}(f,t)=\frac{{\sigma}_{\mathrm{RR}}^{2}\left(t\right)}{1-\sum _{i=1}^{p}{a}_{i}\left(t\right){z}^{-i}{z={e}^{j2\pi f}}_{}.}$$

(7)

Furthermore, the *instantaneous* normalized cross-spectrum (i.e., coherence) can be further computed as $\mathsf{Coh}(f,t)=\frac{{\mathcal{C}}_{uy}\left(f\right)}{\sqrt{{\mathcal{Q}}_{\mathrm{BP}}(f,t){\mathcal{Q}}_{\mathrm{RR}}(f,t)}}$.

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

$${\mathcal{C}}_{uuy}({f}_{1},{f}_{2})=2\mathcal{H}(-{f}_{1},-{f}_{2}){H}_{21}\left({f}_{1}\right){\mathcal{Q}}_{\mathrm{RR}}\left({f}_{1}\right){\mathcal{Q}}_{\mathrm{RR}}\left({f}_{2}\right)=2\mathcal{H}(-{f}_{1},-{f}_{2}){H}_{12}\left({f}_{1}\right){\mathcal{Q}}_{\mathrm{BP}}\left({f}_{1}\right){\mathcal{Q}}_{\mathrm{RR}}\left({f}_{2}\right)$$

(8)

where $\mathcal{H}({f}_{1},{f}_{2})={\sum}_{k=1}^{q}{\sum}_{l=1}^{q}{h}_{kl}{e}^{-j2k\pi {f}_{1}}{e}^{-j2l\pi {f}_{2}}$ denotes the Fourier transform of the 2nd-order kernel coefficients{*h _{kl}*}, and ${\mathcal{Q}}_{\mathrm{BP}}\left(f\right)\text{and}{\mathcal{Q}}_{\mathrm{RR}}\left(f\right)$ denote the power spectra of the BP and R-R series, respectively.

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 $\mathbb{E}\left[y\left(t\right)\right]={\sum}_{i}{\sum}_{j}{h}_{ij}\mathbb{E}\left[x(t-i)u(t-j)\right]={\sum}_{i}{\sum}_{j}{h}_{ij}{C}_{ux}(i-j)={\scriptstyle \frac{1}{2\pi}}\int \mathcal{H}(f,-f){H}_{12}\left(f\right){\mathcal{Q}}_{\mathrm{BP}}\left(f\right)\mathrm{d}f$

Second, we compute the cross third-order cumulant between *u(t)* and *y(t)* (viz. cross bicovariance): ${C}_{uuy}({\tau}_{1},{\tau}_{2})=\mathbb{E}\left\{u(t+{\tau}_{1})u(t+{\tau}_{2})[y\left(t\right)-\mathbb{E}\left[y\left(t\right)\right]]\right\}$. Third, we compute the two-dimensional Fourier transform of *C _{uuy}*(

Let **h**(*t*) denote a vector that contains all of 2nd-order coefficients{*h _{kl}*(

$${\rho}_{t}=\frac{{\mathcal{C}}_{uy}(f,t){\mathcal{C}}_{uy}(f,t)+{\mathcal{C}}_{uuy}({f}_{1},{f}_{2},t)\approx \frac{1}{1+2\mathbf{h}\left(t\right){\mathcal{Q}}_{\mathrm{RR}}(f,t),}}{}$$

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

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

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.

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).

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 **...**

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.

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.

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.

[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. [PMC free article] [PubMed]

[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. [PMC free article] [PubMed]

[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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |