HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Supplemental Document All Versions of this Article: 288/1/H424    most recent 00482.2003v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Barbieri, R. Articles by Brown, E. N. Search for Related Content PubMed PubMed Citation Articles by Barbieri, R. Articles by Brown, E. N.

A point-process model of human heartbeat intervals: new definitions of heart rate and heart rate variability

Neuroscience Statistics Research Laboratory, Department of Anesthesia and Critical Care, Massachusetts General Hospital, and Division of Health Sciences and Technology, Harvard Medical School/Massachusetts Institute of Technology, Boston, Massachusetts

Submitted 23 May 2003 ; accepted in final form 12 September 2004

	ABSTRACT

TOP ABSTRACT THEORY AND METHODS RESULTS DISCUSSION GRANTS REFERENCES

 
Heart rate is a vital sign, whereas heart rate variability is an important quantitative measure of cardiovascular regulation by the autonomic nervous system. Although the design of algorithms to compute heart rate and assess heart rate variability is an active area of research, none of the approaches considers the natural point-process structure of human heartbeats, and none gives instantaneous estimates of heart rate variability. We model the stochastic structure of heartbeat intervals as a history-dependent inverse Gaussian process and derive from it an explicit probability density that gives new definitions of heart rate and heart rate variability: instantaneous R-R interval and heart rate standard deviations. We estimate the time-varying parameters of the inverse Gaussian model by local maximum likelihood and assess model goodness-of-fit by Kolmogorov-Smirnov tests based on the time-rescaling theorem. We illustrate our new definitions in an analysis of human heartbeat intervals from 10 healthy subjects undergoing a tilt-table experiment. Although several studies have identified deterministic, nonlinear dynamical features in human heartbeat intervals, our analysis shows that a highly accurate description of these series at rest and in extreme physiological conditions may be given by an elementary, physiologically based, stochastic model.

tilt table; inverse Gaussian; time-rescaling theorem; Kolmogorov-Smirnov test; autonomic regulation

Heart rate is traditionally estimated as the number of R-wave events (heartbeats) per unit time on the electrocardiogram (ECG) or as the average of the R-R interval reciprocals within a specified time window. Several approaches are used to characterize heart rate variability: elementary statistical measures of R-R interval properties (51), spectral analysis of heart rate or R-R interval time series (1, 3–5, 12, 38, 44, 51), deterministic dynamical systems assessments of heart rate signal properties (28, 29, 41, 49), and approximate entropy measures of R-R interval regularity (40).

Three issues in these approaches used to characterize heart rate variability have yet to be addressed. 1) All methods treat the R-R interval or heart rate series as continuous-valued signals, rather than model them to reflect the point-process structure of the underlying R-wave events. That is, the R-wave events, which mark the electrical impulses from the heart's conduction system that represent ventricular contractions, are a sequence of discrete occurrences in continuous time and, as such, form a point process. 2) None of the approaches assesses model goodness-of-fit, i.e., how well a given model describes the observed R-R interval series as part of its standard analysis framework. 3) None of the approaches provides instantaneous estimates of heart rate variance and R-R interval variance to accompany instantaneous estimates of heart rate and mean R-R interval (heart period). The Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology (51) stated that "the present possibilities of characterizing and quantifying the dynamics of the R-R interval sequence and transient changes in [heart rate variability] are sparse and still under mathematical development." It predicted that "proper assessment of [heart rate variability] dynamics will lead to substantial improvement in our understanding of both the modulations of heart period and their physiological and pathophysiological correlates."

To address these three issues, we model human heartbeat intervals as a history-dependent inverse Gaussian (HDIG) point process and derive from it an explicit probability density that yields new definitions of heart rate and heart rate variability. We estimate the time-varying parameters of the inverse Gaussian model by local maximum likelihood, assess model goodness-of-fit by Kolmogorov-Smirnov (KS) tests derived from the time-rescaling theorem, and use a local likelihood version of Akaike's information criterion (AIC) to guide model selection. We compare our new definitions of heart rate and heart rate variability with current methods in a study of ECG recordings from 10 healthy subjects during a tilt-table experiment.

	THEORY AND METHODS

TOP ABSTRACT THEORY AND METHODS RESULTS DISCUSSION GRANTS REFERENCES

 
Physiology of R-wave events. Our objective is to arrive at precise probabilistic definitions of heart rate and heart rate variability. We begin the derivation of our heartbeat probability model by reviewing the physiology of the R-wave events. We assume that, in an observation interval (0,T], we record 0 < u₁ < u₂ <, ..., < u_k <, ..., < u_K T, K successive R-wave event times from an ECG. As stated above, the R-R intervals are the times between two successive R-wave events.

Each R-wave event is initiated by the synchronous depolarization of the heart's pacemaker cells beginning in the sinoatrial (SA) node of the right atrium and then propagating through its specialized conduction system to the left atrium and to the two ventricles (22). In a manner similar to that whereby a nerve cell generates an action potential, the synchronous depolarization is accomplished everywhere along the conduction. After every depolarization, the transmembrane potentials of these cells return to their resting potentials and begin anew their spontaneous rise toward threshold (22).

The start of this rise marks the start of the waiting time until the next depolarization. The rise of the membrane potential to threshold has been modeled as a Gaussian random walk with drift (47). The probability density of the first passage times for this random-walk process, i.e., the times between threshold crossings (R-R intervals), is well known to be the inverse Gaussian (47, 55). If, as a first approximation, the generation of normal heartbeats (R-wave events) is assumed to be unrelated to any previous or future heartbeat and the membrane potential rise to threshold is modeled as a Gaussian random walk with drift, then the series of elapsed times between beats (the R-R intervals) can be modeled as independent, inverse Gaussian random variables (43).

