On the Standardization of Approximate Entropy: Multidimensional Approximate Entropy Index Evaluated on Short-Term HRV Time Series

,


Introduction
Approximate entropy (ApEn) was introduced by Pincus and coworkers [1] in 1991 as an entropic measurement to quantify the regularity of medical data. Ten years later, Richmann and Moorman [2] introduced sample entropy (SampEn), a variation of ApEn reducing the bias of considering selfcomparisons and being more independent of data length. Both approaches have been widely used to characterize medical disorders or discriminate between healthy and pathological conditions [3][4][5][6][7][8][9]. However, these entropy measures are by definition dependent on two predefined parameter values, namely: m, the embedding dimension (i.e., the length of reconstructed vectors), and r, the tolerance threshold (i.e., similarity value for comparing reconstructed vectors). These parameters have been assigned diverse values in published heart rate variability (HRV) studies (e.g., m = 1 or 2; r = 0.1 to 0.25 times the standard deviation of 2 Complexity the time series) [9][10][11][12][13][14][15][16][17][18][19]. Thus, one main drawback of these entropy indices lies in the comparison between studies, since diverse values of a priori parameters can lead to different physiological interpretations.
Entropy-based methods, and in particular SampEn-based ones, have received great attention in the last ten to fifteen years. Alternatives have been proposed to make this index independent of a priori parameter definition, thus facilitating the standardization of its application. In this regard, multiscale entropy (MSE) and refined MSE approaches were introduced to take into account complexity properties at different scales by applying the coarse-graining technique [15,[20][21][22]. Each scale was characterized by a new time series derived from the original one by averaging and decimation [15]. However, they still required the predefined parameters m and r. Alternatives, such as the use of a density function instead of a tolerance value, was proposed to improve estimation of Renyi entropy, generalizing the Shannon entropy [23,24]. Other entropy-based approaches have been proposed to avoid the need for predefining specific parameters. In this regard, permutation entropy of a given order was calculated from the relative frequencies of all permutations of that order of the analyzed time series [25]. Dispersion entropy was suggested as a more reliable method than sample or permutation entropies. However, the need for defining the tolerance threshold r is not strictly avoided, but this is replaced with the definition of a number of classes [26]. Rankbased methods were investigated to quantify the amount of shuffling by sorting the ranks of the mutual distances between pairs of m-long vectors [27]. Although this approach avoids predefinition of the threshold value, the embedding dimension requires a priori definition. On the other hand, the cumulative histogram was proposed to overcome the limitation of defining a tolerance r value in the computation of sample entropy, yielding the approximate entropy profile [28]. An approach to measure bubble entropy, which was shown to be independent of the embedding dimension m for sufficiently large values of this parameter, was introduced in [29]. Several of the above described approaches showed good performance to characterize pathological versus healthy conditions and led to a reduction in the associated computational load. Nevertheless, they still presented dependence on some parameters, like the embedding dimension or the tolerance value.
Although SampEn was introduced as an improvement of ApEn entropy, advances towards the standardization of ApEn have been investigated as well. Different approaches have been proposed to identify the values of the embedding dimension m or the threshold r used for the characterization of nonlinear dynamics in chaotic time series [30][31][32][33][34][35]. With regard to the embedding dimension, the false nearest neighbor method was used to search for the lowest embedding dimension m that allows phase-space reconstruction. Very different values of m were reported for different models and experimental data [30]. The lowest values of m were found to be 2, 3, and 6 for Hénon, Lorenz, and Mackey-Glass time series, respectively [31]. With regard to the tolerance threshold r, some studies focused on searching for the value, denoted by r max (m), that maximizes ApEn [32,33,36]. The corresponding ApEn-based index, maximum ApEn, was used to show that white noise series was more irregular than crosschirp signal [37]. In addition, this index was reported to better enhance HRV and blood pressure variability differences between supine and upright position with respect to the use of fixed r values, typically 0.1-0.25 times the standard deviation [34]. ApEn profile, based on cumulative histogram computations, was proposed to avoid the shortcomings of a priori definition of the tolerance threshold [38]. HRV indices derived from such ApEn profiles were shown to better separate young versus old individuals as compared to classical ApEn approach. To reduce the influence of the threshold r on the previous method, an adaptive cumulative histogram method was proposed [39]. Total and average cross-approximate entropies, as well as the standard deviation of cross-approximate entropy, distinguished financial time series successfully. It should be noted that although r is generally referred to as a fixed value, the parameter that is fixed is the factor multiplying the standard deviation of the analyzed time series, rather than the threshold r itself. In the literature, identification of the values of r and m used to characterize different time series is still unclear, as highlighted by the highly diverse published values [34,35].
In this study, we aimed to introduce an ApEn-based index that avoids the need for a priori parameter definition and, thus, allows standardization of entropy-based measurements. In particular, is proposed based on summation of maximum ApEn values over different embedding dimensions for each given HRV time series. This multidimensional index does not require the use of the coarse-graining technique and overcomes the variations in data length for different time scales. Furthermore, its computation provides entropic estimation without the need for ad hoc selection of input parameter values. The ability of the proposed index to represent different degrees of randomness in time series is first tested in synthetically generated signals and subsequently used to characterize HRV changes induced by aging and by congestive heart failure (CHF).

Materials and Methods
2.1. Materials. MIX processes: a family of processes combining deterministic and stochastic behavior were studied. The degree of stochasticity was controlled by parameter P. MIX(P) generated a sine for P = 0 (pure deterministic) and became more random as P increased up to 1 (P = 1, pure stochastic):  [41]. Recordings can be downloaded from http://www.physionet.org [42]. The whole 2-hour recordings were analyzed. These databases were available from http://www.physionet.org [42]. The analyzed heartbeat locations of all databases were obtained from an automatic ECG annotator and were subsequently manually reviewed and corrected by experts among Physionet collaborators. In this study, those heartbeat locations were corrected for ectopic beats based on instantaneous heart rate variation [43]. Subsequently, the time intervals between consecutive heartbeats were used to define the RR interval time series. Each RR time series was divided into 300-sample segments, with 50% overlapping, and interpolated at 2 Hz to attenuate the effect of heart rate mean as sampling rate [44]. A 3-hour night period centered on the 300-sample segment with minimum heart rate mean (HRM) was selected in all cases to minimize the influence caused by potentially different daily activities.

Congestive Heart
Estimation of linear and nonlinear HRV indices was performed on each segment. Subsequently, the median value of each index over those segments was computed as the representative value for each of the analyzed time series.
The amount of reconstructed vectors was = − ( − 1) for each embedding dimension m. The distance, computed by the ∞ norm, between each pair of reconstructed vectors, , , was denoted by , . Each distance was compared with a threshold r to compute the number of reconstructed vectors that lie within a hyperspace centered in the reconstructed vector of reference: where is the correlation sum and H is the Heaviside function.
This procedure was repeated with all reconstructed vectors and the probability of a pattern of length m appearing along the time series was denoted by the following.

Maximum Approximate Entropy.
For each embedding dimension m, the threshold r leading to maximum approximate entropy, denoted by r max (m), was determined as in [45] and the maximum value ApEn(m, r max (m)) was calculated [36]. For the sake of comparison with ApEn(2, 0.2), ApEn(2, r max (2)) was calculated.

Multidimensional Approximate Entropy.
The multidimensional entropy index, MApEn(r), was introduced by summing up the approximate entropies for a range of embedding dimensions: where m max = 15. To additionally remove the dependence of the tolerance threshold, the value r max (m) was identified for each embedding dimension m, leading to the definition of .
ApEn, SampEn, and correlation dimension D 2 were computed using the methodology described in [45], considering a range for the tolerance threshold r varying from 0.01 to 3 times the standard deviation of the time series with a resolution of 0.01.

Time-and Frequency-Domain HRV Indices.
Time-and frequency-domain indices were obtained from the raw RR time series and from the modulating signal of the integral pulse frequency modulation (IPFM) model, respectively, using 300 heartbeats as analysis windows. HRM and the square root of the mean squared differences of successive normal heartbeat intervals (RMSSD) were calculated.
Spectral analysis was performed on a modulating signal, m(t), assumed to carry information about the autonomic nervous system (ANS) activity, which was estimated based on the IPFM model compensating for the influence of HRM [46]. Power in the low frequency (P LF , 0.04-0.15 Hz), related to sympathetic and parasympathetic activity, and high frequency bands (P HF , 0.15-0.4 Hz), as a parasympathetic modulation estimation, as well as normalized low frequency power (P LFn = P LF /(P LF +P HF )) was estimated [47].

Statistical Analysis.
The normality of data distributions was tested by the Kolmogorov-Smirnov test. Nonlinear as well as time-and frequency-domain HRV indices were compared between young and elderly subjects and between CHF patients and healthy subjects by t-test or Mann-Whitney U test depending on whether data distributions fulfilled the normality criterion or not. Only paired comparisons (rather than other multiple comparisons) were performed. Spearman correlation between the proposed index and other HRV indices was computed. versus 0.75) all nonlinear indices led to statistically significant differences. In the comparison between MIX(P) processes with P above 0.75 (0.75 versus 1) no statistically significant differences were found for any of the nonlinear indices. ApEn (2, 0.2) and ApEn (2,r max (2)) as well as nonlinear index values were higher for greater levels of randomness until MIX(0.75); however above P = 0.75 an increase in the level of randomness was not captured by ApEn or by any of the nonlinear indices. Regarding the separation of pink and white noise time series, all the analyzed indices were able to separate them. The range of r max (2) values was found to lie in the interval 0.2-0.3 times the standard deviation of the time series, which is commonly used in the literature (see Table  S1 in Supplementary Materials for additional information). The multidimensional index MApEn(r) was able to separate pink and white noise time series for all tested values of the tolerance threshold r (results only shown for r = 0.2, Figure 2(a)).

