## Abstract

The power spectrum of human heart rate (HR) measured over 24 h exhibits “power-law” 1/*f *
^{α}-type spectral behavior with α ≈ 1. We investigate possible nonstationarity in time of the exponent α using maximum likelihood estimation, which allows relatively short data segments to be used. Examination of 24-h HR records from ambulatory normal and congestive heart failure (CHF) subjects indicates that the power-law structure of HR is nonstationary. In addition, α varies with time scale and is different for normal (α ≈ 1) and CHF (α ≈ 1.5) subjects. Simulations suggest that a possible mechanism underlying the observed power-law spectrum may be a switching between values of α near zero (white noise) and near two (Brownian motion). This mechanism generates power-law forms quantitatively similar to CHF subjects when the switching occurs very rapidly and similar to normal subjects when the switching is less rapid.

- power law
- control
- self-similarity
- model
- surrogate data

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 1
*A* shows a 2-h segment of a 24-h R-R interval time series. Figure 1
*B* 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 *
^{α}.

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

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

When one divides a long time series into segments, one is implicitly setting the lowest frequency that can be considered in each segment. If we wish to include frequency *f* in the analysis, our segments must have length at least 1/*f* s. For instance, for frequency scales down to 0.0002 Hz, data segments must be >5,000 s.

The interval *T*
_{s} between samples sets the highest frequency that can be considered (1/2*T*
_{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/*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^{7}; 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.2
*A* and 3
*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. 2
*B* and 3
*B,* where the results are formatted in terms of *T*
_{s} = 1/2*T*
_{s}.

### Assessing Statistical Significance

Figure 2 shows the results of applying the TRA technique to a synthetic signal whose power spectrum (estimated from the FFT) is a perfect power law with α = 1. It can be seen (Fig. 2
*A*) that the estimatesα̂ vary somewhat from segment to segment. Because the synthetic signal is constructed to be stationary, we know that this is sampling variability and results from the fact that each segment is a sample from the perfect power-law signal.

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_{1}, d_{2}, … , d_{10}). Then, we compare each of the surrogate estimates to each other and obtain 10 ∗ 9/2 = 45 measures D_{s} = (ds_{1}, ds_{2}, … , ds_{45}). 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

To illustrate the stationarity-in-time and TRA analyses, we consider three synthetic time series: *1*) a “perfect” stationary 1/*f* time series where α is constant in time and the same at all frequencies; *2*) a time series where α changes in time but is the same at all frequencies at each time; *3*) a time series where α is different for different frequency bands, but is constant in time.

#### α 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.

For the stationarity-in-time analysis, the time series was sampled every 10 s (that is, the 8,192 point time series was used directly) and overlapping segments of length *L* = 256 points were used. For the TRA analysis, the time series was sampled at *T*
_{s}= 10 s, 20 s, … 150 s and divided into segments of length 256 points.

Figure 2 and Table 1 show the results. Sampling fluctuations are visible in Fig. 2
*A,* but these are not qualitatively different than the surrogate data (shown as *x* in Fig. 2
*A*). The histograms showing the spread of α̂ are essentially identical; the interquartile range (IQR) of the distribution of α̂ from the 1/*f* time series (IQR_{d}) and the IQR of the distribution of α̂ from surrogates (IQR_{s}) are roughly the same. The K-S computations similarly show that the distribution of α̂ for the 1/*f* time series is not significantly different than that of the surrogates, which is underlined by a significance level *P* < 0.40. Thus there is no indication of nonstationarity in time here, as would be expected for this type of synthetic signal.

In Fig. 2
*B,* the estimate α̂ is plotted for each sampling time *T*
_{s} for each segment. For fast sampling times (e.g., 10 s) there are many segments in the whole time series and consequently many estimates α̂ (plotted as *x*). For slow sampling times (e.g., 100 s), there are fewer segments and thus fewer estimates for α̂. Note that at a sampling time of 100 s, the time scales comprehended by one estimation range up to 7.1 h. There is no evidence for a time-scale dependence of α, as expected for this type of signal.

#### α 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. 3
*A,* 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).

Figure 3
*B* depicts all the estimates α̂ versus sampling time *T*
_{s}. Here α̂ is approximately constant over a wide range of time scales, although there is no single α in this data set. This suggests that a nonperfect power-law process may appear to have a α̂ that is independent of frequency. But note that the range of α̂ in Fig. 3
*B* is much wider than in Fig. 2
*B,* reflecting the nonstationarity of the signal in Fig.3. This is also underlined by the results from the statistical analysis given in Table 1, *P* < 0.0001.

#### α 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^{16}= 65,536 was synthesized. The power spectrum (Fig.4
*A*) 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. 4
*B*) closely mimicks power spectra seen in studies of HR.

The stationarity-in-time analysis (Fig.5
*A*, *C,* and *D;* Table1) gives no indication of nonstationarity, as is to be expected for this stationary time series. The TRA analysis (Fig. 5
*B*) shows a clear drift in α̂ with time scale.

The power spectrum in Fig. 4
*A* gives a clearer indication of MLE of the structure of α as a function of frequency f than does TRA (Fig. 5
*B*). However, this is misleading. Only those data in Fig. 4
*B* were used in the TRA analysis, and realistic data of any length will always have the sampling fluctuations seen in Fig.4
*B.* The precise straight-line form in Fig. 4
*A* does not incorporate the sampling variability that must be present in any experimentally collected data set. We encourage the reader to perform the following experiment. Without reference to Fig. 4
*A,* locate the frequency at which the bend in Fig. 4
*B* occurs. Compare your estimate to the frequency clearly indicated in 4*A.* We find that an error by a factor of two is typical.

Quantitative estimate of α from the power spectrum requires that linear regression be performed. Figure 6shows the result of such a regression on the power spectrum in Fig.4
*B.* Linear regression was performed over a window from*f*
_{lower} to *f*
_{upper}, and α̂was estimated as the slope of the line as a function of*f*
_{upper}. To facilitate comparison with TRA, we have formatted the results as a function of *T*
_{s} = 1/*f*
_{upper}. Figure 6 shows regressions done with frequency windows of seven (same as used in TRA), six, and five octaves. The shorter windows do not make clearer the sharp transition from α = 2 to α = 0.8.

## 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 8show 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.

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).

The stationarity-in-time analysis is presented for all of the subjects in Fig. 10 and Table2. 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.*

## DISCUSSION

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. 11
*A*) and 32 data points long were concatenated (Fig. 11
*B*). For a model segment length of 256 data points, the estimates α̂ are able to resolve the time series into 1/*f*
^{0} and 1/*f*
^{2}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.

Figure 12
*A*, 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.

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.

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

Our statistical analysis of 24-h HR records indicates that the 1/*f*
^{α} structure in HR fluctuates in time and that α varies with time scale. Analysis of 24-h records as a whole shows α is ∼1, but analysis of short segments (∼1 h) shows that α fluctuates significantly. When looked at another way, concatenating many 1-h segments with widely varying α produces a 24-h signal with α ≈1.

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.

## Acknowledgments

We thank C. K. Peng and A. Goldberger for providing us with the CHF and normal HR data used in this paper.

## Footnotes

Address for reprint requests: D. T. Kaplan, Dept. of Mathematics and Computer Science, Macalester College, 1600 Grand Ave., St. Paul, MN 55105.

This study was partially supported by a “Kurt-Gödel-Fellowship” of the Austrian Ministry for Science and Research, by a fellowship from theÖsterreichische Forschungsgesellschaft, and by grants from the Natural Sciences and Engineering Research Council of Canada, the Fonds de la recherche en santé du Québec, and the Fonds pour la formation de chercheurs et l’aide à la recherche.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. §1734 solely to indicate this fact.

- Copyright © 1999 the American Physiological Society