Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
Published in final edited form as:
PMCID: PMC2952361

Characterizing Nonlinear Heartbeat Dynamics within a Point Process Framework

Zhe Chen, Member, IEEEcorresponding author
Neuroscience Statistics Research Laboratory, Harvard Medical School, Massachusetts General Hospital, Boston, MA 02114 USA, and also with the Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139 USA (ude.tim.tatsoruen@nehcehz)
Emery N. Brown, Fellow, IEEE
Neuroscience Statistics Research Laboratory, Harvard Medical School, Massachusetts General Hospital, Boston, MA 02114 USA, and also with the Harvard-Massachusetts Institute of Technology (MIT) Division of Health Science and Technology, and the Department of Brain and Cognitive Sciences, MIT, Cambridge, MA 02139 USA (ude.tim.tatsoruen@bne)
Riccardo Barbieri, Senior Member, IEEE


Human heartbeat intervals are known to have nonlinear and nonstationary dynamics. In this paper, we propose a model of R–R interval dynamics based on a nonlinear Volterra–Wiener expansion within a point process framework. Inclusion of second-order nonlinearities into the heartbeat model allows us to estimate instantaneous heart rate (HR) and heart rate variability (HRV) indexes, as well as the dynamic bispectrum characterizing higher order statistics of the nonstationary non-Gaussian time series. The proposed point process probability heartbeat interval model was tested with synthetic simulations and two experimental heartbeat interval datasets. Results show that our model is useful in characterizing and tracking the inherent nonlinearity of heartbeat dynamics. As a feature, the fine temporal resolution allows us to compute instantaneous nonlinearity indexes, thus sidestepping the uneven spacing problem. In comparison to other nonlinear modeling approaches, the point process probability model is useful in revealing nonlinear heartbeat dynamics at a fine timescale and with only short duration recordings.

Keywords: Adaptive filters, approximate entropy (ApEn), heart rate variability (HRV), nonlinearity test, point processes, scaling exponent, Volterra series expansion

I. Introduction

The human heartbeat is regulated by the autonomic nervous system, and as a result, heart rate (HR) and heart rate variability (HRV) measurements extracted from the ECG are important quantitative markers of cardiovascular control [1]. A healthy heart is influenced by multiple neural and hormonal inputs that result in variations of the interbeat interval duration. Specifically, various nonlinear neural interactions and integrations occur at the neuron and receptor levels, and underlie the complex output of structures such as the sinoatrial (SA) node in response to changing levels of sympathetic and vagal activities [55]. The complex nature of heartbeat dynamics has been widely considered and discussed in cardiovascular literature. Although detailed physiology behind these complex dynamics has not been completely clarified, several nonlinearity measures of HRV have been pointed out as important quantifiers of complexity of cardiovascular control and have been proved to be of important prognostic value in aging and diseases [4], [25], [26],[46], [58], [60], [62].

Many physiological signals are known to be nonlinear and nonstationary. In biomedical engineering, various nonlinear indexes, such as the Lyapunov exponent, the fractal exponent, or the approximate entropy (ApEn), have been proposed to characterize the nonlinear behavior of the underlying physiological system (e.g., [2]). It has been suggested that such nonlinearity indexes might provide informative indicators for diagnosing cardiovascular or brain diseases. Notably, some difficulties have been often encountered when validating these indexes, such as the presence of noise or artifact, the limited size of data samples, or the low sampling rate of the observed signals. All these issues shall be kept in mind when new statistical indexes are estimated from real signals recorded from a nonlinear system.

In characterizing the nonlinear heartbeat dynamics, both linear and nonlinear system identification methods have been applied to R–R interval series [19], [20], [61]. Examples of higher order characterization for cardiovascular signals, include nonlinear autoregressive (AR) models, Volterra–Wiener series expansion, and Volterra–Laguerre models [2], [32], [33], [36]. Several authors have demonstrated the feasibility and validity of nonlinear AR models, suggesting that future HR dynamics studies should put greater emphasis on nonlinear analysis [19], [20], [31], [61]. However, none of these models have included nonlinear elements in a framework based on a precise statistical characterization of the heartbeat generation process, and all of mentioned studies used either beat series (tachograms) or discretionarily interpolated R–R time series instead of deriving model estimates. In this paper, we apply nonlinear modeling to heartbeat dynamics using a point process paradigm. The point process theory is a powerful statistical tool able to characterize the probabilistic generative mechanism of the heartbeat at each moment in time, thus allowing for estimation of instantaneous HR and HRV measures [7], [8]. Furthermore, inclusion of second-order nonlinear terms to the point process model offers an opportunity to monitor dynamic higher order spectra indexes [39], [40].

The paper is organized as follows. Section II presents some background on nonlinear system identification by Volterra–Wiener series expansion. Section III gives a brief exposition of probabilistic point process model theory for heartbeat intervals, derives the instantaneous HR and HRV indexes, and reviews the adaptive point process filtering algorithm as well as the goodness-of-fit tests. Section IV is devoted to the instantaneous higher order spectral analysis and derivation of the dynamic bispectrum estimate, as well as the nonlinearity test for R–R interval series. Section V describes the synthetic data generated to test the models, as well as two experimental heartbeat datasets. Section VI presents the experimental results on all datasets using the point process models, discussing model selection, nonlinearity assessment, performance comparison, and irregularity characterization. Finally, discussions and conclusion are given in Section VII.

II. Volterra Series for Nonlinear System Identification

The Volterra series expansion, based on the Volterra theorem, is a general method for nonlinear system modeling and identification [36]. In functional analysis, a Volterra series denotes a functional expansion of a dynamic, nonlinear, and time-invariant function. The Volterra series allows for representation of a wide range of nonlinear systems. Because of its generality, Volterra series expansion has been widely used in nonlinear modeling in engineering and physiology [2], [31], [36]. For instance, computational procedures based on a comparison of the prediction power of linear and nonlinear models of the Volterra–Wiener form have been applied to measure the complex dynamics of the heartbeats [6]. However, it shall be pointed out that all of these nonlinear models used only raw R–R intervals without modeling the point process nature of the heartbeats.

Consider a nonlinear single-input and single-output system y = g(x). According to the Volterra series theory, the nonlinear system can be expanded by a (finite or infinite) set of kernel expansion terms


where M is the memory of the nonlinear system. Equation (1) only includes up to the second-order nonlinear term in the Volterra series expansion; however, inclusion of higher order terms is possible. The Volterra kernels {k0,k1,k2,…,} describe the dynamics of the system, each of which is associated with Volterra coefficients at different kernel orders and different time lags. Estimation of the Volterra coefficients is generally performed by computing the coefficients of an orthogonalized series, and then recomputing the coefficients of the original Volterra series. A common method is based on the least squares optimization [36]. In this paper, we apply a point process adaptive filtering approach to recursively estimate the time-varying Volterra coefficients.

III. Heartbeat Interval Point Process Model

