Heart Rate Variability Analyses in Parkinson’s Disease: A Systematic Review and Meta-Analysis

Recent evidence suggests that the vagus nerve and autonomic dysfunction play an important role in the pathogenesis of Parkinson’s disease. Using heart rate variability analysis, the autonomic modulation of cardiac activity can be investigated. This meta-analysis aims to assess if analysis of heart rate variability may indicate decreased parasympathetic tone in patients with Parkinson’s disease. The MEDLINE, EMBASE and Cochrane Central databases were searched on 31 December 2020. Studies were included if they: (1) were published in English, (2) analyzed idiopathic Parkinson’s disease and healthy adult controls, and (3) reported at least one frequency- or time-domain heart rate variability analysis parameter, which represents parasympathetic regulation. We included 47 studies with 2772 subjects. Random-effects meta-analyses revealed significantly decreased effect sizes in Parkinson patients for the high-frequency spectral component (HFms2) and the short-term measurement of the root mean square of successive normal-to-normal interval differences (RMSSD). However, heterogeneity was high, and there was evidence for publication bias regarding HFms2. There is some evidence that a more advanced disease leads to an impaired parasympathetic regulation. In conclusion, short-term measurement of RMSSD is a reliable parameter to assess parasympathetically impaired cardiac modulation in Parkinson patients. The measurement should be performed with a predefined respiratory rate.


Introduction
Parkinson's disease (PD) is one of the most common neurodegenerative disorders. As a multisystem disorder, PD is characterized by motor symptoms and a plethora of nonmotor symptoms [1][2][3]. An increasing number of studies indicate that these are at least partly caused by changes in the gut-brain axis [4], dysbiosis with local inflammation [5], and finally alpha-synuclein accumulation in the enteric nervous system [6]. These pathological proteins may spread via the vagus nerve to the central nervous system [7]. The vagus nerve is the main parasympathetic branch of the autonomic nervous system. Additionally, central dopaminergic cell death can favor central inflammatory processes that affect the vagus nerve directly [8]. An altered function of the autonomic nervous system can cause various nonmotor symptoms, accordingly, many PD patients have symptoms of dysautonomia [9]. Frequently, patients complain of hypersalivation, swallowing difficulties, delayed gastric emptying, constipation, or orthostatic hypotension [10].
Assessing parasympathetic regulation in PD is a promising way to assess symptoms and to improve our knowledge of the course of the disease. Similarly, improved dysautonomia testing may be of additional value in differentiating PD and atypical parkinsonian syndromes. Heart rate variability (HRV) analysis is a simple method to estimate the overall tone of the autonomic nervous system [11]. Whereas heart rate quantifies the number of heartbeats per minute, HRV refers to the fluctuation in the time between successive heartbeats. It reflects the sympathetic and parasympathetic modulation of cardiac activity. In cardiology, there are recommendations for HRV analysis for standards of measurement, physiological interpretation and clinical use [12]. Recommendations for clinical use in patients with neurodegenerative diseases do not exist. In the literature, there are many different approaches to studying HRV in PD.
In general, a distinction can be made in (1) frequency-domain indices, (2) time-domain indices, and (3) nonlinear measurements [11]. Frequency-domain indices are used to describe the power spectral density as a function of frequency using mathematical algorithms. Short-term recordings can distinguish high-frequency (HF), low-frequency and very-low-frequency components. Long-term recordings can additionally distinguish ultralow-frequency components. The HF component is generated mainly by parasympathetic modulation [11]. Time-domain indices quantify HRV observed during monitoring periods ranging from shorter than 1 minute to longer than 24 h. To describe parasympathetic regulation, there are two established parameters. The root mean square of successive normal-to-normal interval differences (RMSSD) in ms reflects the beat-to-beat variance, and the percentage of adjacent normal heartbeat intervals that differ from each other by more than 50 ms (pNN50).
In addition to these mentioned statistical measurements, time-domain recordings can also be converted into nonlinear geometric patterns. However, these geometric patterns mainly describe the total variance and not parasympathetic activity [11], and will therefore not be considered in this review. In this systematic review of parasympathetic modulation in PD, we will focus on the HF spectral component, RMSSD, and pNN50. Lower values of HF, RMSSD, and pNN50 are an indication of reduced parasympathetic regulation.
The objective of our study is to systematically review the literature if a frequencyor time-domain analysis of the HRV may indicate a decreased parasympathetic tone in patients with PD. Our secondary aim is to determine the most suitable HRV method for clinical use in PD.

