Effect of heart rate correction on pre- and post-exercise heart rate variability to predict risk of mortality—an experimental study on the FINCAVAS cohort

The non-linear inverse relationship between RR-intervals and heart rate (HR) contributes significantly to the heart rate variability (HRV) parameters and their performance in mortality prediction. To determine the level of influence HR exerts over HRV parameters' prognostic power, we studied the predictive performance for different HR levels by applying eight correction procedures, multiplying or dividing HRV parameters by the mean RR-interval (RRavg) to the power 0.5–16. Data collected from 1288 patients in The Finnish Cardiovascular Study (FINCAVAS), who satisfied the inclusion criteria, was used for the analyses. HRV parameters (RMSSD, VLF Power and LF Power) were calculated from 2-min segment in the rest phase before exercise and 2-min recovery period immediately after peak exercise. Area under the receiver operating characteristic curve (AUC) was used to determine the predictive performance for each parameter with and without HR corrections in rest and recovery phases. The division of HRV parameters by segment's RRavg to the power 2 (HRVDIV-2) showed the highest predictive performance under the rest phase (RMSSD: 0.67/0.66; VLF Power: 0.70/0.62; LF Power: 0.79/0.65; cardiac mortality/non-cardiac mortality) with minimum correlation to HR (r = −0.15 to 0.15). In the recovery phase, Kaplan-Meier (KM) survival analysis revealed good risk stratification capacity at HRVDIV-2 in both groups (cardiac and non-cardiac mortality). Although higher powers of correction (HRVDIV-4and HRVDIV-8) improved predictive performance during recovery, they induced an increased positive correlation to HR. Thus, we inferred that predictive capacity of HRV during rest and recovery is augmented when its dependence on HR is weakened by applying appropriate correction procedures.


INTRODUCTION
Heart rate (HR) recovery and heart rate variability (HRV) have been used by researchers for assessing the role of autonomic regulation in predicting all-cause and cardiovascular mortality (Freeman et al., 2006). The prognostic capabilities of HR response to exercise and after exercise have been well-documented (Lauer et al., 1996;Cole et al., 1999;Lipinski et al., 2004;Jouven et al., 2005;Kiviniemi et al., 2011) and reviewed by Freeman et al. (2006). Increased sympathetic and decreased parasympathetic activities have been associated with an enhanced risk of sudden death or the vulnerability to ventricular arrhythmias (Lahiri et al., 2008). Subdued time-and frequency-domain HRV indices have been linked with increased risk of mortality in the Framingham cohort (Tsuji et al., 1994), survivors of acute myocardial infarction (MI) (Kleiger et al., 1990;Kiviniemi et al., 2007) and cardiovascular morbidity and mortality (Zuanetti et al., 1996). However, studies determining the prognostic capacity of exercise induced short-term HRV have been sparse and inconsistent. Leino et al. (2010) concluded that none of the HRV indices were good predictors of mortality during peak exercise or recovery phase. In a study by Dewey et al. (2007), a greater short-term HRV during recovery post exercise was associated with an increased risk for all-cause and cardiovascular mortality. This is in contrast to observations made in resting HRV, which implies higher RR-interval variability is associated with better prognosis (Dekker et al., 2000;Leino et al., 2010). Nieminen et al. (2007) justified that the non-linear inverse relationship between RR interval and HR could be the cause for misinterpretation when comparing subjects with different HR levels and this has been concurred by other researchers (Chiu et al., 2003;Virtanen et al., 2007;Bailón et al., 2011). Possible physiological mechanisms involved have also been probed (Perini and Veicteinas, 2003;Goldberger et al., 2006). The non-linear relation between HR and HRV has been addressed by Sacha and Pluta (2008) and correction methods have been suggested to strengthen or weaken the influence of HR (Sacha et al., 2013a;Sacha, 2013). By determining whether decreasing dependence on HR improves the prognostic capacity of HRV, we sought to establish the influence of HR in predicting mortality risk. The aim of our study was to scrutinize these correction techniques and their influence on the predictive capacity of cardiac and non-cardiac mortality in the Finnish Cardiovascular Study (FINCAVAS) cohort.

