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

Authors: Juan Bolea; Raquel Bailón; Esther Pueyo
Journal: Complexity
 (2018)

---

## Abstract

Background. Nonlinear heart rate variability (HRV) indices have extended the description of autonomic nervous system (ANS) regulation of the heart. One of those indices is approximate entropy,ApEn, which has become a commonly used measure of the irregularity of a time series. To calculateApEn, a priori definition of parameters like the threshold on similarity and the embedding dimension is required, which has been shown to be critical for interpretation of the results. Thus, searching for a parameter-freeApEn-based index could be advantageous for standardizing the use and interpretation of this widely applied entropy measurement.Methods. A novel entropy index called multidimensional approximate entropy, M
A
p
E
n
m
a
x, is proposed based on summing the contribution of maximum approximate entropies over a wide range of embedding dimensions while selecting the similarity threshold leading to maximumApEn value in each dimension. Synthetic RR interval time series with varying levels of stochasticity, generated by both MIX(P) processes and white/pink noise, were used to validate the properties of the proposed index. Aging and congestive heart failure (CHF) were characterized from RR interval time series of available databases.Results. In synthetic time series, M
A
p
E
n
m
a
x values were proportional to the level of randomness; i.e., M
A
p
E
n
m
a
x increased for higher values ofP in generated MIX(P) processes and was larger for white than for pink noise. This result was a consequence of all maximum approximate entropy values being increased for higher levels of randomness in all considered embedding dimensions. This is in contrast to the results obtained for approximate entropies computed with a fixed similarity threshold, which presented inconsistent results for different embedding dimensions. Evaluation of the proposed index on available databases revealed that aging was associated with a notable reduction in M
A
p
E
n
m
a
x values. On the other hand, M
A
p
E
n
m
a
x evaluated during the night period was considerably larger in CHF patients than in healthy subjects.Conclusion. A novel parameter-free multidimensional approximate entropy index, M
A
p
E
n
m
a
x, is proposed and tested over synthetic data to confirm its capacity to represent a range of randomness levels in HRV time series. M
A
p
E
n
m
a
x values are reduced in elderly patients, which may correspond to the reported loss of ANS adaptability in this population segment. Increased M
A
p
E
n
m
a
x values measured in CHF patients versus healthy subjects during the night period point to greater irregularity of heart rate dynamics caused by the disease.


---

## Body

## 1. 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 ofApEn reducing the bias of considering self-comparisons 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–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), andr, 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 the time series) [9–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 particularSampEn-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 refinedMSE approaches were introduced to take into account complexity properties at different scales by applying the coarse-graining technique [15, 20–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 parametersm andr. 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 thresholdris not strictly avoided, but this is replaced with the definition of a number of classes [26]. Rank-based methods were investigated to quantify the amount of shuffling by sorting the ranks of the mutual distances between pairs ofm-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 tolerancer 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 dimensionm 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.AlthoughSampEn was introduced as an improvement ofApEn entropy, advances towards the standardization ofApEn have been investigated as well. Different approaches have been proposed to identify the values of the embedding dimensionm or the thresholdr used for the characterization of nonlinear dynamics in chaotic time series [30–35]. With regard to the embedding dimension, the false nearest neighbor method was used to search for the lowest embedding dimensionm that allows phase-space reconstruction. Very different values ofm were reported for different models and experimental data [30]. The lowest values ofm were found to be 2, 3, and 6 for Hénon, Lorenz, and Mackey-Glass time series, respectively [31]. With regard to the tolerance thresholdr, some studies focused on searching for the value, denoted byr
max
(m), that maximizesApEn [32, 33, 36]. The correspondingApEn-based index, maximumApEn, was used to show that white noise series was more irregular than cross-chirp 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 fixedr 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 suchApEn profiles were shown to better separate young versus old individuals as compared to classicalApEnapproach. To reduce the influence of the thresholdr 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 althoughr 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 thresholdr itself. In the literature, identification of the values ofr andm 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 anApEn-basedindex that avoids the need for a priori parameter definition and, thus, allows standardization of entropy-based measurements. In particular, M
A
p
E
n
m
a
x is proposed based on summation of maximumApEn 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 forad hoc selection of input parameter values. The ability of the proposed index M
A
p
E
n
m
a
x 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).


## 2. 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 parameterP. MIX(P) generated a sine forP = 0 (pure deterministic) and became more random asP increased up to 1 (P = 1, pure stochastic):(1)
M
I
X
P
j
=
1
-
Z
j
X
j
+
Z
j
Y
j
,where X
j
=
2
sin
⁡
2
j
π
/
12, Y
j
~ are i.i.d. uniform random variables on [-
3
3], and Z
j
~ are i.i.d. random variables, with Z
j
=
1 having probabilityP and Z
j
=
0 probability 1-P [40]. For eachP value in {
0
,
0.25
,
0.5
,
0.75
,
1
}, 30 processes were generated, each containing 300 samples.Noise time series: (i) White and pink noise (1/f) time series with zero mean and unit variance were studied. A total of 30 white and pink noise processes, respectively, were generated, each containing 300 samples. (ii) Bandpass filtered noise derived from white noise time series filtered in two frequency passbands, representing the low (LF) and high (HF) frequency components, with variable bandwidths. Specifically, three types of bandpass filters with the following passband limits were considered: BW1 (LF = 0.09-0.13 Hz and HF = 0.22-0.28 Hz); BW2 (LF = 0.07-0.15 Hz and HF = 0.19-0.31 Hz); and BW3 (LF = 0.05-0.17 Hz and HF = 0.16-0.34 Hz). In all cases, the power of LF and HF spectral components, P
L
F (0.04-0.15 Hz) and P
H
F (0.15-0.4 Hz), remained invariant. A total of 30 filtered white noise time series per type were generated.


#### 2.1.1. Fantasia Database

