The Autonomic Imbalance of Myocardial Ischemia during Exercise Stress Testing: Insight from Short-Term Heart Rate Variability Analysis

Exercise stress testing (EST) has limited power in diagnosing obstructive coronary artery disease (CAD). The heart rate variability (HRV) analysis might increase the sensitivity of CAD detection. This study aimed to evaluate the correlation between short-term HRV and myocardial ischemia during EST, including the acceleration, maximum, and recovery stages of heart rate (HR). The HRV during EST from 19 healthy (RHC) subjects and 35 patients with CAD (25 patients with insignificant CAD (iCAD), and 10 patients with significant CAD (sCAD)) were compared. As a result, all HRV indices decreased at the maximum stage and no significant differences between iCAD and sCAD were found. The low-frequency power of heart rate signal (LF) of the RHC group recovered relatively quickly from the third to the sixth minutes after maximum HR, compared with that of the sCAD group. The relative changes of most HRV indices between maximum HR and recovery stage were lower in the sCAD group than in the RHC group, especially in LF, the standard deviation of all normal to normal intervals (SDNN), and the standard deviation in the long axis direction of the Poincaré plot analysis (SD2) indices (p < 0.05). The recovery slope of LF was significantly smaller in the sCAD group than in the RHC group (p = 0.02). The result suggests that monitoring short-term HRV during EST provides helpful insight into the cardiovascular autonomic imbalance in patients with significant CAD. The relative change of autonomic tone, especially the delayed sympathetic recovery, could be an additional marker for diagnosing myocardial ischemia.


Introduction
Coronary artery disease (CAD), a narrowing or obstruction of coronary arteries, is the world's most common cause of death [1]. Therefore, early diagnosis of CAD plays an important role in the primary prevention of cardiac death [1]. CAD leads to abnormal electrocardiograms (ECGs) and corresponding arterial waveforms. The abnormal features are enhanced for a CAD patient in exercise. The effects of exercise on the heart can be 2 of 12 monitored via an ECG for a subject running on a treadmill in exercise stress testing (EST). Thus, EST is commonly used as a clinical approach to diagnosing CAD [2]. However, the diagnostic accuracy of EST can be influenced by the patient's age, sex, or clinical characteristics [3]. EST has limited power to rule in or rule out obstructive CAD compared with diagnostic imaging tests [4]. In practice, EST has a low sensitivity of 68% and specificity of 77% for the detection of CAD [5]. The positive results of EST should be treated with caution because of its high false positive rate. According to the 2019 European Society of Cardiology (ESC) Guidelines for diagnosing and managing chronic coronary syndromes, EST is only recommended for risk assessment rather than for a diagnosis of CAD [4]. To clarify the diagnosis of CAD, additional coronary anatomic imaging or computed tomography (CT) is required [3], which is time-consuming or expensive. Therefore, optimizing the accuracy of EST would reduce the financial burden and medical resources.
On the other hand, heart rate variability (HRV) is a relevant marker reflecting the cardiac variation by the sympathetic and vagal components of the autonomic nervous system. Autonomic dysfunction directly correlates with the morbidity and mortality caused by CAD and involves the development of cardiovascular disease or the progression of metabolic disease [6,7]. For example, vagal-mediated HRV indices can be applied for recognizing healthy and diseased states, and have been inversely associated with metabolic diseases, including diabetes, central obesity, dyslipidemia, and hypertension [7]. In the vagal-mediated HRV analysis, long-term HRV was largely applied to predict sudden cardiac death [7]. Short-term HRV is also a method used for enhanced risk assessment in low-to intermediate-risk individuals without known CAD [7,8]. As a dynamic marker while one experiences different loads, HRV appears to be sensitive and responsive to acute stress, including exercise [7,9].
However, the diagnostic potential of exercise-related HRV was still controversial in diagnosing CAD. Bailon et al. identified that HRV indices corrected by mean heart rate (HR) and respiratory frequency could improve the accuracy of EST, ranging from 76% to 95% [10]. Yet, some previous findings reported that the values of HRV indices alone or corrected with HR were insufficient in distinguishing CAD detection [11,12]. Therefore, this study aimed to evaluate the correlation between short-term HRV and myocardial ischemia during EST, including the acceleration, maximum, and recovery stages of heart rate (HR).

