Multifractal Detrended Fluctuation Analysis of Congestive Heart Failure Disease Based on Constructed Heartbeat Sequence

Heart rate variability (HRV) can be used as a common detection method for congestive heart failure (CHF). Existing researches regarding HRV, including both linear indicators and nonlinear characteristics, are mostly based on the <inline-formula> <tex-math notation="LaTeX">$RR$ </tex-math></inline-formula> intervals of the ECG signal. This article proposed a sequence that can reflect the regulation of sympathetic and parasympathetic nerve on heart rate, and on this basis, conducted multifractal detrended fluctuation analysis (MFDFA). We extracted multifractal features to quantitatively compare the complexity of proposed sequence between the healthy and CHF groups. Results showed that abnormal physiological and pathological conditions due to the weakening of autonomic nerve control can reduce the complexity of the heartbeat signal. Estimate the separation performance of all features, the best discrimination is obtained for the area under the mass index spectrum <inline-formula> <tex-math notation="LaTeX">$S1_\tau $ </tex-math></inline-formula> as providing 100% accuracy in separating the Healthy Young and CHF groups, and 90.93% separation accuracy between the Healthy Elderly and CHF groups. This work provide a good basis for the diagnosis of CHF with a novel perspective.


I. INTRODUCTION
Congestive heart failure (CHF) refers to a pathological state due to the decrease in systolic and diastolic functions of the heart, which in turn cannot meet the needs of systemic tissue metabolism [1]. Despite the great progress in diagnosis and treatment, morbidity and mortality are still high, and the prevalence is rising [2]. Around 26 million people worldwide are affected by CHF [3], of which, the prevalence is 1.5%-2.0% in developed countries, and more than 10% in patients >70 years [4]. CHF is a progressive disease, early diagnosis and intervention are extremely important to delay the progress and improve the quality of patients' life.
Heart rate variability (HRV) can be used as a diagnostic basis for abnormal rhythm or sinus arrhythmias in patients with CHF [5], [6]. HRV conventionally refers to the fluctuation of time intervals between adjacent heartbeats, which is generated by the dual modulation of sinus node by the sympathetic and parasympathetic nervous system (SNS & PNS) [7]- [9]. The increase and decrease of heart rate are respectively related to the excitement of the body's The associate editor coordinating the review of this manuscript and approving it for publication was Mohammad Zia Ur Rahman . sympathetic and parasympathetic nerves. The increase in sympathetic tension is often associated with the occurrence of severe arrhythmia, at which time HRV will reduce to a certain extent. Nonetheless, the increase in parasympathetic tone can improve the threshold of arrhythmia, thereby reducing the occurrence of arrhythmia. At this time, HRV appears as a certain increase. Generally, HRV can monitor the tension balance of the sympathetic-parasympathetic nervous system and is a prognostic indicator for evaluating autonomic nervous function and related diseases.
Extensive researches have been carried out regarding HRV, proposing both linear indicators and nonlinear characteristics [10], [11]. Traditional analysis, linear or time-frequency methods, are based on the Fourier transform (FT), which is effective for the analysis of periodic and stationary signals, but can not well reflect the local time-varying characteristics of the signal. Although Wavelet Transform (WT) can focus on any details of the signal, it is not sensitive enough to catch the exact jump time of the frequency components with small amplitude. In distinguishing deterministic ''random'' signals like heartbeat, nonlinear dynamic parameters are more effective than time-frequency analysis methods. Multifractal analysis proved the chaotic characteristics of biomedical signals like EEG, ECG, etc. [12], [13], and has discovered the loss of physiological complexity with CHF [14]. Multifractal is actually a generalization of the fractal system when a single scale exponent is not enough to describe the dynamics of the object, and should be expressed as a continuous exponential spectrum (singular spectrum f (α) ∼ α). Improved formalisms have been developed around multifractal analysis, such as the wavelet transform modulus maxima (WTMM) method and the multifractal detrended fluctuation analysis (MFDFA) method that does not require the modulus maxima procedure. The comparison between these two methods was concerned by Kantelhardt et al. [15], suggesting that the results of MFDFA are slightly more reliable than those of WTMM. Therefore, hereafter we will extract features on the basis of MFDFA to characterize the complexity of heartbeat sequence.
Previous studies on HRV have been based on the RR intervals of ECG signal. Next, we will propose a sequence that helps reflect the alternation of the increase and decrease of heart rate due to the sympathetic and parasympathetic nerve regulation, which will provide powerful diagnostic basis for predicting the physiological and pathological changes.
In this article, we adopted the MFDFA method to analyze the heartbeat signal of CHF patients, and the specific process is as follows: Based on the regulation of heart rate by the autonomic nerve, we first constructed a sequence, that is, the monotonous increase to monotonous decrease amplitude ratios as derived from RR intervals. Then the complexity of the sequence is quantified by multifractal features, including the width of multifractal spectrum, the area under the mass index spectrum and the multifractal spectrum area. Finally, we compared the separation performance of different features.