A random point process is a random element whose values are “point patterns” on a set, where a point pattern is specified as a locally finite counting measure [23]. Specifically in the time domain, a simple 1-D point process consists of series of binary (0 and 1) observations, where the variables 1 marks the occurrence times t [sm epsilon] [0,∞) of the random events. Mathematically, we let N(t) define a continuous-time counting process, and let its differential dN(t) denote a continuous-time indicator function, where dN(t) = 1, when there is an event (such as the ventricular contraction) or dN(t) = 0, otherwise. Point process theory has been widely used in modeling various types of random events (e.g., eruptions of earthquakes, queueing of customers, spiking of neurons, etc.) where the timing of the events are of central interest. Bearing a similar spirit, the point process theory has been used for modeling human heartbeats [7],[8],[16]. The point process framework primarily defines the probability of having a heartbeat event at each moment in time. A parametric formulation of the probability function allows for a systematic, parsimonious estimation of the parameter vector in a recursive way and at any desired time resolution. Instantaneous indexes can then be derived from the parameters in order to quantify important features as related to cardiovascular control dynamics.

A. Heartbeat Interval

Suppose we are given a set of R-wave events {uj}j=1J detected from the ECG, let RRj = ujuj−1 > 0 denote the jth R-R interval, or equivalently, the waiting time until the next R-wave event. By treating the R-wave as discrete events, we may develop a point process probability model in the continuous time domain [7].

Assuming history dependence, the probability distribution of the waiting time tuj until the next R-wave event follows an inverse Gaussian model:


where uj denotes the previous R-wave event occurred before time t, μRR(t) represents the first-moment statistic (mean) of the distribution, and θ > 0 denotes the shape parameter of the inverse Gaussian distribution, whose role is to model the tail shape of the distribution (when θ → ∞, the inverse Gaussian distribution becomes more like a Gaussian distribution). As p(t) indicates the probability of having a beat at time t given that a previous beat has occurred at uj and μRR(t) can be interpreted as signifying the most probable moment when the next beat could occur. By definition, p(t) is characterized at each moment in time, at the beat as well as in-between beats. We can also estimate the second-moment statistic (variance) of the inverse Gaussian distribution as σRR2(t)=μRR3(t)θ. The use of an inverse Gaussian distribution to characterize the R–R intervals' occurrences is motivated by the fact that if the rise of the membrane potential to a threshold initiating the cardiac contraction is modeled as a Gaussian random walk with drift, then the probability density of the times between threshold crossings (the R–R intervals) is indeed the inverse Gaussian distribution [7]. In [16], we have compared heartbeat interval fitting point process models using different probability distributions, and found that the inverse Gaussian model achieved the overall best fitting results. The parameter μRR(t) denotes the instantaneous R–R mean that can be modeled as a generic function of the past (finite) R–R values μRR(t)=g(RRt−1,RRt−2,…,RRth), where RRtj denotes the previous jth R–R interval occurred prior to the present time t. In our previous work [8], [14], [16], the history dependence is defined by expressing the instantaneous mean μRR(t) as a linear combination of present and past R–R intervals (in terms of an AR model), i.e., function g is linear. Here, we propose to include the nonlinear terms of past R–R intervals by defining the instantaneous RR mean as follows:


where <RR>t=1hk=1hRRtk. Here the coefficients a0(t),{ai(t)}, and {bkl(t)} correspond to the time-varying zero-, first-, and second-order Volterra kernel coefficients. The zero order coefficient a0 accounts for the nonzero mean of the R–R series. Equation (2) can be interpreted as a discrete Volterra–Wiener series with degree of nonlinearity d = 2 and memory h = max{p, q} [6]. As μRR(t) is defined in a continuous time fashion, we can obtain an instantaneous R–R mean estimate at a very fine timescale (with an arbitrarily small bin size Δ), which requires no interpolation between the arrival times of two beats. Given the proposed parametric model, the nonlinear indexes of the HR and HRV will be defined as a time-varying function of the parameters ξ(t)=[a0(t),a1(t),…,ap(t),b11(t),…,bqq(t),θ(t)].

B. Instantaneous Indices of HR and HRV

HR is defined as the reciprocal of the R–R interval. For t measured in seconds, a new variable r=c(tuj)−1 (where c = 60 s/min) can be defined in beats per minute (bpm). By the change-of-variables formula, the HR probability p(r)=p(c(tut)−1) is given by


and the mean and the standard deviation of HR r can be derived [7], [8], as given by μHR and standard deviation σHR, respectively



where μ~=c1μRR and θ~=c1θ.

It is known from point process theory [7], [8], [13] that the conditional intensity function (CIF) λ(t) is related to the interevent probability p(t) with a one-to-one relationship


The estimated CIF can be used to evaluate the goodness-of-fit of the proposed probability model for the heartbeat interval point process probability model. The quantity λ(t)Δ yields approximately the probability of observing a beat during the [t,t + Δ) interval in the sense that [23]


where Ht denotes all of available history information (subject to causality) up to time t.

C. Adaptive Point Process Filtering

In order to track the unknown parameters of vector ξ in a nonstationary environment, we can recursively estimate them via adaptive point process filtering [8]. Upon time discretization, we have the following equation updates at discrete-time index k:





where P and W denote the parameter and noise covariance matrices, respectively; Δ =0.005s denotes the time bin size; Δλk = [partial differential]λk/[partial differential]ξk and Δ2λk=2λkξkξkT denote the first- and second-order partial derivatives of the CIF with respect to ξ at time t = kΔ, respectively. The indicator variable nk = 1 if a heart beat occurs in time ((k−1)Δ,kΔ] and 0 otherwise. The point process filtering described in (7) through (10) can be viewed as a “point process analog” of the Kalman filtering (for continuous-valued observations). In (7) and (8), the a priori estimates ξk|k−1 and Pk|k−1 are computed, while in (9) and (10), the a posteriori estimates ξk|k and Pk|k are computed. In (9), [nk − λkΔ] can be viewed as the innovations term computed from the point process filter, and (10) is derived based on a Gaussian approximation of the log posterior. Clearly, as the innovations term is likely to be nonzero even in the absence of a beat, the parameters are always updated at each step.

Once the vector ξ has been estimated within (0, T], one can compute the probability density function (pdf) p(t) as well as the CIF estimate λ(t) [from (6)] in time interval (0, T]. Furthermore, we can compute the cumulative log-likelihood (denoted by L) of the point process observations [13], [23]


where the indicator variable dN(τ) = 1 (or nk = 1), if a beat occurs at time τ (or within the time interval ((k − 1)Δ, kΔ]) and dN(τ) = 0 (or nk = 0), otherwise. The log-likelihood function (11) defines the logarithm of the joint probability of all random events (i.e., beats), and the second equality holds when the bin size Δ is sufficiently small, and [ell]k = nk logλk−λk approximates the instantaneous log-likelihood function [ell](t) = logλ(t)dN(t)/dt−λ(t).

D. Model Selection and Goodness-of-fit Tests

Our method requires the user to predetermine a proper model order {p, q} for the Volterra series expansion. In general, a tradeoff between model complexity and goodness-of-fit arises when a point process model is considered. In practice, the order of the model (2) may be determined based on the akaike information criterion (AIC) (by prefitting a subset of the data using either point process filter or local likelihood method [7], [35]) as well as the KolmogorovSmirnov (KS) statistic in the post hoc analysis. For different values p and q, we can compare the AIC and choose the parameter setup with the minimum AIC value


