Investigation of determinism in heart rate variability - Semantic Scholar

CHAOS

VOLUME 10, NUMBER 2

JUNE 2000

Investigation of determinism in heart rate variability M. E. D. Gomes ˜ o Lineares, Laboratorato´rio de Modelagem, Ana´lise e Controle de Sistemas Na and Departamento de Engenharia Eletroˆnica, Universidade Federal de Minas Gerais, Av. Antoˆnio Carlos 6627, 31270-010 Belo Horizonte, Brazil

A. V. P. Souza ˜ o Lineares, Laboratorato´rio de Modelagem, Ana´lise e Controle de Sistemas Na Universidade Federal de Minas Gerais, Av. Antoˆnio Carlos 6627, 31270-010 Belo Horizonte, Brazil

H. N. Guimara˜es Departamento de Engenharia Ele´trica, Universidade Federal de Minas Gerais, Av. Antoˆnio Carlos 6627, 31270-010 Belo Horizonte, Brazil

L. A. Aguirre ˜ o Lineares, Laboratorato´rio de Modelagem, Ana´lise e Controle de Sistemas Na and Departamento de Engenharia Eletroˆnica, Universidade Federal de Minas Gerais, Av. Antoˆnio Carlos 6627, 31270-010 Belo Horizonte, Brazil

共Received 2 August 1999; accepted for publication 15 March 2000兲 The article searches for the possible presence of determinism in heart rate variability 共HRV兲 signals by using a new approach based on NARMA 共nonlinear autoregressive moving average兲 modeling and free-run prediction. Thirty-three 256-point HRV time series obtained from Wistar rats submitted to different autonomic blockade protocols are considered, and a collection of surrogate data sets are generated from each one of them. These surrogate sequences are assumed to be nondeterministic and therefore they may not be predictable. The original HRV time series and related surrogates are submitted to NARMA modeling and prediction. Special attention has been paid to the problem of stationarity. The results consistently show that the surrogate data sets cannot be predicted better than the trivial predictor—the mean—while most of the HRV control sequences are predictable to a certain degree. This suggests that the normal HRV signals have a deterministic signature. The HRV time series derived from the autonomic blockade segments of the experimental protocols do not show the same predictability performance, albeit the physiological interpretation is not obvious. These results have important implications to the methodology of HRV analysis, indicating that techniques from nonlinear dynamics and deterministic chaos may be applied to elicit more information about the autonomic modulation of the cardiovascular activity. © 2000 American Institute of Physics. 关S1054-1500共00兲01302-1兴

I. INTRODUCTION

Heart rate variability is a signal with a great deal of information regarding the autonomic control of the cardiovascular system. The complexity of the related signals suggests that nonlinear dynamic phenomena may have an important role in cardiac activity. Actually, the literature has recently shown optimistic results when analysis methods from nonlinear dynamics and chaos theory were applied to such signals. The interpretation of such results, however, is much more difficult compared to those obtained with classical methods from the linear systems theory. On the other hand, if the time series has a deterministic signature, greater insight can be gained using techniques borrowed from the field of nonlinear and chaotic dynamics. This work discusses the application of nonlinear modeling and prediction techniques to search for traces of determinism in heart rate variability signals. As a by-product, further issues arise concerning the effect of drugs on the underlying dynamics. 1054-1500/2000/10(2)/398/13/$17.00

The regulation of the cardiac frequency, arterial pressure, and other cardiovascular variables is mediated mainly by the autonomic nervous system, together with intrinsic and extrinsic homeostatic mechanisms 共Levy, 1990; Levy and Schwartz, 1994; Berne and Levy, 1998兲. Complex variations are observed in these variables and the analysis of their variability has become an important noninvasive tool to study the sympatho-vagal interactions in many pathophysiological conditions, such as myocardial infarction and diabetic neuropathy 共Malik and Camm, 1995兲. Heart rate variability 共HRV兲 represents one of the most promising variables for such studies. During the last two decades, a great deal of work has been devoted to understanding the physiological information behind the variability of the cardiac cycle. Among the well established analysis tools, time and frequency domain methods from linear systems theory can provide valuable information for physiological interpretation of HRV and clinical 398

© 2000 American Institute of Physics

Downloaded 20 Jul 2001 to 150.164.35.48. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/chaos/chocr.jsp

Chaos, Vol. 10, No. 2, 2000

use 共Malik and Camm, 1995; Malik, 1998兲. However, nonlinear dynamics are certainly involved in the complex interactions among neural, humoral, and hemodynamic variables. It has been speculated that methods of nonlinear dynamics and deterministic chaos might provide a powerful tool to educe more information for better understanding the mechanisms of cardiovascular control 共Task Force European Soc. Cardiology & North American Soc. Pacing and Electrophysiology, 1996; Mansier et al., 1996; Griffith, 1996; Hoyer et al., 1996; Yambe et al., 1998兲. Since chaotic determinism for cardiac rhythm was first assumed 共Babloyantz and Destexhe, 1988; Goldberger, 1990; Goldberger et al., 1990; Kaplan and Goldberger, 1991兲, the literature has shown ambiguous conclusions 共Kaplan and Glass, 1992; Kanters et al., 1994; Guzzetti et al., 1996; Le Pape et al., 1997; Braun et al., 1998兲. According to the Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology 共1996兲, standards for analyzing HRV signals based on nonlinear systems theory and deterministic chaos ‘‘are lacking and the full scope of these methods cannot be assessed. Advances in technology and the interpretation of the results of nonlinear methods are needed before these methods are ready for physiological and clinical studies.’’ Even though this was stated in 1996, it is still a fact despite all the development in this field. Tools for signal analysis based on linear and nonlinear theory abound in the literature. However, the choice of the methods and the interpretation of the results can be facilitated if one knows, for example, that the underlying dynamics has some degree of determinism. It could be questioned if the chaos theory is really useful if the cardiac system is not manifestly deterministic. In this respect, Govindan et al. 共1998兲 applied surrogate and predictability analysis to study the dynamics of several normal and pathological electrocardiogram 共ECG兲 signals. They found some evidence that suggest deterministic chaos in ECG. Even though ECG signals have very important clinical applications, HRV captures the more intricate relation between the neural and cardiovascular systems. Furthermore, since HRV ‘‘looks’’ random, it is a challenging problem to look for determinism in such signal. It is the objective of this article to search for evidence of deterministic signature in HRV signals. It should be noted that no assumption is made regarding the dynamics of the system. By using surrogate data analysis 共Theiler et al., 1992兲 and NARMA modeling technique 共Aguirre and Billings, 1995b兲, a general framework is built to characterize this important aspect of the HRV time series. If determinism manifests, then appropriate tools, including some from deterministic chaos theory can be used to better understand the cardiac system. The results suggest the presence of determinism in the data in most cases. Moreover, when the autonomic nervous system is pharmacologically blocked at a cardiac level, the deterministic signature in the HRV signal seems to be reduced. This loss of determinism might have important physiological implications in the study of autonomic interactions in cardiac control, and the correct interpretation of this phenomenon is still open to further research. The outline of the article is as follows. Section II de-

Determinism in heart rate variability

399

scribes the data and discusses the test for nonstationarity. Analysis of determinism and surrogate data sets are briefly presented in Sec. III. Section IV describes the framework for NARMA 共nonlinear autoregressive moving average兲 modeling and prediction. The algorithm for testing determinism is described in Sec. V along with tests using a wide variety of well known time series. The results are discussed in Sec. VI and the main conclusions are summarized in Sec. VII. II. THE DATA A. Experimental protocols

1. Experimental animals