Search Strategy and Eligibility Criteria
We conducted our systematic review and meta-analysis in accordance with the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) [13]. This study does not constitute human subject research. Therefore, there was no need for local ethics committee approval. The search strategy was created in consultation with a medical statistician with expertise in systematic review searching. Three electronic databases (MED-LINE, EMBASE, and Cochrane Central) were screened from inception to 31 December 2020, without any restriction on the year of the study. We searched MEDLINE with the Medical Subject Headings "Parkinson disease" (Mesh) AND "heart rate" (Mesh) and EMBASE with the subject headings "Parkinson disease" AND "heart rate variability" and limited the records to humans, English, and articles. Our search in Cochrane Central was conducted using the keywords "Parkinson" AND "heart rate variability" and allowing word variations. The reference lists of included studies or relevant reviews identified through the search were searched to identify any studies missed in the initial search to ensure literature saturation. Only full-text articles were considered for analysis. References were managed using EndNote (version X8; Clarivate Analytics). Titles and abstracts were screened by one reviewer (K.G.H.) and controlled by a second reviewer (T.P.) if the study was excluded.
Studies needed to be original research. We included studies with an observational design, randomized controlled trials, controlled (nonrandomized) clinical trials, and interrupted time series studies if there was at least one reported data point before the intervention. Review articles, book reviews, short communications, and correspondences were excluded.
Two reviewers assessed the full-text articles for eligibility (K.G.H. and T.P.). Different assessments were resolved by consensus. Articles were included if they: (1) were published in English, (2) reported patients with idiopathic PD as a primary study group and healthy adult controls as a secondary study group, and (3) reported at least one HRV parameter, which represents parasympathetic activity (HF spectral component, RMSSD, or pNN50).

Data Extraction
The data extraction form was developed in Microsoft Excel (version 2016; Microsoft Corporation). The data extraction was performed by one member of the research team (K.G.H.) and checked by a second (T.P.). The following data were extracted: study information (author, title, journal, and year of publication), study characteristics (e.g., study setting), participant characteristics (e.g., age, sex, disease severity, and disease duration), condition information (i.e., data sources, total number of participants, and exclusion criteria), and HRV parameters (HF, RMSSD, and pNN50). Data on participant characteristics and HRV parameters were extracted and, if possible, grouped to achieve better evaluability. The results of substudies were extracted separately in order to allow subsequent subgroup analysis based on different subgroup characteristics. If there were incomplete data reported in the publication, we searched supplementary information and documents to locate missing data.

Statistical Analyses
For each study population, the main characteristics were reported as mean ± SD (standard deviation) for continuous variables and the number (%) for categorical variables. If not reported, the sample mean and SD were synthetized and estimated from the sample size, median, range and/or interquartile range [14]. Random-effects meta-analyses by means of the Hartung-Knapp-Sidik-Jonkman method were conducted. If substudies were reported, the raw mean of each study was calculated to avoid an overestimation of the impact. Data were analyzed separately for each HRV parameter. Hedges' g was used as a measure of standardized mean differences (SMD).
We estimated the significance and degree of heterogeneity in the study results by Cochrane's Q, I 2 -statistics, and Tau 2 -statistics. To evaluate the heterogeneity, we searched for outliers. Additionally, we performed influence analyses by means of the leave-one-out method and GOSH plot analyses.
We conducted random-effects subgroup analyses to assess measurement time-dependent between-subgroup differences. Additionally, subgroup analyses were carried out between normal and deep breathing to assess the impact of respiration on HRV. Deep breathing with a predefined respiratory rate standardized the respiration-driven acceleration and deceleration of the heart rate [15].
We considered studies that performed both normal and deep breathing in the same population to assess comparability.
Random-effects meta-regression analyses were performed to examine the relationship between HRV parameters and relevant covariates. We checked for multicollinearity. We considered predictors to be highly correlated if r > 0.8.
Publication bias was assessed through Egger's Test. A p-value of Egger's test < 0.05 is significant, which means that there is substantial asymmetry in the funnel plot. When Egger's test was significant, the trim-and-fill procedure was used to adjust for publication bias.
We did not perform a detailed quality assessment of each included study. First, because mainly no intervention effect was analyzed and second, because the statement of the initiative to strengthening the reporting of observational studies in epidemiology (STROBE statement) should not be used as a tool to assess the methodological quality of cohort studies [16].