PATIENT POPULATION AND FOLLOW-UP
A total of 2212 consecutive patients, who were referred by a physician and willing to undergo exercise stress tests at the Tampere University Hospital, were recruited between 2001 and 2004 for FINCAVAS. Informed consent was obtained from all the participants prior to the interview. Measurements were conducted as stipulated in the Declaration of Helsinki and the study protocol was approved by the Ethical Committee of the Hospital District of Pirkanmaa, Finland. In addition to raw electrocardiograph (ECG), descriptive information, medical history and habitual lifestyle of each patient were recorded. More detailed information regarding the patient population and sample size determination is described elsewhere (Nieminen et al., 2006). Of these, 1288 patients satisfied the inclusion criteria for this study with good quality HRV measurements for at least 2 min during rest phase, immediately prior to exercise, and 2 min during post-exercise recovery immediately after maximum effort.
The follow-up data consisted of information related to causes of death and was collected in 2011. The information for the follow-up was obtained from Causes of Death Register and has been shown to be reliable (Pajunen et al., 2005). The follow-up yielded 66 cardiac deaths and 94 non-cardiac deaths, while the remaining 1128 patients constituted the survival group.

EXERCISE TESTING PROTOCOL
The prognoses of mortality were analyzed using HRV indices obtained from 2 min segments during rest phase before exercise and 2 min recovery immediately after maximal exercise. Resting ECG was measured in the supine position prior to exercise. The exercise stress test was then performed on a bicycle ergometer with electrical brakes and the Mason-Likar modified lead system  (Mason and Likar, 1966) was used for the ECG data acquisition. Initial work load and increments were defined based on patient's age, gender, body mass index (BMI) and physical activity. Starting work load varied between 20 and 30 W and the stepwise increments ranged between 10 and 30 W every minute. ECG and HR were measured continuously during the test. Tests were sign-and symptom-limited with the recommended criteria for termination whereas in the case of post-MI patients, the upper limit for HR was set at 120-130 beats per minute (bpm). The chronotropic response index (CRI), which represents the chronotropic response to exercise, was evaluated as CRI = 100 × (peak HR − resting HR)/(220 − age − resting HR) (Kiviniemi et al., 2011). CRI < 80% was defined as low reserve capacity (Lauer et al., 1996). Measurement during the recovery phase was performed in the sitting position, immediately after exercise.

HRV MEASUREMENT
ECG was recorded at a sampling frequency of 500 Hz using CardioSoft exercise ECG system (Version 4.14, GE Healthcare, Freiburg, Germany) and was analyzed using Modified CASE software (GE Healthcare, Freiburg, Germany). After producing the RR-interval tachogram, the data was preprocessed to remove abnormal intervals and artifacts before they were divided into shorter segments based on the stages of rest and recovery. HRV parameters were determined using the Kubios HRV analysis software (Tarvainen et al., 2014). All intervals were resampled using cubic spline interpolation at 4 Hz. Linear and smoothness prior (smoothing parameter, λ = 500) detrending were performed prior to calculating time-domain parameters. The fast Fourier transform (FFT) spectrum was computed with a window width of 240 samples, which corresponds to the length of 1 min segment with 4 Hz resampling rate. A 50% overlapping window was used for longer segments. Mean RR intervals (RR avg ) were calculated from each segment individually for the HR correction procedure. Post-exercise recovery is marked by sympathetic withdrawal and parasympathetic reactivation. Sympathetic activation and attenuated parasympathetic recovery are significantly associated with adverse prognosis. The parameters included for examination were chosen based on previous HRV studies on mortality prediction and its outcomes. Of the spectral measures, low frequency (0.04-0.15 Hz, LF) power has been found to increase during exercise in normal subjects and reflects both sympathetic and vagal influences (Malliani et al., 1991). In addition, higher log LF power during recovery significantly predicted increased risk of all-cause and cardiovascular mortality (Dewey et al., 2007). Bigger et al. (1993) demonstrated that spectral measures from short segments (2-15 min) correlated significantly with those computed using 24-h periods. Bernardi et al. (1996) indicated that very low frequency (0.0033-0.04 Hz, VLF) power fluctuations were highly dependent on changes in physical activity, rather than preconceived notion of reflecting autonomic tone and thereby, emphasized the importance of activity as a confounding factor. Therefore, VLF power was evaluated due to its independent risk stratification property for all-cause mortality in patients

