INTRODUCTION

Top Abstract Introduction Methods Discussion References

IT HAS BEEN NOTED by numerous investigators (10, 12, 16, 20) that the variance of human heart rate (HR) or R-R intervals tends to increase when examined on longer and longer time scales. When examined using spectral analysis, this increase has the form of a "power law" where the squared amplitude of the Fourier transform at frequency f is proportional to 1/f. Generically, this form is called "one-over-f noise" and indeed it is found empirically that in heart rate   1, as first reported by Kobayashi and Musha (10).

Figure 1A shows a 2-h segment of a 24-h R-R interval time series. Figure 1B depicts the power spectrum P( f ) for the entire 24-h record. One can divide the power spectrum into two main parts. At frequencies above ~0.02 Hz (or, equivalently, time scales <50 s) the HR does not have 1/f structure and in fact shows the two well-known broad peaks: one at the frequency of respiration and one at a time scale of 10 s. At frequencies below ~0.02 Hz the power spectrum is well approximated by a straight line of slope   1. Because the power spectrum is being plotted on log-log axes, this straight line corresponds to a power-law relationship P( f ) ~ 1/f .

View larger version (45K): [in this window] [in a new window]   Fig. 1.   A: 2-h segment of a 24-h R-R interval time series. B: power spectrum of the entire 24-h record, estimated using a fast Fourier transform (FFT.)

Although 1/f noise is observed in a wide diversity of physical, technological, and biological systems, e.g., river flows (4), traffic densities (9), and loudness fluctuations in speech and music (22), there is currently no general explanation for its ubiquity and no generally satisfactory framework for its modeling in specific instances. In cardiac physiology, where there are numerous interacting control systems and adapting mechanisms with a wide range of time scales (3), it is of interest to know why a single exponent  should describe the relationship between fluctuations on time scales from 1 min to several hours.

We are interested here in two questions: Is the exponent  constant, independent of frequency (for f < 0.2 Hz)? Is the exponent  constant over time?

	METHODS

Top Abstract Introduction Methods Discussion References

The fundamental method we employ for assessing possible changes in  with time involves dividing a 24-h R-R interval time series into short segments, computing an estimate in each segment and checking whether is statistically significantly different from one segment to another. However, we confront the trade-off between bias and variance: to track rapid changes in  we need to use short segments, but to get reliable estimates we want segments to be long. Although this trade-off is unavoidable, we optimize our estimates by using an efficient maximum likelihood estimator (MLE) of . This MLE has been found to provide estimates that are superior to linear regression on the logarithmic power spectrum or to Hurst-exponent type estimators (13, 1).

Time-Rescaling Analysis

The interval T_s between samples sets the highest frequency that can be considered (1/2T_s). In studying the 1/f structure of R-R intervals, we would like to be able explicitly to set both the upper and lower frequency bounds. When using a fast Fourier transform (FFT) estimator this is straightforward: explicitly set the frequency window over which linear regression is performed and ignore other frequencies. The MLE technique does not directly permit this, but we can still impose frequency bounds via the segment length and by setting the sampling interval T_s. The method that we have developed, which we term time rescaling analysis (TRA), consists of the following steps: 1) interpolate the raw R-R interval sequence into evenly spaced samples with T_s = 0.5 s; 2) pick a lower and upper frequency of interest (f_lower and f_upper); 3) anti-alias filter and resample the time series with a sampling interval of T_s = 1/2 1/f_upper; 4) divide the resampled time series into segments of length (L) = 1/f_lower (in the work reported here, we make our segments L = 256 points, giving f_lower = f_upper/128 =  f_upper/2⁷; thus we cover 7 octaves of the frequency scale). 5) For each segment, use MLE to compute the estimate .

At this point, we have a time series of : one estimate for each segment. This time series can be examined directly to detect changes in  over time, as in Figs. 2A and 3A.

View larger version (29K): [in this window] [in a new window]   Fig. 2.   Analysis of the simulated perfect 1/f process. A: estimated exponent vs. time for Ts = 10 s for simulated data (solid line) showing mean ± 2 SD (dotted lines) of all surrogate sets and 1 set of surrogates (x marks). B: Time rescaling analysis (TRA) showing exponents vs. sampling time Ts with mean ± 2 SD at each Ts. C: histogram for shown in A. D: histogram for from surrogates shown in A.

