There was no involvement of patient and public in this paper.

Eligibility of studies was based on the use of any indices of HRV to access the non-cardiac, physiological control of heart rate in patients with COPD. Studies not using the conventional measuring techniques for HRV (e.g., tilting technique) based on the guidelines of the European Society of cardiology and the North American Society of Pacing and Electrophysiology were excluded (14). Also, studies involving pharmacological interventions previously reported to influence HRV were excluded. Thus, interventional studies/trials assessing the effect of pharmacological or non-pharmacological treatments on COPD or HRV in COPD patients were excluded.

The title and abstract screening of studies was independently performed by two of the authors (JS and TO). To identify studies based on the inclusion criteria, two authors screen the full-texts articles. All conflicts were resolved via online meetings mediated by a third author (AMA).

The evaluation of the quality of the methods employed in each of the included studies was done by two of the authors (SM and NA) using the modified Newcastle-Ottawa Scale (NOS) (15). The modified NOS included seven domains, each of which was assessed from 0 (high risk of bias) to 3 (low risk of bias), and we averaged the domain scores to get a score between 0 and 3, where a higher number reflects a lower risk of bias. Any discrepancies in the quality rating were settled by discussion with a third author.

According to predefined criteria, we extracted the title, location, design, aims and summary findings of included studies into a predefined table. Further, we extracted the sample sizes (of COPD and healthy controls), gender, age, FEV₁% of COPD patients, the ECG recording equipment use, length of ECG analyzed, sampling frequency of ECG, software used for HRV analysis, HRV indices computed and HRV indices that were reported to be influenced by COPD (Table 1). Meta-analysis was performed to compute the standardized mean differences (SMD) with the 95% confidence interval (CI) of the HRV indices reported. Forest plots of SMDs between HRV of COPD patients and matched controls were computed using the Metan algorithm in Stata/SE16. SMD was used as a measure of the effect size to correct for the heterogeneity in the studies design and analysis (duration of ECG or number of NN-interval analyzed). Random or fixed effect model was used according to the between-study heterogeneity measured by the I² statistic. According to Hopkins et al., SMDs are interpreted as trivial (<0.2), small (0.2–0.6), moderate (0.6–1.2), large (1.2–2.0), very large (2.0–4.0) and extremely large (>4.0) based on the magnitude of differences observed between the groups (39). Where HRV indices were reported as median and interquartile ranges, conversion to mean ± standard deviation was performed as described by Hozo et al. (40). Sensitivity test using the “meta, leave-one-out” algorithm on Stata was performed where computed studies are more than two to test whether any of the studies within the meta-analyzed set is exaggerating the effect size. We did not perform for analysis involving two studies as it will be impossible to analyze one study if the other was excluded. Also, the definition of the various HRV indices available is presented in Supplementary Table S4.

Included studies involve a total of 839 patients clinically diagnosed as having COPD with sample size ranges between 10 and 77. Majority of the studies were case control (14/26) with 9 (33%) conducted in Brazil (Table 1). Of the 27 observational studies, 27% (7/26) had a high risk of bias, whereas the rest scored more than 1.5 in the risk assessment, indicating a low risk of bias (Supplementary Table S3). More details about the included studies presented in the Supplementary Table S1.

Twenty two of the 26 included studies analyzed time domain indices of HRV in COPD patients including the standard deviation of NN intervals (SDNN), SDNN Index (SDNNI), standard deviation of the average NN intervals for each 5 min segment of a 24 h ECG recording (SDANN), root mean square of successive NN interval (RMSSD), coefficient of variation of NN intervals (i.e., SDNN/mean NN), count and percentage of NN intervals that differ by 50% (NN50 and pNN50). Of the 22 studies, 8 compared the indices in COPD patients with healthy controls and were included in the computation for effect sizes.

A single study reported compared HRV non-linear indices between patients with COPD and healthy controls (20). The authors compared the long-term (α2) and short-term (α1) fractal components of the detrended fluctuation analysis (DFA) which accesses the presence of fractal correlation in the NN intervals. The authors reported a significance difference in the short-term fractal components between the groups with lower α1 of the DFA found in COPD patients compared with healthy controls (0.9 ± 0.18 vs. 1.02 ± 0.18). This translates to impairment in fractal correlation in the dynamics of NN intervals in COPD patients compared with healthy controls.