FIGURE 1 | Predictive performance of heart rate variability (HRV) parameters for: (A) cardiac mortality and (B) non-cardiac mortality groups.
Area under the receiver operating characteristics curves (AUC) and correlation coefficients (r), between HRV parameters and HR, for different correction methods during rest and recovery after exercise. AUC > 0.5 indicates that higher heart rate variability (HRV) is associated with better prognosis and AUC < 0.5 indicates higher HRV is associated with worse prognosis.
with acute MI (Bigger et al., 1993). Although high frequency (0.15-0.4 Hz, HF) power has been frequently used to measure parasympathetic tone in resting HRV, interpreting values during recovery after exercise is complicated due to tonic autonomic activity and residual adrenergic activity (Dewey et al., 2007). Goldberger et al. (2006) demonstrated that short-term (as small as 30 s windows) root mean squared difference of successive RR intervals (RMSSD), which represents high frequency variations in HR, is adequate for measuring parasympathetic reactivation in recovery phase.

HR CORRECTION
As described by Sacha et al. (2013a), the HRV dependence on HR is strengthened or weakened by multiplying or dividing the HRV indices by the corresponding segment's RR avg , respectively. In addition to normal determination of HRV indices, eight other classes for the indices were assessed in this study: HRV MUL-0.5 -multiplying HRV indices by RR avg to the power 0.5; HRV MUL-2 -multiplying HRV indices by RR avg to the power 2; HRV MUL-4 -multiplying HRV indices by RR avg to the power 4; HRV DIV-0.5 -dividing HRV indices by RR avg to the power 0.5; HRV DIV-2 -dividing HRV indices by RR avg to the power 2; HRV DIV-4 -dividing HRV indices by RR avg to the power 4; HRV DIV-8 -dividing HRV indices by RR avg to the power 8; and HRV DIV-16 -dividing HRV indices by RR avg to the power 16. With these classes, different levels of dependence/independence to HR were attained and can be considered significant in determining the contribution of HR in prognosis of cardiac and non-cardiac mortalities.

STATISTICAL ANALYSES
The relative risks for cardiac and non-cardiac mortality were assessed for individual characteristics, clinical condition and medication using univariate Cox models. The measure of the predictive power for different HR correction methods for each segment was computed using area under the receiver operating characteristics (ROC) curve. Spearman's rank correlation was performed to determine the degree of correspondence to HR. The cut-off points for Kaplan-Meier (KM) survival analyses were defined from the ROC analyses for each segment. The point of highest overall predictive performance (average of sensitivity and specificity) was chosen as the cut-off to distinguish mortality and survival groups based on HRV observed in the patient population. It has to be noted that these cut-off points were not optimized in order to preserve uniformity during comparisons. The Log-rank chi-square estimates were then used to evaluate the significance of the correction methods based on this classification.

RESULTS
During the follow-up of the patients who satisfied the inclusion criteria, 66 cardiac deaths were recorded, which included 31 sudden cardiac deaths, with a mean follow-up time of 54 months (min: 4.8 days; max: 99.5 months). 94 patients died of noncardiovascular causes between 1.2 and 110.7 months of follow-up (mean: 60.2 months). The baseline characteristics, clinical conditions and medications used by patients who suffered cardiac and non-cardiac deaths are listed in Table 1.
The area under the ROC curve (AUC) for HR was found to be 0.57/0.70 (rest/recovery) for cardiac mortality and 0.53/0.64 for non-cardiac mortality, implying that HR is a better predictor during recovery than during rest phase. The AUC for RMSSD, VLF and LF power, calculated under different correction methods during rest and recovery phases are presented in Figures 1A,B. Correlation with HR (r, presented in Figure 1) indicated increasing dependence or independence of HRV to HR, based on the method of correction used. AUC > 0.5 suggested that higher HRV are indicative of better prognosis. HRV DIV-2 , which revealed minimum correlation to HR, was the best predictor for both outcomes (cardiac and non-cardiac mortality) in the rest phase. However, during recovery, higher standard HRV (i.e., HRV without correction) was associated with worse prognosis (AUC < 0.5), as seen in Figure 1. In addition, similar associations were observed for HRV parameters multiplied by different powers of RR avg (HRV MUL-0.5 , HRV MUL-2, and HRV MUL-4 ). Conversely, after division by higher powers of RR avg (i.e., for HRV DIV-2 , HRV DIV-4 , HRV DIV-8 , and HRV DIV-16 ), higher HRV was associated with better prognosis (AUC > 0.5). Though higher orders of correction resulted in better predictive capacity, it also induced moderate/strong positive correlation to HR (in the case of HRV DIV-4 , HRV DIV-8 , and HRV DIV-16 ).
These results were further corroborated by KM survival analysis. Log-rank estimates at different degrees of correction for both cardiac and non-cardiac mortality are presented in Table 3. HRV MUL-4 and HRV DIV-16 were excluded due to their very strong correlation to HR. Mortality prediction was most significant for HRV DIV-2 in the rest phase. During recovery, the division of HRV by higher powers of RR avg resulted in better risk stratification for cardiac and non-cardiac deaths. Although HRV DIV-4 and HRV DIV-8 exhibited better predictive powers during recovery, the HRV indices exhibited strong positive correlation to HR (r = 0.6  www.frontiersin.org June 2014 | Volume 5 | Article 208 | 7 to 0.9 across both groups, as shown in Figure 1) at these correction levels. On the contrary, HRV DIV-2 was a good predictor of cardiac (p < 0.001) and non-cardiac (p < 0.05) mortality during recovery, with minimum influence of HR (r = −0.15 to 0.15). Figures 2, 3 represent the survival curves for HRV DIV-2 during rest and recovery.

