Possible Caveats of Ultra-Short Heart Rate Variability Reliability: Insights from Recurrence Quantication Analysis

— The heart rate variability (HRV) is the difference between consecutive R-R intervals of heartbeats measured in milliseconds. HRV indices represent the role of sympathetic and parasympathetic autonomic branches. Even though HRV is considered an indirect biomarker of Autonomic Nervous System, there are not yet standardized protocols providing reliable clinical measures. One of the reasons is because HRV techniques requires long recording periods. There are attempts of decreasing the required recording, such as the strategy of ultra-short HVR recording (< one minute), which could make the uti-lization of the technique easier. However, there is little published about its reliability. This work proposes a method to evaluate the reliability of ultra-short HVR based in Poincare map and Recurrence Quantification Analysis, well known methods to assess nonlinear and dynamic information from a system, in order to verify the reliability of the use of ultra-short term HRV. Then, these results was compared with the classical HRV coefficients, such as rMSSD, recorded from subjects in spontaneous breathing and also, in controlled breathing protocols. As a conclusion, using the proposed methods, we were able to show the discrepancy between the segments of interest, both on mean and in variance, explained in the analysis of main components.


I. INTRODUCTION
Heart rate variability (HRV) represents a variation between consecutive R-R intervals of heartbeats measured in milliseconds [1]. The fluctuations of HRV are controlled by homeostatic systems, including baroreceptors, respiration, thermoregulation. Also, the cardiovascular system is influenced by sympathetic and parasympathetic autonomic branches [2]. The quantitative analysis of HRV can be performed in time and frequency, providing indices related to the sympathetic and parasympathetic autonomic axis [3]. Several studies report that clinical disorders like obesity, diabetes, cardiovascular diseases, psychiatric disorders present impairment on HRV indices, indicating an autonomic dysfunction [4][5][6][7][8]. Thus, HRV is considered a non-invasive biomarker of the autonomic nervous system (ANS) suitable for clinical and health applications [5][6][9][10][11]. Although clinical findings have already been reported by groups of researchers [2][3][4][5], HRV assessment is not part of the clinical routine due to the lack of standardized protocols compatible with the clinical setting.
HRV analysis consists of the measure of successive R-R or N-N peaks from electrocardiogram (ECG). This temporal series, also known as tachogram, can be analyzed commonly through frequency or time domain strategies. In the frequency domain, it can be extracted the High Frequency (HF) and Low Frequency (LF) indices from the tachogram, in which HF correspond only to indirect vagal modulation, and LF indices correspond indirectly to a mixed response of both vias (sympathetic and vagal) [12].
Recently, there has been renewed interest in the improvement of methods for HRV recording to refine clinical suitability [13][14][15]. The gold-standard methods (i.e., electrocardiogram, ECG) for R-R intervals recording presents high correlation coefficients to R-R intervals derived from wireless heart rate modules [13], reducing the financial and operational cost of the HRV technique. To obtain reliable HRV data, it is recommended recordings lengths from five minutes to 24 hours [1][2][3]. Recent studies show some HRV indices that present good correlations between gold-standard time length (5-minutes) and one-minute time epochs proposing an ultra-short term HRV recording (<1-minute) [14 -16]. The results report that Root-Mean-Square of the Successive Normal Sinus R-R interval difference (rMSSD) is the adequate parameter for accessing HRV from ultra-short recordings, while other time or frequency domain indices require longer recording periods [14][15][16].
There is little published data on the reliability of ultrashort HRV data applied to clinical sample or experimental conditions, as most studies involve healthy samples or are recorded at rest. The currents evidence suggest that ultra-short-term rMSSD: a) is sensitive to physical-training induced changes in cardiac autonomic tonus of elite athletes [17][18]; b) is associated to mental stress and emotional frustration perceptions after cognitive workload [19][20]; c) can be used efficiently for autonomic evaluation in diabetes mellitus patients [21]. Nevertheless, ultra-short-term rMSSD reliability remains unclear and many investigations have adopted ultra-short analysis without questioning its validity, there is a clear lack of rigorous methods to identify indices for ultra-short recording [22]. So, other methods of HRV indices extraction could be used to better comprehension of R-R intervals dynamic across distinct recording lengths.
Non-linear strategies such as Poincaré Map (PM) and Recurrence Quantification Analysis (RQA) have been applied to assess dynamics related to ECG and HVR. PM has been used to detect abnormal dynamics of cardiac repolarization [23]; to evaluate dynamic autonomic modulation during general anesthesia induction [24]. There are attempts of using RQA for analysis of heart rate variability and respiratory flow series in patients on weaning trials [25]; heart rate dynamics in young patients with diabetes mellitus [26].
RQA is used for quantification of the number of recurrences of a dynamical system and its duration presented by its phase space trajectory [27]. PM is characterized as a discrete dynamical system that represents the continuous periodic flow of another system [28]. Thus, due to the capability of RQA and Poincaré Map to assess nonlinear and dynamic information from a system, we will apply both strategies to verify the reliability of ultra-short term HRV. Then, the results will be compared with the classical HRV coefficients, such as rMSSD, recorded from subjects in spontaneous breathing and also, in controlled breathing protocols.