Physiologically, the autonomic nervous system responds to a stressful situation by increasing sympathetic activity and decreasing parasympathetic activity, leading to a state of alertness (43). Interestingly, common chronic diseases such as cancer (44), obesity and metabolic syndrome (45), heart failure (46), cirrhosis (47), smoking (48), and hypertension (49), are associated with an activation of sympathetic and a decrease in parasympathetic nervous activity. In our study, we found that patients with COPD have a major reduction in all HRV parameters, including the LF band, which indexes the combined actions of the vagal and sympathetic controls of the heart, with a predominance of the sympathetic (50). At the same time, the LF:HF ratio was found to be increased, which indicates a sympathetic dominance (6). Thus, these findings reveal that COPD patients have decreased sympathetic and parasympathetic activity but remains with a predominance of sympathetic activity. In this context, sympathetic hyperactivity in these patients may be explained by several pathophysiological changes due to recurrent hypoxemia, hypercapnia, increased intrathoracic pressure swings due to airway obstruction, increased respiratory effort, systemic inflammation, or the use of betasympathomimetics (35).

However, under acute exacerbation the neural modulation of the cardiac rhythm may be altered, as evidenced in some studies, and interpretation of findings should be in the context of patients’ clinical condition. For instance, Kabbach et al. (51) found higher values of parasympathetic and global HRV, as well as lower values of complexity of the HR signal in AECOPD patients. Similarly, Zamarron et al. (37) revealed higher HF, LF:HF and total HRV power in exacerbated patients compared to stable ones and, Tseng et al. (33) in an emergency department have verified a greater increase in HF and decrease in LF:HF in those patients requiring admission compared to those discharged. Put together, because HRV indices are influenced by varied and diverse physiological pathways, interpretation of results may need to be clinically contextualized.

In terms of reduction in fractal correlates reported in COPD patients, this translates to reduced complexity and increased similarity in the N-N intervals. Physiologically, there is complex interplay between various organ systems in the body in response to an ever-changing environment. These interactions, aimed at maintaining homeostasis translates into higher complexity in physiological time series which may be reduced in patients with chronic disease which may have altered/reduced the ability to maintain a stable state. Indeed, previous studies have shown that reduced complexity in physiological time series may be associated with poorer prognosis in various settings (52–56).

Some studies analyzed link between COPD severity and HRV changes. Incalzi et al. (16) studied patients with stable COPD and found a correlation between 24-h LF: HF and FEV₁ (r = 0.321; p = 0.022). Consistent with these findings, Bedard et al. (18) also observed a correlation between the 24-h LF: HF ratio and FEV₁ (r = 0.342; p = 0.028), which remained significant even after multivariate analysis. These results suggest that greater disease severity (lower FEV₁ values) is accompanied by predominance of vagal modulation, given that a low LF:HF ratio indicates parasympathetic activity dominance (6). This may not be associated with the patients’ clinical condition however, it may be due to the vagal influence on airway narrowing (bronchoconstriction) (51).

Mazzuco et al. (28) showed a significant positive correlation between LF: HF and total gas volume (r = 0.53) and residual volume (r = 0.52), and negative correlation between the HF power and the total gas volume (r = −0.53). The latter association, suggests a relationship between decreased vagal modulation and increased total lung capacity, an index of static hyperinflation which is a marker of disease severity. Corroborating these findings, Corbo et al. (22) found positive associations between parasympathetic HRV parameters (rMSSD and HF) with FRC/TLC ratio. On the other hand, Camillo et al. (8) found that reduced HRV is not necessarily linked to COPD severity, estimated by FEV_1, and BODE index, but rather to the level of physical activity in daily life. Recently, Castello-Simões et al. (38) have demonstrated a dominance of the clinical status (acute or stable) instead of the severity of the disease on brain heart autonomic axis function assessed by HRV indices and heart rate recovery after exercise.

The relationship between COPD and HRV (and autonomic derangement) is further complicated by the increased risk of cardiovascular diseases that have been reported in patients with COPD. Specifically, patients with COPD have been shown to develop arrhythmia (e.g., atrial fibrillation), heart failure, pulmonary circulation, and arterial diseases, as well as ischemic heart disease, all of which may influence the autonomic regulation of the heart and conversely, HRV (57–59). Further, the increased risk of systemic inflammatory response in COPD patients has been linked with arteriosclerosis and myocarditis, which may also derange the autonomic nervous system’s physiological regulation of heart rhythm (60). Indeed, a recent study by Oyelade et al. showed that in patients with a high risk of arrhythmia, heart rate turbulence, which indexes the autonomic regulation of heart rhythm following premature ventricular complexes (PVCs), is more robust and may provide a more accurate alternative to HRV (61).

