The Value of Heart Rhythm Complexity in Identifying High-Risk Pulmonary Hypertension Patients

Pulmonary hypertension (PH) is a fatal disease—even with state-of-the-art medical treatment. Non-invasive clinical tools for risk stratification are still lacking. The aim of this study was to investigate the clinical utility of heart rhythm complexity in risk stratification for PH patients. We prospectively enrolled 54 PH patients, including 20 high-risk patients (group A; defined as WHO functional class IV or class III with severely compromised hemodynamics), and 34 low-risk patients (group B). Both linear and non-linear heart rate variability (HRV) variables, including detrended fluctuation analysis (DFA) and multiscale entropy (MSE), were analyzed. In linear and non-linear HRV analysis, low frequency and high frequency ratio, DFAα1, MSE slope 5, scale 5, and area 6–20 were significantly lower in group A. Among all HRV variables, MSE scale 5 (AUC: 0.758) had the best predictive power to discriminate the two groups. In multivariable analysis, MSE scale 5 (p = 0.010) was the only significantly predictor of severe PH in all HRV variables. In conclusion, the patients with severe PH had worse heart rhythm complexity. MSE parameters, especially scale 5, can help to identify high-risk PH patients.


Introduction
Pulmonary hypertension (PH) is a progressive, complex, and fatal disease. It involves heterogenous etiologies and different mechanisms [1], and eventually leads to right heart failure. The mortality of PH patients is high even after contemporary treatment [2]; however, timely and intensive management can improve outcomes even in high-risk patients. In addition, the dynamic adjustment of PH medications, based on disease status during followup, also plays an important role in PH management [3][4][5]. Therefore, a useful tool for PH risk stratification is urgently needed to guide PH treatment. Several prognostic factors of PH have been verified, including sex, exercise tolerance, right heart hemodynamics, and functional performance [6][7][8], and they have been applied in different prediction models.
In 2015, the European Society of Cardiology (ESC)/European Respiratory Society (ERS) PH guidelines first proposed a dynamic PH risk assessment tool, including a combination of imaging, biologic, hemodynamic, performance status, and clinical conditions [1]. This tool has shown good survival prediction between different risk groups [9,10]; however, it requires right heart hemodynamic measurements, which are invasive and difficult to apply for continuous monitoring of PH severity in clinical practice. Therefore, in this study, we propose a non-invasive and convenient tool for PH risk assessment derived from heart rate variability (HRV), namely, heart rhythm complexity analysis.
Heart rhythm complexity analyzes the complexity of changes in heart rate using non-linear methods, and it has been shown to have better predictive value for the diagnosis of PH and heart failure outcomes [11][12][13] than traditional HRV linear analysis [14]. In our previous study, we found that heart rhythm complexity was decreased in PH patients, and that it was useful to differentiate PH patients from normal populations [13]. However, whether heart rhythm complexity is useful in the risk stratification of PH patients is unknown. Therefore, we designed this study to investigate the clinical application of heart rhythm complexity in the risk stratification of PH patients.