Independence of the R-R intervals would be a reasonable assumption if the effects of sympathetic and parasympathetic inputs from the autonomic nervous system to the SA node were absent or short-lived. However, although these inputs can occur in milliseconds, their effects can last for several seconds. Consequently, an isolated increase in sympathetic input to the SA node may produce shorter R-R intervals, whereas an isolated increase in parasympathetic input may produce longer R-R intervals, for a time. The duration of the increase or decrease in the R-R intervals depends on the cause of the change in parasympathetic and/or sympathetic input, as well as the response of the other components of the cardiovascular control circuitry. In any case, because the effects of these inputs to the SA node persist for several seconds, R-R interval lengths are not independent but, rather, are a function of the recent history of these inputs. Furthermore, because the sympathetic and parasympathetic inputs to the SA node are part of the cardiovascular control circuitry, these inputs are dynamic. That is, they are always changing and, thus, are continuously altering the propensity of the SA node to initiate heartbeats tomaintain homeostasis in physiological and pathological conditions. Therefore, given the effect of the recent history of autonomic inputs to the SA node, our model must take into account the dynamic character of these inputs.

Heart rate probability model. We assume that, given any R-wave event u_k, the waiting time until the next R-wave event or, equivalently, the length of the next R-R interval, obeys an HDIG probability density f(t|H_uk,), where t is any time satisfying t > u_k, H_uk is the history of the R-R intervals up to u_k, and is a vector of model parameters. The model is defined as

(1)

where H_uk = (u_k,w_k,w_k–1,...,w_k–p+1), w_k = u_k – u_{k – 1} is the kth R-R interval, µ(H_uk) = ₀ + _j=1^p _jw_k–j+1 > 0 is the mean, _{p + 1} > 0 is the scale parameter, and = (₀,₁,...,_{p + 1}). This model represents the dependence of the R-R interval length on the recent history of parasympathetic and sympathetic inputs to the SA node by modeling the mean as a linear function of the last p R-R intervals. If we assume that the R-R intervals are independent (i.e., p = 0), then µ(H_uk,) = ₀, f(t|H_uk,) = f(t|u_k,₀,₁), and Eq. 1 simplifies to a renewal inverse Gaussian (RIG) model. The mean and standard deviation of the R-R probability model in Eq. 1 are

(2)

(3)

respectively.

Because our probability density in Eq. 1 characterizes the stochastic properties of the R-R intervals, we use it to formulate precise definitions of heart rate and heart rate variability. Heart rate is often defined as the reciprocal of the R-R intervals (51). Hence, for any t > u_k, t – u_k is the waiting time until the next R-wave event, and we define r = c(t – u_k)^–1 as the heart rate random variable, where c = 6 x 10⁴ ms/min is the constant that converts the R-R interval measurements recorded in milliseconds to heart rate measurements in beats per minute. Therefore, because r is a one-to-one transformation of t – u_k, we use the standard change-of-variables formula from elementary probability theory (43) and derive from the R-R interval probability density in Eq. 1, f(r|H_uk,), the heart rate probability density defined as

(4)

where µ*(H_uk,) = c^–1µ(H_uk,) and *_{p + 1} = c^–1_{p + 1}. The mean, mode, and standard deviation of the heart rate (HR) probability density are, respectively,

(5)

(6)

(7)

Equation 4 defines the stochastic properties of heart rate in terms of a probability density. If p = 0, then Eq. 4 gives the simple reciprocal of the inverse Gaussian probability density (13). To use Eq. 4 in analyses of R-R interval time series, heart rate should be a representative value from the heart rate probability density. Therefore, we define heart rate as the mean (Eq. 5) or the mode (Eq. 6) of f(r|H_uk,), and the heart rate variance as the square of the standard deviation (Eq. 7). Equations 3 and 7 are important steps in constructing the estimates of heart rate variability, because they make explicit the relation between R-R interval variance and heart rate variance. We will use Eqs. 3 and 7 to construct our indexes of heart rate variability in Local maximum-likelihood estimation of instantaneous heart rate and heart rate variability.

To summarize, f(t|H_uk,) (Eq. 1) is a probability model for the R-R intervals that represents the effects of recent sympathetic and parasympathetic input to the SA node (history dependence) through its mean parameter, µ(H_uk,). Through the standard change of variables, the R-R interval probability model is transformed into the heart rate probability density, f(r|H_uk,) (Eq. 4). To account for the continuous dynamics of the sympathetic and parasympathetic inputs to the SA node, we assume that the model parameter is time varying. Hence, by tracking the time course of and continuously evaluating Eqs. 5 and 7 as evolves, we can track instantaneous heart rate and heart rate standard deviation, respectively.

Local maximum-likelihood estimation of instantaneous heart rate and heart rate variability. We use a local maximum-likelihood procedure to estimate the time-varying parameter (34, 53). For the ECG recording period (0,T], let l be the length of the local likelihood observation interval for t [l, T], and let define how much the local likelihood time interval is shifted to compute the next parameter update. Let t^l be the local likelihood interval (t – l,t], and assume that, within t^l, we observe n_t R-wave event times t – l < u₁ < u₂ <,...,< u_nt t. We let u_{t – l:t} = (u₁,...,u_nt). If is time varying, then at time t we define the local maximum-likelihood estimate _t to be the maximum-likelihood estimate of on t^l.

To compute the local maximum-likelihood estimate of , we first define the local joint probability density of u_{t – l:t}. If we condition on the p R-R intervals preceding each R-wave event in u_{t – l:t}, then the local observation interval t^l induces right censoring of the R-R interval measurements, because the n_t + 1st interval is not completely observed. If we take into account the right censoring, the local log likelihood is

(8)

where w(t) is a weighting function for the local likelihood estimation (34). We chose w(t – u) = e^{–(t – u)}, where is a weighting time constant that governs the degree of influence of a previous observation u on the local likelihood at time t. Increasing decreases the influence of a previous observation on the local likelihood and vice versa. Our method for choosing the value of is described in RESULTS.

