Influence of Heart Rate in Non-linear HRV Indices as a Sampling Rate Effect Evaluated on Supine and Standing

The purpose of this study is to characterize and attenuate the influence of mean heart rate (HR) on nonlinear heart rate variability (HRV) indices (correlation dimension, sample, and approximate entropy) as a consequence of being the HR the intrinsic sampling rate of HRV signal. This influence can notably alter nonlinear HRV indices and lead to biased information regarding autonomic nervous system (ANS) modulation. First, a simulation study was carried out to characterize the dependence of nonlinear HRV indices on HR assuming similar ANS modulation. Second, two HR-correction approaches were proposed: one based on regression formulas and another one based on interpolating RR time series. Finally, standard and HR-corrected HRV indices were studied in a body position change database. The simulation study showed the HR-dependence of non-linear indices as a sampling rate effect, as well as the ability of the proposed HR-corrections to attenuate mean HR influence. Analysis in a body position changes database shows that correlation dimension was reduced around 21% in median values in standing with respect to supine position (p < 0.05), concomitant with a 28% increase in mean HR (p < 0.05). After HR-correction, correlation dimension decreased around 18% in standing with respect to supine position, being the decrease still significant. Sample and approximate entropy showed similar trends. HR-corrected nonlinear HRV indices could represent an improvement in their applicability as markers of ANS modulation when mean HR changes.

Introduction

Heart rate (HR) variability (HRV) has been studied as a non-invasive technique to assess autonomic nervous system (ANS) regulation of the heart. Although, HRV analysis is still controversial (Karemaker, 2006), its content has been related to sympathetic and parasympathetic modulation by Task Force of the ESC-NASPE (1996) and Sassi et al. (2015).

During the last decades, HRV analysis has been extended including nonlinear indices based on chaos theory. These methodologies describe ANS in terms of regularity and complexity. Non-linear indices have been studied in a wide range of cardiovascular diseases revealing discriminant power for risk stratification (Maestri et al., 2007). Correlation dimension was used to stratify women who suffered hypotension during spinal anesthesia in cesarean section (Chamchad et al., 2004; Bolea et al., 2014a). Sample and approximate entropy were studied comparing control vs. children with congenital heart malformation due to effects of cyanotic and acyanotic defects (Aletti et al., 2012). Furthermore, the integration of linear and nonlinear HRV indices has been shown relevant to stratify cardiac risk patients (Voss et al., 2009) and to describe pathophysiological mechanisms in the cardiovascular and neural system control (Signorini et al., 2011). Some studies pointed out that HRV complexity changes as a result of sympathetic activation (Porta et al., 2007; Turianikova et al., 2011; Weippert et al., 2013).

However, the physiological interpretation of HRV as a marker of ANS activity may be blurred since several factors could affect how intrinsic pacemaker cells and ANS activity are expressed in HRV (Yaniv et al., 2014a). The nonlinear relationship between temporal and complexity HRV indices with respect to HR has been addressed emphasizing the importance of attenuating this effect (Zaza and Lombardi, 2001; Platisa and Gal, 2006; Monfredi et al., 2014; Yaniv et al., 2014b). Furthermore, different mathematical models have demonstrated a relationship between HRV amplitude and HR correcting it (Chiu et al., 2003; Meste et al., 2005; Bailón et al., 2011; Sacha, 2014; Billman et al., 2015). HR correction effect on HRV analysis was studied to predict risk mortality (Pradhapan et al., 2014). Non linear indices, such as correlation dimension, sample, and approximate entropy, are computed over linearly detrended and normalized series so this effect is already compensated for (Osaka et al., 1993; Pincus et al., 1993; Porta et al., 2007; Voss et al., 2009; Bolea et al., 2014a). Despite this normalization, HR may still influence nonlinear HRV indices due to the fact that HR is the intrinsic sampling rate of HRV signal. This implies that the amount of information captured during the same time interval depends on HR. The dependence of nonlinear indices on data length has already been reported (Havstad and Ehlers, 1989), and some studies have computed nonlinear HRV indices over interpolated RR time series at different sampling rates to increase index reliability (Theiler, 1990; Osaka et al., 1993; Hagerman et al., 1996; Radhakrishna et al., 2000; Javorka et al., 2002; Kim et al., 2005). Our hypothesis is that the influence of HR as sampling rate on nonlinear HRV indices is still noticeable even when the same data length is considered.

The goal of this study is to assess and attenuate the HR influence as sampling rate on nonlinear HRV indices in order to provide insight in their physiological interpretation as markers of ANS modulation. To assess the influence of HR on nonlinear HRV indices, a simulation study is conducted in which changes in ANS modulation are independent of changes on mean HR. Based on simulation results, two approaches are proposed to attenuate this mean HR influence. Finally, HR-corrected nonlinear HRV indices are computed over a body position changes database to study their performance under ANS elicitation.

Materials

Body Position Changes (BPC) Database

This database was developed collaboratively at Harvard Medical School, Massachusetts Institute of Technology, and the Favaloro Foundation Medical School. The whole cohort of short-term recordings comes from two data collecting studies. Further details of this database can be found in Sobh et al. (1995).

First Study

Thirteen male subjects of age 21.6 ± 4.4 years (Mean ± SD; range, 19–38 years) with no history of cardiopulmonary disease participated in a study carried out at Clinical Research Center at the Massachusetts Institute of Technology, USA.

Second Study

It comprises groups of subjects of different ages. Only the young group was included in our work (9 subjects, 26.7 ± 4.7 years; range, 20–35 years).

Thus, from the whole database we selected 22 subjects. Two recordings per subject were acquired containing 7-min electrocardiographic (ECG) and respiration (RP) signals, sampled at 360 Hz. The protocol included postural changes. First, ECG and RP signals were recorded while subjects were in supine position. Then, subjects changed to standing position and after 5 min, to allow reaching hemodynamic equilibrium, ECG and RP signals were recorded in standing position. Subjects were asked to breathe following an irregular sequence of tones (Sobh et al., 1995).

Fantasia Database

Twenty young rigorously-screened healthy subjects underwent 120 min of supine resting while continuous ECG and RP signals were recorded at 250 Hz while watching the movie Fantasia, Disney 1940, to help maintain wakefulness. Further database information is available elsewhere (Iyengar et al., 1996) and can be downloaded from http://www.physionet.org (Goldberger et al., 2000).

Methods

ECG Preprocessing