Heartbeat locations from continuous electrocardiographic (ECG) signals sampled at 250 Hz were available from twenty young (21-34 y.o., mean 25.9) and 20 elderly rigorously screened healthy subjects (68-85 y.o., mean 74.5) undergoing 120 minutes of supine resting while watching the movie Fantasia (Disney, 1940) to help maintain wakefulness. Further database information is available in [41]. Recordings can be downloaded from http://www.physionet.org [42]. The whole 2-hour recordings were analyzed.


#### 2.1.2. Congestive Heart Failure and Normal Sinus Rhythm Databases

For the CHF cohort, heartbeat locations from 44 CHF patients, aged 55.5 ± 11.4 y.o., of the Congestive Heart Failure RR Interval Database and the BIDMC Congestive Heart Failure Database (including patients in all NYHA classes) were considered.For the healthy population, heartbeat locations from 72 healthy subjects, aged 54.6 ± 16 y.o., of the MIT-BIH Normal Sinus Rhythm Database and the Normal Sinus Rhythm RR Interval Database were considered.These databases were available fromhttp://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.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.


### 2.2. Methods

#### 2.2.1. Approximate Entropy

Letx
n
,
n
=
1
,
…
,
N, be the time series of interest, which is normalized by its standard deviation and y
i
m the ith reconstructed vector of embedding dimensionm.(2)
y
i
m
=
x
i
,
x
i
+
1
,
x
i
+
2
,
…
,
x
i
+
m
-
1
TThe amount of reconstructed vectors was N
m
=
N
-
m
-
1 for each embedding dimensionm. The distance, computed by the L
∞ norm, between each pair of reconstructed vectors, y
i
m
,
y
j
m, was denoted by d
i
,
j
m
.Each distance was compared with a thresholdr to compute the number of reconstructed vectors that lie within a hyperspace centered in the reconstructed vector of reference:(3)
C
i
m
r
=
1
N
m
∑
j
=
1
N
m
H
r
-
d
i
,
j
m
,where C
i
m is the correlation sum andH is the Heaviside function.This procedure was repeated with all reconstructed vectors and the probability of a pattern of lengthm appearing along the time series was denoted by the following.(4)
ϕ
m
r
=
1
N
m
∑
i
=
1
N
m
log
⁡
C
i
m
r
(5)
A
p
E
n
m
,
r
=
ϕ
m
r
-
ϕ
m
+
1
rIn this study,m = 2 andr = 0.2 were considered, leading to computation ofApEn(2, 0.2).


#### 2.2.2. Maximum Approximate Entropy

For each embedding dimensionm, the thresholdr leading to maximum approximate entropy, denoted byr
max
(m), was determined as in [45] and the maximum valueApEn(m, r
max
(m)) was calculated [36]. For the sake of comparison withApEn(2, 0.2), ApEn(2, r
max
(
2)) was calculated.


#### 2.2.3. Multidimensional Approximate Entropy

The multidimensional entropy index,MApEn(r), was introduced by summing up the approximate entropies for a range of embedding dimensions:(6)
M
A
p
E
n
r
=
∑
m
=
1
m
=
m
m
a
x
A
p
E
n
m
,
r
,wherem
max = 15. To additionally remove the dependence of the tolerance threshold, the valuer
max
(m) was identified for each embedding dimension m, leading to the definition of M
A
p
E
n
m
a
x.(7)
M
A
p
E
n
m
a
x
=
∑
m
=
1
m
=
m
m
a
x
A
p
E
n
m
,
r
m
a
x
m
ApEn,SampEn, and correlation dimensionD
2 were computed using the methodology described in [45], considering a range for the tolerance thresholdr varying from 0.01 to 3 times the standard deviation of the time series witha resolution of 0.01.


#### 2.2.4. 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].


### 2.3. 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 indexM
A
p
E
n
m
a
x and other HRV indices was computed.


## 2.1. Materials

MIX processes: a family of processes combining deterministic and stochastic behavior were studied. The degree of stochasticity was controlled by parameterP. MIX(P) generated a sine forP = 0 (pure deterministic) and became more random asP increased up to 1 (P = 1, pure stochastic):(1)
M
I
X
P
j
=
1
-
Z
j
X
j
+
Z
j
Y
j
,where X
j
=
2
sin
⁡
2
j
π
/
12, Y
j
~ are i.i.d. uniform random variables on [-
3
3], and Z
j
~ are i.i.d. random variables, with Z
j
=
1 having probabilityP and Z
j
=
0 probability 1-P [40]. For eachP value in {
0
,
0.25
,
0.5
,
0.75
,
1
}, 30 processes were generated, each containing 300 samples.Noise time series: (i) White and pink noise (1/f) time series with zero mean and unit variance were studied. A total of 30 white and pink noise processes, respectively, were generated, each containing 300 samples. (ii) Bandpass filtered noise derived from white noise time series filtered in two frequency passbands, representing the low (LF) and high (HF) frequency components, with variable bandwidths. Specifically, three types of bandpass filters with the following passband limits were considered: BW1 (LF = 0.09-0.13 Hz and HF = 0.22-0.28 Hz); BW2 (LF = 0.07-0.15 Hz and HF = 0.19-0.31 Hz); and BW3 (LF = 0.05-0.17 Hz and HF = 0.16-0.34 Hz). In all cases, the power of LF and HF spectral components, P
L
F (0.04-0.15 Hz) and P
H
F (0.15-0.4 Hz), remained invariant. A total of 30 filtered white noise time series per type were generated.


### 2.1.1. Fantasia Database

Heartbeat locations from continuous electrocardiographic (ECG) signals sampled at 250 Hz were available from twenty young (21-34 y.o., mean 25.9) and 20 elderly rigorously screened healthy subjects (68-85 y.o., mean 74.5) undergoing 120 minutes of supine resting while watching the movie Fantasia (Disney, 1940) to help maintain wakefulness. Further database information is available in [41]. Recordings can be downloaded from http://www.physionet.org [42]. The whole 2-hour recordings were analyzed.


