INTRODUCTION

Top Abstract Introduction Methods Results Discussion References

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 high-risk 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 R-R 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 R-R intervals, which can simply be measured using the skewness of the differences of consecutive R-R intervals.

	METHODS

Top Abstract Introduction Methods Results Discussion References

Data Collection

Stationarity

(1)

i₀ + l  1

(2)

For l we used all possible powers of 2 starting at each consecutive 1,024-beat segment allowing for extension of the window onto the second day of the recording. According to Isliker and Kurths (5), we estimate the probabilities p_k(S₀^all) and p_k(S₀^half)

(3)

for S = S₀^all and S = S₀^half. Then we compare the estimated probability distributions by a ²-test to see whether distributions are equal. The test statistic is defined as

(4)

This test statistic can be expected to have a ²-distribution with m  1 degrees of freedom. To achieve a proper estimate of the probability distribution, the initial bin width was 5 ms in accordance to the limited resolution of our data (±2.5 ms). We used five points as a lower bound for the number of data points in a bin. If necessary, neighboring intervals were grouped.

Surrogate Data Sets

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 distributed R-R₁, ... , R-R_N as given by

(5)

n = 1, ... , N, where the coefficients a_i are uniquely determined by the power spectrum of the data. Because the data never had a Gaussian distribution (see Fig. 2 for example), this model could be rejected without simulations.

(6)

Generating Surrogate Data

Linear modeling (model B). For the stationary segment of our time series {R-R_i} of R-R intervals, we calculate an optimal linear predictor. A linear predictor of order k

(7)

for a given time series {R-R_i} is called optimal if the average value of the square of the prediction error is minimal. The residue r_i is defined by r_i =   R-R_i, i = 1, ... , N. We denote the coefficients of the linear predictor by a₀, ... , a_k. The residues r_i, together with the linear predictor, are then used to produce a surrogate time series {s_n} as follows

(8)

where n = k + 1, ... , N. The start values s_n, (n = 1, ... , k) are R-R_n, (n = 1, ... , k). The term µ_n is chosen randomly from the set {r_i} to destroy all components of the original time series that are not captured by the coefficients a₀, ... , a_k. In this way, we obtain a surrogate time series {s_n}, which has the same autocorrelation coefficients, at least for time lags i < k, as the original time series. The residues are not automatically normally distributed but adapted from our original time series. By repeating this process, we create a set of 100 surrogate data sets, all of which have almost the same linear properties as the original time series. Following a suggestion by Takens (22) to take a rather high order, we set k to 25 for all analyses. These surrogate data sets allow us to check whether the stationary segment of R-R intervals can be regarded as the realization of a linear process (model B).

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

i+k1

(9)

(10)

The correlation integrals C_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 k values. For statistical comparison, x is chosen ({R-R_n}), where ({R-R_n}) denotes the standard deviation of the R-R intervals in the stationary segment.

Time Irreversibility

n+k1

n+k1

It follows for a time-reversible segment that the return map R-R_n+1 versus R-R_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

(11)

(E denotes the expectation value), which can be estimated by

(12)

Because the differences of R-R intervals have zero mean, R₁ is the skewness of the differences of R-R intervals.

	RESULTS

Top Abstract Introduction Methods Results Discussion References

Stationarity

View larger version (52K): [in this window] [in a new window]   Fig. 1.   A: R-R intervals of a 24-h recording of a healthy person (subject 6). Recording started at 5:39 PM. There is a prolonged sleeping period between beats 30,000 and 70,000. B: stationary segment as indicated in A.

View larger version (22K): [in this window] [in a new window]   Fig. 2.   Probability density of R-R intervals from the stationary period of Fig. 1. Thick line shows probability density of R-R intervals of whole segment; dashed line shows probability density of R-R intervals in the first half of this segment. Probability density is estimated by dividing the frequency in a bin by the bin width.

View larger version (17K): [in this window] [in a new window]   Fig. 3.   A: percentage of stationary segments of 1,024 consecutive R-R intervals. Each point represents one subject. Variability between subjects is greater than variability between the first and the second day for same subject. B: percentage of stationary segments of 4,096 consecutive R-R intervals. This percentage is much smaller for all subjects than percentage of stationary segments of 1,024 consecutive R-R intervals.

The bins used for estimating the local probability distributions varied in number depending on the individual R-R interval distributions. With regard to 1,024 consecutive R-R 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 R-R intervals. Figure 4 shows the distribution of the stationary segments for two consecutive 24-h recordings.

View larger version (97K): [in this window] [in a new window]   Fig. 4.   Stationary segments (grey regions) of 1,024 consecutive R-R intervals of a person (subject 6) at 2 consecutive days. Night phase contains more stationary segments in both recordings.

Surrogate Data With Identical Linear Properties

Model A. The probability density of the R-R 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 R-R 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 R-R intervals.

View larger version (21K): [in this window] [in a new window]   Fig. 5.   Correlation integrals from original time series and surrogate data set. Dashed line corresponds to average correlation integral of 100 linear autoregressive models of order 25 with adapted noise (C10, model B). Solid lines correspond to 2-range of correlation integrals. Dotted line shows correlation integral of stationary segment of Fig. 1. Dimension of reconstruction vectors is 10. Correlation integrals were compared at x=156, which is the standard deviation of the stationary segment.

Model C. All return maps of stationary R-R 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₁ is plotted in Fig. 7A. Here, the distribution of ₁ for surrogate data sets (model C) is compared with ₁ 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. 7B). The number of standard deviations  by which ₁ differs from the mean of ₁ of the surrogate data was <2 only in 6 of 42 stationary segments containing 8,192 R-R intervals.

View larger version (30K): [in this window] [in a new window]   Fig. 6.   Two-dimensional return map of stationary segment (A) (Fig. 1) and a surrogate data set (B). Return map of stationary segment shows more structure than return map of surrogate data set. In addition, asymmetry with respect to line of identity is obvious.

View larger version (15K): [in this window] [in a new window]   Fig. 7.   A: test statistic 1 of 100 surrogate data sets produced by iterative scheme (model C) and of stationary segment of R-R intervals. B: histogram of number of standard deviations , test statistic 1 of stationary segment differs from mean of 1 of the corresponding surrogate data sets. Only first stationary segment of 8,192 R-R intervals of each 24-h records is plotted in the histogram.

	DISCUSSION

Top Abstract Introduction Methods Results Discussion References

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

Our study shows that only relatively short segments of 24-h 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., 12-day recordings (14), contain longer stationary segments. This may be of considerable importance for the description of high-dimensional systems.

Rejection of Linear Models

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 N-N intervals >50 ms), SDSD (standard deviation of successive differences), RMSSD (root-mean square of successive differences) (23), which measure the variability of two consecutive R-R intervals, the correlation integral measures the variability of sequences of d consecutive R-R 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 beat-to-beat 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 beat-to-beat variability.

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 the R₁ 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 R-R 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 up-going and down-going 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

Limitations

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 and C. 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.