Because the reliability of the HRV analysis can be compromised by low sampling frequency of ECG recordings (Merri et al., 1990), the ECGs belonging to BPC and Fantasia database were interpolated by cubic splines to a frequency of 1080 and 1000 Hz, respectively. Then, heartbeat times, t(k) where k symbolizes the k^th beat, were estimated using an ECG wavelet-based detector (Martínez et al., 2004). Ectopic beats were identified imposing a time-varying threshold on instantaneous heart rate variations. Then, these ectopic beats are corrected using the IPFM model, as described in Mateo et al. (Mateo and Laguna, 2003).

Non-linear HRV Analysis

Approximate, sample entropy and correlation dimension are methods that exploit the phase-space representation of a time series based on Taken's theorem (Takens, 1981). These nonlinear methods are described in the following and further mathematical details are provided in Appendix.

Approximate and Sample Entropy

SampEn and ApEn are irregularity measurements of the time series (Pincus, 1995). Although both entropies are closely related to each other, SampEn was introduced to overcome the self-pairs-related limitation of ApEn computation. Briefly, patterns of time series values (reconstructed vectors) of a certain length (embedded dimension, m) are compared to the rest of the possible pattern candidates. Those comparisons whose differences are below a threshold (r) are summed up and used to calculate correlation sums. The final entropy value measures the changes produced when increasing the length of the patterns in one unit. The parameters m and r have to be previously defined to estimate the entropy values. In this work parameter values are set to m = 2 and r = 0.15 for SampEn. For ApEn, m = 2 and r is set as the threshold that maximizes approximate entropy (ApEn_max; Yentes et al., 2013). ApEn_max was selected instead of ApEn to avoid the bias introduced by in ApEn when considering self-comparisons (Mayer et al., 2016). Its computation is based on a previously published algorithm (Bolea et al., 2014a).

Correlation Dimension

Correlation dimension, D₂, measures the degree of complexity of the system that generates the time series (Grassberger and Procaccia, 1983). In a previous work, we developed techniques to improve the estimation of correlation dimension (Bolea et al., 2014a). On that study log-log curves (logarithm of correlation sums vs. logarithm of thresholds) were fitted to sigmoid curves, thus increasing the accuracy of maximum slope estimation. Moreover, another estimate of correlation dimension denoted as D_2(max) based on the points that maximize the difference between each pair of sigmoid curves was presented. Both D₂ and D_2(max) were computed by varying m = 1–16 and r = 0.01–3 in steps of 0.01.

Non-linear indices estimation may be compromised when the amplitude value of time series appears discrete in a reduced set of values due to the lack of variation. A pre-processing stage is included and details can be found elsewhere (Bolea et al., 2014a).

Simulation Study

A simulation study is conducted to assess the mathematical relationship between HR and nonlinear HRV indices. The simulation study was carried out based on a HRV representation through the IPFM model. This model assumes that the ANS influence on the sinoatrial node can be represented by a modulating signal, (t) (Mateo and Laguna, 2000). According to this model, when the integral of 1+ (t) reaches a threshold, T, a new heartbeat is generated at time instant t(k). Threshold T represents the inverse mean HR.

Fantasia database was selected to compute modulating signals. Assuming that (t) is causal, band-limited and (t) < 1 then the instantaneous HR can be described as:

Instantaneous heart rate d_HR(t) is obtained from the heartbeat times t(k) based on the IPFM model (Mateo and Laguna, 2003), and sampled at 4 Hz. A time-varying mean heart rate d_HRM(t) is computed by low pass filtering d_HR(t) with a cut-off frequency of 0.03 Hz. The heart rate variability signal is obtained as d_HRV(t) = d_HR(t) − d_HRM(t). Finally, the modulating signal, (t), is approximated by $d_{HRV}(t)/\overline{d_{HRM}(t)}$ (Bailón et al., 2011), that is the HRV signal corrected or normalized by the mean HR.

Spectral analysis was applied to 5-min modulating signals (t) by Welch periodogram. Frequency domain indices were estimated based on spectral bands (LF band from 0.04 to 0.15 Hz and HF band from 0.15 to 0.4 Hz). Respiratory frequency was checked to be within the HF band.

Among all modulating signals, only those which presented one marked peak on each band (LF and HF band) were selected. Spectral indices such as the powers and the frequency peaks were used to generate synthetic modulating signals using an autoregressive moving average technique (ARMA; Orini et al., 2012). A total of one hundred 5-min segments were selected and their spectral indices were used to feed the ARMA model. A total of M = 50 stochastic modulating signals ^j(t) with j = 1,…,M, were simulated for each (t). Figure 1 shows the spectra of 50 stochastic realizations, their median spectrum and the one of the segment-recording they are based on.

FIGURE 1

Figure 1. Spectra derived from 50 stochastic realizations of simulated modulating signals during supine conditions (data simulates subject conditions from Fantasia database) applying ARMA technique fixing LF and HF content (red lines). Average spectrum is shown in circles (blue) and spectrum belonging to real data in dashed line (green).

Then, the IPFM model was applied on each stochastic realization, varying the parameter T_n, where n = 1,…,16, corresponding to T from 0.46 to 1.1 s in 0.04 s steps, to simulate the heartbeat occurrence times, $t{}_{Tn}^j\left(k\right)$, from which simulated 300-sample, are obtained. In this way, simulated RR series are generated where ANS modulation is independent from changes in mean HR. Simulation scheme is illustrated in Figure 2. D₂, SampEn and ApEn_(max) are computed over these simulated RR time series.

FIGURE 2

Figure 2. Simulation data scheme is illustrated. Heartbeat time occurrences are detected from ECG. Based on the IPFM model the ANS modulating signal 𝔐(t) is estimated. Spectral analysis by Welch periodogram is computed on (t) in order to estimate the parameters needed for the simulation, frequency and power of LF and HF components, which are used to construct a new set of modulating signals, ^j(t), through ARMA technique, with M = 50 realizations. Then, IPFM model is used to generate simulated heartbeat occurrences $t{}_{Tn}^j\left(k\right)$ with different values of T from 0.46 to 1.1 s. Simulated RR series are computed from the $t{}_{Tn}^j\left(k\right)$.

Another simulation was done based on the BPC database characteristics. However, since subjects were asked to breathe following an irregular sequence of tones, the HF band does not show a dominant peak. In this case, the low and high frequencies used to feed the ARMA model were placed in the middle of LF and HF band, respectively. Then, modulating signals were simulated from spectral indices derived from supine and upright positions. This extends the analysis of HRV dependence on HR under enhanced sympathetic conditions.

