## Abstract

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 a vital moment-to-moment indicator of cardiovascular integrity measured on every physical examination. Heart rate is also monitored continuously in patients under anesthesia during surgery and in those treated in an intensive care unit, as well as in fetuses during labor. Heart rate variability is an important quantitative marker of cardiovascular regulation by the autonomic nervous system. Its significance was first appreciated over 40 years ago, when it was discovered that fetal distress is associated with appreciable changes in heart rate variability before any change in heart rate (27). Physicians routinely use Holter monitor studies, i.e., 24–72 h of continuous electrocardiography, in which heart rate variability is assessed to diagnose diseases that affect the autonomic nervous system, follow their progression, and measure the efficacy of therapy. Such diseases include diabetes, Guillain-Barre syndrome, multiple sclerosis, Parkinson's disease, and Shy-Drager orthostatic hypotension (17, 19, 32, 38). Loss of heart rate variability is an independent predictor of mortality after an acute myocardial infarction (6, 33, 35), is indicative of ventricular dysfunction in patients with congestive heart failure (41, 42), and is associated with fetal distress during labor (20).

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

#### 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*_{1} < *u*_{2} <, …, < *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 *k*th R-R interval, μ(*H*_{uk}θ) = θ_{0} + ∑_{j=1}^{p} θ_{j}*w*_{k−j+1} > 0 is the mean, θ_{p + 1} > 0 is the scale parameter, and θ = (θ_{0},θ_{1},…,θ_{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},θ) = θ_{0}, *f*(*t*|*H*_{uk},θ) = *f*(*t*|*u*_{k},θ_{0},θ_{1}), 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 × 10^{4} 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*_{1} < *u*_{2} <,…,< *u*_{nt} ≤ *t*. We let *u*_{t − l:t} = (*u*_{1},…,*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}|*H**unt*,θ_{t}), and if *t* > *u*_{nt} and the last interval is censored, then the contribution to the local likelihood is log ∫_{t−unt}^{∞}*f*(*v*|*H*_{unt},θ_{t})d*v*. 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.

The study was conducted at the Massachusetts Institute of Technology (MIT) General Clinical Research Center (GCRC) and was approved by the MIT Institutional Review Board and the GCRC Scientific Advisory Committee. The subjects were five men and five women: age (mean ± SD) 28.7 ± 1.2 yr, height 172.8 ± 4.0 cm and weight 70.6 ± 4.5 kg. The protocol began with each subject lying supine for 5 min; then each subject underwent three types of up-down pairs. The tilt pairs were the following: rapid up tilt, in which the subject was moved from horizontal to vertical in <3 s with a rapid down tilt, in which the subject was moved from vertical to horizontal in <3 s; slow up tilt, in which the subject was moved from horizontal to vertical in ∼1 min with a slow down tilt, in which the subject was moved from vertical to horizontal in ∼1 min; and standing upright, in which the subject stood up immediately supporting his/her own weight with lying supine, in which the subject lay supine immediately from standing supporting his/her own weight. The subject remained in each tilt state for 3 min, and each tilt pair (rapid up-down, slow up-down, and stand up-lie supine) was repeated twice. The order of the tilt pairs was chosen at random for each subject. The protocol lasted 55–75 min (3,300–4,500 s). A single-lead ECG was continuously recorded for each subject during the study, and the R-R intervals were extracted using a curve length-based QRS detection algorithm (58).

## RESULTS

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

Figure 2 illustrates the range of KS plots obtained from the fits of the three models to the heartbeat interval series of the 10 subjects. These include an example of the best HDIG fit (Fig. 2*A*), the best improvement of an HDIG model fit relative to an RIG model fit (Fig. 2*B*), the best RIG fit (Fig. 2*C*), and the worst HDIG fit (Fig. 2*D*). The KS distances (Table 1) can be seen in these plots by computing for each model the maximum vertical distance between the KS plot curve for each model and the 45° line (Fig. 2, thin black diagonal line). As suggested by the KS distances (Table 1), the LA model gave the poorest fits to the heartbeat interval series, inasmuch as their KS plots (Fig. 2*D*) showed near-maximal deviation from the 45° lines, which correspond to an exact fit. With the exception of *subject 7* (Fig. 2*C*), for whom the agreement with the data was excellent, the KS plots for the RIG model (Fig. 2) lie above the 95% confidence bounds for lower quantiles and below the 95% confidence bounds for upper quantiles.

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. 2*B*), 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. 2*D*). 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. 3*A*). The plot of τ_{k} vs. τ_{k − 1} for the LA model (Fig. 3*B*) 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. 3*C*) 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. 3*D*), with a correlation coefficient of 0.05.

