1 INTRODUCTION

Coronavirus disease 2019 (COVID-19), the illness caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), mainly affects the respiratory system and shows wide range of clinical presentations varying from asymptomatic/mild symptoms to severe viral pneumonia with respiratory failure, systemic inflammatory syndrome, and death (Gorbalenya et al., 2020; Wang et al., 2020; Zhou et al., 2020). Although the respiratory tract is the primary target for SARS-CoV-2 virus, cardiovascular involvement has been documented in acute phase of the disease in different studies (Alareedh et al., 2020; Madjid et al., 2020). In recovered patients, palpitations and dyspnea on exertion may prolong and adversely affect health even after full recovery from COVID-19. Extensive diagnostic workups of these patients may not show any reinfection, cardiac or respiratory involvement, and any laboratory abnormalities. This may indicate post-acute sequelae of COVID-19 in which there is altered autonomic function (Goldstein, 2021).

Heart rate variability (HRV) is considered to be a sensitive and noninvasive method for the quantitative assessment of sympathetic and parasympathetic branches of the autonomic nervous system. A heart rate that is variable, and responsive to demands is believed to bestow a survival advantage, whereas lower HRV is associated with higher risk of cardiovascular events and all-cause mortality (Fang et al., 2020; Thayer et al., 2010; Tsuji et al., 1996). HRV has been widely used for decades to quantify risk in a wide variety of both cardiac and noncardiac disorders (Kleiger et al., 2005).

Since recovered patients from COVID-19 experience some kind of dysautonomia, we hypothesized that COVID-19 may cause alterations in cardiac sympathetic and parasympathetic systemin post-acute phase of the disease. Moreover, HRV parameters after recovery from COVID-19 have not been evaluated so far. Therefore, in the present study, we investigated HRV parameters in patients after recovery from COVID-19 who had previously symptomatic and were followed in outpatient setting.

2 METHODS 2.1 Study population

The study population was recruited from the cardiology outpatient clinic of Adıyaman University Hospital between December 2020 and May 2021. The present study complied with the Declaration of Helsinki, and the study protocol was approved by the local ethics committee (Adıyaman University Institutional Ethics Committee).

The study group composed of 50 consecutive subjects with a confirmed history of COVID-19 infection who required no hospitalization or home oxygen treatment and without severe respiratory or other major organ involvements. The control group composed of 50 consecutive healthy subjects without a history of COVID-19 infection and vaccination. All healthy control subjects underwent a nasal swab test to exclude asymptomatic COVID-19 infection. The study group also was also tested for SARS-CoV-2 to exclude repeat infection.

Subjects were excluded if they had active COVID-19 infection, history of COVID-19 infection which required hospitalization, home oxygen treatment or severe respiratory involvement, hypertension, diabetes mellitus, beta-blocker or calcium channel blocker use, moderate-to-severe valvular heart disease, prosthetic heart valves, coronary artery disease, left ventricular dysfunction, atrial fibrillation, frequent atrial or ventricular extrabeats, chronic obstructive pulmonary disease, asthma, obstructive sleep apnea, BMI over 30 kg/m2, renal failure, cerebrovascular disease or thyroid disease, chronic liver disease, and inflammatory and autoimmune disorders.

2.2 Echocardiographic evaluation

All subjects underwent 2-dimensional and pulsed- and tissue-Doppler echocardiographic evaluation to exclude any heart disease that may affect HRV parameters. The following 2-dimensional echocardiographic parameters were measured: left ventricular end-diastolic diameter (LVEDD, mm), left ventricular end-systolic diameter (LVESD, mm), aortic root diameter (mm), left atrium diameter (LA, mm), interventricular septum thickness in diastole (IVST, mm), and posterior wall thickness in diastole (PWT, mm). The pulsed Doppler mitral inflow velocities were obtained from the apical four-chamber view with the sample volume placed just below the mitral leaflet tips and peak transmitral flow velocity in early diastole (E) and late diastole (A) were measured. Using tissue Doppler, the early (Em) diastolic velocities were assessed at the septal and lateral insertion of the annulus of mitral valve, and the average value between two measurements was determined. Afterward, E/A and E/Em ratios were calculated. All measurements were done according to the recommendations of the European Association of Cardiovascular Imaging (Galderisi et al., 2017).

