Nonlinear measures of heart rate time series: influence of posture and controlled breathing

Autonomic Neuroscience: Basic and Clinical 83 (2000) 148–158 www.elsevier.com / locate / autneu

Nonlinear measures of heart rate time series: influence of posture and controlled breathing a a b, R.K.A. Radhakrishna , D. Narayana Dutt , Vikram Kumar Yeragani * b

a Department of ECE, Indian Institute of Science, Bangalore, India Department of Psychiatry, Wayne State University School of Medicine, Detroit, MI, USA

Received 9 February 2000; accepted 27 April 2000

Abstract In this study, we investigated measures of nonlinear dynamics and chaos theory in regards to heart rate variability in 27 normal control subjects in supine and standing postures, and 14 subjects in spontaneous and controlled breathing conditions. We examined minimum embedding dimension (MED), largest Lyapunov exponent (LLE) and measures of nonlinearity (NL) of heart rate time series. MED quantifies the system’s complexity, LLE predictability and NL, a measure of deviation from linear processes. There was a significant decrease in complexity (P,0.00001), a decrease in predictability (P,0.00001) and an increase in nonlinearity (P50.00001) during the change from supine to standing posture. Decrease in MED, and increases in NL score and LLE in standing posture appear to be partly due to an increase in sympathetic activity of the autonomous nervous system in standing posture. An improvement in predictability during controlled breathing appears to be due to the introduction of a periodic component.  2000 Published by Elsevier Science B.V. Keywords: Nonlinear; Chaos; Spectral; Heart rate variability; Posture; Breathing

1. Introduction Heart rate variability (HRV) has been extensively used as a noninvasive tool to study cardiac autonomic function in health as well as in disease (Malik and Camm, 1990; Malliani et al., 1991; Yeragani, 1995). A decrease in HRV has been found to be an important predictor of cardiac events and sudden cardiac death in patients with cardiac illness as well as normal controls (Molgaard et al., 1991; Bigger et al., 1992). Spectral analysis has shown that the short-term HR time series reveal a peak around 0.1 Hz, which is related to the baroreceptor mechanisms, a peak between 0.15 and 0.5 Hz, due to respiratory sinus arrhythmia (RSA), and also a peak below 0.04 Hz, related to peripheral vascular and thermoregulatory mechanisms (Akselrod et al., 1981; Pomeranz et al., 1985; Lindqvist et al., 1990). A relative increase in low frequency (LF: 0.04–0.15 Hz) power and a decrease in high frequency (HF: 0.15–0.5 Hz) power from supine to standing posture are attributed to a predominance of sympathetic activity *Corresponding author, Wayne State University School of Medicine, Flat No. 16, K.C.N. Mansion, Bangalore-560001, India. Tel.: 191-80228-7715; fax: 191-80-341-5694. E-mail address: [email protected] (V.K. Yeragani).

and vagal withdrawal in standing posture (Pagani et al., 1986; Malliani et al., 1991; Yeragani et al., 1993a, 1993b). Goldberger and West (1987) stressed the importance of nonlinear techniques to study HRV. The output from most of the biological systems such as the cardiovascular system appears to be due to the result of internal dynamics as well as input into the system from external factors. Several recent studies have used different nonlinear measures including fractal dimension (FD), correlation dimension (CD), approximate entropy (APEN), measures of nonlinearity (NL), measures derived from symbolic dynamics and 1 /f scaling to characterize HR time series (Kobayashi and Musha, 1982; Katz, 1988; Glenny et al., 1991; Pincus et al., 1991; Yeragani et al., 1993b, 1997, 1998, 2000; Guzzetti et al., 1996; Lombardi et al., 1996; Voss et al., 1996; Ho et al., 1997; Braun et al., 1998; Storella et al., 1998; Kagiyama et al., 1999; Silipo et al., 1999). Nonlinear dynamics can be applied to study systems in which output is not proportional to input. Fractal structure has been described in several biological systems and some of the studies suggest that the condition of health may be associated with a high degree of complexity of these systems. FD measures the space-filling propensity and complexity of the time series (Katz, 1988). APEN quantifies the regularity in time series data and measures the

1566-0702 / 00 / $ – see front matter  2000 Published by Elsevier Science B.V. PII: S1566-0702( 00 )00173-9

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

logarithmic likelihood that runs of patterns that are close remain close on the next incremental comparisons, and the higher the value, the more random the time series (Pincus et al., 1991). We have previously shown that there is a high degree of correlation between FD (as measured by the Katz algorithm) and APEN and the HF power of HR time series, thus reflecting cardiac vagal activity (Yeragani et al., 1993b, 1997, 1998). CD is a measure of system complexity and has been used to study HR time series in health and disease (Kagiyama et al., 1999). Symbolic dynamics is based on coarse graining technique and the time series of HR or HP are converted into symbol sequences using the mean value of the time series and a parameter of deviation from the mean (Voss et al., 1996). Word count, a measure obtained from symbolic dynamics may be effectively used to understand the nonlinear complexity of HR time series (Yeragani et al., 2000). As the HR time series are fractallike, the spectral densities of long-term time-series conform to 1 /f distribution (Kobayashi and Musha, 1982) where the log-power is inversely proportional to logfrequency of the power spectrum. Quantification of nonlinearity (NL) measures basically provide information as to the deviation of the time series from linear surrogates and have not been fully explored in regards to HR time series (Di Garbo et al., 1998). Some of these techniques have underscored the importance of the addition of nonlinear measures to study cardiac autonomic function in health and disease (Pincus et al., 1991; Yeragani et al., 1993b, 1997, 1998, 2000; Guzzetti et al., 1996; Lombardi et al., 1996; Voss et al., 1996; Ho et al., 1997; Braun et al., 1998; Kagiyama et al., 1999; Makikallio et al., 1999). Analysis of time series using methods of nonlinear dynamics consists of estimating dynamical invariants of the reconstructed attractor, such as dimensions, Lyapunov exponents (LE) and degree of nonlinearity (NL). Considering a simple example of Lorenz equations which is expressed as a set of three coupled nonlinear differential equations looks innocuous enough, but the dynamics the equations display have been a source of wonder and intense mathematical study dx(t) ]] 5 s ( y(t) 2 x(t)), dt dy(t) ]] 5 2 x(t)z(t) 1 rx(t) 2 y(t), dt dz(t) ]] 5 x(t)y(t) 2 bz(t). d(t) In the above equations, ‘x’,‘y’ and ‘z’ are the variables, and ‘s ’, ‘r’ and ‘b’ are the constants. Thus one can say that this system can be understood in three dimensional phase space and the minimum embedding dimension (MED) required is 3. The CD which is the actual number of degrees of freedom and gives an estimation of the system’s complexity is always less than this MED and in this case it is 2.14, and also the CD will be fractal. These