Exercise Stress Testing (EST)
The EST follows the Bruce protocol. It consists of seven stages, with a gradual incremental increase in the speed and gradient of the treadmill during the EST. The predicted peak heart rate was calculated as 220-age. Individuals were encouraged to exercise until they experienced limiting symptoms, even if 85% of the maximum predicted heart rate was achieved [13]. Criteria for exercise termination were physical exhaustion or maximum heart rate greater than the age-predicted maximum. During each exercise stage and recovery stage, symptoms, blood pressure, heart rate, and exercise workload in metabolic equivalents (METS) were recorded. Each stage lasted three minutes to allow the patient to habituate to increasing speed and gradient before advancing to the next stage. Following peak exercise, individuals walked for a two-minute cool-down period at 1.5 mph and 2.5% grade. The exercise tests were performed, analyzed, and reported with a standard protocol utilizing a computerized database. The results of the EST were reconfirmed by two cardiologists.

RR Intervals Extraction
All patients underwent short-term electrocardiographic recording, which was performed using the "Medilog" Digital Holter Recorder (Medilog Darwin v2.3 AR4 plus) during EST. RR intervals were extracted from subjects' ECG during EST. Ectopic and other artifactual beats within the RR interval series were corrected using the interpolation method [14]. Drift in RR interval series was minimized using empirical mode decomposition. After the preprocessing, nine segments of the 1 min RR interval series were extracted at different stages before and after the exercise. The 9 stages include the rest period before exercise (Res), the first minute before maximal heart rate (Pre), the first minute when the maximum heart rate was reached (Max), and every 1 min of 6 min in the post-exercise recovery phase (R1~R6).

Analysis of Heart Rate Variability
HRV analysis is suitable for evaluating the short-term autonomic regulation of heart rate under physiologically stable conditions [15]. HRV analysis was applied to the nine one-minute RR interval series of all 54 patients. The HRV metric includes indices from time and frequency domain analysis and non-linear Poincaré plot analysis.
Time domain HRV indices include the root mean square of successive RR differences (RMSSD), the standard deviation of all normal to normal intervals (SDNN), the mean value of the RR interval series (MRRI), and the mean of heart rate (MHR). SDNN reproduces the total variability, and RMSSD reflects parasympathetic activity [7]. Frequency domain HRV indices include the total power from 0.0033 to 0.4 Hz (TP), the very low-frequency power from 0.0033 to 0.04 Hz (VLF), the low-frequency power from 0.04 to 0.15 Hz (LF), the high-frequency power from 0.15 to 0.4 Hz (HF), the ratio of low-frequency power to highfrequency power (LF/HF), the normalized low-frequency power (LF (nu) = LF/(LF + HF)), and the normalized high-frequency power (HF (nu) = HF/(LF + HF)). Index LF referring to the modulation of the RR interval changes corresponds to the sympathetic and parasympathetic activity together. HF modulation of the RR interval changes is primarily regulated by the parasympathetic nerve through the innervations of the heart [7]. Non-linear indices of Poincaré plot analysis include the standard deviation in the short axis direction (SD1), the standard deviation in the long axis direction (SD2), and SD1/SD2. SD1 is the fast beat-tobeat variability in the RR intervals, while SD2 describes the longer-term variability. SD1 reflects mainly the parasympathetic input to the heart, while SD2 reflects the sympathetic and parasympathetic contributions to the heart [7,16].
Further, for each HRV index of a patient, the slope of the linear regression from Max to R6 stages was calculated to quantify the recovery ability of the patient.

Statistical Methods
The HRV data were presented as mean ± standard deviation. The degree of difference for each HRV index between any two groups was performed in terms of p value using the Wilcoxon rank-sum test (which is suitable for small sample comparison based on nonparametric statistics). The statistical significance was set at p < 0.05. All data analyses were performed using Python's open-source statistical package, Pingouin v0. 3

Results
A total of 54 patients were enrolled in the study. They were divided into three groups: RHC (n = 19; mean age: 51.9 ± 8.8 years; male/female, 13/6), sCAD (n = 10; mean age: 54.7 ± 10.0 years; male/female, 8/2), and iCAD (n = 25; mean age: 53.4 ± 8.9 years; male/female, 16/9). There were no significant differences in sex and age among the groups. Figure 1 compares the HRV indices between the RHC and the sCAD groups at all 9 stages. From Res to Max stages, most HRV indices decreased except for the MHR. From R3 to R6 stages, most HRV indices of the RHC group, except for the HF (nu), recovered faster than those of the sCAD group. Note that the mean value of LF in Max stage was nearly negligible. This indicates that the values of LF in Max stage were nearly the same for all subjects in both groups. This is a unique feature to be applied later.

Comparison of HRV between RHC and sCAD Groups in Individual EST Stages
on non-parametric statistics). The statistical significance was set at p < 0.05. All data analyses were performed using Python's open-source statistical package, Pingouin v0.3.3. [17]