Selection and Study Population
We screened 485 unique citations. Of these, 135 were assessed for eligibility, and 47 studies were included in this systematic review and meta-analysis . For further details on study selection, see Figure 1. The study population includes 2772 subjects. Of these, 1566 had PD, and 1206 were healthy controls. The summarized findings of the population are shown in Table 1.
(STROBE statement) should not be used as a tool to assess the methodological quality of cohort studies [16].

Selection and Study Population
We screened 485 unique citations. Of these, 135 were assessed for eligibility, and 47 studies were included in this systematic review and meta-analysis . For further details on study selection, see Figure 1. The study population includes 2772 subjects. Of these, 1566 had PD, and 1206 were healthy controls. The summarized findings of the population are shown in Table 1.

Study Characteristics
All included studies used HRV analysis to assess autonomic function in PD patients compared to controls. Most studies focused on the measurement of HRV under particular conditions regarding time, duration, and method, or they intended to describe a relation between HRV and the severity of PD. An overview of the included studies can be found in Table A1 in the Appendix A.
The duration of the measurements varied from 1 minute to 24 h. Because there is no uniformly accepted limit, we categorized the study into short-term (shorter than one hour) and long-term (longer than one hour) measurements. Short-term measurements were carried out in 34 studies, with a range from 1-20 min. Long-term measurements were carried out in 16 studies, with a range from 6-24 h. Short-term measurements were mainly conducted during daytime and long-term measurements during day-and nighttime.

Random-Effects Meta-Analyses
We conducted a random-effects meta-analysis for each HRV parameter. The results are shown in Table 2(A), raw analyses. 3.5 ± 0.5 5.5 ± 0.8 −0.59 −1.06 −0.12 0.020 64 0. 28 9 Values are given as mean or as mean ± SD, unless otherwise indicated: (A) raw analyses; (B) heterogeneity analyses; ci.lb: confidence interval lower bound; ci.ub: confidence interval upper bound; HF: high frequency components of the power spectral density; HRV: heart rate variability; k: number of considered studies; nu: normalized units; PD: Parkinson's disease; pNN50: number of normal-to-normal intervals differing by more than 50 ms divided by the total number of normal-to-normal intervals; RMSSD: Root mean square of successive normal-to-normal interval differences; SMD: standardized mean differences.

Figure 2.
Forest plot HFms 2 , random-effects meta-analysis. HFms 2 : high-frequency components of the power spectral den sity in ms 2 ; N: number; SD: standard deviation; CI: confidence interval. Rhombus indicates meta-analytically pooled est mate of the 95% confidence interval. The numbers after the authors indicate the included substudies according to th information provided in Appendix A, Table A1.

Figure 2.
Forest plot HFms 2 , random-effects meta-analysis. HFms 2 : high-frequency components of the power spectral density in ms 2 ; N: number; SD: standard deviation; CI: confidence interval. Rhombus indicates meta-analytically pooled estimate of the 95% confidence interval. The numbers after the authors indicate the included substudies according to the information provided in Appendix A, Table A1.

Between-Study Heterogeneity
High values for I 2 were found for all HRV parameters. We evaluated heterogeneity by performing outlier analysis, influence analysis, and GOSH plot analysis for each HRV parameter. Then, we excluded studies that were identified with a potentially high risk of

