Abstract
We present a systematic approach for detecting nonlinear components in heart rate variability (HRV). The analysis is based on twentythree 48h Holter recordings in healthy persons during sinus rhythm. Although many segments of 1,024 RR intervals are stationary, only few stationary segments of 8,192–32,768 RR intervals can be found using a test of Isliker and Kurths (Int. J. Bifurcation Chaos 3:1573–1579, 1993.). By comparing the correlation integrals from these segments and corresponding surrogate data sets, we reject the null hypothesis that these time series are realization of linear processes. On the basis of a test statistic exploring the differences of consecutive RR intervals, we reject the hypothesis that the RR intervals represent a static transformation of a linear process using optimized surrogate data. Furthermore, time irreversibility of the heartbeat data is demonstrated. We interpret these results as a strong evidence for nonlinear components in HRV. Thus RR intervals from healthy persons contain more information than can be extracted by linear analysis in the time and frequency domain.
 RR intervals
 stationarity
 nonlinear dynamics
 surrogate data
 linear models
various linear parameters of heart rate variability have been used for risk stratification of patients following myocardial infarction (1, 10), however, with limited success (4, 13). The search for robust parameters continues in an effort to improve the characterization of highrisk patients (12, 23), since therapy for primary prevention has to consider the implantable cardioverter defibrillator (15). It is our hope that risk stratification can be improved by taking into account all, i.e., also the possible nonlinear components of the heartbeat signal.
The time series of RR intervals between successive sinus beats that can be derived from a Holter recording show in general a complex behavior (see Fig. 1), which suggests the presence of nonlinearities. The objective of our study is to present a systematic approach for detecting nonlinear components in heart rate variability to define ways for their routine evaluation. In refinement of earlier studies (6, 7), we systematically search for stationary data segments to exclude trends that may simulate nonlinear properties. In addition, we use optimized methods to generate surrogate data sets (19), which are fundamental in the proof of nonlinear components.
We demonstrate in this study that heartbeat signals contain more components than can be assessed by spectral analysis or corresponding parameters in the time domain. A characteristic nonlinear feature is the time irreversibility of RR intervals, which can simply be measured using the skewness of the differences of consecutive RR intervals.
METHODS
Data Collection
Holter recordings of 48 h (two consecutive days) were obtained from 23 healthy subjects (mean age: 36.9 yr, range: 21–75 yr, 16 male and 7 female) who maintained their usual daily activities. For electrocardiogram recording and analysis, the ELATEC system in connection with the 2448 Holter recorder (ELA Medical, Munich, Germany) was used, allowing for a precision of ±2.5 ms for the recognition of the onset of QRS complexes. All intervals between 300 and 1,800 ms except for premature beats were entered into the analysis. The average percentage of premature beats was <1% in the stationary segments.
Stationarity
Most statistical measures of a time series are only meaningful if the time series can be regarded as stationary. According to a definition given by Takens (22), a time series {RR
_{n}} is stationary if for each k and for each continuous function gthe average
For l we used all possible powers of 2 starting at each consecutive 1,024beat segment allowing for extension of the window onto the second day of the recording. According to Isliker and Kurths (5), we estimate the probabilitiesp
_{k}(S
Surrogate Data Sets
For modeling the linear components in the heartbeat time series, we consider the RR intervals as output of a linear system.
Linear models.
Model A is the simplest conceivable linear system (called a linear autoregressive Gaussian process), where a Gaussian white input ε_{i} (i = 1, … , N) yields Gaussian distributedRR
_{1}, … ,RR
_{N} as given by
Model B is an alternative model consisting of a linear process (see model A) where the white input ε_{i} does not have a Gaussian distribution. In this case, RR _{n} given byEq. 5 is in general not Gaussian distributed. The ε_{i} can be chosen such that the resulting distribution of RR _{i} fits the measured data. These surrogate data sets are used to reject the null hypothesis that the data are realizations of the linear autoregressive process with arbitrary input.
Model C is a more realistic linear model and is defined as output of a linear autoregressive Gaussian process, where the imaginary
Generating Surrogate Data
Linear modeling (model B).
For the stationary segment of our time series {RR
_{i}} of RR intervals, we calculate an optimal linear predictor. A linear predictor of orderk
Phase randomization (model C).
We compute the discrete Fourier transform of our original data, which consists of an amplitude and a phase at each frequency. We then shuffle the data for destroying all correlations and compute the discrete Fourier transform of the shuffled data. Now we replace each amplitude by the amplitude of the same frequency of the original data. After taking the inverse Fourier transform, we adjust the amplitude of the surrogate data by applying a nonlinear transform to give the surrogate data the same distribution as the original data. This completes the first iteration step and here usually the algorithms for generation of surrogate data stop (3, 8, 24). Following a suggestion by Schreiber and Schmitz (19), we start a second step by calculating the discrete Fourier transform of our surrogate data and again replace the amplitudes and so on. This iteration scheme allows us to produce surrogate data, which have the same distribution as the original data and almost the same power spectrum. This power spectrum is in better concordance to the power spectrum of the original time series than the power spectra of traditionally produced surrogate data (19). The surrogate data, generated in this fashion, are output of a linear Gaussian process (model C).
Correlation Integral
The correlation integral is a measure to describe nonlinear characteristics of a time series. In this context, it is used to detect differences between the original time series and the surrogate data set computed with the method described under linear modeling. For the surrogates we calculate mean and variance of the correlation integrals. In short, the kdimensional vectors Rec_{i} = (RR
_{i},RR
_{i+1}, … ,RR
_{i+k−1}), where i = 1, … N − k + 1, of a time series are called kdimensional reconstruction vectors. For a bounded stationary and infinite time series {RR
_{i}} the correlation integral C
_{k}(x) is defined as the probability that two randomly and independently chosen reconstruction vectors are within distancex. An unbiased and minimal variance estimate for this correlation integral, based on an independent sample of Nreconstruction vectors Rec_{i}, is given by the following formula
The correlation integralsC _{k}(x) will be estimated for an embedding dimension of k = 10, because depending on the individual recording, a significant difference to the correlation integrals of linear models can be observed only at higher kvalues. For statistical comparison, x is chosen ς({RR _{n}}), where ς({RR _{n}}) denotes the standard deviation of the RR intervals in the stationary segment.
Time Irreversibility
We define a stationary time series to be time reversible if the probability of the vector (RR _{n}, … ,RR _{n+k−1}) is equal to the probability of the vector (RR _{n+k−1}, … , RR _{n}) for all k and n (2). The reversibility of linear Gaussian processes has been shown by Tong (25) and obviously holds also for all static transformations of a linear Gaussian process. Thus all realizations of model C are reversible.
It follows for a timereversible segment that the return mapRR
_{n+1} versusRR
_{n} is symmetrical with respect to the line of identity. A quantitative analysis of the irreversibility of a time series is provided by the asymmetry of the distribution of the difference of consecutive intervals with respect to zero. An appropriate scale invariant test statistic, which is frequently used to detect deviations from time reversibility (18, 20) is given by
RESULTS
Stationarity
Figure 1 A shows the sequence of RR intervals from a 24h record, which contains a particularly long sleeping period (beats 30,000–70,000). Typically, the average cycle length of 1,200 ms during this period is interrupted by brief decreases in cycle length (down to 500 ms). These heart rate spikes are known to occur during major body movements during sleep (11). A remarkable long stationary segment of 16,384 heartbeats (marked region) occurs in the sleeping period. This segment is depicted in Fig. 1 B. Stationarity of this segment is demonstrated in Fig.2, where the estimated probability density of the RR intervals of the whole segment as well as the first half is plotted. The two longest stationary segments had a length of 32,768 consecutive RR intervals. They were found in two patients, one in each. The length of stationary segments could not be increased when analyzing 48 h of the recordings. The heartbeat time series of 11 subjects had a single stationary segment of 16,384 RR intervals corresponding to 4–5 h. Whereas many stationary segments of 1,024 consecutive RR intervals were found in all cases, only a few stationary segments of 8,192 consecutive heartbeats were present in each 21 of 23 subjects. Figure 3,A and B, indicates the percentage of stationary nonoverlapping segments of a fixed length of 1,024 and 4,096 consecutive RR intervals, respectively. It is evident from Fig. 3, A and B, that in general the variability between the subjects is greater than the variability between the first and the second day for the same subject.
The bins used for estimating the local probability distributions varied in number depending on the individual RR interval distributions. With regard to 1,024 consecutive RR intervals, for instance, the number of bins ranged between 25 and 50.
The stationary periods usually occurred during the nighttime. From 12:00 to 6:00 AM, 70% of the recording contains stationary segments of 1,024 RR intervals. Figure 4 shows the distribution of the stationary segments for two consecutive 24h recordings.
Surrogate Data With Identical Linear Properties
Model A.
The probability density of the RR intervals is not Gaussian due to the skewness of the estimated probability density (see Fig. 2 for an example). The skewness differs significantly from the skewness of a Gaussian distribution for 86% of all stationary segments (P< 0.05). This implies that modeling the RR intervals with a simple linear autoregressive model (model A) with Gaussian noise with identical variance to the noise in the original time series is in general not justified, since these linear models have a Gaussian distribution of their values.
Model B.
To exclude simple rescaling of the original data, we compared the variances of surrogate data sets and of the original data. The variances of the original data were in the 2ς range of the variances of the corresponding surrogate data sets for all subjects with a stationary segment of 8,192 consecutive RR intervals.
The correlation integral of the stationary segment in Fig. 1 differs from the correlation integral of the surrogate data sets generated according to model B, as shown in Fig.5. Here we calculated the correlation integral for 200 different values of x (equals distance of the reconstruction vectors) for the 100 surrogate data sets.
In 19 of 23 subjects, the hypothesis that the data are a realization ofmodel B could be rejected because at x = ς({RR _{n}}), theZ value was > 4. Z was even >10 in 12 subjects. Only in two persons were the correlation integrals within the 2ς range for all values of x, which leads to aZ value <2 in these cases. In the remaining two subjects, no stationary segment of 8,192 intervals could be found.
Model C.
All return maps of stationary RR intervals have shown asymmetry with respect to the line of identity. An example of this is given in Fig.6. For a quantitative analysis, the distribution of the test statistic R _{1} is plotted in Fig. 7 A. Here, the distribution of R̂ _{1} for surrogate data sets (model C) is compared withR̂ _{1} of the original set. The test shows that the null hypothesis, i.e., the time series are the output of linear Gaussian processes, can be rejected with high significance in almost all cases (Fig. 7 B). The number of standard deviations ς by which R̂ _{1} differs from the mean of R̂ _{1} of the surrogate data was <2 only in 6 of 42 stationary segments containing 8,192 RR intervals.
DISCUSSION
Our analysis is based on the comparison of stationary original data with surrogate data sets generated by optimized methods. Surrogate data were produced using three different models to preserve the linear properties, while all nonlinear components are destroyed. Various test statistics, including correlation integrals and a measure of time irreversibility, unambiguously showed a difference between the surrogate data and the original stationary segments, demonstrating nonlinear components in the latter.
Stationarity
For a reliable estimate of the correlation integral, the time independence of the probability distribution must be guaranteed for the embedding (5). This is best achieved by the test for stationarity suggested by Isliker and Kurths (5), since this test is based on the probability distribution. Other tests for stationarity perform less well in this regard, since they are usually confined to the first and second moments (8, 17). The probability distribution of the onedimensional reconstruction vectors is equal to the probability distribution of the RR intervals. For practical purposes, however, it is not advisable to estimate higher dimensional probability distributions, because the statistical uncertainties increase with higher dimensions due to the scarcity of data. The estimation of the correlation integral of higher dimensions is more robust in this regard, since the correlation integral is a function of scalars and the probability distribution is a function of highdimensional vectors.
Our study shows that only relatively short segments of 24h recordings can be considered stationary as a result of the stationarity test proposed by Isliker and Kurths (5). In particular, during the daytime the underlying biological rhythms may change in a way that characteristic quantities cannot be regarded time invariant. Use of stationary segments may even help to separate different nonlinear dynamic states (6). It remains to be determined whether much longer recordings, e.g., 12day recordings (14), contain longer stationary segments. This may be of considerable importance for the description of highdimensional systems.
Rejection of Linear Models
By comparing the correlation integrals from stationary heartbeat time series and their corresponding surrogate data sets, we can reject the null hypothesis that these time series are realizations of a linear process with nonnormal input. By comparing the incremental distributions of stationary heartbeat series and their corresponding surrogate data sets, we can also reject the null hypothesis that these time series are realizations of a linear Gaussian process with a nonlinear measurement function. Because nonstationarities can be excluded as reason for positively discriminating statistics, we interpret the results as evidence for nonlinear components in our data.
We have no formal proof for nonlinearities, since linear models other than autoregressive models have not been evaluated. However, the models chosen here represent the whole class of linear processes that can be described by linear autoregressive models. These models have no structure other than the autocorrelation coefficients (22) or what is equivalent, the power spectrum. Despite these theoretical limitations, a conclusion can be drawn for practical purposes: The power spectrum does not describe all structures in heartbeat time series. Moreover, consideration of nonlinear analysis may provide more sensitive methods for assessing physiological states (16, 21, 26).
Correlation integral.
In comparison to the clinically established time domain parameters pNN50 (percentage of differences between adjacent NN intervals >50 ms), SDSD (standard deviation of successive differences), RMSSD (rootmean square of successive differences) (23), which measure the variability of two consecutive RR intervals, the correlation integral measures the variability of sequences of d consecutive RR intervals, where d is the dimension of the reconstruction vectors. The fact that the original time series have greater values of the correlation integral than the surrogate data reflects smaller distances between the reconstructed vectors in the original time series. This implies that the beattobeat variability of the linear model is larger than in the actual data. The return map (Fig. 6) of an original data set and a representative surrogate time series illustrate this phenomenon. In particular at short cycle lengths (<900 ms), the surrogate data show higher beattobeat variability.
Thus the correlation integral appears to be a robust measure to detect nonlinear components in the data. This, however, does not hold true for the correlation dimension, which is derived from the correlation integral (7). The computation of a correlation dimension assumes a power law behavior of the correlation integral for a significant scaling region. This could not be found in our data. In addition, there was no convergence of the “slopes” for higher embeddings.
Time irreversibility.
Time irreversibility is known to be a powerful indicator of nonlinearities (2, 18, 20). Here we have shown the time irreversibility of stationary segments with high significance. The fact that theR _{1} statistics yielded zero values for the surrogate data does by itself not imply time reversibility of the surrogate data. However, because of inherent properties of time series created by model C, time reversibility has been shown (25). Thus methods that explore the time direction in the data can be used for a more comprehensive analysis of RR intervals. The widely used parameter pNN50 and other clinically established heart rate variability measures such as RMSSD, SDSD, and NN50count (23) do not differentiate about the relation between upgoing and downgoing changes in heart rate and therefore do not contain any information about time direction. The same holds true for spectral analysis of the heartbeat signal, since spectral analysis of the data in reverse order leads to same results than the original time series.
Comparison to Other Studies
Our findings confirm earlier reports on the use of surrogate data to detect nonlinearities (6, 7, 9). To our knowledge, this is the first application of the improved method for generating surrogate data from stationary heartbeat time series. By the almost perfect match of the surrogate data to our original data in terms of spectrum and histogram, we can exclude trivial falsepositive test statistics for nonlinearities. In conjunction with the selection of stationary segments of up to 5 h, our approach is more rigorous. In addition, our approach is novel in the use of the correlation integral and the time irreversibility as measures for nonlinearities in heart rate variability.
Limitations
As shown in Fig. 3, A and B, stationary segments often comprise only relatively small segments of a 24h recording period. Therefore, this approach leads to results that are restricted to these episodes and are not simply applicable to longer segments. Thus methods should be sought to deal with nonstationary segments. The daily activity of the subjects was not controlled during recording, but this is not achieved in hospitalized patients. A natural standardization occurred during the night time.
Our approach does not show that the process cannot be due to a linear process with nonnormal noise that is passed through a static nonlinearity, i.e., a combination of our models B andC. But in this scenario there are nonlinear components in the data due to the admittedly trivial static nonlinear transformation of the output.
In conclusion, the analysis presented here provides strong evidence for nonlinear components in heart rate variability. We demonstrate that the heart rate variability signal contains components that cannot be assessed by spectral analysis. The search for new parameters comprising these nonlinear components will lead to a more complete description of heart rate variability. The physiological and prognostic value of this nonlinear approach remains to be established.
Acknowledgments
We thank Dr. Thomas Schreiber, University of Wuppertal, and Dr. Rainer Scharf, University of Leipzig, for helpful discussions.
Footnotes

Address for reprint requests: M. Meesmann, Medizinische Klinik der Universität Würzburg, JosefSchneiderStr.2, 97080 Würzburg, Germany.

This study was supported by the Bundesministerium für Bildung und Forschung, Germany, with additional support by St. Jude Medical GmbH Ventritex, Leverkusen, within the project “Nichtlineare EKGAnalysen zur Risikostratifizierung und Therapiebeurteilung von Herzpatienten.”
 Copyright © 1998 the American Physiological Society