II. METHODS AND RESULT A. DATA
There is no obvious limitation in the age of onset of CHF, so both healthy young people and healthy elderly people were selected as the control group. The ECG recordings used herein were from the PhysioNet [16], which is a repository of freely-available medical research data and managed by the MIT Laboratory for Computational Physiology. Recordings of 14 CHF patients (NYHA class 3-4) (aged 22-71) were from the BIDMC Congestive Heart Failure Database [17], which is 20 hours in duration, and sampled at 250 Hz. The original analog electrical signals were made with ambulatory ECG recorders (bandwidth 0.1-40Hz). Recordings of 19 healthy young (aged 21-34) and 15 healthy elderly people (aged 68-85) were from the Fantasia database [18]. They underwent rigorous screening and had ECG taken during 2 hours of continuous supine rest. The signals were digitized at 250 Hz. To ensure accuracy, annotation files of both patients and healthy subjects were detected by an automated detector and additional manual checks.

B. PREPROCESSING
Based on the regulation of heart rate by the autonomic nerve, an amplitude-ratio sequence AR was constructed (Fig. 1). Thereinto, Fig. 1

(a) is a segment of ECG signal where RR(n)
is the time interval between peaks R n and R n+1 ; Fig. 1(b) is the extracted RR(n) sequence; Fig. 1(c) shows the construction process of AR (Eq. 1), where m + i represents the trough position of RR(n), and as the sequence grows, m − i and m + i+1 are the nearest peak and trough, respectively. IA(i) represents the amount of monotonous increase in the heartbeat intervals, indicating the activation of parasympathetic nerve. Correspondingly, DA(i) represents the increased sympathetic regulation. Therefore, AR can reflect the balance of cardiac autonomic nerve. .
(1) Fig. 2 shows the probability density function of (a) the RR interval exhibiting a normal distribution, and (b) the AR sequence presenting a power-law distribution. All subjects in the database show similar results.
Once a certain power-law distribution is found, we can deduce whether this is due to the existence of fractals.

C. METHODS
Based on the MFDFA proposed by Kantelhardt et al. [19], we extracted the multifractal features of AR sequence in i) healthy young people, ii) healthy elderly people and iii) CHF VOLUME 8, 2020 patients, and compared the differences between them. The following are the important steps involved in this method: 1) Calculate the cumulative sum of the de-averaged time series where x ave is the average of non-stationary time series x(i).
2) Use the box of size s to divide Y (i) in order, and discard the parts that cannot be divided, then N s non-overlapping boxes can be obtained. In order to make full use of the series, repeat the above segmentation for inverted Y (i). Ultimately, 2N s boxes will be obtained altogether.
3) Fit the series in each box by least squares and calculate the corresponding variance. For v, v = 1, 2, 3, . . . , N s , and where y v (i) is the fitting polynomial of the series in the v-th box. 4) Calculate the q-th order fluctuation function Here, order q can be any real value except 0. For q = 0, F 0 (s) can be calculated by a logarithmic averaging procedure 5) Vary box size s and repeat the step 2)∼ 4) above. Obviously, F q (s) will increase as s goes up. For long-range correlated series, F q (s) will be a power function of s In general, the value of exponent h(q) may depend on q. For stationary time series, h(2) is identical to the well-known Hurst-exponent H . Thus, h(q) can be called the generalized dimension D q . For q > 0, the box v with large variance F 2 (s, v) may dominate the fluctuation function F q (s). In this case, h(q > 0) basically describes the scaling behavior of boxes with large fluctuations. Accordingly, h(q < 0) represents the scaling behavior of boxes with small fluctuations. While for monofractal time series, h(q) is constant for all values of q. Thereinto, 0.5 < H < 1 indicates a (positively) long-range correlated series; H = 0.5 corresponds to a random walk series; and 0 < H < 0.5 demonstrates that the series is long-range anti-correlated, that is, the future increment is negatively correlated with the past increment, and the series has abrupt jump reversal. 6) Mass index τ (q) can be calculated by the following relationship with h(q) Obviously, for a monofractal time series with unique h(q), τ (q) will grow linearly with respect to q, and for a multifractal time series, τ (q) will depend nonlinearly on q. The greater the curvature of curve τ (q) ∼ q, the more complicated the series is. Singularity strength α and singularity spectrum f (α) can be obtained by Legendre transformation where α represents the fractal dimension of subset series in a box, called the local fractal dimension. Due to the large number of boxes, a sequence of different dimensions can be obtained. And f (α) is used to describe the frequency of the dimension. The width of spectrum α = α max − α min can be used to denote the range of fractal dimensions and can provide the complexity information of the fractal body. A larger α represents a greater complexity.