2.3 HRV analysis

All the study participants underwent 24-hour ambulatory ECG using a 12-channel digital Holter recorder (iH-12Plus Holter System, Biocare, China) during normal daily activities with their normal sleep–wake. The recordings were analyzed by another cardiologist who was blinded to the study population. HRV measurements were obtained using a Holter software program (iH-12Plus Holter System Software, Biocare, China). After computerized primary analysis, all recordings were reviewed and edited manually for careful elimination of ectopic beats and artifacts. Only recordings less than 5% artifacts and ectopies were included for analysis. Both linear and nonlinear HRV parameters were derived from the same time series free from artifacts. Recordings of sufficient time (>22 h) for evaluation were included in the analysis.

The following time domain indices of HRV were derived from the 24-h Holter recordings: (1) The mean of all normal RR intervals in milliseconds (RR, ms); (2) the standard deviation of all normal-to-normal RR intervals (SDNN, ms); (3) the mean of the standard deviation of all normal-to-normal RR intervals for all 5-minute segments (SDNN index, ms); (4) the standard deviation of 5-minute mean RR intervals (SDANN, ms), (5) the percentage of successive normal RR intervals with a difference of more than 50 ms (pNN50, %); (6) the square root of the mean of the squares of the differences between successive normal-to-normal RR intervals (RMSSD, ms) and HRV index (HRVI).

The frequency domain indices of HRV were measured using a fast Fourier transform (FFT) on RR intervals (4 Hz) with a window length of 250 data points and 50% overlap. The following parameters were derived from the 24-hour Holter recordings: total power (TP, <0.4 Hz), high-frequency power (HF, 0.15–0.40 Hz), low-frequency power (LF, 0.04–0.15 Hz). Normalized units of LH and HF (LFnu and HFnu) were also calculated by dividing the power of the LF and HF components by TP minus VLF and multiplying by 100.

Nonlinear indices of HRV were derived from detrended fluctuation analysis (DFA), recurrence plot analysis (RPA), and entropy values. DFA included the short-term fractal coefficient (α1) and the long-term fractal coefficient (α2). These parameters are slopes of a log-log plot, from which α1 was derived by default with a window width of 4–11 beats and α2 within a window width of 11–64 beats. DFA gives information about fractal correlation properties of RR intervals. α1 component values of approximately 1.0 are those that shows greater proximity to normal physiological behavior (Gronwald & Hoos, 2020). RPA parameters included mean line length in recurrence plot (RP) (Lmean, beats), maximum line length in RP (Lmax, beats), percentage of recurrence points in RP (REC, %), determinism or percentage of recurrence points which form diagonal lines in RP (DET, %), and Shannon entropy of line length distribution (ShanEn). RPA is suitable for analysis of non-stationary sequences of RR intervals and higher values correspond to a pattern of low variable state (Tarvainen et al., 2014). Other nonlinear parameters of HRV were approximate entropy (ApEn) and sample entropy (SampEn), both of which refer to the rate at which a complex system produces information, and predict the complexity of a dynamic system (Richman & Moorman, 2000).

2.4 Statistical analyses

Statistical analyses were performed with the use of SPSS software, version 21.0 (SPSS Inc.). The Kolmogorov-Smirnov test was used to evaluate whether the variables were normally distributed. Continuous variables were presented as means ± SD or median with 25th-75th percentiles and categorical variables as frequency. Chi-square test or Fisher's exact test was used for categorical variables where appropriate. For continuous data, Student's t-test for was used normally distributed variables. Non-normally distributed variables were transformed logarithmically. All p values were two-tailed, and values of less than 0.05 were considered to indicate statistical significance.

3 RESULTS

Baseline characteristics of the study group (50 subjects, 29 women and 21 men; mean age 40.8 ± 10.3 years) and the control group (50 subjects, 26 women and 24 men; mean age 38.2 ± 12.0 years) are summarized in Table 1. There were no statistically significant differences between groups with regard to age, gender, and body mass index. The mean duration from COVID-19 to Holter recording was 20.0 ± 11.4 weeks.

