Analysis of Heart Rate Variability by Applying Nonlinear Methods with Different Approaches for Graphical Representation of Results

—There is an open discussion over nonlinear properties of the Heart Rate Variability (HRV) which takes place in most scientific studies nowadays. The HRV analysis is a non-invasive and effective tool that manages to reflect the autonomic nervous system regulation of the heart. The current study presents the results of HRV analysis based on 24-hour Holter ECG signals of healthy and unhealthy subjects. Analysis of heart intervals is performed with the use of original algorithms and software, developed by the authors, to quantify the irregularity of the heart rate. The main aim is the formation of the parametric estimate of patients' health status, based on mathematical methods that are applied on cardiac physiology. The obtained results show that the analysis of Holter recordings by nonlinear methods may be appropriate for diagnostic, forecast and prevention of the pathological cardiac statuses. Different approaches of graphical representation and visualization of these results are used in order to verify this.


INTRODUCTION
Cardiovascular disease is among the top 10 causes of death worldwide.Ischemic heart disease, stroke and chronic obstructive lung have remained the top major causes of death for the past decade.Latest scientific research has shown that cardiovascular diseases can be mathematically analyzed and predicted by cardiac screening for early detection of possible complications to prevent disease or improve quality of life.The diagnostic parameter defined by the electrocardiograms (ECG) is the heart rate variability (HRV) which is calculated by measuring the time between heartbeats.In 1996, the European Society of Cardiology and North American Society of Pacing and Electrophysiology gave recommendations on the clinical usage of the HRV method for evaluating cardiology disease risks as: myocardial infarction (heart attacks), sudden cardiac death and essential hypertension.In addition, diabetes mellitus is associated with a reduced vagal tone and elevated sympathetic activity [1].HRV is reduced in patients with diabetes mellitus, suggesting dysfunction of cardiac autonomic regulation which has been associated with increased risk for pathological cardiac events [2].The main advantage of HRV signals is that it can be calculated real time in non-invasive manner, while all current biomarkers used in clinical practice are discrete and imply blood sample analysis [3].
Methods of HRV analysis are divided in two groups: linear and non-linear methods.Linear methods can be used in timeor frequency-domain for HRV analysis [4].Time-domain parameters are statistical calculations of consecutive RR time intervals, and they are correlated with each other (SDNN, SDANN, pNN50, etc.).Frequency-domain parameters are based on spectral analysis.They are used to evaluate the contribution on HRV of autonomic nervous system (VLF, LF, HF, HF/LF ratio) [5][6][7].The last years, the interest for nonlinear analysis is constantly increasing, because the HRV measurements are non-linear and non-stationary and a considerable part of information is coded in the dynamics of their fluctuation in different time periods.Through implementation of conventional (linear) mathematical methods part of the important characteristics of signal dynamics are missed.Therefore, development and implementation of new non-linear mathematical methods and knowledge, based on the fractal, multifractal and wavelet theory will allow scientists to identify new reasons for HRV fluctuations.This research is conducted on selected control groups of patients with cardiovascular diseases, forming part of an information database of 300 patients of a cardiological clinic; each with 24www.ijacsa.thesai.orghour Holter recordings.Each record contains information of around 100,000 heartbeats.The first problem is based on the fact that fluctuation of the physiological signals possesses hidden information in the form of self-similarity, scale structure, monofractality and multifractality, through the application of these methods [8][9][10][11][12][13].The fractal, multifractal and wavelet-based multifractal analysis of the fluctuations is useful not only for getting the comprehensive information for physiological signals of patients, but also provide a possibility for foresight, prognosis and prevention of the pathological statuses.The second problem is due to the large volume of research information is extremely important the correct determination the trend of the disease of the patients for relatively large period of time -several years, for instance.Such tests are performed periodically to compare the graphic characteristics of images from clinical studies undertaken as a result of treatment and give an idea of the patient's condition and treatment quality.Nonlinear analysis of cardiology data is a relatively new scientific approach to assess the dynamics of heart activity.The above described problems are solved by setting of the following interesting and important objectives: The key components of the proposed approaches are directed at the analysis of HRV by applying nonlinear methods for graphical representation of results, as new direction of the scientific investigations in parallel with research works with the linear methods.

A. Subjects
In this article two groups of signals are analyzed: RR time series of 16 normal subjects and 16 congestive heart failure (CHF) patients.These signals are consisting of around 100 000 data points, corresponding to 24-hour Holter ECG recordings.