Results
A total of 54 patients were enrolled in the study. They were divided into three groups: RHC (n = 19; mean age: 51.9 ± 8.8 years; male/female, 13/6), sCAD (n = 10; mean age: 54.7 ± 10.0 years; male/female, 8/2), and iCAD (n = 25; mean age: 53.4 ± 8.9 years; male/female, 16/9). There were no significant differences in sex and age among the groups. Figure 1 compares the HRV indices between the RHC and the sCAD groups at all 9 stages. From Res to Max stages, most HRV indices decreased except for the MHR. From R3 to R6 stages, most HRV indices of the RHC group, except for the HF (nu), recovered faster than those of the sCAD group. Note that the mean value of LF in Max stage was nearly negligible. This indicates that the values of LF in Max stage were nearly the same for all subjects in both groups. This is a unique feature to be applied later.   Table 1 shows a comparison of HRV indices in different recovery stages between the RHC and the sCAD groups. In the R1 and R2 stages, there was no significant difference between the two groups. In the R3 stage, the VLF, LF, and LF (nu) were lower in the sCAD group than in the RHC group. In the R4 stage, the VLF, LF, LF/HF, TP, LF (nu), rMSSD, and SD1 were lower in the sCAD group. In the R5 stage, the SDNN and SD2 were lower in the sCAD group. In the R6 stage, the LF, SDNN, and SD2 were lower in the sCAD group. Among the HRV indices listed in Table 1, the LF in R4 (35.77 ± 64.06 vs. 3.39 ± 4.3 ms 2 ) exhibited the best performance (p < 0.01), which is labeled in bold in Table 1.  Figure 2 shows the difference in the value of each HRV index between each stage and the Max stage. From the Res to Max stage, the relative reductions of the HRV indices were significantly higher in the RHC groups than in the sCAD group, especially for SDNN, SD2 and total power. From the Max to R4 stage, the relative increase of LF was significantly higher in the RHC group than in the sCAD group. From the Max to R6 stage, the relative increase of HRV was significantly higher in RHC groups, especially for SDNN and SD2. Table 2 shows a comparison of the HRV indices between the RHC and sCAD groups in individual EST stages with respect to Max stage. Only particular stages are listed in which the HRV indices between the two groups are significantly different. For the purpose of comparison, LF, HF, and LF/HF are also listed. Yet, HF and LF/HF showed no statistical difference between the two groups. Among the HRV indices listed in Table 2, the SDNN (34.45 ± 21.28 vs. 11.08 ± 17.89 ms) and SD2 (48.3 ± 30.04 vs. 14.89 ± 25.42 ms) in Res with respect to Max exhibited the best performance (p < 0.01). Table 2. A comparison of HRV indices between RHC and sCAD groups in individual EST stages with respect to Max stage. Only particular stages are listed in which the HRV indices between the two groups are significantly different. For the purpose of comparison, LF, HF, and LF/HF are also listed, which showed no statistical difference between the two groups. The SDNN and SD2 in Res with respect to Max, in bold, exhibited the best performance (p < 0.01).   Table 2 shows a comparison of the HRV indices between the RHC and sCAD groups in individual EST stages with respect to Max stage. Only particular stages are listed in which the HRV indices between the two groups are significantly different. For the purpose of comparison, LF, HF, and LF/HF are also listed. Yet, HF and LF/HF showed no statistical difference between the two groups. Among the HRV indices listed in Table 2, the SDNN (34.45 ± 21.28 vs. 11.08 ± 17.89 ms) and SD2 (48.3 ± 30.04 vs. 14.89 ± 25.42 ms) in Res with respect to Max exhibited the best performance (p < 0.01). According to Table 2, the recoveries of SDNN from Max to R6, SD2 from Max to R6, and LF from Max to R4 in the RHC group were significantly higher than those in the sCAD group. Figure 3 illustrates the corresponding recovery amount of SDNN, SD2, and LF of each patient in the RHC and the sCAD groups. Green and red arrows indicate an increase and decrease in HRV values, respectively. In the RHC group, most patients had a significant increase in SDNN, SD2, and LF indices. Furthermore, the arrow colors of LF tend to be more consistent for all patients than those of SDNN or SD2 in each group.