We use a Newton-Raphson procedure to maximize the local log likelihood in Eq. 8 and compute the local maximum-likelihood estimate of _t (34). Because there is significant overlap between adjacent local likelihood intervals, we start the Newton-Raphson procedure at t with the previous local maximum-likelihood estimate at time t – . Once _t is computed, the interval t^l is shifted to (t – l + ,t + ], and the local maximum-likelihood estimation is repeated. The procedure is continued until t = T.

Given _t, the local maximum-likelihood estimate of at time t, it follows from Eqs. 5 and 7 that the instantaneous estimates of mean R-R, R-R interval standard deviation, mean heart rate, and heart rate standard deviation at time t are

(9)

(10)

(11)

(12)

, respectively.

A key feature of our analysis framework is that we can estimate mean R-R interval (Eq. 9), R-R interval standard deviation (Eq. 10), heart rate (Eq. 11), and heart rate standard deviation (Eq. 12) in continuous time. This is because the HDIG model is defined in continuous time; hence, the heart rate probability model is as well. For any t (l,T], the local likelihood procedure uses all data in the local likelihood interval (t – l,t] to compute _t. If t = u_nt, the contribution to the local likelihood is log f(u_nt – u_{nt –1}|Hunt,_t), and if t > u_nt and the last interval is censored, then the contribution to the local likelihood is log _t–untf(v|H_unt,_t)dv. Whether the last interval is completely observed or censored, it is used in the estimation. This means that the local likelihood function is defined at t, the right end point of the observation interval and the point where the local maximum-likelihood estimate is to be computed. Therefore, because can be estimated in continuous time, all the indexes in Eqs. 9–12 can be computed in continuous time, because all are continuous functions of _t. Equations 9 and 11 are our estimates of the instantaneous mean R-R interval and heart rate, respectively. Equations 10 and 12 are our estimates of instantaneous R-R interval standard deviation and heart rate standard deviation, respectively. We designate Eqs. 10 and 12 our indexes of heart rate variability.

Assessing model goodness-of-fit and choosing the model order and local likelihood interval length. Our R-R interval probability model, together with the local maximum-likelihood method, provides an approach for estimating instantaneous mean R-R interval, heart rate, R-R interval standard deviation, and heart rate standard deviation from a time series of R-R intervals. For this approach, it is also necessary to evaluate model goodness-of-fit, i.e., to determine how well the model describes the series of R-R intervals. Because Eq. 1 defines an explicit point-process model, we can use the time-rescaling theorem (2, 10, 37) to evaluate model goodness-of-fit. To do this, we compute the time-rescaled, or transformed, R-R intervals defined by

(13)

where the u_k values (the R-wave event times) represent a point process observed in (0,T) and (t|H_t,_t) is the conditional intensity function defined as

(14)

evaluated at the local maximum-likelihood estimate _t. The conditional intensity is the history-dependent rate function for a point process that generalizes the rate function for a Poisson process. By the time-rescaling theorem, the _k values are independent, exponential random variables with a unit rate (2, 10, 37). Under the further transformation z_k = 1 – exp(–_k), the z_k values are independent, uniform random variables on the interval (0,1]. Therefore, we can construct a KS test to assess agreement between the z_k values and the uniform probability density (2, 10, 37). Because the transformation from the u_k values to the z_k values is one-to-one, close agreement between the z_k values and the uniform probability density is true if, and only if, there is close agreement between the point-process probability model and the series of R-R intervals. Hence, the time-rescaling theorem provides a direct means of measuring agreement between a point-process probability model and an R-R interval series.

We compute the KS test in terms of the KS plot or the KS distance. The KS plots are constructed and analyzed as described previously (2, 10, 37). The KS distance measures the largest distance between the cumulative distribution function of the R-R intervals transformed to the interval (0,1] and the cumulative distribution function of a uniform distribution on (0,1]. The smaller the KS distance, the closer is the agreement between the original heartbeat interval time series and the proposed model.

Under the time-rescaling theorem, the R-R intervals should be transformed to _k values, which are independent, regardless of the degree of dependence in the original R-R intervals. Although it is difficult to measure independence in a time series, if the _k values are independent, then they should at least be uncorrelated. Therefore, as a further assessment of model goodness-of-fit, we plot _k vs. _{k – 1} for each model. The more unstructured the plot of _k vs. _{k – 1}, the less the correlation and, hence, the more consistent the _k values are with being independent. To assess further the correlation structure and, hence, possible dependence in the _k values that may be present beyond first order, we plot the serial correlation function of the _k values for each subject for 60 lags (1 min). Small values of the serial correlation function at all lags would suggest that the _k values are uncorrelated and would be consistent with their being independent. Approximate independence of these transformed intervals would suggest that the proposed model was highly consistent with the heartbeat interval series.

To conduct our local maximum-likelihood analysis, we set p, the order of the history dependence in Eq. 1, and l, the length of the local likelihood interval. To help choose these two parameters, we computed the approximate AIC for local maximum-likelihood analyses (34).

Tilt-table protocol. We illustrate our methods in a study of the R-R interval time series recorded from 10 healthy subjects in whom cardiovascular and autonomic regulation were studied using a tilt-table protocol (Fig. 1) (25). In a tilt-table study, a subject is first placed horizontally in a supine position, with restraints used to secure him/her at the waist, arms, and hands. The legs are not secured. The subject is then tilted from the horizontal to the vertical position and returned to the horizontal position. The wide range of changes in autonomic input to the heart induced by varying the speed and duration of the postural changes induces a wide range of changes in the R-R intervals and, thereby, makes the tilt-table protocol an excellent paradigm for testing the ability of our new algorithm to track R-R interval and heart rate dynamics.

View larger version (43K): [in this window] [in a new window]   Fig. 1. R-R interval time series from subjects 1–10 in a tilt-table study plotted as a function of time at which each interval was initiated. Study consisted of an initial 5 min of baseline recordings and 3 types of postural changes: rapid up-down tilt, slow up-down tilt, and standing upright and lying supine. Each postural change was repeated twice. Sequence of tilts was performed in random order for each subject. #,*Subjects whose analyses are described in RESULTS (Figs. 5 and 6, respectively).

	RESULTS