where dim(ξ) = p + q2 + 2 denotes the dimensionality of parameter vector ξ in the nonlinear model. Once the order {p, q} is determined, the initial Volterra coefficients will be estimated by the method of least squares [59]: specifically, the coefficients {ai} are optimized by solving a Yule–Walker equation for the linear part using the first 200 sample points, and the coefficients {bij} are estimated by fitting the residual error via least squares. Here, we use a separate estimation instead of a joint estimation procedure for the Volterra coefficients because we like to preserve the interpretation of the linear AR coefficients (such as the stability, which is assured by keeping the roots inside the unit circle). A joint estimation procedure is possible based on orthogonal projection, cross correlation, or least squares [36],[59], but it may destroy the structure described by the linear AR coefficients {ai}, which will be used to estimate the parametric AR spectrum defined later.

The goodness-of-fit of the point process model is based on the KS test [13]. Given a point process specified by J discrete events: 0 <u1 <…<uJ <T, the random variables zj=uj1ujλ(τ)dτ are defined for j = 1, 2, …, J − 1. If the model is correct, then the variables vj = 1−exp(−zj) are independent, uniformly distributed within the region [0, 1], and the variables gj = Φ−1(vj) (where Φ(•) denotes the cumulative distribution function (cdf) of the standard Gaussian distribution) are sampled from an independent standard Gaussian distribution. To compute the KS test, the vj s are sorted from smallest to largest, and plotted against the cdf of the uniform density defined as j − 0.5/J. If the model is correct, the points should lie on the 45° line. The 95% confidence interval lines are defined as y=x±1.36/(J−1)1/2. The KS distance, defined as the maximum distance between the KS plot and the 45° line, is used to measure the lack-of-fit between the model and the data. The autocorrelation function of the gj s: ACF(m)=1Jmj=1Jmgjgj+m, can also be computed. If the gjs are independent, ACF(m) shall be small (around 0 and within the 95% confidence interval 1.96/(J−1)1/2) for any lag m.

V. Quantitative Tools: Spectral Analysis and Nonlinearity Test

A. Instantaneous Higher Order Spectral Analysis

Given the Volterra–Wiener expansion for the instantaneous R–R interval mean {μRR(t)}, we may compute the time-varying parametric (linear) autospectrum


By integrating (12) in each frequency band, we may compute the index within the very low frequency (VLF) (0.01–0.05 Hz), LF (0.05–0.15 Hz), or HF (0.15–0.5 Hz) ranges. In addition, let B(f1,f2,t)=k=1ql=1qbkl(t)ej2kπf1ej2kπf2 denote the Fourier transform of the second-order kernel coefficients {bkl(t)} (all of which together are viewed as discrete samples from a 2-D impulse response function). From (2), it is known that [39], [40]


where C(f1,f2t) denotes the bispectrum (Fourier transform of the third-order moment). Note that we use the approximation “≈” instead of equality “=” in (13), since the equality only strictly holds when the input variables are jointly Gaussian, which is not necessarily true in our case. The bispectrum is an important tool for evaluating the presence of nonlinearity in stationary time series [9], [39], [40]. From (13), we then can estimate the dynamic bispectrum C(f1,f2,t). From the Parseval theorem, we also know that that the sum (or integral) of the square of a function is equal to the sum (or integral) of the square of its transform, namelybkl2=B(f1,f2)2=B(f1,f2)2=B(f1,f2)2 (the second equality follows from the conjugate symmetry property), or bkl=B(f1,f2). Let b(t) denote a vector that contains all of coefficients {bkl(t)}, in light of (13), we may compute an index that quantifies the fractional contribution of the linear terms on the total power as follows:


where | · | denotes either the norm of a vector or the modulus of a complex variable. The spectrum norm defines the area integrated over the frequency range under the spectral density curve. Since the norm units of spectral and bispectral density are the same, their ratio ρ(t) is dimensionless (note that the unit of {bkl}{bkl} is in 1/second, and the unit of norm |Q(f, t)| is in second, thus their product is unitless). As a function of the estimated parameters, this index can also be estimated at each moment in time and is updated at the beat as well as in-between beats. The ratio ρ defined in (14) can be viewed as a dynamic counterpart of the following static power ratio [footnote: As one reviewer pointed out, the ratio defined in (15) hardly reaches close to 0, and the reviewer also suggested to define an alternative ratio index (power spectrum – bispectrum)/power spectrum, which is bounded between 0 and 1]:


where the power spectrum and bispectrum will be calculated by Fourier transform using the observed (nonequally spaced) R–R interval time series (with a stationarity assumption). A lower value of the ratio (15) implies that the fraction of the bispectral power is higher, thus pointing at a more significant nonlinear component in the time series.

It shall be noted that the frequency units appearing in Q(f,t) and B(f1,f2,t) are both cycles/beat, since the autoregression of μRR (t) is conducted on previous beats instead of previous instantaneous {μRR (ti)}. This alternative modeling strategy, however, would require a large number of tags in the linear AR model to compensate for the use of fine timescale. Another way to compute the spectra of interest in the unit of cycles/second is to consider the estimated {μRR(t)}or {μHR(t)} series and compute the power spectrum or bispectrum using a direct method. However, this would require a windowing technique and would not allow for instantaneous estimates. As a consequence of the change from cycles/second to cycles/beat, certain spectral distortion between the spectrum of counts (SOC) and the spectrum of interval (SOI) might be expected [11], [24], [37], [48], especially when the beat-to-beat intervals have a large variance. As this issue could become critical when precise estimate of specific oscillatory frequencies are needed, its effects are less noticeable in total power computations.

B. Nonlinearity Test

In the literature, there are many nonlinearity indexes being proposed for time series analysis, such as the correlation dimension [29], the Lyapunov exponent [2], [29], [51], the time reversibility index [29], [50], and the prediction error [6]. Common methods require the computation of surrogate time series in order to construct a hypothesis test. The standard procedure is to assume Gaussianity and stationarity, and to perform a Fourier transform followed by phase randomization and inverse Fourier transform (such a procedure preserves the first and second-order moment statistics while discarding the phase information). In this paper, we consider a specific established time-domain method [9] as applied to the R–R time series for testing the presence of nonlinearity in the heartbeat intervals.

The test developed in [9] uses a phase scrambled bootstrap technique for testing the presence of nonlinearity of a time series based upon the third-order moment statistics. The basic idea of this method is to compare the estimated third-order moment of the tested series with a set of limits generated from linear stationary phase scrambled bootstrap data: large differences shall indicate nonlinearity or possibly nonstationarity [9]. The null hypothesis assumes that the given time series is linear and stationary. The result of the hypothesis test is either H = 0 (which indicates that the null hypothesis is accepted and P >0.05) or H = 1 (which indicates that the null hypothesis is rejected with 95% confidence). In the considered simulated series and real recordings, we restricted the test to short-term dependence by setting the number of laps M = 8, and a total of 500 bootstrap replications were simulated for every test.

V. Data

In order to test the tracking ability of the nonlinear point process model, and to compare its performance with the standard filter with only linear dynamics, we generate a synthetic heartbeat dataset. Specifically, without postulating a second-order nonlinear system as assumed in our model, we used the chaotic Rössler time series governed by the following differential equations [49] [Footnote: Of note, the heartbeat dynamics reflected in the synthetic set are not directly associated with real physiological generation mechanisms, and it is neither implied that the heartbeat dynamics be chaotic].


