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

**|**HHS Author Manuscripts**|**PMC3059729

Formats

Article sections

- Abstract
- I. Introduction
- II. A Probability Model for Heartbeat Intervals
- III. Frequency Analysis for Differential BRS
- IV. Data and Results
- V. Summary and Discussion
- References

Authors

Related links

Conf Proc IEEE Eng Med Biol Soc. Author manuscript; available in PMC 2011 March 17.

Published in final edited form as:

PMCID: PMC3059729

NIHMSID: NIHMS247729

Z Chen, Senior Member, IEEE, PL Purdon, EN Brown, Fellow, IEEE, and R Barbieri, Senior Member, IEEE^{}

The publisher's final edited version of this article is available at Conf Proc IEEE Eng Med Biol Soc

See other articles in PMC that cite the published article.

Modeling heartbeat variability remains a challenging signal-processing goal in the presence of highly non-stationary cardiovascular control dynamics. We propose a novel differential autoregressive modeling approach within a point process probability framework for analyzing R-R interval and blood pressure variations. We apply the proposed model to both synthetic and experimental heartbeat intervals observed in time-varying conditions. The model is found to be extremely effective in tracking non-stationary heartbeat dynamics, as evidenced by the excellent goodness-of-fit performance. Results further demonstrate the ability of the method to appropriately quantify the non-stationary evolution of baroreflex sensitivity in changing physiological and pharmacological conditions.

Recently, point process probability models have been advocated for characterizing human heartbeat intervals [1, 2, 3, 4, 5]. Unlike previous methods, the point process paradigm allows to estimate *instantaneous* heart rate (HR) and HR variability (HRV), as well as specific cardiovascular/cardiorespiratory functions such as respiratory sinus arrhythymia (RSA) or baroreflex sensitivity (BRS). In order to track the non-stationary heartbeat dynamics, an adaptive point process filtering approach has been proposed to trail the instantaneous model parameters [2, 4]. In this work, we propose a distinctive modeling perspective to tackle the non-stationary nature of the heartbeat intervals as well as potential other physiological covariates (such as blood pressure or respiration measures). The main feature of the proposed algorithm is the inclusion of a linear regression on the mean of the point process probability density by use of an autoregressive integrated moving average (ARIMA) model. Such framework also allows for a newly defined *differential* BRS index. Notably, the new index is conceptually similar to the time-domain “sequence method” [6, 9, 10], which has been often used to measure the sensitivity of the change of blood pressure (BP) relative to the change of R-R interval. However, more suitably in the presence of dynamic changes in a non-stationary environment, our index provides an instantaneous BRS characterization at arbitrarily small time resolutions.

The paper is organized as follows. We first present an overview of the point process modeling paradigm, the new frequency analysis framework associated with the ARIMA model, and the definition of the BRS index. We then apply the proposed model to both synthetic heartbeat series and experimental data recordings consisting of non-stationary heartbeat intervals and systolic BP during a propofol induction of anesthesia protocol [11] and a tilt-table protocol [1]. We demonstrate the noticeable improvement in model goodness-of-fit while fitting highly non-stationary recordings. Additionally, we use these examples to illustrate the ability of our new model to estimate the instantaneous *differential* BRS index. Finally, we conclude the paper with some discussions.

Given a set of R-wave events * {u_{j}}_{j=1}^{J}* detected from the electrocardiogram (ECG), let

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

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

$$p(t)={\left(\frac{1}{2\pi {\sigma}_{\text{RR}}^{2}(t)}\right)}^{\frac{1}{2}}\text{exp}\left(-\frac{{(t-{u}_{j}-{\mu}_{\text{RR}}(t))}^{2}}{2{\sigma}_{\text{RR}}^{2}(t)}\right),$$

In point process theory, the inter-event probability ** p(t)** is related to the conditional intensity function (CIF)

Heart rate (HR) is defined as the reciprocal of the R-R intervals. For time ** t** measured in seconds, the new variable

