## Abstract

The aim of this study was to characterize how different nonlinear methods characterize heart rate and blood pressure dynamics in healthy subjects at rest. The randomized, placebo-controlled crossover study with intravenous terbutaline was designed to induce four different stationary states of cardiovascular regulation system. The R-R interval, systolic arterial blood pressure, and heart rate time series were analyzed with a set of methods including approximate entropy, sample entropy, Lempel-Ziv entropy, symbol dynamic entropy, cross-entropy, correlation dimension, fractal dimensions, and stationarity test. Results indicate that R-R interval and systolic arterial pressure subsystems are mutually connected but have different dynamic properties. In the drug-free state the subsystems share many common features. When the strength of the baroreflex feedback loop is modified with terbutaline, R-R interval and systolic blood pressure lose mutual synchrony and drift toward their inherent state of operation. In this state the R-R interval system is rather complex and irregular, but the blood pressure system is much simpler than in the drug-free state.

- nonlinear dynamics
- complexity
- dimensionality
- entropy

traditional linear analysis methods of heart rate and blood pressure time series data, such as the time- and frequency-domain methods, measure the strength of oscillations in heart rate and blood pressure within a specific frequency range, for example, in the low-frequency (0.04–0.15 Hz) and high-frequency (0.15–0.4 Hz) bands. The spectral powers obtained can be used, for example, to estimate sympathetic and parasympathetic nervous activity and, with cross-spectral approaches, to characterize arterial baroreflex functions. Generally, the linear methods (time- and frequency- domain methods) are useful and have been widely adopted in studies of health and disease because results from linear methods are quite easy to interpret in physiological terms. But they also have limitations, and criticisms have been raised against their use (5).

Multiple feedback loops in cardiovascular regulation systems make rapid adaptations possible under a large variety of physiological and environmental conditions. In analyzing heart rate and blood pressure time series with traditional linear methods, we lose a lot of information on the dynamic patterns used by the cardiovascular regulation systems to adjust heart rate and blood pressure. Linear methods have not been designed to yield information on the systems' inherent dynamic properties. Nonlinear methods of signal analysis can be more useful when characterizing complex dynamics. Thus the idea of using nonlinear statistics in the analysis of heart rate and blood pressure time series data is theoretically very sound and is a challenging objective for both cardiovascular physiologists and system theoreticians. Nonlinear analyses giving information on heart rate and blood pressure dynamics could potentially expand our knowledge on how the human cardiovascular regulation system operates.

Several nonlinear methods have already been used in clinical and experimental cardiovascular studies to characterize heart rate and blood pressure dynamics. The interpretation of nonlinear statistics, however, has been complicated frequently by the fact that the literature available on how nonlinear statistics describe cardiovascular changes induced by various physiological or environmental stimuli is scanty. Therefore, the clinical use of nonlinear methods is still limited, and physiologically relevant and accepted data interpretation requires basic research.

The first aim of this study was to characterize how different nonlinear analysis methods describe heart rate and blood pressure dynamics in healthy subjects at a resting state and how the nonlinear statistics describe cardiovascular changes induced by a potent and selective β_{2}-adrenoceptor agonist, terbutaline. Because terbutaline was used with constant infusion rates, the cardiovascular system could, after a certain transient period, be studied under stationary conditions. Therefore, it was possible to adjust the operation point of the regulation system. By exploring the changes in the nonlinear statistics of the regulation system during these infusion phases, we can gain new insight on the internal structure of the dynamics. The second aim of this study was to offer comparative data on the validity and robustness of the nonlinear methods used in the study and to present the need for simultaneously using several different methods to understand the information obtained adequately. Because we have used most of the nonlinear methods in time series analysis of heart rate and blood pressure in previous studies, we can also validate how well they can characterize inherent properties of the cardiovascular system.

## METHODS

#### Study design.

Six healthy male subjects with a mean age of 24 yr participated in the study. The study followed a randomized, placebo-controlled crossover design. The study protocol was approved by the ethics committee of the Turku University Hospital. Each subject was studied twice; the washout period between the two studies was 2 wk. Each study consisted of two randomized sessions: one session with terbutaline infusion and one placebo session. In each session there were four phases, with each study lasting 60 min. *Phase 1* was without any infusion (baseline), and it was followed by three similar phases (*phases 2–4*) with infusion. The drug infusion concentration was 5.0 mg of terbutaline in 50 ml of physiological saline. The infusion rates were 6 ml/h (10 μg terbutaline/min), 12 ml/h (20 μg terbutaline/min), and 18 ml/h (30 μg terbutaline/min). In the placebo session, physiological saline (50 ml) was used with the same infusion rates. For control purposes, terbutaline plasma concentrations were measured (13).