II. MATERIALS AND METHODS
To evaluate the reliability of ultra-short-term heart rate variability, we performed the HRV analysis by evaluating the tachogram characteristics from two groups, n = 20 subjects under controlled breathing and the same subjects under spontaneous breathing, by using Poincaré Map coefficients SD1 and SD2, through six coefficients from Recurrence Quantification Analysis, that are: Recurrence Rate (RR), Determinism (DET), Entropy (ENT), Maximal diagonal line length (Lmax), Laminarity (L), Trapping Time (TT). Also, it will be obtained the classical rMSSD coefficient from HVR signal.

A. Participants
The sample of this study was composed of 20 participants (all male) aged between 18 and 28 years old (M = 23.65 ± 3.24). The following exclusion criteria were applied to the sample: (a) Do not follow the pace of the breathing pacer (6 breaths / min) satisfactorily or feel discomfort when performing this breathing rhythm; (b) Diagnostic of cardiac dysfunction; (c) Do not drink alcohol (24h before) or caffeine (3h before) drinks before the experiment. This study was approved by a local Ethics Committee, and was conducted according to the Helsinki Declaration.

B. Physiologic Recording and HVR analysis
PolarH7 heart rate monitor (Polar, Finland) was used for the acquisition of RR intervals (iRR, in ms). All acquisitions have occurred in the period between 9 and 12 hours in a well illuminated and quiet room. The resting baseline was made in sitting position, followed by a breathing paced condition (6 breaths/min), both periods with 6-minute length. To reduce the effects of autonomic adaptation to body posture, the participant received the guidelines of the experimental protocol, signed the consent form and performed a test record in the sitting position To evaluate different periods of time, four fragments were randomly selected, in which the first is 6-minute length, the second of 3-minute, third with 2-minute, and last with oneminute. Each fragment was sorted 100 times to minimize the bias of selection of the evaluation periods.

C. Poincaré Map
The Poincaré Maps was generated by scatter plots given by the past R-R intervals (RRt-n) against the present R-R intervals (RRt). This is done qualitatively by studying cluster shapes and quantitatively using cluster deviations SD2 and SD1, which represent the major axis and the minor axis of an ellipse, respectively. SD2 quantifies the point distribution across the line of identity (LOI), and SD1 indicates the point of distribution across the perpendicular line to LOI.

D. Recurrence Quantification Analysis
RQA is based on the evaluation of dot-structures presented on the Recurrence Plot (RP) , through different nonlinear factors [26]. Each dot ( , ) of the RP is calculated according to the recurrence rule, described by , = 1 , if ‖⃗ − ⃗ ‖ ≤ , and , = 0 otherwise; where ‖ . ‖ is the Euclidian distance metric, is the recurrence threshold and ⃗ are vectors pointing to the amplitude coordinates in the phase space ( , + , … , +( −1) ). As the columns and rows of a RP correspond to temporal coordinates of the time series, the elements of the matrix indicate each moment that a relevant state (according to a recurrence threshold ) of a dynamic system repeat along time [28].
To quantify the HVR recurrences for each period, we choose six factors displayed by the Recurrent Plots (RPs) that are [26, 28]: • Recurrence Rate (RR) -quantifies the density of recurrences.

E. Statistical Analysis
All data were submitted to the Kolmogorov-Smirnov normality test for data verification [29]. Since the samples were not considered as following a normal distribution, nonparametric tests were used.
Group differences were assessed using Kruskal-Wallis test, followed by the Sîdák post-hoc test [30]. For two-sample comparisons, we used the Mann-Whitney test. Additionally, it was used Bartlet test to compare differences in variances [31]. Each test was repeated 100 times to evaluate the consistency and time-series stationarity in relation to each epoch 1 min, 2 min, 4 min and 6 min. The significance level for all statistical analyses was established as α = 0.05. Principal Component Analyis (PCA) [32] was applied to explore the differences between RQA coefficients and PM features obtained from rMSSD for the fragments of all lengths.