Patients
We prospectively enrolled 54 Taiwanese patients with PH from a single center, including 35 with pulmonary arterial hypertension (PAH) and 19 with chronic thromboembolic pulmonary hypertension (CTEPH) from May 2012 to April 2018. Based on the ESC guidelines [1], the diagnosis of PH was made when the patient had suspicious clinical symptoms, and with mean pulmonary arterial pressure (mPAP) no less than 25 mmHg in right heart catheterization. The World Health Organization (WHO) recognizes five groups of PH, categorized by etiology or comorbidity. The PAH was in the WHO group 1 and CTEPH was in the WHO group 4. The diagnosis of WHO group 1 was made when the pulmonary artery wedge pressure (PAWP) less than 15 mmHg, and pulmonary vascular resistance (PVR) more than 3 Wood Units, and without the evidence of left heart disease. The diagnosis of WHO group 4 was made when the ventilation-perfusion lung scintigraphy showed filling defects in PH patients with the same hemodynamics criteria in right heart catheterization as in PAH. PAH and CTEPH have similar pathophysiological mechanisms as vascular arteriopathy [15], presenting as elevated pre-capillary vessel pressure and pulmonary vascular resistance [16]. Other types of PH may involve complex disease mechanisms, such as lung disease or heart failure, which may result in patient heterogeneity, and were excluded from this study. Therefore, we only enrolled these two PH subgroups in the present study to avoid the confounding influence of other pathophysiologies.
All patients underwent echocardiography, blood sampling, right heart catheterization, and 24-h ambulatory electrocardiogram Holter recording. A full record of medical history of each patient was documented, including dyslipidemia, diabetes mellitus, hypertension, and coronary artery disease. The prescription of PH specific medication was recorded as well. The diagnosis of PH was confirmed by right heart catheterization, based on the ESC guidelines [1]. Parameters, including mPAP, PAWP, right atrial pressure, cardiac output, and cardiac index, were all recorded during right heart catheterization. Blood sampling was obtained during the right heart catheterization. We tested N-terminal probrain natriuretic peptide (NT-proBNP), hemoglobin, and creatinine. Both echocardiogram and Holter recordings were performed 2 months before or after right heart catheterization. A six-minute-walk-distance (6MWD) test was recorded 3 months before or after right heart catheterization if the patient was tolerable.
The patients were divided into two groups based on PH severity [17]. The highrisk group was defined as (1) WHO functional class IV or (2) WHO functional class III with severely compromised hemodynamics (right atrial pressure: Pra > 15 mmHg or cardiac index < 2.0 L·min −1 ·m 2 ). PH patients in WHO functional class I to II and in WHO functional class III without severely compromised hemodynamics were in the low-risk group [18,19]. There were 20 patients in the high-risk group (group A) and 34 patients in the low-risk group (group B).
This study was approved by the Institutional Review Board of National Taiwan University Hospital (approval numbers NTUH REC No. 201003042R), and all subjects provided written informed consent. All research was performed in accordance with relevant guidelines and regulations. Reporting of this study followed the Strengthening the Reporting of Observational Studies in Epidemiology statement [20].

Echocardiogram
All patients underwent typical transthoracic echocardiography (iE33 x MATRIX Echocardiography System, Philips, Amsterdam, Netherlands). According to the recommendations of the American Society of Echocardiography, tricuspid regurgitation pressure gradient (TRPG) was measured as the peak flow velocity of tricuspid regurgitation (TRV) using a simplified Bernoulli equation: TRPG = 4 × TRV 2 . Left ventricular ejection fraction in M-mode was measured in the parasternal long axis view [21]. The presence of pericardial effusion or not was documented as well.

24-h Holter Recording and Data Processing
All patients received 24-h ambulatory electrocardiogram Holter recording (Zymed DigiTrak Plus 24-Hour Holter Monitor Recorder and DigiTrak XT Holter Recorder 24-Hour, Philips, Amsterdam, Netherlands) and maintained their original daily activity during the examination without specific limitations. The data were automatically processed using an algorithm and then checked by two technicians. The adopted length of RR Intervals for both linear and non-linear HRV analysis was 4-h and the following criteria was met: (1) between 9 a.m. and 6 p.m.; (2) patients were in awake status; and (3) without sudden increases in heart rate of more than 40 bpm within 1 min to avoid the potential influences of sleep and strong physical activities for both linear and nonlinear analysis, while maintaining enough time length for nonlinear analysis. HRV parameters were processed automatically with MATLAB software.
Nonstationarity can significantly compromise the results of complexity analysis especially for the arrhythmias [22]. We identified the QRS complexes by implementing an adaptive threshold, based on the concept of order-statistic filter, which can be effective for wide ranges and variations of heart rate [23]. Then, the detected QRS peaks were visually inspected to avoid automatic misdetections, and the arrhythmic beats, such as atrial premature contractions, and ventricular premature contractions were removed and replaced by the estimated RR using cubic spline interpolation. Only the RR intervals segments with less than 5% removal were used in this study. In addition, to avoid an unwanted effect of external nonstationarity, we used the empirical mode decomposition method to de-trend the RR series for the oscillation longer than 1 h [24].