TABLE 1. Baseline and echocardiographic characteristics of the study population Variables

Recovered COVID-19 group

(n = 50)

Control group

(n = 50)

p-value

Age, years 40.82 ± 10.31 38.24 ± 12.02 .252 Women/ men 29/21 26/24 .688 Duration to Holter recording, weeks 20.0 ± 11.4 — — BMI, kg/m2 25.70 ± 3.86 26.10 ± 4.50 .635 Ectopic beats 63.52 ± 27.24 60.74 ± 25.12 .540 Ejection fraction, % 63.12 ± 4.30 62.80 ± 4.54 .804 LVEDD, mm 47.60 ± 2.62 48.20 ± 2.40 .237 LVESD, mm 27.62 ± 3.90 26.82 ± 3.83 .341 Interventricular septum, mm 9.06 ± 0.93 9.12 ± 1.02 .760 Posterior wall, mm 9.34 ± 1.08 9.20 ±1.00 .784 Aortic diameter, mm 24.50 ± 2.05 25.12 ± 2.17 .141 Left atrial diameter, mm 32.34 ± 3.94 31.86 ± 4.33 .597 E/A ratio 1.57 ± 0.19 1.62 ± 0.16 .233 E/Em ratio 6.36 ± 3.45 6.25 ± 3.33 .251 Abbreviations: BMI, body mass index; LVEDD, left ventricular end-diastolic diameter; LVESD, left ventricular end-systolic diameter.

The two-dimensional and Doppler echocardiographic findings of the study and control groups are shown in Table 1. No differences between groups were found in two-dimensional, pulsed- and tissue-Doppler echocardiographic variables between groups.

Data of the time, frequency and nonlinear domain analysis for the study, and control groups are presented in Table 2. Regarding the time domain parameters, SDNN, SDANN, and SDNN index were significantly lower in the study group when compared to the control group (p < .0001 for all). HRVI was also significantly lower in patients with a past history of COVID-19 than in those without a past history of COVID-19 (Figure 1 and Table 2). pNN50 and RMSSD were depressed in patients with a past history of COVID-19 (p < .0001 for all). RR intervals were not different between groups.

TABLE 2. Time, frequency, and nonlinear HRV parameters of the study population Variables

Recovered

COVID-19 group

(n = 50)

Control group

(n = 50)

p-value RR, ms 757 ± 88 777 ± 78 .227 SDNN, ms 122.40 ± 30.90 161.30 ± 30.80 <.0001 SDANN, ms 113.90 ± 30.20 144.70 ± 33.60 <.0001 SDNN index, ms 50.10 ± 13.40 63.60 ± 14.80 <.0001 pNN50 1.03 ± 0.29 1.23 ± 0.25 <.0001 RMSSD 1.45 ± 0.16 1.62 ± 0.18 <.0001 HRVI 17.20 ± 5.05 20.40 ± 4.95 .030 Total power 3.36 ± 0.24 3.55 ± 0.21 <.0001 VLF power 3.17 ± 0.23 3.33 ± 0.20 .001 LF power 2.71 ± 0.31 2.95 ± 0.28 <.0001 HF power 2.29 ± 0.33 2.62 ± 0.34 <.0001 LFnu 69.60 ± 11.60 67.80 ± 13.90 .482 HFnu 28.10 ± 11.10 34.30 ± 16.60 .033 α1 0.82 ± 0.23 0.95 ± 0.18 .002 α2 0.89 ± 0.14 0.98 ± 0.12 .001 Lmax, beats 420 ± 71 378 ± 90 .009 Lmean, beats 33.55 ± 10.20 27.11 ± 12.43 .004 DET, % 98.67 ± 0.49 96.33 ± 0.84 <.0001 REC, % 36.91 ± 7.07 29.90 ± 7.82 <.0001 ShanEn 3.91 ± 0.42 3.59 ± 0.51 .001 ApEn 0.81 ± 0.19 0.93 ± 0.17 .002 SampEn 0.92 ± 0.19 1.12 ± 0.43 .003 Note Data are presented with mean and SD. Abbreviations: α1, short-term fractal coefficient; α2, long-term fractal coefficient; ApEn, approximate entropy; DET, percentage of recurrence points which form diagonal lines in recurrence plot; HF, high frequency; HFnu, normalized high frequency; HRV, heart rate variability; HRVI, heart rate variability triangular index; LF, low frequency; LFnu, normalized low frequency; Lmax, maximum line length in recurrence plot; Lmean, mean line length in recurrence plot; pNN50, percentage of successive normal RR intervals exceeding 50 milliseconds; REC, percentage of recurrence points in recurrence plot; RMSSD, the square root of the mean of the squares of the differences between successive normal-to-normal RR intervals; RR, mean RR interval; SampEn, sample entropy; SDANN, the standard deviation of 5-minute mean RR intervals; SDNN index, the mean of the standard deviation of all normal-to-normal RR intervals for all 5-minute segments; SDNN, the standard deviation of all normal-to-normal RR intervals; ShanEn, Shannon entropy; VLF, very low frequency.