III. RESULTS AND DISCUSSIONS
In Fig. 1, a fragment of HRV signals based on R-R successive measures from spontaneous breathing and controlled breathing is presented. After that the recurrences from controlled and spontaneous breathing were plotted and presented in Fig. 2. It is possible to notice the amount of differences in the phase space and in the RP of each signal that the attractors can capture. While the attractor and the RP of the signals with controlled breathing present the most regular and deterministic patterns, the attractor and the RP of spontaneous breathing have classic patterns of coupled systems, with tendencies of chaoticity and stochasticity.  Fig. 3 shows examples of Poincare Maps for each period of analysis considered (6 min, 3 min, 2 min and 1 min). It is possible to see how the reduction of the time length of the R-R signals, contributes for a clear the visualization of the Poincare map. This could also impact on the quality of any analysis that is going to be carried out. The distribution of the signals with controlled breathing, despite being more regular, have greater amplitudes, which is also possible to notice in the PM. The following step was to verify if the coefficients obtained by the PMs and the RMSSD from spontaneous and controlled breathing R-R signals were significatively different (Fig. 4). For this analysis, it was compared the data obtained for each period of interest (1 min, 2 min, 4 min and 6 min). Hundred fragments were randomly selected with reposition throughout each record of each participating subject.
Above each plot it is possible to see the proportion in which the 100 comparisons analyzed (from the 100 random selections), showed statistical difference, in each test. For the first analysis of SD1, for example, from 100 comparisons, none was significant between the different time lengths and only 19 comparisons showed any significance between the 4time lengths of interest). Then, it was performed a comparative analysis between mean values and variance from the RQA coefficients for controlled and spontaneous breathing (Fig.5). The first 6 plots exemplify the 100 comparisons made using Kruskal-Wallis (with Sîdák, posthoc) for comparisons between means and Bartlet for comparisons between variances. The title of each plot also shows the number of times, proportionally, each comparison was significant, in the 100 comparisons performed. For both controlled and spontaneous breathing, the Lmax and ENT coefficients were the most relevant in the all comparisons. Finally, it was plotted a Principal Component Analysis using coefficients PM+rMSSD and the coefficients of RQA for R-R signals, for controlled and spontaneous breathing. The first row of the plots shows that using only the coefficients SD1, SD2 and RMSSD as features in the PCA, it is impossible to make a distinction between the four periods of interest. This first result may precisely lead to the misleading conclusion that there are no differences between the four sections. However, when using the RQA coefficients as features of the PCA, the difference in each of the four analyzed sections is evidenced. For the signals of controlled breathing it is possible to observe the large difference in the average of the clusters and their spread, especially between the clusters associated with the 1 min and 6 min periods. With the RQA, it is also possible to visualize a separation by blocks between each period, distinguishing them completely, even in a situation of spontaneous breathing.

IV. CONCLUSIONS
Most standard techniques, that estimates temporal characteristics (or in the frequency domain) of tachogram (R-R signals) to assess heart rate variability are quite promising and, in many situations, efficient. Although, it is necessary to be careful with the conclusions taken from these indices since analyzes bias may occur because of the lack of a standard criterion for a minimum registration time for the estimation of these factors. As the cardiorespiratory system is a non-linear system, with characteristics of stochasticity and chaoticity, it is possible that standard techniques are not able to detect important non-linear fluctuations, which are added up in larger periods, and are extremely important for in HRV analysis. This is evident when applying techniques capable of detecting these non-linear fluctuations such as RQA. Using the RQA, we were able to show the discrepancy between the four segments of interest, both on mean and in variance, explained in the analysis of main components.
As a limitation of this study, it is important to mention that cardiac autonomic tone adaptability to breathing frequency can vary for each individual. This is a relevant bias to be considered when using a controlled breathing rhythm. In slow paced respiratory frequency protocols, the time difference of automatic cardiac tone mobilization between the participants can influence on the replication of short time epochs. Another important limitation of controlled breathing is that the cardiac autonomic tone adaptation to the body posture can generate non-stationary artifacts that impair the replication of ultra-short HRV epochs. Thus, to obtain replicable epochs to ultra-short HRV analysis, it is necessary a period of adaptation to cardiac autonomic tone to the posture.
ACKNOWLEDGMENT This work was supported by PRONEX Program (Programa de Núcleos de Excelência -NENASC Project) of FAPESC-CNPq-MS, Santa Catarina Brazil (process number 56802/2010). RW is a Researcher Fellow from CNPq (Brazilian Council for Scientific and Technologic Development, Brazil) and HMM is supported by CAPES/DS scholarship. We also would like to thank for the Call of the Support of Assistive Technology Projects -TA2016/CNPq, process number: 442563/2016-7