Differentiation of Heart Failure Patients by the Ratio of the Scaling Exponents of Cardiac Interbeat Intervals

Heart failure (HF) is one of the most frequent heart diseases. It is usually characterized with structural and functional cardiac abnormalities followed by dysfunction of autonomic cardiac control. Current methods of heartbeat interval analysis are not capable to differentiate HF patients and some new differentiation of HF patients could be useful in the determination of the direction of their treatment. In this study, we examined potential of the ratio of the short-term and long-term scaling exponents (α1 and α2) to separate HF patients with similar level of reduced cardiac autonomic nervous system control and with no significant difference in age, left ventricular ejection fraction (LVEF) and NYHA class. Thirty-five healthy control subjects and 46 HF patients underwent 20 min of continuous supine resting ECG recording. The interbeat interval time series were analyzed using standardized power spectrum analysis, detrended fluctuation analysis method and standard Poincaré plot (PP) analysis with measures of asymmetry of the PP. Compared with healthy control group, in HF patients linear measures of autonomic cardiac control were statistically significantly reduced (p < 0.05), heart rate asymmetry was preserved (Cup > Cdown, p < 0.01), and long-term scaling exponent α2 was significantly higher. Cluster analysis of the ratio of short- and long-term scaling exponents showed capability of this parameter to separate four clusters of HF patients. Clusters were determined by interplay of presence of short-term and long-term correlations in interbeat intervals. Complementary measure, commonly accepted ratio of the PP descriptors, SD2/SD1, showed tendency toward statistical significance to separate HF patients in obtained clusters. Also, heart rate asymmetry was preserved only in two clusters. Finally, a multiple regression analysis showed that the ratio α1/α2 could be used as an integrated measure of cardiac dynamic with complex physiological background which, besides spectral components as measures of autonomic cardiac control, also involves breathing frequency and mechanical cardiac parameter, left ventricular ejection fraction.