View larger version (32K): [in this window] [in a new window]   Fig. 3.   Analysis of the simulated time-varying 1/f  process. A: estimated exponent vs. time for Ts = 10 s for simulated data (solid line) showing mean ± 2 SD (dotted lines) of all surrogate sets and 1 set of surrogates (x marks). Dashed line indicates exponent  used for generation of simulated data as a function of time. B: TRA showing exponents versus sampling time Ts with mean ± 2 SD at each Ts. C: histogram for shown in A. D: histogram for from surrogates shown in A.

To investigate how  depends on the frequency f_lower, we repeat steps 1-5 for several different values of f_lower, as seen in Figs. 2B and 3B, where the results are formatted in terms of T_s = 1/2T_s.

Assessing Statistical Significance

In considering whether nonstationarity in  is responsible for segment-to-segment variation in , we need to take into account sampling variability. One way to do this is to use synthetic data of the sort generated in Fig. 2 that has a perfect power-law structure. This approach was used in previous work studying the sampling distribution of the MLE estimator (13). However, it has not yet been established that R-R interval data is indeed perfect power-law data. We have therefore used a more general technique for generating synthetic data: the method of surrogate data (7, 8, 19).

Surrogate data are obtained in the following way: 1) take the Fourier transform of the original time series; 2) randomize the phases at each frequency to be uniformly distributed in [0, 2]; 3) take the inverse Fourier transform.

Because the amplitude of the Fourier transform remains unchanged in step 2, the power spectra of the surrogate data and the original data are identical. Different realizations of surrogate data, all with the same power spectrum, can be generated by using different random phases in step 2.

Surrogate data exhibit the following important characteristics. 1) By construction, surrogate data are equivalent to linear filtering of stationary Gaussian white noise. 2) The process that generates surrogate data is stationary, i.e., the coefficients of the filter process remain constant as do the properties of the Gaussian white noise.

Thus surrogate data stem from a linear and stationary process. We note that although surrogate data have the same linear correlations as the original data, any nonlinear correlations in the original data do not appear in the surrogate data. Insofar as nonlinear correlations exist in R-R interval data [a controversial issue that is the subject of current research (6, 7, 14, 17)], the surrogate data will be qualitatively different from the R-R interval data. However, although nonlinearities might possibly be important in the generation of 1/f noise, estimates of  based on spectral analysis or the MLE method used here will be insensitive to the nonlinear correlations in the data. Therefore we claim that surrogate data provide an appropriate statistical technique for estimating the sampling distribution of the MLE estimator. However, without a specific model of how some nonlinearity underlies 1/f noise, we do not know how to construct an empirical test of the validity of this claim.

Using surrogate data, we can estimate the sampling distribution for for a stationary process with the same power spectrum as the 24-h R-R interval data. We use three approaches to investigating whether the distribution from surrogate data differs from that of the original R-R interval data; we term this the stationarity-in-time analysis.

The first approach is to graphically display the histograms of the segment-by-segment estimates for the R-R interval data and a surrogate data set. If the R-R interval data have a nonstationary , the histogram should be broader for the R-R interval data than for the surrogate data. In addition, we plot out a time series of versus time for both the R-R-interval data and one realization of surrogate data.

The second approach is to quantify the width of the distribution of using the interquartile range and compare the width of the estimates for the R-R interval to the estimates for 10 different realizations of surrogate data using a two-sided t-test.

A third approach is to characterize general differences in the distribution using the Kolmogorov-Smirnov (K-S) test (15, 21). The K-S measure d is defined as the maximum value of the absolute difference between two cumulative distribution functions. For each R-R interval time series we obtain a set of estimates . Ten realizations of surrogate data are used, producing 10 sets of from stationary processes. First, we compare the data estimates to each of the 10 surrogate estimates using the K-S statistic and obtain 10 measures D = (d₁, d₂, ... , d₁₀). Then, we compare each of the surrogate estimates to each other and obtain 10 * 9/2 = 45 measures D_s = (ds₁, ds₂, ... , ds₄₅). Finally, we perform Student's t-test to determine whether the two samples D and D_s could have the same mean.

