HEART RATE VARIABILITY AS DETERMINISM WITH JUMP STOCHASTIC PARAMETERS ∗

We use measured heart rate information (RR intervals) to develop a one-dimensional nonlinear map that describes short term deterministic behavior in the data. Our study suggests that there is a stochastic parameter with persistence which causes the heart rate and rhythm system to wander about a bifurcation point. We propose a modified circle map with a jump process noise term as a model which can qualitatively capture such this behavior of low dimensional transient determinism with occasional (stochastically defined) jumps from one deterministic system to another within a one parameter family of deterministic systems.


1.
Introduction. Modeling the behavior of human cardiovascular system is an interesting problem which draws extensive attention from researchers. A fundamental and challenging question is how best to provide a simplified representation of both the deterministic and stochastic aspects of heart dynamics. Suder et al. [1] recorded RR intervals while restricting the paced respiration cycle lengths above 8 s, from which they observe that heart rate variability obeyed a dynamic rule that can be expressed by a one-dimensional, nonlinear circle map. Later, Jason et al. [2] showed that even during spontaneous breathing (with subject at rest), the angle component of RR-interval still has a highly deterministic structure after filtering out low frequency components. Shiau et al. [3] argued that the deterministic characteristic could result from sympathetic activation and thermo-regulation, which are primarily evidenced in the low frequency component of the RR-interval signal. They used a simple nonlinear noise-reduction method [4] to remove the high frequency component and used the next angle map to reconstruct a deterministic attractor.
Our approach is based on Shiau's work, but we use cubic smoothing spline to remove the high frequency component of RR-intervals time series data. This filter approach (via spline) was selected (1) for ease of implementation 1 , but also (2) because the smoothing spline enforces a continuity, consistent with our expectation that heart rate variability should vary continuously . We find the resultant data can 2010 Mathematics Subject Classification. Primary: 15A15, 15A09, 15A23. Key words and phrases. Heart rate variability, electrocardiography, next angle map, circle map, jump process.
The first author is supported by NSF grant DMS-0708083. 1 Smoothing spline requires choice of a single smoothing parameter, while Schreiber's noisereduction methods requires optimizing over a pair of parameter values. Also note that we can obtain similar results, which will be discussed, using Low-Pass filtering by FFT. be easily related to a modified sine circle map, but with an interesting noise behavior which we describe as a stochastic parametric factor with persistence. Here we define 'persistence' as piecewise constant in time, with occasional random jumps. The stochastic behavior that we find is similar in character to that identified in Lerma's recent study on the stochastic aspects of Cardiac Arrhythmias [5], which argued that in the neighborhood of bifurcation points, the fluctuations induced by the stochastic opening and closing of individual ion channels in the cell membrane results in membrane noise that may lead to randomness in the observed dynamics of cardiac rhythm systems. In [6], Kuusela et al. give a simple one-dimensional Langevin-type stochastic difference equation which can model the heart rate fluctuations in a time scale from minutes to hours. The similarity between to our work is that both provide a stochastic model which aims at uncovering the interaction of determinism and stochastic control of cardiac dynamics. However, we focus our study on the low frequency component of RR-intervals data, which relates to sympathetic and vagal activity of heart, to reveal the stochastic jump with persistence around a bifurcation point.
The paper is organized as follows: Section 2 discusses our filtering approach and our decision to focus on the low frequency (LF) component of the Heart Rate Variability signal (HRV) in the phase space reconstruction. Section 3 discusses the phase reconstruction from the ECG data and discusses the complex behaviors observed in the data. In Section 4, we give a reasonable model to simulate these observed behaviors, with conclusions in section 5.
2. Spectral components of HRV. Over the past two decades, the general body of research has recognized a significant relationship between the autonomic nervous system and cardiovascular mortality, including sudden cardiac death. Although cardiac automaticity is intrinsic to various pacemaker tissues, heart rate and the rhythm are largely under control of the autonomic nervous system [7]. Different frequency ranges of the HRV have been related to various physiological phenomenons [7]. Studies of spectral components of short term recordings of HRV show that the efferent vagal activity is a major contributor to the high frequency (HF) component, while the LF component is considered as a marker for sympathetic modulations, with some studies also suggesting that LF reflects both sympathetic and vagal activity [7].
In Suder's experimental method [1], they obtain a one-dimensional, nonlinear deterministic observable from the HF component by restricting the respirationcycle to greater than 8s. Their theoretical foundation is that controlling respiration can induce an increase in the HF signal. After high pass filtering, they identify a one-dimensional deterministic process. Similarly, Janson's work [2] obtains a lowdimensional structure by extracting the HF component of HRV, but without control on respiration. In contrast, Shiau [3] obtains a one-dimensional deterministic process from the LF component under spontaneous breathing conditions. Regardless of the frequency domain being studied, there is evidence showing that low-dimensional deterministic processes describing the autonomic nervous system can be observed, with the low-dimensional models providing insight into heart rate dynamics [8,9]. In this paper, we focus on the LF component while generating a one-dimensional map, which corresponds to sympathetic and vagal activity [7]. A primary reason for choosing to analyze the LF component is that it is easier to reduce the effects of noise introduced by measurement error.

