The Prognostic Value of Non-Linear Analysis of Heart Rate Variability in Patients with Congestive Heart Failure—A Pilot Study of Multiscale Entropy

Aims The influences of nonstationarity and nonlinearity on heart rate time series can be mathematically qualified or quantified by multiscale entropy (MSE). The aim of this study is to investigate the prognostic value of parameters derived from MSE in the patients with systolic heart failure. Methods and Results Patients with systolic heart failure were enrolled in this study. One month after clinical condition being stable, 24-hour Holter electrocardiogram was recording. MSE as well as other standard parameters of heart rate variability (HRV) and detrended fluctuation analysis (DFA) were assessed. A total of 40 heart failure patients with a mea age of 56±16 years were enrolled and followed-up for 684±441 days. There were 25 patients receiving β-blockers treatment. During follow-up period, 6 patients died or received urgent heart transplantation. The short-term exponent of DFA and the slope of MSE between scale 1 to 5 were significantly different between patients with or without β-blockers (p = 0.014 and p = 0.028). Only the area under the MSE curve for scale 6 to 20 (Area6–20) showed the strongest predictive power between survival (n = 34) and mortality (n = 6) groups among all the parameters. The value of Area6–20 21.2 served as a significant predictor of mortality or heart transplant (p = 0.0014). Conclusion The area under the MSE curve for scale 6 to 20 is not relevant to β-blockers and could further warrant independent risk stratification for the prognosis of CHF patients.


Introduction
Congestive heart failure (CHF) remains to be one of the major cardiovascular disorders in the world [1]. Despite its high expenditure in healthcare budgets [2], the mortality rate of CHF patients can be up to 8 times higher than the agematched control population [3]. The present treatment protocols of CHF patients, such as administrating angiotensin converting enzyme inhibitors (ACE-I) and b blockers, have been proven to lower the mortality and hospital admission rate [4]. Nevertheless, the residual risk for mortality and morbidity of CHF remains high even under such treatment protocols [5,6]. Therefore, further novel prognostic predictor is needed to strengthen the treatment strategy in addition to neurohormonal inhibition therapy.
Conventional linear heart rate variability (HRV) analyses, including frequency and time domain analyses, have been reported as prognostic factors for CHF [7,8]. However, heart rate fluctuations have been recognized as complex behaviors originated from nonlinear processes and often with nonstationary property [9][10][11]. Applying linear algorithms to those seemingly irregular and ''patchy'' patterns of heart rate fluctuations [12] may cause the intrinsic computational errors of the linear algorithms [11,13,14]. Properly use of the analyses based on fractals and chaos theory [15][16][17] to qualify or quantify the characteristics of heart rate time series are suggested to serve as a more reliable index of physiological systems in many clinical studies [10,11,18]. As one of such mathematic methods, multiscale entropy (MSE) analysis has focused specifically on characterizing heterogeneous complexity [19]. Such complex structure is ''breakdown'' (loss of information richness) and points to poor prognosis in CHF patients [19,20]. We hypothesized that MSE could yield a prognostic marker which was not relevant to neurohormonal inhibition therapy in CHF patients. The aims of this study were 1) to evaluate the influences of b-blockers on parameters derived from MSE; 2) to assess the prognostic significance of parameters derived from MSE for CHF patients.

Study Population
Patients with manifestation of exertional dyspnea, leg edema and systolic heart failure (LVEF,45% by echocardiography) at the National Taiwan University Hospital were enrolled after giving their inform consents. Baseline information, including age, sex, etiologies for heart failure,diabetes mellitus, hypertension, dyslipidemia (total cholesterol .220 mg/dl), and cardiovascular medication use (b-blockers, ACE-I, angiotensin-II receptor blockers, and spironolactone) was reviewed in medical records and charts. Patients with renal dysfunction (defined by creatine §2.0 mg/dl) were excluded. One month after clinical condition being stable, standard ambulatory 24-hour electrocardiogram (ECG) recorders were placed on all participants. The ECG signals were sampled at 250 Hz and stored in SD memory card for offline analysis on a microcomputer. Subsequently, these patients were followed up and mortality or heart transplantation will be noted as end-point for follow-up. The Ethics Committee of National Taiwan University Hospital approved the study and all patients provided written informed consent.