Linear HRV Analysis
The interpolated normal-to-normal RR intervals were further used to calculate conventional linear HRV based on the recommendations of the North American Society of Pacing Electrophysiology and the European Society of Cardiology [25]. We analyzed time domain and frequency domain parameters. Time domain analysis included mean RR interval (mean RR), standard deviation of RR interval (SDRR), percentage of absolute differences in normal RR intervals greater than 20 ms (pNN 20 ), and percentage of absolute differences in normal RR intervals greater than 50 ms (pNN 50 ), representing autonomic nervous system modulation of heart rhythm. The RR intervals were first linearly interpolated at 4 Hz and fast Fourier transform was carried out on the resampled signals. The summation of the power over a different frequency band, including 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 Hz), were calculated as the frequency domain parameters.

Non-Linear HRV Analysis
For non-linear HRV analysis, we applied the multiscale entropy (MSE) and detrended fluctuation analysis (DFA) to quantify the fractal properties of the signals, such as longterm memory effect and information richness over different time-scales. DFA was used to quantify the correlation property of inter-beat interval dynamics in the time series, while eliminating the external nonstationarity by removing the linear-fitted trends in a different time-scale (box-size) [26]. Initially, the average of the normal-to-normal intervals was removed. The resultant signal was then integrated and then divided into segments of equal samples n. The fluctuation F(n) of the signal in the corresponding time-scale n was calculated by the root-mean-square of the integrated signal after removing the fitted trends in the segments. The procedure was then repeated in a different time-scale from a small box-size (e.g., n = 4) to a large box-size (e.g., n = 100). On a double log graph of F(n) and the corresponding box-size (n), the slope of the line was defined as the α exponent, representing the fractal correlation property of the time series. Both short-range (α1: 4-11 beats) and long-range time scales (α2: 11-64 beats) were calculated [27].
MSE analysis is used to measure the complexity of the finite length time series. Compared to a traditional single scale, entropy estimation only measures the degree of regularity on a single time scale; MSE uses "coarse graining" proceeding multiple time scales and provides information richness over different time-scales as the complexity of the system. To estimate entropy, we calculated sample entropy (SampEn) for each coarse-grained time series, and then plotted this as a function of the scale factor. To quantify the complexity of the heartbeat dynamics, in short and long time scales, we calculated the entropy values of scale 5 (scale 5), the linear-fitted slope of scale 1-5 (slope 5), area under the MSE curve for scale 1-5 (area 1-5), and area under the MSE curve for scale 6-20 (area 6-20) [28].

Statistical Analysis
Continuous variables were expressed as mean ± standard deviation for normally distributed variables, and median (interquartile range, 25th and 75th percentiles) for nonnormally distributed variables. Categorical variables were expressed as absolute and relative frequencies (percentage). Comparisons were made using the independent t-test and the Mann-Whitney U test between two groups. The chi-square test or Fisher's exact test was used to examine differences between proportions. The discrimination abilities of HRV parameters to high-risk PH were assessed using the receiver operating characteristic (ROC) curve analysis. Logistic regression analysis was used to assess associations between variables and high-risk PH. Significant determinants in univariable logistic regression analysis (p < 0.05), including creatinine, PAH group 1, serum creatinine level, plasma NT-proBNP level, mPAP, PVR, LF/HF ratio, DFAα1, slope 1-5, MSE scale 5, and area 6-20, were then tested in multivariable logistic regression analysis with stepwise selection to identify independent factors that could predict high-risk PH. Category-free (continuous) net reclassification improvement (NRI) and integrated discrimination improvement (IDI), were used to evaluate improvements in the accuracy of the prediction after adding a single nonlinear parameter into a model using only linear parameters. NRI is equal to the sum of the increasing probability for survivors and decreasing probability for non-survivors subtracted by the decreasing probability and increasing probability for non-survivors after adopting the updated model. IDI is defined as the average improvement of survival probability for all patients after adopting the updated model [29,30]. The significance of NRI and IDI statistics was based on approximate normal distributions. All statistical analyses were performed using R software 4.0.3 (http://www.r-project.org, accessed on 10 October 2020) and SPSS version 25 for Windows (SPSS Inc., Chicago, IL, USA). The significance level was set at 0.05 (p < 0.05).

Patient Characteristics
The clinical, echocardiographic, and hemodynamic variables of the enrolled patients are listed in Table 1. There were 20 patients in the high-risk group (group A) and 34 patients in the low-risk group (group B). Compared to group B, there were significantly more patients in group A, in WHO group 1, who had pericardial effusion. In addition, group A had higher levels of serum creatinine and NT-proBNP, and higher TRPG than group B. In pulmonary hemodynamic studies, PVR, and mPAP were significantly higher in group A. The PAH specific medication was listed in Table 1.