Experiments were conducted in eleven conscious unrestrained Wistar rats weighing between 210 to 285 g. All animals were housed individually in cages and a standard rat chow and water were available ad libitum. 2. Surgery

Rats were anesthetized with ether for the surgical procedures described below. The left femoral vein was catheterized 共PE10 sterile polyethylene catheter sealed to PE50 under hot air and filled with heparinized saline兲 to enable drug injections. The catheter was closed with a stainless pin when not in use. Silver electrodes were subcutaneously implanted for the measurement of standard lead II electrocardiogram 共ECG兲 signal. The distal end of the catheter and lead wires were tunneled under the skin to exit through a small incision in the scruff of the neck. This surgical procedure allowed the chronically catheterized and implanted rats unrestricted movement to all areas of the cage for the duration of the experiment. After completion of surgical procedures, a period of at least two days was allowed for the animal to recover. Proper placement of the catheter was confirmed at necropsy. 3. Extraction of the HRV signals

The ECG signal was conditioned by an electrocardiographic amplifier and analog-to-digital converted with a sampling frequency of 500 Hz and a resolution of 12 bits. The tachogram, defined as the sequence of RR interbeat intervals, was determined from the digitized ECG by taking the time interval between two successive maxima of R waves of the QRS complexes with an accuracy of 1 ms using an algorithm based on first derivative with an adaptive threshold. The automatic detection of R waves was visually checked to eliminate incorrect detection. This cardiac event series were finally transformed in a regularly sampled time series, with a sampling time of 250 ms, by an appropriate preprocessing technique based on a cubic polynomial interpolation of inverse interval function values. This technique generates a minimum of harmonic components and other artifacts and it ˜es and Santos 共1998兲 and in the is fully described by Guimara references therein. In the present study, we have considered the HRV as the characteristic signal of the oscillations between consecutive instantaneous heart rate around the mean heart rate. All artifact-corrupted data were rejected.


400

Gomes et al.

Chaos, Vol. 10, No. 2, 2000

FIG. 1. Strip plot of a 15 min HRV time series derived from the preprocessing of the tachogram for the 共a兲 control, 共b兲 atropine, and 共c兲 atenolol protocol sections for one representative rat 共rat #4兲.

4. Drugs

6. Data

Muscarinic blockade was performed with 2 mg/kg of atropine methyl sulfate and ␤ 1 -adrenergic blockade with 1 mg/kg of atenolol. All drugs were dissolved in sterileisotonic saline and were administered by intravenous 共i.v.兲 injection.

Finally, one short stationary 256-point segment was carefully chosen from each one of the above time series to compose the data of 33 stationary HRV time series. The tests of stationarity are discussed in the next section. Figure 2 shows one representative example of the stationary HRV time series for the protocols control, atropine, and atenolol. The power spectrum density 共PSD兲 of each HRV signal is shown to illustrate power alterations due to the pharmacological blockade. The PSD shows three important frequency bands, namely: 共a兲 ‘‘high-frequency band’’ or the respiratory sinus arrhythmia related band, which occurs at about the mean breathing frequency of the animal 共1.2 Hz兲; 共b兲 ‘‘lowfrequency band’’ or Mayer wave-like sinus arrhythmia related band, which occurs at about 0.1 Hz spontaneous vasomotor activity; and 共c兲 ‘‘very-low-frequency band,’’ which is the band of less than 0.05 Hz frequencies. Power alterations in the high-frequency band have been attributed to changes in cardiac parasympathetic activity, whereas in the lowfrequency band they have been related to changes in both

5. Experimental protocols

Thirty-three HRV time series of approximately 3600 points were taken from the following experimental protocol sections: 共i兲 five rats numbered #2–#6 submitted to the sequence control 共baseline condition兲, atropine, and atenolol, with 10 min interval between each phase and ECG recordings of 15 min at each phase; 共ii兲 six rats numbered #7–#12 submitted to the sequence control, atenolol, and atropine, with the same intervals between phases and same ECG recording time. Figure 1 shows a strip plot of one representative preprocessed HRV time series.


Chaos, Vol. 10, No. 2, 2000


401

FIG. 2. Stationary HRV time series and respective power spectrum density 共PSD兲 for the 共a兲 control, 共b兲 atropine, and 共c兲 atenolol protocol sections for one representative rate 共rat #4兲 that have shown a deterministic signature in its data. The spectral power in the ‘‘high-frequency band’’ 共about 1.2 Hz兲 is drastically reduced by the application of atropine, an antagonist for muscarinic cholinoceptors that blockade the vagal activity on the heart. The power in the ‘‘low-frequency band’’ 共about 0.1 Hz兲 have been related to changes in both sympathetic and parasympathetic activities. It becomes more prominent during the atropine administration. The application of atenolol, an antagonist for ␤ 1 -adrenoceptors, blockades the sympathetic activity at cardiac level suppressing the power in the lowfrequency band. Power remains in the ‘‘very-low-frequency band,’’ since other physiological factors influence this band.

sympathetic and parasympathetic activities. The physiological elucidation of the very-low-frequency band power is, however, less established, albeit some researchers have related it to thermoregulatory activity, to the renin-angiotensin system and to other systems that involve humoral factors 共Cerutti et al., 1991; Task Force European Soc. Cardiology & North American Soc. Pacing and Electrophysiology, 1996; Persson, 1997兲. B. Test of Stationarity

One requirement for the application of most linear or nonlinear signal processing techniques is stationarity. Changes in the dynamics during the measurement constitute an undesired complication for the analysis. Unfortunately it is not a trivial task to confidently test for stationarity. The literature has recent work on this matter 共Isliker and Kurths, 1993; Manuca and Savit, 1996; Schreiber, 1997兲. While the classical methods look for statistical fluctuations of some parameters in nonoverlapping segments S 1 , . . . , SL of the time series 共Bendat and Piersol, 1986; Priestley, 1988兲, these new methods focus on the dynamics underlying the time series. In this study, two techniques were used to select, from the measurements, a stationary HRV time series of length N⫽256. The techniques discussed below were extensively ˜es 共1996兲. Stationary used in a previous work by Guimara HRV time series selected from these methods have proven adequate for time and frequency analysis and the underlying dynamics extracted from the analysis were consistent with physiological interpretation. However, these tests of stationarity are not foolproof. It is important to note that physiological data are inherently nonstationary, and in the present work the data were obtained from conscious unrestrained rats in-

teracting with an ever-changing environment. The tests of stationarity search for segments that are as close to stationary as possible. Some positively-tested stationary HRV signals used in this work seem to show traces of nonstationarity and they were excluded from the analysis. A thorough discussion appears in Secs. V and VI. The technique proposed by Bendat and Piersol 共1986兲, known as reverse arrangements test, divides the time series y(k), 0⭐k⭐N⫺1, into L segments of length M, where the lth segment is defined as y l (m)⫽y(m⫹lM ) for 0⭐m⭐M ⫺1 and 0⭐l⭐L⫺1. The mean, ¯y l , and the standard deviation, ␴ l , of each segment are calculated. A set 兵 c l 其 is defined, where c l can be either ¯y l or ␴ l . A test function hi j⫽

再

1

if c i ⬎c j

0

otherwise,

共1兲

0⭐i⭐L⫺1, 0⭐ j⭐L⫺1, i⬍ j, is used to determine the total number of reverse arrangements, given by the random variable L⫺2

A⫽

兺

i⫽0

Ai ,

共2兲

