## Abstract

Oscillations of heart rate and blood pressure are related to the activity of the underlying control mechanism. They have been investigated mostly with linear methods in the time and frequency domains. Also, in recent years, many different nonlinear analysis methods have been applied for the evaluation of cardiovascular variability. This review presents the most commonly used nonlinear methods. Physiological understanding is obtained from various results from small animals.

## 1. Introduction

The determination of short-term cardiovascular control by means of studying oscillations at different frequencies in heart rate and blood pressure has become widespread, in human as well as in animal studies. Since the study by Akselrod *et al*. (1981) describing the relationship between spectral components of heart rate variability (HRV) and sympathetic and vagal modulation, a vast number of studies have emerged describing the influence of the autonomic nervous system on HRV (Malliani *et al*. 1991; Eckberg 1997). The power spectral density description of HRV shows two distinct peaks (Aubert & Ramaekers 1999; Mainardi 2009): one in the so-called low-frequency (LF) band (0.04–0.15 Hz in humans) and another in the high-frequency (HF) band (0.15–0.4 Hz in humans). HF fluctuations have been attributed to vagal modulation and LF fluctuations appear to be jointly mediated by sympathetic and vagal influences, probably through the baroreflex mechanism (Aubert & Ramaekers 1999). Similarly, blood pressure (BP) variability (BPV) expresses both a HF and a LF component. The latter is an expression of sympathetic vasomotor tone (Verheyden *et al*. 2007), whereas the former has been suggested to be indirectly influenced by vagal RR modulation (Appel *et al*. 1989; Stauss 2007*a*) and mechanical effects of respiration (Triedman & Saul 1994). Autonomic cardiovascular control remains valid for other species, but frequency limits have to be adapted owing to differences in heart rate in different species (Japundzic *et al*. 1990; Rubini *et al*. 1993; Rohrer *et al*. 1998; Aubert *et al*. 1999; Janssen *et al*. 2000; Just *et al*. 2000; Ramaekers *et al*. 2002*a*,*b*; Aubert *et al*. 2004; Fazan *et al*. 2005; Baudrie *et al*. 2007; Stauss 2007*b*; Laude *et al*. 2008).