### 2.1.2. Congestive Heart Failure and Normal Sinus Rhythm Databases

For the CHF cohort, heartbeat locations from 44 CHF patients, aged 55.5 ± 11.4 y.o., of the Congestive Heart Failure RR Interval Database and the BIDMC Congestive Heart Failure Database (including patients in all NYHA classes) were considered.For the healthy population, heartbeat locations from 72 healthy subjects, aged 54.6 ± 16 y.o., of the MIT-BIH Normal Sinus Rhythm Database and the Normal Sinus Rhythm RR Interval Database were considered.These databases were available fromhttp://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.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.


## 2.1.1. Fantasia Database

Heartbeat locations from continuous electrocardiographic (ECG) signals sampled at 250 Hz were available from twenty young (21-34 y.o., mean 25.9) and 20 elderly rigorously screened healthy subjects (68-85 y.o., mean 74.5) undergoing 120 minutes of supine resting while watching the movie Fantasia (Disney, 1940) to help maintain wakefulness. Further database information is available in [41]. Recordings can be downloaded from http://www.physionet.org [42]. The whole 2-hour recordings were analyzed.


## 2.1.2. Congestive Heart Failure and Normal Sinus Rhythm Databases

For the CHF cohort, heartbeat locations from 44 CHF patients, aged 55.5 ± 11.4 y.o., of the Congestive Heart Failure RR Interval Database and the BIDMC Congestive Heart Failure Database (including patients in all NYHA classes) were considered.For the healthy population, heartbeat locations from 72 healthy subjects, aged 54.6 ± 16 y.o., of the MIT-BIH Normal Sinus Rhythm Database and the Normal Sinus Rhythm RR Interval Database were considered.These databases were available fromhttp://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.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.


## 2.2. Methods

### 2.2.1. Approximate Entropy

Letx
n
,
n
=
1
,
…
,
N, be the time series of interest, which is normalized by its standard deviation and y
i
m the ith reconstructed vector of embedding dimensionm.(2)
y
i
m
=
x
i
,
x
i
+
1
,
x
i
+
2
,
…
,
x
i
+
m
-
1
TThe amount of reconstructed vectors was N
m
=
N
-
m
-
1 for each embedding dimensionm. The distance, computed by the L
∞ norm, between each pair of reconstructed vectors, y
i
m
,
y
j
m, was denoted by d
i
,
j
m
.Each distance was compared with a thresholdr to compute the number of reconstructed vectors that lie within a hyperspace centered in the reconstructed vector of reference:(3)
C
i
m
r
=
1
N
m
∑
j
=
1
N
m
H
r
-
d
i
,
j
m
,where C
i
m is the correlation sum andH is the Heaviside function.This procedure was repeated with all reconstructed vectors and the probability of a pattern of lengthm appearing along the time series was denoted by the following.(4)
ϕ
m
r
=
1
N
m
∑
i
=
1
N
m
log
⁡
C
i
m
r
(5)
A
p
E
n
m
,
r
=
ϕ
m
r
-
ϕ
m
+
1
rIn this study,m = 2 andr = 0.2 were considered, leading to computation ofApEn(2, 0.2).


### 2.2.2. Maximum Approximate Entropy

For each embedding dimensionm, the thresholdr leading to maximum approximate entropy, denoted byr
max
(m), was determined as in [45] and the maximum valueApEn(m, r
max
(m)) was calculated [36]. For the sake of comparison withApEn(2, 0.2), ApEn(2, r
max
(
2)) was calculated.


### 2.2.3. Multidimensional Approximate Entropy

The multidimensional entropy index,MApEn(r), was introduced by summing up the approximate entropies for a range of embedding dimensions:(6)
M
A
p
E
n
r
=
∑
m
=
1
m
=
m
m
a
x
A
p
E
n
m
,
r
,wherem
max = 15. To additionally remove the dependence of the tolerance threshold, the valuer
max
(m) was identified for each embedding dimension m, leading to the definition of M
A
p
E
n
m
a
x.(7)
M
A
p
E
n
m
a
x
=
∑
m
=
1
m
=
m
m
a
x
A
p
E
n
m
,
r
m
a
x
m
ApEn,SampEn, and correlation dimensionD
2 were computed using the methodology described in [45], considering a range for the tolerance thresholdr varying from 0.01 to 3 times the standard deviation of the time series witha resolution of 0.01.


### 2.2.4. 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].


## 2.2.1. Approximate Entropy

Letx
n
,
n
=
1
,
…
,
N, be the time series of interest, which is normalized by its standard deviation and y
i
m the ith reconstructed vector of embedding dimensionm.(2)
y
i
m
=
x
i
,
x
i
+
1
,
x
i
+
2
,
…
,
x
i
+
m
-
1
TThe amount of reconstructed vectors was N
m
=
N
-
m
-
1 for each embedding dimensionm. The distance, computed by the L
∞ norm, between each pair of reconstructed vectors, y
i
m
,
y
j
m, was denoted by d
i
,
j
m
.Each distance was compared with a thresholdr to compute the number of reconstructed vectors that lie within a hyperspace centered in the reconstructed vector of reference:(3)
C
i
m
r
=
1
N
m
∑
j
=
1
N
m
H
r
-
d
i
,
j
m
,where C
i
m is the correlation sum andH is the Heaviside function.This procedure was repeated with all reconstructed vectors and the probability of a pattern of lengthm appearing along the time series was denoted by the following.(4)
ϕ
m
r
=
1
N
m
∑
i
=
1
N
m
log
⁡
C
i
m
r
(5)
A
p
E
n
m
,
r
=
ϕ
m
r
-
ϕ
m
+
1
rIn this study,m = 2 andr = 0.2 were considered, leading to computation ofApEn(2, 0.2).


## 2.2.2. Maximum Approximate Entropy

