Introduction

Fractals are figures and objects that are self-containing in recursive fashion. Fractals can be generated by recursive application of mathematical functions. Any part of the object is a replica of the whole object, but of a smaller size, to infinitesimal levels. The world of fractals is very fascinating. The recursive nature of fractals and how close they resemble natural elements make it an important facet to mathematics and computer graphics. Not only geometrical figures behave as fractals, but also mathematically generated time series behave in a self-similar way. It seems that some processes from the real world (e. g. earthquake and river flooding data, Brownian motions₁) may be regarded as fractal time series. To characterize the entire complex behaviour of a fractal in quantitative terms is a major goal in this area. One of the simplest ways to do so is by fractal dimension (FD) estimation.

For a smooth (i.e., nonfractal) line, an approximate length L(r) is given by the minimum number N of segments of length r needed to cover the line: L(r)=Nr. As r goes to zero, L(r) approaches a finite limit, the length L of the curve. Similarly one can define the area A or the volume V of nonfractal objects as the limit of an integer power law of r:

Figure 1

Where the integer exponent is the Euclidean dimension E of the object.

This definition cannot be used for fractal objects: as r tends to 0, we enter finer and finer details of the fractal and the product Nr ^E may diverge to infinity. However, a real number D exists so that the limit of Nr ^D stays finite. This exponent provides an illustration of what the fractal dimension (FD) is. A commonly used practical way to estimate the fractal dimension is the so-called box counting dimension ₁:

Figure 2

In the frequency domain, fractal time series exhibit power law properties:

Figure 3

Where Psd is the power spectral density, f, and the exponent γ is the so.called power-spectral index ₂,₃. For the values region between FD=1.1 and FD = 1.8 the following relationship between FD and γ is valid:

Figure 4

A time series with random phases and power spectrum following the property (2) may be regarded as a synthetic signal with known fractal dimension. Higuchi ₂ generated synthetic signals with different values of γ(and correspondingly different FD values). This allowed him to propose a reliable method for fractal dimension estimation (FD_H).

A common problem in analysing fractal properties of heart rate variability data emerges from the fact that the HRV signal is not a purely fractal process ₃. Available spectral methods ₄ present (besides other) the drawback that some FD values from real subjects are beyond the linearity region where (2) is valid (see our results below). Besides, the robustness of FD_H in presence of nonfractal components has not been evaluated so far.

It seems plausible to propose estimates of FD based on different approaches. This perhaps can reduce any method-related bias. At the same time, agreement between methods when real data are being evaluated, as well as providing sense-making differences between groups of subjects may further support their use.

The importance of Heart Rate Variability is in providing a non-invasive and reliable means to assess autonomic nervous function. Heart rate is primarily controlled by the balance between parasympathetic and sympathetic nervous system acting on the intrinsic pacemaker discharge frequency of the sinoatrial node. The normal intervals between heart beats, (only a small fraction of the entire ECG signal) contains a great deal of information about intracardiac and extracardiac processes including its active function controlled by the autonomic nervous system (ANS) as before-mentioned, therefore the heart rate variability (HRV) is an objective procedure used in medicine, especially in the domain of disorders of ANS. A robust HRV research finding is that it is the best predictor of sudden death ₄ and for mortality from numerous and diverse medical causes ₄. Its clinical relevance was appreciated when was confirmed that phoetal distress was preceded by alterations in interbeat intervals before any appreciable change occurred in the heart rate itself ₅.

Because HRV can show the status of the autonomic nervous system it can indicate the early degree of deterioration of the ANS in diabetes mellitus ₆, HIV/AIDS, Guillain-Barré syndrome, orthostatic hypotension of the Shy-Drager type in parkinsonism, multiple sclerosis, behavioural disorders, brain injury, as well as in chronic alcoholism, and other possible degenerative neurological conditions ₄ HRV also can be used as a simple tool for monitoring therapeutic effectiveness, and it has become known as a death predictor (7-finlandes}.