The time series were simulated by the Runge–Kutta integration using conditions a = 0.15, b = 0.20, and c = 10.0, with step size of 0.01. A total of 3000 data points were generated, one for every three x-axis values was chosen, and 1000 data points (representing the generated R–R intervals) were finally selected. The simulated deterministic time series is illustrated in Fig. 1 (first panel). The nonlinearity test described earlier indicates the synthetic time series are significantly nonlinear (H = 1 and P < 1e-6).

Figure 1
Synthetic R–R interval series and its estimated indexes μRR (superimposed on R–R series), σRR, μHR, and σHR, using a nonlinear model (10,4).

In order to validate the proposed algorithms' performance as related to real physiological dynamics, we have considered two experimental datasets. The first heartbeat dataset was recorded under the “autonomic blockade assessment of the sympathovagal balance and respiratory sinus arrhythmia” protocol. Detailed description of the experimental data was given in [53]. The recorded R–R interval time series last about 5 min for each epoch. In the drug administered state, after a control recording stage in rest condition, either atropine (ATR, 0.04 mg/kg iv over 5 min, parasympathetic blockade) or propranolol (PROP, 0.2 mg/kg iv over 5 min, sympathetic blockade) was delivered to the subject. In the double blockade (DB) epoch, the inputs from both sympathetic and parasympathetic branches of the autonomic nervous system were suppressed [53]. A total of 17 healthy volunteers participated in the study. Due to space limit, the results of four representative subjects (two from the ATR group and two from the PROP group, both were randomly selected) are listed in Table I. These four subjects have been tested and reported on a previous analysis with a linear predictive model [14], [16]. In Fig. 2, we show the R–R interval series of one representative subject in the control supine condition.

Figure 2
Instantaneous heartbeat statistics computed from a representative subject (subject 11, control, supine, from the autonomic blockade protocol [53]) using a nonlinear model. In the first panel, the estimated μRR (black trace) is superimposed on ...
Results of Model Selection for Synthetic Data.