149

equations are very sensitive to initial conditions. If one starts with two very close initial conditions after some time it is seen that their trajectories start diverging sufficiently, meaning if the exact initial condition is not known, one will not be able to predict the system’s behavior correctly beyond a certain time. This ability of predicting or sensitivity to initial conditions is quantified in terms of LE. If a system can be put in a phase space whose MED is ‘m’, then there will be ‘m’ number of Lyapunov exponents, and all of these are not required to characterize the system’s sensitivity to initial conditions. The largest of them, largest Lyapunov exponent (LLE) is sufficient in most cases. The LLE for Lorenz equations is 1.5, the reciprocal of this is proportional to time which gives an estimation of predictable time for the system. The constants used in Lorenz equations play a vital role in deciding the behavior of the system dynamics. Their interrelation or interaction decides the amount of linearity or nonlinearity in the system dynamics, which is quantified by the degree of nonlinearity (NL). There are many algorithms available in the literature to estimate these invariants, but most of them require a long-term time series to obtain reliable estimates. In most real world situations, acquiring long stationary time series may be a tough task for various reasons. In our present study, we have investigated methods which work well on short-term time series, and result in reliable estimates of dynamical invariants. In this study, we have specifically investigated the effect of postural change from supine to standing, and the effect of controlled breathing on minimum embedding dimension (MED), largest Lyapunov exponent (LLE) and measures of NL of HR time series to understand if these measures provide additional information compared to the frequently used frequency domain measures such as spectral powers in various bands, and how these new measures may relate to the function of the autonomic nervous system.

2. Subjects and methods

2.1. Subjects Twenty-seven normal controls (16 males, 11 females; 30.267.5 years (mean6S.D.)) participated in this study. The studies were part of ongoing projects related to the investigation of cardiac autonomic function in normal controls. These studies were approved by the Institutional Review Boards at the Wayne State University School of Medicine, Detroit, MI and the Wright State University School of Medicine, Dayton, OH, USA. All subjects were healthy and informed consent was obtained prior to their participation in these studies. Electrocardiogram (ECG) was recorded by a HewlettPackard 78173A ECG Monitor in lead II configuration. The signal was recorded onto a PC using a 12-bit A / D

150

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

board, at a sampling rate of 500 Hz. This was done on a PC with custom-developed software (HEM Data Corp., Southfield, MI, USA). Supine data were obtained after the subjects rested for 10 min and the standing data were obtained 2 min after the subjects stood up. We used 256 s of supine as well as standing data for the analyses in this study. We used a peak detection algorithm to identify the R–R intervals (in ms) from the ECG.

2.2. Data processing We used HR time series free of ventricular premature beats and noise. The HR (beats / min: 60,000 / R–R interval in ms) time series were sampled at 4 Hz using the method described by Berger et al. (1986) to obtain instantaneous HR. This step-wise continuous instantaneous HR signal maintains an amplitude equal to the reciprocal of the current R–R interval. This sampling rate would allow an accurate estimation of the power spectrum up to 1 Hz. All 27 subjects had both supine and standing data during spontaneous breathing and 14 of them had data during controlled breathing (12 breaths per min) as well. All the analyses were done using custom developed software on IBM compatible personal computers. The reconstruction of HR time series and the calculation of the MED, LLE and NL, were all computed automatically using a PC with custom designed software according to the following methods.

2.2.1. Time delay embedding and attractor reconstruction In a typical time series analysis scenario, there is a time series from a system whose dynamics and parameters are all unknown. We do not even know how many components or parameters are responsible for this output time series. The best way to start is to estimate the MED, i.e. to determine the number of variables responsible in generating this time series. Again considering the typical example of Lorenz equations, if we have time behavior of one of the variables x(t), which has embedded within it the effect of the other two variables. The task is to reconstruct the attractor in the required phase space with the available observed time series of a single variable. An attractor is a set of points in phase space visited by a single trajectory (path traced by a signal vector) as the time evolves. An attractor is referred to as a ‘strange attractor’ if it has a non-integer dimension and associated with nonperiodic motions. A strange attractor is also characteristic of chaotic systems. The reconstructed attractor should preserve all properties of the original attractor. The theorem proposed by Taken et al. ensures this reconstruction and the procedure is better known as time-delay embedding. If x(n) for n 5 1, 2, . . . , N is the time series, then time-delay vectors in phase-space are formed as:

X(i) 5 [x(i 1 t ), x(i 1 2t ) . . . x(i 1 (m 2 1)t ] i 5 1, 2, . . . , N 2 (m 2 1)t where t is the time delay and m is the embedding dimension. There are many methods to choose the timedelay. We have used the autocorrelation method in our analysis. Fig. 1a and b shows both the original and reconstructed attractors of a Lorenz system in a chaotic regime. Fig. 1c and d correspond to attractors of HR time series in supine and standing posture in three-dimensional phase space for pictorial representation (MED is assumed to be 3 here). Fig. 2a and b show HR time series of a subject in supine and standing postures.

2.2.2. Estimation of minimum embedding dimension ( MED) Proper reconstruction of an attractor is guaranteed if the dimension of phase-space is sufficient to unfold the attractor. It is shown that an embedding dimension of m . 2d 1 1 can achieve this, where d is the dimension of the attractor (Takens, 1980). In most cases of observed time series analysis we neither have knowledge of d or m. There are many different algorithms used to estimate these quantities (Grassberger and Procaccia, 1983; Broomhead and King, 1986; Mees et al., 1987; Theiler, 1987; Kennel et al., 1992), but many of them have the disadvantage of either being too subjective requiring a large number of data points or being computationally very intensive. The method proposed by Liangyue Cao (1997) overcomes these difficulties and is suitable for short-term time series. Additionally, this method gives more reliable estimates of MED even when the dimension is sufficiently large. After time delay embedding, we compute a quantity iXi (m 1 1) 2 Xn(i, m) (m 1 1)i a(i, m) 5 ]]]]]]]] iXi (m) 2 Xn(i, m) (m)i 5 1, 2, . . . , N 2 mt

i (1)

where i . . . i is some measurement of Euclidean distance and is given in this paper by maximum norm defined as: iXk (m) 2 Xl (m)i 5 max ux k 1jt 2 xl 1jtu. 0#j #m21

(2)

Xi (m 1 1) is the ith reconstructed vector in dimension. n(i, m) is an integer in the range 1 # n(i, m) # N 2 dt such that Xn(i, m) (m) is the nearest neighbor of Xi (m) in the dimensional reconstructed phase space. The n(i, m) in the numerator of Eq. (1) is the same as that of the denominator. If Xn(i, m) (m) equals Xi (m), then we take the second nearest neighbor. A quantity E(m) which is the mean value of all a(i, m)s is computed as

O

1 N 2mt E(m) 5 ]] a(i, m). N 2 mt i 51 This averaging removes subjectivity involved with fixing threshold and tolerances of the false nearest neigh-

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

151

Fig. 1. (a, b) The original and reconstructed attractors of Lorenz system; (c, d) the attractors of HR time series in supine and standing postures.

Fig. 2. (a, b) The time series of the heart rate in beats per min in supine and standing postures.

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

152

bor (FNN) method (Kennel et al., 1992). E(m) is dependent only on the dimension m and lag t. To investigate its variation from m to m 1 1, we define E1(m) 5 E(m 1 1) /E(m). It is found that E1(m) stops changing when it is greater than some value m 0 if the series comes from an attractor, then m 0 1 1 is the minimum embedding dimension, which accommodates the attractor efficiently. This technique is applied to variable ‘x’ of the Lorenz system and in Fig. 3a, it can be clearly seen after m 5 3 we do not observe any change in E1(m) indicating that this time series came from a system whose attractor can be embedded in phase space of dimension 3. Deterministic signals are those which can be described mathematically by a formula whereas stochastic signals require statistical description. Another quantity is determined which is useful in distinguishing deterministic signals from stochastic signals and it is given by 1 E*(m) 5 ]] N 2 mt

O ux

N 2mt

i 1mt

2 x n(i, m)1mtu

i 51

and its variation from m to m 1 1 as

E2(m) 5 E*(m 1 1) /E*m where n(i, m) has the same meaning defined earlier (Eq. (1)). For random time series E1(m) will never attain a saturation value as m is increased, but because of limited data samples and practical computations it may be difficult to ascertain whether E1(m) is slowly changing or has stopped changing. In such situation E2(m) will be very useful, since for random data future values are independent of past values, E2(m) will be equal to 1 for any m, whereas for deterministic signals there exists some ms such that E2(m) ± 1. We computed both E1(m) and E2(m). We applied this method on the time series of some of the standard maps and found their MED tallying with the literature. Fig. 3a shows the calculation of MED of the Lorenz attractor. Fig. 3b is the plot of E1(m) from which MED is determined as 14 for supine and 11 for standing conditions.

2.2.3. Largest Lyapunov exponent ( LLE) Lyapunov exponents are another invariant, which could be used to characterize the dynamical system. It quantifies sensitivity of the system to initial conditions. An mdimensional dynamical system has m Lyapunov exponents.

Fig. 3. (a, b, c) The calculation of the MED for Lorenz system and MED and LLE for the HR time series in supine and standing postures.

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