Synthetic
While statistically significant differences were found for both MApEn(0.2) and when comparing different stochastic processes, the former presented decreasing values as randomness increased, while the latter showed opposite results. To analyze these differences, the results of evaluating ApEn for two embedding dimensions, namely, m = 2 and m = 6, are illustrated in Figure 3. The use of a fixed threshold did not always allow ApEn to characterize the randomness level of the MIX(P) processes. As an example, for m = 6, the use of the fixed threshold r = 0.2 led to MIX(0.25) having much larger complexity than MIX(0.5), MIX(0.75), and MIX(1). For m = 2, however, MIX(0.25) presented the lowest complexity. When combining results from different scales, MApEn(0. 2) showed lower values as the randomness level increased, thus not being able to reflect the degree of stochasticity in the time series (Figure 2(a)). However, the use of ApEn(m, r max (m)) when computing served to characterize MIX(P) processes according to their stochasticity level for any value of m, which then rendered suitable to characterize stochasticity levels (Figure 2(b)).
In Figures 4(a) and 4(b) ApEn(m, 0.2) and ApEn(m, r max (m)) are shown for pink and white noise when the embedding dimension m was varied from 1 to 15. ApEn(m, 0.2) values were higher for white noise than for pink noise for low embedding dimension values, while presenting the opposite behavior for large embedding dimensions. The multidimensional index MApEn(0.2) resulted in higher values for white than for pink noise, reflecting the behavior found for the lowest embedding dimensions (Figure 2(a)). On the other hand, ApEn(m, r max (m)) was consistently larger for white than pink noise across all embedding dimensions and, consequently, differences between these noise time series were amplified when computing ( Figure 2(b)) with respect to MApEn(0.2) (Figure 2(a)).
As the statistical significance in the comparison of simulated data highly depends on the number of simulated processes (30 in the comparisons shown in Figure 2(b)), a study was performed to characterize such a dependence for the proposed index when this is compared in all cases shown in Figure 2(b). Results are presented in Table 1 for a number of processes varying from 5 to 30 in 5-process steps. As shown in the table, p-values decreased as the number of simulated processes increased. For a pvalue of 0.05, statistical significance was already achieved for a number of simulated processes as low as 5 in all comparisons indicated as significant in the results presented in Figure 2(b) corresponding to 30 simulated processes.
The median spectrum obtained from each type of bandpass filtered noise time series is illustrated in Figure 5 (2)), ApEn (6,r max (6)), r max (2), and r max (6) values are displayed for each case. * * * * * * * * *   series were found only for by Mann-Whitney U test (Figure 5(b)). Additionally, the dependence of the statistical significance on the number of realizations was studied for RMSSD, P HF , P LF , P LFn , and (see Table  S6 in Supplementary Materials).