To further assess the independence of the transformed intervals from the HDIG model fits, we computed for each subject the correlation function of the τ_{k} values for 60 lags to evaluate the correlation structure beyond *lag 1* (Fig. 4). The largest values of the correlation function were the *lag 4* correlation of 0.15 for *subject 8*, the *lag 3* correlation of 0.09 for *subject 4*, the *lag 4* correlation of 0.09 for *subject 6*, and the *lag 2* correlation of 0.09 for *subject 1*. The remaining correlations for all 10 subjects were within the approximate 95% confidence bounds for the null hypothesis of no correlation (8). These results show that any nonzero correlations in the τ_{k} values were small, consistent with the idea that, for each subject, these transformed R-R intervals for the HDIG model were approximately independent for all subjects. These results, together with the KS distance and KS plot findings, suggest that of the three models, the HDIG model best describes original human heartbeat interval data for all the subjects.

#### Dynamic analysis of human heartbeat intervals.

For each subject, we used the HDIG model to analyze heart rate variability by computing instantaneous estimates of the mean R-R interval (*Eq. 9*), heart rate (*Eq. 11*), R-R interval standard deviation (*Eq. 10*), and heart rate standard deviation (*Eq. 12*) for the entire study. The results of the entire heart rate and heart rate variability analysis for *subject 2* are shown in Fig. 5.

Our heart rate and heart rate variability indexes followed almost instantaneously the changes in the heartbeat interval series in Fig. 1. The estimated instantaneous heart rate for the HDIG model (Fig. 5*B*) was more similar to the reciprocal of the R-R intervals (Fig. 5*A*) across all rest and upright posture periods. The instantaneous heart rate estimate obtained from the RIG or LA model (Fig. 5*C*) were much smoother than the HDIG estimates, even though all three models use the same local likelihood interval of 60 s. We showed only the RIG instantaneous heart rate estimate, inasmuch as the LA and RIG model estimates are algebraically equivalent. This demonstrates that computation of heart rate as the reciprocal of the average of the R-R intervals in a given time period (LA model) is equivalent to estimation in which heartbeat intervals are represented as an RIG model.

The instantaneous R-R interval standard deviation (Fig. 5*D*) and heart rate standard deviation (Fig. 5*E*) 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. 6*A*), a slow tilt (Fig. 6*B*), and an upright position (Fig. 6*C*), 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. 6*A*). 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.

Throughout the slow upward tilt, the heart rate of *subject 4* progressively increased, whereas the mean R-R interval decreased (Fig. 6*B*). Paralleling the increase in heart rate and the decrease in mean R-R interval, the R-R interval standard deviation also decreased. As the tilt was executed, the heart rate standard deviation did not decrease but showed more fluctuations. During the upright posture, the subject’s heart rate and mean R-R interval fluctuated with short paroxysms of increases and decreases. With the gradual decrease in the tilt angle, the subject’s heart rate decreased and his mean R-R interval increased toward their pretilt levels. The increase in the R-R interval standard deviation began in advance of the downward tilt and continued upward for the balance of the rest period, whereas the heart rate standard deviation initially began its rise in advance of the downward tilt. However, with the decrease in the tilt angle, the heart rate standard deviation initially decreased.

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. 6*C*). In a manner analogous to that observed for the rapid tilt (Fig. 6*A*), 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.

In the frequency domain analysis, we computed LF/HF as the ratio of the spectral power of the heart rate series at 0.04–0.15 Hz (LF) to that at 0.15–0.4 Hz (HF). LF/HF for the standard analysis was derived by spectral analysis of heart rate time series computed by cubic spline interpolation and resampling at 3 Hz (3). LF/HF for our heart rate series was derived by spectral analysis of our series of mean heart rate. Because our indexes were computed at 200 Hz, i.e., we updated our local likelihood estimates every 5 ms, we downsampled each one at 3 Hz for a more accurate, direct comparison.

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. 7*E*). 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. 7*B*). 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 × (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.

These results suggest that summaries comparable to SDNN and LF/HF analyses can be performed with heart rate series derived from our HDIG model. These results, coupled with the goodness-of-fit analyses, which demonstrated that our models offer an accurate description of the stochastic structure in the heartbeat interval series, suggest that the time domain summaries and LF/HF derived from our methods may be a more accurate description of these quantities.

## DISCUSSION

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

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

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.

- Copyright © 2005 by the American Physiological Society