where A i ⫽⌺ L⫺1 j⫽i⫹1 h i j . The time series is assumed stationary if A L;(1⫺ ␣ /2) ⬍A⬍A L;( ␣ /2) , where ␣ ⫽0.05 共confidence interval of 5%兲. In the second technique, the time series y(k), 0⭐K ⭐N⫺1, is divided into two nonoverlapping segments of equal length 共128 points兲. Autoregressive models of the same order are estimated for the two segments. The order is determined by Akaike criteria over the analysis of the first segment of the HRV time series. The poles and respective confidence elliptical regions associated with the first segment are plotted in the complex plane. Next, only the poles associated


402

Gomes et al.

Chaos, Vol. 10, No. 2, 2000

FIG. 3. Test of stationarity of two HRV time series: 共a兲 possibly stationary time series since 共b兲 the poles of AR models of two nonoverlapping segments are confined in a confidence region; 共d兲 in this case some poles are outside the confidence region, suggesting that the HRV time series shown in 共c兲 is nonstationary.

with the second segment are plotted in the same plane. If all the poles fall within the confidence regions, the time series is ˜es, assumed to be stationary 共Pagani et al., 1986; Guimara 1996兲. Figures 3共a兲, 3共b兲 show the result of this test applied to a possibly stationary HRV time series. The poles of the AR models estimated for two segments of 128 points are located inside the confidence regions, suggesting that the linear dynamics of these segments are similar. On the other hand, Fig. 3共d兲 shows some poles outside the confidence regions, suggesting that the HRV time series of Fig. 3共c兲 is not stationary. III. DETECTING DETERMINISM

A number of test statistics for the detection of determinism have been used in the literature. The method proposed by Kaplan and Glass 共1992兲 is adequate to quantify determinism in densely sampled data. Also the technique of false nearest-neighbors by Kennel et al. 共1992兲 can be used. Kaplan 共1994兲 proposed a method to test for determinism based on consistency with a continuous dynamical mapping. He argued that the technique could be applied to short data sets from moderate-dimensional noisy systems or lowdimensional systems with large Lyapunov exponents. Mansier et al. 共1996兲 applied a nonlinear prediction method 共Sugihara and May, 1990兲 and surrogate data sets to address the question if the HRV series is the output of a deterministic dynamical system. They showed that prediction index is better for the experimental series than for its surrogate series and suggested that these differences are thus in favor of a nonlinear deterministic system for the HRV times series under study. Chon et al. 共1997兲 proposed a method to test chaotic determinism based on fitting a nonlinear autoregressive model to the time series followed by the analysis of the characteristic exponents of the model over the observed probabil-

ity distribution of states for the system. They applied the method to HRV time series of 4096 points and their preliminary analysis suggested a nonchaotic deterministic component in these signals. Another approach was proposed by Barahona and Poon 共1996兲, which is similar to the method used in this work. They found nonlinear determinism in HRV signals of length 1000 points. This article discusses the similarities and peculiarities of their method in Sec. IV. The predictor used in the methods presented above usually require a time series with a large number of points, which is difficult to obtain with biomedical signals when stationarity is required. The approach presented in this article seems more appropriate for HRV signals, since NARMA modeling can be used with a very modest number of points and the surrogate data sets can be easily generated. A. Surrogate data sets

The method of surrogate data analysis 共Theiler et al., 1992兲 is being used in this work to test for determinism. The main idea behind this approach is as follows. Take the original times series and generate an artificial sequence that is assumed to be a realization of a Gaussian linear process or a class that is somewhat wider. Usually this new sequence retains some properties of the original signal, such as power spectrum. This ‘‘surrogate data’’ is, by definition, random and it represents the null hypothesis that should be rejected if the original signal has a deterministic signature. To discriminate the deterministic aspect of the time series, a prediction index is used after submitting the original time series and the surrogate data sets to nonlinear modeling and prediction. This method can be used to exclude certain classes of dynamics, but it does not allow for a more general characterization of the time series. In this study, phase randomized surrogate data and amplitude-adjusted surrogate data 共Theiler et al., 1992兲 are


Chaos, Vol. 10, No. 2, 2000


403

FIG. 4. The original HRV time series 共upper trace兲 along with its phase-randomized surrogate 共center兲 and its amplitude-adjusted surrogate 共lower trace兲.

generated from the HRV time series. This type of surrogates were adequate during the validation of the algorithm and there was not enough evidence in the HRV data to suggest the use of other types of surrogates. This, of course, does not exclude the idea of using different surrogates 共Kaplan and Glass, 1993兲 or even design a special one for HRV signals. An up-to-date discussion can be found in Schreiber 共1999兲. The phase randomized surrogate sets y s (k) are generated by first calculating the discrete time Fourier transform Y (n) of the HRV sequence y(k), where 0⭐n⭐N⫺1 and 0⭐k ⭐N⫺1. The phase information is randomized and then transformed back to the time domain by 1 y s共 k 兲 ⫽ N

N⫺1

兺

n⫽0

e j ␣ n 兩 Y 共 n 兲兩 e j2 ␲ nk/N ,

共3兲

where 0⭐ ␣ n ⬍2 ␲ are independent uniform random numbers. Therefore, this new sequence has the same secondorder properties as the original HRV signal, but it is otherwise random. The amplitude-adjusted surrogate sets are obtained by first generating a set of Gaussian random numbers x(k), 0 ⭐k⭐N⫺1. Then, the HRV times series y(k) is rescaled according to the Gaussian distribution so that the ranks of both time series agree; that is, if y(k) is the kth smallest of all the y ⬘ s, then x(k) will be the kth smallest of all the x ⬘ s. Next this new sequence is phase randomized by the previous algorithm. Finally, the surrogate set is created by again rescaling the last sequence by the histogram of the original HRV time series. Figure 4 shows one of the measured HRV signal 共upper trace兲, the corresponding phase randomized surrogate sequence 共center兲 and the amplitude-adjusted surrogate sequence 共lower trace兲. A visual impression from Fig. 4 is that complexity remains in the surrogate data.

IV. NARMA MODELING AND PREDICTION A. Model identification from data

Consider the nonlinear autoregressive moving average model 共NARMA兲共Leontaritis and Billings, 1985兲 y 共 k 兲 ⫽F l 关 y 共 k⫺1 兲 , . . . ,y 共 k⫺n y 兲 ,e 共 k 兲 , . . . ,e 共 k⫺n e 兲兴 , 共4兲 where n y and n e are the maximum lags considered for the process and noise terms, respectively, y(k) is a time series, e(k) accounts for uncertainties, possible noise, or unmodeled dynamics, and F l 关 • 兴 is some nonlinear function of y(k) and e(k) with degree of nonlinearity l苸Z⫹ . In this paper, the map F l 关 • 兴 is a polynomial of degree l. In order to estimate the parameters of this map, Eq. 共4兲 has to be expressed in prediction error form as y 共 k 兲 ⫽ ␺ T 共 k⫺1 兲 ␪ˆ ⫹ ␰ 共 k 兲 ,

共5兲

where ␰ (k) is white noise and

␺ T 共 k⫺1 兲 ⫽ 关 ␺ Ty 共 k⫺1 兲 ␺ Ty ␰ 共 k⫺1 兲 ␺ T␰ 共 k⫺1 兲兴 ␪ˆ ⫽ 关 ␪ˆ Ty ␪ˆ Ty ␰ ␪ˆ ␰T 兴 T

共6兲

and where ␺(k⫺1) is a vector which contains linear and nonlinear combinations of terms of the form y(k⫺i) and ␰ (k⫺ j) up to and including time k⫺1. It should be noted that the model output regressors are linear and nonlinear combinations of degree up to l of the variables 关 y(k ⫺1) . . . y(k⫺n y ) 兴 . Hence the model dynamical order is n y . The indication k⫺1 in Eq. 共5兲 is just a reminder that the latest information used by the model corresponds to time k ⫺1, thus emphasizing the predictor character of the model. The independent variables of model 共5兲, called regressors, for the sake of presentation were grouped in three sets, namely ␺y (k⫺1) which includes linear and nonlinear combinations of output terms only 共this is the deterministic part