Real Data.
One subject from the young cohort and another one from the elderly cohort were discarded for analysis, as their RR time series were not suitable for ectopic beats correction due to the fact that their variability values reached the maximum allowed HR variation, thus precluding distinction between normal and ectopic beats [43].
ApEn(2, 0.2), ApEn(2, r max (2)), MApEn(0.2), , SampEn(2, 0.2), and D 2 values showed a trend for higher values in young than in elderly subjects, but only showed statistically significantly higher values in young subjects (see Figures 6 and 7). On the other hand, all analyzed nonlinear HRV indices, with the exception of MApEn(0.2), were statistically significantly larger in CHF patients than in healthy subjects (Figures 6 and 7). HRM and RMSSD values are displayed in Figure 9 for the studied groups. No statistically significant differences were found for HRM when comparing young and elderly subjects, although a tendency to lower HRM values in the elderly group was observed. HRM was statistically significantly higher in CHF patients as compared to healthy subjects. On the other hand, the RMSSD index showed statistically significant differences between both pairs of groups, with elderly subjects and CHF patients presenting reduced RMSSD values.
Values of the frequency-domain HRV indices P LF , P HF , and P LFn are shown in Figure 10 for the four studied groups. Significant differences were found for P LF and P HF in the two comparisons. P LFn only showed statistically significant differences for the comparison between CHF patients and healthy subjects. Figure 11 illustrates the relationship between and RMSSD (a) and P HF (b). Correlation coefficients computed over the four analyzed groups are shown in Table 2.