The presence of a positive Lyapunov exponent indicates chaos. It also quantifies the amount of instability or predictability of the system. A fully deterministic system will have a zero Lyapunov exponent since it is fully predictable, whereas a random system will have a large positive exponent indicating no predictability. In most applications it is sufficient to compute only the largest Lyapunov exponent (LLE) instead of all Lyapunov exponents. The LLEs quantify predictability of the system, i.e. how far in time the prediction can be done which can assess the model used for analysis. The reciprocal of LLE is proportional to prediction time, i.e. time over which predictions are reasonable. A fully deterministic system has zero LLE suggesting that it is predictable over infinite time (e.g. y(t) 5 A.sin(wt)), whereas a random or stochastic system will have large positive LLE indicating it is predictable only over a very short time. The LLE corresponding to Lorenz is 1.5. There are many algorithms available to estimate LLE and the Lyapunov spectrum (Sano and Sawada, 1985; Wolf et al., 1985; Sato et al., 1987; Zeng et al., 1991). Most of them are unreliable when operated on small data sets. In our present work we use a method proposed by Rosenstien et al. (1993) which is robust against data length. After reconstructing the attractor, this algorithm looks for nearest neighbors of each phase point Xi on the trajectory. Distance between two neighboring points at instant n 5 0 is defined by d i (0) 5 min Xj iXj 2 Xi i where i . . . i is the Euclidean norm. This algorithm imposes constraint that nearest neighbors are temporally separated at least by the mean period of the time series. The LLE is then estimated as the mean rate of separation of nearest neighbor’s, i.e. we can write d j (i) ¯ Cj e

l 1 (iDt )

(3)

where Cj is the initial separation. Taking the logarithm on both sides of (3) we obtain ln (d j (i)) 5 ln Cj 1 l1 (iDt).

(4)

Eq. (4) represents a set of approximately parallel lines (for j 5 1, 2 . . . , M), where the slope is roughly proportional to the LLE. In practice, the Lyapunov exponent is easily and accurately estimated using a least-squares fit to the average line defined by 1 y(n) 5 ]kln d i (n)l Dt where kl denotes the average over all values of i. This last averaging step is the main feature that allows an accurate evaluation of l even when we have a short and noisy data set. In this algorithm, the points track the divergence of all pairs of nearby points and then the average determines the divergence rate. Applying this algorithm on two HRV time

153

series shown in Fig. 2a and b, we found their values to be 0.1103 for supine and 0.1708 for standing conditions (Fig. 3c).

2.2.4. Tests for nonlinearity The erratic fluctuations that are observed in an experimental time series owe their dynamical variation to a mix of various influences: chaos, nonchaotic but still nonlinear determinism, linear correlation, and noise, both in the dynamics and in the measuring set-up. This emphasizes the need for estimating nonlinear structure in the time series. In our present work we investigate nonlinear structure present with HRV time series using two methods. We check whether nonlinear time correlations are present among the time series values. The two protocols used are based on number of local extremas and their broken length. Since extremas are related to singularities in complex time plane, this is related to system dynamics. A larger value for ‘S’ (NL score) means the system is more nonlinear or it is deviating from a linear process which shares many properties of the system under consideration like mean, standard deviation, autocorrelation and power spectrum, by a greater extent. With our HR time series it is noted that the ‘S’ value is notably higher in standing than in supine posture indicating the time series is more nonlinear in standing posture than supine posture. This is indicative of change in process behavior from one posture to the other. Both methods are based on the analysis of the extrema (local maxima or minima) as proposed by Di Garbo et al. (1998). 2.2.5. Nonlinearity test based on extrema of a time series Snet and Spsc are significance values indicative of nonlinearity when two different protocols are used. Snet stands for significance value when a protocol based on number of extremas is used whereas Spsc is the significance value when a protocol based on the length of a line connecting these extremas is used. It has been shown that the dynamical behavior of the real time solution of an ordinary differential equation (ODE) is strongly connected to its analytic properties in the complex time plane, and in particular to the distribution of the singularities nearest to the real axis (Ramani et al., 1989). The second consideration arises from a general property of a stochastic process, which states that given a mean square differentiable stochastic process, the expected number of its extrema for unit time is contained in the joint density function of x(t), x9(t) and x0(t) (Soong, 1973). These theoretical and numerical results suggest that the sequence of extrema of a time series contain dynamical information of the process generating them. Both methods statistically discriminate measures evaluated based on extremas for original and surrogate data sets. Two types of surrogates are considered in our analysis,

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

154

Gaussian scaled (GS) and phase randomized (PR) surrogates. The following is the procedure used in generating Gaussian scaled surrogates. 1. 2. 3. 4. 5.

Histogram transformation. Fast Fourier transformation (FFT). Phase randomization. Inverse FFT. Inverse histogram transformation.

In the case of PR surrogates, steps 1 and 5 are eliminated. We will give a brief description of the two techniques used in quantifying nonlinear structure in the time series.

2.2.6. Pattern of singularities in the complex time plane ( PSC) algorithm The steps involved in quantifying nonlinear correlations with the PSC method are as follows. 1. Determine the couples hs t j , t j for j 5 1, 2, . . . nj corresponding to local maxima and time, at which it occurred. 2. Determine the length of the broken line joining these extremas ]]]]]]]] n21

L5

œ O hs j 51

s t j 11 2 s t jd 2 1st j 11 2 t jd 2 j.

3. n number of surrogates are generated and L for each surrogate is computed. 4. Determine mean Ls and standard deviation ss of these quantities. 5. Determine the measure of significance as proposed by Theiler et al. (1992)

uL 2 Lsu Spsc 5 ]]. ss

2.2.7. Number of extrema for unit time ( NET) The following is the protocol of the NET method. 1. The number of extrema N0 for unit time, T 0 of the given time series is determined and used as discriminating statistics. 2. n numbers of surrogate data sets are generated and number of extremas for each surrogate set Ni (i 5 1, . . . , n) are computed. 3. The average number of extrema for unit time Ns and their standard deviation ss are determined and they are statistically discriminated by computing the significance: uN0 2 Nsu Snet 5 ]]]. ss