For each embedding dimensionm, the thresholdr leading to maximum approximate entropy, denoted byr
max
(m), was determined as in [45] and the maximum valueApEn(m, r
max
(m)) was calculated [36]. For the sake of comparison withApEn(2, 0.2), ApEn(2, r
max
(
2)) was calculated.


## 2.2.3. Multidimensional Approximate Entropy

The multidimensional entropy index,MApEn(r), was introduced by summing up the approximate entropies for a range of embedding dimensions:(6)
M
A
p
E
n
r
=
∑
m
=
1
m
=
m
m
a
x
A
p
E
n
m
,
r
,wherem
max = 15. To additionally remove the dependence of the tolerance threshold, the valuer
max
(m) was identified for each embedding dimension m, leading to the definition of M
A
p
E
n
m
a
x.(7)
M
A
p
E
n
m
a
x
=
∑
m
=
1
m
=
m
m
a
x
A
p
E
n
m
,
r
m
a
x
m
ApEn,SampEn, and correlation dimensionD
2 were computed using the methodology described in [45], considering a range for the tolerance thresholdr varying from 0.01 to 3 times the standard deviation of the time series witha resolution of 0.01.


## 2.2.4. 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].


## 2.3. 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 indexM
A
p
E
n
m
a
x and other HRV indices was computed.


## 3. Results

### 3.1. Synthetic Data

Figures1 and 2 show the values of the analyzed nonlinear HRV indices computed over MIX(P) processes and pink and white noise time series. In the comparison of MIX(P) processes with randomness degree belowP = 0.75 (0 versus 0.25; 0.25 versus 0.5; and 0.5 versus 0.75) all nonlinear indices led to statistically significant differences. In the comparison between MIX(P) processes withP above 0.75 (0.75 versus 1) no statistically significant differences were found for any of the nonlinear indices.ApEn(2, 0.2)andApEn(2,r
max
(
2)) as well as nonlinear index values were higher for greater levels of randomness until MIX(0.75); however aboveP = 0.75 an increase in the level of randomness was not captured byApEn 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 ofr
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 indexMApEn(r) was able to separate pink and white noise time series for all tested values of the tolerance thresholdr(results only shown forr= 0.2, Figure 2(a)).Figure 1
Nonlinear indices: (a)ApEn(2,0.2); (b)SampEn(2,0.2); (c)ApEn(2,r
max
(
2)); and (d)D
2 evaluated over MIX(P) processes with varying stochastic levels and pink and white noise time series. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) when comparing MIX(P) processes or pink versus white noise time series. Data are shown as median and interquartile range over 30 realizations.
(a)
(b)
(c)
(d)Figure 2
(a)MApEn(0.2)and (b) M
A
p
E
n
m
a
x evaluated over MIX(P) processes with varying stochastic levels as well as white and pink noise. (∗) indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) when comparing MIX(P) processes or pink versus white noise time series. Data are shown as median and interquartile range over 30 realizations.
(a)
(b)While statistically significant differences were found for bothMApEn(0.2) and M
A
p
E
n
m
a
x 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 evaluatingApEn for two embedding dimensions, namely,m= 2 andm = 6, are illustrated in Figure 3. The use of a fixed threshold did not always allowApEn to characterize the randomness level of the MIX(P) processes. As an example, form = 6, the use of the fixed thresholdr = 0.2 led to MIX(0.25) having much larger complexity than MIX(0.5), MIX(0.75), and MIX(1). Form = 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 ofApEn(m, r
max
(m)) when computing M
A
p
E
n
m
a
x served to characterize MIX(P) processes according to their stochasticity level for any value ofm, which then rendered M
A
p
E
n
m
a
x suitable to characterize stochasticity levels (Figure 2(b)).Figure 3
Approximate entropy values computed over MIX(P) processes considering a set of thresholds and evaluated form = 2 (a) andm= 6 (b) for a particular stochastic realization. Dashed line indicatesr = 0.2 in both panels.ApEn(2,r
max
(
2)),ApEn(6,r
max
(
6)),r
max
(
2), andr
max
(
6)values are displayed for each case.
(a)
(b)In Figures4(a) and 4(b)
ApEn(m, 0.2)andApEn(m, r
max
(m)) are shown for pink and white noise when the embedding dimensionm 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 indexMApEn(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 M
A
p
E
n
m
a
x (Figure 2(b)) with respect toMApEn(0.2)(Figure 2(a)).Figure 4
ApEn computed across embedding dimensions for pink noise (grey) and white noise (black) using different thresholds (a) r(m) = 0.2 and (b) r(m) = r
m
a
x(m). ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between pink and white noise for a given embedding dimension. Data are shown as median and interquartile range over 30 realizations.
(a)
(b)As the statistical significance in the comparison of simulated data highly depends on the number of simulated processes (30 in the comparisons shown in Figure2(b)), a study was performed to characterize such a dependence for the proposed index M
A
p
E
n
m
a
x 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 ap-value 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.Table 1
P-values from Mann-Whitney U test comparing the index M
A
p
E
n
m
a
x over different synthetic time series when varying the number N of simulated stochastic processes.
M
A
p
E
n
m
a
x
N = 5
N = 10
N = 15
N = 20
N = 25
N = 30
MIX(0) vs. MIX(0.25)
0.0079
6.4e-05
6.9e-07
8.0e-09
9.7e-11
1.2e-12
MIX(0.25) vs. MIX(0.5)
0.0079
1.8e-04
3.4e-06
6.8e-08
1.4e-09
3.02e-11
MIX(0.5) vs. MIX(0.75)
0.0079
0.001
4.02e-05
1.2e-06
4.9e-08
4.9e-09
MIX(0.75) vs. MIX(1)
0.55
0.43
0.80
0.32
0.17
0.25
Pink vs. White noise
0.0079
1.8e-04
3.4e-06
6.8e-08
1.4e-09
3.02e-11The median spectrum obtained from each type of bandpass filtered noise time series is illustrated in Figure5(a). Time- and frequency-domain indices, as well as the proposed M
A
p
E
n
m
a
x index, were estimated over these filtered noise time series (see Figure S4 in Supplementary Materials for time- and frequency-domain indices). Statistically significant differences when comparing the different filtered noise time series were found only for M
A
p
E
n
m
a
x by Mann-Whitney U test (Figure 5(b)). Additionally, the dependence of the statistical significance on the number of realizations was studied forRMSSD,P
HF
, P
LF
, P
LFn
, and M
A
p
E
n
m
a
x (see Table S6 in Supplementary Materials).Figure 5
(a) Median spectrum and (b)M
A
p
E
n
m
a
x values obtained per type of filtered white noise time series. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05). Data are shown as median and interquartile range over 30 realizations.
(a)
(b)