Searching for Standardizing Approximate Entropy.
Nonlinear HRV analysis has extended the description of ANS regulation of cardiac electrical activity. The use of nonlinear HRV indices has provided new physiological insights into ventricular dysfunction, ventricular tachycardia, obstetrical complications under anesthesia, mental disorders, and aging, among others [3][4][5][6][7][8][9]. Nevertheless, the physiological interpretation of some common nonlinear HRV indices could be biased by the a priori selection of parameter values intrinsic   (2)). In this study, the index has been proposed based on summing the maximal approximate entropy values over all analyzed embedding dimensions. Although the computation of requires setting a maximum embedding dimension value, the relative contribution of ApEn(m,r max (m)) values to the final estimation largely decreases as m increases, thus rendering practically independent of such selection. In fact, in this study ApEn(m,r max (m)) was found to present a remarkable decay for m varying from 1 to less than 10, while it remained approximately constant or with a very slow decay from that point onwards. Thus, a value of 15 or any other larger value for the maximum embedding dimension is expected to provide essentially the same outcomes. As a consequence, the proposed index, , can be seen as an entropybased index not requiring a priori parameter definition for its computation. The algorithm described in [45] for D 2 estimation was used in this study to compute the analyzed nonlinear indices. By using this algorithm, nonlinear HRV indices could be computed in less than a second for a 300sample time series, whereas the computational time for the sequential algorithm over the same time series was 1086 seconds (18 mins), in both cases using a Windows5 7 based PC, Intel5 Core6 i7 3.5GHz, 16Gb RAM with Matlab5 R2015a.
Other SampEn-based indices measuring time series irregularity such as rank entropy, permutation entropy, and bubble entropy have been proposed as methods aimed at providing entropy definitions independent of a priori selected parameter values [25,27,29]. The multiscale extension of these SampEn-based methods, by considering the coarse-graining derived series, served as a basis for this study to formulate a multidimensional approach, overcoming the coarse-graining limitations for short-term HRV data analysis.

Analysis of Synthetic Data.
The proposed index has been tested for characterization of synthetic time series. Statistically significantly larger values of were found for white than for pink noise. This behavior was also reproduced by all other tested nonlinear HRV indices including MApEn(0.2), although such differences were found to be attenuated when using a fixed tolerance threshold as compared to using the one leading to maximum entropy.
The evaluation of the proposed marker over synthetic data generated by MIX(P) processes revealed that the greater the stochastic level, the higher the value of , which is in agreement with the results shown by other nonlinear HRV indices, like D 2 , ApEn, and SampEn, describing complexity of the time series. On the contrary, the use of a fixed threshold r, e.g., with value 0.2, in the computation of MApEn(r) was not able to reflect the level of randomness of the MIX(P) time series.
An additional simulation study was carried out by analyzing three types of filtered white noise time series whose timeand frequency-domain indices were statistically similar. The analysis showed that provides complementary information to that from time-and frequency-domain indices. According to the simulation results, was sensitive to changes in spectrum shape involving changes in HRV regularity. The narrower the spectral components, the lower the value is, according to a more regular signal.
In the comparisons of simulated populations, the number of stochastic processes is expected to influence the statistical significance. Results have been presented to confirm the decrease in the p-values when the proposed index was compared over different synthetic populations with an increasingly higher number of simulated processes. The number of processes required to attain a given significance was derived for each performed statistical comparison. Our results, in line with previous studies published in the literature [48], point out the importance of reporting the number of considered processes along with the p-value of any statistical comparison over synthetically generated data.