We note that surrogate data were developed originally (19) to detect nonlinearity in possibly chaotic data. This use involves employing a statistic that is sensitive to nonlinear structure. Here, we use a statistic that is based on the autocorrelation function and hence models the signal as linear. Because of this, differences between the surrogates and the test data do not point directly to nonlinearity of the test data. Instead, they reflect other violations of the null hypothesis associated with surrogate data, in particular the stationarity of surrogate data.

Illustration of the Techniques Using Synthetic Data

Stationary in time, independent of frequency. A perfect 1/f time series of length 65,536 sample points was generated using a spectral synthesis method, and then a subsample of length 8,192 was used as our test data as in Ref. 13. A length of 8,192 is roughly the same as the 24-h R-R interval time series when sampled every 10 s.

Changes in time, independent of frequency. The ability of MLE to track nonstationarity was tested by estimating the exponent  of test data simulating time-varying behavior of exponent . The data set was generated by concatenating time series, 256 data points long, with time-varying  in the range 0.5-1.8. As seen in Fig. 3A, the estimates (solid line) track closely the  used in the simulation (dashed line). The surrogates show no such pattern. The histogram of is much broader for the nonstationary data than for the surrogate data (which has more or less the same width as in the stationary case illustrated in Fig. 2).

Stationary in time, but changes with frequency. The ability of TRA to discern variations in  at different frequencies is illustrated by application to a time series with frequency-varying exponent . A time series of length 2¹⁶ = 65,536 was synthesized. The power spectrum (Fig. 4A) shows that  = 2 for f < 0.0017 and  = 0.8 for f  0.0017. For the analysis, a segment of length 8,192 was extracted. The power spectrum of this segment (Fig. 4B) closely mimicks power spectra seen in studies of HR.

View larger version (52K): [in this window] [in a new window]   Fig. 4.   A: Power spectrum of a perfect frequency-varying 1/f  process with  = 0.8 for time scales smaller than 10 min, i.e. frequencies higher than f = 1.7 · 103,  = 2 for larger time scales, and of length n = 216 (64 K). B: power spectrum of a segment of length n = 213 (8 K) of the signal in A.

View larger version (29K): [in this window] [in a new window]   Fig. 5.   Analysis of the simulated frequency-varying 1/f  process shown in Fig. 4. A: estimated exponent vs. time for Ts = 10 s for simulated data (solid line) showing mean ± 2 SD (dotted lines) of all surrogate sets and 1 set of surrogates (x marks). B: TRA showing exponents vs. sampling time Ts with mean ± 2 SD at each Ts. C: histogram for shown in A. D: histogram for from surrogates shown in A.

View larger version (16K): [in this window] [in a new window]   Fig. 6.   Exponents from linear regression on a moving window of the power spectrum in Fig. 4B vs. frequency fupper of the window (formatted in terms of sampling time Ts). Three different lengths of the regression window are shown, covering a frequency range of 7 (solid line), 6 (dashed line), and 5 octaves (dash-dotted line) data points.

	ANALYSIS OF HEART RATE

We applied the stationarity-in-time and TRA analysis to 20 recordings of R-R intervals over 24 h from which abnormal beats had been deleted. These time series were provided to us by C. K. Peng and A. L. Goldberger (see Ref. 12). There were 9 healthy subjects and 11 congestive heart failure (CHF) patients.

Figures 7 and 8 show the stationarity-in-time and TRA analysis for a representative normal subject and a representative CHF patient. In both cases, nonstationarity in  is suggested by the broad histograms of for the R-R interval data compared with the surrogate data. The CHF patient shows clear sustained deviations in from the mean value that last for several hours.

View larger version (31K): [in this window] [in a new window]   Fig. 7.   Results of the analysis of case no. 16,786, a normal subject. A: estimated exponent vs. time for Ts = 10 s for R-R interval data (solid line) showing mean ± 2 SD (dotted lines) of all surrogate sets and 1 set of surrogates (x marks). B: TRA showing exponents vs. sampling time Ts with mean ± 2 SD at each Ts. C: histogram for shown in A. D: histogram for from surrogates shown in A.

View larger version (30K): [in this window] [in a new window]   Fig. 8.   Results of the analysis of case no. 9,706, a congestive heart failure (CHF) subject. A: estimated exponent vs. time for Ts = 10 s for R-R interval data (solid line) showing mean ± 2 SD (dotted lines) of all surrogate sets and 1 set of surrogates (x marks). B: TRA showing exponents vs. sampling time Ts with mean ± 2 SD at each Ts. C: histogram for shown in A. D: histogram for from surrogates shown in A.