Data pre-processing
Each digitalized 24-hour ECG data was annotated by an automated algorithm and the annotated file was then carefully inspected and corrected by technicians for extracting the RR intervals. The ectopic beats (including atrial or ventricular premature beats) were interpolated by its adjacent RR intervals. A four-hour period of RR intervals in daytime (between 9AM-5PM) was selected from each recording to avoid confounding effects on nonlinear or linear analysis caused by different sleep stage or diurnal rhythm [21,22]. Only subjects consisted of more than 80% of qualified normal sinus beats were included for further analysis (the typical RR-interval tracings and the corresponding recurrent plots for survival and mortality groups ( Figure S1) as well as all of the RR-interval data are provided as supplementary materials (Materials S1)).

Time and frequency domain analysis
Standard deviation of normal RR intervals (SDNN) and percentage of absolute differences in normal RR intervals greater than 50 ms (pNN 50 ) were calculated to represent the total variance and vagal modulation of the HR. In addition, the spectrum analysis was carried out in accordance with the recommendations of the European Society of Cardiology and the North American Society of Pacing Electrophysiology [23]. The spectral density of each frequency band-high frequency (HF) (0.15-0.4 Hz), low frequency (LF) (0.04-0.15 Hz), and very low frequency (VLF) (0.003-0.04) were computed by average power spectrum.

Nonlinear methods
Nonlinear analysis enables the researchers to probe the fundamental characteristics of the signals. However, unwanted inferences such as noise and nonstationarity may introduce spurious features to the signals [24,25]. The underlying mechanisms of the irregular and unpredicted behavior of the signals can be misinterpreted and the reliability of the results of analysis can be compromised. Two methods had been chosen for their ability to evaluate the main properties of the signals [14,26].

Detrended fluctuation analysis (DFA)
DFA is a modified root-mean-square analysis used to evaluate the fractal correlation beneath the heart rate fluctuation originated from the interacted regulatory mechanisms. The algorithm has been described in detail elsewhere [14]. To briefly introduce this method, at first, it eliminates the environmental inferences by removing the linear-fitted ''local'' trend over different time scales (''box sizes'') in an integrated time series. Next, the root-meansquare fluctuation of this integrated and detrended time series is calculated. This procedure is repeated over different time scales and then the slope of the curve (a exponent) can be estimated on the log-log plot of fluctuations versus box sizes.
In addition, a crossover phenomenon of a exponent in heart rate dynamics between short (4-11 beats) and long (11-64 beats) time scales has been proposed. The short-term (a 1 ) as well as longterm (a 2 ) fractal correlation exponents were calculated to provide better understanding of the fractal correlation property in physiological system [14].

MSE analysis
Instead of simply using single time scale to estimate the complex pattern (irregularity) of a time series, MSE extended this concept to evaluate the complexity of physiological signals on multiple time scales. It comprises of two steps: 1) coarse-graining the signals into different time scales; 2) quantify the degree of irregularity in each coarse-grained time series using sample entropy (SpEn) [27]. Finally, the entropy is calculated as a function of scale, providing a measure of information richness embedded in different time scales. In addition, it has shown that different features of small and large scales in different groups of subjects may assist the clinical categorization [19] and thus three different parameters were derived from the MSE profile: the summations of quantitative values of scale 1-5 (Area 5 ) or scale 6-20 (Area 6-20 ) which represent the complexity exhibit in short and long time scales, respectively; and the linear-fitted slope of the first 5 scales (Slope 5 ) ( Figure 1). Although MSE was successfully applied in physiological signals [18,19,24], nonstationary artifacts especially trends can compromise the estimation of entropy-based analysis by increasing the standard deviation of the data. Hence, detrending process was used to attenuate the spurious influence caused by nonstationarity [24]. In this study, the empirical mode decomposition (EMD) method was adopted as an adaptive filter to eliminate the oscillations slower than VLF range in the original R-R interval signals [28]. The data subsequently evaluated by the MSE analysis after detrending. This algorithm, instead of removing trend with a priori mathematical formulas such as linear or polynomial functions [14,29], could evaluate the hidden dynamics of heart beat fluctuations better [9,10,29].