Aging and Congestive Heart
Failure Assessment by Nonlinear HRV Analysis. The effect of aging was evaluated in healthy subjects during supine resting conditions, with subjects being awake. Lower values were found in the elderly cohort with respect to the young one for all analyzed nonlinear indices, with differences being statistically significant for our proposed index. Our results are in agreement with previous studies where aging has been reported as a cause for decreased complexity and irregularity of the beatto-beat RR series [8,20,38,41,49]. This reduction has been associated with an impairment to adapt against external or internal perturbations [41,50]. HRM, RMSSD, and P HF values were also decreased, while P LFn values were increased in the elderly cohort, pointing to a potential enhancement in the sympathetic modulation of sinoatrial node activity.
The impact of CHF was assessed by comparing all our analyzed indices between failing patients and healthy controls during the night period while they were sleeping. According to the results obtained for all the studied nonlinear HRV indices, including the proposed index, statistically significantly higher complexity and irregularity were found in CHF patients with respect to healthy subjects. Similar results were found when restricting the age of the individuals to the range covering from 55 to 75 years old, pointing out that age has no relevant effect on the comparison between CHF patients and healthy subjects for the compared cohorts. Similarly, NYHA class had no impact on the distinction between CHF patients and healthy subjects, as demonstrated by restricting the analysis to CHF patients in NYHA classes III-IV, which led to similar results to those found for the whole CHF population (see Figure  S2 in Supplementary Materials for additional details on the performed comparisons).
Spectral HRV analysis showed a significant decrease in P LFn for CHF patients with respect to healthy subjects. In previous studies, CHF has been reported to be associated with an increase in the sympathetic tone and a decreased peripheral response to adrenergic input [51,52]. Studies considering pharmacological blockades have demonstrated that sympathetic blockade leads to increases in ApEn and decreases in P LF and P LFn , while sympathetic activation decreases the values of all nonlinear indices and increases P LFn [53,54]. Thus, an increase in P LFn as well as a decrease in irregularity/complexity indices could be expected in CHF patients when compared to healthy subjects. However, in our data both a decrease in P LFn and an increase in irregularity/complexity indices were observed in CHF with respect to controls. A decrease in the low frequency content of HRV has been related to the progression of heart failure in CHF patients with advanced disease [52]. This apparently contradictory results may be explained by a decreased responsiveness to sympathetic modulation [55]. CHF has been characterized by ApEn in the literature, but results have been found to be dependent on the applied method [56]. For instance, the use of a pointprocess paradigm for describing the heartbeat generation led to greater irregularity and/or complexity values in healthy subjects as compared to CHF patients, whereas the use of raw RR interval time series showed the opposite tendency [56]. MSE, which extends the classical SampEn definition to a time-scale representation, has been applied to describe heart rate dynamics in CHF [20]. Costa et al. reported that healthy subjects presented greater MSE values than CHF patients, but in particular the elderly healthy cohort showed lower MSE values in the first scale than the CHF group. Although a relationship between MSE-derived indices and sympathetic and parasympathetic modulation has been reported [57], the methodological differences comparing MSE-derived indices and do not allow translation of the physiological interpretation to the present study.
In our study CHF patients and healthy subjects showed diverse ApEn values, being greater in one or the other group depending on the tolerance r and the embedding dimension m values. Therefore, the different results obtained in the present study and in [21] could be due to the method itself (ApEn here and SampEn in [21]) and/or to the values used for the tolerance threshold (r max (m) in this study and fixed thresholds in [21]). ApEn(m, r max (m)) values were found to be larger in CHF patients than in healthy subjects for all considered embedding dimensions, contrary to the lower MSE values in CHF patients reported for all scales in [21]. Additionally, to discard the effect of the segment length of analysis in underlying differences between the results presented here and results in [21], our proposed multidimensional index was also tested for longer segments of HRV analysis (up to 4000 samples rather than 300 samples). Results confirmed that presented essentially the same discriminative capacity between healthy and CHF groups for the different analyzed segment lengths, with higher values corresponding to CHF patients (see Figure  S3 in Supplementary Materials). This led us to conclude that the selected segment length for HRV analysis was not a factor explaining the higher values found in CHF patients.
Of note, the CHF patients included in the present study were enrolled in a long-term study to evaluate the efficacy of drugs, in particular milrinone and digoxin, and these drugs may have had an effect on autonomic modulation of heart rate. In addition, no information was available in relation to other concomitant disorders that CHF patients might suffer, such as obstructive sleep apneas, particularly considering that these were evaluated during the night period [58].
Finally, our proposed index, even if not showing remarkably higher discriminative capacity than other traditional HRV indices in separating young versus elderly subjects and CHF patients versus healthy subjects, could be capturing different information from that of those other indices and, consequently, it might prove useful as a complementary tool to characterize the effects of diseases or conditions. This is supported by the simulation study considering filtered white noise time series whose time-and frequency-domain indices were invariant. Other studies have also reported that the combination of entropy-based indices with time-and frequency-domain ones, when analyzed in patient populations similar to those of the present study (CHF patients versus healthy subjects), provided better classification results than by using only time-and frequency-domain indices [59,60]. The added value of in clinical and experimental applications should be the focus of future research studies.

