## Abstract

We studied heart rate variability in rats by power scaling spectral analysis (PSSA), autoregressive modeling (AR), and detrended fluctuation analysis (DFA), assessed stability by coefficient of variation between consecutive 6-h epochs, and then compared cross-correlation among techniques. These same parameters were checked from baseline conditions through acute and chronic disease states (streptozotocin-induced diabetes) followed by therapeutic intervention (insulin). Cross-correlation between methods over the entire time period was *r* = 0.94 (DFA and PSSA), *r* = 0.81 (DFA and AR), and *r* = 0.77 (AR and PSSA). Under baseline conditions the scaling parameter measured by DFA and PSSA and the high-frequency (HF) component measured by AR fluctuated around an average value, but these fluctuations were different for the three methods. After diabetes induction, a strong correlation was found between the HF power and the short-term scaling parameter. Despite their differences in methodology, DFA and PSSA assess changes in parasympathetic tone as detected by autoregressive modeling.

- heart rate variability
- autoregressive modeling
- detrended fluctuation analysis
- circasemiseptan rhythm
- power scaling spectral analysis

assessment of cardiac autonomic tone by analysis of heart rate variability (HRV) has been used to investigate such conditions as congestive heart failure, myocardial infarction, hypertension, neurological injury, diabetic neuropathy, prematurity of birth, pharmacodynamics, and mental stress and the effect of exercise. HRV has been assessed by using parametric (autoregressive modeling, AR) methods to determine the power spectral density (PSD) as well as by scaling techniques. The advantages of parametric methods for calculating the PSD are the following: smoother spectral components that can be distinguished independently of preselected frequency, easy postprocessing of the spectrum with an automatic calculation of low- and high-frequency (HF) power components and easy identification of the central frequency of each component, and an accurate estimation of PSD even on a small number of samples on which the signal is supposed to maintain stationarity (23). Nonstationarity in typical heart rate time series severely limits the range of frequencies that can be studied by conventional frequency-domain analytic methods. The frequency-domain analysis is used to measure the HF power to estimate parasympathetic nervous system (PNS) activity. A scaling technique in the frequency domain, power spectrum scaling analysis (PSSA), has also been used (4). PSSA estimates the scaling parameter as the slope of the linear regression of the log-log plot of the PSD. Only the linear region of the plot is selected for regression. PSD is calculated with a nonparametric method [fast Fourier transform (FFT)]. Because the method cannot identify the underlying structure of physiological fluctuations if there are trends due to external environmental influences, a scaling technique called detrended fluctuation analysis (DFA) was developed (25). DFA permits the detection of correlations embedded in the seemingly nonstationary time series, and this avoids the spurious detection of apparent long-term correlations that are an artifact of nonstationarities (17). In the present study, the scaling methods are used to investigate short-term correlations in the heartbeat expressing short-term memory between beats.

Although HRV has been studied in clinical and animal studies with a variety of techniques [6 studies (Refs. 2,7, 8, 13, 28,29) compared FFT to autoregressive power spectral analysis (AR) with a time period of 5–10 min; 2 studies (Refs.1, 6) compared FFT to AR with a 24-h epoch], there have been very few studies carrying out simultaneous comparisons. One study (15) compared spectral powers and scaling parameters with a 2-h epoch; two studies (17, 34) compared the scaling parameters with spectral powers with a short-term epoch; and two studies (22, 27) compared the two with a 24-h epoch. Yet no multiple HRV measurement techniques have ever been simultaneously investigated in instrumented rats under baseline conditions, disease states, and therapeutic interventions comparing spectral powers with scaling methods.

Cardiac autonomic neuropathy is a common complication in insulin-dependent diabetes mellitus (IDDM) and is clearly displayed by a decrease in HRV (3). Streptozotocin (STZ) is widely used to induce diabetes in animals as a model of IDDM (18). STZ impairs cardiac vagal tone in rats within 5 days of injection (21). There is a decrease in time-domain HRV (9) as well as frequency-domain HRV (10) in rats with chronic STZ-induced diabetes. This decrease in frequency-domain HRV has also been found in chronic STZ-induced diabetic Yucatan pigs (24). However, insulin treatment restores normal cardiovascular function during the early stages of diabetes in rats as measured by baroreceptor reflex control (5) and time-domain HRV (33).