TOP ABSTRACT THEORY AND METHODS RESULTS DISCUSSION GRANTS REFERENCES

 
We performed our analysis in four parts. 1) We determined which HDIG models best described the heartbeat interval time series of the 10 tilt-table subjects in function of the order of the autoregression p, the local likelihood interval l, and the weighting coefficient . 2) We compared the fits of the best HDIG models with those of a local average (LA) model and an RIG model using goodness-of-fit assessments based on the time-rescaling theorem. 3) We illustrated the use of our indexes of heart rate variability to characterize the dynamic properties of heartbeat intervals in the 10 tilt-table study subjects. 4) We compared standard time domain measures of heart rate variability with our indexes of heart rate variability, and we compared spectral measures of heart rate variability computed from standard interpolated heart rate estimates with these spectral measures computed from our HDIG estimates of heart rate.

The R-R interval series for the 10 subjects during the tilt protocols showed a broad range of patterns (Fig. 1). The R-R intervals shortened with tilt periods and lengthened with rest periods, making these periods readily apparent in the R-R interval series of subjects 1, 2, 5, 7, and 10. In contrast, subjects 3, 4, 6, and 8 showed minimal changes in R-R intervals across all tilt and rest periods. Subjects 2 and 9 had a wide range of R-R interval lengths in the tilt and rest periods. Although in the heartbeat interval series of subject 2 the tilt and rest periods were apparent in the R-R interval series, they were not apparent in the R-R interval series of subject 9.

Model selection and model goodness-of-fit. To help choose a model order p, a local likelihood interval l, and the weighting function time constant , we fit the HDIG model with orders p = 2, 4, 6, 8, 10, and 12, l = 15–180 s in increments of 15 s, and weighting coefficients = 0, 0.01, 0.02, 0.05, and 0.1 to the R-R interval time series for all 10 subjects with the shift parameter = 5 ms. The AIC and KS plots from this preliminary analysis suggested that p = 2 for subjects 2 and 7, p = 4 for the remaining eight subjects, l = 60 s, and = 0.01.

For each subject, we compared the fit of the HDIG model with the fit of the RIG model (p = 0, l = 60) and also with heart rate estimates computed as the average of the inverse of the R-R intervals in the same 60-s local likelihood interval. We designated the latter the LA model.

To compare for each subject the fit of the best HDIG model with the RIG and LA models, we applied the KS test, based on the time-rescaling theorem, using the conditional intensity functions estimated for each of the three models (Table 1). For all 10 subjects, the LA model was in poor agreement with the data in terms of KS distance (Table 1). For all subjects, the RIG model showed a substantial improvement (smaller KS distance) compared with the LA model, and for subjects 7–9, the RIG model KS distances were close to those of the best HDIG fits. The HDIG model had the smallest KS distances for all 10 subjects. For subjects 3–5, the HDIG fits were nearly perfect, in that the KS distances were within the 95% confidence bounds (Fig. 2).

View larger version (26K): [in this window] [in a new window]   Fig. 2. Kolmogorov-Smirnov (KS) plots of time-rescaled R-R intervals (Eq. 13) derived from history-dependent inverse Gaussian (HDIG, bold solid line), renewal inverse Gaussian (RIG, dashed line), and local average (LA, dotted line) model fits for 4 subjects. A proposed model was considered consistent with R-R interval data if its KS plot fell entirely within 95% confidence bounds, defined by parallel thin solid diagonal lines in A–D. A: best HDIG model fit; KS plot of HDIG model fit entirely within 95% confidence bounds. B: best improvement of an HDIG model fit relative to an RIG model fit (dashed line). C: best RIG model fit. D: worst HDIG model fit. KS plot of subject 3 resembled plots of subjects 4 (A) and 5 (B); KS plots of subjects 1, 6, 8, 9, and 10 resembled plot of subject 7 (C).

In contrast, the KS plots of the HDIG model (Fig. 2, bold solid curves) were entirely within the 95% confidence bounds for subjects 4 and 5 (Fig. 2B), just above the upper 95% bounds for the 0.50–0.70 quantile for subject 7, and below the 95% confidence bounds for lower quantiles for subject 2 (Fig. 2D). The KS plot of subject 3 resembled those of subjects 4 and 5, whereas those of subjects 1, 6, 8, 9, and 10 resembled that of subject 7.

By the time-rescaling theorem, the KS plot analyses suggest that the _k values, transformed heartbeat intervals from the HDIG model, are consistent with an exponential probability density with mean waiting time of 1 or equivalently, the z_k values are consistent with a uniform probability density on the interval (0,1] (see THEORY AND METHODS). However, in addition, if this model describes the data accurately, then by the time-rescaling theorem, these transformed intervals should also be independent. Hence, a plot of _k vs. _{k – 1} for the best-fitting model should be uncorrelated, and its scatterplot should resemble a random point cloud. The results of the correlation analysis of _k values for subject 10, which are identical to those from subjects 1–9, are shown in Fig. 3. The adjacent R-R intervals are highly correlated, with a correlation coefficient of 0.96 (Fig. 3A). The plot of _k vs. _{k – 1} for the LA model (Fig. 3B) has a persistently high correlation of 0.89. It is not surprising that the LA model's scatterplot closely resembles the original R-R interval plot, inasmuch as this model estimates the conditional intensity function (Eq. 12) as a simple, locally constant function. The _k values from the RIG model (Fig. 3C) also remain highly correlated, with a correlation coefficient of 0.70, despite improved KS plots relative to the LA model. For the HDIG model, however, the plot of _k vs. _{k – 1} is an almost perfect point cloud (Fig. 3D), with a correlation coefficient of 0.05.

View larger version (26K): [in this window] [in a new window]   Fig. 3. Scatter-plot analysis of time-rescaled R-R intervals for subject 10. A: plot of R-Rk vs. R-Rk – 1 shows strong correlation of 0.96. Correlation of time-rescaled R-R intervals, k values, was 0.89 for LA model (B), 0.70 for RIG model (C), and 0.05 for HDIG model (D). Inset: enlargement of scatter plot in D showing strongly random structure in transformed times for HDIG model and suggesting that they were uncorrelated.