Conclusions
A multidimensional approximate entropy index, , is proposed as a parameter-free entropy-based measurement. By considering maximum approximate entropy values across embedding dimensions, the proposed index increases the estimation robustness. This has been tested by evaluating the methodology on synthetic time series, including MIX(P) processes as well as pink and white noise, confirming that consistently characterized the degree of randomness in the series, whereas ApEn computed by using fixed thresholds produces inconsistent results depending on the analyzed embedding dimension. provided additional information in the analysis of the effects of aging. Specifically, was found to be higher in elderly than young subjects, thus capturing heart rate complexity changes possibly reflecting loss of ANS adaptability in the elderly population. In addition, was found to be increased in CHF patients as compared to healthy subjects during the night period, indicating greater HR complexity due to CHF. In conclusion, can provide robust information on nonlinear HR dynamics with the advantage of being independent of preselected parameter values.

Data Availability
Synthetic data can be found in the following link: https://drive .google.com/open?id=1wCjF0dJmmXqdwn7Qx371GHXrN N0WyjP. The databases, Fantasia Database, Congestive Heart Failure RR Interval Database, BIDMC Congestive Heart Failure Database, and Normal Sinus Rhythm, used for nonlinear HRV analyses can be downloaded from http://www.physionet.org.

Disclosure
CIBER is a center of the Instituto de Salud Carlos III in assistance from the European Regional Development Fund. The computation was performed by the ICTS "NANBIO-SIS", more specifically by the High Performance Computing Unit of the CIBER in Bioengineering, Biomaterials & Nanomedicine (CIBER-BBN) at the University of Zaragoza.

Supplementary Materials
Supplementary Materials are provided. Three figures and five tables that can be found in separate files. Table S1.
(2) values computed on synthetic and real data. Mann-Whitney U test was used to evaluate statistical differences between consecutive row data. Data are shown in terms of median (1st|3rd) quartiles. Figure S1. Approximate entropy values computed for thresholds r set at 0.1 and 0.15 across embedding dimensions in young (black) and elderly (grey) groups shown in the upper panels, and in healthy subjects (dotted grey) and CHF patients (dashed black) during the night period in the lower panels. * indicates statistically significant differences. Data are shown as median and interquartile range. Figure S2. Results of nonlinear indices computed on two subsets from healthy and CHF databases during night period, restricting age range covering from 55 to 75 years old and restricting CHF patients in NYHA classes III -IV. * indicates statistically significant differences. Data are shown as median and interquartile range. Figure S3.
MApEn max values computed for young (solid black) versus elderly (solid grey), and healthy subjects (dotted black) versus CHF patients (dashed grey) in time series of 4000 samples.
* indicates statistically significant differences ( -values < 0.05 (dark), p-values < 0.07 (grey)). Data are shown as median and interquartile range. Table S2. P-values and AUC (Area Under Roc Curve) computed on synthetic data for the nonlinear indices. Mann-Whitney U test was used to evaluate statistical differences. Table S3. P-values computed in the comparison between ApEn(m, 0.2) and ApEn(m, (m)) across embedding dimensions. Mann-Whitney U test was used to evaluate statistical differences. Data are shown in terms of median (IQR). Table S4. P-values and AUC (Area Under Roc Curve) computed on HRV data for nonlinear indices. Mann-Whitney U test was used to evaluate statistical differences. Table S5. P-values and AUC (Area Under Roc Curve) computed on HRV data for temporal and spectral indices. Mann-Whitney U test was used to evaluate statistical differences. Figure S4. (a) RMSSD, (b) P HF , (c) P LF , and (d) P LFn values obtained per type of filtered white noise time series. * indicates statistically significant differences by Mann-Whitney U test ( -values < 0.05). Data are shown as median and interquartile range over 30 realizations. Table  S6. P-values from Mann-Whitney U test comparing the indices RMSSD, P HF , P LF , P LFn , and MApEn max over filtered white noise time series by considering BW1, BW2, and BW3, varying the number N of simulated processes. P-values < 0.05, marked in bold, are considered statistically significant. (Supplementary Materials)