Towards The Estimation Of The Fractal Dimension Of Heart Rate Variability Data
J Hernández Cáceres, H Foyaca Sibat, R Hong, M Sautié, A Namugowa
Keywords
fractal dimension, heart rate variability
Citation
J Hernández Cáceres, H Foyaca Sibat, R Hong, M Sautié, A Namugowa. Towards The Estimation Of The Fractal Dimension Of Heart Rate Variability Data. The Internet Journal of Cardiovascular Research. 2004 Volume 2 Number 1.
Abstract
With the aim of determining the Fractal Dimension (FD) of Heart Rate Variability (HRV) signals, three indices (Ix) were evaluated using HRV data. The relationship between each (Ix) and FD was established via a Numerical Experiments (NE) approach. Signals with known FD, were evaluated using different Ix. Once the functional relationship between FD and Ix was found, an estimate of the FD, based on a particular Ix can be proposed (FD_{Ix}).
The following Ix's were analysed: A) Fractal Dimension estimated by Higuchi's method (FD_{H}); B) Long Range Slope (LRS), obtained from time detrended data, as well as C) the Autoregressive Dimension Index (ARDI) a measure recently proposed by Enzmann et al. Our results showed that, for fractal signals, a functional dependence between FD and each of the other two Ix's was found in all the cases , explaining more than 96% of the total variance for 1.51.95 FD values region. Adding nonfractal components to a pure fractal signal increased FD_{Ix}. The three indices behaved similarly, i.e. 57 % FD_{Ix} increase in the presence of about 30% of non fractal contribution.
A comparison between FD_{Ix} of old
Introduction
Fractals are figures and objects that are selfcontaining 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 selfsimilar way. It seems that some processes from the real world (e. g. earthquake and river flooding data, Brownian motions_{1}) 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
Where the integer exponent is the Euclidean dimension E of the object.
This definition cannot be used for fractal objects: as
In the frequency domain, fractal time series exhibit power law properties:
Where
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 _{2} 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 _{3}. Available spectral methods _{4} 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 methodrelated bias. At the same time, agreement between methods when real data are being evaluated, as well as providing sensemaking differences between groups of subjects may further support their use.
The importance of Heart Rate Variability is in providing a noninvasive 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 beforementioned, 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 _{4} and for mortality from numerous and diverse medical causes _{4}. 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 _{5}.
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 _{6}, HIV/AIDS, GuillainBarré syndrome, orthostatic hypotension of the ShyDrager type in parkinsonism, multiple sclerosis, behavioural disorders, brain injury, as well as in chronic alcoholism, and other possible degenerative neurological conditions _{4} HRV also can be used as a simple tool for monitoring therapeutic effectiveness, and it has become known as a death predictor (7finlandes}.
Besides the welldocumented relationship between HRV parameters and ANS, there are ample evidences about the presence of fractal properties in the HRV signal _{3},_{7},_{8},_{10},_{11},_{12}. 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 nonlinear nonparametric 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 abovementioned indices to find:

The functional relationship between FDH and power law index of synthetic signals whose length was comparable to that of 2hour 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 nonfractal 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
Simulated data.
For generating a time series with a known FD value, the methodology introduced by Higuchi_{2} was followed. For that purpose, Inverse Fourier Transformed (IFT) was applied to vectors with random phases and spectral densities following (4).
For further details about the procedure for generating fractal time series Higuchi's article_{2} is advised.
By conveniently changing the power law index, a set of simulated signals with FD values ranging from FD = 1.5 to FD = 1.95 was submitted to further analysis.
Non purely fractal synthetic signals.
For estimating the influence of non fractal sources into the different FD, a fractal signal with FD = 1.85 was mixed at different proportions by adding to a non fractal signal whose spectrum was not equal to zero only at one frequency value (F=3000), and equal to zero elsewhere. The contribution of the nonfractal signal was measured in proportions of total variance. It ranged from 0% to 29% (see table I)
HRV data
The database included 10 recordings of intervals between ECG Rwaves (IRR), corresponding to 5 young subjects aged between 2134 years old (Y), and 5 elderly, 6881 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_{11}.
Each trace corresponded to an IRR signal obtained from 2 hours of continuous ECG recording in a supine position, and contained at least 4000 heartbeat counts. The authors of (_{11}) provide a further description about the data.
The following indices were estimated from both simulated HRV traces.
Fractal Dimension Following Higuchi's method FD.
Higuchi _{2} proposed a method to calculate the fractal dimension of selfsimilar 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+[(Nk)/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
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)
LRS is one of the parameters estimated during detrended fluctuation analysis _{14}. The overall idea of the method is to estimate the roughness of a time series that was previously corrected for any local linear trend. In our version the time series to be analysed is divided into boxes of equal length, n. In each box of length n, a least squares line (or polynomial curve of order k) is fit to the data (representing the trend in that box). Next, the time series is detrended by subtracting the local trend in each box. Then the detrended time series is integrated The rootmeansquare fluctuation of this integrated and detrended time series is calculated and denoted as F(n) .
This computation is repeated over all time scales (box sizes), from n = minbox to n = maxbox, to characterize the relationship between F(n) , the average fluctuation, and n, the box size. Typically, F(n) will increase with box size n. A linear relationship on a loglog plot indicates the presence of power law (fractal) scaling. Under such conditions, the fluctuations can be characterized by a scaling exponent, i.e., the slope of the line relating log[F(n) ] to log[n]. It has been shown that in HRV _{16}data instead of a continuous line, we observe two lines with a breakpoint. One line correspond to short duration boxes (e.g. n<5), whereas the second corresponds to boxes with n>7. The slope of the long range segment (LRS) was considered in our computations. There is one difference between our procedure and that of _{13}. In the original version, the signal is integrated prior to be detrended. However, our version proved to have a high predictive power in predicting haemodynamic instability among haemodialysis patients _{14}.
Autoregressive Dimension Index (ARDI)
Enzmann et al_{15}. firstly introduced this index as a measure for evaluating HRV in patients under haemodialysis treatment. The rationale of the method is to apply a nonlinear identification approach to nonstationary data. Nonlinear identification has proven to be adequate for the analysis of short duration (180500 data points) time series whose dynamical nature is unknown _{17}. 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 _{17}. 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 _{18},_{19}. 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 nonlinear autoregressive model was fit:
In = f (I_{n1}, I_{n2}...I_{nr}) + γ_{n} (5)
Where I_{n1}, I_{n2}...I_{nr} are the nth, (n1)th... RR intervals in the sequence.
f is a multivariate nonlinear 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 nonlinear autoregressive model. The function F is estimated nonparametrically _{20}.
According to this method, the estimate of f at an an arbitrary point (Z_{t1}, Z_{t2},...Z_{tp}) of the state space is obtained as a weighted average of the data.
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 _{21}.
After computing the optimal order r for each segment of the whole trace it is possible to compute ARDI as the proportion of rvalues equal or higher than 15 corresponds to ARDI.
ARDI =
Estimation of fractal dimension via index evaluation
For estimating the fractal dimension on the basis of a particular index (FD_{IX}), NE was carried out. FD_{LRS}, and FD_{ARDI} were considered, whereas FD_{H} was taken as the “golden rule”.
At least 10 simulated signals with FD ranging from FD=1.1 to FD=2.0 were submitted to the estimation of each index.
The dependence of FD respect to each Ix was fit to a polynomial
FD =P_{l} (Ix) (7)
Where l is the degree of the polynomial P.
Accordingly, the FD of a given time series was obtained via estimating each index and applying (7). For FD_{PLI} a 3^{rd} degree polynomial was used linear regression was used for estimating FD_{LRS}, and 5th degree polynomial was used for FD_{ARDI}.
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_{22} 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_{0}) of no differences was tested.
Results
Numerical experiments.
We first tried to reproduce the numerical experiments performed by Higuchi_{2} 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 FDvalues region:

FD = 2.136361.363636*(LRS) (n=9, r = 0.96) (9)

FD =1.917338 +0.0098492*(ARDI)  0,0007324*(ARDI)^2 + 0.0000156*(ARDI)^3 –1.105*107*(ARDI)^4 4.02*109*(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
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 nonfractal 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 nonfractal component to the total variance increases the estimated FD value in about 57%. According to different authors, real HRV data contain about 20% of nonfractal contribution top the overall variance_{3}. 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 %_{4}. 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
Figure 9
(In table III, a summary is provided for the statistical processing of the data).
Figure 10
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. _{11}
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 nonlinear 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_{11}.
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 2min 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_{17}, 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 PomeauManneville map_{23}. This map is a deterministic nonlinear low dimensional processes. However, its power spectrum is undistinguishable from that of a typical fractal with “1/f” noise. The nonlinear 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 selforganized criticality may bring us to the understanding of physiological regulatory mechanisms not encased in the classical framework of the feedback theory scenario_{24}, _{25}. However, for these views being productive, we need a simple conviction: that fractallike processes take place in the human body. Our work has been an attempt in that direction.
Acknowledgements
We thank Ms Allison Swanby for discussing the paper as well as for redaction corrections. Especially we thank an unknown referee that managed to be strict and encouraging at the same time.