#### Data acquisition.

Electrocardiogram (ECG; cardioscope Olli Monitor 432, Kone Instrument Division), continuous noninvasive blood pressure from the left middle finger (Finapres 2300 NIBP Monitor, Ohmeda), and pulmonary air flow (M901, Medikro Oy) were measured and analog-to-digitally converted with a sampling rate of 200 Hz. The R-R interval (RRI) time series were generated by detecting the R peaks on the ECG signal, and the systolic arterial pressure (SAP) time series were generated by finding the maximal pressures on the blood pressure signal in each RRI. The pulmonary flow signal was used to check whether there were significant changes in the minute ventilation. Data acquisition was initiated 15 min after the start of each phase to have stabilized hemodynamics. During measurements, the breathing of subjects was paced (breathing rate was 0.25 Hz) with a sound signal.

#### Data analysis.

A set of different nonlinear methods was used to estimate the complexity [approximate entropy (ApEn), sample entropy (SampEn), and Lempel-Ziv entropy (LZEn)], fractal nature, dimensionality [correlation dimension (CorrDim)], oscillation patterns [symbol dynamic (SymDyn) entropy and number of forbidden words] and stationarity (StatAv) of RRI and SAP time series. All methods are suitable for analysis of time series of limited length. Mathematical descriptions of the methods are provided in the . In general, we used 800 data points, with the exception of a very few cases in which heart rate during the 10-min recordings included <800 beats. When the data series included >800 points, the section for the analysis was positioned in the middle of the time series. We excluded such nonlinear methods as deviation quantities of Poincare plots (12), which merely measure variability of the time series. All data analyses were done using special software (WinCPRS, Absolute Aliens Ay).

The RRI and SAP time series were analyzed with the same methods and equal parameters. Because all methods are independent of the units and absolute variability, the results of RRI and SAP time series analysis can be compared directly. Before entropy calculations, linear trends were removed. LZEn was calculated using the mean value as a threshold level. The results of nonlinear analysis are collected in Fig.1 for RRI and in Fig.2 for SAP. The differences between placebo and terbutaline treatments were analyzed by *t*-test for independent samples.

## RESULTS

Heart rate (expressed in terms of RRI) was constant during placebo infusions, whereas terbutaline significantly increased heart rate (Fig.1
*A*), as we reported previously (13). The changes of SAP (Fig. 2
*A*) were similar during placebo and terbutaline treatments. The finger cuff of the Finapres device probably altered the finger circulation observed as a small increase in the mean blood pressure. It should be noted that, although heart rate was almost doubled by terbutaline in *phase 4* compared with the baseline level values, the corresponding SAP changes were rather small.

The entropy values of SAP are systematically lower than those of RRI measured with ApEn, SampEn, or LZEn (Figs. 1, *B–D*, and2, *B–D*). Thus, in terms of different approaches of entropy, the SAP signals are more regular than the corresponding RRI signals. All entropy quantities for SAP time series also decreased with increased terbutaline effects (Figs. 1, *B–D*, and 2,*B–D*). The result shows that terbutaline simplifies the basic character of SAP regulation systems. In our previous work (13), we showed that terbutaline dose-dependently decreases the spectral powers of RRI but not those of SAP. The decreased entropy in SAP time series indicates that the complexity of the signal decreases but not the spectral powers.

The values of SampEn for RRI and SAP are always higher than the values of ApEn. The result is expected because SampEn does not include the comparison of the template vector with itself. Accordingly, the differences of SampEn values between placebo and terbutaline treatments, particularly in *phases 3* and *4*, are higher than the corresponding differences in ApEn values. We can assume that SampEn can be more sensitive to changes in the complexity of the signals. The differences between SampEn and ApEn values in our study, however, were not large enough to significantly improve statistical significance between placebo and terbutaline treatments.