404

Gomes et al.

Chaos, Vol. 10, No. 2, 2000

of the model兲, and 关 ␺ Ty ␰ (k⫺1) ␺ ␰T (k⫺1) 兴 T which includes nonlinear cross terms of output and residuals and linear and nonlinear combinations of residual terms 关 ␰ (k⫺ j) 兴 only. It is vital to realize that in making predictions only the deterministic part of the model will be used, that is yˆ (k) ⫽ ␺y ␪ˆ y . The stochastic part of the model is used only during parameter estimation as a means of drastically reducing parameter bias induced by noise in the observations. This is a standard procedure in system identification 共Ljung, 1987兲. If model 共5兲 is taken over a set of data with N observations, Eq. 共5兲 can be used to define a set of equations 共constraints兲 that involve the unknown vector of parameters ␪. Such a set of equations can be written in matrix form as y⫽⌿ ␪ˆ ⫹ ␰,

共7兲

where ⌿苸R(N⫺n y )⫻(n p ⫹n ␰ ) is known as the regressor matrix. Each column of this matrix is generated by taking each regressor in ⌿T (k⫺1) over the data. The vector of residuals, ␰, is defined as the difference between the vector of measured data, y, and the one-step-ahead prediction ⌿ ␪ˆ . The parameter vector ␪ can be estimated by orthogonal leastsquares techniques 共Chen et al., 1989兲. Hence parameter estimation is performed for a linear-in-the-parameters model which is closely related to Eq. 共5兲 and which is represented as n p ⫹n ␰

y 共 k 兲⫽

兺

i⫽1

g iw i共 k 兲 ⫹ ␰ 共 k 兲 ,

共8兲

where n p ⫹n ␰ is the number of 共process plus noise兲 terms in n p ⫹n ␰ are constant parameters and the regresthe model, 兵 g 其 i⫽1 n p ⫹n ␰ sors 兵 w i (k) 其 i⫽1 are constructed from the original monomials to be orthogonal over the data records. Finally, parameters of the model in Eq. 共5兲 can be calculated from the n p ⫹n ␰ . 兵 g i 其 i⫽1 A criterion for selecting the most important terms in the model can be devised as a by-product of the orthogonal parameter estimation procedure. The maximum mean-squared prediction error 共MSPE兲 of the signal is by definition y 2 (k), where the overbar indicates time averaging. Hence squaring Eq. 共8兲 and taking time average gives rise to squared terms of the form w 2i (k) and a great number of cross terms. Because of orthogonality only w 2i (k) and ␰ 2 (k) will not vanish. Therefore, the reduction in the MSPE due to the inclusion of the ith term, g i w i (k), in the auxiliary model of Eq. 共8兲 is (1/N)g 2i w 2i (k). Expressing this reduction as a percentage of the signal total MSPE yields the error reduction ratio 共ERR兲关 ERR兴 i ⬟

g 2i w 2i 共 k 兲 y 2共 k 兲

⫻100,

i⫽1,2, . . . ,n p ⫹n ␰ .

共9兲

The error reduction ratio, defined in Eq. 共9兲, yields a convenient way of ranking all possible candidate terms in decreasing order of importance. Hence those terms with large values of ERR are selected to form a parsimonious model. In practice Eq. 共9兲 is unsuitable for computation because it requires knowing the parameters g i that have to be estimated. In fact, a compact and very efficient algorithm can be implemented to rank the candidate regressors in descending order of im-

portance without needing prior knowledge of the parameters 共Korenberg et al., 1988兲. This procedure has two major advantages, namely 共i兲 it reduces inaccuracies due to numerical ill-conditioning, and 共ii兲 aids in selecting the structure of the final model. The choice as to where to stop including terms into the model has to be decided by other means, for instance using information criteria. A similar model representation, but with a different structure selection procedure, has been used to analyze cardiovascular signals 共Celka et al., 1999兲. B. Data prediction

As mentioned before, the main idea of testing for determinism using predictions and surrogates consists of using some kind of predictor to make forecasts over real and surrogate data. Because the surrogates are designed to be random, unless the predictions over the data are significantly different from the surrogate counterparts there will be no evidence of determinism in the data. For reasons to be discussed below, it is important to notice that the predictions over the real data need not be better but must be qualitatively different from those obtained over the surrogates. The idea of using predictions to assess data characteristics is not new. Sugihara and May 共1990兲 used predictions as a means to distinguish low-dimensional chaos from noise 共high-dimensional chaos兲. Casdagli 共1991兲 used model predictions to quantify the degree of determinism and nonlinearity in several sets of data. Schreiber 共1999兲 reviews prediction-based methods to detect nonlinear dynamics in data sequences. More recently, predictions have been used as a tool to detect determinism in physiological data 共Cao, 1997; Govindan et al., 1998兲. A common feature to most works is that the prediction used is based on local linear prediction techniques 共Farmer and Sidorowich, 1987兲 which typically require generous amounts of data in order to perform well. In this article a global nonlinear predictor of the form ˆ T (k⫺1) ␪ˆ is used to make forecasts over the meaˆy (k)⫽ ␺ y sured data and respective surrogates. In this predictor, ˆ T (k⫺1) contains previously predicted values of the signal, ␺ therefore y(k) is sometimes referred to as the free-run prediction of the model or as to the model predicted output. Among the advantages of this approach the following can be listed: 共i兲 data requirements are much less exacting than for other methods 共typically a few hundred observations are sufficient兲; 共ii兲 the predictors are very compact and can be written in closed form; and 共iii兲 predictors are globally nonlinear. Some important differences between local linear and globally nonlinear model structures have been commented in Kadtke et al. 共1993兲 and Aguirre 共1994兲. For each one of the 33 HRV 256-point time series, 16 random phase surrogates and 16 amplitude-adjusted surrogates were generated. NARMA models are then fitted to the first 200 points of each one of these 1089 sequences 共a total of 33 original HRV sequences and 1056 surrogates兲. The details of the procedure will be given in the next section. In order to compare model performance over the observed data and the surrogates, the following index was used:


Chaos, Vol. 10, No. 2, 2000

RMSE⫽

L 冑⌺ k⫽1 共 y 共 k 兲 ⫺yˆ 共 k 兲兲 2 , L ¯ 共 k 兲兲 2 冑⌺ k⫽1 共 y 共 k 兲 ⫺y


共10兲

where yˆ (k) is the free-run predicted signal, ¯y is the average value of the measured signal y(k) over an estimation data window of length L. The interpretation of Eq. 共10兲 is straightforward. This index compares the model-predicted output performance with the trivial mean value predictor. Values smaller than unity indicate improvement compared to the ‘‘trivial’’ predictor and suggest a greater degree of predictability of the time series. Therefore, for each training window 共of measured data or surrogates兲 a different predictor was estimated. Such predictors are, in a sense, optimal for the given number of model terms and maximum lags 共see below兲. It is also important to notice that for each training window, criterion 共9兲 was used to select the most important regressors. The performance of each predictor on the measured data and respective surrogates were then compared. The results are summarized in Sec. VI. A work that shares a number of common features with the procedure put forward in this paper has been published by Barahona and Poon 共1996兲. The model representation 共5兲 is the same and the index used to quantify predictor performance is very similar to Eq. 共10兲. The main differences of the procedures are the following. Firstly, the objective is different. Barahona and Poon 共1996兲 seek to establish if there are signs of nonlinearity in the data. To do this they fit linear and nonlinear models to the data and to surrogates. Unless the nonlinear models perform significantly better than the linear ones, the null hypothesis of a linear Gaussian process is not rejected. In this article we are concerned in finding determinism, which could turn out to be either linear or nonlinear. Secondly, Barahona and Poon 共1996兲 use an information criterion to select the structure of the fitted models. The order of the terms considered is predefined and is always the same. In this article we also use an information criterion to decide when to stop including terms in the model, but criterion 共9兲 is used to rank the candidate terms in order of importance. This is known to yield significantly better models. Thirdly, because the structure of our models is carefully chosen, free-run predictions are, in the vast majority of the cases, stable. This enables us to use such predictions to compare predictor performance in a way which is much more demanding and revealing. On the other hand 共presumably because the models structure is somewhat inadequate兲, Barahona and Poon 共1996兲 use one-step-ahead predictions to compare model performance. Finally, in this article, an autoregressive model is fitted to the data and a moving average is fitted to the residuals in order to ensure that ␰ (k) are white and therefore significantly reduce parameter bias. V. THE ALGORITHM FOR TESTING DETERMINISM