Statistical analysis
For the independence of different nominal variables between groups, the chi-square test or Fisher exact test were performed. The continuous variables were represented as mean value 6 SD and the normality of those variables was evaluated by using the Shapiro-Wilk test. Then, the Mann-Whitney U test or Student's t test was applied to the between-group comparison accordingly while the Wilcoxon sign test or Student's paired t test was calculated for the intra-group comparison. The receiver operating characteristic curve (ROC) was constructed by the sensitivity and specificity of the continuous variable in predicting the end-point. Area under the ROC curve (AUCs) gave an estimate of the overall discriminate ability. Furthermore, the most predictive indexes will be selected to seek the optimal cut point within the 30 th to 70 th percentile for all patients in 5 th -percentile step. The maximal hazards ratio and independent correlation of variables with event status (mortality) was determined by Cox regression analysis. Then, Kaplan-Meier event probability curves and log rank analysis of the dichotomized groups were obtained. The statistical significance was set at p,0.05.

Characteristics of Patients
A total of 40 heart failure patients (30 males and 10 females) with a mean age of 56616 years were enrolled and followed-up for 6846441 days. Twenty-five patients received b-blockers (either carvedilol or metoprolol). Carvedilol was titrated from 3.25 mg per day and metoprolol was titrated from 12.5 mg per day to maximal tolerable doses. Six patients died or received heart transplantation during follow-up period of this study. The demographic and clinical data were showed in Table 1. No clinical variable was significantly different between these two groups (with or without b-blocker therapy).

Effect of b-blockers on autonomic activities and its fractal properties
While the conventional HRV measurements showed no significant different between these two groups, DFAa 1 showed significantly higher values (p = 0.014) in patients with b-blocker treatment ( Table 2).

Effect of b-blockers on the dynamical complexity assessed by the MSE analysis
There was no significant difference when any single sample entropy value of scales 1 to 20 was compared between patient groups with or without b-blockers. However, the value of Slope 5 was not only significantly lower (p = 0.028) in patients without bblocker therapy but also exhibit a negative value ( Table 2).

Discussion
Although the effects of b-blockers on HRV indices have been extensively studied [30][31][32], the effects of b-blockers on MSE are still unclear. The present study was the first study to assess the relationship between b-blockers and MSE in CHF patients. The main findings of this preliminary study were that the b-blocker therapy may change the short-term complexity (the slope of MSE between scale 1 to 5). But the most significant predictor of mortality or heart transplantation was the long-term complexity (the area under the MSE curve for scale 6 to 20). Therefore, we found an alternative prognostic predictor of CHF in addition to neurohormonal inhibition therapy by assessing the nonlinear characteristics of the heart rate fluctuations.