View larger version (24K): [in this window] [in a new window]   Fig. 4. Correlation functions of time-rescaled R-R intervals, k values, from HDIG model fit for lags 1–60 for subjects 1–10. Horizontal parallel lines are approximate 95% confidence bounds computed as ±2(Kj) (8), where Kj is the number of R-R intervals for subject j, for j = 1,...,10. Correlations outside these bounds were statistically significant. Except for large correlations for subjects 1, 4, 6, and 8, correlations are within confidence bounds, suggesting that k values were approximately uncorrelated for each subject. This lack of correlation was consistent with k values being independent and suggested, by the time-rescaling theorem, that the HDIG model agreed closely with the R-R interval series of each subject.

View larger version (36K): [in this window] [in a new window]   Fig. 5. Analysis of instantaneous heart rate and heart rate variability for subject 2 for the entire tilt-table study. A: reciprocal R-R intervals. B: HDIG model instantaneous heart rate (HR) series estimate. C: LA/RIG model instantaneous heart rate series estimate. D: HDIG model instantaneous R-R interval standard deviation (sd) estimate. E: HDIG instantaneous heart rate standard deviation estimate. HDIG heart rate estimates closely resembled reciprocal of R-R intervals. LA and RIG heart rate estimates were identical, although these approaches differed in how they model the R-R interval series. R-R interval and heart rate standard deviations show that there can be low and high variability at rest and in upright posture. bpm, Beats per minute.

The instantaneous R-R interval standard deviation (Fig. 5D) and heart rate standard deviation (Fig. 5E) estimates showed dynamics different from those shown by the instantaneous heart rate estimates. The R-R interval standard deviation estimate showed that the variation in R-R intervals can be high or low at rest and in the upright posture. Similarly, the heart rate standard deviation showed that the heart rate variance can also be high or low at rest and in the upright posture. The heart rate standard deviation showed finer fluctuations on a second-to-second basis than the R-R interval standard deviation. Because the instantaneous R-R standard deviation and the instantaneous heart rate standard deviation are second-moment estimates, they were, as expected, sensitive to sharp changes in the R-R series length, such as occurred with an abrupt change due to a very short or long interval from an ectopic beat (Fig. 5, D and E). Whereas the effect of an ectopic beat on the instantaneous heart rate estimate was short-lived, its effect on the R-R interval and the heart rate standard deviations lasted for a longer period. These results show that it is possible to accompany an instantaneous estimate of heart rate with an instantaneous estimate of heart rate variability computed as instantaneous R-R interval or heart rate standard deviations.

A detailed summary of the subjects is given in supplemental data for this article, which may be found at http://ajpheart.physiology.org/cgi/content/full/00483.2003/DC1.

For each postural change, our four measures of heart rate dynamics gave different pattern changes. To illustrate these, we show the instantaneous mean R-R interval, instantaneous heart rate, and instantaneous R-R interval and heart rate standard deviations from subjects 1, 4, and 5 during the 200 s before and after a rapid tilt (Fig. 6A), a slow tilt (Fig. 6B), and an upright position (Fig. 6C), respectively. Before all three postural changes, all four measures were approximately constant. At the start of the rapid tilt, the heart rate of subject 1 increased immediately and, after a small transient decline, continued to increase, whereas the mean R-R interval decreased immediately and, after a small transient increase, continued to decrease (Fig. 6A). With the rapid tilt, the R-R interval and heart rate standard deviations increased and then declined to near-constant values for the balance of the tilt.

View larger version (29K): [in this window] [in a new window]   Fig. 6. HDIG model analysis of heart rate and heart rate variability. A: 1st rapid up-down tilt at 880–1,550 s from subject 1. B: 1st slow up-down tilt at 2,180–2,900 s from subject 4. C: 2nd upright standing-to-supine postural change at 2,471–3,080 s from subject 5. 1, R-R intervals; 2, HDIG model instantaneous heart rate series estimate; 3, HDIG model instantaneous R-R interval standard deviation estimate; 4, HDIG instantaneous heart rate standard deviation estimate. Vertical dashed lines demarcate onset and offset of postural change. Double vertical lines in B demarcate transitions from horizontal (vertical) to vertical (horizontal) for slow up-down tilt.

When subject 5 stood upright, heart rate (mean R-R interval) increased (decreased) rapidly and continued to fluctuate while the subject was in the upright posture (Fig. 6C). In a manner analogous to that observed for the rapid tilt (Fig. 6A), the R-R interval standard deviation increased initially and then declined. The heart rate standard deviation also rose very rapidly when the subject stood upright and declined to a constant level within 50 s. This pattern of dynamics in the heart rate standard deviation closely resembled the pattern of changes in cardiac output observed when a subject stands upright (46, 50).

For each of these three positional changes, the heart rate and mean R-R interval returned to their values before the positional change, yet the R-R interval standard deviations were above their values before the positional change, suggesting that the physiological states of the subjects after the positional changes were different from those before the positional changes. A summary analysis of all the subjects for the rapid tilts, slow tilts, and upright-standing intervals is given in the supplemental data for this article (see above). The patterns shown here in the dynamics of these four measures were also observed in the other subjects (see supplemental data, Fig. S1).

Comparison with standard heart rate and heart rate variability indexes. We have focused on the use of our paradigm to derive instantaneous measures of heart rate and heart rate variability. Here we compare our indexes of heart rate variability with a standard time domain measure, the standard deviation of the R-R intervals (SDNN). We also compared estimates of the low-frequency (LF)-to-high-frequency (HF) index (LF/HF) computed from our heart rate series with this index computed from interpolated heart rate series (3). For these analyses, we used, from each subject, 12 (6 rest and 6 upright) 120-s epochs of his/her heartbeat interval series. For the value of each index during the rest (upright) period, the rest (upright) average for the subject was computed in the 120-s epoch and then averaged across the 10 subjects.