B. Methods
Detrended Fluctuation Analysis (DFA) DFA is a technique for detecting correlations in time series [14].These functions are able to estimate several scaling exponents from the RR time series being analyzed.The scaling exponents characterize short or long-term fluctuations.The main steps of the DFA algorithm [15] are as follows: Step 1: The RR interval time series is integrated using (1): . Where:  y(k) is the k th value of the integrated series;  RR j is the j th inter beat interval;  __ RR is the average inter beat interval over the entire series.
Step 2: The integrated time series is divided into boxes of equal length n.
Step 3: In each box of length n, a least square line is fitted to the RR interval data and y n (k) denote these regression lines.
Step 4: The integrated series y(k) is detrended by subtracting the local trend in each box.The root-mean-square fluctuation of this integrated and detrended series is calculated using (2): .
Where: F(n) is a fluctuation function of box size n.
Step 5: The procedure is repeated for different time scales.The relationship on a double-log graph between fluctuations F(n) and the time scales n can be approximately evaluated by a linear model that provides the scaling exponent α given in (3): .α n ≈ F(n) (3) The parameter α is depended by the correlation properties of the signal.By changing the parameter n can be studied how to change the fluctuations of the signal.Linear behavior of the dependence F (n) is an indicator of the presence of a scalable behavior of the signal.The slope of the straight line is used to determine the value of the parameter α.For uncorrelated signals, the value of this parameter is within the range (0, 0.5), where α > 0.5-it is an indication for the presence of correlation.When α = 1, the signal is 1 / f -noise, while α = 1.5usuallyBrownian motion.
In the case of RR time series, DFA shows typically two ranges of scale invariance, which are quantified by two separate scaling exponents, α 1 and α 2, reflecting the short-term and long-term correlation [14].The short-term fluctuation is characterized by the slope α 1 obtained from the (log n, log F(n)) graph within range 4≤n≥11 and the slope α 2 obtained from the range 12≤n≥64.

Rescaled adjusted range Statistics plot (R/S)
The rescaled range is a statistical measure of the variability of a time series introduced by British hydrologist Harold Hurst [16].The Hurst exponent is one closely associated method with the R/S .This exponent is a measure that has been widely used to evaluate the self-similarity and correlation properties of fractional Brownian noise, the time series produced by a fractional Gaussian process.The self-similarity means that the www.ijacsa.thesai.orgstatistical properties (all moments) of a stochastic process do not change for all aggregation levels.The main properties of self-similar processes include:  Slowly decaying variancethe variance of the sample is decreased more slowly than the reciprocal of the sample size.
 Long-range dependence -the process is called a stationary process with longrange dependence if its autocorrelation function is non-summable.The speed of decay of autocorrelations is more hyperbolic than exponential.
 Hurst effectit expresses the degree of self-similarity.
Based on the Hurst exponent value, the following classifications of time series can be realized:  H=0.5 indicates a random series;  0<H<0.5thedata in the signal are anti-correlated;  0.5<H<1the data in the signal are long-range correlated.
The R/S method for the time series X(n) is asymptotically given by a power law:  R(n) is the range which is the difference between the minimum and maximum accumulated values;  S(n) is the standard deviation estimated from the observed data X(n);  H is the Hurst exponent.
To estimate the Hurst exponent is plotted R(n)/S(n) versus n in log-log axes.The slope of the regression line approximates the Hurst exponent [17].