The procedure followed in this work basically tries to assess if the performance of predictors is significantly better for the measured data when compared to the surrogates. As

405

the latter are random by design, predictions should not be significantly better than the data average value, which can be considered a trivial predictor. If the performance of a predictor over measured data is significantly better from the average performance of predictions over the respective surrogates, then there is some evidence of determinism. The proposed algorithm consists basically of three stages: 共i兲 generation of surrogate data sets, 共ii兲 NARMA modeling, and 共iii兲 prediction and hypothesis testing. This section presents the complete procedure of the algorithm and illustrates how the algorithm performs with some well known systems.

A. Procedure

The following steps summarize the complete procedure. A few remarks are contained in each step to help choose the necessary parameters. • Step 1. Given a time series y(k) N1 , generate a surrogate data set. Phase randomized surrogate data and amplitudeadjusted surrogate data were considered in this work. Different surrogates may be used to model a class that is wider than the Gaussian linear process 共Theiler et al., 1992兲, and an ad hoc surrogate can be designed in specific cases 共Kaplan and Glass, 1993兲. Many realizations are recommended to validate the statistical test. In this article, 16 surrogates were generated for each class, giving a total of 32 for each time series to be tested. • Step 2. Define n y , n e , and l, the maximum lag for the process terms, the maximum lag for the noise terms, and the degree of nonlinearity, respectively. Usually, n y and n e were set to 5 and 2, respectively. If some information about the system is known, it is possible to adjust these variables to better values. For example, smaller n y works better for a logistic map. The degree of nonlinearity was set to 2 for all tests. This has shown to give more stable models with the data at hand. Furthermore, for l⬎2 no significant improvement in model performance was observed. • Step 3. Generate the candidate terms of the regressor matrix. In this article these terms were chosen as all possible linear and nonlinear combinations of degree up to l of the variables 关 y(k⫺1)•••y(k⫺n y ) 兴 and the linear combinations of the noise terms 关 ␰ (k⫺1)••• ␰ (k⫺n e ) 兴 . • Step 4. Divide the times series into two parts, namely (i) the identification window, and (ii) the prediction window. The time series used in this work are 256 points long: 200 points are considered for the identification window and the remaining 56 points for the prediction window. • Step 5. Obtain a family to NARMA models with number of terms ranging from n min to n max . The terms of each model are chosen among all possible candidates by using the ERR criterion. The models are nested since a model with n p terms has the n p ⫺1 terms of the previous model plus a new term selected using ERR. All parameters, on the other hand, are totally different for each model. In this article the following values were used n min⫽2 and n max⫽20.


406

Chaos, Vol. 10, No. 2, 2000

Gomes et al.

FIG. 5. Each line corresponds to one well-known time series. The first column shows the original time series, the second column one of the 32 surrogates 共random phase兲, and the third column presents the prediction performance, given by the prediction error RMSE, of the original time series 共solid line兲 and the interval of two standard deviation around the mean RMSE of the surrogates 共dotted line兲. Paradigm 1 is verified in 共a兲, 共b兲, and 共c兲. The RMSE curves for the original time series were different from the mean RMSE of the surrogates, showing better performance over some prediction time 关the scale of the ordinate in 共a兲 and 共c兲 starts at zero since the model fits very well and the RMSE approaches zero during some prediction time兴. According to the algorithm, the null hypothesis of randomness could be rejected, suggesting the presence of determinism in these time series. This is in agreement of the underlying dynamics. Paradigm 2 is verified in 共d兲 and 共e兲. The RMSE curve of the time series is most of the time inside the interval of 2 standard deviation, and no significant differences can be observed between the original and the surrogate series. Models cannot explain the dynamics of the time series better than the randomness in the surrogates and no determinism can be inferred.

• Step 6. From the (n max⫺nmin⫹1) NARMA models found in Step 5, select the best model by using an information criterion (e.g., Akaike). • Step 7. Use the model chosen in Step 6 to perform free-run predictions over the data prediction window. • Step 8. Calculate the RMSE for the original time series and for the surrogates. Compute the mean RMSE and the standard deviation for the surrogates. To illustrate the computation of the RMSE vector, recall that the prediction window used in this work has 56 points. A moving window of 16 points is used to calculate each RMSE value in the vector. For example, RMSE共1兲 takes 关 y(201)•••y(216) 兴 of the time series 共or equivalently the 16 first points of the prediction window兲 and 关 yˆ (201)•••yˆ (216) 兴 , the corresponding 16 predicted values of the prediction window. RMSE共2兲 takes 关 y(202)•••y(217) 兴 and 关 yˆ (202)•••yˆ (217) 兴 , and so on. Then, the last index of the RMSE vector is 41. • Step 9. Plot the RMSE vector of the original time series and the mean RMSE vector of the surrogates and compare the prediction performance. Note that no distinction between the two types of surrogates is made when comparing results. B. Validation

Before using the algorithm to investigate the presence of determinism in HRV time series, it was tested using a wide variety of well-known short time series (N⫽256). Model

identification was carried out with 200 points, and the remaining 56 were compared with 56 free-run predicted points according to the error criterion given by the RMSE. The first set of time series includes three deterministic systems: two continuous dynamical systems—Ro¨ssler and Lorentz—and one discrete map—the logistic equation. Since the time series taken from these systems are stationary and deterministic, it is expected that their prediction performance are superior to their respective surrogates, given a RMSE below unity for some prediction time. This scenario is referred to as the first paradigm in this article, namely: if the original time series is deterministic, the RMSE should be smaller than the mean RMSE of the surrogates 共usually it is approximately 1兲. The first paradigm regards the rejection of the null hypothesis of randomness and it is possible to suggest the presence of a deterministic signature in the data. The Ro¨ssler’s equation is defined 共Ro¨ssler, 1976兲 by a set of equations x˙ ⫽⫺(y⫹z), y˙ ⫽x⫹ ␣ y, and z˙ ⫽ ␣ ⫹z(x ⫺ ␮ ), with ␣ ⫽0.2 and ␮ ⫽5.7. The sampling time was chosen according to Aguirre 共1995兲. The time series was taken from the x-variable sampled at T s ⫽0.08. The Lorenz equations are given by x˙ ⫽ ␴ (y⫺x), y˙ ⫽ ␳ x ⫺y⫺xz, and z˙ ⫽xy⫺ ␤ z, with ␴ ⫽10, ␤ ⫽8/3, and ␳ ⫽28. This system was simulated and the time series was taken from the x-component sampled at T s ⫽0.04. The logistic equation y(k)⫽␭y(k⫺1) 关 1⫺y(k⫺1) 兴 was simulated for 11 values of the parameter ␭, ranging from 2.99 to 3.99 with increments of 0.1.


