Associations between heart rate asymmetry expression and asymmetric detrended fluctuation analysis results

Abstract The relation between recently established asymmetry in Asymmetric Detrended Fluctuation Analysis (ADFA) and Heart Rate Asymmetry is studied. It is found that the ADFA asymmetric exponents are related both to the overall variability and to its asymmetric components at all studied time scales. We find that the asymmetry in scaling exponents, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha ^{+}<\alpha ^{-}$$\end{document}α+<α- is associated with both variance-based and runs-based types of asymmetry. This observation suggests that the physiological mechanisms of both types are similar, even though their origins and mathematical methods are very different. Graphical abstract The graphical abstract demonstrates strong, nonlinear association between the expression of Heart Rate Asymmetry measured using relative descriptors and the Asymmetric Detrended Fluctuation Analysis results. It is clear that there is a strong relation between the two theoretically disparate approaches to signal analysis. The technique to demonstrate the association is loess fit.


Introduction
Asymmetry in the RR intervals time series has recently found interest in the heart rate variability research [1][2][3][4][5]. Starting from the most generic approaches to time-irreversibility, which in fact spans an enormous number of physiological phenomena and heart rate variability (HRV) measures [5][6][7] through specialized predictive methods like Phase Rectified Signal Averaging (PRSA) [8,9] or Heart Rate Turbulence (HRT) [10][11][12] down to phenomena like Heart Rate Asymmetry (HRA) along with its mathematical tooling.
A recent addition to the literature is the Asymmetric Detrended Fluctuation Analysis, which was developed mainly for studying asymmetric correlations in economic time series [13,14], but when applied to the RR intervals time series, revealed an asymmetric physiological effect [15,16].
HRA includes asymmetric effects defined and established in the variance of the RR intervals time series, in its structure and complexity. In all of these areas asymmetry is prevalent, consistent and unidirectional [1,2,[17][18][19][20]. Thus, we hypothesize that HRA and the above-mentioned ADFA results should be connected, and our aim in this paper is to relate the + < − asymmetric effect established in ADFA to HRA.

Asymmetric detrended fluctuation analysis
Detrended fluctuation analysis is one of the most often used methods for analyzing the time series of RR intervals. The main information it provides is the scaling properties of noise left over after detrending the time series. The details of the method may be found in [21][22][23], and they will be reviewed below very quickly so as to establish the notation.
Let us define the RR intervals time series as the distance between successive R-waves in an electrocardiogram [24,25] in the following way: Let us also define a derivative, summed and mean-subtracted time series: where RR stands for the mean of the whole time series and y(k) defines the so-called box of length n within which a trend is found and subtracted. Since our aim is to identify rising and falling trends we only use first order (linear) polynomials where y n (k) is a line fit with slope a n and intercept b n to the specific box and k is the error term we are interested in studying. A mean-square root function is defined as: and calculated over all the scales which are available in the studied time series -in practice n changes from 4 to N/4, where N is the length of the time series. The values of log 10 (F(n)) are plotted against log 10 (n) and if the resulting plot follows a straight line, the existence of power law is concluded in the scaling of the mean-square root function, i.e., where is the scaling exponent. The values of this exponent are interpreted as signifying the presence of negative long-range correlations ( 0 < < 0.5 ), white noise ( = 0.5 ), positive long-range correlation ( 0.5 < < 1 ), 1/f noise ( = 1 ), long-range correlations not following the power law ( 1 < < 1.5 ) or the consistency of the detrended time series with random walk ( = 1.5 ) [21][22][23].
The fact that local trends in the boxes can be linear makes it possible to define two different mean-square root functions, depending on the sign of the tangent of the fitted line.
Defining M as the overall number of segments and M + n as the number of segments in which the trend is increasing (or, using Eq. (3), a n > 0 ) we define: and correspondingly for the decreasing trends: where the outer summation (over j) goes over all boxes, (±) selects boxes with increasing (Eq. (6)) or decreasing trends (Eq. (7)) exclusively (compare [14,16]). If the above functions are presented on a doubly logarithmic scale with n, two scaling exponents: may be defined, provided that the dependence is linear.
In [13] the local version of ADFA was developed in which the scaling exponents ± are calculated in a window moving along the analyzed time series.
In [15] we applied the local version of ADFA to the time series of RR intervals (moving window length 100) and the result was that there is a highly statistically significant asymmetric effect with + < − . In [16] we systematically analyzed this effect and found that it was present in both global and local versions of ADFA with the use of windows of length 100 through 1000 in 30-min ECG recordings, but it was much weaker in the global version.

Variance-based heart rate variability descriptors
The variance of the time series (1) is defined in the following way [1,24,26]: where N is the length of the time series. Variance SDNN 2 can be partitioned into short-term and long-term variability in the following way [1,24,26,27]: The reasons for calling SD1 2 and SD2 2 the short-term and long-term variability and the details on their calculations are explained in detail in [1,24,26].

Variance-based heart rate asymmetry descriptors
The source of variance-based HRA descriptors is the decomposition of variance-based HRV descriptors, such as SDNN 2 , SD1 2 and SD2 2 into parts which only depend on decelerations or accelerations.
In [1] it was shown that SD1 2 can be partitioned into two parts dependent separately on decelerations and accelerations in the following way: In [17] it was demonstrated that long-term variability may be partitioned in the following way: and the full variance may be partitioned in the following way: For numerical and algorithmic details of the above see [1,17].
The respective parts of variance can be normalized in order to minimize interpersonal variability [1,17]. For short-term variance: and there is: For long-term variance: where: And finally, for total variance; where: The above descriptors, when applied to the RR intervals time series, reveal a strong asymmetry of this object. First of all, the contribution of heart rate decelerations to short-term variability is greater than that of accelerations, i.e., SD1 2 d > SD1 2 a . In long-term variability this is reversed with SD2 2 d < SD2 2 a , and this is also true in the case of total variability with SDNN 2 d < SDNN 2 a . If the normalized contributions defined above (14), (16), (19) are taken into account, the asymmetry in short-term, long-term and total variability may be respectively expressed as C1 d > 0.5 , C2 d < 0.5 and C d < 0.5.

The runs method
A run is an uninterrupted sequence of RR intervals which constantly shortens (heart rate accelerates) or constantly lengthens (heart rate decelerates) or constantly does not change, and which is preceded and followed by a different type of run. Figure 1 shows the partitioning of a segment of an RR intervals time series into disjoint accelerating and decelerating runs. It can be easily noted that this partitioning is unambiguous [2]. A detailed definition of runs may be found in [2].

Number of runs
The most natural descriptor is the number of runs of a specific type. So, by DR i we will denote the number of deceleration runs of length i, and AR i will mean the number of acceleration runs of length i.

Entropic descriptors of runs of decelerations and accelerations
In [2] the Shannon entropy [28] associated with the distribution of decelerating and accelerating runs was partitioned into parts depending only on decelerations or only on accelerations in the following way (dropping the device-dependent neutral runs): The mathematical details of the above formulas may be found in [2]. In [2,18] it was found that the runs of accelerations are longer in terms of the number of beats than the runs of decelerations. In [2] it was found that H DR < H AR . Both these effects are highly statistically significant.
The runs method turned out to have a significant predictive value for long-term survival in patients after myocardial infarction [18] and in patients who underwent clinically indicated exercise test [19].
Runs have also been found to be very useful in studying and predicting sleep apnea [20,29] as well as selecting patients who will respond to the proper treatment of obstructive sleep apnea. It has also been applied to the diagnosis of late sepsis in neonates [30].

Methods
We used 388 stationary 30-min ECG recordings from healthy young subjects, age range 20-40 years, 233 women. The study was performed at rest in the supine position, and the subjects were kept quiet in a neutral environment. The subjects were allowed to breathe spontaneously during the whole study. The 30-min recording was taken after a preceding 15-min period used for cardiovascular adaptation. The ECG curves were sampled at 1600 Hz with the use of the analog-digital converter (Porti 5, TMSI, Holland). The libRASCH/RASCHlab (v. 0.6.1, www. libra sh. org, Raphael Schneider, Germany) [31] software was used for post-processing and automatic classification into beats of sinus, ventricular and supraventricular origin as well as artifacts. The automatic classification was reviewed by a trained technician who corrected any wrong classifications of the beats. To obtain the asymmetric descriptors from the annotated RR intervals time series we used in-house, free GPL3 software written in Python, HRAexplorer, which can be reviewed and downloaded at https:// github. com/ jarop is/ HRAEx plorer. An interactive online version of this software in the R programming language may be found at https:// hraex plorer. com/. The RR intervals time series were carefully filtered for each technique -the specifics of dealing with ectopic beats are described in [25] and [2]. The above-mentioned software uses all these filtering techniques.
Since both theory and the Shapiro-Wilk test reject the normality of all HRA descriptors, the non-parametric paired Wilcoxon test was used to establish asymmetric relations. The binomial test was used to establish the departure from 0.5 of recordings exhibiting HRA, which would be the case in symmetric (e.g., shuffled) data.
In the present paper we used the local version of ADFA with jumping window of length 100 beats, since in a previous study [16] it was the shortest window in which asymmetry was clear and consistent with longer windows. By jumping windows we mean disjoint windows of length 100 fully covering the analyzed recording. This can be contrasted with a sliding window which means a window sliding along the recording, moving by either one beat or by a time unit -for details see [16,32]. To calculate + and − in-house, free, GPL3 software written in Python with the use of Cython (https:// github. com/ kosmo 76/ adfa) was used. Any ectopic beats were linearly interpolated according to [33].
The time series of + and − obtained for each recording were summarized by medians for the purpose of comparisons. Since the mean values of + and − did not have normal distribution (Shapiro-Wilk test), the non-parametric paired Wilcoxon test was used for their comparison.
The association analyses between ADFA asymmetric exponents and the descriptors of HRA were carried out with the use of the non-parametric Spearman correlation test. To check whether HRA entails + < − we built the univariate logistic regression models for 2 × 2 contingency tables tabulating the existence of HRA and + < − .
All statistical calculations were carried out with the use of the R statistical language and its libraries.

Presence of asymmetry in ADFA
The median values of the scaling exponents were + = 0.951 (IQR (0.825, 1.052)), − = 0.993 (IQR (0.885, 1.099)), the p value is < 0.0001 . These results should be analyzed in view of the possible values that ± may take (see discussion after formula (5)). If the median values above were different by a larger amount, this would mean a total flip in the properties of the noise after detrending.
The + < − was present in 278 cases, which is 74% of the entire group, the binomial test for this value to be consistent with 50% gives the p-value < 0.0001 Fig. 1 An example of the partitioning of an RR intervals time series into monotonic runs. The runs are marked as DRn (Deceleration Run of length n) and ARn (Acceleration run of length n); the Nn symbols stand for neutral runs which may break the deceleration/acceleration runs. Full black circles denote beginnings of decelerations runs and full gray circles mark the beginnings of accelerations runs -these can be thought of as reference points for the respective runs

Presence of asymmetry in the variance-based descriptors
The presence of asymmetry in the variance-based descriptors was established by comparing the deceleration-and acceleration-based parameters as well as checking whether globally the proportion of subjects with a specific type of asymmetry was different from the theoretically expected value of 0.5 if there is no asymmetry [1,17].
Short-term asymmetry The parameters of the distribution of SD1 d , SD1 a and C1 d may be found in Table 1. The number of cases in which short-term asymmetry can be established is 291, which is 75% of the entire group, the binomial test for this value to be consistent with 50% gives the p-value < 0.0001.

Long-term asymmetry
The parameters of the distribution of SD2 d , SD2 a and C2 d may be found in Table 2. Long-term asymmetry is present in 269, which is 69.3% of the entire group, the binomial test for this value to be consistent with 50% gives the p-value < 0.0001.
Total asymmetry The parameters of the distribution of SDNN d , SDNN a and C d may be found in Table 3.

Presence of asymmetry in the runs based descriptors
The summary of the acceleration and deceleration runs distributions may be found in Fig. 2 and in Table 4. It can be concluded that the asymmetric effect, i.e., runs of accelerations being longer than those of decelerations, can be observed for runs of lengths 4 through 12, with lengths 1-3 and over 12 being not statistically significant or too few to carry out the statistical tests. Thus, the runs-based descriptors demonstrate the presence of HRA in the studied group. Table 5 demonstrates the distributions of runs entropy for H DR and H AR shows the comparison between the two types.
The direct comparison between the entropic runs summaries yields p < 0.0001 , so there is a highly statistically significant asymmetric effect.

Associations between variance-based descriptors and local asymmetric scaling exponents
The associations between variance-based HRA descriptors and ADFA exponents were studied with the use of the Spearman correlation as well as loess type regression to visualize the type of association. The results are presented in Table 6. From the table above it can be seen that the asymmetric scaling exponents are significantly associated with both HRV and HRA magnitudes. The associations between asymmetric contributions of HRA to short-term, long-term and total variability and the asymmetric exponents are significant for decelerations and not significant for accelerations. Therefore it is necessary to study this in more detail. This is undertaken in the next section.

Associations between the presence of asymmetry in variance-based descriptors and the presence of asymmetry in local asymmetric scaling exponents
To answer the question whether or not the asymmetry observable in variance-based descriptors is related to the asymmetry observable in the scaling exponents is actually present rather than there being just an association between ± and the variance of time series we first build a 2 × 2 contingency table between the two categorical variables. Then for each type of asymmetry we build a logistic regression model to study the strength and significance of the association, i.e., the predictor in each case is the presence of asymmetry expressed according to the standard definition, e.g., using the inequalities from points 1.3, 1.4 and 1.5. The presence of asymmetry is coded as 1.
The relations between various types of HRA and ADFA are summarized in Table 7. The results of the   logistic model of the above asymmetry as predictor of + < − are presented in Table 8. It can be concluded that variance-based heart rate asymmetry and asymmetry in local ADFA are related.

Associations between runs-based HRA descriptors and ADFA
Tables 9 and 10 demonstrate the correlations between deceleration / acceleration runs and ADFA asymmetric exponents. The associations are quite strong for runs of length greater than 3. This may mean that there is both an association between heart rate asymmetry and ADFA and an association between the presence of various length runs and ADFA. Thus it is necessary to study the cooccurrence of both types of asymmetry.

Associations between the presence of asymmetry in runs-based descriptors and the presence of asymmetry in local asymmetric scaling exponents
We use the same method as above for runs of length greater than 3, since for these runs the heart rate asymmetry is clearly observable. The results of the logistic regression models for all the analyzed run lengths can be found in Table 12.  From both Tables 11 and 12 it can be concluded that the presence of asymmetry in monotonic runs predicts the presence of asymmetry in ADFA.

The association between asymmetry in the entropic HRA descriptors and ADFA
The associations between ± and the quantities summarizing all the above runs, i.e., H DR and H AR are presented in Table 13.
As we did for the other descriptors, let us build the 2 × 2 contingency table (Table 14) and the logistic model for the co-occurrence of the two types of asymmetry.
The coefficient of the model is 1.856 with p < 0.0001 , so again the association is strong.   Table 7 Associations between occurrence of HRA and the asymmetry in ADFA  Clearly, in all the types of descriptors, the presence of asymmetry in HRA entails the presence of asymmetry in ADFA.

Discussion
In the present paper we have shown that HRA in both its versions presented here, i.e., variance based and runs based, is associated with the asymmetry observable in ADFA, i.e., + < − . The study was carried out in a group of healthy young people from whom stationary 30-min-long recordings were obtained.
In the studied group we observed clear and highly statistically significant heart rate asymmetry at all time scales, i.e., short-term ( C1d > 0.5 ), long-term ( C2d < 0.5 ) and total ( Cd < 0.5).
Runs-based asymmetry is also clearly visible and significant for all observable runs of length > 3 as well as in the summary runs-based descriptors, i.e., H DR < H AR .
This group also exhibits a clear and highly statistically significant ADFA asymmetry, defined as + < − .
We have found that both scaling exponents are highly statistically significantly associated with the heart rate variability of the analyzed recordings, i.e., they were associated with total ( SDNN 2 ), long-term ( SD2 2 ) and short-term ( SD1 2 ) variability. For this reason they were also associated with the unnormalized HRA Poincaré plot descriptors like SD1 2 d∕a , SD2 2 d∕a and SDNN 2 d∕a . The associations between the scaling exponents and the magnitude of the relative contributions to asymmetry ( C1 d , C2 d , C d ) were significant for + , but not significant for − . At this point it is probably impossible to interpret this result.
The associations of ± with the magnitudes of the runsbased HRA descriptors are also highly statistically significant.
However, the strongest associations between ADFA and HRA in the RR intervals time series are observable in the agreement in between the two types. From the numbers presented in Section 3.4 on the predictive power of HRA for ADFA asymmetry, it is almost certain that both approaches describe the same physiological phenomenon, even though they use totally different methods, assumptions and even different language (HRA is based on statistical considerations, and ADFA on studying long-range correlations in dynamic systems).
An important limitation of this study should be raised at this point. Since the recordings used in the present study are quite short, the ± cannot be considered measures of  Table 11 Associations between occurrence of HRA in runs of length 4-10 and the asymmetry in ADFA   fractality and they will be strongly related with the variance of the recording. Therefore, the results involving variance-based HRA descriptors are weakened by this observation. However, the relation with the normalized descriptors is much more robust to this effect since the dependence on variance is largely eliminated by normalizing by variance-based parameters. Additionally, the results relating the incidence of asymmetry in ADFA with HRA should be fully resistant to this effect. The explanation of HRA has not yet been established. In the crudest approximation of the Autonomic Nervous System (ANS) it can be said that the parasympathetic branch of the ANS is responsible for decelerations and the sympathetic branch is responsible for accelerations. In this picture it might be tempting to ascribe the deceleration-based descriptors ( SD1 d , SD2 d , SDNN d , C1 d , C2 d , C d , DRx as well as + ) to the parasympathetic branch and the rest to the sympathetic branch. Yet, in reality, both branches can lead to both accelerations and decelerations through activation or deactivation. It would be more prudent to state that HRA reflects the interaction between both branches as well as describing the patterns of accelerations and decelerations. ± in the approach assumed in ADFA are influenced by both accelerations and decelerations and reflect the differences in noise left over after removing linear trends. Thus, their asymmetric behavior most possibly reflects the different interactions, both long and short term, present during decelerating and accelerating trends.
Possibly the best way in which the two types of asymmetry can be related is through the monotonic runs. It is fair to hypothesize that, since falling and rising trends consist of individual runs, rising trends in the RR intervals time series should be dominated by decelerating runs and falling trends by accelerating runs. Since, as we have shown in this and other papers, the dynamics of these runs differ, the long-range correlations reflected by the ADFA scaling exponents, should also differ.
The above hypothesis is strictly mathematical. As far as physiological relation is concerned,we can safely hypothesize that both measures are linked through a common, underlying physiological process. One such process which immediately comes to mind is respiratory sinus arryhythmia, which is a strong driver of heart rate variability. However, as demonstrated in [34], the link between HRA and RSA is by no means clear. The authors find no link between period variability asymmetry and respiratory sinus arrhythmia in healthy and chronic heart failure individuals. This means that further, physiologically oriented studies are necessary.
Since this is a strictly observational study we refrain from physiological explanation for HRA and its relation with ADFA as this would call for a full paper, we would, however, like to note that a few interesting and convincing physiological mechanisms have been identified in [35]. Another point that should be raised at this point is that we use a specific approach to HRA, namely the variance-based and runs-based descriptors. There are other approaches, like using predictability markers [36] which could also shed light on the analyzed problem.
To sum up, ADFA asymmetry is associated with HRA and one is an almost perfect predictor of the other.

Conclusion
ADFA, which is a new method for studying asymmetry in time series has been applied to the time series of RR intervals. It has been found that it is associated to other approaches to asymmetry in this time series. Since ADFA has unique properties, like the ability to study long-range correlations or scaling behavior of the time series, it is a promising direction in the HRA analysis, even though it does have some interpretational difficulties.
Funding This study was partially funded from the sources of the project 'Predicting adverse clinical outcomes in patients with implanted defibrillating devices' (TEAM/2009-4/4), which is supported by the Foundation for Polish Science -TEAM program co-financed by the European Union within the European Regional Development Fund.

Declarations
Ethics approval This study was accepted by the Bioethics Committee of the Poznan University of Medical Science.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. J. Piskorski is a theoretical physicist and medical biologist interested mainly in heart rate asymmetry, heart rate variability, signal processing and statistical methods in biomedical research. He is employed at the University of Zielona Gora and often consults for biomedical and pharmaceutical companies. D. Mieszkowski is a theoretical physicist, mainly interested in applications of theoretical physics to biological systems and in industrial processes. His main area of interest is Asymmetric Detrended Fluctuation Analysis. Dr. Mieszkowski works in industry for the automation company OEM Automatic Ltd. He spends his free time fishing, travelling and having fun with his family.

S.
Żurek works at the Institute of Physics of the University of Zielona Góra, Poland. His main areas of interest include the application of methods of computational physics and high-performance computing to complex problems in physics and medicine. Lately, the most explored topics concern analysis of biomedical time series, specifically the heart rate variability and complexity analysis of RR intervals. In his previous work, he contributed to the field of polymer translocation through pores and was working in the EPR spectroscopy developing some artificial intelligence methods that helped in the labeling of experimental data. He also has some experience in molecular dynamics simulations. B. Biczuk is PhD candidate at the Institute of Physics of the University of Zielona Gora, Poland. His main areas of interest include application of methods of digital signal processing to biomedical signals and study of quantum optical systems.
S. Jurga is a medical doctor with a specialization in neurology, head of the neurology department of the University of Zielona Góra hospital. His main scientific interest lies within cerebral vascular diseases, neurosonology and sleep medicine. He is interested in applying mathematical concepts to studying these areas and is eager to apply them to his clinical data.
T. Krauze works at the Department of Cardiology-Intensive Therapy and Internal Diseases, Poznan University of Medical Sciences as a research/technical worker. In his scientific research is focused mainly on autonomic modulation of the cardiovascular system, noninvasive measure of hemodynamic parameters and arterial stiffness.
A. Wykrętowicz is the head of the Department of Cardiology, Poznan University of Medical Sciences. His interests are diverse, including Heart Rate Variability and Heart Rate Asymmetry. The Department headed by prof. Wykretowicz is responsible for the most number of HRA/HRV papers in Poland and possibly also in the world. P. Guzik is a cardiologist working at the Department of Cardiology-Intensive Therapy, Poznan University of Medical Sciences, Poznan, Poland. His scientific interests are related to cardiology, including cardiac arrhythmias (e.g., atrial fibrillation), electrocardiology, and cardiovascular time series analysis. Together with Jaroslaw Piskorski, they have described the Heart Rate Asymmetry phenomenon related to unequal input of heart rate accelerations and decelerations to Heart Rate Variability. In clinical practice, he consults patients with cardiovascular diseases, e.g., acute myocardial infarction, heart failure, arrhythmias, or cerebral stroke. Several interdisciplinary papers related to bioengineering, medical physics, and signal analysis are among his publications.