Besides the well-documented relationship between HRV parameters and ANS, there are ample evidences about the presence of fractal properties in the HRV signal ₃,₇,₈,₁₀,₁₁,₁₂. Even when the research on this field is not recent, literature about HRV fractal properties sometimes is contradictory. Besides, hitherto we do not have a physiological counterpart for the putative fractal properties changes observed in different groups of patients. It is our opinion that we need reliable ways of confirming and characterizing the fractal properties of HRV. In this paper, we rely on a numerical experiments approach, where different indices are calculated from synthetic signals whose fractal dimension is known. The original data are “corrupted” with the presence of nonfractal components in order to explore the robustness of these indices for FD. Assessment.

In this work we studied three different indices for evaluating the FD in a set of HRV traces:

• The FDH according to the definition provided by Higuchi 2.

• The Long Range Slope (LRS) index, a measure of self similarity obtained from time detrended signals 13,14.

• The Autoregressive Dimensional Index (ARDI), obtained from the application of non-linear non-parametric identification analysis to fractal time series 15.

Even when related to each other, each index is based on a different theoretical grounding.

Prior to the FD estimation from HRV data, numerical experiments were carried out with each of the above-mentioned indices to find:

• The functional relationship between FDH and power law index of synthetic signals whose length was comparable to that of 2-hour HRV data. This was done for confirming that FDH is a reliable FD estimate that can be applied to real HRV signals.

• The other two Ix's were compared to FDH for demonstrating the usefulness in FD estimation.

• The robustness of FD estimates obtained using each method with respect to the presence of non-fractal components in the time series.

We concluded that all the three indices had similar performance for estimating the FD of a time series.

The separate use of each of these indices, as well as their combined use seems to be recommendable for the study of fractal time series with unknown FD. The FD was estimated from HRV traces downloaded from the “Fantasia” database. A set of ten HRV recordings corresponding to 5 young (Y) and 5 elderly (O) clinically healthy male humans was used. We obtained that FD values estimated using each of the three indices, were lower in the group of elderly subjects, when compared to their young counterparts. This is in agreement with data from literature and points to the utility of all the three indices.

Material And Method

Data Sets

HRV data

The database included 10 recordings of intervals between ECG R-waves (I-RR), corresponding to 5 young subjects aged between 21-34 years old (Y), and 5 elderly, 68-81 year old subjects (O). All subjects were rigorously screened and a healthy condition was certified. ASCII files with individual recordings (O1.txt, O2.txt, O5.txt, Y1.txt, Y5.txt) were downloaded from the “Fantasia” database, freely available at “www.physionet.org”. Details about the “Physionet” website, as well as about the possibility to use these data for research purposes are described in Iyengar's manuscript₁₁.

Each trace corresponded to an I-RR signal obtained from 2 hours of continuous ECG recording in a supine position, and contained at least 4000 heartbeat counts. The authors of (₁₁) provide a further description about the data.

Signal Analysis

The following indices were estimated from both simulated HRV traces.

Fractal Dimension Following Higuchi's method FD.

Higuchi ₂ proposed a method to calculate the fractal dimension of self-similar curves in terms of the slope of the straight line that fits the length of the curve versus the time interval (the lag) in a double logarithmic plot. The method consist of considering a finite set of data taken at a regular interval v(1), v(2), ..., v(N). From this series it may be constructed a new series v(m,k), defined as

v(m), v(m+k), v(m+2k),...,v(m+[(N-k)/k].k); with m=1,2,...,k. (3)

Where [] denotes Gauss' notation, that is, the bigger integer, and m,k are integers that indicate the initial time and the time interval respectively. The length of the curve v(m,k) is defined as

Figure 5

Then, the length of the curve for the time interval k is given by mean(L(k)), the average over k sets Lm(k). Finally, if mean(L(k)) is proportional to k^D , then the curve is fractal with fractal dimension D. We call this estimate FD_H

Long Range Slope (LRS)

Autoregressive Dimension Index (ARDI)