Chaos, Vol. 10, No. 2, 2000


407

FIG. 6. 共a兲 Time series from a logistic map with nonstationarity in the prediction window 共change of parameter ␭ in 15 points兲. The nonstationarity part of the prediction window can be identified by observing that the prediction performance becomes worse. 共b兲 Time series from a logistic map with nonstationarity in both the identification and prediction windows 共change of parameter ␭ in 15 points in each window兲. The prediction performance is worse in the nonstationary part of the prediction window. The nonstationarity in the identification window affects the overall performance of the predictor model.

Figures 5共a兲, 5共b兲, 5共c兲 show the time series, one example of its surrogates, and the prediction performance 共RMSE兲 for each one of these dynamical systems. It can be verified that the algorithm worked well with all of them, since the results were coherent with the deterministic characteristic of the data. The prediction performance observed in these cases is the scenario of paradigm 1. Now, two stochastic processes are considered: a uniformly distributed random sequence with zero mean and unit variance, and a Gaussian random process. Since these time series are stationary and stochastic, it is expected that their prediction performance would be similar to their respective surrogates. This means that each predicted time series of length 56 is around its mean, giving a RMSE around one, for both the original and surrogates time series. This is the second paradigm considered in this article: if the original time series is stochastic, the graphical representation of the RMSE should be around 共‘‘around’’ means plus–minus one standard deviation兲 the RMSE for the surrogates 共usually it is approximately 1兲. In this case, the null hypothesis of randomness cannot be rejected and it is not possible to advocate for a deterministic signature in the data. The algorithm was tested using 11 realizations of each process. The Gaussian random process can be represented by a randomly forced linear dynamical system x˙ ⫽Ax⫹Be and y⫽cx, where A and B are matrices of constants, c is a con-

stant vector, and e is a Gaussian distributed time series. In this article the following values were used: A⫽ 关 ⫺11; ⫺0.50兴 , B⫽ 关 0.5;0.5兴 , and c⫽ 关 10兴 . Figures 5共d兲, 5共e兲 show the time series, one example of its surrogates, and one representative example of the prediction performance 共RMSE兲 for each one of these two stochastic processes. The results were coherent with the randomness characteristic of the data. No significant differences could be observed between the prediction performance of the original and the surrogate time series in 100% of all realizations of these processes, which is the scenario of paradigm 2. Therefore, no traces of determinism could be verified by the test. To demonstrate the robustness of the algorithm in the presence of noise, the time series from the dynamical systems were corrupted with about 10% additive uncorrelated noise. The prediction performance suggested the presence of determinism in the data and demonstrated that paradigm 1 is still verified. Tests using the HRV signals, to be discussed thoroughly in the next section, revealed that a given predictor over measured data can be significantly worse than the performance over surrogates. A bad performance could be caused by the fact that the dynamics underlying the identification window are different from those underlying the prediction window. This is a common situation when the data are nonstationary and the models are nonlinear, because such models are quite

FIG. 7. This figure shows four cases that suggest determinism in the data. The prediction index is consistently better for the original HRV time series 共solid line兲 than the prediction index RMSE computed for the surrogates 共the mean of the index and an interval of one standard deviation is shown on both sides of the prediction curve for the surrogates—dotted line兲. These are typical scenarios of paradigm 1 in Table I. Observe how atenolol reduced the deterministic signature in the HRV signal of the rat #4. Notes: ctr-control protocol section; atatenolol protocol section; atr-atropine protocol section; ctr04-number identifies the rat.


408

Gomes et al.

Chaos, Vol. 10, No. 2, 2000

FIG. 8. This figure presents three examples where it is not possible to advocate for determinism in the data. The prediction index curve for the original HRV time series 共solid line兲 is most of the time inside the interval of one standard deviation from the mean of the surrogates 共dotted line兲, and no longer significant differences can be observed between the original and the surrogate series. Models cannot explain the dynamics of the HRV signal better than the randomness in the surrogates. These are typical scenarios of paradigm 2 in Table I. Notes: ctr-control protocol section; at-atenolol protocol section; atr-atropine protocol section; ctr02-number identifies the rat.

sensitive to changes in the underlying dynamics. To test such a situation, nonstationarity was introduced in the dynamics of a logistic map. A few nonstationary time series were simulated by introducing changes in the parameter ␭ in the prediction window and in both prediction and identification windows. The nonstationarity caused the nonlinear models to predict worse in the nonstationary parts of the prediction window. Figure 6 shows two examples of this. In 共a兲, change in ␭ was made in the prediction window and in 共b兲, change in ␭ was made in both prediction and identification windows. This scenario is referred to as the third paradigm in this article, namely: if the original time series is nonstationary, the RMSE can be erratic, changing from smaller to bigger than the mean RMSE of the surrogates. When paradigm 3 is present, we suspect of nonstationarity in the data even though this may not be the case in general. In any case, the present algorithm cannot be applied, since, as it was stated earlier, the method relies on stationarity. Furthermore, if the data are really nonstationary, the null hypothesis can incorrectly be rejected because the surrogates are stationary, even if the data are random 共Timmer, 1998兲. VI. TEST OF DETERMINISM IN HRV SIGNALS

The 33 HRV 256-point time series and respective 16 random phase surrogates and 16 amplitude-adjusted surrogates were submitted to the test of determinism, according to the procedure described in the previous section. The maximum lag for the process terms, n y , was set to 5 and the maximum lag for the noise terms, n e , was set to 2. The degree of nonlinearity was set to 2 for all time series. A degree of 3 did not give better performance and models became more unstable, due to overparametrization 共Aguirre and Billings, 1995a兲. Given n y , n e , and the degree of nonlinearity, 23 candidate terms were made available for the application of the ERR criterion 共9兲, in order to rank them in order of importance. n min and n max were set to 2 and 20, respectively. Then, 19 models were identified using an identification window of 200 points. The model with the minimum Akaike’s index was automatically selected as the predictor model. It is interesting to note that the models for the surrogates had on average more linear terms 共78.03%兲 than the models for the HRV time series 共60.61%兲. This observation agrees with the fact that the surrogates used in this work are real-

izations of a Gaussian linear process and reinforces the idea that the HRV signals have nonlinear components. The predictor was used to simulate the remaining 56 points of the time series and the RMSE index vector was calculated. The results are shown in Figs. 7 and 8. Note that the RMSE index for the surrogates is usually around 1 for all cases. An interesting feature seen in Fig. 7 is that the prediction performance over the data and surrogates tend to become similar after 16 points, which is equivalent to 4 s 共recall that the sampling time is 250 ms兲. It should be borne in mind that all forecasts are model free-run predictions. Because of the nature of the data, time correlations fade after 6 to 8 s. Consequently, the analysis should focus mainly on the first 4 s of forecasts which corresponds to, approximately, 22 heartbeats. Paradigm 1 discussed above is the scenario of Fig. 7, when measurement predictions are better than surrogates. As the predictor is moved along the complete data set, the performance deteriorates indicating that the underlying dynamics are moving apart from the dynamics underlying the iden-