Electrocardiogram data.
From an open source repository (physionet.org), we obtained RR-intervals data taken from volunteers who were supine and asked to breathe at a fixed rate of 0.25 Hz for 10 min. (Data originally from [10].) The twodimensional embedding of the raw data, shown in Figure 1(a), has no apparent low dimensional structure. To identify a low-dimensional model, we apply the following sequence of processing steps: First, we filter out the HF component by applying a cubic smoothing spline interpolation to the RR-intervals data, as illustrated in Figure 1(b), an alternative filtering technique to that of [4].  As second step, for each filtered data point, we compute an angular coordinate representation φ n measured relative to the centroid of the data set, RR : We then construct a time delay embedding of the angular coordinates to produce the "next angle map." Figure 2(a) shows the time delay embedding representation of the next angle map (1) data, which appears to be reasonably well described by a one dimensional curve, but with data lying nearly on, but both above and below, the identity line on that graph. We find that there exist two typical cases for the portion of the trajectory that lies near the intersection between identity line and the data: • Case 1: The trajectory evolves along a lower branch, below the identity line, resulting in decreasing values; • Case 2: The trajectory evolves along an upper branch, which intersects the identity line with two fixed points, generating increasing values. Figure 2(b)) illustrates these two behavior. We observe from our data that the trajectory appears to persist along a branch (either upper or lower) for several iterations before "jumping" to the other branch. Of special importance is that although   (1) indicates that a onedimensional representation may be a reasonable approximation.
(b) Cobweb along typical data trajectories. Highlighting three trajectory segments near the line φ n+1 = φ n . We note that the trajectory sometimes travels following a lower branch (below the line and decreasing -case 1) while at other times, it follows an upper branch (above the line and increasing-case 2).
the full data set may appear as a "cloud" as if it were noisy observations, any particular trajectory follows a smooth (seemingly deterministic) path for several iterates. The noisy cloud appearance results from overlaying numerous deterministic partial trajectories, each with slightly different parameters. We describe this behavior as "persistence," and view it as parametric noise. Describing this phenomena as noise is not meant to imply that the underlying process is truly stochastic, but is simply a recognition that even if it is deterministic, we are not modeling the complex control system behavior, which would be affected by the body condition, circumstance, emotion,etcetera. . We label '1' to the case when the trajectory follows the lower branch, with red patches. The corresponding interval of RR n is concave down above the average line or concave up and below average; and '2' to the case when the trajectory follows the upper branch, with blue patches. The corresponding interval of RR n is downward cavity when below average or upward concavity when above average. (Bottom) Series data of φ. Decreasing φ n relates to the patch '1'; increasing φ n relates to the patch '2'. For both graphs, the data is discrete, with the curve drawn for clarity of illustration.
For convenience, the following discussion uses '1' as shorthand notation for case 1 trajectories, and '2' for case 2. In 3, we plot the sequence data of RR n and φ(n). From (1), we can understand the relationship between these data for each monotone interval of φ(n) : Case 1 -φ(n + 1) < φ(n), or φ(n) is decreasing, with data below the identity line of Fig. 1(b) . The corresponding interval for RR n would be the one where RR n is concave down above the average line or concave up and below average (See Fig. 3 for the red patches).
Case 2 -φ(n + 1) > φ(n), or φ(n) is increasing, data above the identity line in Fig. 1(b). The corresponding interval of RR n is downward cavity when below average or upward concavity when above average (See Fig. 3 for the blue patches).
These observed behaviors lead to a natural symbolic labeling and transition graph representation using symbols '1' and '2.' In this context, the "persistence" of the system to stay on a particular branched for several iterates would be evidenced in the transition diagram ( Figure 4) as p ii > 1/2, with larger values indicating greater persistence. We remark that our data and case description admits some ambiguity that might appear in longer dataset: (1) it is possible that the trajectory following some branch may experience a parameter 'jump', but land on a new branch that is the same case as before the jump, and that jump would not result in a transition to the other node, (2) the trajectory might jump to a Case 2 branch (with a stable fixed point) while the system state is above the fixed point, in which case the trajectory would decrease toward the stable fixed point, with this behavior essentially indistinguishable from a Case 1 trajectory. Figure 4. Transition relationship between upper and lower branch trajectories, as observed in the data. When p ii is large, the system will tend to "persist" along a either the upper or lower branch.
In this paper, we are not trying to provide a physiological explanation for this behavior. Rather, we simply have observed the phenomena. In the next section, we build a model based on the circle map [11,12] which can mimic the features discussed above through representation as a stochastic process with an interesting form of 'noise' in that process. 4. The circle map. The circle map is a generic term describing a family of dynamical systems whose state space can be interpreted as angles of a circle. A simple example of circle map using the modulo function and is given by φ n+1 = f (φ n ) = ω + φ n (mod T ), where ω and T are constant. A second example is the sine circle map, which was introduced by Kolmogorov and well studied by Vladimir Arnold [11,12], is defined by where ω and k are constants. Figure 5 plots the Sine Circle Map with ω = 5 and k = 1.25 as the red curve, while the plotted points (green) are from the data. The circle map provides a prototypical model for systems that are controlled by two pacemakers, with the term k sin φ describing the effect of the nonlinear oscillator coupling [12], which is why it has been previously studied as a model for the heart [11]. When k = 1.2 trajectories would mimic the Case 1 trajectories of our data. As k is increased, the map gets closer to the identity line, yields trajectories with slow passage between the identity line and the curve [13], a ghost of the fixed point associated to the saddle node bifurcation that occurs at approximately k = 1.283. When k = 1.35, the curve has crossed through the identity line. Two fixed points are created, one is stable, the other is unstable. Trajectories on this map would mimic Case 2 behavior during the transient approach toward the stable fixed point.
In order to simulate similar behaviors of our ECG data using the circle map, we proceed as follows: We assume that k is piecewise constant (the persistence) with changes to k occurring as Poisson arrivals. When k changes, the new value for k is chosen as an i.i.d. variable from a Gaussian distribution: where 0 ≤ p c ≤ 1, p n ∼ U [0, 1], and n ∼ N [0, 1]. Critical value p c governs how often k n changes its value. (The Gaussian scale parameters as well as the choice p c = 1/17 were selected heuristically to provide a good visual match to the observed data trajectory.) A visualization of such k n is given in Figure 6(b). Whenever k n < 1.283 (below the red horizontal line in Figure 6(b)) the dynamics are similar to case 1 of our data. For k n > 1.283, the behavior is like case 2. Figure 7(a) shows an example of a typical trajectory generated by our model. The cobweb around the crossover section is plotted in Figure 7  . Whenever k n is above the red horizontal line, the sine circle map is behaves similarly to case 2 of our data; whenever k n is below the red horizontal line, the sine circle map is behaving similarly to case 1 of our data; segments, illustrative reasonable qualitative agreement with our data. We note that long sequences of '2' behavior results in the nearly horizontal part of φ n in Figure  7(a). We similarly give a graph representation for the modified circle map in Figure  8. For our model, we can compute the transition probabilities as follows: Define a k-refresh as the condition that k n = k n−1 , the situation where the system has jumped to a new parameter value. Then the conditional probability that refresh results in an "upper branch" (Case 2) behavior is given by is independent of system state and is easily computed from the normal distribution function. For the particular parameters described above, this yields p u ≈ 0.361. The transition from '1' to '1' is due to either no refresh of parameter or a refresh, but with new k n < 1.283, so that The other transition probabilities can be easily computed in a similar fashion.