Non-linear Indices Dependence on HR as Sampling Rate

The methodology used to compute nonlinear HRV indices, considered in this study, is applied over linearly detrended and normalized RR time series. The detrending ensures that mean HR values are removed from the series whereas the normalization eliminates the influence of mean HR on HRV amplitude. Despite this fact, the effect of mean HR as sampling rate might still be present on them. In this section this effect is investigated on the simulation study, where changes in mean HR are independent from changes in ANS modulation. First, a mathematical relationship between nonlinear HRV indices and HRM is assessed by two regression formulas; then, a HR-correction is proposed based on these formulas. Second, interpolation of RR series is proposed to attenuate the sampling rate influence by mean HR on nonlinear HRV indices.

Regression Formulas

In order to explore the relationship between nonlinear HRV indices and HR the following regression models were proposed.

$\begin{array}{rc}\text{X}=\beta +\alpha RR,(\text{Linear model}),& (2)\end{array}$

$\begin{array}{rc}\text{X}=\beta (R{R}^{\alpha }),(\text{Parabolic model}),& (3)\end{array}$

where X ∈ {D₂, SampEn, ApEn}, α and β are regression coefficients.

Based on the former models HR-correction formulas were obtained by projecting each nonlinear index onto a standard level of RR = 0.5 s, hence:

$\begin{array}{rc}\mathrm{\text{Linear}}:{\text{X}}_{\text{C}1}=\text{X}+\xi (0.5-RR),& (4)\end{array}$

$\begin{array}{rc}\mathrm{\text{Parabolic}}:{\text{X}}_{\text{C}2}=\text{X}{\left(\frac{0.5}{RR}\right)}^{\xi },& (5)\end{array}$

where ξ is the correction factor.

Transformation of X_C1 or X_C2 and RR into linear relationship was used to compute Pearson correlation coefficient ρ. Then, optimization was assessed by total least squares providing correction factors by the Golden Cut Search satisfying ρ(ξ) = 0.

Correction factors were computed on each stochastic realization. Thus, subject-specific correction was defined considering the correction factors of the 50 stochastic realizations for each modulating signal and computing the median of the HR-corrected indices.

Furthermore, a unique correction parameter was computed considering all stochastic realizations for all modulating signals. The transformation and optimization technique described above was applied to median values for each nonlinear index, thus defining a median correction approach to obtain X^ξ_C1 and X^ξ_C2.

Interpolation

RR series are unevenly sampled being the HR its sampling rate. This implies that the number of data information for the same time interval is dependent on HR. On the other hand, it is known that estimation of nonlinear indices such as correlation dimension, sample, and approximate entropy are dependent on data length (Havstad and Ehlers, 1989). Therefore, interpolating RR time series at the same sampling rate may alleviate the influence of mean HR on nonlinear HRV indices since it allows using the same number of data for the same time interval. Interpolation at 2, 4, and 8 Hz were studied (X_I2, X_I4, and X_I8 respectively).

Statistical Analysis

Kolmogorov–Smirnov test was used to test the normality of data distributions. Mann–Whitney U-test was used when necessary; otherwise paired T-test was applied. Furthermore, Pearson correlation was used to assess linear correlation between corrected nonlinear HRV indices and RR. p < 0.05 are considered as statistically significant.

Bland–Altman plots were used to analyse the agreement of subject-specific vs. median correction formulas. The intra-classes coefficient (ICC) was computed by SPSS for Windows, Version 15.0. Chicago, SPSS Inc.

Results

Non-linear HRV Indices and Mean HR Reflect Body Position-Induced Changes

Non-linear HRV indices (D₂, SampEn, and ApEn_max) were computed for the BPC database considering 300-sample segments. All of them were found significantly higher in supine than in standing position (see Figure 3B). Mean HR was also significantly higher in supine than in standing position (Figure 3A), which might explain the statistical differences observed in the computed nonlinear HRV indices.

FIGURE 3

Figure 3. BPC database (22 subjects) was analyzed in which supine and standing positions were compared. (A,B) Illustrate mean RR and uncorrected nonlinear HRV indices while (C) shows HR-corrected nonlinear HRV indices. Asterisks (*) indicates p < 0.05 by Mann–Whitney U-test between supine (Sup) and standing (Std) positions.

Relationship between Non-Linear HRV Indices and RR in the Simulation Study

The relationship between nonlinear HRV indices and RR is assessed in the simulation study where RR is changed without changes in ANS modulation. Non-linear HRV indices computed from simulated data are illustrated in Figure 4 (median values shown as blue circles). The correlation of nonlinear indices and mean HR was evaluated by Pearson correlation coefficient finding high correlation, for a wide range of median index values being 3.5–5.02 for D_2(max), 0.42–1.02 for SampEn, and 0.72–1.24 for ApEn_max (see Table 1).

FIGURE 4

Figure 4. Non-linear HRV indices computed from simulation study varying mean heart period, distribution of all uncorrected nonlinear HRV indices are shown in (blue) circles, corrected by linear regression in (brown) down triangles, corrected by parabolic regression in (green) up triangles, both former cases only subject-specific correction value of non-linear indices for each RR were shown. Finally, HR-corrected non-linear indices by interpolating RR time series at 2, 4, and 8 Hz in (orange) diamonds. Data corresponds to all simulations derived from Fantasia database. Distributions are represented by median and interquartile range.

TABLE 1

Table 1. Pearson correlation factor ρ and p-values of non-linear indices and RR obtained from simulation study.

HR-Corrected Non-linear Indices by Regression Formulas

Regression formulas were applied to each simulated modulating signal (subject-specific approach) providing corrected indices with minimal mean HR correlation. The obtained HR-corrected nonlinear indices are shown in Figure 4 (median values considering all segments, in triangles right and left for linear and parabolic regressions, respectively). The application of correction formulas alleviated the correlation between nonlinear indices and mean HR. Furthermore, the range covered by them was highly reduced (3.81–3.95 for D_2(max), 0.57–0.61 for SampEn, and 0.8–0.9 for ApEn_max), see Table 1.

A set of correction factors (median approach) was obtained by considering the median of all nonlinear index values for each heart rate and computing global correction parameters (Table 1). To evaluate the agreement between subject-specific vs. median correction approaches, the ICC was analyzed, being above 0.8 for all HR-corrected nonlinear indices and for both proposed regression formulas. The Bland–Altman plot in Figure 5 illustrates the difference between both approaches in which D_2(max) is shown as an example.