Enzmann et al₁₅. firstly introduced this index as a measure for evaluating HRV in patients under haemodialysis treatment. The rationale of the method is to apply a non-linear identification approach to nonstationary data. Nonlinear identification has proven to be adequate for the analysis of short duration (180-500 data points) time series whose dynamical nature is unknown ₁₇. The method allows to separate the deterministic and stochastic components of a stochastic nonlinear system. In particular, most of the known classical chaotic attractors could be reproduced by this method ₁₇. It can be viewed either as an extension of classical linear autoregressive estimation to the nonlinear case or either as an extension of classical chaos theory approach to the case when the nonlinear system is fed by innovation noise ₁₈,₁₉. From all the wealth of information that is obtained by this method we selected a single parameter, namely, the optimal order of the autoregressive model.

For ARDI estimation a recording of duration N = 5000 was divided into 25 non overlapping segments 200 data points long each. To each segment the following non-linear autoregressive model was fit:

In = f (I_n-1, I_n-2...I_n-r) + γ_n (5)

Where I_n-1, I_n-2...I_n-r are the n-th, (n-1)-th... R-R intervals in the sequence.

f is a multivariate non-linear function relating the nth interval to the k preceding intervals in the sequence.{γ_n} corresponds to a random, independent, identically distributed variable. The parameter r is the order of the non-linear autoregressive model. The function F is estimated non-parametrically ₂₀.

According to this method, the estimate of f at an an arbitrary point (Z_t-1, Z_t-2,...Z_t-p) of the state space is obtained as a weighted average of the data.

Figure 6

The bandwidth parameter h determines the weight of each neighbouring point in the phase space. In particular, if his too large we have just averaging, whereas for a too short h noise will be incorporated into the deterministic function. A cross minimal validation error criterion has been used for selecting the bandwidth parameter (17,18,19 , Caceres, los 3 papers, peter folklore). The determination of the optimal order of an autoregressive model is difficult task even for parametric models. The introduction of likelihood criteria is an attempt to penalize the good prediction at expenses of too many parameters. For that purpose a cross validation criterion has also being used (17, 18, 19). For a description of the use of cross validation in kernel; nonlinear autorregression the reader is refered to ₂₁.

After computing the optimal order r for each segment of the whole trace it is possible to compute ARDI as the proportion of r-values equal or higher than 15 corresponds to ARDI.

ARDI = Number of segments with optimal order higher than 15 Total number of segments.

Estimation of fractal dimension via index evaluation

Statistical Analysis

Mean and standard deviation for FD_h, FD_LRS, FD_ARDI were obtained for both the Y and O groups. Conventional analysis under Gaussianity assumption was carried out.

Nonparametric permutation analysis is recommended for small samples with unknown distribution₂₂ For each index, all the 10 individual values from the Y and O traces were allocated according to all possible permutations into 2 groups with 5 individuals each. The difference between the mean values for each randomly generated pair of groups was computed. An empirical distribution for the differences between the means was generated. Accordingly, the probability for the mean difference between Y and O the corresponding null hypothesis (H₀) of no differences was tested.

Results

Numerical experiments.

We first tried to reproduce the numerical experiments performed by Higuchi₂ for the case when the duration of the segment is relatively short (N=5000). We obtained that in the FD values region from 1.1 to 1.9 the values corresponded nicely to the power law index values (γ) in accordance with the theoretically expected relationship (2). Our regression yielded:

FD_H=2.48- 0.46γ(n=87, r=0.97); (8)

Thus we regard this result as supporting the validity of FD_H as an estimate for FD.

Previous results from our laboratory suggested that FD values for HRV data are expected in the 1.7 – 1.9 values region. As the result of the NE carried out, the following relationships were found in that FD-values region:

• FD = 2.13636-1.363636*(LRS) (n=9, r =- 0.96) (9)

• FD =1.917338 +0.0098492*(ARDI) - 0,0007324*(ARDI)^2 + 0.0000156*(ARDI)^3 –1.105*10-7*(ARDI)^4 4.02*10-9*(ARDI)^5 (n= 20, 95 % of the variance explained) (10)

The relation between ARDI and FD is illustrated in fig (1). FD_LRS, and FD_ARDI are obtained from (8) and (9) respectively.

Figure 7