The TRA analysis shows that for the CHF subject increases from 1 at T_s near 20 s to 1.5 at T_s near 100 s. This pattern is sustained in almost all of the records. Figure 9 shows the TRA analysis for all 20 records. The CHF patients typically show a higher than the normal patients at large time scales (corresponding to low frequencies).

View larger version (22K): [in this window] [in a new window]   Fig. 9.   Estimated exponents vs. sampling time Ts for 24-h R-R interval data from 11 CHF patients (dotted lines) and 9 normal subjects (solid lines). Mean for each subject at each Ts is plotted.

The stationarity-in-time analysis is presented for all of the subjects in Fig. 10 and Table 2. In almost all subjects, the IQR of the distribution of is roughly twice as large in the R-R interval data compared with the width expected purely due to sampling variability (as indicated by surrogate data). The K-S computations similarly show that the distribution of for the R-R-interval data is different than that of the surrogates and point to this difference being statistically significant at a level P < 0.01 in all except CHF subject 9,778. 

View larger version (14K): [in this window] [in a new window]   Fig. 10.   Interquartile range (IQR) of vs. patient index no. for sampling time Ts = 10 s. A: IQR of from normals (solid line) and mean IQR ± 2 SD of all corresponding surrogate sets (dashed and dotted lines). B: IQR of from CHF patients (solid line) and mean IQR ± 2 SD of all corresponding surrogate sets (dashed and dotted lines).

View this table: [in this window] [in a new window]   Table 2.   Results of t-test using K-S statistic describing differences between distributions of data estimates  and surrogate estimates for sampling time Ts = 10 s

	DISCUSSION

Top Abstract Introduction Methods Discussion References

R-R interval time series show in almost all cases "mixed" process behavior, i.e., exponent  exhibits time and frequency dependence. Similar results were also found by Meesmann et al. (11), who showed that in the HR of young healthy subjects, periods of 1/f noise had intermittent subperiods of white noise during the night. Ichimaru and Katayama (5) suggest that the exponent  may be different for each sleep stage. Di Rienzo et al. (2) found that measured via power-spectral regression varied with frequency f.

One possibility suggested by Fig. 3 is that nonstationarity and the observed 1/f pattern in HR are related. In that figure,   1 at large T_s, although  is not constant. For instance, at T_s = 140 s, estimates are being made using segments that span 10 h, during which time  is varying substantially, yet at these long segment lengths does not vary much from segment to segment.

Although there is no general theoretical explanation for systems with  = 1, there are simple theoretical models for  = 0 (white noise) and  = 2 (Brownian motion, the cumulative sum of white noise). Is it possible that a system that switches between white noise and Brownian motion might produce  = 1?

From our stationarity-in-time analysis of R-R interval data, there is no reason to suspect that the periods of white noise or Brownian motion might last for as long as one analysis segment (2,560 s). If this were the case, we would expect to see values of near 0 or 2, but we do not. Instead, we hypothesize that the switching might occur much faster, over say one-quarter or one-half of the analysis segment (320 s or 640 s).

In terms of control systems, Brownian motion corresponds to control being turned off, whereas white noise variability indicates that control is on: the system may be knocked off its set point by an outside disturbance, but quickly returns to it, only to be knocked off again. But why should HR control be turned off?

HR is one of several effector variables involved in regulating the cardiovascular system. By "effector variable" we mean a feature of the system that is used to bring a regulated quantity to the appropriate level or range of levels. Other effector variables are peripheral resistence, blood volume, and vascular compliance. Examples of regulated quantities are blood pressure and blood gas levels that are directly sensed or cardiac output measured indirectly via blood gas levels and blood pressure. Note that HR is generally considered not to be a regulated quantity; it is a means to an end rather than an end in itself. Although there are baroreceptors and chemoreceptors to measure blood pressure and blood gas levels, there is no analogous receptor for HR. This suggests that HR might be allowed to drift (Brownian motion). However, drifts in HR might be punctuated by periods of resetting when control systems attempt to use HR as an effector variable and therefore HR is locked in to the appropriate value for the quantity that is being regulated.