Effect of b-blocker therapy on linear and nonlinear properties of HRV
All linear HRV measurement showed no difference between patients with or without b-blocker treatment in the present study. While the conventional HRV measurements showed no significant different between these two groups, nonlinear indices, DFA a1 and the value of Slope 5 , were significantly lower in patients without b-blocker therapy.     Other researchers have proposed that the restoration of autonomic function can be assessed by HRV indices [30,31,33]. The discrepancy between our finding and previous studies could be due to several factors. The first, the b-blockers were titrated according to patients' tolerance in our study. Therefore, the duration and dosage of b-blockers were variable. The second, the ECG was recorded 1 month after clinical condition being stable in our study. However, significant changes of most linear HRV parameters are found after 12 weeks of treatments. The short-term fractal scaling correlation index, DFAa 1 , was markedly higher in b-blocker group. The reversal of DFAa 1 is also reported after administrating b-blocker in patients with CHF [33]. Our data showed similar results.
The new nonlinear method, MSE, allows us to evaluate the information richness in heart beat time series over different time scales. Although the underlying mechanisms responsible for the quantitative feature of the MSE in different time scales are still unclear, previous studies have shown that the complexity decreased significantly during aging and further deteriorated in patients with CHF [19]. We applied three different parameters, Area 5 , Area 6-20 , and Slope 5 , to compare patients with or without b-blocker treatment. The summation of entropy values at different time scales may give the quantitative estimation of information richness over certain time scales. That is, Area 5 and Area 6-20 can probe the complexity structure of the heart rate dynamics in short (e.g., 1 to 5 heart beats) and longer (e.g., 6 to 20 heart beats) time scales, respectively. The Slope 5 also outlined the structure of heart rate dynamics in short time scale [19]. The negative value of Slope 5 indicates random-like patterns in short time scales. Therefore, the significant difference between b-blocker and non b-blocker groups might due to dysfunction of the short-term regulatory mechanisms coincided with the results assessed by DFAa 1 .

Complexity analysis as a prognosis predictor
Applying Fourier-based method to nonlinear or nonstationary signals may result in inaccurate estimations [9,29] and compromise the sensitivity of linear HRV measurements. Recently, DFAa 1 has been proposed to be a better predictor for CHF patients [34]. In the present study, SDNN and DFAa1 failed to predict the prognosis of CHF patients. Makikallio et al. noticed that HRV indexes such as SDNN is less sensitive in CHF patient with NYHA.III [35]. Sixteen patients (40%) were classified as NYHA III-IV in our study and five out of them were in mortality group. Moreover, administration of beta-blockade has shown a reversal effect in either linear and nonlinear parameters such as SDNN, HF, LF and DFAa 1 [30,32,33]. Half of the patients with poor outcome were treated with beta blocker which may potentially influence the parameters. SDNN and DFAa 1 may be, therefore, insensitive to predict their prognosis. Area 5 and Area 6-20 derived from MSE were markedly lower in the mortality group. This phenomenon was in agreement with those found by Costa et al. in MUSIC study [19]. Although the underlying mechanisms were still unclear, this preliminary study provided a new insight for the prognosis of CHF by probing the dynamical complexity on the system level. It could potentially offer an alternative marker for the outcome of CHF in addition to neurohormonal inhibition therapy.

Limitations of study
First, our study had small sample size and no placebo-controlled group. Second, all ECG data were recorded in normal ''freerunning'' conditions with possible confounding factors (e.g., physical activities, different breathing patterns, and so on). The additive or dynamical noises may affect the properties of the signals. Although we examined and interpreted the results of the analysis cautiously and the differences characteristics of the patients were unlikely due to the noise, we did not assess the features or level of noise for more detailed information that may benefit the exploration of the underlying deterministic rules. Finally, some parameters related to possible physiological mechanisms of MSE were not collected, such as baroreflex sensitivity, catecholamine levels, and chemoreflex activities.
In conclusions, the area under the MSE curve for scale 6 to 20 is not relevant to b-blockers and could further warrant independent risk stratification for the prognosis of CHF patients. Figure S1 The 4-hour time tracings of RR intervals of survived patient (A) and patient who died after 1 year (B) and the return map traits of three-dimensional RR time series reconstruction (x-axis for RRn, y-axis for RR n+1 and z-axis for RR n+2 ) of survived patient (C) and expired patients (D). Note that the trait of the map in survival patient was similar to that in expired patient. (TIFF) Materials S1 The RR intervals of each patient was output in a one-column Ascii file and categorized into two groups according to their outcomes after 2.5 years follow up. All of the files were packed into a compressed RAR file. (RAR)