In the time domain analysis, the mean heart rate and mean R-R interval computed by simple averages were identical to the average values derived from our HDIG model (Table 2). In particular, both showed the same magnitude of increase in heart rate and decrease in R-R interval with rest compared with upright posture. We computed the SDNN as the standard deviation of the R-R interval series in each 120-s epoch (Table 2). On average, during rest, the SDNN was 57 ms and decreased to 49 ms in the upright posture, whereas our R-R interval standard deviation was, on average, 37 ms and decreased to an average of 27 ms in the upright posture.

Higher values of LF/HF (Table 2) were associated with the upright epochs for the standard and new heart rate variability indexes. This is consistent with an increase in the sympathovagal balance. LF/HF computed from our new indexes increased from 2.98 to 3.90, whereas the standard index increased from 2.28 to 5.14 after posture change.

The reasons for these differences were apparent on examining the spectra from a rest and an upright period such as those from subject 5 (Fig. 7). The spectra computed from the HDIG heart rate estimates were significantly different from those computed from the interpolated series (Fig. 7, A and C), as assessed by KS plot analysis (Fig. 7E). There are two reasons for these differences. 1) The spectra computed from the HDIG heart rate series had less power in LF and HF than the standard spectra. During the rest epochs, LF/HF was close to 1 in both cases (0.87 and 0.126, respectively). On the other hand, in the upright posture, there was, for the HDIG and interpolated series, an increase in LF power and a sharp decrease in HF power (Fig. 7B). With change of posture, LF/HF for the HDIG heart rate series and standard series, increased (3.23 and 3.87, respectively). However, the lower LF power content of the HDIG heart rate series led to a smaller increase in LF/HF than in the standard series. Of the 120-KS plot analyses [10 subjects x (6 rest + 6 tilt)] comparing spectra from standard interpolated heart rate series and our HDIG heart rate series, 117 showed the differences between the spectra to be statistically significant (P < 0.05). 2) The spectra of the HDIG heart rate time series had power at frequencies higher than those seen in spectra computed from the heart rate series after interpolation. The heart rate series interpolated from the original heartbeat interval series can contain no frequencies beyond the approximate cutoff frequency defined by the mean R-R interval, i.e., _cutoff 0.5(mean R-R)^–1 (Fig. 7, B and D). For example, if the R-R interval is 800 ms, then _cutoff = 0.62 Hz. In contrast, our analysis framework represents the heartbeat interval series in continuous time using point-process models; hence, it places no restriction on the spectral content of the heartbeat interval series.

View larger version (17K): [in this window] [in a new window]   Fig. 7. Spectral analysis of rest (120 s before upright posture to start of upright posture) and upright posture (60 s after upright posture to 180 s after upright posture) segments from 1st slow tilt from subject 5. A–D: spectrum of rest [0–0.5 Hz (A) and 0.5–1.5 Hz (B)] and upright posture [0–0.5 Hz (C) and 0.5–1.5 Hz (D)] segments of interpolated reciprocal R-R beat series (dashed curve) and HDIG instantaneous heart rate estimates (solid curve). Low-frequency (LF) range is 0.04–0.15 Hz, high-frequency (HF) range is 0.15–0.40 Hz, and very-high-frequency (VHF) range is 0.5–1.5 Hz. At rest (upright), low- to high-frequency ratio (LF/HF) was 1.26 (3.87) for interpolated heart rate series and 0.87 (3.23) for HDIG heart rate series. E: KS plots comparing spectral power distribution of HDIG heart rate estimate (x-axis) and interpolated reciprocal R-R beat series (y-axis) for rest and tilt segments. Because both KS plots fell outside 95% confidence bounds, the 2 spectra differ when the subject is at rest and in the upright posture (21).

	DISCUSSION

TOP ABSTRACT THEORY AND METHODS RESULTS DISCUSSION GRANTS REFERENCES

 
To study heart rate and heart rate variability, we have developed a statistical model of human heartbeat intervals built on three statistical and physiological properties of these intervals: 1) a description of the point-process nature of the R-R intervals, 2) the dependence of the interval length on the recent state of autonomic inputs to the SA node, and 3) the dynamic character of these inputs. We represented these characteristics using an HDIG model to describe the point-process structure of R-R intervals and used local maximum likelihood to estimate the model's time-varying parameters. The inverse Gaussian model describes the first passage to threshold of the membrane voltages of the heart's pacemaker cells, the model's autoregressive structure describes the dependence of the R-R interval lengths on the recent history of the autonomic inputs to the SA node, and the model's time-varying parameters capture the dynamic character of these SA node inputs.

From the HDIG model, we derived an explicit probability density for heart rate and formulated precise definitions for mean R-R interval (Eq. 2), R-R interval standard deviation (Eq. 3), heart rate (Eq. 5), and heart rate standard deviation (Eq. 7), as well as their interrelations. We propose the instantaneous estimate of heart rate (Eq. 11) as the new definition of heart rate and the instantaneous estimates of R-R interval standard deviation (Eq. 10) and instantaneous heart rate standard deviation (Eq. 12) as new estimates of heart rate variability. All these measures of heartbeat dynamics are computed from the same probability framework. This analysis, which to our knowledge is the first report of instantaneous estimates of heart rate standard deviation and instantaneous estimates of R-R interval standard deviation to accompany instantaneous estimates of heart rate, provides an approach to constructing the type of dynamic algorithms for tracking heart rate variability called for by the Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology (51).

Our goodness-of-fit and model selection analysis showed that an order 2 or 4 HDIG model and a local maximum-likelihood interval of 60 s gave accurate descriptions of the heartbeat interval series of 10 subjects in a tilt-table study during multiple transitions between rest and tilt. Our analysis showed that the HDIG model gave a substantially better fit than the commonly used LA and RIG models. Although several studies have found nonlinear dynamical characteristics in human heartbeat intervals, our analysis provides compelling evidence that these data may be accurately described as a stochastic process with nonstationary parameters during multiple transitions between rest and extreme physiological conditions.