It should be noted that during the last step of terbutaline infusion there is a trend for an increase, not a decrease, for both ApEn and SampEn of RRI. In ApEn and SampEn values for SAP, however, we can detect only a monotonically decreasing trend. Thus RRI and SAP time series data behave differently in terms of ApEn and SampEn, which indicates that these systems have different underlying dynamic structures. In LZEn values, however, we cannot detect the increasing trend in RRI. Obviously, the binary sequence in LZEn ignores a lot of information, and the LZEn method is too rough to catch the trend: when the absolute variability of the signal is low, a binary sequence cannot catch fine details. Because ApEn and SampEn use the original data values of the time series, ApEn and SampEn methods are less prone to lose discreteness. The advantage of LZEn is that the method is nonparametric and fast to compute. However, the method requires that the variability of the signal under study be relatively large.

Terbutaline affects CorrDim of RRI and SAP (Figs. 1
*E* and2
*E*) differently. CorrDim of RRI decreases already in*phase 2*, is almost at the same level in *phase 3*, and returns back to the original level in *phase 4*. On the other hand, CorrDim of SAP does not change in *phase 2*, drops dramatically in *phase 3*, and remains low in *phase 4*. The RRI regulation system can restore the original dimensionality at the highest rates of terbutaline infusions. The phenomenon parallels with the findings on ApEn and SampEn values, which also show a paradoxical increasing trend in *phase 4*. If the explanation for the findings is the true changes in the dimensionality of the systems (attractor in phase space), CorrDim can be regarded as a more sensitive measure, because it is more directly related to system structure. CorrDim values for SAP remain low at the maximal terbutaline infusion rates. The phenomenon is found also in the corresponding entropy values. The lowest CorrDim values can be found already in*phase 3*. This phenomenon could be interpreted as a sign of saturation. Terbutaline is able to maximally simplify the SAP regulation system already in *phase 3*.

Fractal dimension calculated by curve lengths (FD-L) of both RRI and SAP decreases slightly with increasing terbutaline dosage (Figs.1
*F* and 2
*F*). The actual decreases are too small to characterize changes in the RRI or SAP regulation systems. The methods based on fractal dispersion analysis (FD-DA) are also unsuitable to characterize terbutaline-induced changes in the RRI and SAP regulation systems. The most prominent changes are found during placebo infusions, whereas the measures are unaffected by terbutaline infusions.

Our data on RRI and SAP symbol dynamics reveal that the responses of RRI and SAP to increasing doses of terbutaline differ greatly from one another. There are no significant differences in the baseline level of SymDyn entropy of RRI and SAP (Figs. 1
*H* and 2
*H*). However, the percentage of forbidden words (Figs. 1
*I* and2
*I*) is much higher for SAP than for RRI. The result shows that there are fewer patterns in SAP regulation than in RRI regulation. This finding is consistent with results on ApEn, SampEn, and LZEn, although the number of forbidden patterns measures quite different features of regulation systems.

With the increased terbutaline effects, the SymDyn entropy of RRI does not change but the percentage of forbidden words (patterns) increases dramatically. This finding can be interpreted to mean that there are only a few dominating patterns in the RRI regulation system at rest and a quite large set of patterns that are used rarely. The condition was explored by studying the distribution of the words; we could detect only a few high peaks and a widespread low-level background. Terbutaline totally abolishes the most seldom-used words, and the percentage of forbidden words increases dramatically. However, SymDyn entropy does not change, because the seldom-used words (patterns) contribute only nonsignificantly to the total SymDyn entropy and the dominant words contributing to the SymDyn entropy are unaffected by terbutaline treatment.

Terbutaline has no effects on the number of forbidden words in the SAP time series data, as shown in Fig. 2
*I*. This result is due to the fact that in the SAP time series there are only a few patterns and these patterns do not disappear with terbutaline treatment. The significant decreases in SymDyn entropy (Fig. 2
*H*) also require explanation. Possibly, there are changes in the few dominant patterns contributing to the total SymDyn entropy. Apparently, the system dynamics concentrates only on a very few patterns that result in more regular behavior and in decreases in SymDyn entropy.

StatAv analysis reveals changes only in RRI signals (Fig.1
*J*), whereas StatAv of SAP signals is unchanged by terbutaline (Fig. 2
*J*). StatAv measures the baseline shifts of the signals coarsely. Therefore, it is difficult to explain why the StatAv of RRI signal is increased by terbutaline. The explanation may be the fact that the result also depends on the segment length, which was 20 data points in our analyses. In theory, higher StatAv values indicate increased irregularity. However, StatAv, unlike entropy measures, is insensitive to rapid changes in the signal. As symbol dynamic analyses show, the RRI regulation system utilizes only a few patterns in *phase 4*, and the patterns utilized were the dominant patterns of the baseline. The significant decreases in the number of patterns cannot explain the increases of irregularity. However, at the same time, the regulation systems become more complex in terms of CorrDim, SampEn, and ApEn. ApEn and SampEn characterize the conditional probability of the patterns found on the signal, whereas they do not measure how frequently patterns evolve. SymDyn entropy, on the other hand, measures only the relative probability of the patterns. Thus the most obvious explanation for the increases of StatAv is that the few utilized patterns emerge more randomly in the time scale of 20 beats because of the more complex dynamics.