Predictors of Interest: HRV Analysis
In linear HRV analysis, group A had significantly lower LF/HF ratio compared to group B. Other linear parameters were comparable between the two groups ( Table 2). In non-linear HRV analysis, group A had significantly lower DFAα1, slope 1-5, scale 5 and area 6-20 compared to group B ( Table 2). The entropies over different time scales in group A and group B were shown in Figure 1.

The Effect of Adding Heart Rhythm Complexity to the Linear HRV Parameters to Identify High-Risk PH Patients
In both NRI and IDI models, the MSE scale 5 significantly improved the discrimination power of all linear HRV parameters, including mean RR, SDRR, VLF, LF, HF, and LF/HF ratio. Area 6-20 significantly improved the discrimination power of mean RR, VLF, and HF in both NRI and IDI models, and SDRR, LF, and LF/HF ratio in IDI model. DFAα1 significantly improved the discrimination power of SDRR, VLF, LF, and HF in both the NRI and IDI models, and mean RR in the IDI model (Table 4).

Discussion
The main finding of this study was that heart rhythm complexity was significantly depressed in high-risk PH patients. In addition, adding heart rhythm complexity predictors to traditional linear HRV parameters improved the power to predict high-risk PH patients. This is the first study to demonstrate an association between heart rhythm complexity and severity of PH, and the better performance of heart rhythm complexity in identifying high-risk PH patients than traditional HRV parameters.
PH is a critical disease, which needs an early diagnosis and timely management. Patients classified as being at high risk according to the 2015 ESC/ERS PH guidelines have a worse prognosis compared to patients at low risk. Sitbon et al. demonstrated that poor functional status was associated with poor outcomes. In their study, PH patients in WHO functional class IV and those in class III with severely compromised hemodynamics had the worst outcomes [19]. Previous studies have demonstrated that early interventions including both pharmacological and multidisciplinary team care can improve the outcomes of PH patients, even those with severe disease and poor functional status [5]. Therefore, identifying high-risk patients is essential for the management of PH. Several survival prediction models have been proposed for PH patients; however, they are complex and difficult to use [31]. Currently, the 2015 ESC/ERS PH guidelines advocate assessing the risk of PH by using a combination of several different tools, and this method is widely used in daily practice [1]. However, risk assessment requires invasive right heart catheterization, which is difficult to apply in frequent monitoring during follow-up. Therefore, there is still a strong unmet need for an easy-to-use tool to allow for both timely and continuous monitoring of disease status to improve the clinical care of PH patients.
HRV is a useful non-invasive tool, which has been studied in many diseases, including coronary artery disease, heart failure, and even pulmonary hypertension [32][33][34]. It has been correlated with autonomic dysfunction and used as an outcome predictor. Porte et al. demonstrated that heart rate complexity parameters decreased due to sympathetic activation during postural change [35]. Another study showed that sympathetic activation during senescence was associated with impaired heart rate complexity [36]. These evidences supported that the usefulness of heart rate complexity in monitoring sympathovagal balance. Pulmonary hypertension is characterized by increased pulmonary artery pressure leading to right ventricular failure [37]. The serum norepinephrine increased in patients with PAH similar to those with congestive heart failure as the indicator of cardiac sympathetic activation [38]. Furthermore, sympathetic activation has also been correlated with the severity of PH [39,40]. Several studies also showing that measuring autonomic system regulation using HRV could be a predictor of disease severity and long-term outcomes in PH [41][42][43][44]. Since that, overactivation of sympathetic systems is likely to be one of the major reasons explaining the worse HRV and complexity in severe PH patients. Bienias et al. demonstrated that patients with arterial or chronic thromboembolic PH had significantly impaired heart rate turbulence, a linear HRV parameter [45]. Recently, Peng et al. proposed the heart rhythm complexity derived from two non-linear parameters of HRV, DFA, and MSE, to evaluate complexity change in the biological systems [26,28,46]. Heart rhythm complexity was shown to have better efficacy and predictive power for various diseases than traditional HRV [14,47].
Heart rhythm complexity measures the complexity of changes in the R-R interval, which contains detailed information derived from heart rate dynamics. Once a biological system has become diseased, the complexity breaks down, and non-linear HRV analysis, can detect subtle changes at an early stage [48]. In a retrospective study, abnormal DFAα1 in asymptomatic heart failure patients was associated with the onset of heart failure years in advance of the first clinical event [49,50]. Tsai et al. recently demonstrated that heart rhythm complexity had a better prognostic value for cardiovascular events in patients undergoing peritoneal dialysis compared with linear HRV analysis [47]. In recent years, heart rhythm complexity was extensively studied in many fields, including left heart failure [51], post-infarction myocardial function [52], patients undergoing dialysis [12,47,53], severity of abdominal aorta calcification [54], primary aldosteronism [24], stroke [55], and PH [56]. These studies support the importance of heart rhythm complexity in clinical practice and its potential role in disease risk stratification. In the present study, we demonstrated that heart rhythm complexity parameters, especially MSE scale 5, were significantly associated with PH disease severity and could be used in PH risk stratification. To the best of our knowledge, this is the first study to apply heart rhythm complexity to the prediction of PH disease severity. Although the improvement of the complexity can be attribute to not only the enhanced complexity characteristics but the magnitude of HRV [57], combining different parameters of MSE can give us better information related to the "quality" (complexity) or the "quantity" (magnitude of HRV) changes [58]. Furthermore, model-free complexity can assess the embedded space with variable scales grouped by the K-nearest-neighbor to avoid coarse-graining that may introduce bias due to the fixed dimensions as well as the aliasing filter effect [57,59]. Recently, the local version of the sample entropy was proposed to eliminate the nonstationary effect on the results of complexity analysis [60]. The research by using those new methods warrants further study. In addition, the cardiopulmonary coupling is another important issue in the HRV analysis focusing on the interaction between cardiovascular and cardiorespiratory systems. The cardiorespiratory coupling between the systems is thought to be with each other in a nonlinear way [61]. MSE has been used to evaluate the cardiorespiratory coupling and the asynchrony. Platiša et al. demonstrated that primary alterations in the regularity of cardiac rhythm resulted in changes in the regularity of the respiratory rhythm in heart failure patients [62]. However, there were limited studies investigating the cardiorespiratory coupling in PH patients. Further studies may be needed to integrate the role of cardiorespiratory coupling in PH patients.
Compared with heart rhythm complexity, linear HRV parameters, including SDRR, SDRR index, VLF, LF/HF ratio, and heart rate turbulence have been widely studied to assess PH [45,63]. Recent studies have also demonstrated an association between impaired linear HRV parameter, SDRR, and PH disease severity markers, including impaired WHO functional status, decreased 6MWD, impaired tricuspid annular plane systolic excursion, right ventricular systolic function, higher TRPG, and NT-proBNP level [64][65][66]. In this study, we first demonstrated a better association between heart rhythm complexity and PH disease severity compared to traditional HRV analysis. Second, the discrimination power of linear HRV for PH disease severity improved significantly after combining heart rhythm complexity parameters. The combination of linear and non-linear HRV parameters to form a new predictive model may have further improved its risk stratification ability and outcome prediction.
There are several limitations to this study. First, this is a pilot study. The number of cases was small, and further studies are needed to validate the results. In addition, model-free complexity analysis or entropy with local characteristics can preserve more information instead of a one-fit-all algorithm. Those methods should be included for a large-scale study to probe the underlying pathophysiological mechanisms related to the changes of the complexity of the PH patients. Second, we only enrolled PH patients in WHO group 1 and group 4, and future studies should enroll different groups of PH patients to investigate the potential application of HRV in these patients. Third, this pilot study is a cross-sectional design and lacks clinical long-term follow-up data. A prospective cohort study with clinical end-point follow-up is needed to confirm the utility of heart rhythm complexity on clinical outcome predictions.

Conclusions
This study demonstrated that high-risk PH patients had worse heart rhythm complexity. MSE scale 5 had the best discrimination power to predict high-risk PH patients. Moreover, adding MSE scale 5, area 6-20 or DFAα1 to linear HRV parameters significantly improved the predictive power for high-risk PH patients. Heart rhythm complexity can potentially be used as (i) an indicator of PH disease severity, and (ii) to stratify the risk of PH.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.