2.3. Statistical analysis We used paired t-tests to compare linear and nonlinear measures of HRV between supine and standing postures and two-way ANOVA for repeated measures to compare spontaneous and controlled breathing conditions during supine and standing postures. Pearson’s product moment correlations were used to examine the relationship between linear and nonlinear measures of interest. All tests were two-tailed and a probability value of 0.05 was accepted as significant.

3. Results Table 1 shows the results of supine and standing nonlinear and linear measures of HR. As expected, there

Table 1 Supine and standing nonlinear and linear measures of HR (mean6S.D.) (N527)a

Mean HR (bpm) MED LLE Snet Spsc Total power (0–0.5 Hz) VLF power (0–0.04 Hz) LF power (0.04–0.15 Hz) Relative LF power HF power (0.15–0.5 Hz) Relative HF power LF / HF ratio

Supine

Standing

t

P

61.667.8 11.8961.81 0.14260.025 7.2362.73 6.77610.02 2.5760.86

79.5612.6 9.9362.69 0.18860.039 12.1464.02 10.0265.14 3.6360.74

9.44 5.78 7.10 5.12 2.73 6.64

,0.00001 ,0.00001 ,0.00001 0.00001 0.011 ,0.00001

1.4361.0

2.7160.72

5.99

,0.00001

1.5360.86

2.8660.92

7.48

,0.00001

37.6613.2 1.0961.05

48.2614.2 1.0660.94

3.0 0.18

0.007 0.86

26.2613.2 2.061.6

9.366.5 7.966.9

7.0 4.4

,0.00001 0.0002

a Spectral power units are in Ln of beats per min squared. MED, minimum embedding dimension; LLE, largest Lyapunov exponent; Snet and Spsc , nonlinearity measures (NL).

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

155

Fig. 4. The significant change of HR, MED, LLE and NL scores from supine to standing posture.

was a significant increase in HR, total power, VLF and LF powers, and LF / HF ratio, and a decrease in relative HF power in standing posture. There was a highly significant decrease in MED while there were significant increases in LLE, Snet and Spsc in standing posture. Snet showed a more significant difference compared to Spsc scores (Fig. 4). Table 2 shows the results of nonlinear measures during spontaneous and controlled breathing. A two-way ANOVA showed a significant increase of MED during controlled breathing in supine posture (F 59.82; df51,13; P50.008) and a significant decrease of LLE during controlled breathing, especially during supine posture (F 512.5; df5 1,13; P50.004). Measures of NL did not significantly change during controlled breathing (Fig. 5). There was also a significant increase of HF power during controlled breathing (F 541.2; df51,13; P,0.00001). The correlation analysis has not shown any significant relationship between age (in this age group) and any of the nonlinear measures. There was a negative correlation

between LLE and MED in supine as well as standing postures (r5 20.62 and 20.94, respectively). LLE correlated significantly with supine VLF power (r50.41; P5 0.03), standing total power (r50.51; P50.005), and standing VLF power (r50.42; P50.03) and LF power (r50.52; P50.005). There were no significant correlations between LF / HF ratio and LLE, MED or NL measures in supine or standing postures.

4. Discussion

4.1. Spectral analysis We found that there was a significant increase in HR, relative LF power, an increase in LF / HF ratios, and a decrease in relative HF power, which is in agreement with previous reports (Pagani et al., 1986; Malik and Camm,

Table 2 Supine and standing nonlinear measures during spontaneous and controlled breathing (CB) (mean6S.D.) (N514)a

Mean HR (bpm) MED LLE Snet Spsc HF power (0.15–0.5 Hz)

Supine

Supine CB

Standing

Standing CB

60.166.2 12.6361.28 0.12960.019 7.2263.27 7.7864.25 0.7461.04

60.065.1 14.3661.01 0.10160.02 7.1362.32 7.3563.40 1.6460.97

73.269.2 11.061.92 0.15960.013 11.1862.62 10.0962.93 0.7160.89

76.269.0 11.2161.19 0.15660.017 12.9162.70 11.8064.32 1.3260.81

a Spectral power is in Ln of beats per min squared. MED, minimum embedding dimension; LLE, largest Lyapunov exponent; Snet and Spsc , nonlinearity measures (NL).

156

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

Fig. 5. The significant change of MED and LLE from supine to supine controlled breathing condition.

1990; Yeragani et al., 1993a, 1993b). This is in line with decreased cardiac vagal activity and a relative increase in sympathetic function in standing posture. Though it is believed that LF / HF ratios reflect sympathovagal ratios, this issue has been disputed by some investigators (Cacioppo et al., 1994).

4.2. Nonlinear measures The degree of freedom as indicated by MED, decreased from supine to standing posture, whereas the time series become more nonlinear as shown by increased NL scores. The decrease in MED is significantly correlated with an increase in LLE in both postures. Thus the increased LLE and decreased MED in standing posture appear to be partly due to the dominance of the sympathetic system and vagal withdrawal. The MED computed were in the range of 8–15, indicating a high dimensional chaos. This is similar to the values of ,20 used in the study of Kagiyama et al. (1999). The chaotic nature of the signal was verified while computing MED. Among the two types of surrogates the PR surrogates were expected to perform better since measurements were based on number of local maximas and their broken length. Other recent studies have also

