Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Biosystems. Author manuscript; available in PMC 2009 July 1.
Published in final edited form as:
PMCID: PMC2561955

Application of Dynamic Point Process Models to Cardiovascular Control


The development of statistical models that accurately describe the stochastic structure of biological signals is a fast growing area in quantitative research. In developing a novel statistical paradigm based on Bayes’ theorem applied to point processes, we are focusing our recent research on characterizing the physiological mechanisms involved in cardiovascular control. Results from a Tilt Table study point at our statistical framework as a valid model for the heart beat, as generated from complex mechanisms underlying cardiovascular control. The point process analysis provides new quantitative indices that could have important implications for research studies of cardiovascular and autonomic regulation and for monitoring of heart rate and heart rate variability measures in clinical settings.

Keywords: Heart Rate Variability, Point Processes, Adaptive Filters, State Estimation

1. Introduction

Heart rate (HR) and heart rate variability (HRV) are important dynamic measures of the state of the cardiovascular system and the autonomic nervous system (Stauss, 2003; Task Force, 1996). Heart rate is traditionally estimated as the average of the reciprocal of the RR intervals within a specified time window, or as the number of R –wave events (heart beats) per unit time on the electrocardiogram (ECG). The R –wave events mark the electrical impulses from the heart’s conduction system that represent ventricular contractions. Hence, they are a sequence of discrete occurrences in continuous time, and as such, form a point process. Rather than modeling them to reflect the point process structure of the heart beats, most current methods either treat the heart beat RR interval series as continuous-valued signals, or convert them into continuous-valued, evenly spaced measurements for analysis by interpolation of either the RR intervals or their reciprocals. We have recently derived new definitions of HR and HRV based on an explicit point process Bayesian probability model for heart rate under the assumption that the stochastic properties of the R-R intervals are governed by an inverse Gaussian renewal model. We can estimate the time-varying inverse Gaussian parameters by either local maximum likelihood (Barbieri et al., 2005) or by adaptive point process estimation (Barbieri et al., 2006), and assess model goodness-of-fit by Kolmogorov-Smirnov tests based on the time-rescaling theorem. These models give a more physiologically sound representation of the stochastic structure in heart beat generation than those provided by current definitions and analysis methods. In particular, the adaptive filter algorithm can compute updates in an on-line fashion and at any desired temporal resolution, and it may be at the core of a new device to monitor heart beat dynamics in clinical setting such as the intensive care unit, the operating room and during labor and delivery (Fig. 1). We here show the application of our adaptive paradigm to data from ten healthy subjects during postural changes.

Figure 1
From ECG noninvasive recordings, RR interval peak can be detected, and the adaptive filter algorithm can compute instantaneous updates of heart rate and heart rate variability indices in an on-line fashion and at any desired temporal resolution. ...

2. Methods

In this section, we present the heart beat interval and the heart rate probability models, the heart beat interval model parameters, the point process adaptive filtering algorithm to derive instantaneous estimates of heart rate and heart rate variability, and the goodness-of-fit test to evaluate how well these estimates describe the stochastic structure of the R – wave events extracted from an ECG.

2.1. Point Process Probability Model of Heart Beat Intervals

Each R –wave event is initiated by a coordinated depolarization of the heart’s pacemaker cells that begins in the sino-atrial (SA) node and propagates throughout the cardiac muscle. Deterministic models of this integrate (rise of the transmembrane potential)-and-fire (depolarization) mechanism are used regularly to simulate heart beats or R –wave events (De Boer et al., 1985; Berger et al., 1986). An elementary, stochastic integrate-and-fire model is the Gaussian random walk model with drift, and the probability density of the first passage times for this random walk process, i.e., the times between threshold crossings (RR intervals), is well-known to be the inverse Gaussian (Tuckwell, 1988; Chhikara and Folks, 1989). Therefore, we assume that given any R –wave event uk, the waiting time until the next R –wave event, or equivalently, the length of the next RR interval, obeys the following history-dependent inverse Gaussian (HDIG) probability density:

equation M1

where 0 < u1 < u2 <,…,< uk <,…<uKT are the K successive R –wave event times from an ECG in an observation interval (0, T], t is any time satisfying t >uk >, Hk = {uk, wk, wk−1,…,wk−p+1}is the history of the RR intervals up to uk, wk = ukuk−1 = is the kth RR interval, equation M2 is the mean, θp+1 > 0 is the scale parameter, and θ =(θ0, θ1, …, θp+1) is the vector of model parameters. The autoregression on the mean allows for consideration of the effect of the recent history of the sympathetic and parasympathetic inputs to the SA node.

The mean and standard deviation of the RR probability model in (1) are, respectively,

μRR = μ(Hk, θ)
equation M4

Heart rate is often defined as the reciprocal of the RR intervals, thus we define r = c(tuk)−1 as the heart rate random variable and use the standard change-of-variables formula from elementary probability theory to derive the mean and standard deviation of the heart rate probability density defined as

equation M5
equation M6

where equation M7.

To track the non-stationary behavior in heart beat dynamics that occurs due to changes in state under both physiological and pathological conditions, we assume that the parameter θ is time-varying, and we model the time-varying behavior of θ using a state space model. To define the state model and the observation model, we choose J large, and divide (0, T] into J intervals of equal width Δ = T / J, so that there is at most one spike per interval. The adaptive parameter estimates will be updated at jΔ for j = 1,…, J.

From the heart beat probability model in (1) we define the associated conditional intensity function as

equation M8

The conditional intensity function provides a canonical characterization of a point process that gives a history-dependent generalization of the rate function of a Poisson process (Brown et al., 2003).

Once the state model and the observation process model are defined (Barbieri et al., 2006), it follows from (Barbieri et al., 2004; Eden et al., 2004) that the point process adaptive filter algorithm for this system is (One-Step Prediction)

θj|j−1 = θj−1|j−1

(One-Step Prediction Variance)


(Posterior Mode)

θj|j = θj|j−1Wj|j−1(∇logλj)[nj − λjΔ]

(Posterior Variance)

equation M12

where λj =λ (jΔ | Hj, θj|j−1) and [nabla]([nabla]2) denotes the first (second) derivative of the indicated function with respect to θ for j =1,…,J. The notation θ j|k defines the state at time, jΔ given the observations from (0, kΔ].

Given θj|j, the point process adaptive filter estimate of θ at time jΔ, it follows from (2), (3), (4) and (5), that the instantaneous estimates of mean RR, RR interval standard deviation, mean heart rate and heart rate standard deviation at time jΔ are respectively (Mean R-R Interval)

μRR(jΔ) = μ(Hj, θj|j)

(R-R Interval SD)

equation M14

(Mean Heart Rate)

equation M15

(Heart Rate SD)

equation M16

Equation (13) provides a new algorithm for computing an instantaneous estimate of heart rate whereas (12) and (14) are new indices of heart rate variability.

2.1 Model Goodness-of-Fit

To evaluate goodness-of-fit for the heart beat probability model-point process adaptive filter algorithm, i.e., determine how well this model describes the sequence of ECG R -wave events, we use the Kolmogorov-Smirnov test based on the time-rescaling theorem for point processes (Brown et al., 2002, Barbieri et al., 2005; Barbieri et al., 2006). Close agreement between the uniform transformation of the ordered observations (empirical quantiles) and the ordered observations from a uniform probability density (model quantiles) is true if and only if there is close agreement between the point process probability model and the series of R-R intervals. Graphically, the better the model describes the data, the more the KS plot lies along the 45° line and inside the 95% confidence bounds.

To assess the correlation structure and hence, possible dependence that may be present beyond first order, we also perform a normal transformation of the observations and compute the serial correlation function for 60 lags (~ 1 min). Small values of the serial correlation function at all lags would suggest that the series of the normal transformation of the observations are independent. Approximate independence suggests that the original model is highly consistent with the R-R interval series.