The RRI/SAP cross-entropy (Fig. 1
*K*) in *phases 3*and *4* decreases statistically significantly from baseline values. The result indicates that certain regulation patterns are more likely to occur in both signals during the last two phases of terbutaline infusions than before drug administration.

## DISCUSSION

In this study we used a set of different nonlinear methods to characterize heart rate and blood pressure dynamics in healthy volunteers before and after terbutaline administration. Heart rate and blood pressure time series can be regarded to represent complex system dynamics, because both heart rate and blood pressure are regulated via neural and hormonal feedback mechanisms including various types of baro-, volume- and chemoreceptors located in the heart, great arteries and veins, lungs, and brain (11). For example, the parasympathetic and sympathetic nervous systems regulate the beat-to-beat control of blood pressure and heart rate via sinoaortic baroreflex actions (15). The nonlinear analysis methods are superior to the traditional linear methods in characterizing complex system dynamics. Thus nonlinear analyses of heart rate and blood pressure time series could not only complement the traditional time- and frequency-domain analyses but also yield essential information on human heart rate and blood pressure dynamics.

Terbutaline is a potent and selective β_{2}-adrenoceptor agonist that interferes strongly with the regulation of arterial blood pressure, heart rate, and arterial baroreflex functions. Effects of terbutaline on the cardiovascular regulation are mediated mainly via activation of vascular, neuronal, and cardiac β_{2}-adrenoceptors. The activation dilates blood vessels particularly in the skeletal muscle, interferes with the release of various neurotransmitters and hormones, increases heart rate, and, at low doses, induces a reflex activation of arterial baroreflex. High terbutaline doses drastically decrease cardiac parasympathetic modulation and dampen arterial baroreceptor functions. Thus the three-step intravenous terbutaline infusions used in this study induce large dose-dependent changes in the cardiovascular regulation of healthy volunteers.

We used several different measures of entropy to estimate the complexity of our RRI time series data. The significant decrease of ApEn, SampEn, and LZEn could be well explained by the terbutaline-induced decrease in the parasympathetic modulation of heart rate. It is well known that cardiac respiratory sinus arrhythmia is related to cardiac parasympathetic modulation, and sinus arrhythmia explains a major fraction of human heart rate variability at rest (8, 18). SymDyn entropy, however, was surprisingly not decreased by terbutaline, whereas the another symbol dynamic measure, the percentage of forbidden words, was drastically increased by terbutaline. The symbol dynamic indices measure the probability of certain dynamic patterns within the RRI time series data occurring or not occurring. The large increases of the percentage of forbidden words could be explained by the uncoupling of arterial baroreceptor functions (13) and decreased cardiac parasympathetic modulation. We discussed in detail in results why terbutaline does not decrease SymDyn entropy of RRI. The paradoxical increases of ApEn and SampEn by the last terbutaline infusion step are hard to explain on the basis of available data. The effect could reflect the loss of selectivity of terbutaline action. Terbutaline at high concentrations may act nonselectively and activate other β-adrenoceptor subtypes and possibly activate even dopamine or serotonin receptor subtypes. The loss of selectivity could produce more complex RRI time series. High doses of terbutaline may also strongly interfere with central nervous system functions and peripheral cholinergic nerve transmission. Thus high terbutaline doses could explain the increases of RRI SampEn and ApEn. On the other hand, the effect may reflect physiological adaptation to the long-lasting and strong terbutaline effects: the attempts of the regulation systems to normalize heart rate dynamics (17). If this explanation is valid, the relationships between RRI entropy and cardiac parasympathetic modulation are complex.

Terbutaline significantly decreased the complexity of SAP time series data expressed in ApEn, SampEn, or LZEn. In contrast to RRI time series data, terbutaline also significantly decreased SymDyn entropy. The effect can be well explained by terbutaline-induced decreases in cardiac parasympathetic modulation and uncoupling of arterial baroreceptor functions (13). However, terbutaline had no effect on the percentage of forbidden words. Thus, in symbol dynamic terms, the dynamics of blood pressure differs greatly from the dynamics of blood pressure. The mathematical background of this interesting finding is explained and discussed in detail in results.