### 3.2. 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),
M
A
p
E
n
m
a
x
, SampEn(2, 0.2), andD
2 values showed a trend for higher values in young than in elderly subjects, but only M
A
p
E
n
m
a
x 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 ofMApEn(0.2), were statistically significantly larger in CHF patients than in healthy subjects (Figures 6 and 7).Figure 6
Regularity and complexity indices: (a)ApEn(2,0.2); (b)SampEn(2,0.2); (c)ApEn(2,r
max
(
2)); and (d)D
2 comparing young versus elderly, and healthy subjects versus CHF patients. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups. Data are shown as median and interquartile range.
(a)
(b)
(c)
(d)Figure 7
(a)MApEn(0.2)and (b) M
A
p
E
n
m
a
x values computed for young (solid black) versus elderly (solid grey), and healthy subjects (dotted) versus CHF patients (dashed). ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups. Data are shown as median and interquartile range.
(a)
(b)Figure8 presentsApEn(m,0.2) andApEn(m, r
max
(m)), as contributors toMApEn(0.2) and M
A
p
E
n
m
a
x (Figures 6(a) and 6(b)) separately for each of the analyzed cohorts. The indexApEn(m, 0.2) was found to have higher values in young subjects than in elderly subjects for low embedding dimensions (m <= 3), even if not achieving statistical significance in all cases (Figure 8(a)). Similarly, CHF patients showed higherApEn(m, 0.2) values than healthy subjects for low values of the embedding dimensionm (Figure 8(c)). For large embedding dimensions, opposite results were found in the two comparisons. Similar results were found forApEn(m, 0.1) andApEn(m, 0.15) (see Figure S1 in Supplementary Materials, where additional results are presented for other commonly usedr values). RegardingApEn(m, r
max
(m)), statistically significantly higher values were found in young versus elderly and in CHF patients versus healthy subjects for most values of the embedding dimensionm, while for the remaining values ofmthe compared values could not be statistically distinguished.Figure 8
Approximate entropy values computed for thresholds set at 0.2 andr
max
(m) across embedding dimensions in young (black) and elderly (grey) groups shown in (a) and (b), and in healthy subjects (dotted grey) and CHF patients (dashed black) during the night period in (c) and (d). ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups for a given embedding dimension. Data are shown as median and interquartile range.
(a)
(b)
(c)
(d)HRM andRMSSD values are displayed in Figure 9 for the studied groups. No statistically significant differences were found forHRM when comparing young and elderly subjects, although a tendency to lowerHRM values in the elderly group was observed.HRM was statistically significantly higher in CHF patients as compared to healthy subjects. On the other hand, theRMSSD index showed statistically significant differences between both pairs of groups, with elderly subjects and CHF patients presenting reducedRMSSD values.Figure 9
(a)HRM and (b)RMSSD values in young versus elderly and CHF patients versus healthy subjects during the night period. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups.
(a)
(b)Values of the frequency-domain HRV indicesP
LF,P
HF, andP
LFn are shown in Figure 10 for the four studied groups. Significant differences were found forP
LF andP
HF in the two comparisons.P
LFn only showed statistically significant differences for the comparison between CHF patients and healthy subjects.Figure 10
Frequency-domain indices: (a)P
LF; (b)P
HF; and (c)P
LFn computed comparing young and elderly, and CHF patients and healthy subjects during the night period. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups.
(a)
(b)
(c)Figure11 illustrates the relationship between M
A
p
E
n
m
a
x andRMSSD (a) andP
HF (b). Correlation coefficients computed over the four analyzed groups are shown in Table 2.Table 2
Spearman correlation coefficients betweenM
A
p
E
n
m
a
x and RMSSD and between M
A
p
E
n
m
a
x and P
HF.
M
A
p
E
n
m
a
x vs. RMSSD
M
A
p
E
n
m
a
x vs. P
HF
ρ
p-value
ρ
p-value
Young
0.02
0.93
0.21
0.37
Elderly
-0.04
0.86
0.08
0.73
Healthy
-0.13
0.26
0.023
0.85
CHF
-0.57
0.0001
-0.65
0.0001Figure 11
Scatter plots showing the relationship betweenM
A
p
E
n
m
a
x andRMSSD (a) and between M
A
p
E
n
m
a
x andP
HF (b).
(a)
(b)


## 3.1. Synthetic Data