$${\mu}_{\text{HR}}={\tilde{\mu}}^{-1}+{\tilde{\theta}}^{-1},\text{}{\sigma}_{\text{HR}}=\sqrt{(2\tilde{\mu}+\tilde{\theta})/\tilde{\mu}{\tilde{\theta}}^{2}},$$

(1)

where = *c*^{−1} μ_{RR} and = *c*^{−1} θ. Essentially, the instantaneous indices of HR and HRV are characterized by the mean **μ**_{HR} and standard deviation **σ**_{HR}, respectively.

In general, the heartbeat interval can be modeled by Wiener-Volterra series expansion, using previous R-R intervals and other physiological covariates (e.g., blood pressure or respiration) as the input variables. Although inclusion of nonlinear or bilinear terms in modeling is possible [3, 5], in this paper we focus on the investigation on linear modeling approaches.

We can use the bivariate AR model to model the instantaneous mean **μ**_{RR}** (t)** [3]:

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

(2)

where ** a_{0}** compensates the nonzero mean effect of the R-R measurements,

In time-series (linear) modeling, it is a common practice to “detrend” a time series by taking differences if the series exhibits undesired non-stationary features. The autoregressive integrated moving average (ARIMA) process may provide a suitable framework to achieve such goal from a modeling point of view [12]. Simply, the original time series is applied by a difference operator (one or more times) until the non-stationary trends are not observed in the final series of interest. Motivated by this idea, we here define the “increment of R-R series” *{δRR _{t−i}}*

$${\mu}_{\text{RR}}(t)={\mathit{\text{RR}}}_{t-1}+{\displaystyle \sum _{i=1}^{p}{a}_{i}(t)\delta {\mathit{\text{RR}}}_{t-i}+}{\displaystyle \sum _{j=1}^{p}{b}_{j}(t)\delta {\mathit{\text{BP}}}_{t-j}}.$$

(3)

Note that the ** a_{0}(t)** term in (2) has been replaced by

Let **ξ** = ** [{a_{i}}_{i=1}^{p},{b_{j}}_{j=1}^{p},θ]^{T}** denote the vector that contains all unknown parameters from the new model (3), we can recursively estimate them via adaptive point process filtering [2]:

$$\begin{array}{cc}\hfill & {\xi}_{k|k-1}={\xi}_{k-1|k-1}\hfill \\ \hfill & {P}_{k|k-1}={P}_{k-1|k-1}+W\hfill \\ \hfill & {\xi}_{k|k}={\xi}_{k|k-1}+{P}_{k|k-1}(\nabla \text{log}{\lambda}_{k})[{n}_{k}-{\lambda}_{k}\mathrm{\Delta}]\hfill \\ {P}_{k|k}=\hfill & {\left[{P}_{k|k-1}^{-1}+\nabla {\lambda}_{k}\nabla {\lambda}_{k}^{T}\frac{\mathrm{\Delta}}{{\lambda}_{k}}-{\nabla}^{2}\text{log}{\lambda}_{k}[{n}_{k}-{\lambda}_{k}\mathrm{\Delta}]\right]}^{-1}\hfill \end{array}$$

where ** P** and

The goodness-of-fit of the point process model is tested based on the *time-rescaling theorem* [1, 2]. Given a point process specified by ** J** discrete events:

In addition, the autocorrelation function of the ** g_{j}**s, defined by $\mathit{\text{ACF}}\phantom{\rule{thinmathspace}{0ex}}(m)=\frac{1}{J-m}{\displaystyle {\sum}_{j=1}^{J-m}{g}_{j}{g}_{j+m}}$
, is computed. If the

In light of equation (3), we can derive the instantaneous BP**→**RR transfer function in the feedback loop, which essentially relates to the *differential* BRS

$$H(f,t)=\frac{{\displaystyle {\sum}_{i=1}^{q}{b}_{i}(t){z}^{-i}{|}_{z={e}^{j2\pi {f}_{2}}}}}{1-{\displaystyle {\sum}_{i=1}^{p}{a}_{i}(t){z}^{-i}{|}_{z={e}^{j2\pi {f}_{1}}}}}.$$

(4)