In general, our entropy measurements indicate that the regulation of SAP is much simpler than the regulation of heart rate. The time courses of ApEn and SampEn of RRI after terbutaline administration also differ from the time courses of SAP. This effect was studied more thoroughly by CorrDim.

CorrDim computations of time series data estimate the nature of regulation dynamics. As discussed in results, our CorrDim results are in good conformity with the entropy results of RRI ApEn and SampEn and SAP ApEn and SampEn. Thus the physiological backgrounds explaining the findings are likely to be at least partly similar. For example, the abrupt decreases of CorrDim and SampEn of SAP by terbutaline could well reflect the uncoupling of arterial baroreflex functions and large decreases in cardiac parasympathetic activity. Arterial baroreflex actions may represent a “dimension” in blood pressure time series data, a dimension abolished by terbutaline.

If we compare the results obtained with the methods measuring entropy, we can conclude that all entropy quantities (ApEn, SampEn, and LZEn) revealed the main changes of the system due to terbutaline treatment. From the mathematical point of view, SampEn should be more sensitive than ApEn in detecting changes in entropy. In this study, however, the sensitivity of methods was very similar. Obviously, the better sensitivity of SampEn is expressed in analyzing very short time series. However, one should very careful in drawing any conclusions on system features when the time series include only a few data points. The LZEn results confirmed the results obtained with the ApEn and SampEn methods, because the computational approach of LZEn differs totally from that of ApEn and SampEn. The low absolute variability of heart rate and blood pressure time series, however, can limit the usefulness of the LZEn method. CorrDim, symbol dynamics quantities (SymDyn entropy and forbidden words), and stationarity analysis (StatAv) yield truly independent information on the cardiovascular regulation system. The methods based on fractal image analyses were disappointing: neither FD-L nor FD-DA methods could show any relevant changes. Methods based on image analysis are obviously not useful to reconstruct the dynamics of the underlying system; they perhaps merely reflect changes in the variability of the signals. There are some reports in the literature that these methods could be useful in characterizing cardiovascular regulations systems, but our results do not support those suggestions. We were unable to find any relevant physiological explanation for the results obtained with these methods, and therefore we must conclude that they have severe limitations in time series analysis of RRI and SAP.

Our results show that the regulation system consists of two mutually connected heart rate and blood pressure subsystems both of which have unique features. The system theoretical relations between the essential parts of the heart rate and arterial blood pressure regulation systems are schematically presented in Fig. 3. When the baroreflex is damped, RRI can oscillate or merely wander (there is some variation in RRI because entropies are clearly positive but the amplitude of the oscillation is much lower than normal) without any control of the pressure systems. However, arterial baroreflex still controls SAP by adjusting peripheral vascular resistance. This subsystem has a rather simple dynamics. However, the subsystem can obviously regulate SAP. Although the RRI regulation systems work without baroreflex control, cardiac output still affects peripheral vascular resistance. Thus one should also measure cardiac output and total peripheral vascular resistance to explore human blood pressure dynamics more thoroughly.

Before terbutaline administration, the subsystems are tightly linked together and they share many similar dynamic properties. When the mutual interaction has been weakened or cut off completely by terbutaline, the regulation systems drift to a state in which blood pressure does not affect heart rate. The more independent RRI and SAP subsystems show new modes of function. In terms of correlation dimension and entropy, the new RRI state is similar to the state before terbutaline administration. However, the new and old states differ in terms of stationarity. The new SAP state, however, is rather simple in its coherent dynamic patterns. Thus the regulation of heart rate also remains complex after baroreflex failure. The simple inherent dynamics of SAP regulation systems may reflect the restrictions required to keep arterial blood pressure within narrow physiological limits. On the other hand, the pressure part of the system under normal conditions can control the RRI subsystem easily because it is working close to its open-loop condition.

Linear cross-spectral analysis of RRI and SAP time series yields information on arterial baroreflex actions. In this study we analyzed RRI-SAP interactions using a cross-entropy approach. Cross-entropy measures the probability of finding similar patterns in RRI and SAP time series. Terbutaline significantly decreased cross-entropy. Thus the pattern architecture of the time series becomes more similar after terbutaline administration. Terbutaline, however, is known to dampen arterial baroreflex actions and to desynchronize RRI/SAP behavior with respect to time. Possibly, when baroreflex actions are weakened, other external factors can affect RRI and SAP regulation more powerfully and similarly. We believe that this external force is respiration, which can affect both subsystems directly.