FIGURE 5

Figure 5. Bland-Altman plots and intra-class coefficient (ICC) illustrate the agreement between subject-specific and median correction approaches computed on Fantasia database to correct HR effect for both proposed regressions for corrected D_2(max) index: linear (left panel) and parabolic (right panel).

HR-Corrected Non-Linear Indices by Interpolation

The nonlinear indices were computed from simulated RR time series resampled at 2, 4, and 8 Hz. As shown in Figure 4, the corrected nonlinear index values obtained applying regression formulas projected onto RR = 0.5 s and by interpolating RR time series at 2 Hz have similar median values. Despite the fact that Pearson correlation factor computed between HR-corrected nonlinear HRV indices by interpolation and mean HR is still significant their range is much reduced [3.84–3.86 for D_2(max), 0.588–0.592 for SampEn, and 0.824–0.836 for ApEn_max] being negligible compared to the range of uncorrected nonlinear ones.

Application to BPC Database

The proposed HR-corrections were evaluated in the BPC database. The results shown in Figure 3C illustrate the differences found between supine and standing conditions. Median and interquartile range of uncorrected and HR-corrected nonlinear HRV indices are provided in Table 2.

TABLE 2

Table 2. Uncorrected Non-linear HRV indices and HR corrected by proposed approaches.

In a first study, the value of the median correction factor ξ extracted from the simulation study was used. It is worth noting that after linear correction there was no significant difference in SampEn and ApEn_max between supine and standing positions, while parabolic correction only reduced differences below significance for SampEn.

In a second study, simulation of each recording's characteristics was computed to apply subject-specific correction, derived independently from supine, and standing recordings. HR-corrected D_2(max) was found statistically significantly different for linear and parabolic regression formulas whereas ApEn_max was only significant for parabolic.

Finally, nonlinear HRV indices were computed on RR time series interpolated at 2, 4, and 8 Hz. We can conclude that the higher the interpolation order, the lower the nonlinear HRV values. In all cases HR-corrected nonlinear indices calculated by interpolation showed statistical differences between positions regardless of the interpolation order used (results of 4 and 8 Hz not included in the manuscript) being their range notably reduced.

Discussion

HRV analysis has been widely used as a non-invasive technique to assess and quantify cardiac ANS modulation (Task Force of the ESC-NASPE, 1996; Voss et al., 2009; Sassi et al., 2015). However, HRV analysis is still under investigation due to HRV characteristics that could lead to physiological misinterpretations (Osaka et al., 1993; Chiu et al., 2003). ANS modulation is linked to ANS tone (HR mean) and, as a consequence, an increase in the sympathetic activity and a decrease in the vagal tone are related to an increment in the HR and a reduction on its variability (Chiu et al., 2003; Kazmi et al., 2016). This implies that there is a physiological correlation between HRV and HR. However, we have demonstrated that there exists also a methodological influence between nonlinear HRV indices and HR due to the fact that HR is the intrinsic sampling rate of HRV. To do this, we have conducted a simulation study. ANS modulating signals were generated as realizations of a stochastic process (Orini et al., 2012). Then, heart beat occurrences are generated using an IPFM model, which is based on action potential generation in SA node cells, and has been proven appropriate to describe the genesis of HRV (Mateo and Laguna, 2000). This simple model allows keeping the ANS modulation constant for different mean HR values, which is not possible in the reality. This simulation allowed assessing the nonlinear dependence of the standard deviation of normal beats on mean HR as it has been pointed out in previous studies (Zaza and Lombardi, 2001; Monfredi et al., 2014; Yaniv et al., 2014b; data not included in the manuscript). Two approaches have been proposed to attenuate the effect shown in the simulation study: regression formulas and interpolation.

Regression Formulas

Regression formulas are commonly used to characterize the relationship between two magnitudes such as ventricular repolarization and heart rate (Pueyo et al., 2004; Smetana et al., 2004; Baumert et al., 2010). The relationship between correlation dimension, sample, and approximate entropy, computed over simulated 300-sample RR series, and mean HR was studied by linear and parabolic regression. Although, the use of other regression models different from linear or parabolic ones may suppose an improvement, it would be unjustified since coefficients of determination R² ≥ 0.9 were obtained for all cases for these models. Then, a correction was proposed based on regression formulas derived for each simulated case, the so-called subject-specific correction, minimizing nonlinear HRV indices correlation to mean HR. A correction based on regression formulas derived from median parameters was proposed as an extension to be applied to other databases. ICC values >0.8 were found when evaluating subject-specific vs. median correction approaches for all nonlinear HRV indices, suggesting the usage of either of the approaches, see Figure 5.

Interpolation

Simulated RR time series were interpolated at 2, 4, and 8 Hz. The higher the interpolation rate, the lower the nonlinear index values. The addition of new data, resulting from interpolation, can be interpreted in terms of entropy as an increase in signal regularity being in concordance with a previous work in which electroencephalogram complexity through correlation dimension was evaluated varying the sampling rate (Jing and Takigawa, 2000). In our study, interpolation was used as a technique to alleviate the dependence of nonlinear HRV indices on mean HR as sampling rate effect since it allows estimating nonlinear indices over the same time interval and the same number of points. Sampling rate should be above the maximum HR. HR-correction nonlinear HRV indices computed by interpolating at 2 Hz and by regression formulas presented similar values and range. In some studies, RR time series were interpolated to increase the number of data points to increase reliability of nonlinear measurements, thus compensating mean HR effect on them. The used sampling frequency varies including 2, 4, and 8 Hz, or even 20 KHz (Osaka et al., 1993; Hagerman et al., 1996; Kim et al., 2005). However, since nonlinear HRV indices estimates are strongly dependent on the selected sampling rate, results should be compared with caution.

Despite the dependence of nonlinear HRV indices on mean HR revealed in the simulation study, no HR-correction of nonlinear HRV indices is applied in most of the studies found in the literature, where mean HR values are even not provided in some cases (Penttilä et al., 2003; Platisa and Gal, 2006; Melillo et al., 2011; Kunz et al., 2012; Moon et al., 2013; Weippert et al., 2013). The application of nonlinear indices without HR correction should be restricted to HR steady-state group conditions, as for example in Weippert et al. (2013).