Reduced HRV has a powerful prognostic significance for cardiac events, including angina pectoris, myocardial infarction, coronary heart disease death, and congestive heart failure (62). In populations without known cardiovascular disease, low HRV is associated with a 32–45% increased risk of a first cardiovascular event (63). Similarly, a previous systematic review and meta-analysis showed reduced HRV indices in patients with cirrhosis reported overall HRV which may be linked with poorer prognosis (47). As shown here, COPD patients have an overall decrease in HRV parameters. This result may provide a physiological basis for the increased cardiovascular morbidity and mortality observed in COPD patients (64). In addition, COPD is also independently associated with cardiac arrhythmias, including atrial fibrillation/atrial flutter, non-sustained, and sustained ventricular tachycardia (65). This is another evidence suggesting that COPD patients are more susceptible to electrical disturbances and abnormalities of heart rhythm, which may predispose cardiac autonomic dysfunction. Additionally, HRV may have valuable prognostic value for other important COPD-related events, such as the risk of hospitalizations for exacerbation (33) and sudden death (66).

Considering the multifactorial basis of cardiac autonomic dysfunction, and the complex aetiology and clinical characteristics of COPD with multiple phenotypes and endotypes (e.g., frequent exacerbator, non-hyperinflator, emphysematous, overlap comorbidities, physical frailty) (67), the translation of the findings into real practice should be cautious. Traditionally, studies have established the cardiac autonomic profile at rest. However, a dynamic assessment and a reduced responsiveness to a stressor stimulus (physical exercise, postural change, and respiratory manoeuvres) may provide in-depth information and characterize numerous pathophysiological states.

This review has various limitations some of which are inherent to meta-analysis (68). Firstly, most of the included studies were cross-sectional, and the different measurement conditions of HRV parameters as well as heterogeneity of the populations studied possibly led to high inter- and intra-operator variability which is reflected in the high between-studies heterogeneity presented in the forest plots. However, the sensitivity test based on leave-one-out analysis showed not exaggerated effect. Further, while there are similarities between the study populations, they are not clinically identical. Also, although we found that HRV is reduced in COPD according to most studies, whether COPD severity and prognosis is linked with magnitude of HRV changes is still unclear. In addition, most of the studies determine severity criteria based on spirometry measurements, which can be a limitation for clinical applicability. An interesting perspective for future studies would be to associate HRV changes with the ABCD classification proposed by GOLD, which is based on information about the burden of symptoms and risk of exacerbation and has greater clinical applicability (69). In addition, the health status of controls was not sufficiently detailed in most of the studies, which may have influenced HRV. This may also have attenuated the magnitude of differences in HRV parameters when comparing COPD patients and controls. However, most of the included studies were homogeneous according to funnel plots and we computed the standardized mean difference (SMD) which is robust to some of these differences. Additionally, the final number of patients included in the meta-analysis is low, which may make it difficult to generalize the findings. Finally, the inclusion criteria allowed only non-classical HRV measurement algorithms, which may limit the scope of the results. Currently, wearable systems have been used to facilitate the capture of biosignals in clinical practice (60). However, the validity of these methods still needs to be extensively proven. In this review, included studies used mixture of analysis conditions, whereby some studies conducted HRV assessments in exercise conditions (17, 21) during non-invasive ventilation stimulus (19), and autonomic modulation (8), which deserves attention when generalizing these findings.

As a final point, it is important to highlight that in 2020 Catai et al. (70) have published a 30-item checklist to HRV procedures as crucial to dissemination of the HRV findings. The authors described the most common shortcomings and recommendations for solutions in HRV, which included the need for complete characterization of individuals and preparation for signal recordings, environmental condition, data collection method and standard devices, recording length and software, data processing and analysis details. Although it was not the focus of our study, it remains as an overall limitation in the interpretation of HRV studies and deserves to be considered in future studies for HRV to become useful in monitoring COPD patients.