3. Analysis of Heart Beat Interval Dynamics

To illustrate the application of our point process adaptive filter algorithm to actual data, we have analyzed heart beat series from ten healthy subjects performing a tilt table study.

The protocol begins with subjects lying supine for approximately 10 minutes, after which, each subject undergoes 3 types of up-down tilt pairs. The tilt pairs are: rapid up (down) tilt in which the tilt table moved from horizontal (vertical) to vertical (horizontal) in less than 3 s; slow up (down) tilt in which the tilt table moved from horizontal (vertical) to vertical (horizontal) in approximately 1 min; and stand-up (supine) in which the subject stood up immediately supporting his or her own weight and then lied supine immediately from having been standing supporting his own weight (Heldt et al, 2002, 2003).

The wide range of changes in heart rate that can be induced in either short or long time periods make the tilt table protocol and excellent paradigm for testing the ability of our new algorithm to track RR interval and heart rate dynamics. Because the subject's legs are not supported, weight bearing from standing and muscle tone in the legs, which normally help return blood to the heart, are diminished. As a consequence, to maintain perfusion, there are compensatory changes in cardiovascular and autonomic activity reflected as changes in heart rate and blood pressure. Similarly, by tilting the subject back from the vertical to the horizontal, a different set of compensatory changes in cardiovascular and autonomic activity can be induced. The degree of compensatory cardiovascular and autonomic activity can be affected by changing the rates at which the postural changes are induced and/or the length of time the subject is maintained in the supine or upright posture.

Following a preliminary partial autocorrelation analysis of the RR interval series, we implemented the point process adaptive filter algorithm assuming an AR(8) HDIG model to analyze these data. The initial parameter values were computed by local maximum likelihood analysis (Barbieri et al, 2005) using the first 60 s. We chose the state covariance matrix Wε to be the diagonal matrix with values of 6×10−2 for θ0, 3×10−1 for equation M17, and 1×10−6for the autoregressive parameters by applying the same preliminary local maximum likelihood analysis.

Figure 2 shows results from one subject. In this recording session two rapid tilts were initiated at around 1000 s and 2900 s, slow tilts at 350 s and 2450 s, whereas the subjects stood up at around 1560 s and 2010 s. Our analysis provides time-varying estimates of mean RR (Fig. 2A), heart rate (Fig. 2B), RR standard deviation (Fig. 2C) and heart rate standard deviation (Fig. 2D) computed using (1114). While the mean estimates show dynamics consistent with the original beat series, and with the hemodynamic changes elicited by changes in posture, the variability instantaneous indices show very fine fluctuations not directly correlated with posture changes, particularly with significantly faster dynamics than the respective log likelihood instantaneous indices previously estimated (Barbieri et al, 2005). Most importantly, the RR and heart rate standard deviations show that time-varying estimates of mean RR and heart rate alone do not capture all the dynamic features of heart beat interval.

Figure 2
Instantaneous time varying estimates of A. the mean RR, B. heart rate, C. RR standard deviation, and D. heart rate standard deviation computed for Subject 1 as defined in Eqs. (1114). E. Kolmogorov-Smirnov (KS) plots of ...

The goodness-of-fit analysis resulted in the KS plot in Fig. 2E, and the autocorrelation plot in Fig. 2F. The KS plot for the model fit does not lie entirely within the 95% confidence bounds. Nevertheless, the curve does not significantly deviate from the diagonal, keeping all values relatively close to the confidence bounds. Looking at the autocorrelation function of the transformed times, only three lags (lag 1, 2, and 9) have autocorrelation estimates that fall significantly outside the 95% confidence bounds, suggesting that its transformed times are nearly independent. In summary, the goodness-of-fit analysis suggests that the AR(8) HDIG model with time-varying parameters describes well the stochastic structure in this subject’s heart beat series, but also invites for exploration of more refined algorithms and more complex models that could explain more exhaustively the mechanisms behind heart beat dynamics.