Classical nonlinear HRV indices evaluated in the BPC database showed around 21, 34, and 21% of reduction in median values from standing with respect to supine position for D_2(max), SampEn, and ApEn_max respectively while mean HR increased by around 28%. Changes in these indices may reflect changes in mean HR as well as additional changes in ANS modulation, as suggested in a previous study (Platisa and Gal, 2006; Yaniv et al., 2014a).

In the BPC database, HR-corrected nonlinear indices were computed under supine and standing conditions and D_2(max) was found to be significantly different for all regression approaches while ApEn_max only for subject-specific using parabolic regression. Linear and parabolic regression formulas were selected to be suitable for the three indices under simulation conditions, although coefficients of determination were slightly lower for ApEn_max and SampEn than for D_2(max).

On the other hand, all nonlinear HRV indices were found still significantly different when corrected by interpolation. It was found a statistically significant reduction in standing with respect to supine position of 18, 30, and 12% for D_2(max), SampEn, and ApEn_max respectively, mostly reflecting ANS modulation changes while mean HR effect was attenuated. HR-corrected nonlinear index ranges, calculated as the difference of median values for supine and standing positions, were found reduced when compared to uncorrected nonlinear HRV index ranges.

Although, regression formulas were studied as HR-correction approach, their suitability depends on simulation requirements. Possible mismatches of simulated data with respect to real data difficult their application and therefore, we propose to compute nonlinear indices over interpolated RR series to attenuate mean HR effect, since no simulation is required, it saves computational time and still differentiates between both positions. This correction may lead to better understanding complexity and regularity under ANS changes unbiased by mean HR as natural sampling rate of RR time series.

Note that, although HR-correction attenuates the effect of mean HR as sampling rate, HR-corrected nonlinear HRV indices may be still correlated with mean HR due to their physiological dependence. After HR-correction, nonlinear HRV indices are capable of capturing information about ANS modulation in response to body position changes.

HR-corrected nonlinear HRV indices addressed in this study, pointed out a reduction in the complexity of the underlying system and an increase in the HRV series regularity caused by an increase of the sympathetic activity, when changing from supine to standing position, being in agreement with previous works with similar conditions, considering tilt table test or even exercise (Osaka et al., 1993; Kamen et al., 2000; Radhakrishna et al., 2000; Javorka et al., 2002; Porta et al., 2007; Bolea et al., 2014b). Nevertheless, these results and their physiological interpretation are limited by the low number of subjects of study and further studies are needed.

Conclusion

In this work, changes in nonlinear HRV indices were studied under different sympathetic conditions where mean HR also changed. It is studied to what extend changes in nonlinear HRV indices are explained by HR ones. Correlation dimension, approximate and sample entropy dependence on mean HR as sampling rate is explored. A simulation study was carried out emulating ANS modulation no linked to mean HR. Simulation results showed that heart rate affects nonlinear indices as it is the intrinsic sampling rate of HRV even when considering the same data length. Two HR-correction methodologies, regression formulas and interpolation, were proposed. Their evaluation on a BPC database revealed a reduction of all studied HR-corrected nonlinear HRV indices in supine and standing positions. After HR-correction, nonlinear HRV indices are capturing changes in the sympathetic modulation by body position-induced changes. HR-correction by interpolation was found suitable to attenuate nonlinear HRV indices effect on mean HR and its application could represent an improvement in their applicability extending it in such cases of non-steady mean HR.

Author Contributions

All authors equally contributed to the conception of the work, revising it critically for important intellectual content, final approval of the version to be published, and agreement to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. EP and RB were supervisors of the research and MO gave methodological support. In addition, JB was responsible for drafting the work.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was funded under projects TEC2013-42140-R, TIN2013-41998-R, and TIN2014-53567-R by MINECO (Spain) and by BSICOS Group (T96) from Government of Aragón, European Social Fund (EU) and by European Research Council (ERC-2014-StG 638284) to (EP). 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 “NANBIOSIS,” more specifically by the High Performance Computing Unit of the CIBER in Bioengineering, Biomaterials, and Nanomedicine (CIBER-BBN) at the University of Zaragoza.

References

Aletti, F., Ferrario, M., Almas de Jesus, T. B., Stirbulov, R., Borghi Silva, A., Cerutti, S., et al. (2012). Heart rate variability in children with cyanotic and acyanotic congenital heart disease: analysis by spectral and non-linear indices. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2012, 4189–4192. doi: 10.1109/embc.2012.6346890

PubMed Abstract | CrossRef Full Text | Google Scholar

Bailón, R., Laouini, G., Grao, C., Orini, M., Laguna, P., and Meste, O. (2011). The integral pulse frequency modulation model with time-varying threshold: application to heart rate variability analysis during exercise stress testing. EEE Trans. Biomed. Eng. 58, 642–652. doi: 10.1109/TBME.2010.2095011

PubMed Abstract | CrossRef Full Text | Google Scholar

Baumert, M., Seeck, A., Faber, R., Nalivaiko, E., and Voss, A. (2010). Longitudinal changes in QT interval variability and rate adaptation in pregnancies with normal and abnormal uterine perfusion. Hypertens. Res. 33, 555–560. doi: 10.1038/hr.2010.30

PubMed Abstract | CrossRef Full Text | Google Scholar

Billman, G. E., Huikuri, H. V., Sacha, J., and Trimmel, K. (2015). An introduction to heart rate variability: methodological considerations and clinical applications. Front. Physiol. 6:55. doi: 10.3389/fphys.2015.00055

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolea, J., Laguna, P., Remartínez, J. M., Rovira, E., Navarro, A., and Bailón, R. (2014a). Methodological framework for estimating the correlation dimension in HRV signals. Comput. Math. Methods Med. 2014:129248. doi: 10.1155/2014/129248

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolea, J., Pueyo, E., Laguna, P., and Bailón, R. (2014b). Non-linear HRV indices under autonomic nervous system blockade. Proc. Ann. Int. Conf. IEEE Eng. Med. Biol. Soc. 2014, 3252–3255. doi: 10.1109/embc.2014.6944316

PubMed Abstract | CrossRef Full Text | Google Scholar