D. RESULTS
The multifractal features of AR sequence in i) healthy young people, ii) healthy elderly people and iii) CHF patients were extracted based on MFDFA method. Thereinto, the sequence length of AR was uniformly intercepted as N = 1249 in accordance with the minimum of all subjects. Research on the MFDFA of nonstationary time series [19] shows that,  systematic deviations from the scaling behavior in Eq. 7 will occur for s < 10. Additionally, fluctuation function F q (s) will become statistically unreliable for large sizes s > N /4. Thus, in this article, it is an adequate range for s varying from 10 to 200. Moreover, we chose q ∈ [−5, 5] to ensure the effectiveness and stability. In order to select the adequate order of the polynomial fitting function, we applied the MFDFA, varying the order from 1 to 3. Fig. 3 shows the results of generalized dimension D q ∼ q. The curves are just slightly different from each other. In terms of simplicity of calculation, we adopted detrending polynomial fitting functions of order 1 in the MFDFA. Fig. 4 is the results of MFDFA for the healthy and CHF groups. Fig. 4a shows the change of generalized dimension D q , singularity strength α and singularity spectrum f (α) with order q. α decrease monotonically with q, indicating that AR is multifractal characterized in all groups. Fig. 4b gives   FIGURE 5. The multifractal features for one healthy young subject. Plot of (a) the area under the mass index spectrum S1 τ , and (b) the area under the multifractal spectrum S1 f . the mass index curve τ (q) ∼ q. It can be concluded from Eq. 8 that, the linear dependence of τ (q) on q indicates a monofractal characteristic, and that more complex AR will correspond to greater curvature. Fig. 4c exhibits the relation between singularity spectrum f (α) and singularity strength α, where α is a typical quantitative indicator of multifractal degree. From the three subplots showed in Fig. 4, conclusion about the degree of complexity that healthy young > healthy elderly > CHF patients can be obtained.
Next, we extracted four features from Fig. 4 that can characterize the degree of multifractality, and performed quantitative comparison between different groups: 1) The width of multifractal spectrum α (Fig. 4c).
2) Area under the mass index spectrum S1 τ (Fig. 5a), which is enclosed by the curve τ (q) ∼ q (blue solid line) and its end-to-end connection (black dotted line). S1 τ contains the two-dimensional space information of mass index spectrum.
3) The classical multifractal spectrum area S(f ) = f (α)dα. Integration of f (α) and α contains all dimension information of AR sequence. 4) Area under the multifractal spectrum S1 f (Fig. 5b) proposed by Wang et al. [20]. S1 f is enclosed by curve f (α) ∼ α (blue solid line) and its end-to-end connection (black dotted line). S1 f effectively avoids the deviation caused by the movement of curve position, and can better reflect the nonlinear nature of AR sequence.
The higher the nonlinear complexity of AR, the greater the features above. Table 1 lists the statistical results of four features, suggesting that the complexity of AR sequence is more in healthy groups than those with CHF groups, and that, data of young people are more complex than the elderly. Significant difference between the healthy and CHF groups (p < 0.01) can be obtained in all features (Fig. 6), among which, S f (Fig. 6c) can effectively distinguish between healthy young people and healthy elderly people (p < 0.05). VOLUME 8, 2020 FIGURE 6. Comparison of α, S1 τ , S f and S1 f between the healthy young, healthy elderly and CHF groups. Thereinto, ' * * ' represents p < 0.01, and ' * ' represents p < 0.05.
For data with small samples, we adopt the following measures to estimate the separation performance between healthy and CHF groups.
1) The ratio η where µ 1 , µ 2 and σ 1 , σ 2 are respectively the means and standard deviations of the features of the two groups being compared. For a good discrimination, difference between the means is much greater than the standard deviation, which will result in η 1. 2) Statistical distance d 2 , which represents the ''distance'' from the means of the ''boundary'' that minimizes incorrect classification for both groups.
Similarly, the larger the d 2 the better the discrimination provided by the method.
3) The total separation rate TSR which ranges in [0,1]. Thereinto, L denotes the length of the overlap of standard deviations between the healthy and CHF groups. TSR = 1 means a clear separation and that the feature can separate subjects into two nonoverlapping groups. Table 2 shows the performance of α, S1 τ , S f and S1 f in the discrimination between healthy and CHF groups. It can be obtained that S1 τ performs better, and provides 100% accuracy in separating the Young and CHF groups, and 90.93% accuracy in separating the Elderly and CHF groups.