Heart rate variability triangular index in the control and study groups

Regarding the frequency domain parameters, TP, VLF, LF, and HF were significantly depressed in the study group in comparison with the control group. However, when LF and HF were converted into their normalized forms, HF continued to be depressed in the study group, while LF was not. There was no correlation between any HRV parameter and the duration from COVID-19 to Holter recording.

Nonlinear parameters of HRV were statistically significantly different between patients with a past history of COVID-19 and without a past history of COVID-19. α1 and α2 components of DFA were depressed in the study group when compared with the control group. By contrast, RPA parameters including Lmax, Lmean, DET, REC, and Shannon entropy increased in the study group in comparison with the control group. ApEn and SampEn were depressed in the study group compared to the control group.

4 DISCUSSION

The main findings of our study are as follows: (1) SDNN, SDANN, SDNNi, and HRVI, all of which are time domain parameters, are decreased and indicate altered autonomic function (2) reduced power in most spectral bands further implies altered autonomic function in this population, (4) RMSSD, pNN50, HF, and HFnu, all of which are a specific expression of vagal activity, are depressed, thus confirming decreased parasympathetic activity in this population, (5) nonlinear domain measures including those which extracted from DFA, RPA, and entropy-based analysis are all altered in this population. To the best of our knowledge, this is the first study to have investigated all three domains of HRV in a cohort of patients who recovered from COVID-19.

HRV was studied in the hospitalized COVID-19 patients with severe or critical infections as well as in outpatient settings. Kaliyaperumal et al compared 5 min HRV parameters in COVID-19 patients with that of healthy population and found that frequency parameters of HRV, which are LF and HF, were significantly decreased among COVID-19 patients as compared to healthy individuals. In addition, they found that rMSSD was significantly higher in the COVID-19 group compared to healthy individuals (Kaliyaperumal et al., 2021). In line with above findings, we found reduced LF and HF bands after recovery from COVID-19. However, when normalized units were used, only HF continued to be decreased while LF was not in our study. The presentation of LF and HF in normalized units highlights controlled and balanced behavior of the two branches of the autonomic nervous system. Furthermore, the normalization seems to decrease the effects of the changes in TP on the values of LF and HF (Zhou et al., 2020). In our study, both HF and HFnu are reduced, while LFnu are similar between groups. This may be due to the fact that HF power reflects primarily parasympathetic influences and LF power has been shown to reflect both sympathetic and parasympathetic influences (Kleiger et al., 2005). In another study, Hasty et al evaluated SDNN using 7 min short-segment HRV data together with CRP levels in patients presented with hypoxic respiratory failure requiring high-flow nasal cannula or mechanical ventilation in intensive care unit and found that SDNN substantially decreased in this population and the decreases in HRV preceded elevations in CRP levels (Hasty et al., 2020). In our study, SDNN continued to be depressed even after acute COVID-19 with a mean duration of 20 weeks after infection.