Figure 3 better illustrates the different pattern changes of our four measures of heart rate dynamics as elicited by each postural change. Here we show, from top to bottom respectively, the instantaneous estimates of mean RR, heart rate, RR standard deviation and heart rate standard deviation from subjects 1, 4, and 5 during approximately 200 s before and after a rapid tilt (Fig. 3A), a slow tilt (Fig. 3B), and an upright position (Fig. 3C). Of relevance in this new adaptive analysis is the increase in both RR and heart rate standard deviations immediately after the rapid tilt and after standing up. The relatively constant level of heart rate standard deviation after the slow tilt, when compared to the sharp increase immediately after rapid tilt and standing up, suggests new patterns of faster dynamics which may be useful to distinguish physiological differences in these states.

Figure 3
Heart rate and heart rate variability estimates during the three gravitational transients. A. Rapid tilt initiated at 1000 s for Subject 1, B. Slow tilt initiated at 2353 s for Subject 4, C. Stand up initiated at 2644 s for Subject 5. Transitions are ...

4. Analysis of Postural Changes

The results in Fig. 3 suggest that our heart rate variability indices help provide a more detailed description of heart beat interval dynamics for an individual in each of the three types of postural change. To illustrate how our measures can be used to help characterize heart rate variability for the population, we analyzed for all 10 subjects one of the periods in the study where the most dynamics in the heart beat intervals would be expected. That is the period prior to the postural change while the subject was in the supine position to the period immediately after the subject assumed an upright posture by a rapid tilt, slow tilt or standing up. For each type of postural change, we define this period of maximal dynamics as beginning 60 seconds before onset of the postural change and extending to 180 s after the onset (Fig. 4). For each subject, we compute each of our four indices during this dynamic period and summarize the results as the average and standard deviation. We term the time interval 60 s after the onset of the postural change as the transient period. We report the four series of RR, heart rate, RR interval standard deviation and heart rate standard deviation respectively as ΔRR, ΔHR, ΔRRV and ΔHRV because we subtract from each series its mean level during the 60s prior to the postural change.

Figure 4
Postural Changes. Summary analysis of point process adaptive heart rate (A.RR interval mean and B. heart rate mean) and heart rate variability (C.RR interval standard deviation and D. heart rate standard deviation) for the 10 tilt ...

Because we subtract the respective mean levels, all four estimates fluctuate during this period about zero (Fig. 4). At the initiation of the rapid tilt (Fig. 4A and 4B, top graph), ΔRRHR) decreases (increases) abruptly and remains at an almost stationary level for the balance of the tilt. In the vertical position there is a wider range of fluctuations, reflected in the standard deviation range (gray curves), in both ΔRR and ΔHR compared to rest. Both ΔRRV and ΔHRV (Fig. 4C and 4D, top graph) show initial abrupt increases with the onset of the rapid tilt. As the subject stabilizes in the vertical position at around 200 s, ΔRRV fluctuations also stabilize below rest levels, whereas the ΔHRV mean level goes back to rest level. Of note, both of these indices show large fluctuations about their mean levels, both in the transient period and in the period immediately following the transient period.

During the slow tilt there is a gradual decrease (increase) in ΔRRHR) and the fluctuations around the estimates are consistent with the fact that the stress of the postural change is slowly increasing (Fig. 4A and 4B, middle graph). The ΔRRHR) decreases (increases) slowly through the transient period, is well below (above) the baseline of zero by the end of the transient period and stays well below (above) zero for the balance of this postural change. The ΔRRV and the ΔHRV (Fig. 4C and 4D, middle graph) show a gradual increase through the first part of the slow tilt with an increased level of fluctuations about the time course of each mean. At approximately thirty seconds into the transient period, both show an increase in level and fluctuations about this level. Compared with the preceding rest period, ΔRRV stabilizes to levels markedly below rest levels, whereas ΔHRV goes back to initial levels, with both measures showing greater fluctuations than during the rest period.