Figures1 and 2 show the values of the analyzed nonlinear HRV indices computed over MIX(P) processes and pink and white noise time series. In the comparison of MIX(P) processes with randomness degree belowP = 0.75 (0 versus 0.25; 0.25 versus 0.5; and 0.5 versus 0.75) all nonlinear indices led to statistically significant differences. In the comparison between MIX(P) processes withP above 0.75 (0.75 versus 1) no statistically significant differences were found for any of the nonlinear indices.ApEn(2, 0.2)andApEn(2,r
max
(
2)) as well as nonlinear index values were higher for greater levels of randomness until MIX(0.75); however aboveP = 0.75 an increase in the level of randomness was not captured byApEn 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 ofr
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 indexMApEn(r) was able to separate pink and white noise time series for all tested values of the tolerance thresholdr(results only shown forr= 0.2, Figure 2(a)).Figure 1
Nonlinear indices: (a)ApEn(2,0.2); (b)SampEn(2,0.2); (c)ApEn(2,r
max
(
2)); and (d)D
2 evaluated over MIX(P) processes with varying stochastic levels and pink and white noise time series. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) when comparing MIX(P) processes or pink versus white noise time series. Data are shown as median and interquartile range over 30 realizations.
(a)
(b)
(c)
(d)Figure 2
(a)MApEn(0.2)and (b) M
A
p
E
n
m
a
x evaluated over MIX(P) processes with varying stochastic levels as well as white and pink noise. (∗) indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) when comparing MIX(P) processes or pink versus white noise time series. Data are shown as median and interquartile range over 30 realizations.
(a)
(b)While statistically significant differences were found for bothMApEn(0.2) and M
A
p
E
n
m
a
x 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 evaluatingApEn for two embedding dimensions, namely,m= 2 andm = 6, are illustrated in Figure 3. The use of a fixed threshold did not always allowApEn to characterize the randomness level of the MIX(P) processes. As an example, form = 6, the use of the fixed thresholdr = 0.2 led to MIX(0.25) having much larger complexity than MIX(0.5), MIX(0.75), and MIX(1). Form = 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 ofApEn(m, r
max
(m)) when computing M
A
p
E
n
m
a
x served to characterize MIX(P) processes according to their stochasticity level for any value ofm, which then rendered M
A
p
E
n
m
a
x suitable to characterize stochasticity levels (Figure 2(b)).Figure 3
Approximate entropy values computed over MIX(P) processes considering a set of thresholds and evaluated form = 2 (a) andm= 6 (b) for a particular stochastic realization. Dashed line indicatesr = 0.2 in both panels.ApEn(2,r
max
(
2)),ApEn(6,r
max
(
6)),r
max
(
2), andr
max
(
6)values are displayed for each case.
(a)
(b)In Figures4(a) and 4(b)
ApEn(m, 0.2)andApEn(m, r
max
(m)) are shown for pink and white noise when the embedding dimensionm 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 indexMApEn(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 M
A
p
E
n
m
a
x (Figure 2(b)) with respect toMApEn(0.2)(Figure 2(a)).Figure 4
ApEn computed across embedding dimensions for pink noise (grey) and white noise (black) using different thresholds (a) r(m) = 0.2 and (b) r(m) = r
m
a
x(m). ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between pink and white noise for a given embedding dimension. Data are shown as median and interquartile range over 30 realizations.
(a)
(b)As the statistical significance in the comparison of simulated data highly depends on the number of simulated processes (30 in the comparisons shown in Figure2(b)), a study was performed to characterize such a dependence for the proposed index M
A
p
E
n
m
a
x 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 ap-value 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.Table 1
P-values from Mann-Whitney U test comparing the index M
A
p
E
n
m
a
x over different synthetic time series when varying the number N of simulated stochastic processes.
M
A
p
E
n
m
a
x
N = 5
N = 10
N = 15
N = 20
N = 25
N = 30
MIX(0) vs. MIX(0.25)
0.0079
6.4e-05
6.9e-07
8.0e-09
9.7e-11
1.2e-12
MIX(0.25) vs. MIX(0.5)
0.0079
1.8e-04
3.4e-06
6.8e-08
1.4e-09
3.02e-11
MIX(0.5) vs. MIX(0.75)
0.0079
0.001
4.02e-05
1.2e-06
4.9e-08
4.9e-09
MIX(0.75) vs. MIX(1)
0.55
0.43
0.80
0.32
0.17
0.25
Pink vs. White noise
0.0079
1.8e-04
3.4e-06
6.8e-08
1.4e-09
3.02e-11The median spectrum obtained from each type of bandpass filtered noise time series is illustrated in Figure5(a). Time- and frequency-domain indices, as well as the proposed M
A
p
E
n
m
a
x index, were estimated over these filtered noise time series (see Figure S4 in Supplementary Materials for time- and frequency-domain indices). Statistically significant differences when comparing the different filtered noise time series were found only for M
A
p
E
n
m
a
x by Mann-Whitney U test (Figure 5(b)). Additionally, the dependence of the statistical significance on the number of realizations was studied forRMSSD,P
HF
, P
LF
, P
LFn
, and M
A
p
E
n
m
a
x (see Table S6 in Supplementary Materials).Figure 5
(a) Median spectrum and (b)M
A
p
E
n
m
a
x values obtained per type of filtered white noise time series. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05). Data are shown as median and interquartile range over 30 realizations.
(a)
(b)