reported that nonlinearity could be detected with HR measures and could be quantified (Storella et al., 1998; Silipo et al., 1999). LLE describes the rate of exponential divergence of trajectories and the sensitive dependence on the initial condition. The decrease in predictability (increase in LLE) from supine to standing posture can be due to an increased nonlinear interaction of variables. This again may be related to increased sympathetic activity and changes in peripheral vascular mechanisms in standing posture, as there is a positive correlation between standing VLF and LF powers, and LLE. However, it should be noted that these correlations are about 0.5 or lower, suggesting that measures such as LLE yield additional information from the spectral measures. The correlation analysis has shown that age is not significantly related to any of these nonlinear measures in the age group studied. It is also important to note that LLE did not significantly correlate with HF power either in supine or standing postures and thus may not be significantly related to respiratory sinus arrhythmia. Other measures such as FD using the method of Katz (1988) and APEN using the technique of Pincus et al. (1991), correlate significantly with HF power and thus reflect RSA which is related to cardiac vagal function (Yeragani et al., 1993a, 1993b, 1998). Thus LLE, NL and MED measures appear to give additional information about HR time series. We have also found that Word count obtained from long-term HR time series using the technique of Voss and co-workers also does not correlate significantly with other measures yielding unique information (Yeragani et al., 2000). However, we cannot reliably obtain Word count on time series of 256 s duration. Kagiyama et al. (1999) have recently reported the effect of tilt from supine posture on nonlinear measures of CD and LLE in normal controls and patients with essential hypertension. They have found that CD was significantly lower in patients compared to controls but LLE values were not. The LLE values were not significantly different during tilt compared to supine posture in their study. On the contrary, we have shown highly significant increases in LLE in standing posture. This discrepancy is difficult to explain because spectral measures in their study also revealed an increase in LF / HF ratios and a decrease in relative HF power similar to our results. It is difficult to believe that standing posture vs. tilt may have been responsible for this difference. Kagiyama et al. (1999) used R–R intervals and also 8 Hz sampling rates for the chaos analysis and 1 Hz for spectral analysis, whereas we used a sampling rate of 4 Hz for all the analyses. They have used 500 points for chaos analysis and 256 data points for spectral analysis while we used 1024 data points for either analysis. The mean age of controls in our study was 30 whereas in the study of Kagiyama et al. it was 40 years. It should also be noted that the methodology to calculate LLE is different between the two studies and

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

future studies using different techniques to calculate these nonlinear indices may help address these issues.

4.3. Controlled breathing During spontaneous breathing, the subject does not pay any attention to his or her breathing and during controlled breathing, the subject breathed at a particular frequency for the duration. However, the subjects were asked to breathe normally. All subjects were conscious during the studies. Controlled breathing in supine posture resulted in an increase in MED compared to normal breathing inferring involvement of extra system parameters. At 12 per min breathing rate, there is an increase in HF power, which is likely due to an increase in cardiac vagal activity (Weise and Heydenreich, 1989) and this might have contributed toward an increase in MED. Also the decrease in MED in standing posture may have been partly due to decreased vagal activity. However, the correlations between MED and HF power were not significant and this measure should be evaluated further in future studies in regards to its significance. The improvement in predictability (decreased LLE) with controlled breathing in supine posture can be attributed to a more rhythmic nature of the signal. However, Kanters et al. (1997) did not find any significant differences in CD, a measure of complexity in normal controls between spontaneous vs. 12 / min breathing.

157

nature, suggesting a combination of stochastic and nonlinear models may be more appropriate to analyze them. Some nonlinear measures correlated modestly with linear measures such as spectral LF power, whereas, some invariants are independent measures. The tools used to estimate dynamical invariants are robust with data length. The protocol used for NL measures confirms its reliability for low dimensional time series. The estimated invariants of the attractor are significant and infer information about underneath dynamics, which may be applicable to the study of cardiac autonomic function in health and disease. These nonlinear measures appear to be valuable additions to the more commonly used linear measures in different studies investigating cardiac autonomic function in various disorders as some studies have already suggested decreased magnitude of nonlinear measures in cardiac disease and hypertension (Guzzetti et al., 1996; Lombardi et al., 1996; Ho et al., 1997; Kagiyama et al., 1999; Makikallio et al., 1999). Some recent studies have applied waveletbased time series analysis, which can characterize nonstationary behavior and elucidate such phase interactions (Ivanov et al., 1996). There is also evidence to suggest that there is a rhythmic structure of the chaotic component of the HRV over a 24-h period (Curione et al., 1998). Thus future studies should address the additional value of these measures to discriminate the states of health and illness and should investigate the amount of NL during awake as well as sleep conditions.

5. Conclusions The present study of HRV time series using nonlinear techniques has shown highly significant differences between supine and standing postures, and thus could give additional insight into underlying dynamics of HRV and in the investigations of cardiac autonomic function.

6. Perspectives The nonlinear dynamics play an important role in analyzing physiological time series. The complexity of heart rate dynamics encourages a researcher to use nonlinear models rather than linear or stochastic models. These models can throw light on heart dynamics and result in better understanding. The nonlinear measures used in this paper are based on local extremas and do not depend on mean and standard deviation of time series, i.e. it is providing extra information on system dynamics. The time taken to compute them is also not very large. But, all nonlinear dynamical measures do not enjoy this luxury and most of the measures are computationally intensive. Also, the role of noise in the time series becomes very crucial. Large noise content can lead to diverged inference. We observed that HR time series are not of low dimensional chaotic systems, but of high dimensional