DISCUSSION
The HRV indices computed from RR-interval measurements correlated with HR as a result of the non-linear relationship between the RR-interval and instantaneous HR (Chiu et al., 2003;Nieminen et al., 2007;Virtanen et al., 2007;Bailón et al., 2011). Higher variability during rest and lower variability during recovery were associated with better prognosis, and this corresponds to observations made by Dewey et al. (2007). Our results indicate that the predictive capacity of HRV at rest was highest when the correlation to HR was minimum (HRV DIV-2 , r = −0.15 to 0.15), suggesting that exclusion of HR influence on resting HRV improved prognostic capacity for cardiac and non-cardiac mortality. Since HR is a poor predictor at rest, removal of HR influence perchance resulted in improved prognostic capacity. On the contrary, HR during recovery phase exhibited significant risk stratification for both outcomes. Thus, increasing HRV's dependence on HR enhanced its predictive capacity (observed in HRV DIV-4 and HRV DIV-8 ).
However, higher degrees of correction produced moderate/strong positive correlation to HR, similar to observations made by Sacha et al. (2013c), Sacha (2014). To attain true independence, the correction technique that yields HRV least influenced by HR, needs to be identified. In our study, HRV DIV-2 demonstrated improvement in predictability of mortality risk during recovery phase with minimum dependence on HR. However, conclusive evidence could not be established to distinguish between cardiac and non-cardiac related deaths. This is in contrast to findings by Sacha et al. (2013b), who suggested that increasing the HRV dependence on HR resulted in greater predictive ability for cardiac death and increasing its independence indicated greater predictive power for non-cardiac death. One possible explanation could be that the study population analyzed by Sacha and coworkers comprised only post-MI patients whereas the current study included more heterogeneous patient material.
This study suffered certain limitations. First, the risk factors for individual, clinical conditions and medication were not modeled to determine their contribution toward mortality prediction. By including these variables to the analyses, a more definite conclusion on the cause of mortality could have been established. Second, the patients were not controlled for the type of medication prescribed. The effects of beta blockers and nitrates have been known to affect HR, which could have an effect on the results of HR correction. However, the purpose of the current study was to evaluate the effects of HR correction methods in mortality prediction and therefore, these issues need to be considered in future studies.

CONCLUSION
The findings of this study indicate that the predictive power of HRV parameters for both cardiac and non-cardiac mortality is augmented when its dependence on HR is weakened during rest and recovery. In addition, when HR is a good predictor, increasing HRV's dependence on HR further enhances the risk stratification for both modes of death.

AUTHOR CONTRIBUTIONS
The study was conceptualized by Tuomo Nieminen, Kjell Nikus, Terho Lehtimäki, Mika Kähönen, and Jari Viik. Data acquisition and analysis was performed by Paruthi Pradhapan, Mika P. Tarvainen, Rami Lehtinen, and Jari Viik. All authors contributed equally in drafting and revising the manuscript.