III. DISCUSSION
The multifractal phenomenon in heartbeat dynamics is an inherent characteristic of its control mechanism, which does TABLE 2. Separation Performance of α, S1 τ , S f and S1 f Between the Healthy (Young, Y; Elderly, E) and Patient (CHF, C) Groups.
not simply change with other factors, such as external stimuli, intensity of physical activity, and changes in body posture. In fact, aging and disease can lead to a decrease in the body's autonomic control, which in turn affect the degree of nonlinearity in the electrocardiographic system. However, there is no definite mathematical equation to describe such a complex nonlinear dynamic system, so it is very necessary to simulate its complexity through experiments.
In this article, we first proposed a sequence that is derived from ECG signals and can reflect the regulation of autonomic nerve to heart rate. Next, for subjects with different physiological and pathological conditions, we used MFDFA to analyze complexity for the proposed sequence. The following conclusions can be obtained: 1) AR sequence, which is derived from RR intervals, can reflect the alternating changes of heart rate due to the regulation of sympathetic and parasympathetic nerves. Different from a typical normal probability distribution showed by RR sequence, AR presents a power-law distribution with long tail. The emergence of power-law distribution allows us to speculate whether this is due to the existence of fractals.
2) The multifractal analysis of AR sequence confirms the conclusion that aging and disease are often accompanied by a decrease in HRV. Additionally, quantitative results of the complexity of cardiac dynamics for each group are shown as, healthy young > healthy elderly > CHF patients. The decrease in the nonlinear complexity of the electrocardiographic system may be caused by changes in the self-similar structure of the heart tissue.
3) Features extracted from MFDFA, α, S1 τ , S f and S1 f , are all sensitive to the difference between the healthy and CHF groups. Thereinto, the best discrimination is obtained for S1 τ as providing 100% accuracy in separating the Young and CHF groups, and 90.93% separation accuracy between the Elderly and CHF groups.
Existing studies associated with heart failure have exhibited high discriminating performance. Hu et al. compared the discrimination of different indicators for genders. The AUC (area under the receiver operating characteristic curve) corresponding to heart rate deceleration capacity, acceleration capacity, left ventricular ejection fraction, triangle index and SDNN, are: Male, 0.88, 0.84, 0.98, 0.77, 0.78; Female, 0.97, 0.92, 0.95, 0.80, 0.81, respectively [21]. In flexible analytic wavelet transform framework, Kumar et al. used accumulated entropies to detect CHF based on short-term HRV signals, and achieved an accuracy of greater than 97.71% [22]. Additionally, according to the previous literatures, only two mentioned the use of MFDFA to analyze the complexity of RR sequence in CHF patients [23], [24], and the results showed only differences in the width of the multifractal spectrum α. Therefore, we believe that this article has great advantages and contributions in the research of CHF detection based on MFDFA.
In this article, we provide a good basis for the diagnosis of CHF in different ages. Besides, the presentation of AR sequence is noteworthy as it relates the regulation of cardiac autonomic nerve to the multifractal complexity of body's heartbeat signal, which provides a novel perspective for the detection of CHF disease.