Cardiac Autonomic Dysfunction and Incidence of de novo Atrial Fibrillation: Heart Rate Variability vs. Heart Rate Complexity

Background The REACT DX registry evaluates standard therapies to episodes of long-lasting atrial tachyarrhythmias and assesses the quality of sensing and stability of the lead and the implantable cardioverter-defibrillator (ICD) (BIOTRONIK Lumax VR-T DX and successors) over at least a 1-year follow-up period. Objective To study the association between the risk of de novo device-detected atrial fibrillation (AF), the autonomic perturbations before the onset of paroxysmal AF and a 7-days heart rate variability (7dHRV) 1 month after ICD implantation. Methods The registry consists of 234 patients implanted with an ICD, including 10 with de novo long-lasting atrial tachyarrhythmias with no prior history of AF. The patients were matched via the propensity-score methodology as well as for properties directly influencing the ECGs recorded using GE CardioMem CM 3000. Heart rate variability (HRV) analysis was performed using standard parameters from time- and frequency-domains, and from non-linear dynamics. Results No linear HRV was associated with an increased risk of AF (p = n.s.). The only significant approach was derived from symbolic dynamics with the parameter “forbidden words” which distinguished both groups on all 7 days of measurements (p < 0.05), thereby quantifying the heart rate complexity (HRC) as drastically lower in the de novo AF group. Conclusion Cardiac autonomic dysfunction denoted by low HRC may be associated with higher AF incidence. For patients with mild to moderate heart failure, standard HRV parameters are not appropriate to quantify cardiac autonomic perturbations before the onset of AF. Further studies are needed to determine the individual risk for AF that would enable interventions to restore autonomic balance in the general population.