HRVI is a geometrical measure of time domain HRV and expresses overall HRV, while other time domain parameters such as RR intervals, SDNN, SDANN, RMSDD, pNN50 are statistical measures of HRV. Decreased HRVI reflects autonomic dysfunction, but does not distinguish between specific changes in sympathetic and parasympathetic activity. Major advantage of HRVI is that it is less affected by noise and artifacts (Hämmerle et al., 2020). The prognostic value of HRVI for predicting all-cause mortality and arrhythmic events is robust in heart failure and atrial fibrillation (Hämmerle et al., 2020; Wijbenga et al., 1998). We showed in the present study that HRVI was significantly depressed in patients after recovery from COVID-19 and further reflects some degree of autonomic dysfunction in this group of patients.

Although time and frequency domain parameters of HRV measures variability on various time scales, nonlinear HRV analysis has been adopted for better understanding the characteristics and complexity of the beat to beat variability. Nonlinear methods in HRV do not tell the amplitude of the variability, but rather the qualitative characteristics of dynamics of the signal (Aubert et al., 2003). Nonlinear methods used in our study were DFA, RPA, and entropy-based analysis. The reduction in α1 and α2 in DFA observed in recovered COVID-19 patients may suggest the loss or disorganization of the properties of short-term and long-term fractal correlations of the heart rate dynamics toward more random dynamics. RPA including Lmax, Lmean, DET, REC, and Shannon entropy showed significantly higher values for patients with a history of COVID-19 than in those without a history of COVID-19, translating into less variation and therefore higher autonomic physiological impairment. Entropy-based analysis also measures the complexity or irregularity of HRV. Several different measures of entropy have been used to quantify the heart rate dynamics. Approximate entropy quantifies the unpredictability of fluctuations in a time series. The sample entropy is a measure similar to the approximate entropy. The accuracy of sample entropy is higher than approximate entropy, thus reducing bias (Richman & Moorman, 2000). In our study, both approximate and sample entropies were depressed in patients with a past history of COVID. This shows low variability in this group of patients.

The pathophysiological mechanisms through which COVID-19 causes autonomic dysfunction after full recovery remain unclear. However, there are several pathophysiological mechanisms that may be associated with decreased HRV after COVID-19. Firstly, various conditions such as pain, fever, anorexia, nocturnal sweating, emotional stress, prolonged bed rest, or sleep disorders that occur in clusters at different stages of the COVID-19 may activate sympathetic nervous system and these negative effects may take longer to diminish. Secondly, SARS-CoV-2 could invade the brainstem and alter functions of medullary centers resulting in increased central sympathetic discharges. Recent research on SARS-CoV-2 has shown that this virus can invade tissues by binding to the angiotensin-converting enzyme 2 (ACE2) receptor on glial cells, neurons, and brain microvascular endothelial cells (Xia & Lazartigues, 2010). The olfactory nerve and transsynaptic transfer are also potential pathways for entry of the SARS-CoV-2 into the brain, as the olfactory epithelium, neuronal membrane, and cytoplasm all highly express ACE2 (Xiao et al., 2013). The entry proteins ACE2 may have a role in the binding of the virus to the olfactory nerves followed by travel into the brain (Barrantes, 2020; Stefano et al., 2020). Thirdly, SARS-CoV-2 might infect and destroy, via toxin-mediated or immune-mediated, extracardiac postganglionic neurons of sympathetic nervous system, secondarily increasing cardiac sympathetic discharge. However, there is no published evidence supporting this mechanism.

This study has several limitations. Sample size of the study population was rather small. We did not include those who healed after critically ill during acute phase with an intention to better represent the majority of the post-COVID-19 population because a recent report from WHO says that 80% of infections are mild or asymptomatic with no mortality, 15% are severe disease with no mortality and 5% are critical disease (Madjid et al., 2020). Graphic representations could not be given. Long-term follow-ups are lacking and certainly needed to determine whether lower HRV persists or not in this patient population.

In conclusion, the present data confirm that abnormalities in the three domains of HRV which indicate some degree of dysfunction in cardiac autonomic nervous system are present in subjects with a past history of COVID-19. In addition, our findings indicate an impairment of cardiac parasympathetic function in this population. HRV may be a useful tool to investigate the development of preclinical autonomic dysfunction in these patients. We believe that this preliminary research can serve a starting point for future research in this direction.