The Recovery Slope of HRV of the RHC, iCAD, and sCAD Groups
To quantify the recovery intensity of each HRV index, we calculated the slope of the linear regression from the Max to the R6 stages of the index for each subject in the RHC, sCAD, and iCAD groups. Table 3 shows the population characteristics of the 54 patients and their recovery slope (mean ± standard error) of each HRV index with a p value between any two groups during the recovery phase. For comparison, the slopes of LF, SDNN, and SD2 were significantly lower in the sCAD group than in the RHC group. Among the three indices, the slope of LF exhibited the best performance (3.93 ± 5.95 vs. 0.85 ± 0.68 ms 2 /min, p = 0.02) to separate the sCAD group from the RHC group. On the contrary, although the slope of LF still exhibited the best performance (3.82 ± 3.82 vs. 0.85 ± 0.68 ms 2 /min, p = 0.06) between the iCAD and sCAD groups, there was no significant difference (p = 0.06). Similarly, there was also no significant difference (p = 0.57) between the iCAD and RHC groups. and LF from Max to R4 in the RHC group were significantly higher than those in the sCAD group. Figure 3 illustrates the corresponding recovery amount of SDNN, SD2, and LF of each patient in the RHC and the sCAD groups. Green and red arrows indicate an increase and decrease in HRV values, respectively. In the RHC group, most patients had a significant increase in SDNN, SD2, and LF indices. Furthermore, the arrow colors of LF tend to be more consistent for all patients than those of SDNN or SD2 in each group.  Table 2.

Discussion
EST increases the load of the heart via physical exercise and is a risk assessment method for screening asymptomatic CAD patients [1]. However, the accuracy of EST was lower for diagnosing CAD and could be influenced by the patient's gender, age, and underlying disease [3,5]. A meta-analytic review in 2012 reported a positive likelihood ratio  Table 3. Population characteristics and the recovery slope (means ± standard deviation) of each HRV index with a p value between any two groups during the recovery phase. Bold highlights the best performances. The symbol * indicates a significant difference at the 0.05 level (p < 0.05).