## 3.2. 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),
M
A
p
E
n
m
a
x
, SampEn(2, 0.2), andD
2 values showed a trend for higher values in young than in elderly subjects, but only M
A
p
E
n
m
a
x 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 ofMApEn(0.2), were statistically significantly larger in CHF patients than in healthy subjects (Figures 6 and 7).Figure 6
Regularity and complexity indices: (a)ApEn(2,0.2); (b)SampEn(2,0.2); (c)ApEn(2,r
max
(
2)); and (d)D
2 comparing young versus elderly, and healthy subjects versus CHF patients. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups. Data are shown as median and interquartile range.
(a)
(b)
(c)
(d)Figure 7
(a)MApEn(0.2)and (b) M
A
p
E
n
m
a
x values computed for young (solid black) versus elderly (solid grey), and healthy subjects (dotted) versus CHF patients (dashed). ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups. Data are shown as median and interquartile range.
(a)
(b)Figure8 presentsApEn(m,0.2) andApEn(m, r
max
(m)), as contributors toMApEn(0.2) and M
A
p
E
n
m
a
x (Figures 6(a) and 6(b)) separately for each of the analyzed cohorts. The indexApEn(m, 0.2) was found to have higher values in young subjects than in elderly subjects for low embedding dimensions (m <= 3), even if not achieving statistical significance in all cases (Figure 8(a)). Similarly, CHF patients showed higherApEn(m, 0.2) values than healthy subjects for low values of the embedding dimensionm (Figure 8(c)). For large embedding dimensions, opposite results were found in the two comparisons. Similar results were found forApEn(m, 0.1) andApEn(m, 0.15) (see Figure S1 in Supplementary Materials, where additional results are presented for other commonly usedr values). RegardingApEn(m, r
max
(m)), statistically significantly higher values were found in young versus elderly and in CHF patients versus healthy subjects for most values of the embedding dimensionm, while for the remaining values ofmthe compared values could not be statistically distinguished.Figure 8
Approximate entropy values computed for thresholds set at 0.2 andr
max
(m) across embedding dimensions in young (black) and elderly (grey) groups shown in (a) and (b), and in healthy subjects (dotted grey) and CHF patients (dashed black) during the night period in (c) and (d). ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups for a given embedding dimension. Data are shown as median and interquartile range.
(a)
(b)
(c)
(d)HRM andRMSSD values are displayed in Figure 9 for the studied groups. No statistically significant differences were found forHRM when comparing young and elderly subjects, although a tendency to lowerHRM values in the elderly group was observed.HRM was statistically significantly higher in CHF patients as compared to healthy subjects. On the other hand, theRMSSD index showed statistically significant differences between both pairs of groups, with elderly subjects and CHF patients presenting reducedRMSSD values.Figure 9
(a)HRM and (b)RMSSD values in young versus elderly and CHF patients versus healthy subjects during the night period. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups.
(a)
(b)Values of the frequency-domain HRV indicesP
LF,P
HF, andP
LFn are shown in Figure 10 for the four studied groups. Significant differences were found forP
LF andP
HF in the two comparisons.P
LFn only showed statistically significant differences for the comparison between CHF patients and healthy subjects.Figure 10
Frequency-domain indices: (a)P
LF; (b)P
HF; and (c)P
LFn computed comparing young and elderly, and CHF patients and healthy subjects during the night period. ∗ indicates statistically significant differences by Mann-Whitney U test (p-values < 0.05) between groups.
(a)
(b)
(c)Figure11 illustrates the relationship between M
A
p
E
n
m
a
x andRMSSD (a) andP
HF (b). Correlation coefficients computed over the four analyzed groups are shown in Table 2.Table 2
Spearman correlation coefficients betweenM
A
p
E
n
m
a
x and RMSSD and between M
A
p
E
n
m
a
x and P
HF.
M
A
p
E
n
m
a
x vs. RMSSD
M
A
p
E
n
m
a
x vs. P
HF
ρ
p-value
ρ
p-value
Young
0.02
0.93
0.21
0.37
Elderly
-0.04
0.86
0.08
0.73
Healthy
-0.13
0.26
0.023
0.85
CHF
-0.57
0.0001
-0.65
0.0001Figure 11
Scatter plots showing the relationship betweenM
A
p
E
n
m
a
x andRMSSD (a) and between M
A
p
E
n
m
a
x andP
HF (b).
(a)
(b)


## 4. Discussion

### 4.1. 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–9]. Nevertheless, the physiological interpretation of some common nonlinear HRV indices could be biased by the a priori selection of parameter values intrinsic to their definitions. Some investigations have addressed the automatic search for those parameter values based on specific time series characteristics [32, 35, 38]. Indices likeTotalApEn and maximum approximate entropy are still dependent on the definition of the embedding dimensionm. The values ofm vary from one study to another even though there is a general tendency to setm = 2 to estimateApEn(2, r
max
(
2)).In this study, the index M
A
p
E
n
m
a
x has been proposed based on summing the maximal approximate entropy values over all analyzed embedding dimensions. Although the computation of M
A
p
E
n
m
a
x requires setting a maximum embedding dimension value, the relative contribution ofApEn(m,r
max
(m)) values to the final M
A
p
E
n
m
a
x estimation largely decreases asm increases, thus rendering M
A
p
E
n
m
a
x practically independent of such selection. In fact, in this studyApEn(m,r
max
(m)) was found to present a remarkable decay form 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, M
A
p
E
n
m
a
x, can be seen as an entropy-based index not requiring a priori parameter definition for its computation. The algorithm described in [45] forD
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 300-sample 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 Windows® 7 based PC, Intel® Core™ i7 3.5GHz, 16Gb RAM with Matlab® R2015a.OtherSampEn-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 theseSampEn-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.


### 4.2. Analysis of Synthetic Data