The second heartbeat dataset, which was retrieved from a public source: Physionet ( [28], consists of R–R time series recorded from 12 congestive heart failure (CHF) patients (from BIDMC-CHF Database) and 16 healthy subjects (from MIT-BIH Normal Sinus Rhythm Database). Each R–R time series was artifact-free (upon human's visual inspection and artifact rejection) and lasted about 50 min (small segments of the original over 20 h recordings). In Fig. 3, we show the R–R interval series from one representative healthy subject. Since these recordings have longer durations, they have been deemed as particularly suitable for studying complex heartbeat interval dynamics [42], [46].

Figure 3
Instantaneous indexes computed from a representative healthy subject (recording no. 16483, MIT-BIH Normal Sinus Rhythm Database from Physionet [28]) using a nonlinear model. In the first panel, the estimated μRR (black trace) is superimposed on ...

V1. Results

A. Model Selection and Goodness-of-Fit Tests

Using the synthetic dataset, we have conducted several analyses to assess model order selection for both linear and nonlinear models. The results of AIC and KS statistics are shown in Table I, which are computed from fitting all simulated data points. As seen from Table I, according to AIC, the best fit is given by the nonlinear model (8,4), followed by (10,4), whereas according to the KS statistic (smaller KS distance), the best fit is given by the nonlinear model (10,4). Overall, it is important to notice that for the same level of model complexity, the nonlinear model generally achieves a better KS statistic than the linear model, but only when the predictive power from the linear part is sufficient—this can be seen from the relatively poorer performances of nonlinear models (4,4) and (5,9) in Table I. In this analysis, we selected the nonlinear model (10,4) as the optimal nonlinear model (with estimated instantaneous indexes shown in Fig. 1), for which the KS plot and autocorrelation plot for fitting the synthetic heartbeat data are shown in Fig. 4. As a comparison, the KS plot from the linear AR(14) model is also shown (see Fig. 4, left panel). It is worth noting that we have also simulated a linear Gaussian AR model for the R–R time series and have compared the performance between the linear and nonlinear predictive models—it was found that goodness-of- fit performance by the linear model is generally better than the one by a nonlinear model with the same model complexity (data not shown).

Figure 4
(Left) KS plots obtained from a nonlinear (10,4) model (solid line) and a linear AR(14) model (dotted line) for fitting the synthetic heartbeat interval data. (Right) Autocorrelation plot from the nonlinear (10,4) model. The dashed lines in all plots ...

For the two experimental datasets, we also conducted a preliminary model selection analysis (based on the AIC using the first 5-min recordings). Specifically, for testing the linear model alone, AIC analysis indicated p = 8 as the optimal linear order in almost all cases. In order to keep the number of unknown parameters relatively small, while the size of parameters from both linear and nonlinear models remain approximately the same (for fair comparison), we set p = 8 for linear modeling and empirically set p = 4 ~ 6 and q = 2 for nonlinear modeling. The fitting results are summarized in Tables II and III. Some representative tracking results are shown in Figs. 2 and and3,3, and their respective KS plots are illustrated in Figs. 5 and and6.6. In general, the results related to model fitting improvement vary among subjects in both datasets, and in some subjects we also found that the nonlinear model did not improve the goodness-of-fit compared with the linear model with equal model complexity.

Figure 5
KS plot comparison between linear (top) and nonlinear (bottom) modeling (subject 11, control, supine).
Figure 6
KS plot comparison for the healthy control subject 16483 (top linear versus bottom nonlinear).
Results From Selected Four Subjects in First Experimental Dataset.
Results From 28 Subjects (12 CHF and 16 Healthy) in Second Experimental Dataset

B. Nonlinearity

Model selection and goodness-of-fit tests on synthetic data validate the nonlinear quantification as evaluated by the point process framework, demonstrating that indexes, such as the ρ and ratio (defined in (14) and (15), respectively) are able to proportionally and correctly discern between series with linear and nonlinear prevailing dynamics. Based on the model with the best fit (i.e., the nonlinear model (10,4) in Table I), the ratio and the mean ρ value are computed as 0.40 and 0.49, respectively, for the simulated R–R series generated from the Rössler equations. For other simulated R–R time series with a linear Gaussian AR model, these values' estimates typically range from 0.91 to 0.99 (from various Monte Carlo simulations, data not shown).

In the first experimental dataset, the nonlinearity test showed that the level of nonlinearity varies from different postures and pharmacological conditions. For instance, all R–R time series in the control upright condition failed to reach significance in the nonlinearity test (see Table II). In addition, a higher presence of nonlinearity was observed when injecting ATR (control→ATR and PROP→DB), where parasympathetic modulation is absent, in contrast, lower nonlinearity (or higher linearity) was observed when injecting PROP (control→PROP and ATR→DB), where vagal activity is absent. Computation of the ratio and mean ρ statistics indicated that they typically had greater values in supine than in upright condition, suggesting that nonlinear interactions (LF: 0.01–0.15 Hz) become more prevalent due to the increase of cardiac sympathetic nerve activity and the reduction of vagal nerve activity. We did not observe consistent changes during the ATR or PROP administration.

In the second experimental dataset, from the results of the nonlinearity test (see Table III), it appeared that 15 out of 16 R–R time series from the healthy subjects showed significant nonlinearity (P <0.05), whereas in the CHF group, 5 out of 12 R–R time series failed to reach significance (test level 5%). Our test result confirms that the heartbeat dynamics from healthy subjects are more nonlinear to some degree. The fact that a lower degree of nonlinearity was found in the CHF patients suggests that pathological conditions might reduce the nonlinearity in the heartbeat interval series, which is also consistent with previous finding that a healthy heartbeat presents more pronounced nonlinear dynamics [5], [26], [27], [46], [62].

The nonlinearity effect from the second experimental dataset can also be observed in the computed time-averaging ρ index (within the LF range, 0.05–0.15 Hz). Generally, a time series with higher nonlinear dynamics would result in a lower ρ value, since the coefficients in the second-order Volterra terms would have relatively greater values [see (14)]. Note that the index ρ can be computed in a dynamic fashion, as the instantaneous estimate is obtained at each single time step. This is arguably more accurate than the batch estimate [ratio, defined in (15)], since the data are likely to be nonstationary. When the time series is close to be purely linear, both ρ and ratio will be close to 1. Statistical test (rank-sum test) show that these two indexes reveal significance differences between the CHF and healthy groups (P <0.002), with a higher level of nonlinearity in the healthy group (median ρ 0.9258 and median ratio 0.9564) than the CHF group (median ρ 0.9663 and median ratio 0.9995). It seems that the dynamic ρ index could detect higher level of nonlinearity possibly due to its more versatile nature than the ratio index. Particularly in the CHF group, the ratio indexes across all subjects seem to saturate to the level of pure linearity. Results in Table III present all of estimated statistics for each subject being analyzed. Fig. 7 shows two examples of tracking nonlinearity for one healthy control and one CHF subject.

Figure 7
Two tracking examples of dynamic nonlinearity index ρ(t) for one healthy control and one CHF subject. Note that the y-axis range of the second and fourth panels are different, where numbers in the boxes indicate the mean value averaged over time. ...

In comparing the two dynamic indexes in the second and fourth panels of Fig. 7, the more complex heartbeat structure in the healthy subject is revealed by substantial dynamic variations in the ρ index, associated with frequent marked increases in nonlinearity (note the drops of the ρ index below 0.85 at around 1200, 1750, 2050, and 2600 s). In contrast, a much less variable (and more uniformly linear) index is observed in the CHF patient. Clearly, the instantaneous index is useful in revealing different dynamic signatures of nonlinearity between healthy and CHF subjects that could not have been observed by using any other method.

C. Performance Comparison

The performance comparison between the proposed nonlinear Volterra modeling and the standard linear modeling was measured by the KS distance: the smaller the KS distance, the better the model fit. The comparative results are summarized in Tables II and III for the first and second experimental datasets, respectively. In the first experimental dataset, it is observed that nonlinear Volterra modeling generally improves the model fit, especially when the R–R time series exhibit more nonlinearity (i.e., smaller p-values in the nonlinearity test). Overall, among a total of 17 × 6 = 102 epochs, 87% (89/102) of fittings improved in terms of the KS statistic. The same observation can also be found in the second experimental dataset (improvement were found in 71% (20/28) subjects), although neither linear nor nonlinear models being tested, thus far has fully passed the KS test (i.e., all points within the 95% confidence bounds).

The imperfect fitting performance in both experimental datasets confirms that modeling real heartbeat dynamics (in a probabilistic sense) remains a challenging task. The lack of fit in some experimental heartbeat series may be due to the fact that the choice of the inverse Gaussian distribution is insufficient for characterizing the highly complex dynamics involved in these nonstationary heartbeat time series (which sometimes involve dramatic transient changes, or possibly, a sudden switch to a different physiological state), or it may be also due to the fact that we have not included cardiovascular covariates in this analysis. Meanwhile, increasing the model memory and including covariates might improve goodness-of-fit, but it might not necessarily improve the AIC statistics. Determining an optimal tradeoff between complexity and performance remains an issue that needs to be standardized by further research. In our experiences, choosing a proper probabilistic model and informative covariates is more crucial and effective than increasing the model order in improving the goodness-of-fit. To this extent, further investigation will be required to improve our model.

D. Quantification of Self-Similarity Via Scaling Exponent

Nonlinearity is often related to the complexity (regularity or predictability) of the random time series. For heartbeat time series, many nonlinearity measures, such as the ApEn [44], [45], sample entropy [34], [47], multiscale entropy [21], [57], and Poincaré plot [12], have been proposed to study the irregularity of the heartbeat [22], [52]. Specifically, complex dynamics have been observed in heartbeat interval series from healthy subjects [27], [46], and there have been growing interests in developing nonlinearity indexes able to characterize the irregularity of heartbeat dynamics in both healthy and pathological conditions [4], [30], [42]. Research effort has been largely devoted to characterizing such nonlinear behavior at different timescales using relatively short recordings [3]. In time-series analysis, detrended fluctuation analysis (DFA) is a method for determining the statistical self-affinity of a nonstationary signal [41], [42]. Essentially, DFA constructs a trend based on polynomial fitting to extract and quantify fluctuations at different time scales, which is useful for detecting long-range correlations in time series. Hence, it may reveal the fractal structure of time series that often appears to be a long-memory process with power-law decaying autocorrelation function. Specifically, let x(t) denote the time series of length N, whose fluctuation are to be studied, an integrated series y(k) is computed as follows:


where x denotes the sample average of {x(t)}. The resultant y(k) series is then divided into n-length subsequences without overlap, and for each subsequence, a linear regression yn is fitted against k. The root mean square error fluctuation between y(k) and the local linear trend yn(k) is defined as follows:


and the power-law behavior en [proportional, variant] nα is characterized by DFA [Footnote: As an outcome, DFA computes the scaling α-exponent that is similar to the Hurst exponent. Notably α = 0.5 indicates that x(t) is uncorrelated white noise (or y(k) is a random walk), α < 0.5 implies the signal x(t) is anticorrelated (i.e., negative correlation), α = 1 implies 1/f noise and long-range correlation, and α = 1.5 indicates Brownian noise (i.e., the integration of white noise)].

The advantage of DFA over other conventional methods (e.g.,spectral analysis) for estimating the fractal exponent α is that it is able to detect long-range correlation for nonstationary time series and also to avoid false detection of long-range correlation due to artifact [42]. Several authors have tried to apply DFA to the raw heartbeat interval analysis for characterizing its fractal-like scaling properties [38], [42], [43]. Here, it is our intention to investigate the estimate of this method computed from different timescales using the second experimental dataset. Specifically, we computed the DFA α-exponent, using both the original R–R time series and the estimated μRR(t), for each subject. The fine temporal resolution (5 ms) enabled us to reveal the fractal structure of the evenly spaced continuous-time signals without using any interpolation technique. We also computed the two scale DFA α1 and α2 exponents (assuming α1 ≥ α2) using the estimated μRR(t) (we did not compute two-scaled DFA exponents using raw R–R series due to the small number of sample points). The purpose of computing two-scale α-exponent is to investigate if there is a presence of the so-called “crossover” phenomenon reported in [42]. Our analysis showed that the healthy control subjects have a clear crossover point in scaling (i.e., the data points can be fitted better with two straight lines with two different slope parameters α1 andα2), whereas the crossover point is less obvious in the CHF subjects. An example of this finding on four subjects (two from each group) is illustrated in Fig.8.

Figure 8
Representative plots of log F (n) versus log n for computing the slope parameter α in DFA. (a) Two healthy control subjects. (b) Two CHF subjects. Note that in the healthy subjects, there appears a “crossover” phenomenon in scaling, ...

To evaluate the statistical differences between the healthy and CHF groups for the indexes computed from our method and DFA, we also conducted a nonparametric test under the null hypothesis that the medians of two sample groups are equal. The results are summarized in Table IV. As seen from the table, on average the CHF patients have lower HRV and greater HR. The insignificance of the α-exponent computed from the R–R time series might be due to the insufficiency of samples, in contrast, the α-exponent estimated from μRR(t) seems to be more accurate in characterizing the group difference (using 24 h recordings, it was reported in [42] that the scaling exponent statistic computed from raw R–R time series is 1.24 ± 0.22 for the CHF group and 1.00 ± 0.11 for the healthy group). Also, we found that using the μRR(t) estimates, the α1 andα2 exponents both show significant differences between healthy and CHF groups, with α2 appearing slightly more discriminating (P <1e-3) than α1 (P = 0.042). Furthermore, the gap between the α1 andα2 exponents is noticeably bigger in the healthy subjects, confirming the presence of the crossover phenomenon. Hence, the HR, HRV, ρ-index, and scaling α-exponent statistics can serve as useful metrics to distinguish the healthy and pathological conditions given relatively short heartbeat recordings. As an illustration, Fig. 9 presents a few scatter plots of selected mean estimated statistics between 12 CHF and 16 healthy subjects.

Figure 9
Scatter plots of some estimated statistics for 12 CHF (triangle) and 16 healthy (circle) subjects.
Mean±STD Group Statistics of Statistical Indexes Computed From Second Experimental Dataset.

E. Connection to Other Methods

As evidenced from our data analysis, the method presented here provides some new perspectives to characterize and measure the nonlinear dynamics of the heartbeat interval. Our quantification can also be combined with other nonlinearity indexes (e.g., the DFA) previously proposed in the field. Hence, our method can serve as a complementary tool for assessing the nonlinearity of heartbeat R–R time series. For instance, we can apply any desired nonlinearity index, such as the popular ApEn [44], to the estimated instantaneous μRR(t) series (see Table IV). ApEn is a statistical measure, which indicates the degree of randomness and regularity as applied to time series [10].As an example, we computed the ApEn(m, r) statistic for both the R–R interval series and the estimated μRR series. For the R–R interval series, we set r = 0.15 ~ 0.25, and m is either 1 or 2. The ApEn values in Table IV reflect the group mean ± STD statistics, whereas each subject's statistic is averaged from the estimates computed on a number of 5-min nonoverlapping segments. For the estimated μRR series, we set r = 0.1 ~ 0.15 and used the estimate at 100 ms timescale, so each 5-min segment contains 3000 samples for evaluating ApEn. As seen from Table IV, no statistical difference was found between two groups while using the raw R–R series. Varying the combination of parameters m and r did not change the qualitative statement about the group difference, although higher values in the mean ApEn statistic (computed with the raw R–R series) from the CHF group indicate a more erratic behavior in R–R interval series due to the loss of circadian variation [10]. In contrast, significant differences were found between the two groups (CHF versus healthy; rank-sum test, P <0.01) when using the estimated μRR series, interestingly, the mean ApEn statistics computed with the estimated μRR series are all higher in the healthy group. This suggests that the use of our instantaneous indexes examines different levels of irregularity and enhances the discriminating power of the ApEn metric between the pathologic and control groups.

IV. Discussion and Conclusion

We have presented a method for characterizing nonlinear dynamics of the human heartbeat within a point process paradigm. Unlike other nonlinear modeling methods developed in the literature, our point process probability model computes time-varying nonlinear indexes simultaneously with instantaneous HR and HRV statistics. Based on the second-order Volterra–Wiener series expansion, we have devised an adaptive point process filter to track the kernel coefficients and estimate the instantaneous parametric autospectrum and bispectrum, as well as the dynamic power ratio. It is noteworthy that it is also possible to incorporate a physiological covariate (such as respiration or blood pressure measures) into (2) of the point process model and produce further instantaneous indexes from their dynamic cross spectrum and cross bispectrum [17], [18]. Our model and method proposed here can be viewed as a further extension of previous models [8], [14], [16], which expands the horizon of modeling human heartbeat intervals.

Unlike other paradigms for estimating nonlinearity indexes developed in the literature [21], [26], [44], [57], our method is formulated within a probabilistic framework specifically developed for point process observations (the R–R intervals). Moreover, most other nonlinearity indexes are derived from nonparametric models, whereas our model is purely parametric and the analytically derived indexes can be evaluated in a dynamic and instantaneous fashion. We believe these strengths will enable our method as a useful tool for assessing nonlinear dynamics of heartbeat intervals in a nonstationary environment. Meanwhile, just like other approaches, our method also has caveats: besides the increased computational complexity, our model also requires sensible initial estimates of ξ and W, which might need to be reestimated from time to time in a dramatically changed environment.

Timescale is an important issue to evaluate the nonlinear nature of a physiological signal. It has been shown in [57] that different dynamical systems can exhibit similar nonlinearity signatures depending on the sampling time or sampling interval, and that similar systems can showdifferent degrees of nonlinearity by varying the timescale. This is naturally anticipated, since a purely irregular or nonlinear time series would have a specific range of dependence or correlation statistic in time [Footnote: For example, in the study of the chaotic Rössler time series used in Section V, it was observed that the ApEn statistic at different sampling rates varies quite a bit (range from 0.18 to 0.58, for the setup of m = 2 and r = 0.2), and the value does not monotonically change according to change in the sampling rate]. The point process framework provides us a reliable tool to examine the unevenly spaced heartbeat intervals at very high temporal resolutions, without resorting to any interpolation method. By performing a proper “up-sampling” using the estimate from our point process model, information can be discovered in a way that the original observed data cannot reveal. Furthermore, unlike other methods that might require large sample size (while directly operating on the raw R–R intervals), the point process method is potentially useful to examine short recordings of the physiological signals of interest. Certainly, the estimate at fine temporal resolution from the point process method is achieved at the cost of increasing computational complexity and requires a tradeoff approach by the practitioner.

The nonlinearity test provides a quantitative measure of the regularity of a tested time series [5], [9], [56]. In the study of the autonomic blockade protocol, it was found that the level of nonlinearity varies from different postures and pharmacological conditions, which essentially influence the sympathetic/parasympathetic balance in the autonomic nervous system. In comparing the healthy and CHF subjects, the heartbeat exhibits lower nonlinear dynamics in the pathological condition, which was confirmed by both the nonlinearity tests and the relative linear/nonlinear power ratio. These quantitative nonlinearity indexes can reveal statistical differences between groups with different cardiac conditions. We have also applied two well-established nonlinearity indexes: α-exponent and ApEn, to the estimated μRR series. Results under specific hypothesis testing reveal significant group differences between the healthy and CHF groups. Importantly, we have showed that by changing the timescale, our method can reveal different nonlinearity signatures. Further interpretation on the effect of timescale change on the nonlinearity indexes will be the subject of future investigation.

Finally, to conclude the paper, the probabilistic point process framework provides a new characterization for human heartbeat interval that allows us to estimate instantaneous indexes of HR and HRV, as well as indexes derived from the (linear) autospectrum and (nonlinear) bispectrum. Our experimental results in both synthetic and experimental heartbeat data, have demonstrated that our proposed point process model is useful in characterizing the inherent nonlinearity of the heartbeat dynamics. In the near future, we are planning to pursue this direction further in order to validate our new measures on more experimental studies, and investigate their potential use in real-time monitoring for clinical practice.


The authors would like to thank Dr. J. B. Schwartz (University of California, San Francisco) and Dr. G. B. Stanley (Georgia Institute of Technology and Emory University) for providing the experimental data used in this study. Recordings were performed during the tenure of a research fellowship from the American Heart Association, California Afliate (G. B. Stanley). The work was performed at the University of California, San Francisco, General Clinical Research Center. They also thank two anonymous reviewers for their valuable comments that helped improve the presentation of the paper.

This work was supported by the National Institutes of Health under Grant R01-HL084502, Grant R01-DA015644, Grant DP1-OD003646, Grant GM-26691, and Grant AG-9550, and by the Division of Research under Grant RR-79.


This work was presented at Proceedings of the IEEE Engineering in Medicine and Biology Conference (EMBC), 2008, Vancouver, BC, Canada.

Color versions of one or more of the figures in this paper are available online at


[1] Task Force of the European Society of Cardiology and the North American Society of Pacing Electrophysiology Heart rate variability. Circulation. 1996;93(5):1043–1065. [PubMed]
[2] Akay M, editor. Nonlinear Biomedical Signal Processing, Volume II: Dynamic Analysis and Modeling. Wiley-IEEE Press; New York: 2000.
[3] Armoundas AA, Ju K, Iyengar N, Kanters JK, Saul PJ, Cohen RJ, Chon KH. A stochastic nonlinear autoregressive algorithm reflects nonlinear dynamics of heart-rate fluctuation. Ann. Biomed. Eng. 2002;30:192–201. [PubMed]
[4] Atyabi F, Livari MA, Kaviani K, Tabar MRR. Two statistical methods for resolving healthy individuals and those with congestive heart failure based on extended self-similarity and a recursive method. J. Biol. Phys. 2006;32:489–495. [PMC free article] [PubMed]
[5] Baillie RT, Cecen AA, Erkal C. Normal heartbeat series are nonchaotic, nonlinear, and multifractal: New evidence from semiparametric and parametric tests. Chaos. 2009;19:028503-1–028503-5. [PubMed]
[6] Barahona M, Poon C-S. Detection of nonlinear dynamics in short, noisy time series. Nature. 1996;381:215–217.
[7] Barbieri R, Matten EC, Alabi AA, Brown EN. A point-process model of human heartbeat intervals: New definitions of heart rate and heart rate variability. Amer. J. Physiol. Heart Cicr. Physiol. 2005;288:424–435. [PubMed]
[8] Barbieri R, Brown EN. Analysis of heart beat dynamics by point process adaptive filtering. IEEE Trans. Biomed. Engin. 2006 Jan;53(1):4–12. [PubMed]
[9] Barnett AG, Wolff RC. A time-domain test for some types of nonlinearity. IEEE Trans. Signal Process. 2005 Jan;53(1):26–33.
[10] Beckers F, Ramaekers D, Aubert AE. Approximate entropy of heart rate variability: Validation of methods and application in heart failure. Cardiovasc. Eng. 2001;1(4):177–182.
[11] Brennan M, Palaniswami M, Kamen P. Distortion properties of the interval spectrum of IPFM generated heartbeats for heart rate variability analysis. IEEE Trans. Biomed. Eng. 2001 Nov;48(11):1251–1264. [PubMed]
[12] Brennan M, Palaniswami M, Kamen P. Do existing measures of Poincare' plot geometry reflect nonlinear features of heart rate variability? IEEE Trans. Biomed. Engr. 2001 Nov;48(11):1342–1347. [PubMed]
[13] Brown EN, Barbieri R, Eden UT, Frank LM. Likelihood methods for neural data analysis. In: Feng J, editor. Computational Neuroscience: A Comprehensive Approach. CRCPress; London, U.K.: 2003. pp. 253–286.
[14] Chen Z, Brown EN, Barbieri R. Proc. ICASSP'2008. Las Vegas, NV: A study of probabilistic models for characterizing human heart beat dynamics in autonomic blockade control; pp. 481–484. [PMC free article] [PubMed]
[15] Chen Z, Brown EN, Barbieri R. Characterizing nonlinear heartbeat dynamics within a point process framework. IEEE Annu. Conf.Eng. Med. Biol. (EMBC2008); Vancouver, Canada. pp. 2781–2784. [PMC free article] [PubMed]
[16] Chen Z, Brown EN, Barbieri R. Assessment of autonomic control and respiratory sinus arrhythmia using point process models of human heart beat dynamics. IEEE Trans. Biomed. Engr. 2009 Jul;56(7):1791–1802. [PMC free article] [PubMed]
[17] Chen Z, Purdon PL, Pierce ET, Harrell G, Brown EN, Barbieri R. Proc. ICASSP'2009. Taipei, Taiwan: Assessment of baroreflex control of heart rate during general anesthesia using a point process method; pp. 333–336. [PMC free article] [PubMed]
[18] Chen Z, Purdon PL, Pierce ET, Harrell G, Walsh J, Salazar AF, Tavares CL, Brown EN, Barbieri R. Linear and nonlinear quantification of respiratory sinus arrhythmia during propofol general anesthesia. IEEE Annu. Conf. Eng. Med. Biol. (EMBC2009); Minneapolis, MN. pp. 5336–5339. [PMC free article] [PubMed]
[19] Chon KH, Mullen TJ, Cohen RJ. A dual-input nonlinear system analysis of autonomic modulation of heart rate. IEEE Trans. Biomed. Engr. 1996 May;43(5):530–544. [PubMed]
[20] Christini DJ, Bennett FM, Lutchen KR, Ahmed HM, Hausdofi JM, Oriol N. Application of linear and nonlinear time series modeling to heart rate dynamics analysis. IEEE Trans. Biomed. Eng. 1995 Apr;42(4):411–415. [PubMed]
[21] Costa M, Goldberger AL, Peng C-K. Multiscale entropy analysis of complex physiological time series. Phys. Rev. Lett. 2002;89(6):068102-1–068102-4. [PubMed]
[22] Cysarz D, Lange S, Matthiessen PF, Leeuwen P. Regular heartbeat dynamics are associated with cardiac health. Amer. J. Physiol. (Regul. Integr. Comp. Physiol.) 2007;292(1):R368–372. [PubMed]
[23] Daley D, Vere-Jones D. An Introduction to the Theory of Point Processes, vol. I: Elementary Theory and Methods. 2nd Springer-Verlag; New York: 2003.
[24] de-Boer RW, Karemaker JM, Strackee J. Comparing spectra of a series of point events particularly of heart rate variability data. IEEE Trans. Biomed. Eng. 1984 Apr;BME-31:384–387. [PubMed]
[25] Glass L. Synchronization and rhythmic processes in physiology. Nature. 2001;410:277–284. [PubMed]
[26] Goldberger AL, Peng C-K, Lipsitz LA. What is physiologic complexity and how does it change with aging and disease? Neurobiol. Aging. 2002;23:23–26. [PubMed]
[27] Goldberger AL, Amaral LAN, Hausdorff JYM, Ivanov P, Peng C-K, Stanley HE. Fractal dynamics in physiology: Alterations with disease and aging. Proc. Nat. Acad. Sci. USA. 2002;99:2466–2472. [PubMed]
[28] Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov P, Mark RG, Mietus JE, Moody GB, Peng CK, Stanley HE. PhysioBank, physiotoolkit, and physioNet: Components of a new research resource for complex physiologic signals. Circulation. 2000;101(23):215–220. [PubMed]
[29] Heath R. Nonlinear Dynamics Techniques and Applications in Psychology. Lawrence Erlbaum Associates; Mahwah, NJ: 2000.
[30] Ivanov PC, Amaral LA, Goldberger AL, Havlin S, Rosenblum MG, Struzik ZR, Stanley HE. Multifractality in human heartbeat dynamics. Nature. 1999;399:461–465. [PubMed]
[31] Jo JA, Blasi A, Valladares EM, Juarez R, Baydur A, Khoo MCK. A nonlinear model of cardiac autonomic control in obstructive sleep apnea syndrome. Ann. Biomed. Eng. 2007;35(8):1425–1443. [PubMed]
[32] Korenberg MJ. Parallel cascade identification and kernel estimation for nonlinear systems. Ann. Biomed. Eng. 1991;19(4):429–455. [PubMed]
[33] Korenberg MJ, Paarmann LD. Orthogonal approaches to timeseries analysis and system identification. IEEE Signal Proc.Mag. 1991 Jul;8(3):29–43.
[34] Lake DE, Richman JS, Griffin MP, Moorman JR. Sample entropy analysis of neonatal heart rate variability. Amer. J. Physiol. (Regul. Integr. Comp. Physiol.) 2002;283(3):R789–R797. [PubMed]
[35] Loader C. Local Regression and Likelihood. Springer-Verlag; New York: 1999.
[36] Marmarelis PZ. Identification of nonlinear biological systems using Laguerre expansions of kernels. Ann. Biomed. Eng. 1993;21:573–589. [PubMed]
[37] Mateo J, Laguna P. Improved heart rate variability signal analysis from the beat occurrence times according to the IPFM model. IEEETrans. Biomed. Eng. 2000 Aug;47(8):985–996. [PubMed]
[38] Morren G, Lemmerling P, Daniëls H, Naulaers G, Van Huffel S. Sensitivity of detrended ßuctuation analysis applied to heart rate variability of preterm newborns. IEEE 27th Annu. Conf. Eng. Med.Biol. (EMBC2005); Shanghai, China. pp. 319–322. [PubMed]
[39] Nikias CL, Mendel JM. Signal processing with higher-order spectra. IEEE Signal Proc. Mag. 1993 Jul;10(3):10–37.
[40] Nikias C, Petropulu AP. Higher Order Spectra Analysis: A Non- Linear Signal Processing Framework. Prentice-Hall; Englewood Cliffs, NJ: 1993.
[41] Peng C-K, Buldyrev SV, Havlin S, Simons M, Stanley HE, Goldberger AL. Mosaic organization of DNA nucleotides. Phys. Rev. E. 1994;49:1685–1689. [PubMed]
[42] Peng C-K, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos. 1995;5:82–87. [PubMed]
[43] Perfetto JC, Ruiz A, Attellis CD'. Detrended fluctuation analysis (DFA) and R-R interval variability:Anewlinear segmentation algorithm. Proc. Comput. Cardiol. (CinC2006) :629–632.
[44] Pincus SM. Approximate entropy as a measure of system complexity. Proc. Nat. Acad. Sci. USA. 1991;88:2297–2301. [PubMed]
[45] Pincus SM, Goldberger AL. Physiological time-series analysis: What does regularity quantify? Amer. J. Physiol. (Heart Circ. Physiol.) 1994;266:H1643–H1656. [PubMed]
[46] Poon C-S, Merrill CK. Decrease of cardiac chaos in congestive heart failure. Nature. 1997;389:492–495. [PubMed]
[47] Richman JS, Moorman JR. Physiological time series analysis using approximate entropy and sample entropy. Amer. J. Physiol. (HeartCirc. Physiol.) 2000;278(6):H2039–H2049. [PubMed]
[48] Rompelman O, Snijders J, Van Spronsen CJ. The measurement of heart rate variability spectra with the help of a personal computer. IEEETrans. Biomed. Eng. 1982 Jul;BME-29(7):503–510. [PubMed]
[49] Rössler OE. An equation for continuous chaos. Phys. Lett. 1976;35A:397–398.
[50] Schmitz A, Schreiber T. Testing for nonlinearity in unevenly sampled time series. Phys. Rev. E. 1999;59:4044–4047.
[51] Shibaa M, Kikuchia A, NMiaob T, Haraa K, Sunagawaa S, Yoshidaa S, Takagia K, Unnoc N. Nonlinear analyses of heart rate variability in monochorionic and dichorionic twin fetuses. Gynecol. Obstet. Invest. 2008;65(2):73–80. [PubMed]
[52] Signorini MG, Ferrario M, Marchetti M, Marseglia A. Nonlinear analysis of heart rate variability signal for the characterization of cardiac heart failure patients. IEEE Annu. Conf. Eng. Med. Biol.(EMBC2006); New York. pp. 3431–3434. [PubMed]
[53] Stanley GB, Verotta D, Craft N, Siegel RA, Schwartz JB. Age and autonomic effects on interrelationships between lung volume and heart rate. Amer. J. Physiol. (Heart Circ. Physiol.) 1996;270(5):H1833–H1840. [PubMed]
[54] Sugihara G, Allan W, Sobel D, Allan KD. Nonlinear control of heart rate variability in human infants. Proc. Nat. Acad. Sci. USA. 1996;93:2608–2613. [PubMed]
[55] Sunagawa K, Kawada T, Nakahara T. Dynamic nonlinear vagosympathetic interaction in regulating heart rate. Heart Vessels. 1998;13:157–174. [PubMed]
[56] Theiler J, Eubank S, Longtin A, Galdrikian B, Farmer JD. Testing for nonlinearity in time series: The method of surrogate data. PhysicaD. 1992;58(1–4):77–94.
[57] Thuraisingham RA, Gottwald GA. On multiscale entropy analysis for physiological data. Physica A. 2006;366:323–332.
[58] Tulppo MP, Kiviniemi AM, Hautala AJ, Kallio M, Seppanen T, Makikallio TH, Huikuri HV. Physiological background of the loss of fractal heart rate dynamics. Circulation. 2005;112(3):314–319. [PubMed]
[59] Westwick DT, Kearney RE. Identification of Nonlinear Physiological Systems. Wiley; New York: 2003. Explicit least-squares methods; pp. 169–206.
[60] Wu G-Q, Arzeno NM, Shen L-L, Tang D-K, Zheng D-A, Zhao N-Q, Eckberg DL, Poon C-S. Chaotic signatures of heart rate variability and its power spectrum in health, aging and heart failure. PLoS ONE. 2009;4(2):e4323. [PMC free article] [PubMed]
[61] Zhang Y, Wang H, Ju KH, Jan K-M, Chon KH. Nonlinear analysis of separate contributions of autonomous nervous systems to heart rate variability using principal dynamic modes. IEEE Trans. Biomed.Eng. 2004 Mar;51(2):255–262. [PubMed]
[62] Special issues on nonlinearity on heart rate. Chaos. 2009;19