Note that the new form of transfer function is reminiscent of the “sequence method” [6, 9, 10] previously used to define the BRS. However, the sequence method employed a time-domain batch method, while ours is a frequency-domain instantaneous method. We further assume that the effect of time discretization is small such that the linear dependence between two variables remains unchanged, i.e. $\frac{y(t)}{x(t)}\mathbf{\approx}\frac{y(t+\delta )-y(t)}{x(t+\delta )-x(t)}$. As expected, due to the difference operator used in the time domain, the new AR coefficients ** {a_{i},b_{i}}** are more sensitive to the “incremental changes”; consequently, in the frequency domain,

We simulated the heartbeat R-R series (unit: ms) using a stationary linear Gaussian model plus a linear increasing trend along time (for simplicity, we use only a univariate AR model in simulation):

$$({\mathit{\text{RR}}}_{t-1}-1000)=\alpha t+{\displaystyle \sum _{i=1}^{8}{a}_{i}({\mathit{\text{RR}}}_{t-i}-1000)+n(t)}$$

(5)

where **α**=0.001, and ** n(t)** is a Gaussian noise process with zero mean and variance of 100. As seen from the top panel of Fig. 1, the simulated R-R series has a linear increasing trend (from ~800 ms to ~1200 ms) during a 6-min interval. The overall R-R series has a global mean of 1000 ms and standard deviance of 110. The AR coefficients

We analyzed the experimental recordings from a general anesthesia experiment. The experimental protocol was approved by the Massachusetts General Hospital (MGH) Department of Anesthesia and Critical Clinical Practices Committee, the MGH Human Research Committee and the MGH General Clinical Research Center, and has been reported previously [3, 11]. Briefly, subjects were delivered several concentration levels of propofol anesthetic drug (0.0, 1.0, …, 5.0 **µ**g/ml) administered by anesthesiologists. Each epoch lasted about 15 minutes.

In examining the data, we have observed that upon propofol administration, the physiological measurements exhibit highly non-stationary dynamics (as opposed to the awake baseline steady state), which brings in the challenge of correctly modeling these data while simultaneously assessing cardiovascular control. Here, we consider only two experimental epochs during induction of anesthesia (level-1 propofol concentration) to illustrate our model. The experimental data recordings from two subjects are shown in Fig. 2. The model parameters in (3) are initialized by fitting a bivariate AR model with a least-squares method on the first 100 samples. The model order ** p** is chosen from 4 to 8 and the final optimum order is chosen according to the Akaike information criterion (AIC). The estimated instantaneous HR, HRV, and BRS indices from two subjects are shown in Figs. 3 and and4,4, along with the KS plot and autocorrelation plot. As seen, the KS plots are almost aligned with 45° line, indicating a nearly perfect model fit.

Two experimental R-R interval series and SBP series (left panels) from two recording epochs during level-1 propofol concentration. Their corresponding incremental series are also shown at the right panels. In subject 13, R-R series has a clear increasing **...**

The estimated instantaneous indices from the data of subject 13 and the associated KS and autocorrelation plots. The KS statistic is improved from 0.0681 (old model, red curve) to 0.0152 (new model, blue curve).