References Akselrod, S., Gordon, D., Ubel, F.A., Shannon, D.C., Barger, A.C., Cohen, R.J., 1981. Power spectrum analysis of heart rate fluctuation: a quantitative probe of beat-to-beat cardiovascular control. Science 213, 220–222. Berger, R.D., Akselrod, S., Gordon, D., Cohen, R., 1986. An efficient algorithm for spectral analysis of heart rate variability. IEEE. Trans. Biomed. Eng. 33, 900–904. Bigger, J.T., Fleiss, J.L., Steinman, R., Rolnitzky, L.M., Kleiger, R.E., Rottman, J.N., 1992. Frequency domain measures of heart period variability and mortality after myocardial infarction. Circulation 85, 164–171. Braun, C., Kowallik, P., Freking, A., Hadeler, D., Kniffki, K.D., Meesmann, M., 1998. Demonstration of nonlinear components in heart rate variability of healthy persons. Am. J. Physiol. 275, H1577–H1584. Broomhead, D.S., King, G.P., 1986. Extracting qualitative dynamics from experimental data. Physica D 20, 217–236. Cacioppo, J.T., Berntson, G.G., Ginkley, P.F., Quigley, K.S., Uchino, B.N., Fieldstone, A., 1994. Autonomic cardiac control. II. Noninvasive indices and basal response as revealed by autonomic blockades. Psychophysiology 31, 586–598. Curione, M., Bernardini, F., Cedrone, L., Proietti, E., Danese, C., Pellegrino, A.M., De Ross, R., Di Siena, G., Vacca, K., Cammarota, C., Cugini, P., 1998. The chaotic component of human heart rate variability shows a circadian periodicity as documented by the correlation dimension of the time-qualified sinusal R–R intervals. Clin. Ter. 149, 409–412. Di Garbo, A., Balocchi, R., Chillemi, S., 1998. Nonlinearity tests using the extrema of a time series. Int. J. Bifurcation Chaos 8, 1831–1838.

158

R.K. A. Radhakrishna et al. / Autonomic Neuroscience: Basic and Clinical 83 (2000) 148 – 158

Glenny, R.W., Robertson, H.T., Yamashiro, S., Bassingthwaighte, J.B., 1991. Application of fractal analysis to physiology. J. Appl. Physiol. 70, 2351–2367. Goldberger, A.L., West, B.J., 1987. Fractals in physiology and medicine. Yale J. Biol. Med. 60, 421–435. Grassberger, P., Procaccia, I., 1983. Measuring the strangeness of strange attractors. Physica D 9, 189–208. Guzzetti, S., Signorini, M.G., Cogliati, C., Mezzetti, S., Porta, A., Cerutti, S., Malliani, A., 1996. Non-linear dynamics and chaotic indices in heart rate variability of normal subjects and heart-transplanted patients. Cardiovasc. Res. 31, 441–446. Ho, K.K., Moody, G.B., Peng, C.K., Mietus, J.E., Larson, M.G., Levy, D., Goldberger, A.L., 1997. Predicting survival in heart failure case and control subjects by use of fully automated methods for deriving nonlinear and conventional indices of heart rate dynamics. Circulation 96, 842–848. Ivanov, P.C., Rosenblum, M.G., Peng, C.K., Meitus, J., Havlin, S., Stanley, H.E., Goldberg, A.L., 1996. Scaling behavior of heart beat intervals obtained by wavelet-based time series analysis. Nature 383, 323–327. Kagiyama, S., Tsukashima, A., Abe, I., Fujishima, S., Ohmori, S., Onaka, U., Ohya, Y., Fujii, K., Tsushihashi, T., Fujishima, M., 1999. Chaos and spectral analyses of heart rate variability during head-up tilting in essential hypertension. J. Auton. Nerv. Syst. 76, 153–158. Kanters, J.K., Hojgaard, M.V., Agner, E., Holstein-Rathlou, N.H., 1997. Influence of forced respiration on nonlinear dynamics in heart rate variability. Am. J. Physiol. 272, R1149–1154. Katz, M.J., 1988. Fractals and the analysis of waveforms. Comp. Biol. Med. 18, 145–156. Kennel, M.B., Brown, R., Abarbanel, H.D.I., 1992. Determining minimum embedding dimension using a geometrical construction. Phys. Rev. A 45, 3403–3411. Kobayashi, M., Musha, T., 1982. 1 /f fluctuation of heart beat period. IEEE. Trans. Biomed. Eng. 29, 456–457. Cao, L., 1997. Practical method for determining the minimum embedding dimension of a scalar time series. Physica D 110, 43–50. Lindqvist, A., Jalonen, J., Parviainen, P., Antilla, K., Laitinen, L.A., 1990. Effect of posture on thermally stimulated cardiovascular oscillations. Cardiovasc. Res. 24, 373–380. Lombardi, F., Sandrone, G., Mortara, A., Torzillo, D., Rovere, M.T.L., Signorini, M.G., Cerutti, S., Malliani, A., 1996. Linear and nonlinear dynamics of heart rate variability after acute myocardial infarction with normal and reduced left ventricular ejection fraction. Am. J. Cardiol. 77, 1283–1288. Makikallio, T.H., Koistinen, J., Jordaens, L., Tulppo, M.P., Wood, N., Golasarsky, B., Peng, C.K., Goldberger, A.L., Huikuri, H.V., 1999. Heart rate dynamics before spontaneous onset of ventricular fibrillation in patients with healed myocardial infarcts. Am. J. Cardiol. 83, 880–884. Malik, M., Camm, A.J., 1990. Heart rate variability. Clin. Cardiol. 13, 570–576. Malliani, A., Pagani, M., Lombardi, F., Cerutti, S., 1991. Cardiovascular neural regulation explored in the frequency domain. Circulation 84, 482–492. Mees, A.I., Rapp, P.E., Jennings, L.S., 1987. Singular value decomposition and embedding dimension. Phys. Rev. A 37, 340–346. Molgaard, H., Sorensen, K.E., Bjerregard, P., 1991. Attenuated 24-h heart rate variability in apparently healthy subjects, subsequently suffering sudden cardiac death. Clin. Auton. Res. 1, 233–237. Pagani, M., Lombardi, F., Guzzetti, S., Rimoldi, O., Furlan, R., Pizzinelli, P., Sandrone, G., Malfatto, G., Dell’Orto, S., Piccaluga, E., Turiel, M., Basellli, G., Cerutti, S., Malliani, A., 1986. Power spectral analysis of heart rate and arterial pressure variabilities as a marker of sympathovagal interaction in man and conscious dog. Circ. Res. 59, 178–193. Pincus, S.M., Gladstone, I.M., Ehrenkranz, R.A., 1991. A regulatory statistic for medical data analysis. J. Clin. Monit. 7, 335–345.