Chamchad, D., Arkoosh, V. A., Horrow, J. C., Buxbaum, J. L., Izrailtyan, I., Nakhamchik, L., et al. (2004). Using heart rate variability to stratify risk of obstetric patients undergoing spinal anesthesia. Anesth. Analg. 99, 1818–1821. doi: 10.1213/01.ane.0000140953.40059.e6

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiu, H. W., Wang, T. H., Huang, L. C., Tso, H. W., and Kao, T. (2003). The influence of mean heart rate on measures of heart rate variability as markers of autonomic function: a model study. Med. Eng. Phys. 25, 475–481. doi: 10.1016/S1350-4533(03)00019-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberger, A. L., Amaral, L. A. N., Glass, L., Hausdorff, J. M., Ivanov, P. C., Mark, R. G., et al. (2000). PhysioBank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals. Circulation 101, e215–e220. doi: 10.1161/01.CIR.101.23.e215

PubMed Abstract | CrossRef Full Text | Google Scholar

Grassberger, P., and Procaccia, I. (1983). Characterization of strange attractors. Phys. Rev. Lett. 50, 346–349. doi: 10.1103/PhysRevLett.50.346

CrossRef Full Text | Google Scholar

Hagerman, I., Berglund, M., Lorin, M., Nowak, J., and Sylvén, C. (1996). Chaos-related deterministic regulation of heart rate variability in time- and frequency domains: effects of autonomic blockade and exercise. Cardiovasc. Res. 31, 410–418. doi: 10.1016/0008-6363(95)00084-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Havstad, J. W., and Ehlers, C. L. (1989). Attractor dimension of nonstationary dynamical systems from small data sets. Phys. Rev. A 39, 845–853. doi: 10.1103/PhysRevA.39.845

PubMed Abstract | CrossRef Full Text | Google Scholar

Iyengar, N., Peng, C. K., Morin, R., Goldberger, A. L., and Lipsitz, L. A. (1996). Age-related alterations in the fractal scaling of cardiac interbeat interval dynamics. Am. J. Physiol. 271, R1078–R1084.

PubMed Abstract | Google Scholar

Javorka, M., Zila, I., Balhárek, T., and Javorka, K. (2002). Heart rate recovery after exercise: relations to heart rate variability and complexity. Braz. J. Med. Biol. Res. 35, 991–1000. doi: 10.1590/S0100-879X2002000800018

PubMed Abstract | CrossRef Full Text | Google Scholar

Jing, H., and Takigawa, M. (2000). Low sampling rate induces high correlation dimension on electroencephalograms from healthy subjects. Psychiatry Clin. Neurosci. 54, 407–412. doi: 10.1046/j.1440-1819.2000.00729.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamen, P. W., Krum, H., and Tonkin, A. M. (2000). The correlation dimension of heart rate variability reflects cardiac autonomic activity. Ann. Noninvasive Electrocardiol. 2, 206–214. doi: 10.1111/j.1542-474X.1997.tb00328.x

CrossRef Full Text | Google Scholar

Karemaker, J. M. (2006). Why measure cardiovascular variability at all? J. Appl. Physiol. 101:684. doi: 10.1152/japplphysiol.00643.2006

CrossRef Full Text

Kazmi, S. Z. H., Zhang, H., Aziz, W., Monfredi, O., Abbas, S. A., Shah, S. A., et al. (2016). Inverse correlation between heart rate variability and heart rate demonstrated by linear and nonlinear analysis. PLoS ONE 11:e0157557. doi: 10.1371/journal.pone.0157557

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, W. S., Yoon, Y. Z., Bae, J. H., and Soh, K. S. (2005). Nonlinear characteristics of heart rate time series: influence of three recumbent positions in patients with mild or severe coronary artery disease. Physiol. Meas. 26, 517–529. doi: 10.1088/0967-3334/26/4/016

PubMed Abstract | CrossRef Full Text | Google Scholar

Kunz, V. C., Borges, E. N., Coelho, R. C., Gubolino, L. A., Martins, L. E., and Silva, E. (2012). Linear and nonlinear analysis of heart rate variability in healthy subjects and after acute myocardial infarction in patients. Braz. J. Med. Biol. Res. 45, 450–458. doi: 10.1590/S0100-879X2012007500025

PubMed Abstract | CrossRef Full Text | Google Scholar

Maestri, R., Pinna, G. D., Accardo, A., Allegrini, P., Balocchi, R., D'Addio, G., et al. (2007). Nonlinear indices of heart rate variability in chronic heart failure patients: Redundancy and comparative clinical value. J. Cardiovasc. Electrophysiol. 18, 425–433. doi: 10.1111/j.1540-8167.2007.00728.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Martínez, J. P., Almeida, R., Olmos, S., Rocha, A. P., and Laguna, P. (2004). A wavelet-based ecg delineator evaluation on standard databases. IEEE Trans. Biomed. Eng. 51, 570–581. doi: 10.1109/TBME.2003.821031

PubMed Abstract | CrossRef Full Text | Google Scholar

Mateo, J., and Laguna, P. (2000). Improved heart rate variability signal analysis from the beat occurrence times according to the IPFM model. IEEE Trans. Biomed. Eng. 47, 985–996. doi: 10.1109/10.855925

PubMed Abstract | CrossRef Full Text | Google Scholar

Mateo, J., and Laguna, P. (2003). Analysis of heart rate variability in the presence of ectopic beats using the heart timing signal. IEEE Trans. Biomed. Eng. 50, 334–343. doi: 10.1109/tbme.2003.808831

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayer, C., Bachler, M., Holzinger, A., Stein, P., and Wassertheurer, S. (2016). The effect of threshold values and weighting factors on the association between entropy measures and mortality after myocardial infarction in the cardiac arrhythmia suppression trial (CAST). Entropy 18:129. doi: 10.3390/e18040129

CrossRef Full Text | Google Scholar

Melillo, P., Bracale, M., and Pecchia, L. (2011). Nonlinear heart rate variability features for real-life stress detection. Case study: students under stress due to university examination. Biomed. Eng. Online 10:96. doi: 10.1186/1475-925X-10-96

PubMed Abstract | CrossRef Full Text | Google Scholar

Merri, M., Farden, D., Mottley, J., and Titlebaum, E. (1990). Sampling frequency of the electrocardiogram for spectral analysis of the heart rate variability. IEEE Trans. Biomed. Eng. 37, 99–106. doi: 10.1109/10.43621

PubMed Abstract | CrossRef Full Text | Google Scholar

Meste, O., Khaddoumi, B., Blain, G., and Bermon, S. (2005). Time-varying analysis methods and models for the respiratory and cardiac system coupling in graded exercise. IEEE Trans. Biomed. Eng. 52, 1921–1930. doi: 10.1109/TBME.2005.856257