TABLE I. The lines in this table refer to the protocol sections used to collect the HRV time series: control, atenolol, and atropine. Column ‘‘Paradigm 1’’ presents the percentage number of time series in each protocol section that show strong evidence of determinism 共Fig. 7兲. The prediction performance of the original HRV time series is superior to the prediction performance of the surrogates. The hypothesis of randomness in the data can be rejected. Autonomic blockade by atenolol and atropine seems to reduce the predictability of the HRV signals. The column ‘‘Paradigm 2’’ shows the percentage of the time series that appears to be random. It is not possible to advocate for a deterministic signature in the HRV signal since the performance of the prediction over measured data is similar to the performance over surrogates 共Fig. 8兲. The null hypothesis of randomness cannot be rejected in this case. Because a few HRV time series seems to be random, it is possible to suggest that the deterministic dynamics is predominant in the HRV time series. The column ‘‘Paradigm 3’’ shows the percentage number of time series that were left out of the analysis, since they may be nonstationary, even though they passed the stationarity tests discussed in the Sec. II B 共Fig. 9兲. Considering the valid HRV time series 共excluding those from paradigm 3兲, the column ‘‘Suggest Determinism’’ shows the percentage of the valid time series that suggest the presence of deterministic signature in the HRV signals. Protocol

Paradigm 1

Paradigm 2

Paradigm 3

Suggest Determinism

Control Atenolol Atropine

55% 36% 18%

9% 46% 64%

36% 18% 18%

86% 44% 22%


Chaos, Vol. 10, No. 2, 2000


409

FIG. 9. This figure shows two examples of paradigm 3. The model is not consistent in predicting the underlying dynamics. This might suggest that the time series was not stationary or that a good model could not be fitted to the data. Notes: ctr-control protocol section; atatenolol protocol section; ctr06-number identifies the rat.

tification window. These results suggesting determinism consistently appeared in 86% of the cases computed for control HRV sequences, that is, when the data were collected before the administration of any drugs 共baseline conditions兲. However, the results obtained for the HRV signals labeled atenolol and atropine, that is, date collected after the administration of autonomic blockade drugs suggest a loss of determinism. Similar results were also observed by Mansier et al. 共1996兲. After atropine injection, they observed a reduction of a predictability index using local linear predictors. The reduction of deterministic signature brings the time series closer to a stochastic signal. These results suggest that drugs 共atenolol and atropine兲 affect the cardiovascular system by reducing the degree of determinism in HRV signal. The efferent activity of the autonomic nervous system at a cardiac level is substantially blocked and it seems that the dynamics of the cardiovascular control system loses some determinism. This might have an important physiological interpretation but more investigation should be carried out to better understand this phenomenon. Paradigm 2 discussed above is the scenario of Fig. 8. These results do not suggest determinism in the data. It is not possible to advocate for deterministic signature in the HRV signal since the performance of the prediction over measured data is similar to the performance over surrogates. Only a small number of cases were classified according to this paradigm, suggesting that the determinism is prevalent in HRV signal. A few HRV time series were classified according to paradigm 3, as presented in Table I 共see Fig. 9 for two examples兲. Because paradigm 3 could be induced by nonstationarity, such time series were left out of the analysis, even though they passed the stationarity tests discussed above. More investigation is being carried out to elucidate more information about this paradigm. Table I presents a summary of the results. Recall that 33 HRV time series were analyzed: 11 for control, 11 for atenolol, and 11 for atropine protocol sections. For each protocol section, the results were classified according to the three paradigms discussed above. The column ‘‘Paradigm 1’’ presents the percentage number of time series in each protocol section that show strong evidence of determinism 共Fig. 7兲. Note, as pointed out, how the autonomic drugs influence the predictability of the HRV signals. Smaller number of HRV time series showing determinism were found for atenolol and atropine protocol sections. The column ‘‘Paradigm 2’’ shows the percentage of the time series in each protocol section that appears to be random 共Fig. 8兲. Because

a few time series seem to be random, it is possible to suggest that the deterministic dynamics is predominant in the HRV time series. ‘‘Paradigm 3’’ is included in the third column of Table I to inform how many time series were left out of the test, since they may be inappropriate by the possible presence of nonstationarity. VII. CONCLUSION

This article was concerned with the analysis of deterministic signature in heart rate variability 共HRV兲 time series. No assumption was made regarding the dynamics of the system. By using NARMA modeling and surrogate data analysis, traces of determinism in the dynamics of the HRV signals under study have been found. The work considered 33 HRV time series obtained from Wistar rats submitted to pharmacological autonomic blockade. The results have shown that the prediction index is good 共86% of the cases兲 before drug injection, which suggest the presence of determinism in the normal HRV signals. After pharmacological intervention, the prediction performance was no longer very different between the original and the surrogate series, which suggests that a deterministic signature is prevalent in the dynamics of a normal cardiac function. The apparent reduction of determinism for the cases, in which autonomic blockade drugs were administered, might have an interesting, but not obvious, physiological interpretation and more investigation in this direction is under way. ACKNOWLEDGMENT

M.E.D.G., A.V.P.S., and L.A.A. gratefully acknowledge financial support from CNPq 共Brazil兲. Aguirre, L. A. 共1994兲. ‘‘Some remarks on structure selection for nonlinear models,’’ Int. J. Bifurcation Chaos Appl. Sci. Eng. 4, 1707–1714. Aguirre, L. A. 共1995兲. ‘‘A nonlinear correlation function for selecting the delay time in dynamical reconstructions,’’ Phys. Lett. A 203, 88–94. Aguirre, L. A., and Billings, S. A. 共1995a兲. ‘‘Dynamical effects of overparametrization in nonlinear models,’’ Physica D 80, 26–40. Aguirre, L. A., and Billings, S. A. 共1995b兲. ‘‘Retrieving dynamical invariants from chaotic data using narmax models,’’ Int. J. Bifurcation Chaos Appl. Sci. Eng. 5, 449–474. Babloyantz, A., and Destexhe, A. 共1988兲. ‘‘Is the normal heart a periodic oscillator?,’’ Biol. Cybern. 58, 203–206. Barahona, M., and Poon, C. 共1966兲. ‘‘Detection of nonlinear dynamics in short, noisy time series,’’ Nature 共London兲 381, 215–217. Bendat, J. S., and Piersol, A. G. 共1986兲. Random Data: Analysis and Measurement Procedures 共Wiley, New York兲. Berne, R. M., and Levy, M. N. 共1998兲. Physiology 共Mosby, St. Louis兲. Braun, C., Kowallik, P., Freking, A., Hadeler, D., Kniffki, K., and Meesmann, M. 共1998兲. ‘‘Demonstration of nonlinear components in heart rate variability of healthy persons,’’ Am. J. Physiol. 275, H1577–H1584.


410

Chaos, Vol. 10, No. 2, 2000