Pomeranz, B., Macaulay, R.J.B., Caudill, M.A., Kutz, I., Adam, D., Gordon, D., Kilborn, K.M., Barger, A.C., Shannon, D.C., Cohen, R.J., Benson, H., 1985. Assessment of autonomic function in humans by heart rate spectral analysis. Am. J. Physiol. 248, H151–H153. Ramani, A., Grammaticos, B., Bountis, T., 1989. The Painleve property and singularity analysis of integrable and non-integrable systems. Phys. Rep. 180, 161–245. Rosenstien, M., Collins, J.J., De Luca, C.J., 1993. A practical method for calculating largest Lyapunov exponents from small data sets. Physica D 65, 117–134. Sano, M., Sawada, Y., 1985. Measurement of the Lyapunov spectrum from a chaotic time series. Phys. Rev. Lett. 55, 1082–1085. Sato, S., Sano, M., Sawada, Y., 1987. Practical methods of measuring the generalized dimension and largest Lyapunov exponents in highdimensional chaotic systems. Prog. Theor. Phys. 77, 1–5. Silipo, R., Deco, G., Vergassola, R., Gremigni, C., 1999. A characterization of HRV’s nonlinear hidden dynamics by means of Markov models. IEEE. Trans. Biomed. Eng. 46, 978–985. Soong, T.T., 1973. In: Bellman, R. (Ed.), Random Differential Equations in Science and Engineering. Academic Press. Storella, R.J., Wood, H.W., Mills, K.M., Kanters, J.K., Hojgaards, M.V., Holstein-Rathlou, N.H., 1998. Approximate entropy and point correlation dimension of heart rate variability in healthy subjects. Integr. Physiol. Behav. Sci. 33, 315–320. Takens, F., 1980. Detecting strange attractors in turbulence. In: Rand, D., Young, L. (Eds.). Dynamical Systems and Turbulence. Lecture Notes on Mathematics, Vol. 898. Springer-Verlag, Berlin, pp. 366–381. Theiler, J., 1987. Efficient algorithm for estimating the correlation dimension from a set of discrete points. Phys. Rev. A 36, 4456–4462. Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., Farmer, J.D., 1992. Testing for nonlinearity in time series: the method of surrogate data. Physica D 58, 77–94. Voss, A., Kurths, J., Kleiner, H.J., Witt, A., Wessel, N., Separin, P., Osterzeil, K.J., Schurath, R., Dietz, R., 1996. The application of methods of nonlinear dynamics for the improved and predictive recognition of patients threatened by sudden cardiac death. Cardiovasc. Res. 31, 419–433. Weise, F., Heydenreich, F., 1989. Effects of modified respiratory rhythm on heart rate variability during orthostatic load. Biomed. Biochim. Acta 8, 549–556. Wolf, A., Swift, J., Swinney, H., Vastano, J., 1985. Determining Lyapunov exponents from a time series. Physica D 16, 285–317. Yeragani, V.K., 1995. Heart rate and blood pressure variability: implications for psychiatric research. Neuropsychobiology 32, 182–191. Yeragani, V.K., Nadella, R., Hinze, B., Yeragani, S., Jampala, V.C., 2000. Nonlinear measures of heart period variability: decreased measures of symbolic dynamics in patients with panic disorder. Depression Anxiety, In press. Yeragani, V.K., Pohl, R., Berger, R., Balon, R., Ramesh, C., Glitz, D., Srinivasan, K., Weinberg, P., 1993a. Decreased heart rate variability in panic disorder patients: a study of power spectral analysis of heart rate. Psychiatry Res. 46, 89–103. Yeragani, V.K., Srinivasan, K., Vempati, S., Pohl, R., Balon, R., 1993b. Fractal dimension of heart rate time series: an effective measure of autonomic function. J. Appl. Physiol. 75, 2429–2438. Yeragani, V.K., Sobolewski, E., Kay, J., Jampala, V.C., Igel, G., 1997. Effect of age on long-term heart rate variability. Cardiovasc. Res. 35, 35–42. Yeragani, V.K., Sobolewski, E., Jampala, V.C., Kay, J., Yeragani, S., Igel, G., 1998. Fractal dimension and approximate entropy of heart period and heart rate: awake vs. sleep differences and methodological issues. Clin. Sci. 95, 295–301. Zeng, X., Eykholt, R., Pielke, R.A., 1991. Estimating the Lyapunovexponent spectrum from short time series of low precision. Phys. Rev. Lett. 66, 3229–3232.