Heart failure (HF) is one of the most frequent heart diseases. It is usually characterized with structural and functional cardiac abnormalities followed by dysfunction of autonomic cardiac control. Current methods of heartbeat interval analysis are not capable to differentiate HF patients and some new differentiation of HF patients could be useful in the determination of the direction of their treatment. In this study, we examined potential of the ratio of the short-term and long-term scaling exponents (α 1 and α 2 ) to separate HF patients with similar level of reduced cardiac autonomic nervous system control and with no significant difference in age, left ventricular ejection fraction (LVEF) and NYHA class. Thirty-five healthy control subjects and 46 HF patients underwent 20 min of continuous supine resting ECG recording. The interbeat interval time series were analyzed using standardized power spectrum analysis, detrended fluctuation analysis method and standard Poincaré plot (PP) analysis with measures of asymmetry of the PP. Compared with healthy control group, in HF patients linear measures of autonomic cardiac control were statistically significantly reduced (p < 0.05), heart rate asymmetry was preserved (C up > C down , p < 0.01), and long-term scaling exponent α 2 was significantly higher. Cluster analysis of the ratio of short-and longterm scaling exponents showed capability of this parameter to separate four clusters of HF patients. Clusters were determined by interplay of presence of short-term and longterm correlations in interbeat intervals. Complementary measure, commonly accepted ratio of the PP descriptors, SD2/SD1, showed tendency toward statistical significance to separate HF patients in obtained clusters. Also, heart rate asymmetry was preserved only in two clusters. Finally, a multiple regression analysis showed that the ratio α 1 /α 2 could be used as an integrated measure of cardiac dynamic with complex physiological background which, besides spectral components as measures of autonomic cardiac control, also involves breathing frequency and mechanical cardiac parameter, left ventricular ejection fraction.
Keywords: heart failure, scaling exponents, detrended fluctuation analysis, Poincaré plot, autonomic cardiac control, asymmetry INTRODUCTION Real-life experience warns us that patients do not have the same clinical response to a single treatment regimen and that traditionally accepted clinical parameters are not good enough predictors of the success of the performed therapy. This is also the case in various areas of the heart failure (HF) treatment. For instance, despite the availability of advanced imaging techniques and strict clinical, echocardiographic and electrocardiographic criteria for the selection of patients with HF and an indication for cardiac resynchronization therapy device implantation, still 20-40% of these patients do not have positive response to this therapy (Dhesi et al., 2017). Therefore, in this field it is necessary to find new ways of separating HF patients, with the intent of defining a group of them who would benefit most from the cardiac resynchronization therapy.
Both intrinsic heart rate and its modulation are primarily determined by alterations in the autonomic tone. The autonomic tone is the general activity rate of the autonomic nervous system (ANS) and it is considered to refer to its long-term activity. It is accepted that measurements of linear time domain and frequency domain variables of heart rate variability (HRV) are simple and practical tools to assess autonomic function (Akselrod et al., 1981; Task Force of the ESC and the NASPE, 1996). However, introduction of methods derived from non-linear dynamics for analyzing HRV have shown new complementary insights into dynamic and structural nature of HRV signals (Eke et al., 2002(Eke et al., , 2012Sassi et al., 2015;Platiša et al., 2016). Since multiple regulatory mechanisms of cardiac control operate across different time scales and have shown scaling behavior, classical linear signal analysis methods are often unsuitable to quantify complex HRV content which is far from simple periodicity. One of the featured groups of non-linear measures in quantifying complex physiological dynamics is a group of fractal measures which are used to assess self-affinity of heartbeat fluctuations over multiple time scales. Fractal organization of healthy sinus rhythm dynamics represents complex output from linear and non-linear processes, usually with non-stationary properties. Long-term scaling properties of interbeat time series was first quantified by the scaling exponent (β) as a slope of the regression line of the power versus frequency relation on loglog graph (Kobayashi and Musha, 1982;Saul et al., 1987). In general, the power law scaling exponent is typically calculated in the frequency domain as the β or in the time domain as the Hurst exponent (H) (Bassingthwaighte et al., 1994;Eke et al., 2000). The technique of detrended fluctuation analysis (DFA), based on modified root mean square analysis of a random walk, was proposed to assess the intrinsic correlation properties of a complex cardiac system where scaling exponent (α) of approximately 1 indicates fractal-like behavior of healthy heart rate dynamics (Peng et al., 1995). The obtained exponent is similar to the Hurst exponent, except that DFA may also be applied to non-stationary signals. With this method, the presence of correlations in the fluctuations of heartbeat intervals can be quantified by short-and long-term scaling exponents (α 1 and α 2 ) in two distinct linear regions that determine range of the shortand long-term correlation properties.
Beside monofractal complexity, multifractal analysis introduced by Ivanov et al. (1999) revealed new informations about physiological complexity of HRV signals. Multifractal complexity arises from a large number of local scaling exponents which physiological background was explained with involvement of coupled cascades of feedback loops in healthy cardiac system. More, Amaral et al. (2001) showed significant impact of parasympathetic control on the multifractal properties on HRV where atropine administration resulted with a marked loss of multifractality. Reduced multifractality also was found in pathological state, in patients with congestive HF (Ivanov et al., 1999). Further, Silva et al. (2014) found that loss of multifractality may indicate an impairment of the left ventricular ejection fraction (LVEF) in patients after acute myocardial infarction. One of the later papers of Ivanov et al. (2009) showed that physiological or physical systems with similar 1/f scaling behavior can differ in various levels of complexity which depend on the nature of control mechanisms. The necessity for developing new methods to detect the network between individual organ systems as well as between network of coupled control mechanisms, in response to changes in (patho)physiological conditions, has been recognized in the proposed interdisciplinary area known as a Network Physiology (Lin et al., 2016;Ivanov et al., 2017). This new field focuses on understanding physiological functions with new theoretical framework and analytic formalism.
Fractal view of physiology has become the basis in understanding and controlling physiological networks where both homeostatic and allometric control mechanisms existed (West, 2009). Homeostatic control has a negative feedback character which is local and rapid while allometric control can take into account long-range interactions in complex phenomena (West, 2010). Hence, in recent modeling the dynamics and control of complex physiological phenomena the fractional calculus is more frequently applied (Bogdan et al., 2013). Furthermore, there have been efforts of reconstructing the network between physiological processes when some physiological signals are not observed yet capturing their fractional behavior through fractional order derivatives (Gupta et al., 2018).
In this study, we analyzed heterogeneity of HF patients through application of various methods of HRV analysis. More, we examined potential of the ratio of short-and long-term scaling exponents to differentiate subgroups of HF patients with different short-term and long-term correlation properties. We also calculated standard linear HRV measures with twofold aim. The first one was to show that our data are in line with previously reported results (their reduction and alterations in HF patients) and the second one was to help us to reveal the physiological background of the obtained results.

Subjects
The group of HF patients comprised 46 patients (nine females) with a mean age of 59 ± 2 years (range 37-78 years). All patients had HF with reduced LVEF, on average below 30%. They all had sinus rhythm without any supraventricular or ventricular rhythm disorders, including supraventricular and ventricular extrasystoles. Most of them were in functional class NYHA II (36 patients) and 10 patients had worse functional capacity (NYHA III). Majority of HF patients in this study were receiving optimal therapy for HF, including β-blockers, angiotensinconverting enzyme (ACE) inhibitors, and aldosterone blockers. An individual therapeutic approach had always been applied to the patients, which means that patients were not exclusively receiving in the guidelines recommended, but maximum tolerated doses of drugs, and rarely, in the presence of certain contraindications, some of previously mentioned group of drugs was left out of therapy. Control healthy group was formed from 35 volunteers (17 females) with a mean age 44 ± 2 years (range 35-59 years). All subjects were non-smokers, without medical history. Ethic Committee of the Faculty of Medicine, Belgrade University approved this study .

Experimental Protocol and Data Acquisition
Experiments were performed in the morning between 7:00 and 10:00 a.m. in a quiet setting at the Pacemaker Center of the Clinical Center of Serbia and at the Laboratory for Biosignals, Institute of Biophysics, Faculty of Medicine, Belgrade University. Twenty minutes electrocardiogram (ECG) and respiratory signal data were recorded with sampling rate 1 kHz using Biopac and AcqKnowledge 3.9.1 software (BIOPAC System, Inc., Santa Barbara, CA, United States). The ECG data were collected using the ECG 100C electrocardiogram amplifier module. The classic 1 channel ECG for the measurement of Lead I based on three electrodes placed on left and right shoulder and the right leg was used. The RSP 100C respiratory pneumogram amplifier module with TSD 201 transducer attached to the belt (adjustable nylon strap) was used to measure abdominal expansion and contraction. Transducer was placed on the abdomen, at the point of minimum circumference (maximum expiration). Subjects were relaxed before measurement, and they were supine and breathed with spontaneous breathing frequency. Interbeat (RR) intervals and interbreath intervals were extracted from recorded signals using the Pick Peaks tool from Origin 6.0 (Microcal Origin, Northampton, MA, United States). Breathing frequency (BF) was obtained as a reciprocal value of the mean interbreath interval.

HRV Analysis in Time Domain
A few standard HRV parameters in time domain were determined from time series of RR intervals with our Matlab program: (1) standard deviation of the RR intervals (SDNN), (2) root mean square difference between successive RR intervals, and (3) the percent of RR intervals which were longer by more than 50 ms than the immediately following RR interval (pNN50) (Task Force of the ESC and the NASPE, 1996).

Frequency Domain Analysis
Heart rate variability analysis in frequency domain was performed using Origin 6.0 (Microcal Origin, Northampton, MA, United States) (Platiša and Gal, 2006;Kapidžić et al., 2014). RR interval series were resampled using the mean RR-value for each subject. Power spectrum densities were obtained using FFT with Hanning window (Microcal Origin, Northampton, MA, United States). Absolute values of spectral components were determined carrying out an integration of the power spectrum in the range of total power (TP, 0.0033-0.4 Hz), very low frequency (VLF, 0.0033-0.04 Hz), low frequency (LF, 0.04-0.15 Hz) and high frequency (HF, 0.15-0.4 Hz) (Task Force of the ESC and the NASPE, 1996). The power of RR variability in VLF range is related to long-term regulation mechanisms related to the thermoregulation, to the renin-angiotensin system and to other humoral factors (Task Force of the ESC and the NASPE, 1996). The physiological interpretation of LF spectral component is still controversial because both sympathetic and parasympathetic contributions are involved in this measure (Task Force of the ESC and the NASPE, 1996; Billman, 2011), while HF spectral power is generally accepted as a marker of parasympathetic activation (Akselrod et al., 1981;Billman, 2011). In order to achieve a normal distribution of the data, natural logarithm of spectral powers in absolute units was taken. Relative units of spectral powers were calculated by dividing each spectral component with the sum of all three spectral powers.

Detrended Fluctuation Analysis (DFA)
A modification of the random walk model analysis has been used to quantify the fractal-like scaling properties of RR interval time series (Peng et al., 1995). The root-mean-square fluctuations of the integrated and linear detrended data F(n) were measured in observation windows of varying sizes n and then plotted against the size of window on a log-log scale. The power-law behavior was quantified as the slope of the linear regression line, log F(n) ∼ α log n. This slope is defined as the scaling exponent α. In this study detrended fluctuation function F(n) was calculated by using algorithm from PhysioNet 1 . The short-term scaling exponent α 1 was calculated over the window size n = (4-16) and the long-term scaling exponent α 2 was calculated over the window size n ≥ 16. Scaling exponents were estimated with standard errors and the coefficient of determination (R 2 ) was calculated in OriginPro 8 (OriginLab Corporation, Northampton, MA, United States). An α < 0.5 characterizes signal with anticorrelations (with stronger anticorrelations when α is closer to 0). If α = 0.5 there are no correlations and signal represents (Gaussian) white noise; if α ≈ 1 represents 1/f noise and if α = 1.5 the signal is random walk (Brownian motion) (Peng et al., 1995;Bashan et al., 2008).

Asymmetry of the Poincaré Plot
A typical HRV Poincaré plot (PP) represents scatter graph of RR i+1 = f (RR i ). The two standard parameters SD1 and SD2, called the PP descriptors, describe distribution of points around two diagonals. It is accepted that SD1 describes instant heartbeat intervals variability and quantifies short-term HRV, while SD2 quantifies long-term HRV. Pearson correlation coefficient of PP, noted as r, measured the association between all pairs (RR i , RR i+1 ) in the time series of RR intervals of one subject. Beside excellent visualization capability and short-and long-term HRV information quantification, the standard PP technique revealed asymmetry as one of unexploited physiological phenomena in resting healthy people (Piskorski and Guzik, 2007). Guzik et al. (2006) extended standard descriptor SD1 (dispersion of the PP across the line of identity, y = x) to two new finer descriptors SD1 up and SD1 down . Shortly, pattern of heart rate during acceleration is different to the pattern of deceleration, i.e., in deceleration contribution of the points above the line of identity (RR i < RR i+1 ) is higher than that of the points below the line (RR i > RR i+1 ). They introduced two variables C up and C down : which determine the relative contribution of SD1 up and SD1 down to SD1. Analysis of the PP of RR intervals was done with our Matlab program (The MathWorks Inc., United States). Relation between C up (C down ) and Bauer's deceleration (acceleration) capacity (Bauer et al., 2006) is explained in the Supplementary Material.

Statistics
Normal distribution of data was examined by the Shapiro-Wilk test (appropriate for a small sample size of up to 50 subjects). We used natural logarithm of the spectral powers to obtain their normalized values. The K-means cluster analysis with Euclidean distance measure was performed for continuous variable -the ratio of the scaling exponents α 1 /α 2 . One-way ANOVA was applied to find significant difference in mean values of each variable or parameter: (1) between control and HF group; (2) between four clusters in HF group with Bonferroni post hoc test; and (3) between control and each cluster of HF with Bonferroni post hoc test. Multiple regression analysis was applied to find which variables and parameters predict the ratio of the scaling exponents in HF patients. Statistical analyses were performed using IBM Statistical Package for the Social Sciences (SPSS) version 19.0. Data are presented as mean ± standard errors. P < 0.05 was used as statistically significant. Table 1 shows clinically relevant data of HF subjects and descriptive statistic results for all linear and non-linear HRV measures in healthy control group and HF patients. It can be observed that, as it was expected, patients with HF had statistically significantly reduced values of linear measures of autonomic cardiac control compared with healthy control group. However, some quantifiers of HRV properties were not statistically different between groups (C up and C down , as well as short-term scaling exponent α 1 ). Also, younger control subjects had higher heart rate than HF patients (p < 0.01).

RESULTS
. Data are presented as mean values ± standard errors. LVEF, left ventricular ejection fraction; SDNN, standard deviation of RR intervals; RMSSD, root mean square of successive differences of RR intervals; pNN50, percentage of consecutive RR intervals that deviate from one another by more than 50 ms. SD1, SD2, and r (Pearson PP), C up , and C down are parameters of Poincaré plot. RR, mean interbeat interval; VLF, very low frequency spectral component; LF, low frequency spectral component; HF, high frequency spectra component; TP, total power of spectral power density; BF, breathing frequency; α 1 , short-term scaling exponent; α 2 , long-term scaling exponent.
By applied cluster analysis we found that the ratio of the shortand the long-term scaling exponents, α 1 /α 2 , was significantly capable to differentiate four clusters of HF subjects. All curves of detrended fluctuation functions F(n) plotted versus segment size n on the log-log graph are presented on Figure 1 for healthy subjects and on Figure 2 for the four clusters of HF patients. Values of estimated scaling exponents for each patient in each cluster of HF group could be found in the Tables 2-5. Table 6 reports results of cluster analysis, with the mean and standard errors of all linear and non-linear HRV indexes for each cluster. Averaged values of short-term scaling exponent (α 1 ), long-term scaling exponent (α 2 ), as well as their ratio for each cluster and control group of healthy subjects are shown on Figure 3. It can be seen that in HF patients with reduced autonomic cardiac control, interplay of different correlation properties in heart rhythm over short-term and long-term scales determines four independent subgroups. In comparison of clusters with control group we found that there is no significant difference in α 1 neither in the ratio α 1 /α 2 between control and 2nd cluster, while α 2 was significantly lower in control group than in the first two clusters (Figure 3). Comparison of heart rate dynamics indices between patients with HF showed that there is no significant difference between time HRV measures neither between absolute values of spectral components. However, we found significant difference between subgroups in relative values of spectral powers as well as in their comparison with control group (Figure 4). Results of statistical analysis could be found in the legend of the Figure 4. Also, there is no significant difference with respect to age, NYHA, LVEF, or breathing frequency between clusters patients ( Table 6).
Poincaré plot descriptors SD1 and SD2 showed that they were not able to separate patients in obtained clusters, but there was a tendency toward statistical significance for their ratio SD2/SD1 (Table 6 and Figure 5). Asymmetry variables C up and C down , determined from the Poincaré plots, indicated that HF patients could be separated only into two subgroups: one with dominant deceleration mechanism and the second one, where there was no statistical difference between deceleration and acceleration pattern of regulatory mechanisms (Figure 6). We also found that there is no statistical difference between control group and each cluster in C up and C down .  Values are mean ± standard error. R 2 , coefficient of determination (the goodness of linear fit). Values are mean ± standard error. R 2 , coefficient of determination.
A multiple regression analysis was performed to determine predictors of the ratio α 1 /α 2 . We found that relative spectral powers (rHF and rVLF), the LVEF, normalized total spectral power, and breathing frequency statistically significantly predicted α 1 /α 2 in HF subjects, F(5,40) = 20.966, p < 0.01, R 2 = 0.724 (Table 7). This result indicates complex physiological background of the ratio of scaling exponents which comprised relative cardiac vagal and central autonomic control, mechanical efficiency of the left ventricle, total variability as well as modulatory effect of the breathing process.

DISCUSSION
In the last few decades DFA method, i.e., short-term and longterm scaling exponents separately, showed a greater prediction potential in several cardiac diseases than any parameter from HRV analyses. In this study, we found that in HF patients short-and long-term correlation properties quantified by scaling exponents could gradually change in the opposite directions. Values are mean ± standard error. R 2 , coefficient of determination. Values are mean ± standard error. R 2 , coefficient of determination.
We introduced new parameter, the ratio of short-and longterm correlation properties of heart interbeat fluctuations, which could be used as an integrative parameter of regulatory cardiac mechanisms in HF patients. This parameter is capable to differentiate four clusters which could not be simply classified by some clinical parameters (LVEF or NYHA) or solely by linear and non-linear measures of autonomic cardiac control. The finding for the first cluster could be a good example of previously recognized interactions between short-term and long-term cardiovascular control mechanisms under specific pathological conditions (Hoyer et al., 2007). Namely, Hoyer et al. (2007) showed that cardiovascular system incorporates dominant short-and long-term control mechanisms which are in optimal adjustment in normal healthy conditions. Under pathological conditions, realized with increased collapse of the short feedback loop, the inverse association between their randomness has been recognized. As a system inherent type of readjustment, short-term randomness increased and long-term randomness decreased. We found in the first cluster of HF patients that the ratio of the scaling exponents had the lowest value; α 1 was significantly reduced and significantly lower than α 2 . The shortterm scaling exponent has shown greater clinical discrimination of various cardiac diseases. A reduced α 1 indicates loss of fractal organization in cardiac interbeat intervals and it is a good predictor of mortality in post-infarction patients (Huikuri et al., 2000;Tapanainen et al., 2002), specific risk factor for cardiac death in the elderly and it has been proposed as a strong predictor for HF patients (Ho et al., 1997;Mäkikallio et al., 2001). In this state of cardiac system with reduced ANS activity, loss of short-term correlations in heartbeat intervals could result from reduced capacity for cardiovascular adaptation, vagal tone (high percentage of rHF) and/or alterations in breathing pattern independent on the breathing frequency. This is an interesting finding because there are plenty of data on factors that affect the breathing frequency in HF, such as the effects of some drugs like sedatives, or LVEF, NYHA, some lung diseases, but we have not found out enough about alterations in breathing patterns so far (Forleo et al., 2015). In our earlier papers we showed that fractal organization of interbeat intervals dynamic could also be altered with some physiological processes like standing, exercise, and recovery (Platiša and Gal, 2008). Additionally, in our study with voluntary breathing at different breathing frequencies we found reduction of α 1 with increase in breathing frequency (Platiša and Gal, 2010). Perakakis et al. (2009) also showed that breathing frequencies may bias evaluation of short-term fractal scaling properties. Compared with control, in patients from this cluster fractal long-term correlations become stronger (α 2 increased). Some previous studies have reported change of fractal scaling behavior with aging, where both scaling exponents increased with aging (Ryan et al., 1994;Lipsitz, 1995;Pikkujämsä et al., 1999). That type of correlations is usually related to the rigidness of regulatory mechanisms with reduced control ability. According Guzik's variables of asymmetry these patients also preserved asymmetry in short HRV, as it is observed in healthy subjects, which indicates that dominant deceleration mechanism had a different pattern from acceleration mechanism/s (Guzik et al., 2006). Relative spectral powers of very low frequency (rVLF) and of high frequency (rHF) regions had similar percentage values and we assumed that regulatory mechanisms in this range of frequencies are equally involved in the state of reduced autonomic cardiac control as the result of readjustment in system control mechanisms.
In the second cluster we did not find statistically significant crossover between scaling exponents (fractal correlations are very similar in both ranges). This cluster is very similar to the control group, but with significantly higher long-term scaling exponent α 2 . However, comparison of distributions of relative spectral powers revealed similar rHF, but with a significant reorganization of rLF and rVLF and with domination of slower regulatory mechanism from domain of very low frequencies in these HF patients. More, asymmetry was not reached, and the pattern of short-term deceleration and acceleration mechanisms was not different. We only may speculate that patients from second cluster were in state with some developed compensatory mechanisms of cardiac control.
unchanged. This type of interactions between short-term and long-term randomness under HF as pathological condition, also may be another type of system inherent readjustment, here in the opposite direction compared with the first cluster. In the literature, it can be found that the increase in α 1 is usually related to vagal withdrawal (Tulppo et al., 2001;Castiglioni et al., 2011) but in this HF states probably with preserved feedback loops of regulatory mechanisms. Tendency to decreased α 2 in the fourth cluster may be related to sinus node dysfunction (Shin et al., 2011) or the β-receptor blockade (Castiglioni et al., 2011). It is believed that complexity of interbeat interval fluctuations is generated solely by ANS activity and that quantified linear and non-linear parameters of HRV are quantifiers of modulatory mechanisms of ANS. However, recently published data revealed that interbeat intervals reflect intrinsic complexity with origin in sinoatrial node cells (Ponard et al., 2007;Yaniv et al., 2013) which usually is integrated in the whole complexity of the heart rhythm. Hence, Yaniv et al. (2013) concluded that HRV is determined by the intrinsic properties of cells in the sinoatrial node and the competing influences of the two branches of the autonomic neural input. With these findings, the importance of assessments of long-term scaling exponent has increased. Hotta et al. (2005) showed that α 2 is a more powerful measure for predicting cardiac morbidity and mortality. FIGURE 4 | The distribution of relative spectral components in clusters determined by the ratio of the scaling exponents. Values are given as mean plus standard errors. Relative VLF spectral component (rVLF) was significantly different between control group and two clusters (2nd and 3rd), p < 0.01. Relative LF spectral component (rLF) was significantly different between control group and three clusters (1st, 2nd, and 3rd), p < 0.01. Relative HF spectral component (rHF) was significantly different between control group and two clusters (1st and 3rd), p < 0.05, and with tendency toward statistical significance with the 4th cluster (p = 0.08). Significant and toward significant differences of comparison between clusters of heart failure patients are for rVLF (1st vs. 2nd, p = 0.07 and 1st vs. 3rd, p = 0.08), for rLF (1st vs. 4th, p < 0.01, and 2nd vs. 4th, p = 0.01), and for rHF (1st vs. 2nd, 1st vs. 3rd, and 1st vs. 4th, p < 0.01).
In the study of Shin et al. (2011) long-term scaling exponent was the only parameter which was capable of discriminating differences in heart rate dynamics between patients with sinus node dysfunction. They concluded that reduced value of α 2 is a robust measure and could be an adjunctive tool for improvement of diagnostic performance in detection of sinus node dysfunction. The question is whether this finding can be applied to our patients, because they were mostly treated with highly selective beta blockers, whereas in the previously mentioned study, nonselective propranolol was used in blocking the beta receptors. Also, in healthy subjects in a much lesser extent, such alterations of scaling exponents are characterized as age-related degradation of integrated physiological regulatory systems (Iyengar et al., 1996). Our patients probably had superposition of both effects which are more pronounced in the fourth cluster. The reason for decrease of α 2 may be related to some other physiological background ( Figure 2D). In the sixth and eighth patient, even with statistically approved linear regression analysis, a new scaling regime which is probably related to some slower regulatory mechanism could be identified. Unfortunately, our time series of 20 min length with approximately 2 10 points was not long enough to detect this regime in all patients of the fourth cluster.
It can be observed that the third and fourth cluster are different with respect to the distribution of asymmetry variables. While there was asymmetry in the third cluster, in the fourth cluster short-term asymmetry has ceased to exist. In order to  clarify all these findings, future studies with longer recordings and/or pharmacological recognition need to be done. We also noticed that even younger control subjects had higher heart rate than HF patients. If it is known that the main pathophysiological feature of HF is the imbalance of ANS with increased sympathetic activity, then it is expected that patients with HF will have faster resting rate than control subjects (Hori and Okamoto, 2012). However, all patients with HF in this study were receiving β-blockers, and some of them were treated with additional antiarrhythmics, such as amiodarone or digoxin. Thus, in all these patients there was iatrogenic decrease of heart rate. Also, in previous study we found that HF patients with sinus rhythm and without ventricular extrasystoles have a significantly lower heart rate compared to HF patients with arrhythmias, either ventricular or supraventricular (Radovanović et al., 2018). Finally, it should be added that some researchers suggest the existence of selectivity of effects of aging on autonomic function in healthy subjects; more precisely, they believe that sympathetic function remains unchanged with increasing age (Parashar et al., 2016).
Results of multiple regression analysis showed that the ratio of the scaling exponents was significantly predicted with the small alterations in five independent measures which determined four clusters of HF patients. Comparison of standardized regression coefficient β values showed that relative spectral components, rHF and rVLF, had strong negative relationship with the ratio of the scaling exponents. Relative strength of relationships of LVEF and lnTP with this ratio was much weaker, while with BF it was the weakest. Recognized heterogeneity in HF patients points to the necessity of introducing new approaches in the analysis of cardiac dynamics which will comprise the interactions with the other coupled physiological systems (Ivanov et al., 2017;Kuhnhold et al., 2017). In this pathological condition, analyses of the heart-brain interactions are of special importance because of recognition and quantification of neuroplasticity changes in the dynamics of the brain stem integrators.

Limitations
The significant limitation of this study is a statistically significant age difference between healthy controls and HF patients. It is known that aging, as well as diseases, is accompanied by significant cardiovascular modifications, both structural and functional, although there are studies that indicate that fractal temporal organization of cardiac dynamics does not break down with healthy aging (Ferrari, 2002;Schmitt et al., 2009). What is sure is that with aging the sympathetic activity increases, the renin-angiotensin-aldosterone system activity decreases, respiratory sinus arrhythmia and HRV are reduced, as well as effectiveness of cardiovascular and cardiopulmonary reflexes (Ferrari, 2002;Voss et al., 2015). Some of these changes are also a characteristic of HF and, therefore, a large problem is separating truly HF dependent alterations from those arising from aging. However, we could not age-match our HF patients and healthy controls, because it was not feasible to select subjects without any medical history, even without any risk factor for cardiovascular disease development, among those who were closer in age to the HF patients. In our study, patients were treated with a combination of medications from groups of drugs recommended for the treatment of HF. The differences among patients in doses of given drugs remain an inevitable limitation of this research, since the prescribed dose of the drug was determined by the patient's comorbidities, which were numerous in our patients.
After this study it is necessary to examine the clinical potential of the ratio of short-and long-term scaling exponents, which was the basis for separating of HF patients. We want to determine whether patients of one of these clusters have greater benefit from the cardiac resynchronization therapy. Results of a multiple regression analysis, more precisely, the fact that the examined ratio is at the center of neural, cardiac and respiratory influences, encourages us that this parameter could contribute to a better selection of candidates for this therapy and also be the parameter on which we will rely on in the course of device programming optimization during follow-up period. In the future, this ratio could potentially be used by a device algorithm that would automatically optimize its function. This will surely be the topic of our next research.

CONCLUSION
Our findings showed that in the integral cardiac control quantified by the ratio of short-term and long-term scaling exponents, beside neural cardiac control, mechanical properties of the heart and the modulatory effect of the breathing frequency are significantly involved. This ratio is able of differentiating four clusters of HF patients in sinus rhythm which do not differ in cardiac autonomic control, age, LVEF and NYHA class.

ETHICS STATEMENT
This study was carried out in accordance with the recommendations of the Declaration of Helsinki. The protocol was approved by the Ethic Committee of the Faculty of Medicine, Belgrade University (Ref. Numb. 29/III-5). All subjects gave written informed consent in accordance with the Declaration of Helsinki.

AUTHOR CONTRIBUTIONS
MP, NR, and SP conceptualized the study. MP and NR designed the experiments, performed the study, wrote the manuscript, and collected and analyzed the data. AK developed program for calculation Poincaré plot descriptors and variables of asymmetry. GM and SP conducted clinical examination and recruitment of patients. All authors discussed, edited, and approved the final and revised version of the manuscript.