For the stand up postural change ΔRRHR) shows an abrupt change in the transient period of more than 300 s for the RR interval (20 beats per minute for heart rate) within fifteen seconds (Fig. 4A and 4B, bottom graph). Both indices return almost to baseline by thirty seconds, but then slowly decrease (increase) again for the balance of this postural maneuver. As observed for the rapid and slow tilt series, the variability about this level is increased both during the transient period and following the transient period. The ΔRRV increases in the transient period, then declining below rest level prior to the postural change (Fig. 4C, bottom graph). The fluctuations about this level, both during and after the transient period, are substantially greater than when the subjects were at rest. The ΔHRV change in the transient period mirrors the change reported for ΔHR (Fig. 4D, bottom graph). It shows a very abrupt increase reaching a maximum at around fifteen to twenty seconds into the transient period and declines to a level that is still above baseline by the end of the transient period. Also here, as observed for the other two maneuvers, ΔRRV stabilizes to levels markedly below rest level, whereas ΔHRV goes back to the initial level. This trend, common to all posture changes, confirms on one side ΔRRV as an index of autonomic activation as defined by standard HRV analysis (Task Force, 1996), and points at ΔHRV as a new measure proportional to the transient hemodynamic variations.

In summary, our results show that together our instantaneous HR and HRV estimates provide a different signature for the transient period for each of the three types of postural change (Fig. 4). First, the time courses of the mean R–R intervals and the heart rate series follow the well documented patterns of cardiovascular response to orthostatic stress (Borst et al., 1982; Ewing et al., 1980; Heldt et al., 2002, 2003; Sprangers et al., 2003; Tanaka et al., 1995). In contrast, our heart rate variability series revealed a new pattern of temporal dynamics. In particular, the sharp increase in heart rate standard deviation immediately after standing up, that persisted for at least 20 seconds before returning to baseline (Fig. 4D), may be associated not only with autonomic changes (Stauss, 2003), but also with heart rate modulation due to an increase in venous return as leg muscles contract when the upright posture is assumed (Heldt et al., 2002; Sprangers et al., 2003; Tanaka et al., 1995). Hence, our heart rate standard deviation index may offer a potentially non-invasive correlate of changes in cardiac output and/or total peripheral resistance (Sprangers et al., 2003; Taylor et al., 1995).

5. Conclusion

We have presented a point process algorithm to study heart beat dynamics as applied to recordings from healthy subjects during a Tilt table protocol. The algorithm is based on the inverse Gaussian probability model because this model can be derived as a stochastic version of the widely-applied deterministic integrate-and-fire models used to simulate heart beats. We have used an autoregressive model to define the mean of the inverse Gaussian and to represent the dependence of the RR interval length on the recent state of the autonomic inputs to the SA node. The point process adaptive filter was used to estimate the time-varying nature of the inputs to the SA node. Because the point process heart beat model is defined in continuous time, the algorithm computes estimates of mean RR, RR standard deviation, mean heart rate and heart rate standard deviation that can be updated at any desired temporal resolution. With changes in autonomic state, our analysis shows different behaviors between HR and HRV dynamics, as well as distinctive signatures reflecting the different hemodynamic changes elicited by three kinds of posture variation.

Future work will focus on more complex representations of the probability function, incorporating the point process framework into linear and nonlinear models of cardiovascular control and autonomic regulation, with inclusion of other cardiovascular variables such as arterial blood pressure, central venous pressure and respiration.