Between-Study Heterogeneity
High values for I 2 were found for all HRV parameters. We evaluated heterogeneity by performing outlier analysis, influence analysis, and GOSH plot analysis for each HRV parameter. Then, we excluded studies that were identified with a potentially high risk of bias. Again, a random-effects meta-analysis was conducted. The results are shown in Table 2(B), heterogeneity analyses. Taken together, in PD patients, there were still significantly lower values for HFms 2 . Additionally, there were significant effect sizes regarding RMSSD and pNN50. I 2 decreased considerably. However, there was still a moderate degree of heterogeneity, indicating different patient populations, interventions, or measurement methods.

Meta-Regression
First, we considered patient age, sex, disease duration and severity of disease d mined by Hoehn and Yahr stage, Unified Parkinson's Disease Rating Scale (UPDRS) and UPDRSIII subscore separately. A longer disease duration significantly contribut lower HFms 2 (p = 0.033), but only in the long-term measurement. A higher Hoehn Yahr stage was associated with lower RMSSD (p = 0.0274), but only in the short-term m urement.

Meta-Regression
First, we considered patient age, sex, disease duration and severity of disease determined by Hoehn and Yahr stage, Unified Parkinson's Disease Rating Scale (UPDRS) score and UPDRSIII subscore separately. A longer disease duration significantly contributed to lower HFms 2 (p = 0.033), but only in the long-term measurement. A higher Hoehn and Yahr stage was associated with lower RMSSD (p = 0.0274), but only in the short-term measurement.
Second, we conducted multiple meta-regression analyses. The UPDRS score and UPDRSIII subscore are with a substantial degree of accuracy significantly linearly predicted from the Hoehn and Yahr stage (r p = 0.92 and r p = 0.81, respectively). Accordingly, as covariates, we included patient age, sex, disease duration and Hoehn and Yahr stage. Multiple meta-regressions revealed no significant contribution of any covariate to one of the HRV parameters.