In conclusion, the analysis of heart rate and blood pressure time series with nonlinear methods reveals essential features of human heart rate and blood pressure dynamics. However, because of the complexity of cardiovascular regulation systems we must use a set of different methods to obtain the information required. The heart rate and arterial blood pressure regulating systems share many common dynamic features when the systems are linked to each other by baroreflex actions. However, after baroreflex dampening, the systems show highly different inherent dynamic structures. The heart rate regulation systems are complex and irregular, whereas the blood pressure regulation systems are rather simple.

## Acknowledgments

We thank Dr. Jarmo Hietarinta (University of Turku, Department of Physics) for valuable comments and suggestions.

## Appendix

#### Approximate entropy.

ApEn measures the complexity of signals: ApEn is low in regular signals and high in irregular signals. ApEn also determines the probability of finding specific patterns in the time series (2, 4, 19,21). When ApEn is calculated, the template patterns are constructed from the signal itself, and no a priori knowledge of the systems is needed.

In calculating ApEn, we first construct from the original time series (RRI or SAP in this study), {*x*(*i*)},*i* = 1 … *N*, vectors of the pseudo-phase space using time lags
where *m* is the embedding dimension and τ is the time lag. The optimal τ value can be found by autocorrelation or mutual information method (7). For RRI and SAP time series, a τ of 1 is usually the proper choice. The distance of the vectors **u**(*i*) is determined as a maximal distance*d* of the corresponding components in the vector (we have now set τ = 1)
The second index in *u* corresponds to the vector component. The probability of finding another vector within the distance *r* from the template vector**u**(*i*) is given by
If we use the logarithmic probability
approximate entropy is defined by
This expression can be interpreted as the (logarithmic) conditional probability that vectors, which are close in *m*points, are also close if lengthened to *m* + 1 points.

ApEn depends on three parameters: *m*,* r*, and*N*. In RRI and SAP time series, the embedding dimension*m* cannot be big, and *m* was set to be 2. ApEn depends strongly on the filter parameter* r*, and normally it has a maximum at some value of *r*. It is useful to set the filter parameter *r* to be some percentage of the standard deviation (SD) of the time series because then ApEn does not depend (or depends only slightly) on the absolute variability or the unit of the signal. A commonly used choice is *r* = 20% of SD. If*r* is very small, noise in the signal affects the ApEn values. When very large *r* values are used, essential properties of the signal remain undetected. The number of data points*N* is not critical when the number exceeds 800 points, because ApEn converges asymptotically at the final level when*N* increases. ApEn is very sensitive to linear trends: in distance calculations we compare the values of the coordinates directly. Thus regular signals with small linear trends can be interpreted to be much more irregular than they actually are, because of high entropy values. Therefore, we removed the obvious linear trends from the time series before computing ApEn.

#### Sample entropy.

The basic idea of the SampEn is very similar to the ApEn, but there is a small computational difference (23). In ApEn, the comparison between the template vector and the rest of the vectors also includes comparison with itself. This guarantees that probabilities*C*
(*r*) are never zero. Thus it is always possible to take a logarithm of probabilities. Because template comparisons with itself lower ApEn values, the signals are interpreted to be more regular than they actually are.

With SampEn, the vector comparison with itself is removed
Now we can determine
and finally
SampEn measures complexity of the signal in the same manner as ApEn. However, the dependence on the parameters *N *and*r* is different. SampEn decreases monotonically when*r* increases. In theory, SampEn does not depend on*N*. In analyzing time series including <200 data points, however, the confidence interval of the results is unacceptably large. When* r* and *N* are large, SampEn and ApEn give the same results.

#### Cross-entropies.

Entropy can also be calculated between two different signals, for example, between RRI and SAP signals. This mutual entropy characterizes the probability to find similar patterns within the signals (20,22, 23). Cross-ApEn describes the synchronicity of patterns; it does not measure the time synchronization of the signals. If Cross-ApEn is low the signals have common features in the pattern architecture, and if Cross-ApEn is high the signals differ in their pattern architecture.