Because longer-term physiological fluctuations may be due to endocrine systems and metabolic processes (17), and no long-term sequential HRV studies comparing PNS activity and scaling measures have ever been carried out using healthy rats made diabetic by STZ and then checking for the effect of insulin on these HRV parameters, we investigated the short- and long-term effect in healthy rats made diabetic by STZ, the effect of insulin on the scaling parameter as well as on the autonomic nervous system response, and the relations between them. The main question to be investigated was whether in disease states the short-term scaling parameter could serve as an index for changes in parasympathetic vagal tone as detected by classic spectral analysis.

## METHODS

### Animals

Three male Sabra rats weighing 300–350 g were housed separately in plastic cages and maintained on a 12:12-h light-dark cycle with food and water available ad libitum. The project adhered to the principles of laboratory animal care published by the National Institutes of Health (NIH Publication No. 85-23, Revised 1985).

### Compounds

A solution of 0.2 mg/ml streptozotocin (*N*-[methylnitrosocarbamoyl]-d-glucosamine; Sigma, St. Louis, MO) in physiological solution (0.9% sodium chloride injection USP; Teva Medical, Ashdod, Israel) was freshly prepared. Two international units of injectable insulin (Actrapid HM-biosynthetic human insulin; Novo Nordisk, Bagsvaerd, Denmark) were dissolved in 0.3 ml of physiological solution.

### Drug Treatment Protocol

After 1 wk of baseline recording, rats were injected with STZ solution. Three months after diabetic onset, the rats were treated with insulin injections.

The diabetic state was confirmed by monitoring blood glucose levels before and after diabetic onset and insulin treatment. Blood samples were obtained by tail prick, and blood glucose concentration was measured with a commercial monitor (Glucometer Elite XL; Bayer, Tarrytown, NY). Glucose concentrations were determined to be 100 mg/dl for baseline (“healthy”) condition, ≥300 mg/dl for the diabetic state, and <170 mg/dl after insulin treatment.

### Telemetry System