The proposed indexM
A
p
E
n
m
a
x has been tested for characterization of synthetic time series. Statistically significantly larger values of M
A
p
E
n
m
a
x were found for white than for pink noise. This behavior was also reproduced by all other tested nonlinear HRV indices includingMApEn(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 M
A
p
E
n
m
a
x, which is in agreement with the results shown by other nonlinear HRV indices, likeD
2,ApEn, andSampEn, describing complexity of the time series. On the contrary, the use of a fixed thresholdr, 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 time- and frequency-domain indices were statistically similar. The analysis showed thatM
A
p
E
n
m
a
x provides complementary information to that from time- and frequency-domain indices. According to the simulation results, M
A
p
E
n
m
a
x was sensitive to changes in spectrum shape involving changes in HRV regularity. The narrower the spectral components, the lower the M
A
p
E
n
m
a
x 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 thep-values when the proposed index M
A
p
E
n
m
a
x 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 thep-value of any statistical comparison over synthetically generated data.


### 4.3. 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 proposedM
A
p
E
n
m
a
x index. Our results are in agreement with previous studies where aging has been reported as a cause for decreased complexity and irregularity of the beat-to-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, andP
HF values were also decreased, whileP
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 proposedM
A
p
E
n
m
a
x 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 inP
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 inApEn and decreases inP
LF andP
LFn, while sympathetic activation decreases the values of all nonlinear indices and increasesP
LFn [53, 54]. Thus, an increase inP
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 inP
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 byApEn in the literature, but results have been found to be dependent on the applied method [56]. For instance, the use of a point-process 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 classicalSampEn 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 greaterMSE values than CHF patients, but in particular the elderly healthy cohort showed lowerMSE values in the first scale than the CHF group. Although a relationship betweenMSE-derived indices and sympathetic and parasympathetic modulation has been reported [57], the methodological differences comparingMSE-derived indices and M
A
p
E
n
m
a
x do not allow translation of the physiological interpretation to the present study.In our study CHF patients and healthy subjects showed diverseApEn values, being greater in one or the other group depending on the tolerancer and the embedding dimensionm values. Therefore, the different results obtained in the present study and in [21] could be due to the method itself (ApEnhere andSampEn 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 lowerMSE 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 M
A
p
E
n
m
a
x was also tested for longer segments of HRV analysis (up to 4000 samples rather than 300 samples). Results confirmed that M
A
p
E
n
m
a
x 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 M
A
p
E
n
m
a
x 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 M
A
p
E
n
m
a
x in clinical and experimental applications should be the focus of future research studies.


## 4.1. 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–9]. Nevertheless, the physiological interpretation of some common nonlinear HRV indices could be biased by the a priori selection of parameter values intrinsic to their definitions. Some investigations have addressed the automatic search for those parameter values based on specific time series characteristics [32, 35, 38]. Indices likeTotalApEn and maximum approximate entropy are still dependent on the definition of the embedding dimensionm. The values ofm vary from one study to another even though there is a general tendency to setm = 2 to estimateApEn(2, r
max
(
2)).In this study, the index M
A
p
E
n
m
a
x has been proposed based on summing the maximal approximate entropy values over all analyzed embedding dimensions. Although the computation of M
A
p
E
n
m
a
x requires setting a maximum embedding dimension value, the relative contribution ofApEn(m,r
max
(m)) values to the final M
A
p
E
n
m
a
x estimation largely decreases asm increases, thus rendering M
A
p
E
n
m
a
x practically independent of such selection. In fact, in this studyApEn(m,r
max
(m)) was found to present a remarkable decay form 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, M
A
p
E
n
m
a
x, can be seen as an entropy-based index not requiring a priori parameter definition for its computation. The algorithm described in [45] forD
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 300-sample 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 Windows® 7 based PC, Intel® Core™ i7 3.5GHz, 16Gb RAM with Matlab® R2015a.OtherSampEn-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 theseSampEn-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.


## 4.2. Analysis of Synthetic Data

The proposed indexM
A
p
E
n
m
a
x has been tested for characterization of synthetic time series. Statistically significantly larger values of M
A
p
E
n
m
a
x were found for white than for pink noise. This behavior was also reproduced by all other tested nonlinear HRV indices includingMApEn(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 M
A
p
E
n
m
a
x, which is in agreement with the results shown by other nonlinear HRV indices, likeD
2,ApEn, andSampEn, describing complexity of the time series. On the contrary, the use of a fixed thresholdr, 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 time- and frequency-domain indices were statistically similar. The analysis showed thatM
A
p
E
n
m
a
x provides complementary information to that from time- and frequency-domain indices. According to the simulation results, M
A
p
E
n
m
a
x was sensitive to changes in spectrum shape involving changes in HRV regularity. The narrower the spectral components, the lower the M
A
p
E
n
m
a
x 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 thep-values when the proposed index M
A
p
E
n
m
a
x 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 thep-value of any statistical comparison over synthetically generated data.


## 4.3. 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 proposedM
A
p
E
n
m
a
x index. Our results are in agreement with previous studies where aging has been reported as a cause for decreased complexity and irregularity of the beat-to-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, andP
HF values were also decreased, whileP
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 proposedM
A
p
E
n
m
a
x 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 inP
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 inApEn and decreases inP
LF andP
LFn, while sympathetic activation decreases the values of all nonlinear indices and increasesP
LFn [53, 54]. Thus, an increase inP
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 inP
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 byApEn in the literature, but results have been found to be dependent on the applied method [56]. For instance, the use of a point-process 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 classicalSampEn 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 greaterMSE values than CHF patients, but in particular the elderly healthy cohort showed lowerMSE values in the first scale than the CHF group. Although a relationship betweenMSE-derived indices and sympathetic and parasympathetic modulation has been reported [57], the methodological differences comparingMSE-derived indices and M
A
p
E
n
m
a
x do not allow translation of the physiological interpretation to the present study.In our study CHF patients and healthy subjects showed diverseApEn values, being greater in one or the other group depending on the tolerancer and the embedding dimensionm values. Therefore, the different results obtained in the present study and in [21] could be due to the method itself (ApEnhere andSampEn 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 lowerMSE 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 M
A
p
E
n
m
a
x was also tested for longer segments of HRV analysis (up to 4000 samples rather than 300 samples). Results confirmed that M
A
p
E
n
m
a
x 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 M
A
p
E
n
m
a
x 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 M
A
p
E
n
m
a
x in clinical and experimental applications should be the focus of future research studies.


## 5. Conclusions

A multidimensional approximate entropy index,M
A
p
E
n
m
a
x, 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 M
A
p
E
n
m
a
x consistently characterized the degree of randomness in the series, whereasApEn computed by using fixed thresholds produces inconsistent results depending on the analyzed embedding dimension. M
A
p
E
n
m
a
x provided additional information in the analysis of the effects of aging. Specifically, M
A
p
E
n
m
a
x 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, M
A
p
E
n
m
a
x 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, M
A
p
E
n
m
a
x can provide robust information on nonlinear HR dynamics with the advantage of being independent of preselected parameter values.



---

*Source: 4953273-2018-11-05.xml*