We are grateful to Roger G. Mark and Thomas Heldt, Harvard-Massachusetts Institute of Technology Division of Health Sciences and Technology, for kindly providing the tilt-table data analyzed in this study. This work was supported by NIH Grants R01-HL084502, R01-DA015644 and DP1-OD003646.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


  • Barbieri R, Brown EN. Analysis of heart dynamics by point process adaptive filtering. IEEE Trans Biomed Eng. 2006;53(1):4–12. [PubMed]
  • Barbieri R, Frank LM, Nguyen DP, Quirk MC, Solo V, Wilson MA, Brown EN. Dynamic analyses of information encoding in neural ensembles. Neural Computation. 2004;16:277–307. [PubMed]
  • Barbieri R, Matten EC, Alabi AA, Brown EN. A point process model of human heart rate intervals: new definitions of heart rate and heart rate variability. Am J Physiol Heart Circ Physiol. 2005;288:H424–H435. [PubMed]
  • Berger RD, Akselrod S, Gordon D, Cohen RJ. An efficient algorithm for spectral analysis of heart rate variability. IEEE Trans Biomed Eng. 1986;33:900–904. [PubMed]
  • Borst C, Wieling W, van Brederode JF, Hond A, de Rijk LG, Dunning AJ. Mechanisms of initial heart rate response to postural change. Am J Physiol Heart Circ Physiol. 1982;243:H676–H681. [PubMed]
  • Brown EN, Barbieri R, Ventura V, Kass RE, Frank LM. The time-rescaling theorem and its application to neural spike train data analysis. Neural Computation. 2002;14:325–346. [PubMed]
  • Brown EN, Barbieri R, Eden UT, Frank LM. Likelihood methods for neural data analysis. In: Feng J, editor. Computational Neuroscience: A Comprehensive Approach, ch. 9. London: Chapman and Hall; 2003. pp. 253–286.
  • Chhikara RS, Folks JL. The inverse Gaussian distribution: theory, methodology, and applications. New York: Dekker Inc.; 1989.
  • De Boer RW, Karemaker JM, Strackee J. Spectrum of a series of point event, generated by the integral pulse frequency modulation model. Med Biol Eng Comput. 1985;23:138–142. [PubMed]
  • Eden UT, Frank LM, Barbieri R, Solo V, Wilson MA, Brown EN. Dynamic analysis of neural encoding by point process adaptive filtering. Neural Computation. 2004;16:971–998. [PubMed]
  • Ewing DJ, Hume L, Campbell IW, Murray A, Neilson JM, Clarke BF. Autonomic mechanisms in the initial heart rate response to standing. J Appl Physiol. 1980;49:809–814. [PubMed]
  • Heldt T, Shim EB, Kamm RD, Mark RG. Computational modeling of cardiovascular response to orthostatic stress. J Appl Physiol. 2002;92:1239–1254. [PubMed]
  • Heldt T, Oefinger MB, Hoshiyama M, Mark RG. Circulatory Response to Passive and Active Changes in Posture. Computers in Cardiology. 2003;30:263–266.
  • Sprangers RL, Wesseling KH, Imholz AL, Imholz BP, Wieling W. Initial blood pressure fall on stand up and exercise explained by changes in total peripheral resistance. J Appl Physiol. 1991;70:523–530. [PubMed]
  • Stauss HM. Heart rate variability. American Journal of Physiology - Regulatory, Integrative and Comparative Physiology. 2003;285(5):R927–R931. [PubMed]
  • Tanaka H, Sjoberg BJ, Thulesius O. Cardiac output and blood pressure during active and passive standing. Clin Physiol. 1996;16(2):157–170. [PubMed]
  • Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Heart rate variability: standards of measurement, physiological interpretation and clinical use. Circulation. 1996;93:1043–1065. [PubMed]
  • Taylor JA, Hayano J, Seals DR. Lesser vagal withdrawal during isometric exercise with age. J Appl Physiol. 1995;79(3):805–811. [PubMed]
  • Tuckwell HC. Vol 2 Nonlinear and Stochastic Theories. New York: Cambridge University Press; 1988. Introduction to Theoretical Neurobiology.