When calculating the cross-entropies, the vectors that are compared are taken in pairs from different time series. If we have two time series {*x*(*i*)} and {*y*(*i*)},* i* = 1… *N*, the pseudo-phase vectors are constructed as before
and the vector distance is now
With this definition of distance, we can again compute ApEn and SampEn. However, because normally two time series represent different quantities, as for RRI and SAP time series, they must first be normalized. The most natural normalization method is to subtract the mean value from each data series and then divide it by its SD. Because we compare patterns this normalization is valid.

In our time series analyses, the ApEn approach is not always useful because it is quite possible that no vectors from the second time series can be found that are within the distance *r* to the template vector taken from the first time series leading to the terms of log(0). With the SampEn approach, however, we need only one close vector, and in most cases Cross-SampEn entropy can be calculated. Cross-SampEn gives the same result regardless of the order of the signals, but Cross-ApEn depends on the order. We used Cross-SampEn with the same parameter values as in the case of normal SampEn (*m* = 2, *N* = 1,000,*r* = 20% of SD).

#### Modified CorrDim.

If the dynamic system has a well-defined attractor (the asymptotic path of the system) in the phase space, the dimension of this attractor can reveal the nature of the system. The CorrDim computation is one method to characterize the dimension of this attractor. It should be noted that chaotic systems can have fractal attractors, and therefore the dimension does not have to be an integer. The existence of attractor is not always clear, and noninteger CorrDim does not prove that system is chaotic.

CorrDim can be computed by first calculating the number of vectors within the distance *r* from a given vector in the pseudo-phase space (6, 10, 14). Vectors are constructed as before. The vector number can be expressed as
where *D* is the euclidean distance
When we calculate the average over all vectors, we get the correlation integral
Finally, CorrDim is determined as limit value
These limits cannot be calculated in practice. Instead, one fits a regression line over a reasonable linear region in the log-log plot, and the slope is interpreted as CorrDim.

The computation requires that the embedding dimension *m* is at least 2 × *D*, where *D* is the dimension of the system (*D* estimates the number of independent variables) and the number of data points is at least 10^{m}. If one assumes that the blood pressure regulation system includes, for example, four variables, the time series required to perform adequate calculations should be very long. It also quite difficult to find the range where the relation between log *r* and log*C ^{m}
*(

*r*) is linear because of noise or nonstationarities of the signal. However, a useful estimate of CorrDim, called modified CorrDim, can be computed by using the fixed embedding dimension

*m*= 20 and the time series length including ∼1,000 data points and searching the slope of the regression line on the log-log plot over the range 0.01 <

*C*(

^{m}*r*) < 0.1 (26). We have also used the requirement that the correlation coefficient of the regression line must be >0.8.

#### Fractal dimension by dispersion analysis.

The dimension analysis can also be performed by studying the time series curve itself without trying to reconstruct the dynamic system. This kind of approach is close to image analysis. We can assume that complex dynamic behavior can be seen as apparently unpredictable changes in the signal. Therefore, the fractal analysis of the signal curve can give indirect information about the underlying system.

The method of dispersion analysis is based on the standard deviation
where {x(*i*)} is the time series of *N*data points (1, 26). In the next step, a new time series is constructed that consists of the mean values of two adjacent data points; this new time series has now *N*/2 data points and its standard deviation is SD(2). By continuing this procedure with group sizes of 4, 16, 32, and so on (we must stop when the number of points is <4), we get the series of deviations SD(*m*). By plotting log SD(*m*) as a function of log*m*, the points should be positioned on a straight line. The fractal dimension (FD-DA) is calculated from the relation FD-DA = 1 − slope of the line. FD-DA can have values between 1.0 and 1.5, where a value of 1 corresponds to a constant signal and a value of 1.5 corresponds to a maximal fractal structure (= totally random signal). If the time series contains noise, the SD(*m*) values are not exactly on the line in the log-log plot but the slope can be found by fitting the regression line on the data. Also in this case, we have set an additional condition that the correlation coefficient must be >0.8 for valid estimation of FD-DA.

It should be noted that FD-DA does not depend on the actual unit (and the absolute variability) of the time series because the scaling of the signal is transformed to the offset term in the log-log plot. Unfortunately, FD-DA is quite sensitive to noise, measurement inaccuracies, and nonlinear trends, which restricts the use of this method.

#### Fractal dimension by curve lengths.