Wavelet Transform Modulus Maxima
One of the most popular tools in wavelet-based multifractal analysis is the Wavelet Transform Modulus Maxima (WTMM) method.This method is based on wavelet analysis that is called -mathematical microscope‖ due to its ability to maintain good resolution at different scales [18].The WTMM method is a powerful tool for statistical description of non-stationary signals, because the wavelet functions are localized in time and frequency.
Physiological signals, such as RR time series, can be efficiently represented by decomposition at different frequencies.The conventional method for this approach is a Fourier analysis.This analysis works well for stationary time series, but not for non-stationary signals, when the frequency content changes over time.Scalograms is based on the wavelet transform simultaneously to provide both time and frequency information and is important for investigating the signal.The WTMM method is used for analysis the multifractal scaling properties of fractal signals.This method uses continuous wavelet transform to detect singularities of a signal.WTMM is based on Wavelet analysis (continuous wavelet transform, skeleton construction) and Multifractal formalism (partition function calculation, scaling exponential function estimation, multifractal spectrum estimation).The method consists in the following basic steps [19,20].
Step 1: Calculation of the Continuous Wavelet Transform (CWT) In many previous researches, the wavelet decomposition has been used in the medical domain [21][22][23][24][25][26][27] and in other domains [28][29][30].A wavelet is simply a finite energy function with a zero mean value.The wavelet transform is defined by the continuous time correlation between the time series and the particular wavelet of scaling parameter τ and shift parameter α: Where the analyzing wavelet is a zero average function with local support, centered around zero.The family of wavelet vectors is obtained by the translations and dilatations of the -mother‖ wavelet: The wavelet transform has a time frequency resolution which depends on the scale α.
Step 2: Calculation of the local maxima of the modulus of the CWT The modulus maxima (largest wavelet transform coefficients) are found at each scale α as the supreme of the computed wavelet transforms such that: .0 τ α) W(τ(    Step 3: Calculation of the partition function Z(q,α) based on wavelets, where α is the dilatation and q is a scale factor The originality of the WTMM method is in the calculation of the partition function Z(q,α) from the maxima lines.The space-scale partitioning given by the wavelet tiling or skeleton defines the particular partition function: Where α is the dilatation and q is a scale factor.This partition function effectively computes the moments of the absolute values of the wavelet resonance coefficients W(τ,α).
Step 4: Calculation of the decay scaling exponent τ(q) The scaling exponent τ(q) is the Legendre Transform of the multifractal spectrum f(α) for self-similar time series and relates the fractal dimensions to the order q of the partition function Z(q,α).The slope of the double-logarithmic plot allows the computation of the decay scaling exponent.This slope can be obtained using linear regression: If the scaling exponent is everywhere convex, that indicates multifractal behaviour of the signal.In case of monofractal behaviour, the scaling function is a line.
Step 5: Estimation of the spectrum of singularities Multifractal formalism uses multifractal spectrum for the detailed fractal analysis of the signal.The Multifractal spectrum function shows the scope of all fractal measures.The Multifractal spectrum function is calculated from scaling function via Legendre transform by formulas: The spectrum width on multifractality degree is ∆α=α maxα min , this quantity is a measure of the range of fractal exponents in the time series, so if ∆α is large, the signal is multifractal.

Poincaré plot
The Poincaré plot analysis is a graphical nonlinear method to assess the dynamic of HRV [31,32].The method provides summary information as well as detailed beat-to-beat information on the behaviour of the heart.It is a graphical representation of temporal correlations within the RR intervals derived from ECG signal.The Poincaré plot is known as a return maps or scatter plots, where each RR interval from time series RR = {RR1, RR2, …, RRn, RRn+1} is plotted against next RR interval.The Poincaré plot parameters used in this paper are SD1, SD2 and SD1/SD2 ratio.SD2 is defined as the standard deviation of the projection of the Poincaré plot on the line of identify (y=x) and SD1 is the standard deviation of projection of the Poincaré plot on the line perpendicular to the line of identify.These parameters can be defined by (11), (12) and (13) Where:  i = 1, 2, 3, …, n and n is the number of points in the Poincaré plot;  var(d) is the variance of d; . Parameter SD1 has been correlated with high frequency power, while SD2 has been correlated with both low and high frequency powers.The ratio SD1/SD2 is associated with the randomness of the HRV signal.It has been suggested that the ratio SD1/SD2, which is a measure of the randomness in HRV time series, has the strongest association with mortality in adults.

III. RESULTS AND DISCUSSION
The analyzed data for describing two types of signals are combined in two groups: records of 16 CHF patients and 16 normal subjects.The results of the R/S method applied to the studied signals to determine the value of the Hurst exponent are shown in Fig. 3 (a) and Fig. 3 (b).The obtained results show that RR time series are correlated, i.e. they are fractal time series.For normal subject, the value of Hurst exponent is high due to the variation being chaotic, and for CHF patient this value decreases because the RR variation is low.
The wavelet decomposition of RR sig.nals can be used to provide a visual representation of the fractal structure of the investigated signals (Fig. 4 (a), Fig. 4 (b)).The wavelet coefficients are presented in their absolute values and coloured in accordance with colour bar.Dark colours correspond to lower absolute wavelet coefficient values.Light colours indicate higher absolute wavelet coefficient values corresponding to large heartbeat fluctuations.The wavelet analysis uncovers hierarchical scale invariance and reveals a self-similar fractal structure in the healthy subjects and a loss of this fractal structure in the unhealthy (CHF) patients.
On Fig. 5 (a) and Fig. 5 (b) show the variations of scaling exponents τ(q) over all the values of q for the normal subject and CHF patient.The constantly changing curvature of the τ(q) curves for the normal subject suggests multifractal behaviour.In contrast, τ(q) is a straight line for CHF patient, indicating monofractal behaviour.Fig. 6 (a) illustrates the multifractal spectrum of RR time series of normal subject and Fig. 6 (b) for CHF patient.The curve of normal subject has multifractal behaviour due to the wide range of local values of the Hölder exponent α (∆α=α maxα min ).The range of values of the Hölder exponent α for RR time series of healthy subject is ∆α=0.79177,and for RR time series of CHF patient is ∆α=0.30351.The interval ∆α, corresponding to the RR signal for CHF patient is more than two times lower than the signal for a normal subject.
The results of the Poincaré plot analysis of RR time series of healthy subject and CHF patient are presented in Fig. 7 (a) and Fig. 7 (b).The Poincaré plot for healthy subject is a cloud of points in the shape of an ellipse (‗comet' shaped plot).On the other hand, points for CHF patient are a cloud of points in the shape of a circle (‗complex' shaped plot).The geometry of these plots has been shown to distinguish between healthy and unhealthy subjects.The obtained Poincaré plot parameters are directly related to the physiology of the heart.The parameter SD1 is the length of the semi-minor axis of the ellipse and it is www.ijacsa.thesai.org the reflection of short term variability of heart rate.The parameter SD2, the length of the semi-major axis, is the measure of long term variability.Since the values of the SD1 and SD2 parameters depend on the RR intervals, the ratio of SD1 to SD2 is used to make comparison among Poincaré plots from different subjects.The values of SD1 and SD2 are higher in normal subjects comparatively to the congestive heart failure subjects.
The modern measuring devices (ECG Holters) and computational methods for data analysis allow the researcher and medical expert derive various summary parameters, analyze recurring patterns, transform the data and much more [33][34][35][36][37].
The values of the investigated parameters (mean ± standard deviation) are reported in Table 1.The described methods are realized by the Matlab software developed in the research project for nonlinear analysis of ECG signals.
There are some limitations in the using of the described analytical system.First, in the study are include restricted number of investigated subjects for analysis of RR time series (16 normal subjects and 16 CHF patients).Secondly, the study is oriented only at nonlinear methods for analysis of RR time series.These limitations are included because this paper is prepared not for showing the medical treatment results, but for demonstration the conceptual developments of the software environment and database as a new and modern tools to analyze, quantify and display the results of nonlinear analysis of HRV signals.
Concerning practical application of the above described nonlinear methods of graphical representation, the graphical user interface of the software environment, database and appropriate signal processing algorithms were developed and implemented in MATLAB 6.0.The database includes background parameters, for example subject's gender, age and level of physical activity as well as disease history.The variations are presented by measurement the time of day and describing the circadian rhythm of the subject.The graphical representation of the circadian rhythm for every subject is stored in database, which enable evaluation of the results in a relatively long period of time for the post-analysis of the prescribed treatment, for instance.Moreover, the circadian rhythm of the physiological states can be assessed with this database's results.

IV. CONCLUSION
The paper presents a part of the latest scientific investigations and the results of application of nonlinear mathematical methods, based on the fractal theory for selected group of patients.The used nonlinear graphical methods for HRV analysis are an effective tool to visualize HRV fluctuations.The results of the cardiac data and the selected methods of analysis managed to give detailed information about the physiological status of the patients and to develop a work towards a framework for prognosis and prevention of pathology status in case of cardiovascular disease.This proves the importance of the proposed study.
The use of graphic-oriented methods for analysis and visualization of cardiac diseases facilitates decision-making by health professionals, especially in case of large amounts of information such as the analysis of diurnal heart rhythms.Furthermore, analysis of trends of patients with cardiovascular diseases is more easily run by periodically comparing the plots of about 100.000 entries heart rate, rather than to make comparisons of numerical values or ECG data.Graphic images of the clinical studies for each patient are automatically stored in a database with easy access and the ability to compare the periodic 24-hour recordings.The study of patient data jointly with the application and comparison of graphics obtained from several nonlinear mathematical methods is recommended.

Fig. 1 (
Fig. 1(a) and Fig. 1(b) distribute the RR interval signals of normal subject and CHF patient.The graphs of RR time data are highly nonstationary (the mean, variance and other variation in the time statistical parameters).The fluctuations of heart-beat time series are larger in healthy subject compared to CHF patient.Fig. 2(a) and Fig. 2(b) illustrate the values of scaling exponents and the slope of the line F(n) on double logarithmic plot obtained by using the DFA method for the investigation of signals.The results indicate a significant difference between patients with CHF and healthy controls in short-and long-time scales.Healthy subjects typically show the fractal behavior of heartbeat dynamics while patients with CHF show an alteration in fractal correlation properties.

Fig. 6 .
Fig. 6. (a).Multifractal spectrum for RR series of healthy subject Fig. 6. (b).Multifractal spectrum for RR series of CHF patient

TABLE I .
PARAMETERS FOR HEALTHY SUBJECTS AND CHF PATIENTS