Time domain (51), frequency domain (1, 3, 6, 12, 44, 51), dynamical systems (28, 29, 41, 49), and entropy methods (40) for heart rate variability analysis have been previously reported. Most of these methods must convert R-R interval data to continuous-valued, evenly spaced measurements for analysis by interpolating the R-R intervals or the reciprocal of the R-R intervals. To apply the frequency domain methods appropriately, data must be collected under stationary conditions. Because of the requirement for stationarity, it is more challenging to track changes in the temporal dynamics of heartbeat intervals. The several methods designed to address nonstationarity apply time-varying methods to R-R interval (24, 30, 31, 36) or heart rate (4, 54) time series obtained by interpolating the R-R interval or the reciprocal of the R-R interval series, respectively. On the other hand, methods that do not interpolate typically do not consider updates at temporal resolutions smaller than two consecutive R-R intervals (5, 39, 56). These methods include time-varying spectral analysis (4, 5, 30, 36), wavelet analysis (31, 39, 54), complex demodulation (24), and coarse-graining analysis (56).

Our approach, based on the point-process method for neural spike train data analysis (2, 9, 10), begins an analysis with 1 min of data and computes instantaneous estimates of the heart rate variability, because it can compute sequential updates at any desired time resolution. The key feature of our model construction that makes this computation possible is definition of the point process in continuous time. Because nonstationarity is a physiological property of the R-R interval series, we model this feature as a stochastic process with time-varying parameters. We obviate the need for interpolation, because our heart rate probability model (Eq. 4) defines a precise probabilistic prescription for conversion of the point-process R-R intervals to the continuous-valued functions for heart rate (Eq. 5) and heart rate standard deviation (Eq. 7) in continuous time. Because we use local maximum-likelihood estimation to compute time-varying estimates, our analysis is equally valid for stationary and nonstationary conditions.

Furthermore, because the conditional intensity function gives a canonical characterization of a point process, we use the estimated conditional intensity function to construct goodness-of-fit tests by the time-rescaling theorem. Our goodness-of-fit tests allow us to compare different point-process models of R-R intervals and to help select the width of the local likelihood estimation interval as well as the order of the history dependence. Although we have not systematically studied other waiting models, analyses using log normal and gamma renewal models gave results similar to those from the RIG model. Although other methods give important characterizations of human heartbeat interval dynamics, our approach models the R-wave events as point processes, reports an overall goodness-of-fit assessment as a standard output, and computes instantaneous mean R-R, R-R standard deviation, mean heart rate, and heart rate standard deviation in one framework.

Despite the apparent similarities between our R-R interval model in Eq. 1 and a standard stationary autoregression model, there are important differences. A stationary autoregressive model of order p is typically a model of a continuous-valued Gaussian process sampled at equally spaced time points (8). Because the underlying process is assumed to be stationary, the model coefficients and the error variance are constant in time. In contrast to the Gaussian process in a stationary autoregressive model, which is a continuous-valued process sampled in discrete time, a point process is a discrete-valued process that occurs in continuous time (9, 14). The autoregression for our point-process model in Eq. 1 is constructed by assuming that the R-R intervals obey an inverse Gaussian model and that the mean of the current R-R interval can be expressed as a linear function of the p previous R-R intervals. Because our autoregression is on the R-R intervals, the time points at which these coefficients are computed are not evenly spaced. This type of history dependence makes our model an extension of the class of Wold point-process models (14). Moreover, our model is assumed to be nonstationary; therefore, the autoregession coefficients (Eq. 8) and the process variance (Eq. 10) are time varying.

Our analysis of the heartbeat interval series from the subjects in the tilt-table study showed that 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 (7, 16, 25, 26, 46, 50). Together, the instantaneous heart rate, instantaneous mean R-R interval, R-R standard deviation, and heart rate standard deviation provided a different signature for each of the three types of postural changes (Fig. 6; see Fig. S1 in supplemental data). The most interesting of these was the transient increase in heart rate standard deviation when a subject stood upright (Fig. 6; see Fig. S1 in supplemental data), which closely resembles the pattern and time course of change in cardiac output reported when a subject stood upright (46, 50). The standard indexes of heart rate variability, SDNN, and LF/HF computed from our heart rate estimates differ from those computed from reciprocal R-R interval series. This was most evident in our analysis of LF/HF. Although our new spectral estimates and the standard estimates appear similar and both show the expected increase in LF/HF with tilt, the distribution of their power was significantly different, in that our new estimates had power beyond the standard cutoff frequency.

Our methods suggest several new research directions and applications. 1) Studies of our methods in other protocols are needed to establish the physiological relevance of our new time and frequency domain measures of heart rate dynamics. 2) Our methods suggest new algorithms for simulating human heartbeat intervals using our HDIG model and standard point-process simulation methods (14) as well as a method for identifying and replacing ectopic beats. 3) Our methods may provide a means to characterize normal and pathological conditions in terms of heart rate variability. These analyses could be used to diagnose disease and follow its progression. 4) We are studying the application of our algorithms to real-time monitoring of fetal heart rate variability during antenatal and intralabor studies, as well as real-time monitoring of patients in the intensive care unit and in the operating room. 5) Studying cardiovascular control by combining point-process filtering (11, 15) algorithms with extensions of these methods to state-space (45) representations of the autonomic nervous system may help clarify the relation between stochastic and deterministic models of human heartbeat dynamics.

	GRANTS

TOP ABSTRACT THEORY AND METHODS RESULTS DISCUSSION GRANTS REFERENCES

 
This study was supported in part by National Institutes of Health Grants MH-59733, MH-61637, and DA-015644 and National Science Foundation Grant IBN-0081458 to E. N. Brown and by the PASTEUR Educational Program at Harvard Medical School and the Doris Duke Charitable Foundation Clinical Research Fellowship Program to E. C. Matten.

	ACKNOWLEDGMENTS

 