Figure 1: Relationship between FD and ARDI. Each point corresponds to the Fractal Dimension evaluated from a simulated signal whose ARDI had been previously estimated.

The degree of concordance between different FD estimates was high when purely fractal signals were evaluated, thus for FD_ARDI and FD_PLI the correlation coefficient was r = 0.96 (n=10).

Addition of a non-fractal signal to a purely fractal signal was a cause for FD increase with all the three indices. As shown in table 1, adding up to about 30% of a non-fractal component to the total variance increases the estimated FD value in about 5-7%. According to different authors, real HRV data contain about 20% of non-fractal contribution top the overall variance₃. At the same time, it seems plausible to expect that most of the clinical conditions related to changes in the spectral peaks (LF, and HF) imply changes in the range from the “normal” value of 22% to lower values around 18 %₄. Following the results shown in table I, it means that changes in nonfractal spectral peaks may not account for more than a 1% change in FD values estimated by any of these indices. Thus, our results suggest that all the three proposed FD estimates might be pertinent for HRV fractal dimension evaluation.

Figure 8

Table I: Percent change in different FD as results of adding a non fractal component to a purely fractal signal (FD = 1.85)

HRV data: In table II the estimates of FD according to each index are shown for each of the HRV recording analysed.

Figure 9

Table II: FD values for the 10 recordings from the “Fantasia” database. Recording's codes correspond to the original file names as presented at the website.

(In table III, a summary is provided for the statistical processing of the data).

Figure 10

Table III: Means standard deviations and type I error probabilities for the groups O and Y. In the case of permutations the probability P0 that the difference between the means of the both groups is zero has been computed.

According to the three indices considered FD is decreased in the group of elderly subjects. Statically significant differences were obtained using either conventional or nonparametric permutation analysis. These results agree with ample literature evidence mainly using other estimation methods. ₁₁

Discussion

Here three indices for time domain estimation of FD in HRV data were proposed. Higuchi had previously proposed one of these indices. The second index (LRS) is a modification of the original DFA method. This index exhibits good predictive capability and has been used by our group for few years. The third index appeared from our early attempts to characterize HRV by a non-linear identification approach.

As our result revealed, the two indices are functionally related to fractal dimension when purely fractal signals were analysed. The robustness of all three indices was similar. However, we do agree that more specific numerical experiments are needed for a more complete characterization of these indices' robustness.

All the three indices documented in a similar way a result with ample literature report: the decrease with age of FD from HRV recordings₁₁.

A further step in our research might be the possibility to use these estimates for assessing FD from short duration traces, as well as from traces sampled al lower frequencies. Our preliminary results (not shown) point to a preservation of the predictive power of these estimates even for 2-min duration traces, whereas only one of the indices (FD_H) loses its clinical predictability as the sampling frequency is reduced from 1000 to 100 Hz.

Recent experience with the advent of chaos theory shows that we must be extremely cautious when trying to apply results of mathematical theories into real objects. An avalanche of chaos “demystification” followed the early “epidemics” of chaos finding in almost any area of physiology₁₇, and the sad sensation that our understanding of different physiological mechanisms remains obscure had become a reason for pessimism.

In our opinion, the sole use of power spectrum measures for assessing the putative fractality of a time series whose mechanism is unknown may appear misleading. Even with mathematical objects this can take place. One example is the known intermittent Pomeau-Manneville map₂₃. This map is a deterministic non-linear low dimensional processes. However, its power spectrum is undistinguishable from that of a typical fractal with “1/f” noise. The non-linear identification approach (from which ARDI is derived) is capable of distinguishing the both. Thus we may expect that these indices will provide not only new ways for assessing HRV fractal properties, but also for understanding some of HRV underlying dynamics.

The fractal properties of HRV may open avenues for radically new promising views to the heart beat regulation. Concepts as self-organized criticality may bring us to the understanding of physiological regulatory mechanisms not encased in the classical framework of the feedback theory scenario₂₄, ₂₅. However, for these views being productive, we need a simple conviction: that fractal-like processes take place in the human body. Our work has been an attempt in that direction.

Acknowledgements