To minimize stress involved with data collection, a telemetric monitoring system was implanted in the peritoneal cavity under anesthesia (50 mg/kg ketamine and 10 mg/kg xylazine ip). Preventive antibiotic treatment was given subcutaneously 30 min before the operation. The telemetric device is a small (2 × 0.5") two-lead electrocardiogram (ECG) radio frequency transmitter (TCA-F40; DSI). The leads were tunneled subcutaneously to their positions at the right acrotrapezoidal muscle and the left gluteus muscle. Recording and experiments took place after a recovery period of 1 wk.

### Signal Acquisition System

Data were acquired with a telemetry system that contained the implantable radio frequency transmitter and a receiver (RCA-1020; DSI) located below the cage. Analog signals were transmitted to a computer (586 Pentium, 133 MHz) and digitized at a sampling rate of 600 Hz by an analog-digital converter (PCL 818 HG; ICPC). The R-R interval data were obtained online from the continuous ECG records by a threshold peak algorithm. The R peak was identified by crossing a dynamic threshold ranging between ±20% of the peak initially determined by the user. The following R peak was determined after a specified (35 ms) lag time to prevent a mistaken T peak identification. The threshold midrange was then set with the new detected value. R-R interval data expressed in milliseconds were recorded on continuous 1-h-length files.

### Data Analysis

#### Autoregressive preprocessing.

##### artifact elimination.

Epochs of 120,000- or 1,200,000-ms duration were chosen for calculation. Artifact elimination was performed by replacing an R-R interval sample that exceeded both the previous sample R-R interval and the epoch-average R-R interval by more than a prescribed value (50 and 80 ms, respectively) with the average R-R interval. Higher than 1% correction of a single epoch was determined as the exclusion criterion. In general, approximately <2% of the data were omitted (26). To reduce the effect of a slow nonperiodic variation of the parameters, detrending was performed by calculating a smoothed moving polynomial and subtracting it from the original signal (32).

##### AUTOREGRESSIVE POWER SPECTRUM ANALYSIS.

Autoregressive (AR) power spectrum analysis is based on time series parametric modeling. Parametric modeling for power spectrum estimation consists of an appropriate choice of a model, estimation of the parameters of the model, and then substitution of these estimated values into the theoretical power spectrum density (PSD) expressions. An AR time series model is chosen, the parameters (*a*) of the model are estimated, and PSD is calculated from these parameters. The model output, which is the estimated R-R interval at a time point*n*, is given by
where *p* is the order of the model and denotes the number of terms in the time series model and α_{k} is the *k*th coefficient of the model. The input driving process is assumed to be a white noise sequence of zero mean and variance of ς^{2}. The autocorrelation function, R(*k*), of the time series is first calculated. The Yule-Walker equations, which describe the relationship between the autocorrelation function and the AR parameters are then derived.
where ς^{2} is the variance of the white noise input of the AR model and *j* is an order index. To solve these equations for the AR parameters, the Levinson-Durbin algorithm is implemented. The Levinson-Durbin algorithm is initialized by
and the recursion for *k* = 2,3, …*p* is given by
where *i* is an order index. The Levinson-Durbin algorithm provides the coefficients of the autoregressive model of the order *p*. The power spectrum of an output of a autoregressive model with a white noise input of power spectrum of ς^{2}Δ*t* is related to the AR parameters by
where Δ*t* is the sampling interval of the original R-R signal calculated as the mean R-R interval during the epoch and*f* is the frequency. The model order *p* selection is generally chosen by several criteria such as the Akaike information criterion or Parzen's criterion autoregressive transfer function. The 16th-order AR model with a 301-point moving polynomial was found to be the most suitable and serves as a good compromise between a power spectrum that is too smooth and one with too many peaks (19, 30, 31). The logarithmic values of the distinct HF peak (range 1.35–2.65 Hz) served as markers (index) for changes in PNS activity.

#### Detrended fluctuation analysis.

DFA is a modified root-mean-square analysis of a random walk. In DFA, the numerical value of the scaling exponent α is indicative of the type and degree of correlation present in the heart rate data and can be thought of as an indicator of the “roughness” of the time series data (25).

The correlations in a time series between the steps*u _{i}
* and u

_{i}

_{+n}are defined by the autocorrelation function Equation 1where

*N*is the length of the series and

*n*is the distance between the steps. In systems in which there are long-range correlations, the nature of the correlations is of a power law Equation 2where γ is the autocorrelation exponent. By integration of a time series with long-range correlations, a self-affine process is formed. The scaling properties of the integrated series are connected to the long-range correlation properties of the original time series. If we consider then the profile of the time series, the accumulated sum

*Y*= Σ

_{n}*u*, its fluctuations F(

_{i}*n*) as a function of different scales

*n*are related to C(

*n*) and increase by a power law Equation 3where α is referred to as the scaling exponent or the self-similar parameter and is equivalent to the Hurst exponent in the range 0 < α < 1. However, the values α can take are not limited to that range.

The power spectrum of the original time series decreases by a power law
Equation 4The relation between α and β is
Equation 5The case of α = ½ still represents a fractal, but one that has no correlation, like the steps of a random walker. The fluctuations of a random walk (the standard deviation of its position) acts as a scaling law, *y* ∼ *t*½, with α = ½, so it is a fractal (11). When the data are correlated α > ½, and when the data are anticorrelated α < ½.

#### Calculations of short-time scale scaling exponent.

##### dfa.

The DFA is used to calculate the root mean square fluctuation of the integrated and detrended time series for different time scales and is actually a modified root-mean-square analysis of a random walk. The DFA (order *m*) steps are as follows: *1*) A profile is created
Equation 6
*2*) The profile is divided into nonoverlapping windows of size *l*, and then the local polynomial trend of order*m* in each window is calculated. This is achieved by the least-squares fit method. *3)* The variance is calculated around the local trend in each window, to eliminate trends of order*m − 1* in the original time series; the variances from all windows of size *l* are then averaged over number of windows, and the square root of the result is the fluctuation function F(*l*).