PubMed Abstract | CrossRef Full Text | Google Scholar

Monfredi, O., Lyashkov, A. E., Johnsen, A., Inada, S., Schneider, H., Wang, R., et al. (2014). Biophysical characterization of the underappreciated and important relationship between heart rate variability and heart rate. Hypertension 64, 1334–1343. doi: 10.1161/HYPERTENSIONAHA.114.03782

PubMed Abstract | CrossRef Full Text | Google Scholar

Moon, E., Lee, S. H., Kim, D. H., and Hwang, B. (2013). Comparative study of heart rate variability in patients with schizophrenia, bipolar disorder, post-traumatic stress disorder, or major depressive disorder. Clin. Psychopharmacol. Neurosci. 11, 137–143. doi: 10.9758/cpn.2013.11.3.137

PubMed Abstract | CrossRef Full Text | Google Scholar

Orini, M., Bailón, R., Mainardi, L. T., and Laguna, P. (2012). Synthesis of HRV signals characterized by predetermined time-frequency structure by means of time-varying ARMA models. Biomed. Signal Process. Control 7, 141–150. doi: 10.1016/j.bspc.2011.05.003

CrossRef Full Text | Google Scholar

Osaka, M., Saitoh, H., Atarashi, H., and Hayakawa, H. (1993). Correlation dimension of heart rate variability: a new index of human autonomic function. Front. Med. Biol. Eng. 5, 289–300.

PubMed Abstract | Google Scholar

Penttilä, J., Helminen, A., Jartti, T., Kuusela, T., Huikuri, H. V., Tulppo, M. P., et al. (2003). Effect of cardiac vagal outflow on complexity and fractal correlation properties of heart rate dynamics. Auton. Autacoid Pharmacol. 23, 173–179. doi: 10.1046/j.1474-8673.2003.00293.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pincus, S. (1995). Approximate entropy (ApEn) as a complexity measure. Chaos 5, 110–117. doi: 10.1063/1.166092

PubMed Abstract | CrossRef Full Text | Google Scholar

Pincus, S. M., Cummins, T. R., and Haddad, G. G. (1993). Heart rate control in normal and aborted-SIDS infants. Am. J. Physiol. 264, R638–R646.

PubMed Abstract | Google Scholar

Platisa, M. M., and Gal, V. (2006). Reflection of heart rate regulation on linear and nonlinear heart rate variability measures. Physiol. Meas. 27, 145–154. doi: 10.1088/0967-3334/27/2/005

PubMed Abstract | CrossRef Full Text | Google Scholar

Porta, A., Gnecchi-Ruscone, T., Tobaldini, E., Guzzetti, S., Furlan, R., and Montano, N. (2007). Progressive decrease of heart period variability entropy-based complexity during graded head-up tilt. J. Appl. Physiol. 103, 1143–1149. doi: 10.1152/japplphysiol.00293.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Pradhapan, P., Tarvainen, M. P., Nieminen, T., Lehtinen, R., Nikus, K., Lehtimäki, T., et al. (2014). Effect of heart rate correction on pre- and post-exercise heart rate variability to predict risk of mortality – an experimental study on the FINCAVAS cohort. Front. Physiol. 5:208. doi: 10.3389/fphys.2014.00208

PubMed Abstract | CrossRef Full Text | Google Scholar

Pueyo, E., Smetana, P., Caminal, P., de Luna, A. B., Malik, M., and Laguna, P. (2004). Characterization of QT interval adaptation to RR interval changes and its use as a risk-stratifier of arrhythmic mortality in amiodarone-treated survivors of acute myocardial infarction. IEEE Trans. Biomed. Eng. 51, 1511–1520. doi: 10.1109/TBME.2004.828050

PubMed Abstract | CrossRef Full Text | Google Scholar

Radhakrishna, R. K., Dutt, D. N., and Yeragani, V. K. (2000). Nonlinear measures of heart rate time series: influence of posture and controlled breathing. Auton. Neurosci. 83, 148–158.

PubMed Abstract | Google Scholar

Sacha, J. (2014). Interplay between heart rate and its variability: a prognostic game. Front. Physiol. 5:347. doi: 10.3389/fphys.2014.00347

PubMed Abstract | CrossRef Full Text | Google Scholar

Sassi, R., Cerutti, S., Lombardi, F., Malik, M., Huikuri, H. V., Peng, C. K., et al. (2015). Advances in heart rate variability signal analysis: joint position statement by the e-Cardiology ESC Working Group and the European Heart Rhythm Association co-endorsed by the Asia Pacific Heart Rhythm Society. Europace 17, 1341–1353. doi: 10.1093/europace/euv015

PubMed Abstract | CrossRef Full Text | Google Scholar

Signorini, M. G., Ferrario, M., Cerutti, S., and Magenes, G. (2011). Advances in monitoring cardiovascular signals. Contribution of nonlinear signal processing. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2011, 6568–6571. doi: 10.1109/IEMBS.2011.6091620

PubMed Abstract | CrossRef Full Text | Google Scholar

Smetana, P., Pueyo, E., Hnatkova, K., Batchvarov, V., Laguna, P., and Malik, M. (2004). Individual patterns of dynamic QT/RR relationship in survivors of acute myocardial infarction and their relationship to antiarrhythmic efficacy of amiodarone. J. Cardiovasc. Electrophysiol. 15, 1147–1154. doi: 10.1046/j.1540-8167.2004.04076.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobh, J. F., Risk, M., Barbieri, R., and Saul, J. P. (1995). Database for ECG, arterial blood pressure, and respiration signal analysis: feature extraction, spectral estimation, and parameter quantification. Conf. Proc. IEEE Eng. Med. Biol. Soc. 955–956.

Google Scholar

Takens, F. (1981). “Detecting strange attractors in turbulence,” in Dynamical Systems and Turbulence, Warwick 1980, eds D. Rand and L.-S. Young (Berlin, Heidelberg: Springer), 366–381.

Task Force of the ESC-NASPE (1996). Heart rate variability: standards of measurement, physiological interpretation, and clinical use. Circulation 93, 1043–1065. doi: 10.1161/01.CIR.93.5.1043

PubMed Abstract | CrossRef Full Text

Theiler, J. (1990). Estimating fractal dimension. J. Opt. Soc. Am. A 7, 1055–1073. doi: 10.1364/JOSAA.7.001055

CrossRef Full Text | Google Scholar

