|Home | About | Journals | Submit | Contact Us | Français|
The pacemaking system of the heart is complex; a healthy heart constantly integrates and responds to extracardiac signals, resulting in highly complex heart rate patterns with a great deal of variability. In the laboratory and in some pathological or age-related states, however, dynamics can show reduced complexity that is more readily described and modeled. Reduced heart rate complexity has both clinical and dynamical significance – it may provide warning of impending illness or clues about the dynamics of the heart’s pacemaking system. In this paper, we describe uniquely simple and interesting heart rate dynamics that we have observed in premature human infants -- reversible transitions to large-amplitude periodic oscillations -- and we show that the appearance and disappearance of these periodic oscillations can be described by a simple mathematical model, a Hopf bifurcation.
The normal, healthy heart beats at a variable rate with extraordinarily complex fluctuations across a wide range of time scales (1–9). This phenomenon is commonly referred to as heart rate variability (HRV). Multiple mechanisms contribute to HRV; it is a composite reflection of autonomic outflow (the combination of sympathetic and parasympathetic nervous system activation), neuroendocrine influences, and autonomic responsiveness (the ability of the cardiovascular system to respond to changes in autonomic outflow) (2). It has been demonstrated that this healthy complex variability deteriorates as a consequence of disease and aging (1,2,5,7). Complexity in heart rate may be easy to see but it is hard to define, so there are many different measures of complexity, having varying degrees of utility in different contexts (10–22).
Regularity, on the other hand, is often easy to see and easy to define precisely, especially when that regularity has the form of periodic cycles. Two periodic cycles of heart rate have been extensively studied. Respiratory sinus arrhythmia, the coupling of heart rate to breathing, produces fluctuations in heart rate with a period of a few seconds in adults, or about one second in infants in a neonatal intensive care unit (NICU) (23,24). Another cycle of heart rate is correlated with a cycle of blood pressure called Mayer waves (25,26); the period of these is about ten seconds in adults or 12–14 seconds in NICU patients (23). In these infants, who typically have interbeat intervals near 350 ms (and heart rate around 170 beats per minute), the magnitude of the changes of interbeat intervals due to respiratory sinus arrhythmia is about 2–4 ms, while that associated with Mayer waves is about 10 ms. Combined, they can sum to give periodic fluctuations of magnitude around 5% of the interbeat interval. Other, lower frequency components elevate the range of fluctuations of heart rates further (6), but these are less cyclical.
We describe here in infants a different and previously uncharacterized heart rate cycle with large decelerations. A deceleration is a decrease in heart rate followed by a return to the base rate. In previous work, we have found that heart rate decelerations in the setting of otherwise low heart rate variability often preceded acute neonatal illness (27–33). In those studies, we used statistical measures based on heart rate distributions to make predictive models for clinical use.
For this new work, we devised a wavelet-based detector of decelerations in heart rate records and applied it to a large clinical database. We find that large decelerations are common, similar in shape among infants, and can appear in clusters in which they are sometimes remarkably periodic for epochs as long as two days. Their period is nearly constant at about 15 seconds, and their sizes may increase the interbeat interval from 350 to 550 ms, or reduce the heart rate from 170 to 110 beats per minute. Thus these cycles of decelerations are much too long and much too large to represent respiratory sinus arrhythmia, and, though the period is comparable, they are much larger than reported Mayer waves in infants (23). Indeed, these cycles do not resemble any well-characterized heart rate oscillation, and we propose that they constitute a previously unrecognized form of low-dimensional heart rate dynamics. We give a mathematical interpretation of clusters of periodic decelerations as resulting from a Hopf bifurcation, a common way for a steady state to change to a cycle.
This phenomenon is interesting from both dynamical and clinical perspectives. It is dynamically interesting because it shows that the immature control system of the heart rate can go into a previously uncharacterized oscillatory mode. These observations therefore can provide a new point of contact with mathematical models of the control system (34–36). It is clinically interesting because the occurrence of clusters of decelerations is correlated with risk of neonatal sepsis1.
We give here our methods of computerized detection and measurement of decelerations. We collected heart rate (HR) data from neonatal intensive care units, and extracted interbeat (RR) intervals as previously described (27,28,31–33). (All the rhythms were sinus in origin, and there were no variations in the P-wave morphology. Also there was little variation in the PR intervals; in fact the PR intervals got a little shorter as the RR intervals became longer.)
Some examples of neonatal RR interval time series are shown in Figure 1. Panel (a) is a complex time series of a healthy NICU infant; panel (b) is a time series showing reduced variability of a NICU infant with sepsis (bloodstream infection); panel (c) is a time series with decelerations (each spike in the RR interval time series is a deceleration, or a slowing of the HR followed by return to a base rate); and panel (d) is a time series with periodic decelerations. Compared to case (a), case (d) plainly represents reduced-dimensional dynamics. The period of the decelerations is nearly constant at about 45 beats, or 15 seconds, and, while the decelerations have a variety of heights (which can be anywhere from 20 to 300 ms), they have a common shape.
To investigate the dynamics, we developed a neonatal HR deceleration detector. We studied the HR records of 479 very low birth weight (< 1500 g) neonates from NICUs at the University of Virginia, the University of Alabama at Birmingham, and Wake Forest University, with a total of 513,193 half-hour HR records collected between August 2005 and May 2006. To search these records, we constructed a pattern-matching algorithm inspired by wavelet theory (37) – we defined a template function that gives an adequate description of the shape of a deceleration, and swept this function through the RR interval time series to identify points at which the HR signal matched the template within a specified tolerance.
To detect decelerations in the signals that also contain fluctuations that have the appearance of random noise, we chose a function χb (n − n0), which describes the observed shape of typical decelerations, and we treated the RR interval (time interval between beats) signal as if it were a sum of that deceleration function plus some remainder:
Here RR(n) is the time between beat n and beat n+1, χb (n − n0) is a function describing a deceleration of width b centered about the point n0, and G(n) represents the remainder of the signal (Fig. 2). The value a represents the height or amplitude of the deceleration and is to be determined. We get an estimate of this height from the formula,
The numerator is the projection of a portion of the signal RR(n) onto the waveform of the deceleration:
To employ this scheme effectively, an appropriate template function, χ, is needed. Fig. 3 shows several decelerations from the record shown in Fig. 1d superimposed atop one another. Evidently, decelerations are relatively symmetric in shape and, compared with a Gaussian, are narrower at the peak and wider in the wings. A function that meets these criteria and subsequently has proven to successfully locate decelerations of similar nature throughout signals from many infants is:
where the amplitude a is calculated as described by Equation (2), and b is the width parameter. The goodness of fit of this function was high; when compared with visually identified decelerations the median R-squared value was 0.93. In addition, the χ(n) function consistently provided a reasonable estimate of the height of a deceleration relative to its baseline. Thus, χ(n) was used as a template for a deceleration detector2.
We swept this template function χb (n − n0) through the signal RR(n) with width parameter b ranging from eight beats to 100 beats and calculated the correlation coefficient a(n0, b) at each scale and translation. This process generated a surface of a(n0, b) values like that shown in Fig. 4. We then found the points (n0, b) at which a(n0, b) has a local maximum, that is, where the correlation between the deceleration waveform and the signal was locally strongest.
Having isolated these (n0, b) points, we added several more stipulations in order to identify decelerations. First, if there was more than one local maximum near a particular translation, we picked the largest of the maxima (highest coefficient a) so that one deceleration was not erroneously represented by two waveforms. Next, to avoid spurious correlations, we demanded that the wavelet χb (n − n0) fit the original signal with an R-squared value of at least 0.75.
We found that decelerations may be isolated, or, less commonly, they may occur in clusters lasting up to two days. (In a separate, clinically-oriented paper, we will show that clusters of decelerations are correlated with impending sepsis.) Fig. 5a shows a half-hour excerpt of such a cluster, the whole of which lasted approximately one day. Within such extended clusters of decelerations, we sometimes found shorter intervals of time, lasting up to several hours, in which the decelerations showed remarkable periodicity (Figure 5c) – we refer to these as periodic sequences of decelerations. While decelerations are common, clusters of decelerations are uncommon3, and periodic sequences of decelerations lasting 10 minutes or more without interruption (we refer to these as long periodic sequences of decelerations) are rare; they were seen in six cases out of the population of several hundred patients. In each of these cases, the period of the long periodic sequences of decelerations was about 15 seconds. We focus on periodic and long periodic sequences of decelerations for the remainder of this report.
A cluster of periodic decelerations having a common shape suggests low-dimensional dynamics. Moreover, the clusters began and ended abruptly. Interestingly, small amplitude oscillations of similar period were often present in the vicinity of these clusters (Figures 5c, ,7c).7c). All of these behaviors are consistent with a noisy Hopf bifurcation; in fact, a Hopf bifurcation is the most general way a system can change from stable to oscillatory. Such bifurcations occur for example in laser systems, oscillatory chemical reactions, predator-prey dynamics, and in the Hodgkin-Huxley model of the giant squid axon (38–42). The effect of noise on a system having a Hopf bifurcation was described by Wiesenfeld and his collaborators (43,44).
The theory of Hopf bifurcations is well known, and we summarize just enough of it here to show how we do our simulations. One assumes that the pacemaking system of the heart has feedback loops that can be modeled by a large number of dynamical variables (denoted u=[u1…un]) governed by an equally large number of differential equations4,
These governing equations are assumed to contain many parameters p=p1…pm, which may vary with time in the following way: on short time-scales, they may have small, rapid, noisy fluctuations but generally they do not have large changes; on long time-scales, they have substantial slow variation. The functions f (u;p), governing the rate of change of u, are not presumed to be known. However, it is assumed that these functions have Taylor expansions that converge in the range of interest, and that these functions have a zero (steady state) which can be taken to occur at u=0. For some values of the parameters p, the steady state is assumed to be stable.
If the Taylor expansions are truncated at first degree, one obtains a set of linear equations
The eigenvalues of the matrix M(p) associated with these linear equations must either be real or occur in complex conjugate pairs. When all real parts of these eigenvalues are negative, the steady state is stable. The steady state can go unstable if, as the parameters p change, one real eigenvalue, or a pair of complex-conjugate eigenvalues, crosses the imaginary axis. Hopf theory examines the latter case, and provides two powerful theorems. (1) “Center Manifold Theorem”: In the state space of dynamical variables u, there is a two-dimensional surface (called a center manifold) which is an invariant surface and an attractor. That means: (a) if u(t) lies initially on this two-dimensional surface, it stays on this surface; (b) if u(t) lies initially off the surface, it moves toward the surface in the manner of exponential decay. Furthermore, the surface is analytic (Taylor-expandable) and the evolution in the surface can be described by two differential equations. (2) “Normal Form Theorem”: there exists an analytic change of variables to new coordinates (x, y) such that the governing pair of equations can be reduced to a standard “normal” form:
In polar coordinates, that normal form is
The new parameters μ, a, b, and ω depend on the parameters p. If, as the parameters p vary with time, μ changes from negative to positive, the stable steady state goes unstable, and one of two things happens: (1) if a<0, a stable cycle is created (this is a “soft” Hopf bifurcation, also called “supercritical”); (2) if a>0, an unstable cycle is destroyed; in this case, if b<0, the radius r(t) can go quickly to a large value while the angle θ(t) increases steadily, so the system “jumps” to a large cycle of frequency ω (this is a “hard” Hopf bifurcation, also called “subcritical”). Finally, in both cases (1) and (2), if μ remains just below the bifurcation point, small oscillations of frequency ω will occur in the presence of noise – these are called “noisy precursors”.
We can use this theory to simulate the heart rate if we only postulate that the time between beats RR(t) is some function of one of these unknown variables x(t)=r(t)cos(θ(t)). Since we know the shape of the decelerations, we can use a Fourier representation of that shape to get RR(r(t), θ(t)):
To generate a simulation, we integrated the equations (7), adding random noise to both variables after every time step, Δt. This noise term was multiplied by a coefficient, . The factor ξ controls the strength of the noise and, for our simulations, generally assumed a value between .01 and 0.1. The factor ensures that the statistical properties of the noise fluctuations are independent of the step size5 Δt. The integration was carried out for various values of parameters μ, a, and b (though a was always kept positive and b negative in order to keep the bifurcation “hard”) and for cases in which these parameters fluctuated in time; more specifically, in our model, we set μ initially negative, and allowed it to vary near zero. This integration process yielded output like that shown in Figure 6a; as μ goes from negative to positive, the system goes from low variability to noisy precursors to large oscillations. As μ then decreases, the system crosses another critical point, in this case at μ=−1/18, and returns to its low variability state. The hysteresis shown here is characteristic of a hard Hopf bifurcation. This output in (x(t), y(t)) was then transformed into RR(t) as defined by Equation (9), using only the first 11 terms in this series, as these are sufficient to provide a good approximation to the deceleration shape. This transformation converted the output shown in Figure 6a into that shown in Figure 6b.
Thus, by allowing the parameter μ to vary in the Hopf model and using a Fourier representation to transform the output into an RR interval series, it is possible to produce clusters of periodic decelerations similar to those observed in the clinical data (Fig. 7). Furthermore, by specifying the times at which decelerations begin and end, and forcing the parameter μ to cross its critical points accordingly, one can reasonably capture the essential elements of the HR data including the noisy precursors, abrupt starts and stops, and repeating cycles of decelerations (Figure 7c).
Another way to show the cycle is in Cartesian coordinates (x, y)= (r cos θ, r sin θ), or, better, in coordinates (RR, y). Part of the simulation shown in Fig. 7c is shown in these variables in Fig. 7d.
We conclude that this two-dimensional dynamical model, based on minimal assumptions, gives a credible simulation of periodic neonatal HR decelerations.
Using a new wavelet-based deceleration detector, we have described an apparently new type of periodic behavior in the heart rates of neonates. While oscillatory patterns have been observed previously in the human heart rate in the form of Mayer waves or Cheyne-Stokes dynamics (2,45), the oscillations reported here are up to ten times larger in amplitude than are seen in those cases. We interpret the behavior as a noisy Hopf bifurcation, and we have shown that quantitative output of this model resembles the observations. If this interpretation is correct, then we are reporting here for the first time the observation of a previously unknown bifurcation in the dynamics of the control system of the human heart rate.
Decelerations like those we observe, and in particular long periodic sequences of decelerations, have been seen only in neonates. Therefore one may speculate that the bifurcation producing them is connected with the immature control system of these infants. An obvious candidate for further study is the feedback loops between heart rate, blood pressure, respiration, and oxygen saturation. On the observational side, we are now collecting large databases consisting of waveforms of EKG, oxygen saturation, respiration, temperature and in some cases blood pressure, for some hundreds of infants in two NICU’s. On the mathematical side, we are examining models of the physiological control loops governing heart rate. We know that some of these models (developed from data on adults) exhibit dynamical bifurcations (35,46), and we will see whether appropriate modifications of such models give behavior that is consistent with our data.
We emphasize that since long periodic sequences of decelerations occurred spontaneously in the clinical setting, and were not induced by controlled means, they appear to be a natural mode of low-dimensional dynamics in the human cardiac pacemaking system near the time of birth.
This work was supported by NIGMS 64640. A.A. Flower was supported by National Heart, Lung, and Blood Institute Basic Cardiovascular Research Training Grant 5T32-HL-07284. JBD is supported by NSF grant PHY-0457070. Computing support was provided by the SciClone Cluster Project at the College of William and Mary. We thank Dr. Gene Tracy for his helpful suggestions.
1Sepsis (bloodstream infection) is one of the greatest threats to infants in a NICU, but diagnosis of neonatal sepsis using current clinical and laboratory methods is difficult, and has a high rate of false negatives.
2We do not attach any theoretical significance to this function; rather, it has been reached empirically as a function that represents the data.
3Large decelerations are common: at least one occurs in about 1/3 of our 20–30 minute records. If we define a “cluster of large decelerations” as a collection of at least four decelerations, each having height of at least 100 ms, then such clusters are uncommon; of the approximately 450,000 half-hour neonatal heart records studied, less than 700 contained four or more decelerations of 100 ms height or more.
4Hopf theory can also describe differential equations containing time-delays, or difference equations describing beat-to-beat models.
5The point here is that the width parameter of the distribution resulting from a random walk is proportional to the size of each step times the square root of the number of steps.
Required conflict of interest statement: Medical Predictive Science Corporation of Charlottesville, Virginia, has a license to market technology related to heart rate characteristics monitoring of newborn infants, supplied partial funding for this study, and aided in collection of data. They played no role in study design, analysis, and interpretation of data, or of the writing of the report or the decision to submit the paper for publication. Drs. Lake and Moorman have an equity share in this company.