Discussion
PD patients exhibited a strong decrease in parasympathetically modulated HRV parameters compared to healthy controls. In our meta-analysis, PD was associated with lower values of HFms 2 and lower values of RMSSD during short-term measurements. No decrease was seen regarding HFnu. The normalized unit of HF is the relative value of HF in proportion to the total power minus the very-low-frequency component. Total power includes all frequencies and depends on the regulation of both branches of the autonomic nervous system. Therefore, normalization tends to minimize the parasympathetic influence. Accordingly, it must be emphasized at this point that a determination of HFnu does not allow any valid statement on parasympathetic modulation, because it reflects the overall autonomic tone.
As shown in the meta-regression analyses, there is some evidence that more advanced PD leads to an increasing reduction in these parasympathetically modulated HRV parameters. This outcome corresponds to the assumption of increasing damage to the vagus nerve in the course of the disease [61,65]. Especially regarding HFms 2 , disease duration significantly contributes to the effect size. Additionally, the Hoehn and Yahr stage significantly contributed to the effect size of RMSSD. However, these correlations could not be consistently ascertained for all HRV parameters. Clinical parameters regarding the duration and severity of the disease have not been reported in many studies. The UPDRSIII subscore was only reported in 16 of 47 studies. Only two studies used the MDS-UPDRSIII subscore [29,61], although it was published in 2007 and should be used preferably [66]. Insufficient data reporting precluded studies from the meta-regression model.
It is known that age and sex have an important impact on HRV [67][68][69][70][71]. Aging causes neuronal and structural modifications of the cardiorespiratory system leading to a decreased vagal tone [68,70]. The influence of patient age could not be confirmed by our results. Generally, it is assumed that females show greater vagal tone than males [71]. This differential autonomic tone indicates age-and sex-related predisposition to cardiovascular diseases. We could not find evidence that in PD, HRV parameters depended significantly on sex. This result may be due to progressive sex-independent vagal damage in the course of PD.
Heterogeneity of the studies was high. Besides a limited number of studies, this indicates different patient populations and diverse study settings. Especially regarding frequency-domain measurements, many studies had to be classified as outliers or influence studies. High heterogeneity makes valid regression analysis difficult and complicates the investigation of relevant covariates regarding disease severity. This issue could be improved through a standardized basic measurement in addition to a study-specific measurement and through reporting of defined clinical parameters. This would allow for improved subgroup analysis.
The importance of a subgroup analysis becomes clear when considering the measurement time. Considering measurement time, it is evident that RMSSD is significantly reduced in PD patients only during short-term measurements. Theoretically, a 24-h recording seems to be advantageous because the total variance of HRV increases. However, the prerequisite for this would be that the mechanisms responsible for heart rate modulation remain unchanged during the recording period to ensure a largely stable heart rate and respiratory rate. Otherwise, the results of HRV analysis may be more due to external influences than to autonomic regulation. In particular, physiological mechanisms of heart rate modulation cannot be considered stationary over a 24-h period [72]. Therefore, short-term measurements under stationary conditions, especially regarding physical activity, position, and temperature, seem to be advantageous, which is supported by our meta-analysis.
To ensure literature saturation patients with a leucine-rich repeat kinase 2 mutation (LRRK2) associated PD were not primarily excluded. Additionally, most of the patients were not explicitly genetically tested and therefore it cannot be excluded that LRRK2-associated PD patients were included in other study populations. Autonomic dysfunctions are associated with genetic forms of synucleinopathies like LRRK2-associated PD [73]. Within this meta-analysis we identified two studies, which examine HRV in LRRK2-associated PD patients [26,60]. These studies showed higher parasympathetically modulated HRV parameters in LRRK2-associated PD patients compared to idiopathic PD. Therefore, it can be assumed that the described effect sizes are rather underestimated.
It is well known that respiration has a major impact on heart rate [74]. The heart rate increases during inspiration and decreases during expiration, which is called respiratory sinus arrhythmia [15]. Measurements conducted with a predefined deep breathing rate enable HRV analyses to reflect the vagal tone rather than differences in respiration [72]. The obtained results of our subgroup analyses regarding normal and deep breathing should be interpreted with caution. The results are based on the values of a very limited number of studies. Especially regarding RMSSD and pNN50, the results are obtained from just one study [27]. Therefore, it is strongly recommended to perform HRV analysis with a predefined respiratory rate. Otherwise, the significance of the analysis is very limited. This should be considered in further studies.
In addition to heterogeneity and meta-regression analyses, the results have to be evaluated under consideration of publication bias. Rather unexpectedly, there was no significant adjusted effect size regarding HFms 2 . However, a subgroup analyses revealed that there is evidence of publication bias only in the short-term measurement of HFms 2 , and not in the long-term measurement.

Conclusions
HRV analysis is an easy tool to measure vagal influence on the heart rate, and indirectly, to assess autonomic dysfunction in PD. PD patients showed decreased parasympathetically modulated HRV parameters compared to healthy controls. This finding is in line with the assumption of vagal damage in the course of the disease. After considering heterogeneity analysis, subgroup analysis and evaluation of publication bias, we recommend establishing a short-term measurement of RMSSD as a basic measurement of HRV in future studies. At the very least, this parameter should be reported as a basic measurement in addition to study-specific measurements. To standardize the impact of respiration, the measurement should be performed with a predefined respiratory rate of around six respiratory cycles per minute. To enable further assessment of disease progression and in the attempt to distinguish PD and atypical parkinsonian syndromes, we recommend specifying defined clinical parameters of the patient collective. In particular, patient age, sex, disease duration and severity of disease determined by Hoehn and Yahr stage are easy to collect even for non-neurologists. A specification of the MDS-UPDRSIII would further improve the quality of future studies.

Data Availability Statement:
The data presented in this study are available from the corresponding author on reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest. Values are given as the mean, unless otherwise indicated: #: values are given as the median; na: nonreported data; *: not considered in further analyses because data are not plausible or reported in rarely used units; HC: healthy controls; HF: high-frequency components of the power spectral density; nu: normalized units; PD: Parkinson's disease; pNN50: number of normal-to-normal intervals differing by more than 50 ms divided by the total number of normal-to-normal intervals; RMSSD: root mean square of successive normal-to-normal interval difference.