This method also characterizes the fractal dimension of the signal curve. This is done by counting the number of sticks of various lengths needed to follow the curve. When the length of the stick is decreased, the shape of the curve can be followed better and, naturally, more sticks are needed (3, 9). If the curve has fractal nature the number of sticks increases exponentially. If*N*(*L*) is the number of the sticks of the length*L*, the fractal dimension is determined as a limit value
for a straight line FD-L = 1. In practice, FD-L is determined as a slope of the regression line in the log-log plot. FD-L depends on the unit of the signal, and the most natural normalization is to subtract the mean value and divide by the SD. This also guarantees that FD-L does not depend on the absolute variability.

#### SymDyn.

In this method, the original time series is replaced with a new coarser time series that ignores a lot of information but preserves the most essential dynamic features. The time series is converted to a sequence of a few symbols (18, 24, 25). In this study we generated four symbols (*A*, *B*, *C*, *D*) by using the mean and SD of the signal to create value bands for each symbol. The value bands for *A*, *B*, *C*, and *D* were the following
The time series were converted to the symbol sequence. In the next phase, the symbol sequence was interpreted as a sequence of words. Each word consisted of three symbols. Because we used three-symbol words, the number of different words is 4 × 4 × 4 = 64. Each word corresponds to a certain functional pattern in the original time series. Each word has a different probability of appearing in the time series, because the heart rate and blood pressure dynamics favor some patterns and other patterns are rare. To study the word distribution, we have to first index the words, for example, the word could be AAA, BAA, CAA, DAA, ABA, BBA, CBA, DBA, ACA, … CDD, DDD. The distribution pattern can yield interesting information on the system dynamics. However, the word distribution patterns are frequently difficult to interpret because of the lack of data in the literature. Therefore, we measured the complexity of the distribution by using the Shannon entropy
where *p _{i}
* is the probability of each word and

*n*is the number of different words. It is a common practice to normalize the entropy to the unit bits/word. Another measurement to characterize the word distribution is the percentage of the forbidden words, i.e., those whose probability is lower than 0.001.

#### Lempel-Ziv entropy.

In Kolmogorov entropy, the complexity of the sequence is equal to the length in bits of the shortest program or algorithm producing the sequence. If the program is short, the complexity of signal is low, and vice versa. This approach is frequently impossible to calculate without special restrictions.

The symbol sequence consisting of only two different symbols actually presents a binary sequence, and we can use the symbols “0” and “1”. This kind of sequence can always be constructed with two operations: the insertion of the symbol and copying of the subqueue. The Lempel-Ziv complexity can be determined using this information sequence. For example, if we have the infinite sequence “0000…”, the sequence can be constructed by inserting symbol “0” and then copying it and the complexity is 2. The sequence “010101…” is constructed by inserting “0” and “1” and then copying that subqueue. The complexity in this case is 3. In general, we can determine in this way the complexity of any arbitrary sequence (16, 27).

The LZEn of the totally random sequence of length *n*consisting of two different symbols with equal probabilities is
If we divide the complexity of the sequence by the complexity*b*(*n*) of the random sequence, we get the normalized LZEn, which does not depend on the length of the sequence when *n* is large, and 0 ≤ LZEn ≤ 1. The binary sequence is simple to construct: the data values below the mean have the symbol “0” and the values above the mean have the symbol “1”. LZEn can be calculated for arbitrarily short sequences. However, the normalized LZEn is independent of the length of the sequence only when the sequence is longer than 800.

#### Stationarity.

If the operation point of the system producing the signal is constant, the signal can be regarded as stationary. However, the assumptions on the state of operation point are always vague in analyzing heart rate and blood pressure time series data. In practice, the signals are regarded as stationary when there are no obvious changes in the baseline, amplitude distribution, spectral density, or autocorrelation functions do not depend on time. Each of the conditions determines the stationarity in different terms. The applications determine which of the conditions are the most relevant.

We tested the stationarity of signals by checking whether there are any large changes on the baseline of the signal (18). In this method, the signal was divided into subepochs where the mean value was calculated. The ratio of the standard deviation of subepoch means and the standard deviation of the whole time series is a measure of stationarity, StatAv. The length of the subepochs, of course, affects StatAv. We found empirically that 20-data point subepochs were suitable for RR interval and SAP time series data analyses.

## Footnotes

This work was supported by the Academy of Finland.

Address for reprint requests and other correspondence: T. Kuusela, Dept. of Physics, Univ. of Turku, Vesilinnantie 5, FIN-20014 Turku, Finland (E-mail: tom.kuusela{at}utu.fi).

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

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.10.1152/ajpheart.00559.2001

- Copyright © 2002 the American Physiological Society