Cao, L. 共1997兲. Strong evidence of determinism in multi-channel physiological data 共preprint兲. Celka, P., Vesin, J. M., Vetter, R., Grueter, R., Thoret, G., Pruvot, E., Duplain, M., and Scherrer, U. 共1999兲. ‘‘Parsimonious modeling of biomedical signals and systems: Application to the cardiovascular system,’’ in Nonlinear Biomedical Signal Processing, edited by M. Akay 共IEEE, New York兲. Cerutti, C., Gustin, M. P., Paultre, C. Z., Lo, M., Julien, C., Vincent, M., and Sassard, J. 共1991兲. ‘‘Autonomic nervous system and cardiovascular variability in rats: a spectral analysis approach,’’ Am. J. Physiol. 261, H1292–H1299. Chen, S., Billings, S. A., and Luo, W. 共1989兲. ‘‘Orthogonal least squares methods and their application to nonlinear system identification,’’ Int. J. Control 50, 1873–1896. Chon, K. H., Kanters, J. K., Cohen, R. J., and Holstein-Rathlou, N. 共1997兲. ‘‘Detection of chaotic determinism in times series from randomly forced maps,’’ Physica D 99, 471–486. Farmer, J. D., and Sidorowich, J. J. 共1987兲. ‘‘Predicting chaotic time series,’’ Phys. Rev. Lett. 59, 845–848. Goldberger, A. L. 共1990兲. ‘‘Nonlinear dynamics, fractals, and chaos: Applications to cardiac electrophysiology,’’ Ann. Biomed. Eng. 18, 195–198. Goldberger, A. L., Rigney, D. R., and West, B. J. 共1990兲. ‘‘Chaos and fractals in human physiology,’’ Sci. Am. 262, 42–49. Govindan, R. B., Narayanan, K., and Gopinathan, M. S. 共1998兲. ‘‘On the evidence of deterministic chaos in ECG: Surrogate and predictability analysis,’’ Chaos 8, 495–502. Griffith, T. M. 共1996兲. ‘‘Temporal chaos in the microcirculation,’’ Cardiovasc. Res. 31, 342–358. ˜es, H. N. 共1996兲. ‘‘Ana´lise da Variabilidade da Frequëˆncia Guimara Cardıáca-Me´todos e Implicações Fisiolo´gicas,’’ Ph.D. thesis, Instituto de Cieˆncias Biolo´gicas-Universidade Federal de Minas Gerais, Belo Horizonte, Brazil. ˜es, H. N., and Santos, R. A. S. 共1998兲. ‘‘A comparative analysis of Guimara preprocessing techniques of cardiac event series for the study of heart rhythm variability using simulated signals,’’ Braz. J. Med. Biol. Res. 31, 421–430. Guzzetti, S., Signorini, M. G., Cogliati, C., Mezzetti, S., Porta, A., Cerutti, S., and Malliani, A. 共1996兲. ‘‘Nonlinear dynamics and chaotic indices in heart rate variability of normal subjects and heart-transparented patients,’’ Cardiovasc. Res. 31, 441–446. Hoyer, D., Schmidt, K., Zwiener, U., and Bauer, R. 共1996兲. ‘‘Characterization of complex heart rate dynamics and their pharmacological disorders by nonlinear prediction and special data transformations,’’ Cardiovasc. Res. 31, 434–440. Isliker, H., and Kurths, J. 共1993兲. ‘‘A test for stationarity: Finding parts in a time series apt for correlation dimension estimates,’’ Int. J. Bifurcation Chaos Appl. Sci. Eng. 3, 1573. Kadtke, J. B., Brush, J., and Holzfuss, J. 共1993兲. ‘‘Global dynamical equations and Lyapunov exponents from noisy chaotic time series,’’ Int. J. Bifurcation Chaos Appl. Sci. Eng. 3, 607–616. Kanters, J. K., Holstein-Rathlou, N. H., and Agner, E. 共1994兲. ‘‘Lack of evidence for low-dimensional chaos in heart rate variability,’’ J. Cardiovasc. Electrophysiol. 5, 591–601. Kaplan, D. T. 共1994兲. ‘‘Exceptional events as evidence for determinism,’’ Physica D 73, 38–48. Kaplan, D. T., and Glass, L. 共1992兲. ‘‘Direct test for determinism in a time series,’’ Phys. Rev. Lett. 68, 427. Kaplan, D. T., and Glass, L. 共1993兲. ‘‘Coarse-grained embeddings of time series: random walks, Gaussian random processes, and deterministic chaos,’’ Physica D 64, 431–454. Kaplan, D. T., and Goldberger, A. L. 共1991兲. ‘‘Chaos in cardiology,’’ J. Cardiovasc. Electrophysiol. 2, 342–354.

Gomes et al. Kennel, M. B., Brown, R., and Abarbanel, H. D. I. 共1992兲. ‘‘Determining embedding dimensions for phase-space reconstruction using geometrical construction,’’ Phys. Rev. A 45, 3403–3411. Korenberg, M. J., Billings, S. A., Liu, Y. P., and McIlroy, P. J. 共1988兲. ‘‘Orthogonal parameter estimation algorithm for nonlinear stochastic systems,’’ Int. J. Control 48, 193–210. Le Pape, G., Giacomini, H., Swynghedauw, B., and Mansier, P. 共1997兲. ‘‘A statistical analysis of sequences of cardiac interbeat intervals does not support the chaos hypothesis,’’ J. Theor. Biol. 184, 123–131. Leontaritis, I. J., and Billings, S. A. 共1985兲. ‘‘Input–output parametric models for nonlinear systems. Part II: Stochastic nonlinear sytems,’’ Int. J. Control 41, 329–344. Levy, M. N. 共1990兲. ‘‘Autonomic interactions in cardiac control,’’ Ann. N. Y. Acad. Sci. 601, 209–221. Levy, M. N., and Schwartz, P. J. 共1994兲. Vagal Control of the Heart: Experimental Basis and Clinical Implications 共Futura, Armonk, NY兲. Ljung, L. 共1984兲. System Identification. Theory for the User 共Prentice Hall, Englewood Cliffs, NJ兲. Malik, M. 共1998兲. ‘‘Heart rate variability,’’ Curr. Opin. Cardiol. 13, 36–44. Malik, M., and Camm, A. J. 共1995兲. Heart Rate Variability 共Futura, Armonk, NY兲. Mansier, P., Clairambault, J., Charlotte, N., Me´digue, C., Vermeiren, C., LePape, G., Carre´, F., Gounaropoulou, A., and Swynghedauw, B. 共1996兲. ‘‘Linear and nonlinear analyses of heart rate variability: a minireview,’’ Cardiovasc. Res. 31, 371–379. Manuca, R., and Savit, R. 共1996兲. ‘‘Stationarity and nonstationarity in time series analysis,’’ Physica D 99, 134–161. Pagani, M., Lombardi, F., Fuzzetti, S., Rimoldi, O., Furlan, R., Pizzinelli, P., Sandrone, G., Malfatto, G., Dell’Orto, S., Piccaluga, E., Turiel, M., Baselli, G., Cerutti, S., and Malliani, A. 共1986兲. ‘‘Power spectral analysis of heart rate and arterial pressure variabilities as a marker of sympatho-vagal interaction in man and conscious dog,’’ Circ. Res. 59, 178–193. Persson, P. B. 共1997兲. ‘‘Spectrum analysis of cardiovascular time series,’’ Am. J. Physiol. 42, R1201–R1210. Priestley, M. B. 共1998兲. Nonlinear and Nonstationary Time Series Analysis 共Academic, London兲. Ro¨ssler, O. E. 共1976兲. ‘‘An equation for continuous chaos,’’ Phys. Lett. 57A, 397–398. Schreiber, T. 共1997兲. ‘‘Detecting and analyzing nonstationarity in a time series using nonlinear cross predictions.’’ Phys. Rev. Lett. 78, 843–846. Schreiber, T. 共1999兲. ‘‘Interdisciplinary application of nonlinear time series methods,’’ Phys. Rep. 308, 1–64. Sugihara, G., and May, R. M. 共1990兲. ‘‘Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series,’’ Nature 共London兲 344, 734–741. Task Force European Soc. Cardiology & North American Soc. Pacing and Electrophysiology 共1996兲. ‘‘Heart rate variability—standards of measurement, physiological interpretation, and clinical use,’’ Circulation 93, 1043–1065. Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., and Farmer, D. D. 共1992兲. ‘‘Testing for nonlinearity in time series: the method of surrogate data,’’ Physica D 58, 77–94. Timmer, J. 共1998兲. ‘‘Power of surrogate data testing with respect to nonstationarity,’’ Phys. Rev. E 58, 5153–5156. Yambe, T., Nanka, S., Kobayashi, S., Tanaka, A., Yoshizawa, M., Abe, K., Tabayashi, K., Takeda, H., and Nitta, S. 共1998兲. ‘‘Origin of chaos in the circulation: Open loop analysis with an artificial heart,’’ ASAIO J. 44, M700–M703.