INTRODUCTION
Atrial fibrillation (AF), the most common clinical arrhythmia, is associated with increased risk of stroke, heart failure (HF) and possibly dementia (Mozaffarian et al., 2016). There is growing evidence that cardiovascular disease (CVD) risk factors and CVD itself explains only 50% of AF occurrences (Agarwal et al., 2017), as the understanding of AF pathophysiology is still incompletely understood (Benjamin et al., 2009). Cardiac autonomic dysfunction has long been suspected in the development of AF (Coumel, 1994), therefore the direct observation of heart rate variability (HRV) immediately preceding the onset of arrhythmia will likely permit the documentation of the underlying mechanism(s). However, the role of the autonomic nervous system should be taken into consideration for optimal individual antiarrhythmic treatment.
Twenty years ago, the Task Force of the European Society of Cardiology (ESC) and the North American Society of Pacing and Electrophysiology published the HRV standards of measurement (Malik et al., 1996). The contemporary literature has a plethora of studies applying HRV methodologies and their successes in clinical applications. However, the identification of high-risk patients has been rather limited. Recently, the e-Cardiology ESC Working Group and the European Heart Rhythm Association co-endorsed by the Asia Pacific Heart Rhythm Society wrote a joint position statement about advances in HRV signal analysis (Sassi et al., 2015). They presented a critical review of newly developed HRV methodologies developed after publication of the initial Task Force HRV overview  and their applications in different physiological and clinical studies.
In a letter  in response to this review, we provided several potential explanations for the limited success that merits additional attention: (1) the importance of monitoring respiration for the interpretation of standard HRV analysis, (2) the need to address its complexities using improved signal processing method (Benjamin et al., 2009;Wessel et al., 2009Wessel et al., , 2016Sidorenko et al., 2016) for patients with mild to moderate HF, standard HRV parameters may not be suited due to the alternating sinus rhythm phenomena (Voss et al., 2006;Costa et al., 2017). Therefore, non-linear methods for the description of heart rate complexity may be more appropriate (Kurths et al., 1995;Voss et al., 1996;Wessel et al., 2000;Brindle et al., 2016).
To capture the parameters from these domains, it is necessary to focus on their physiological interpretation. Spontaneous fluctuations of cardiovascular signals were already described more than a 100 years ago (Ludwig, 1847;Koepchen and Thurau, 1959;Wolf et al., 1978) with fluctuations in both heart rate and blood pressure representing oscillations around fixed values and expressing several influences, e.g., respiration and different self-regulating rhythms.
Short-term heart rate regulation is mainly accomplished by neural sympathetic-and parasympathetic-mediated cardiac baroreflexes and peripheral vessel resistance whereas long-term regulation is achieved by hormonal pathways as well as other systems like the renin-angiotensin-system (Berntson et al., 1997). HRV measurements have proven to be independent predictors of sudden cardiac death after acute myocardial infarction, chronic HF or dilated cardiomyopathy (Kleiger et al., 1987;Malik et al., 1996;Tsuji et al., 1996;Szabó et al., 1997). Moreover, it has been shown that short-term HRV analysis yields a prognostic value in risk stratification independent of clinical and functional variables (La Rovere et al., 2003). However, the underlying regulatory mechanisms are still poorly understood. Further work includes the relevance of premature ventricular ectopic beats, which are associated with an increased risk of sudden cardiac death, as well as their sophisticated neuro-cardio-respiratory interactions (Schmidt et al., 1999;Schulte-Frohlinde et al., 2001, and different means to analyze HRV as stationary epochs (Bernaola-Galván et al., 2001;Camargo et al., 2014).
The main objective of the REACT DX registry was to evaluate different standard therapies to episodes of long-lasting devicedetected atrial tachyarrhythmias (duration >6 min at a frequency >190 beats per minute) and assess the quality of sensing and stability of the lead and implantable cardioverter-defibrillator (ICD) (BIOTRONIK Lumax VR-T DX and successor) over the pre-defined 1-year follow-up period. In addition, cardiac autonomic perturbations before the occurrence of de novo AF were investigated. Therefore, the 7-day continuous ECGs recordings 1 month after ICD implantation were analyzed in all patients. Based on standard time-and frequency-domain parameters as well as from non-linear dynamics (Wessel et al., 2007), the study also analyzed the association between the risk of de novo AF and 7-day HRV (7dHRV).

Cohort
Between November 2012 and June 2016, the multicentre REACT DX registry enrolled 234 patients before or ≤1 month after ICD implantation from 14 sites in Germany. Follow-up duration was at least 12 months and up to 24 months (18.8 ± 4.9 months).
The main inclusion criteria were (i) indications and contraindications for ICD implantation according to national and international guidelines, (ii) be available for follow-up visits on a regular basis at an approved investigational centre, (iii) have a first implantation of an ICD-system, (iv) sign the informed consent form, (v) implantation of a Lumax 740 VR-T DX or successor and (vi) have sufficient coverage of mobile phone network.
The main exclusion criteria were (i) presence of permanent atrial tachyarrhythmia, (ii) indication for cardiac resynchronization therapy, (iii) life expectancy of less than 6 months, (iv) expected cardiac surgery within 6 months after enrolment, (v) age less than 18 years and (vi) enrolment in another cardiac clinical investigation.
Device-detected episodes of atrial tachyarrhythmia were evaluated by the investigators and by a clinical event committee consisting of two experienced physicians.
The REACT DX registry was approved by the ethic committee. All patients signed informed consent. The registry was conducted according to GCP and the Declaration of Helsinki, and was registered in Deutschen Register Klinischer Studien (registration number DRKS00010898).

Pre-processing
The analysis of HRV is often difficult due to many artifacts of the arrhythmias signals. While occasional ectopic beats are treated successfully by most pre-processing methods, more complex arrhythmias or arrhythmias which are similar to normal fluctuations, may remain undetected. We therefore developed a method for data pre-processing which is described in detail elsewhere (Wessel et al., 2007).
A detailed motivation, physiological introduction and the purpose of the analysis of HRV analysis can be reviewed in the Task Force HRV (Malik et al., 1996). It is important TABLE 1 | Parameters available for use in matching via the propensity score methodology: Parameters indicated as: "Fixed matching category" are identical for each matched pair, "Allowed for propensity score" were included into the regression to calculate the propensity score which was used to create the match, "Completeness" indicates the availability of this parameter for the cohort (Checkmarks/crosses highlight accepted/rejected parameters).

Parameter
Fixed matching category Allowed for propensity score Propensity score model coefficients

Completeness
Missing value process to address each artifact with the appropriate tools without influencing the results and the simple exclusion of ventricular premature beats may lead to erroneous HRV parameters (Lippman et al., 1994). An appropriate method to handle most of these problems is through adaptive filtering, which has been described in detail elsewhere (Wessel et al., 2007). The main advantage of this method is the spontaneous adaptation to variability changes, which enables a more reliable removal of artifacts, ventricular premature beats and VPCs. This filtering algorithm consists of three sub-procedures: (i) the removal of obvious recognition errors (ii) the adaptive percent-filter (iii) the adaptive controlling filter It must be mentioned our group has analyzed several thousand human and animal time series (Wessel et al., 2007). In those analyses, small variations of the controlling parameters did not significantly influence the filtering results. A MATLAB implementation of the pre-processing algorithm is available from tocsy.agnld.uni-potsdam.de.

Time-and Frequency-Domain Parameters
Standard methods of HRV analysis include time-and frequencydomain parameters that are linear methods. Time-domain parameters are based on simple statistical methods derived from the RR-intervals as well as the differences between them. Mean heart rate is the simplest parameter, but the standard deviation over the entire time series (sdNN) is the most prominent HRV measure for estimating overall HRV. A list of selected parameters with a short explanation is given in Table A1 and has been described before (Wessel et al., 2007).

Non-linear Analysis of Heart Rate Variability
Heart rate and blood pressure variability reflects the complex interactions of many different controlled loops of the cardiovascular system. In relation to the complexity of the sinus node activity modulation system, a predominantly nonlinear behavior has to be assumed. Thus, the detailed description and classification of dynamic changes using time and frequency measures is often not sufficient. We have previously shown that symbolic dynamic is an efficient approach to analyze the dynamic aspects of HRV (Kurths et al., 1995;Voss et al., 1996). The first step in this analysis is the transformation of the time series into symbol sequences, with symbols given an alphabet letter. Although some detailed information is lost in this process, the broad dynamic behavior can be analyzed (Engbert et al., 1997). Wackerbauer et al. (1994) used the methodology of symbolic dynamics for the analysis of the logistic map where a generic partition is known. However, for physiological time series analysis, a more pragmatic approach is necessary. The transformations into symbols should be chosen in a context-dependent manner. For this reason, we have developed a series of complexity measures based on such context-dependent transformations which have a close connection to physiological phenomena and are relatively easy to interpret (cf. appendix A).
All statistical analyses were performed using the R System. In the case of missing data values, a decision was made on a case to case basis as to if imputation of the missing data was possible or if the variable should be removed. A thorough verification of all patient data was conducted.

Matching
To investigate autonomic perturbations visible through HRV measures, we matched a control group from those that did not register an AF episode during the follow-up period on two levels. To reduce the possible number of subjects, the first matching was based on criteria directly impacting the morphology and interpretation of the ECGs ( Table 1). This specific matching was based on the propensity-score methodology (Imbens and Rubin, 2015) with the score used to represent the propensity of the patient to suffer an AF episode before the end of the follow-up period given the demographic and medical data known at enrolment. This matching was meant to account for any previously available information regarding outcomes. Any differences found during

Propensity Score
A linear logistic regression model was selected to separate those who suffered from AF episodes before implantation and those who did not until at least after the follow-up period. This model only used demographic data and variables from the medical history captured during enrolment (variables implying previous AFs were excluded, Table 1). The remaining variables were verified for completeness and, in the case of missing values, they were included using a k-nearest neighbor approach (Templ et al., 2011). The variables that contained more than 25% missing values were excluded. All ordinal factors were translated to orthogonal polynomial contrasts (Chambers and Hastie, 1991) to allow the incorporation of the order relationship into the model. We performed repeated runs of glmnet cross-validation (Friedman et al., 2010) [α = 0.5, each using a rose resampled (Menardi and Torelli, 2014) variant of the data to account for imbalances] to select the appropriate λ. Finally, we fitted the model to the original, but weighted according to the outcome variable imbalance, dataset at the calculated λ. The selected variables, as used in the model and their coefficients, are presented in Table 1.
The response value of this model is the propensity score used for the secondary matching process, the matching results are presented in Table 2.
All long-term ECGs were recorded for 7 days, 1 month after ICD implantation using GE Healthcare CardioMem CM 3000 devices. HRV analyses were performed separately for each of the 7 days. All HRV parameters described above were calculated on a 24-h basis as well as from windowed analyses. For the latter, the time-and frequency-domain parameters were calculated from successive 5-min windows and then averaged over 24 h. The window length for the non-linear HRV parameter was 30 min (Wessel et al., 2007).

RESULTS
The results presented in Table 3 were derived from the windowed analysis. No HRV parameter from time-and frequency-domain showed a consistent statistically significant differences between the de novo AF and the control group for all 7 days (p = n.s.). The only successful approach was based on symbolic dynamics. "Polvar10" showed significant differences on day 1 and 2 whereas the Shannon entropy H k was significant on 5 days. The parameter "forbidden words, " however, significantly distinguished both groups on all 7 days of measurement as performed independently (p < 0.05; the probability for false statistical results is less than 1 to a billion: <0.05ˆ7). This measure quantifies the heart rate complexity which is drastically lower in the de novo AF group. The averaged effects size over all 7 days of "forbidden words" is almost one, values greater than 0.8 are interpreted as strong effects. I.e., the mean values "forbidden words" for the de novo AF and the control group differ about one standard deviation. To emphasize the strong effect between de novo AF and the control The results of the Mann-Whitney-U-test are presented in the columns Day 1 to 7. No time-and frequency-domain parameters demonstrated consistent significant differences between the de novo AF and the control group. In contrast, the symbolic dynamics parameter Forbword and H k demonstrated significant differences for several days (Bold font: p < 0.05).
Frontiers in Physiology | www.frontiersin.org FIGURE 1 | Heart rate variability (HRV) in long-term ECG recordings. (A) Time series is from a de novo AF patient with a high sympathetic tone. Heart rate is strongly increased (86 bpm) and HRV is clearly decreased (day time sdNN = 12 ms), (B) time series is also from a de novo AF patient but with normal heart rate and HRV (65 bpm, day time sdNN = 98 ms) and (C) time series is from a control patient with normal heart rate and HRV (64 bpm, day time sdNN = 105 ms).
groups, we performed the following Monte-Carlo experiment. From each subject we randomly chose one measurement (1 day out of 7) and compared the "forbidden words" values between both groups. Altogether we have seven one day measurements from 20 subjects (10 de novo AF + 10 controls), which offer 7 20 possibilities for such a random selection. To quantify the effect size, we performed a leave one out classification, i.e., a logistic regression model was trained on 19 subjects and the tested while one of the subjects was left out. This leave one out classification was performed 10000 times. In this way we got an overall accuracy of 0.72 (95% CI: 0.71-0.73, p < 0.001), a sensitivity of 0.64, a specificity of 0.8, a positive predictive value of 0.76 and a negative predictive value of 0.70. All values differ significantly high from a mere random effect. Figure 1 demonstrates the differences between HRV in long-term ECG recordings. It is well known that cardiac autonomic dysfunction, characterized by a high sympathetic tone as observed in Figure 1A and documented via an increased heart rate, is associated with a strongly decreased HRV and a strongly increased risk of AF occurrence (Jons et al., 2010;Nortamo et al., 2018). In this case, not only the HRV but also the heart rate complexity is low (FORBWORD = 12). Figures 1B, 2 show one example of a de novo AF patient with "normal" heart rate and HRV. A decrease in complexity is not readily visible at first glance, however, heart rate complexity is decreased due to the occurrences of alternating sinus rhythm patterns (Costa et al., 2017). The time series in Figures 1C, 3 are from a control patient with normal heart rate, HRV and heart rate complexity. The patterns of the two other patients are undistinguishable using only standard HRV parameters. Only non-linear parameters are able to show a strong decrease in HRC and an increase in FORBWORD due to alternating sinus rhythm patterns present in the times series of Figure 2, which are unremarkable for the time series in Figure 3.

DISCUSSION
The objective of this investigation was to study the association between the risk of de novo AF and 7dHRV based of standard time-and frequency-domain parameters as well as from nonlinear dynamics. Therefore, the long-term ECGs of 20 of the 234 patients with mild to moderate HF included in the REACT DX registry were analyzed. For all patients, the implantation of an ICD device was indicated suggesting all were high-risk patients. The main aim of the present study was to detect the patients with a high risk for de novo AF based on HRV parameters. This was a complex process for several reasons: (i) we did not have the simple task of separating patients from a healthy group as done previously (Costa et al., 2017), (ii) a high amount of ectopy was expected, (iii) patients often had high heart rates and low FIGURE 2 | Heart rate complexity of the de novo AF patient with normal heart rate and HRV from Figure 1B. Panel (A) shows the original time series from Figure 1 demonstrating a high HRV over the day. Panel (B) presents a time breakdown of the gray area from panel (A) and reveals the occurrence of alternating sinus rhythm patterns. The HRC is strongly decreased, therefore FORBWORD is increased to 8 which is comparable to the HRC value presented in Figure 1A. Panel (C) shows the ECG from the gray area from panel (B) showing sinus rhythm was maintained, (cf. Costa et al., 2017).
HRV, (iv) our patient match included different disease states, (v) we did not exclude inflammation and fibrosis, (vi) patient were prescribed different medication regimens and (vii) patients may have had different circadian HRV profiles. Therefore, standard HRV parameters failed to predict cardiac autonomic dysfunction and no parameter from time-and frequency-domain showed significant differences between the de novo AF and the control groups. In contrast, the only significant approach was derived from symbolic dynamics with the parameter "forbidden words" which distinguished both groups on all 7 days of measurements independently (p < 0.05), thereby quantifying the HRC as drastically lower in the de novo AF group.
It is known that cardiac autonomic dysfunction, as evidenced by severely diminished HRV, increases the risk of AF occurrence (Jons et al., 2010). However, these perturbations may not be detected by standard HRV analysis, but appear to be a marker of higher AF incidence (Table 2 and Figure 1). A decreased in heart rate complexity may be accompanied with a decreased in HRV. Researchers (Costa et al., 2017) recently referred to the alternating sinus rhythm pattern phenomenon as "heart rate fragmentation." Applying standard HRV methods to such time series with alternating rhythms may have lead to false or incomplete interpretations as reported by (Bettoni and Zimmermann, 2002) who found a vagal predominance before the onset of AF, and (Vikman et al., 1999) who detected an altered complexity. Others (Costa et al., 2017) also reported abnormal patterns in Poincaré plots and other maps (Woo et al., 1992;Brouwer et al., 1996;Huikuri et al., 1996;Domitrovich and Stein, 2003;Stein et al., 2005Stein et al., , 2008Gladuli et al., 2011;Makowiec et al., 2015). The pathophysiology of alternating sinus rhythm patterns remains to be determined. Different mechanisms may be involved such as the sinus node exit block, a very subtle atrial bigeminy and the sinus node parasystole and perturbations of internal pacemaker "clocks" in the SA node (Costa et al., 2017). Costa et al. (2017) speculated that heart rate fragmentation would be of high interest if it was an initial event leading to arrhythmias such as AF or other tachyarrhythmias. The results of our study confirmed this hypothesis. Costa et al. (2017) raised the additional question whether abnormalities in breathing dynamics could be responsible for the fragmentation. An indirect indication of cardiorespiratory coupling involvement is that deep breathing may lead to alternating patterns in heart rate (Löwit, 1879). Following earlier work on cardio-respiratory synchronization changes with sleep and aging (Bartsch et al., 2012), we  introduced the analysis of cardiorespiratory coordination (CRC) during sleep. They found, by using the advanced analysis technique of the coordigram, not only the occurrence of CRC was significantly more frequent during respiratory sleep disturbances than during normal respiration but also more frequent after these events. A further investigation of CRC (Krause et al., 2017) showed the new phenomenon of heartbeat-initiated inspiration at the end of an apnea, which provides new impulses to current approaches to obstructive sleep apnea characterization (Ivanov et al., 1996;Mietus et al., 2000).
The major limitation of our study is the low number of de novo AF events as only 10 patients were included. AF is a powerful risk factor for stroke, independently increasing risk about five-fold throughout all age groups (Mozaffarian et al., 2016). Because AF is often asymptomatic and likely frequently undetected, the risk of stroke attributed to AF may be substantially underestimated. In the setting of AF, important risk factors for stroke include advancing age, hypertension, HF, diabetes mellitus, previous stroke or TIA, vascular disease and female sex (Schmitt et al., 2009;Mozaffarian et al., 2016). Additional biomarkers such as high levels of troponin and BNP also increase the risk of stroke independent of the well-established clinical characteristics (Mozaffarian et al., 2016). Thus, further studies are needed to determine the individual risk for AF based on HRV, heart rate complexity, clinical characteristics and additional biomarkers.

CONCLUSION
Cardiac autonomic dysfunction denoted by low heart rate complexity may be associated with a higher AF incidence.
For patients with mild to moderate HF, standard HRV parameters are not appropriate to quantify cardiac autonomic perturbations before the first occurrence of AF. Further studies are needed to determine the individual risk for AF that would enable interventions to restore autonomic balance in the general population.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by DRKS00010898. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
NW, TH, RB, and AS designed and directed the project. TH, KR, RB, AS, DW, SK, CK, and NK performed the register. NW, AG, KB, JFK, and JK analyzed the data. NW, JK, KB, AG, and KR wrote the manuscript. All authors contributed to the article and approved the submitted version. FUNDING NW, AG and KB were supported by the DFG grant WE 2834/11-1. NW, JFK and KB were supported by the BMBF grant in ZHerz. JK acknowledges support from Russian Ministry of Science and Education "Digital Biodesign and Personalised Healthcare". The study was sponsored by BIOTRONIK SE & Co. KG (Berlin, Germany).

APPENDIX A Symbolic Dynamics
While comparing different types of symbol transformations, we found that the use of four symbols, as explained in Eq. (1), was appropriate for our analysis. The time series x 1 , x 2 , x 3 ..., x N is transformed into the symbol sequence s 1 , s 2 , s 3 ..., s N , s i ∈ A on the basis of the alphabet A = {0,1,2,3}.
The transformation into symbols refers to three given levels where µ denotes the mean beat-to-beat interval and a is a special parameter that we have chosen 0.05. We tested several values of a from 0.02 to 0.08, however the resulting symbol sequences differed not significantly (Wessel et al., 2007). Because there are several quantities that characterize such symbol strings, we analyzed the frequency distribution of words of length 3 (i.e., substrings which consist of three adjacent symbols leading to a maximum of 64 different words), which is a compromise between retaining important dynamical information and having a robust statistical estimate of the probability distribution. We considered 3 measures of complexity: (i) The Shannon entropy H k calculated from the distribution p of words is the classic measure for the complexity in time series: where W k is the set of all words of length k. Larger values of Shannon entropy refer to higher complexity in the corresponding tachograms and lower values to lower ones. (ii) The "forbidden words" in the distribution of words of length 3 which are the number of words that (almost) never occur (with probability of less than 0.1%). A high number of forbidden words reflects a rather regular behavior in the time series and, if the time series is highly complex in the Shannonian sense, only a few forbidden words are to be found. (iii) The parameter "plvar10" characterizing short phases of low variability from successive symbols of another simplified alphabet B and consisting only of symbols "0" and "1." In our study, "0" stands for a small difference between two successive RR-intervals (the resolution of the defibrillators) whereas "1" represents cases when the difference between two successive RR-intervals exceeds a certain specific limit. s n = 1 : |x n − x n−1 | ≥ 10 msec 0 : |x n − x n−1 | <10 msec.
Words consisting of a unique type of symbols (either all "0" or all "1") were counted. To obtain a statistically robust estimate of the word distribution, we chose words of length six, defining a maximum of 64 different words. "Plvar10" represents the probability of the word "000000" occurrence and thus detects even intermittently decreased HRV. BBI stands for the filtered beat-to-beat intervals (NN-intervals).