The fluctuation function is then plotted on a log-log scale, and the scaling exponent is calculated from the slope of the line. In some cases, there exists crossover from a short-scale scaling exponent to a long-scale scaling exponent, meaning that the function has a certain slope in the short-scale regime but a different slope in the long-scale regime. Two separate scaling exponents are then considered. The short-scale exponent, α_{1}, is the slope up to the crossover point, and the long-scale scaling exponent, α_{2}, is the slope from the crossover point (16). The scaling exponents are used as monofractal measures of the heartbeat.

DFA of order *m* measures the scaling exponent accurately in the range 0.5 < α < *m* + 0.5. When a scaling exponent α < ½ is suspected, extra integration is applied to the time series to increase the value of the scaling exponent by 1. The fluctuations are then F̃(*l*) ∼*l*
^{α+1}. DFA can then reliably calculate the new scaling exponent. We subtract 1 from the value of the new scaling exponent to obtain the original scaling exponent. In such a way, values of α < ½ can now be estimated correctly, including negative values.

We summed the data of each 20-min epoch and then applied the DFA of order 2 on the summed series. We then plottedF̃(*l*)/*l* in log-log scale and calculated the slope of the line up to a window of size 10, where there was the crossover point. This scaling exponent is denoted as α
.

##### PSSA.

The power spectrum S(f) is calculated with FFT of the time series where The square root of the power spectrum is then plotted into a log-log scale, and again the crossover point is determined. The slopes are then calculated and named for small frequencies and for high frequencies.

We plotted
in log-log scale and calculated β
from*f* = 0.1 beat^{−1} to *f* = 0.3 beat^{−1} to avoid the peak in the power spectrum associated with breathing. β
corresponds with α
as the high frequencies correspond with small windows. α
was then calculated with *Eq.5
*.

### Baseline and Disease State Measurements and Significance of Analysis

The analysis was done for each rat independently of the others. To compare between the methods, analysis of multiple observations was performed for the same rats with the actual values of α, β, or HF power for 20-min epochs. The values were cross-correlated between the different measures of the same data for the same rat (over the entire period of time), and therefore the sample size is considered not as the number of animals but rather the number of epochs.

The stability of the short-term scaling parameter and of the HF power was investigated during 24 consecutive 6-h epochs of HRV in healthy instrumented rats. Coefficients of variation between consecutive 6-h epochs were used to ascertain stability of measurement. Cross correlations among the three techniques were also evaluated.

Because of the similarity of the results obtained with different methods in different animals, an average value was used to describe the values at each steady condition of the rats (healthy, acute, chronic, insulin). The total epochs of the healthy state formed one average, to which each daily epoch in the disease states was compared. Hence, each state holds a distribution of values (corresponding to the 20-min epochs).

ANOVA tests were performed to distinguish the insulin treatment from the healthy state as well as from the acute and chronic diabetic states. Two-sample *t*-tests were performed for the baseline condition as well as all the diabetic days and the nadirs and peaks of oscillation.

## RESULTS

After the initial week of the baseline period, the rats were given STZ. After a 10-day waiting period, an acute effect of STZ-induced diabetic neuropathy was observed with all three techniques. The instrumented rats were kept in cages for another 3-mo period to evaluate the effects of chronic diabetes. At the end of this 3-mo period, HRV was investigated with the three techniques. After this measurement, insulin was given and the three HRV techniques were again applied. As parameters of measurement, we used only the HF power associated with respiration and the short-term scaling parameter.

Figure 1
*A* shows a typical example of an R-R interval file, and its PSD calculated with the fast Fourier transform is shown in Fig. 1
*B*. Figure2 shows the coefficient of variation of the three methods in healthy rats during a 5-day period. The least amount of variation was in DFA. Thus, although DFA had the lowest amount of noise, it was the least sensitive, with the HF power having >2.5 times higher levels. To investigate whether there were differences in the use of 2-min epochs vs. 20-min epochs, we compared AR plots of both epochs with a moving polynomial spline and found that the two patterns were highly similar (not shown). DFA and PSSA were performed over 20-min epochs only. In the healthy baseline stage, the HRV exhibited uncorrelated behavior in the short time scales (averaged scaling parameter close to 0.5), meaning that the heart beats were not influenced by previous beats (i.e., no short-term memory).

Figure 3, *A–C*, shows the correlation between measurements obtained by the three analysis methods over the entire time period (healthy, acute/chronic diabetic state induced by STZ, insulin). A comparison between the short-term scaling parameters as obtained by the two different scaling techniques, PSSA and DFA, produced the highest correlation (0.94, shared variance = 88%). High correlation was found between the short-term scaling parameter and the logarithm of the HF power. The linear correlation of the scaling parameter as measured by DFA and the logarithm of the HF power yielded (−0.85, shared variance = 72%), and between the scaling parameter as measured by PSSA and the logarithm of the HF power it yielded (−0.8, shared variance = 64%). The DFA and PSSA techniques are used to measure the same quantity (scaling parameter), so we expect strong correlation between the results of the two. However, the HF power is a different quantity that serves as a marker for PNS activity. Our results show that the effect of the disease on both quantities is correlated.

Figure 4 shows the percentage of change of the logarithm of the HF power as detected by AR and of the short-term scaling parameter by two different scaling techniques, through the diabetic progression and after insulin treatment. Because the measured parameters have both positive and negative values, a change of the values from positive to negative results in >100% decrease. From baseline to the STZ acute stage, the most sensitive method was AR with −170% change, followed by PSSA with −163%; the least sensitive method was DFA with −137% change (*P*< 0.0001 for all). From baseline to 13- to 23-day chronic STZ, the most sensitive method was AR with −67% change, followed by DFA with −41%; the least sensitive method was PSSA with −33% change (*P* < 0.0001 for all). From baseline to 3-mo chronic STZ the most sensitive method was the AR with −46% change (*P* < 0.0001). Although large differences were observed by the PSSA and DFA methods (−85% and −65% change, respectively), the significance of these results was less profound (*P*> 0.1 for both). From baseline to insulin treatment, the most sensitive method was PSSA with 47% increase, followed by DFA with 41% (*P* < 0.0001 for both); the least sensitive method was AR with −3% change (*P* < 0.05). After insulin treatment a significant increase of values was observed by all the three methods (*P* < 0.0001) from the far chronic diabetic state. A full recovery to control values was detected by both DFA and AR methods. More profound increase was detected by the PSSA method, where the values were found to be significantly higher than the control values (*P* < 0.05).

Figure 5 shows the evolution of the scaling parameter before and after STZ induction of diabetes. The division for the different states of the disease was based on this type of picture, which repeated itself for different rats. The α-values at the healthy state were not statistically different from 0.5, according to the single distribution *t*-test, therefore representing fractals that have no correlations. After the induction of STZ, HRV showed an anticorrelated behavior (scaling parameter < 0.5) in the short time scales. This shows that the following heartbeats act in an inverse manner to the previous beats. Significant change from the healthy condition (*P* < 0.05) was first observed by the DFA method, 3 days after STZ administration. Such significant change was observed by the other methods (AR and PSSA) 1 day later.

Figure 6 shows strong correlation between the daily averaged scaling parameters for different rats.

Figure 7 shows an interesting phenomenon of an apparent 96-h cycle (circasemiseptan rhythm) in the scaling parameter of the HRV starting ∼10 days after STZ induction. The same pattern was found with both scaling techniques and by the HF power.

## DISCUSSION

This is the first study to investigate the stability and cross-correlation of an indicator of PNS status and the short-term scaling parameter of HRV from baseline conditions through the evolution of a disease, followed by therapeutic intervention. This study measures and compares two different parameters, the HF power (1.5–2.65 Hz) and the short-term scaling parameter. HF power has been established in numerous investigations as an indirect marker related to parasympathetic tone and respiratory sinus arrhythmia (23). Using conventional AR spectral analysis, we found that in different stages of STZ-induced diabetes there were marked changes in PNS activity. This framework enabled us to assess the impact of changes in PNS activity on the short-term memory of HRV. Despite their differences in methodology and frequency, there was a high correlation (−0.85) between the parameters over the entire time period.

This suggests that a relationship exists between respiration and the short-term correlations, a possibility also raised by Toweill et al. (34). Because respiration is a periodic activity that modulates the heart rate, it affects the short-term correlations between heartbeats and introduces a characteristic time scale. Surprisingly, the increased breathing activity was connected with a stronger anticorrelation behavior of the short-term memory of the heartbeat.

The fractal measures of HRV complement the conventional measures because they can uncover “hidden” information in the time series data. DFA was developed in an attempt to distinguish between two different classes of fluctuations in physiological data sets, those that arise from changes in environmental conditions having little to do with the intrinsic dynamics of the system itself and those that arise from complex nonlinear dynamic interactions inherent to the system (25). Whereas in a study using techniques of HRV to predict imminent ventricular fibrillation DFA performed better than any other measure of HRV in differentiating between patients with ventricular fibrillation and controls (22), in another study predicting survival in patients with congestive heart failure, DFA was of borderline predictive significance when multivariate models were used (15).

We found that the scaling parameter obtained by DFA had the least amount of noise in baseline (having the lowest coefficient of variation), and this scaling parameter was accordingly used to define the status of the rat. The scaling parameters obtained by PSSA and DFA, β^{PSSA} and α^{DFA}, are linearly related (14). Indeed, our findings show a high correlation between α
and β
.

We found that, after STZ administration, the scaling parameter decreased for 10 days and then reached a minimal value. After that acute chronic phase, a milder chronic state was obtained. Throughout the chronic period, although the HF power and the scaling parameter values fluctuated, they were significantly lower than the baseline condition. Previously, such a significant decrease in the HF power was recorded after 12–18 wk in STZ-induced diabetic rats (10) and also at age 8–9 mo in the WBN/kob spontaneous diabetic rat model (12).

Schaan et al. (33) showed that insulin infusion increased time domain parameters of HRV in STZ diabetic rats during the early stages of diabetes. We found that insulin administration (via the intramuscular route) increased HF power in 3-mo chronic STZ diabetic rats and restored it to baseline (healthy) levels. More profound increase was detected by both the short-term scaling methods, where in the case of PSSA values even significantly exceeded the control (healthy) levels. In all three methods, however, remarkable changes were detected before and after insulin administration. The effect of insulin administration in the healthy state was not included in the original design of the study.

Whether this study can be extrapolated to human studies remains to be seen. The STZ model is not exactly analogous to human diabetes. A common feature of STZ-induced diabetes is glucose insensitivity (20), which does not occur in the human disease. Finally, the intriguing phenomenon of a 96-h cycle (circasemiseptan rhythm) in the HF power and scaling parameter (Fig. 7) may represent an infradian rhythm.

Our findings show that the short-term scaling parameter is related to the HF power and is proportional to its logarithm. Therefore, we may quantify the changes in the complex behavior of HRV in relation to the disease state and as an index for autonomic nervous system functioning. The changes in parasympathetic vagal tone can be assessed by changes in the short-term scaling parameter.

## Acknowledgments

A. Hoffman is affiliated with the David R. Bloom Center of Pharmacy. This work is part of the PhD dissertation of I. Perlstein.

## Footnotes

Address for reprint requests and other correspondence: A. Hoffman, Dept. of Pharmaceutics, School of Pharmacy, Hebrew Univ. of Jerusalem, PO Box 12065, Jerusalem 91120, Israel (E-mail:ahoffman{at}cc.huji.ac.il).

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. Section 1734 solely to indicate this fact.May 30, 2002;10.1152/ajpheart.00519.2001

- Copyright © 2002 the American Physiological Society