5.
Conclusion and discussion. In summary, we apply cubic smoothing spline to the RR-intervals data to remove the HF component, and then we extract the angles coordinate. A delay embedding of the data indicates that a low-dimensional representation is reasonable, but analysis of structure in the time series indicates that a better representation can be achieved by viewing the heart rate as driven at a parameter point near a saddle node bifurcation. Stochastic fluctuations (with persistence) cause the system to wander back and forth through that bifurcation. To model this behavior which derives from combination of deterministic and stochastic factors, we propose a stochastically perturbed modified circle map. The model   shows qualitative agreement with the two-dimensional embedding of angles from heart rate data in the sense of simulating the stochastic behaviors. We remark that our choice to model the stochastic behavior of parameter k as piecewise constant with Poisson arrival of jumps was based on the simplicity of the approach. It is certainly reasonable to also consider the case where k experiences small fluctuations (rather than constant behavior between jumps), in which case a standard jumpdiffusion model might be appropriate. However, our noisy data did not appear to be sufficient to resolve that low level diffusion.
In this paper, we do not address the question why this modified circle is "good" to model this heart rate data, leaving as open question the issues of what physiological effects might generate such behavior. Moreover, the question of "goodness," from a mathematical framework, requires that we have some way to quantify deviation of model and system. The degree to which a "toy model" might be representative of a more complicated system is a fundamental issue from dynamical systems that is not easily resolved. For example, comparisons between dynamical systems based on L p spaces may fail to be a good judgment if the systems turn out to be chaotic. Mostly Conjugacy [15] appears to be a promising method to address this problem, because it compares dynamical systems in a way of judging the quality of "matching" by looking at their topological difference (homeomorphic defect). Additionally, our choice of specific model parameters was based primarily on a qualitative assessment. Using tools from the theory of commuters [15], tailored to this stochastic setting, may allow for a reasoned way to choose parameters to best match system dynamics. Comparing the models and the heart rate data by using Mostly Conjugacy method will be our next work on this topic. We suggest that applying these parameter estimation techniques to data collected during regular physical examination of heart could provide a new method for automated "change detection" with respect to cardiac function, where one might hope that detection of such changes in dynamics might have clinical relevance.
Acknowledgments. This work was accomplished with funding provided by National Science Foundation grant DMS-0708083 . We would also like to acknowledge recommendations of the reviewers which have significantly improved this manuscript.
Appendix A. Test for determinism using surrogate data. A primary concern when conducting analysis of processed data is that observed structures may be an artifact of the processing, with no relevance to the true phenomena under study. To show that the low dimensional structure of successive angles does not simply result from the low-pass filtering, we compare our observations with surrogate data created by linear Gaussian stochastic process which preserves both the spectrum and the histogram of the empirical RR-interval data [16]. We process that surrogate data using the same filtering methods, as described in Section 3. Figure 9 plots the Angle map based on data from the experimentally measured data [10] compared with processed surrogate data.  In the "eyeball" metric, treating the data as a cloud of points in R 2 , we observe little qualitative difference between measured and surrogate data. However, as a dynamic process, we see distinguishing features. As described in Section 4, in our data, we observe that trajectories (in delay coordinates shown in Figure 9) that are near the line φ(n) = φ(n + 1) that either track upward (when above the diagonal) or downward (when below) with trajectories only infrequently "crossing" the line. Examining the surrogate data, generated from a Gaussian stochastic process with the matching spectrum and distribution as the data, we find frequent crossings of that main diagonal after processing the surrogate data, as represented by the illustration in Figure 10 . In particular, if focus on that portion of trajectories lying near the diagonal (in delay space) where φ n+1 = φ n , and project data onto that line (as illustrated Figure 10(R)). We test as follow: We take as null hypothesis that the time series of RR-interval is generated by a linear stochastic process. Let sampled data c be the projection of the selected (select the points around near the main diagonal in phase space) onto the diagonal. We then generate 10 surrogate data sets {s1, s2, ...s10} for the RR interval data, where we model as a guassian process with the same power spectrum and distribution as the empirical RR-interval data, using method iAAF T [17]. Let cc = [s1 s2 ... s10] bet the concatenation of all 10 surrogate data sets. Then cc can be viewed as the "average" behavior of surrogate data sets. We perform a Kolmogorov −Smirnov test on c and cc. The test results in the rejection of the null hypothesis at the 5% significance level, with p − value = 3.4037e − 004.
Our test is meant to establish that the low dimensional structure of RR-interval data does not come from a Gaussian stochastic process. In some sense, this test seems to argue for greater structure to the underlying process, reasonably captured by our model of deterministic transient behavior, though the formal result is simply rejection of that null hypothesis.