We are grateful to Roger G. Mark and Thomas Heldt (Harvard-MIT Division of Health Sciences and Technology) for kindly providing the tilt-table data analyzed in this study.

Part of this research was performed while E. N. Brown was on sabbatical at the Laboratory for Information and Decision Systems in the Department of Electrical Engineering and Computer Science at MIT.

	FOOTNOTES

 

Address for reprint requests and other correspondence: R. Barbieri, Neuroscience Statistics Research Laboratory, Dept. of Anesthesia and Critical Care, Massachusetts General Hospital, 55 Fruit St., Clinics 3, Boston, MA 02114-2696 (E-mail: barbieri{at}neurostat.mgh.harvard.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 THEORY AND METHODS RESULTS DISCUSSION GRANTS REFERENCES

 

1. Akselrod S, Gordon D, Ubel FA, Shannon DC, Berger AC, and Cohen RJ. Power spectrum analysis of heart rate fluctuation: a quantitative probe of beat-to-beat cardiovascular control. Science 213: 220–222, 1981.
2. Barbieri R, Quirk MC, Frank LM, Wilson MA, and Brown EN. Construction and analysis of non-Poisson stimulus-response models of neural spiking activity. J Neurosci Methods 105: 25–37, 2001.
3. Barbieri R, Triedman JK, and Saul JP. Heart rate control and mechanical cardiopulmonary coupling to assess central volume: a systems analysis. Am J Physiol Regul Integr Comp Physiol 283: R1210–R1220, 2002.
4. Barbieri R, Waldmann RA, Di Virgilio V, Triedman JK, Bianchi AM, Cerutti S, and Saul JP. Continuous quantification of baroreflex and respiratory control of heart rate by use of bivariate autoregressive techniques. Ann Noninvasive Electrocardiol 1: 264–277, 1996.
5. Bianchi AM, Mainardi L, Petrucci E, Signorini MG, Mainardi M, and Cerutti S. Time-variant power spectrum analysis for the detection of transient episodes in HRV signal. IEEE Trans Biomed Eng 40: 136–144, 1993.
6. Bigger JT, Fleiss JL, Steinman RC, Rolnitzky LM, Kleiger RE, and Rottman JN. Frequency domain measures of heart period variability and mortality after myocardial infarction. Circulation 85: 164–171, 1992.
7. Borst C, Wieling W, van Brederode JF, Hond A, de Rijk LG, and Dunning AJ. Mechanisms of initial heart rate response to postural change. Am J Physiol Heart Circ Physiol 243: H676–H681, 1982.
8. Box GEP, Jenkins GM, and Reinsel GC. Time Series Analysis, Forecasting and Control (3rd ed.). Englewood Cliffs, NJ: Prentice-Hall, 1994.
9. Brown EN, Barbieri R, Eden UT, and Frank LM. Likelihood methods for neural data analysis. In: Computational Neuroscience: A Comprehensive Approach, edited by Feng J. London: Chapman & Hall/CRC, 2003.
10. Brown EN, Barbieri R, Ventura V, Kass RE, and Frank LM. The time-rescaling theorem and its application to neural spike train data analysis. Neural Comput 14: 325–346, 2002.
11. Brown EN, Nguyen DP, Frank LM, Wilson MA, and Solo V. An analysis of neural receptive field plasticity by point process adaptive filtering. Proc Natl Acad Sci USA 98: 12261–12266, 2001.
12. Cerutti S, Baselli G, Bianchi AM, Mainardi LT, Signorini MG, and Malliani A. Cardiovascular variability signals: from signal processing to modeling complex physiological interactions. Automedica 16: 45–69, 1994.
13. Chikara RS and Folks JL. The Inverse Gaussian Distribution: Theory, Methodology, and Applications. New York: Dekker, 1989.
14. Daley D and Vere-Jones D. An Introduction to the Theory of Point Process (2nd ed.). New York: Springer-Verlag, 2003.
15. Eden UT, Frank LM, Barbieri R, Solo V, and Brown EN. Dynamic analysis of neural encoding by point process adaptive filtering. Neural Comput 16: 971–998, 2004.
16. Ewing DJ, Hume L, Campbell IW, Murray A, Neilson JM, and Clarke BF. Autonomic mechanisms in the initial heart rate response to standing. J Appl Physiol 49: 809–814, 1980.
17. Ewing DJ, Neilson JM, and Travis P. New method for assessing cardiac parasympathetic activity using 24-hour electrocardiograms. Br Heart J 52: 396–402, 1984.
18. Freeman R, Lirofonis V, Farquhar WB, and Risk M. Limb venous compliance in patients with idiopathic orthostatic intolerance and postural tachycardia. J Appl Physiol 93: 636–644, 2002.
19. Freeman R, Saul JP, Roberts MS, Berger RD, Broadbridge C, and Cohen RJ. Spectral analysis of heart rate in diabetic autonomic neuropathy. Arch Neurol 48: 185–190, 1991.
20. Freeman RK. Problems with intrapartum fetal heart rate monitoring interpretation and patient management. Obstet Gynecol 100: 813–826, 2002.
21. Fuller WA. Introduction to Statistical Time Series. New York: Wiley, 1976.
22. Guyton AC. Textbook of Medical Physiology (8th ed.). Philadelphia, PA: Harcourt Brace, 1991.
23. Halliwill JR, Lawler LA, Eickhoff TJ, Joyner MJ, and Mulvagh SL. Reflex responses to regional venous pooling during lower body negative pressure in humans. J Appl Physiol 84: 454–458, 1998.
24. Hayano J, Taylor JA, Yamada A, Mukai S, Hori R, Asakawa T, Yokoyama K, Watanabe Y, Takata K, and Fujinami T. Continuous assessment of hemodynamic control by complex demodulation of cardiovascular variability. Am J Physiol Heart Circ Physiol 264: H1229–H1238, 1993.
25. Heldt T, Oefinger MB, Hoshiyama M, and Mark RG. Cir