Discussion
EST increases the load of the heart via physical exercise and is a risk assessment method for screening asymptomatic CAD patients [1]. However, the accuracy of EST was lower for diagnosing CAD and could be influenced by the patient's gender, age, and underlying disease [3,5]. A meta-analytic review in 2012 reported a positive likelihood ratio of 3.57 and a negative likelihood ratio of 0.34 for EST [3]. Therefore, EST is only recommended for risk assessment rather than for diagnosing CAD, due to its poor diagnostic accuracy [4]. Although invasive coronary angiography (CAG) is the gold standard for detecting CAD, it may be accompanied by a possible risk of arrhythmia, myocardial infarction, stroke, and some artery injury [1]. Therefore, adding a non-invasive imaging study, energy defect of Hilbert-Huang transform (HHT), or HRV with EST could help to increase the diagnostic accuracy of CAD [4,18].
Physical exercise is associated with parasympathetic withdrawal and sympathetic increase, resulting in a heart rate increase. The analysis of HRV has been used to predict cardiovascular mortality, but is not usually used in predicting CAD due to inconsistent results [12,[19][20][21]. After exercise cessation, HR and HRV demonstrate a time-dependent recovery and return to pre-exercise levels. The initial increasing HR during exercise is due to parasympathetic withdrawal, and recovery of the HR immediately after exercise is mediated by parasympathetic reactivation in response to the activity in arterial baroreceptors [22,23]. The heart rate recovery (HRR), especially in the first 30 s after exercise, reflects the balance of reactivating parasympathetic drive and withdrawal of sympathetic drive [19]. The imbalance of autonomic function, especially the delayed reactivation of the parasympathetic nervous system, plays an important role in HHR and predicts cardiovascular mortality [19][20][21]24]. The last phase of HRR has been attributed to parasympathetic reactivation [25]. Slow HRR of the early post-exercise period reflects the autonomic imbalance and could predict all-cause mortality, sudden death, or silent ischemia [19,22,23] indepen-dent of exercise workload, myocardial perfusion defects, or HR changes during exercise [26]. However, the possible mechanisms mediating post-exercise cardio-deceleration are not well known. They could be related to autonomic denervation of the myocardium, a higher pain threshold during EST, higher β-endorphin, or anti-inflammatory cytokines [22,23]. Slow HRR could predict microvascular or macrovascular myocardial ischemia, but the relation between HRR and the severity of coronary stenosis is still uncertain [22]. HRR in the supine position shortly after exercise may have a negative predictive value for CAD detection [27]. In our study, the sCAD group had a lower maximum HR and slower HRR during the R1 and R2 stage compared with the RHC group. During the recovery stage of EST in the sCAD group, a slower phase of cardio-deceleration was observed, likely mediated by both prolonged parasympathetic reactivation and sympathetic withdrawal.
The complexity of HRV is moderated by the dynamic interaction between the sympathetic and parasympathetic nervous systems [16]. Compared with HRR, the complexity of HRV is not yet fully understood in supposed correspondence with different HRV powers [28]. By applying these frequency range differences in HRV analysis, the individual contribution of parasympathetic and sympathetic systems could be identified by HRV analysis. Typically, LF-nu increases during low-moderate intensity exercise and decreases during higher-intensity exercise, while HF-nu demonstrates the opposite response [25]. High-intensity exercise training can chronically lead to a shift from vagal to sympathetic cardiac modulation, and the progressive sympathetic predominance at peak training load may predict performance during exercise [7]. A transition from parasympathetic to sympathetic predominance in cardiovascular autonomic modulation could be influenced by different intensity during exercise, duration of exercise, and acute post-exercise recovery [7]. Regular exercise would activate the parasympathetic activity and has been shown to increase resting HRV [28]. However, the role of HRV recovery after acute EST still remains unclear. Furthermore, the variation of the HF component is largely contributed to by respiratory sinus arrhythmia (RSA). Vagal-mediated RSA may fluctuate greatly with cardio-acceleration during inspiration and cardio-deceleration during expiration [7]. As the workload of the EST increases, decreasing HRV is noted when respiratory frequency increases [7]. Some patients may experience irregular respiratory rate during angina, thereby rendering difficult the proper interpretation of HRV data [10,12].
In this study, patients with significant CAD may have undergone a lower intensity or shorter duration of exercise, for they usually complained of angina earlier than the expected intensity or duration was reached. This tended to influence their heart rates and the corresponding HRV inconsistently and individually. In the Max stage, although the heart rates and exercise intensities were different among patients, the values of some HRV indices such as LF were found to remain the same, nearly negligible, for all subjects in the RHC and the sCAD groups, as shown in Figure 1. Thus, monitoring the HRV recovery after the Max stage may be an alternative method for assessing heart load recovery.
Among HRV indices, the HF power may fluctuate substantially, and LF/HF also demonstrates inconsistent responses to different intensities of exercise [25]. As they have insufficient power in distinguishing coronary stenosis or vessel number, LF, HF, or VLF components of HRV have been neglected in recent HRV studies [11,29]. The LF of HRV was significantly lower in patients with CAD [11] and could be another index for distinguishing myocardial ischemia. In a previous study, myocardial ischemia induced sympathetic reflex and delayed recovery of LF/HF during the early recovery period of EST [30].
In this study, the delayed reactivating parasympathetic tone of the sCAD group was observed in the HF index during the R3 and R4 stages. Withdrawal of sympathetic tone with prolonged lower LF index was found in patients with significant CAD during the whole recovery stage. Reduced HRV and HRR were risk factors of CAD after multi-variable regression analysis [11,31]. However, we found that delayed HRV recovery, especially in regards to the low recovery slope of LF, occurred in patients with significant CAD. In the recovery stage of EST, the recovery amounts of LF were significantly lower in the sCAD group than in the RHC group. Additionally, the recovery slope of LF power was also sig-nificantly lower in the sCAD group. Therefore, post-exercise HRV, especially in regards to the recovery slope of LF, could be an important tool for investigating autonomic imbalance in patients with significant CAD. Although there was no significant difference between the sCAD and iCAD groups, the mean values of the recovery slopes of LF in the sCAD, iCAD, and RHC groups are relatively small, medium, and large, respectively, as expected. The result implies that the slope of LF can be partially applied to reflect the severity of coronary stenosis. In future studies, combining machine learning techniques with multiple HRV indices may increase the accuracy in quantifying the severity of coronary stenosis.
This study possessed the limitations of being performed at a single institution and in a small patient population. In addition, the gender imbalance of participants may affect the results [32]. The difference in exercise parameters between the three groups is another limitation of this retrospective data. However, the autonomic imbalance was more significant during the post-exercise period irrespective of different exercise stress influences. In the future, the relation between decreasing sympathetic activity and CAD should be assessed in a larger number of subjects with pathophysiologic data. In conclusion, although it is possible that these limitations influenced the results, we conclude that, in this preliminary study, the decreasing sympathetic activity, especially as it relates to the recovery slope of LF, is useful for easily detecting myocardial ischemia in high-risk patients and could increase the accuracy of EST without unnecessary invasive treatments.

Conclusions
The main finding of the present study indicated that short-term HRV could reflect the cardiovascular autonomic imbalance in patients with significant CAD during EST. Moreover, the impaired response of HRV after exercise, especially the prolonged withdrawal of sympathetic tone, was found. A relatively slow recovery slope of LF after exercise could be an additional marker for diagnosing myocardial ischemia.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The anonymized data are available from the author upon reasonable request.