Turianikova, Z., Javorka, K., Baumert, M., Calkovska, A., and Javorka, M. (2011). The effect of orthostatic stress on multiscale entropy of heart rate and blood pressure. Physiol. Meas. 32, 1425–1437. doi: 10.1088/0967-3334/32/9/006

PubMed Abstract | CrossRef Full Text | Google Scholar

Voss, A., Schulz, S., Schroeder, R., Baumert, M., and Caminal, P. (2009). Methods derived from nonlinear dynamics for analysing heart rate variability. Philos. Trans. A Math. Phys. Eng. Sci. 367, 277–296. doi: 10.1098/rsta.2008.0232

PubMed Abstract | CrossRef Full Text | Google Scholar

Weippert, M., Behrens, K., Rieger, A., Stoll, R., and Kreuzfeld, S. (2013). Heart rate variability and blood pressure during dynamic and static exercise at similar heart rate levels. PLoS ONE 8:e83690. doi: 10.1371/journal.pone.0083690

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaniv, Y., Ahmet, I., Liu, J., Lyashkov, A. E., Guiriba, T.-R., Okamoto, Y., et al. (2014a). Synchronization of sinoatrial node pacemaker cell clocks and its autonomic modulation impart complexity to heart beating intervals. Heart Rhythm 11, 1210–1219. doi: 10.1016/j.hrthm.2014.03.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaniv, Y., Lyashkov, A. E., Sirenko, S., Okamoto, Y., Guiriba, T.-R., Ziman, B. D., et al. (2014b). Stochasticity intrinsic to coupled-clock mechanisms underlies beat-to-beat variability of spontaneous action potential firing in sinoatrial node pacemaker cells. J. Mol. Cell. Cardiol. 77, 1–10. doi: 10.1016/j.yjmcc.2014.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Yentes, J. M., Hunt, N., Schmid, K. K., Kaipust, J. P., McGrath, D., and Stergiou, N. (2013). The appropriate use of approximate entropy and sample entropy with short data sets. Ann. Biomed. Eng. 41, 349–365. doi: 10.1007/s10439-012-0668-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaza, A., and Lombardi, F. (2001). Autonomic indexes based on the analysis of heart rate variability: a view from the sinus node. Cardiovasc. Res. 50, 434–442. doi: 10.1016/S0008-6363(01)00240-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix

Non-Linear Method Calculations

Let x(k), k = 1,…,N be the time series of interest, N being the total number of samples. A set of m-dimensional vectors, y_m(i), called reconstructed vectors, can be generated:

$\begin{array}{ll}{y}_i^m={\left[x(i),x(i+1),x(i+2),\dots ,x(i+(m-1))\right]}^T& (\text{A}1)\end{array}$

Then, the amount of reconstructed vectors is N_m = N − (m−1) for each m-embedded dimension. The distance between each pair of reconstructed vectors, $y{}_i^m$, $y{}_j^m$ is denoted as $d_{i,j}^m$. The norm used is L_∞. Then, each distance is compared with a threshold, r, in order to compute how many reconstructed vectors are within a hyper-space centered in the reconstructed vector of reference. The embedded dimension was set to m = 2 in this manuscript for the three following methods. Let see it mathematically on each particular case.

Approximate Entropy (ApEn)

$\begin{array}{ll}{C}_i^m(r)=\frac{1}{N_m}{\sum }_{j=1}^{N_m}H\left(r-d_{i,j}^m\right)& (\text{A}2)\end{array}$

is the correlation sum where H is the Heaviside function.

$\begin{array}{rc}H(x)=\left\{\begin{array}{c}1x\ge 0\\ 0x<0\end{array}\right\}& (\text{A}3)\end{array}$

$\begin{array}{rc}{\varphi }^m(r)=\frac{1}{N_m}{\sum }_{i=1}^{N_m}log\left(C_i^m(r)\right)& (\text{A}4)\end{array}$

$\begin{array}{rc}ApEn(N,m,r)={\varphi }^m(r)-{\varphi }^{m+1}(r)& (\text{A}5)\end{array}$

In addition, due to the intrinsic characteristics of ApEn, there is a threshold, r, for each time series at which ApEn is maximum defining ApEn_max. Computation of ApEn_max is explained in the last subsection of this Appendix.

Sample Entropy (SampEn)

Self-comparisons, $d_{i,i}^m$, bias ApEn estimation being more notable when data length becomes shorter. Thus, SampEn was introduced to solve this bias and to be data length independent. In this case,

$\begin{array}{ll}{C}_i^m(r)=\frac{1}{N_m-1}{\sum }_{j=1,j\ne i}^{N_m-1}H\left(r-d_{i,j}^m\right)& (\text{A}6)\end{array}$

is the correlation sum not considering self-comparisons.

$\begin{array}{rc}{\varphi }^m(r)=\frac{1}{N_m}{\sum }_{i=1}^{N_m}{C}_i^m(r)& (\text{A}7)\end{array}$

$\begin{array}{ll}SampEn(N,m,r)=-\text{In}({\varphi }^m(r)/{\varphi }^{m+1}(r))& (\text{A}8)\end{array}$

SampEn is also commonly calculated at the same range of r as for ApEn, spanning from 0.1 to 0.2 times standard deviation of the data.

Correlation Dimension (D2)

$\begin{array}{ll}{C}_m(r)=\frac{1}{N_m\left(N_m-1\right)}{\displaystyle {\sum }_{i,j=1}^{N_m}H\left(r-d_{i,j}^m\right)}& (\text{A}9)\end{array}$

For deterministic systems, C_m(r) decreases monotonically to 0 as r approaches 0, and it is expected that C_m(r) is well approximated by $C_m(r)\approx r^{D_2^m}$. For correlation dimension estimation the threshold, r, is varied from r = 0.01–3 in steps of 0.01.

Thus, $D_2^m$ can be defined as:

$\begin{array}{ll}{D}_2^m={\displaystyle \underset{r\to 0}{lim}}\frac{log{C}_m(r)}{log(r)}& (\text{A}10)\end{array}$

For increasing m, $D_2^m$ values tend to saturate to a value D₂ which constitutes the final correlation dimension estimate. Embedded dimension is varied from m = 1–16. In Bolea et al. (2014a), a fast computation algorithm was presented and it is further detailed elsewhere. The algorithm considers self-comparisons in the correlation sums allowing ApEn computation in a middle-step. Therefore, ApEn_max is computed at the same time as correlation dimension