Finally, we considered R-R intervals and SBP recordings, from a protocol where each subject underwent a few cycles of “rest” and “tilt” posture conditions (details of the “*tilt-table*” protocol can be found in [1]). For demonstration purpose, we analyzed two representative subjects (#1 and #4). As an illustration, Figure 5 shows the estimated instantaneous indices from Subject #1 (only 5 cycles are shown). As seen, the point process filter tracks the instantaneous differential BRS with fluctuating trend similar to R-R intervals, i.e. decreasing during tilt and increasing back to baseline level as the subject returned to supine position. For statistical analysis, we averaged the differential BRS estimates (HF range) along the “rest” and “tilt” epochs (values are shown in Table I). Notably, comparing the tilt epochs with the rest epochs, the mean differential BRS is lower, in agreement with consolidated previous findings in the literature [7, 8]. As further comparison, we also computed the BRS using the sequence methods (one value for each epoch, as the sequence method cannot perform an instantaneous assessment). The slope of the linear interrelation between SBP and the next R-R interval is calculated when the correlation coefficient is at least 0.7 (violation of this condition will inevitably make the slope estimate inaccurate). As seen from the table, the significant decrease in differential BRS from rest to tilt are confirmed by both methods, with the point process measure providing a more accurate representation in terms of both timescale and methodology.

We propose a novel differential autoregressive modeling approach within a point process probability framework for analyzing R-R interval and blood pressure variations. In testing the proposed model as applied to synthetic and experimental data in highly non-stationary conditions, the model is found to be greatly effective in tracking non-stationary heartbeat dynamics, as evidenced by excellent goodness-of-fit performance. In addition to producing optimal KS statistics, the new model also includes the computational advantage of allowing for lower autoregressive orders, as the incremental series exhibit more stationary features. The differential BRS assessment proves as a valid instantaneous characterization at arbitrarily small time resolutions Overall, the new formulation allows for more precise dynamic measures in a highly non-stationary environment, and provides one further advancement towards potential realtime indicators for ambulatory monitoring and instantaneous assessment of autonomic control in clinical practice.

The authors thank L. Citi, K. Habeeb, G. Harrell, R. Merhar, E. T. Pierce, A. Salazar, C. Tavares and J. Walsh (Massachusetts General Hospital) for assistance with the anesthesia data used in this paper. We also thank Roger G. Mark and Thomas Heldt (MIT) for sharing the tilt-table data in this study.

^{*}This work was supported by NIH Grants R01-HL084502, K25-NS05758, DP2-OD006454, DP1-OD003646 and CRC UL1 RR025758.

1. Barbieri R, Matten EC, Alabi AA, Brown EN. A point-process 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]

2. 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]

3. Chen Z, Purdon PL, Pierce ET, Harrell G, Brown EN, Barbieri R. Assessment of baroreflex control of heart rate during general anesthesia using a point process method; Proc. ICASSP’09; Taiwan: 2009. pp. 333–336. [PMC free article] [PubMed]

4. Chen Z, Brown EN, Barbieri R. Assessment of autonomic control and respiratory sinus arrhythmia using point process models of human heart beat dynamics. IEEE Trans. Biomed. Eng. 2009;vol. 56:1791–1802. [PMC free article] [PubMed]

5. Chen Z, Brown EN, Barbieri R. Characterizing nonlinear heartbeat dynamics within a point process framework. IEEE Trans. Biomed. Eng. 2010;vol. 57(no. 6):1335–1347. [PMC free article] [PubMed]

6. Hughson RL, Quintin L, Annat G, Yamamoto Y, Gharib C. Spontaneous baroreflex by sequence and power spectral methods in humans. Clin. Physiol. 1993;vol. 13(no. 6):663–676. [PubMed]

7. Malliani A, Pagani M, Lombardi F, Cerutti S. Cardiovascular neural regulation explored in the frequency domain. Circulation. 1991;vol. 84:1482–1492. [PubMed]

8. Montano N, Ruscone T, Porta A, et al. Power spectrum analysis of heart rate variability to as- sess the changes in sympathovagal balance during graded orthostatic tilt. Circulation. 1994;vol. 90:1826–1831. [PubMed]

9. Oka H, Mochio S, Yoshioka M, Morita M, Inoue K. Evaluation of baroreflex sensitivity by the sequence method using blood pressure oscillations and R-R interval changes during deep respiration. Eur. Neurol. 2003;vol. 50:230–243. [PubMed]

10. Parati G, DiRienzo M, Mancia G. Dynamic modulation of baroreflex sensitivity in health and disease. Ann. NY Acad. Sci. 2001;vol. 940:469–487. [PubMed]

11. Purdon PL, et al. Simultaneous electroencephalography and functional magnetic resonance imaging of general anesthesia. Ann. NY Acad. Sci. 2009;vol. 1157:61–70. [PMC free article] [PubMed]

12. Vu KM. The ARIMA and VARIMA Time Series: Their Modelings, Analyses and Applications. 2007

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