NIH Public Access
`Author Manuscript
`IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
`Published in final edited form as:
`IEEE Trans Biomed Eng. 2010 June ; 57(6): 1335–1347. doi:10.1109/TBME.2010.2041002.
`
`Characterizing Nonlinear Heartbeat Dynamics within a Point
`Process Framework
`
`Zhe Chen [Member, IEEE]
`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 (zhechen@neurostat.mit.edu)
`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 (enb@neurostat.mit.edu)
`Riccardo Barbieri [Senior Member, IEEE]
`Neuroscience Statistics Research Laboratory, Harvard Medical School, Massachusetts General
`Hospital, Boston,MA02114 USA, and also with the Massachusetts Institute of Technology,
`Cambridge, MA 02139 USA (barbieri@neurostat.mit.edu)
`
`Abstract
`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
`
`Correspondence to: Zhe Chen.
`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 http://ieeexplore.ieee.org.
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`APPLE 1034
`
`1
`
`

`

`Chen et al.
`
`Page 2
`
`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.
`
`IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`2
`
`

`

`Chen et al.
`
`Page 3
`
`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
`
`(1)
`
`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 ∊ [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.
`
`IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`3
`
`

`

`Chen et al.
`
`A. Heartbeat Interval
`
`Page 4
`
` detected from the ECG, let RRj = uj −
`Suppose we are given a set of R-wave events
`uj−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 t−uj 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
`. 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,…,RRt−h), where RRt−j 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:
`
`(2)
`
`. Here the coefficients a0(t),{ai(t)}, and {bkl(t)} correspond to
`where
`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)].
`
`IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`4
`
`

`

`Chen et al.
`
`Page 5
`
`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(t−uj)−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(t−ut)−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
`
`(3)
`
`(4)
`
`(5)
`
`where
`
` and
`
`.
`
`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]
`
`(6)
`
`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:
`
`(7)
`
`(8)
`
`IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`5
`
`

`

`Chen et al.
`
`Page 6
`
`(9)
`
`(10)
`
`where P and W denote the parameter and noise covariance matrices, respectively; Δ =0.005s
` denote the first- and second-
`denotes the time bin size; Δλk = ∂λk/∂ξk and
`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 ) of the point process observations
`[13], [23]
`
`(11)
`
`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 ℓk = nk logλk−λk approximates the
`instantaneous log-likelihood function ℓ(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 Kolmogorov–Smirnov
`(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
`
`IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`6
`
`

`

`Chen et al.
`
`Page 7
`
`{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
` 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:
`
`, 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
`
`(12)
`
`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,
` denote the Fourier transform of the second-
`let
`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]
`
`(13)
`
`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,
` (the second equality follows from
`
`IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`7
`
`

`

`Chen et al.
`
`Page 8
`
`. Let b(t) denote a vector that contains
`the conjugate symmetry property), or
`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:
`
`(14)
`
`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]:
`
`(15)
`
`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.
`
` are both cycles/
`It shall be noted that the frequency units appearing in
` and
`beat, since the autoregression of μRR (t) is conducted on previous beats instead of previous
`instantaneous {μRR (t–i)}. 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.
`
`IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`8
`
`

`

`Chen et al.
`
`Page 9
`
`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.
`
`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].
`
`V. Data
`
`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).
`
`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.
`
`The second heartbeat dataset, which was retrieved from a public source: Physionet
`(http://www.physionet.org/) [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
`
`IEEE Trans Biomed Eng. Author manuscript; available in PMC 2011 June 1.
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`NIH-PA Author Manuscript
`
`9
`
`

`

`Chen et al.
`
`Page 10
`
`healthy subject. Since these recordings have longer durations, they have been deemed as
`particularly suitable for studying complex heartbeat interval dynamics [42], [46].
`
`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).
`
`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 3, and their respective KS plots are illustrated in Figs. 5 and 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.
`
`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 (