A different approach to studying cardiovascular control mechanisms is through the analysis of nonlinear dynamics (Kurths *et al*. 1995; Mansier *et al*. 1996; Voss *et al*. 1997, 2007, 2009; Trzebski *et al*. 1998, 2001; Wagner & Persson 1998; Beckers 2002; Beckers *et al*. 2006*a*,*b*,*c*, 2007; Wessel *et al*. 2007). For example, indices of the 1/*f* slope are used for the long-term analysis of RR intervals (Lombardi 2000), whereas a scaling exponent for short-duration recordings (Peng *et al*. 1995*b*) provides information on the presence or absence of fractal correlation properties. However, their physiological background or relationship to autonomic modulation is still a matter of debate (Dabiré *et al*. 1998; Beckers *et al*. 2006*b*), and no up-to-date standards for nonlinear indices have been published (Task Force 1996). Nonlinear indices seem rather to reflect overall, integrated control of heart rate (Lombardi 2000). Nevertheless, the importance of nonlinear behaviour of cardiovascular control was underlined by Yamamoto & Hughson (1991). These authors have shown that, for persons in the supine (awake) position, the contribution of nonlinear fluctuations to the total was 85.5±4.4 per cent.

Nonlinear control of heart rate (HR) and BP poses potential physiological advantages in the possibility of adapting quickly and more subtly to changes in physiological needs. The basic idea behind nonlinear control mechanisms is that, with a rather small amount of energy, a different endpoint can be reached. Analysis methods derived from nonlinear system dynamics have opened up a new approach for studying and understanding the characteristics of cardiovascular dynamics. Nonlinear analysis methods differ from the conventional methods because they are not designed to assess the magnitude of variability but, rather, the quality, scaling and correlation properties of the signals.

The purpose of this paper is to review some recent developments in nonlinear methodology applied to cardiovascular control mechanisms in small animals (mostly rats and mice). An outline of different methods used to identify nonlinear behaviour and results from the literature will be discussed.

## 2. Non-linear control

The complexity of a system can be estimated by calculating its behaviour in phase space. The dimension of the phase space (correlation dimension, CD; Grassberger & Procaccia 1983; Bogaert *et al*. 2001), the sensitivity to initial conditions (Lyapunov exponent, LE; Rosenstein *et al*. 1993) and the system entropy (approximate entropy, ApEn; Pincus 1991; Beckers *et al*. 2001), together with numerical noise titration (NL) (Barahona & Poon 1996; Poon & Barahona 2001), are some of the indices that describe nonlinear complexity behaviour.

Another nonlinear property of a system is its scaling behaviour. It has been shown that normal spectra of HRV exhibit a 1/*f*^{b} behaviour (Zbilut *et al*. 1989; −*b*=slope of log(amplitude) versus log(*f*)). This 1/*f* behaviour is a property of scaling or fractal behaviour and is affected by various pathological conditions (Butler *et al*. 1993) and ageing (Pikkujamsa *et al*. 1999; Beckers *et al*. 2006*c*). A fractal has the essential characteristic that its statistical properties are self-similar at different time scales. Detrended fluctuation analysis quantifies fractal-like correlation properties of the time series (Peng *et al*. 1995*a*). The short- and long-term scaling exponents (Beckers 2002; Raab *et al*. 2006) can be obtained using the method of detrended fluctuation analysis. The spatial occupancy of the time series can be quantified using the fractal dimension (FD; Peng *et al*. 1995*a*). Although, from a mathematical point of view, these scaling indices are related to each other, all of them investigate a different aspect of scaling property of the time series.

Little attention has been devoted to the physiological significance of these fluctuations. Especially with the increased use of nonlinear dynamics in clinical studies (Goldberger *et al*. 1990; Guzzetti *et al*. 1996, 2000; Voss *et al*. 1998; Yeragani *et al*. 2004; Kim *et al*. 2005; Wu *et al*. 2005; Bar *et al*. 2007), it is important to understand the mechanisms underlying the generation of these fluctuations. However, it is becoming clear that altered nonlinear behaviour can be found in several pathologies (Wu *et al*. 2005; Porta *et al*. 2007*a*,*b*), although the exact contributions or importance of the nonlinear indices in these cases also need to be further investigated.

It has been established that nonlinear cardiovascular fluctuations are generated by the autonomic nervous system (Goldberger *et al*. 1990; Mansier *et al*. 1996; Porta *et al*. 1998; Lombardi 2000; Beckers *et al*. 2001, 2006*a*,*b*,*c*, 2007; Beckers 2002). Suppressing the activity of one branch of the autonomic nervous system using pharmacological blockade should alter the nonlinear characteristics of HR and BP pressure fluctuations in order to compensate for the disturbance of the equilibrated system. Using different types of autonomic blocking agents in small animals allows for the examination of the involvement of the sympathetic and/or vagal nervous system in the generation of nonlinear cardiovascular fluctuations.

The present review in small animals (mostly mice and rats) focuses on: (i) the influence of several autonomic blocking agents on nonlinear HRV and BPV indices and (ii) the physiological response mechanisms contributing to these nonlinear fluctuations.

The influence of different pharmacological interventions in rats will be shown on nonlinear complexity measures: CD; the largest LE; ApEn; NL and nonlinear scaling indices FD; 1/*f* slope; and long- (DFAα2) and short-term (DFAα1) scaling exponents using the detrended fluctuation algorithm. Besides the above-mentioned methods, graphical display methods are also used, the initial example being the Poincaré plot. These nonlinear indices can be compared with standard time- and frequency-domain HRV and BPV parameters.

## 3. Initial rat experiments: graphical methods

Nonlinear dynamics methods for the study of cardiovascular control in rats were studied by Dabiré *et al*. (1998). They used a graphical nonlinear method (Zbilut *et al*. 1998; Marwan *et al*. 2002): a recurrence plot to determine specific nonlinear indices. This method looks for recurrent values in a time series. Some parameters are related to the highest LE. The authors suggested an involvement of the α-adrenoceptor sympathetic nerves in the generation of nonlinear fluctuations, based on observations of the recurrence plot. Nonlinear HRV indices decreased after the administration of atropine and prazosin. Nonlinear BPV indices decreased after the administration of prazosin. No changes were found after β-adrenoceptor blockade. These results were later confirmed by the same group (Mestivier *et al*. 2001) in hypertensive rats. Skinner *et al*. (2000) showed reduced dimensionality of HRV in rats after ketamine anaesthesia in combination with haemorrhage. It was stated that this occurred when the nonlinear dynamics within the autonomic nervous system self-organize to adapt a physiological function. Gonzalez *et al* (2000), using recurrence quantification analysis, described the effects of α-adrenoceptor blockade and showed an increase in nonlinear BPV. Vandeput *et al*. (2007) reached exactly the same conclusion using a new nonlinear analysis technique, NL. Atropine produced a decrease in nonlinear HRV, but not in BPV. In their study, the respiratory component of HR and BP fluctuations was mathematically suppressed (frequencies>0.9 Hz).

## 4. Nonlinear analysis methods of cardiovascular control

The nonlinear characteristics of HRV and BPV can be computed using different methods. These can be divided into two categories: (i) indices that describe the scaling behaviour of the nonlinear system (FD, 1/*f*, DFAα1 and DFAα2) and (ii) indices that describe the complexity of the system (CD, LE, ApEn and NL).

The validation of nonlinear methods is not an easy task. Beckers and colleagues (Beckers *et al*. 2001, 2006*a*; Beckers 2002) created several test signals for each nonlinear method to quantify mathematical outcome. Physiological quantification was used to identify changes in autonomic nervous modulation of HR and BP.

### (a) Scaling

#### (i) Fractal dimension

The basic property of a fractal object is self-similarity or scale invariance: the details of the structure are similar when zooming at different resolutions. In general, it can be estimated from the minimal number of cubes while varying *r*. It is also known as a box-counting dimension (Mandelbrot 1982; Goldberger *et al*. 1990; Mansier *et al*. 1996; Ishizawa *et al*. 2000). For example, when the tachogram is recorded for 5, 50 and 500 min, fast erratic fluctuations seem to vary in a similar manner to the slower fluctuations (self-similarity at different time scales; Grassberger & Procaccia 1983).

FD computation is based on the algorithm of Katz (1988). It describes the planar extent of the time series. The higher the FD, the more self-similar the signal will be.

In his study, Beckers *et al*. (2006*b*) showed a decrease in scaling properties of HR variations after α-adrenoceptor blockade. By contrast, BPV was not affected (Beckers *et al*. 2006*b*). β-adrenoceptor and cholinergic blockade demonstrated no explicit changes.

#### (ii) 1/*f* slope

The 1/*f* slope of the log(power)−log(frequency) plot is obtained from the linear regression from 10^{−4} to 10^{−2} Hz. In case the recording length is insufficiently long to cover the 10^{−4} Hz range, an interpolation has to be started at the lowest frequency possible. The 1/*f* slope is calculated on the spectra obtained by fast Fourier transformation on the whole signal, without windowing. The plots have an uneven density that may give greater weight to data in the higher frequency range. Therefore, a logarithmic interpolation of the log–log plot has to be used, giving a balanced number of points for linear interpolation. A slope of −1 is an indication of scaling behaviour.

The only influence on the 1/*f* slope was in BPV and after β-adrenoceptor blockade a decrease from −1.12±0.19 to −1.69±0.26 (Beckers *et al*. 2006*b*).

#### (iii) Detrended fluctuation analysis

Detrended fluctuation analysis quantifies fractal-like correlation properties of the time series and uncovers short-range and long-range correlations. The root-mean-square fluctuation of the integrated and detrended data is measured within observation windows of various sizes and then plotted against window size on a log–log scale (Pikkujamsa *et al*. 1999). The scaling exponent DFAα indicates the slope of this line, which relates log(fluctuation) to log(window size). Both the short-term (4–11 beats) DFAα1 and the long-term (more than 11 beats) DFAα2 scaling exponents can be calculated. The scaling exponent can be seen as a self-similarity parameter, which is characteristic of a fractal. The values of α around 1 are an indication of scaling behaviour. It was shown from HRV data in rats that preprocessing does not change nonlinear DFA parameters (Gomes *et al*. 2002).

The α-adrenoceptor system has an inhibitory action on the nonlinear scaling properties of BP variations (Beckers *et al*. 2006*b*). The effects of β-adrenoceptor on BPV are minor: DFAα2 was closer to 1. These observations are in line with well-known physiological behaviour: α-adrenoceptor receptors are mostly concentrated in the vascular system and, to a lesser extent, in the heart. The sympathetic nerve catecholamines bind preferentially to the β-adrenoceptors in the heart, whereas, in the blood vessels, they bind preferentially to α-adrenoceptors.

### (b) Complexity

#### (i) Correlation dimension

The CD determines an order of the system, i.e. the number of dimensions needed to model the dynamics of the system under consideration (Grassberger & Procaccia 1983; Bogaert *et al*. 2001; Raab *et al*. 2006). In the presence of nonlinear complex behaviour, an attractor in the phase space characterizes the dynamics of the system and its complexity can be quantified in terms of the properties of the attractor.

The time delay for the reconstruction of the attractor can be calculated for each recording separately with the method of the autocorrelation function (Fojt & Holcik 1998). The embedding dimension varies between 2 and 30 (Guzzetti *et al*. 1996). When a finite value is found for the CD of a time series, correlations are present in the signal. To conclude whether these correlations are linear or nonlinear, a surrogate time series needs to be calculated from the signal (Theiler *et al*. 1992). The surrogate time series has a random phase, but the same power spectrum as the original signal. A significant difference between the CD of the surrogate and the original time series indicates that there are nonlinear correlations present in the HRV signal. The significance level is calculated as the difference between the CD of the surrogate data (CDsurr) and the CD of the original data (CDdata) divided by the standard deviation of the surrogate dataset (SDsurr) (in the case of normal distribution). This is called the S-value (S=[CDsurr−CDdata]/SDsurr). A value of S>2 indicates that the measure reflects nonlinear correlations within the two time series. In the case of S<2, no significant difference is found between the time series and the signal is not chaotic.

Zwiener *et al*. (1996*b*) computed the CD and the LE in rabbits using autonomic blockade. However, the blockade was administered immediately after anaesthesia, which influences cardiovascular control and decreased autonomic modulation (Aubert *et al*. 1999). Despite this, autonomic cholinergic blockade decreased the nonlinear content of HRV even further. Beckers *et al*. (2006*b*) could not find any significant changes for CD in rats.

#### (ii) Lyapunov exponent

The trajectories of chaotic signals in the phase space follow typical patterns. Closely spaced trajectories converge and diverge exponentially, relative to each other. LEs measure the average rate of convergence/divergence of these neighbouring trajectories. For dynamical systems, sensitivity to initial conditions is quantified by the LEs. A positive LE can be considered as a definition of nonlinear complexity, provided the system is known to be deterministic (Rosenstein *et al*. 1993). Larger values of the LE indicate more complex behaviour.

The largest LE can be calculated based on the algorithm of Rosenstein *et al*. (1993), which allows the calculation of this parameter on short datasets.

Zwiener *et al*. (1996*a*), in piglets, reported a decrease in complexity measures of HRV after cholinergic blockade, whereas the LE of the BPV remained unchanged. Beckers *et al*. (2006*b*) showed similar results for the LE in rats after cholinergic blockade.

#### (iii) Approximate entropy

Entropy refers to system randomness, regularity and predictability and allows systems to be classified by rate of information loss or generation. The classification of dynamical systems via entropy stems from the theoretical work of Kolmogorov in the 1950s.

ApEn quantifies the entropy, or regularity, of the system (Beckers *et al*. 2001). ApEn more specifically measures the likelihood that runs of patterns will remain close for subsequent incremental comparisons of the pattern length. The probabilistic nature of this parameter implies that it will be most useful to uncover subtle abnormalities or alterations in long-term data that are otherwise not apparent.

As ApEn is influenced by noise present in cardiovascular data, less biased entropy estimators have been proposed as sample entropy (Richman & Moorman 2000) and corrected conditional entropy (Porta *et al*. 1998), showing an improved accuracy.

In most studies, ApEn is calculated according to the formula of Pincus (1991) with fixed input variables *m*=2 and *r*=15 per cent of the standard deviation of the series, as suggested by Goldberger *et al*. (*m* being the length of compared runs and *r* acts as a filter; Goldberger 1990). Higher values of ApEn indicate a more complex structure in the time series.

ApEn of HRV proved to be most sensitive to all blockade interventions: it increased as well after β-adrenoceptor as after α-adrenoceptor and cholinergic blockade (Beckers *et al*. 2006*b*).

#### (iv) Numerical noise titration

The method of numerical NL is an analytical technique that provides a sufficient and robust numerical test to detect chaos, and it gives a relative measure of chaotic intensity, even in the presence of significant noise contamination (Barahona & Poon 1996; Poon & Barahona 2001; Vandeput *et al*. 2007).

White (or linearly correlated) noise of increasing standard deviation (σ) is added to the data until its nonlinearity goes undetected—within a prescribed level of statistical confidence—by a particular indicator at a limiting value of σ=noise limit (NL) value. As an indicator for the nonlinearity, the Volterra–Wiener nonlinear identification method is used. A detailed version of the numerical NL technique can be found in Vandeput *et al*. (2007).

The condition NL>0 indicates chaos and the value of NL gives an estimate of its relative intensity. Conversely, if NL=0, then it may be inferred either that the series is not chaotic or that the chaotic component is already neutralized by the background noise in the data. Therefore, the condition NL>0 provides a simple sufficient test for chaos.

In a recent study in rats (Vandeput *et al*. 2007), the only influence on NL was in BPV with a significant (*p*<0.01) increase after α-adrenoceptor blockade, while β-adrenoceptor and cholinergic blockade demonstrated no explicit changes. Similar results were found in rats with DFA (Beckers *et al*. 2006*b*), suggesting that both techniques describe the same nonlinear signal properties.

## 5. Concluding remarks

Nonlinear dynamics methods are a powerful addition in the tools for studying cardiovascular oscillations. Although the relationship with the branches of the ANS has been shown (Mansier *et al*. 1996), these methods are still not widely used. The major reasons are as follows: (i) their lack of visual representation compared with the linear methods of frequency analysis, (ii) many methods are still under development and need to be validated, and (iii) often conflicting results are presented. Some of the conflicting findings in the literature may be due to the use of selective or non-selective blocking agents and the use of anaesthetized animals, to specific interventions (hypoxia and haemorrhage) or to specific mathematical techniques. Therefore, there is still a need for new well-controlled experiments, specifically aimed at the different types of receptors.

Numerous methods, derived from nonlinear dynamics, are available for processing physiological data. No single measure seems to be more appropriate than the others for physiological or clinical practice. There is a vast amount of data proving that both branches of the autonomic nervous system can generate, or, at least, are involved in the generation of, nonlinear fluctuations. The many different methods that are used to quantify these oscillations describe different properties of the signal. This means that some indices can be closely linked to sympathetic outflow, while others may be more closely related to vagal outflow. In some situations, linear methods are more appropriate than nonlinear, i.e. for stable stationary systems. In other situations, i.e. which show a broad band spectrum without pronounced peaks indicating dynamics from multivariate sources, nonlinear methods are more appropriate. Both methods are more complementary than opposing each other.

Nonlinear system analysis is still a field under development, therefore one may expect still more parameters to be helpful in physiological research or clinical practice. Nonlinear dynamics methods for HRV and BPV, eventually combined with other diagnostic tests (such as controlled breathing or perturbation or inhibition of specific cardiovascular mechanisms), may provide a wealth of unique information on cardiovascular regulation.

## Acknowledgments

This work was partially funded by ESA-PRODEX grants from the Belgian Federal Office of Scientific Affairs. F.B. (2002–2005) and B.V. are supported by ESA–PRODEX grants from the Belgian Federal Office of Scientific Affairs. J.L. is supported by a bilateral agreement between Belgium and China from the Belgian Federal office of Scientific Affairs.

## Footnotes

One contribution of 15 to a Theme Issue ‘Addressing the complexity of cardiovascular regulation’.

- © 2009 The Royal Society