In addition, control systems acting on very long time scales (e.g., the renal blood-volume control system) might indirectly affect HR via other feedback loops. This could cause slow, long-term drifts in HR. All told, HR might be seen as episodes of Brownian motion interrupted by episodes of white noise for short time scales, all on a background of slow trends.

To model this hypothesized structure for HR, we concatenated short segments of Gaussian white noise and Brownian motion. For each segment, we made a random choice between white noise and Brownian motion. As a background, we added to each segment a linear trend with a random direction for each segment. These trends were matched at their endpoints. Note that none of these three processes individually produce 1/f noise: white noise has  = 0, Brownian motion has  = 2, and the piecewise linear trends have  = 2 for low frequencies and  = 0 for high frequencies.

Figure 11 shows the histograms of the estimated exponents (sampling time T_s = 10 s) for data produced by our model process where model segments 256 (Fig. 11A) and 32 data points long were concatenated (Fig. 11B). For a model segment length of 256 data points, the estimates are able to resolve the time series into 1/f⁰ and 1/f² segments, expressed by two main peaks at  = 0 and  = 2 in the histogram. For model segments of length 32,  is not able to resolve the time-varying character of exponent . The histograms of from surrogate data are substantially narrower than for the model data, providing the same kind of evidence for nonstationarity seen in the R-R interval data.

View larger version (26K): [in this window] [in a new window]   Fig. 11.   Histograms from the model nonstationary process (see DISCUSSION). A1: results from data generated by concatenation of model segments 256 data points long. A2: results from its surrogates. B1: results from data generated by concatenation of model segments 32 data points long. B2: results from its surrogates.

Figure 12A, which depicts the TRA analysis of vs. T_s shows a pattern similar to that seen in the R-R interval data. We point out the similarity of the short model segment (32 points) TRA to that found in normal subjects, and the long model segment (256 points) TRA to that found in CHF patients.

View larger version (11K): [in this window] [in a new window]   Fig. 12.   TRA analysis on the model nonstationary process. vs. sampling time Ts for 4 nonstationary test data sets generated by concatenating white and Brownian motion segments with different segment lengths, 256 (solid line), 128 (dash-dotted line), 64 (dashed line), and 32 (dotted line) data points. A: TRA results from time series with local random trends added. B: TRA results from time series without local random trends.

Figure 12B shows the effect of deleting the local linear trends part of the model. For the short model segment ("healthy") data, there is a loss of verisimilitude.

Visual comparision (Fig. 13) of a randomly selected segment of R-R interval data from a normal subject and a synthetic signal from the nonstationary segments-with-trend model of  suggests that this simple model produces realistic time series with realistic power spectra.

View larger version (43K): [in this window] [in a new window]   Fig. 13.   A1: a 2-h segment of an R-R interval time series sampled at Ts = 10 s. A2: 2-h segment from the nonstationary model of . B1: power spectrum of a 24-h R-R interval time series showing 1/f-like spectral behavior. B2: power spectrum of 24-h simulated time series showing 1/f-like spectral behavior.

This simple model of the 1/f structure of HR suggests that the power-law structure of HR may be nonstationary over fairly short intervals, similar to the 5-min segments commonly used in short-term HR variability analysis (18). It also points to a possible difference between cardiovascular control in normal and CHF subjects.

In addition, the stationarity-in-time analysis raises the question of whether it is possible to obtain more information about the 1/f behavior of HR data by analyzing longer time series. The answer might be no. The range of exponents encountered indicates that the observed 1/f pattern might be a property of the interaction of the estimators with nonstationarity in the exponent , as we have seen in the simulations.

Perspectives

We speculatively carry this process to an extreme: using computer simulations we show that by concatenating very short segments with  = 0 or  = 2 one produces 1-h segments where  is close to 1 but varies significantly on either side of it and 24-h patterns with  quite close to 1. We cannot confirm this rapid-switching model by direct measurement of  on very short segments: the statistical methods have too high a variance on the short segments and we would need to use even shorter segments to find out when the hypothetical switching occurs.

It may be possible to assess the rapid-switching model by measurements other than . For instance, the amplitudes of the well-known 10-s and respiratory HR spectral peaks may vary as cardiac control switches from an  = 0 to an  = 2 pattern. At this point in the research, we see the value of the model as a way of supporting and organizing the hypothesis that the 24-h 1/f structure may be an epiphenomenon and is the statistical result of a conceptually much different process.