Recurrence Quantitative Analysis of Wavelet-Based Surrogate Data for Nonlinearity Testing in Heart Rate Variability

Exploring the presence of nonlinearity through surrogate data testing provides insights into the nature of physical and biological systems like those obtained from heart rate variability (HRV). Short-term HRV time series are of great clinical interest to study autonomic impairments manifested in chronic diseases such as the end stage renal disease (ESRD) and the response of patients to treatment with hemodialysis (HD). In contrast to Iterative Amplitude Adjusted Fourier Transform (IAAFT), the Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) surrogates preserve nonstationary behavior in time series, a common characteristic of HRV. We aimed to test synthetic data and HRV time series for the existence of nonlinearity. Recurrence Quantitative Analysis (RQA) indices were used as discriminative statistics in IAAFT and PWIAAFT surrogates of linear stationary and nonstationary processes. HRV time series of healthy subjects and 29 ESRD patients before and after HD were tested in this setting during an active standing test. Contrary to PWIAAFT, linear nonstationary time series may be erroneously regarded as nonlinear according to the IAAFT surrogates. Here, a lower proportion of HRV time series was classified as nonlinear with PWIAAFT, compared to IAAFT, confirming that the nonstationarity condition influences the testing of nonlinear behavior in HRV. A contribution of nonlinearity was found in the HRV data of healthy individuals. A lower proportion of nonlinear time series was also found in ESRD patients, but statistical significance was not found. Although this proportion tends to be lower in ESRD patients, as much as 60% of time series proved to be nonlinear in healthy subjects. Given the important contribution of nonlinearity in HRV data, a nonlinear point of view is required to achieve a broader understanding of cardiovascular physiology.

Introduction

Exploring the presence of nonlinearity in data provides insights into the nature of physical and biological systems (Schreiber and Schmitz, 2000). With this aim, surrogate data testing as proposed by Theiler et al. (1992), is applied by the creation of several versions of time series that no longer involve nonlinearity despite preserving statistical properties. A nonlinear statistical measure is then applied to the surrogates and the original time series; any deviation of this measure in the surrogates as compared to the one obtained from such original series is used to discriminate this series from the null hypothesis, which is met by the linear surrogates. Other authors have also proposed to improve the null hypothesis to address more specific types of behavior (Lancaster et al., 2018). Fourier transform-based surrogates generate constrained realizations that virtually preserve the same correlation function of the original data. The statistical null hypothesis of the Iterative Amplitude Adjusted Fourier Transform (IAAFT) technique is that the data represent a stationary linear Gaussian process, measured through an invertible, time-independent instantaneous measurement function (Schreiber and Schmitz, 1996; Lancaster et al., 2018). This technique has been extensively used in physical and biological systems, such as in the analysis of heart rate variability (HRV) (Porta et al., 2015; Faes et al., 2019).

HRV refers to the instantaneous changes in heart rate, measured as the time interval of consecutive R waves in electrocardiography (ECG) recordings, and it is a powerful and simple tool for the study of cardiovascular physiology (No authors listed, 1996). HRV is used in clinical studies and, through decades, it has been considered in various medicine applications (Sassi et al., 2015). Short-term HRV recordings are analyzed for the study of the autonomic nervous system and its influence on the cardiovascular system (No authors listed, 1996; Sassi et al., 2015). The approach of IAAFT surrogates has been used in short-term HRV time series (Porta et al., 2015; Faes et al., 2019), as well as other Fourier transform-based surrogates (Yamamoto et al., 1993; Porta et al., 2000; Lerma et al., 2003; Faes et al., 2004), both offering different results regarding the presence of nonlinear behavior in HRV data.

One of the main disadvantages of Fourier transform-based surrogates is that the original time series must be limited to stationary processes (Borgnat and Flandrin, 2009), while HRV is often nonstationary (Porta et al., 2004). Some other approaches for surrogate data testing consider, as well, the condition of nonstationarity in their null hypothesis (Keylock, 2006; Faes et al., 2009; Lucio et al., 2012). One of them is Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) (Keylock, 2007), which conserves the nonstationary behavior in the surrogate data in a controlled fashion. The analysis of nonlinearity in these settings has led to the increasing application of nonlinear tools in which time series are not required be neither very long nor nonstationary; this is the case of the recurrence plots (RP) (Marwan et al., 2002, 2007). Recurrence quantitative analysis (RQA) is used to quantify diverse nonlinear behaviors in the RP and is widely used in physiological time series, such as electroencephalography (EEG) (Ouyang et al., 2008; Heunis et al., 2018; Pitsik et al., 2020) and ECG (Marwan et al., 2002; Naschitz et al., 2003; Gonzalez et al., 2013; Martín-González et al., 2018).

To study the influence of the autonomic nervous system on cardiovascular dynamics by HRV analysis, the active standing test is used for eliciting controlled parasympathetic predominance at supine position and a sympathetic influence during active standing (Carnethon et al., 2002). Also, it has been of great clinical interest to study autonomic impairments manifested in chronic diseases such as the end stage renal disease (ESRD) (Echeverria et al., 2017; Bokhari et al., 2018). ESRD patients treated with hemodialysis (HD) are in fact subjected to a significant physiological stress during each HD session (Kooman et al., 2018), which involves a sympathetic “challenge” on a regular basis (Lerma et al., 2015) and thus becomes a robust model for the study of autonomic impairments. According to changes indicated by RQA indices during an active standing test, the cardiovascular dynamics associated with both ESRD and HD are consistent with the loss of access to some dynamic physiological conditions (Gonzalez et al., 2013; Calderon-Juarez et al., 2020). Furthermore, recent reports of the correlation between the mean duration of the cardiac cycle (meanNN) with RQA indices in HRV time series (Calderon-Juarez et al., 2020; Robles-Cabrera et al., 2021), suggest that the meanNN parameter, which reflects changes in the cardiac activity required to address different hemodynamic challenges, may influence the nonlinear dynamics of HRV as well. Yet, it has not been fully demonstrated whether the RQA indices in short-term HRV time series exhibit nonlinear dynamics by using surrogate data approach, in particular, considering the nonstationary behavior of these time series.

The purpose of this work was to assess RQA indices as discriminative nonlinear statistics using IAAFT and PWIAAFT surrogates applied to short-term HRV time series of healthy subjects and ESRD patients (before and after treatment with HD) collected during an active standing test.

Materials and Methods

Synthetic Data

Given that, a priori, the type of dynamical behavior of HRV time series was unknown, linear synthetic signals were first tested to ensure that a proper confirmation of the null hypothesis in stationary and nonstationary settings was achieved by the combination of RQA and the algorithms used for the generation of surrogate data; thereby preventing misinterpretations of false detections of nonlinearity in HRV data. The following second order autoregressive (AR2) processes were used as controls for surrogate data testing. These are the same used by Keylock (2007): (a) AR2s, AR2 process with broad energy spectrum, considered as stationary Eq. (1); (b) AR2ns, AR2 process with a peaked energy spectrum Eq. (2). Using these processes, one thousand and eight hundred (1800) values were obtained and the first 1500 values were discarded from both AR2 processes, avoiding transient changes at the beginning of the time series. The remaining 300 values were thus considered to evaluate the linear null hypothesis in synthetic time series, as those of short segments of HRV data, in which near 300 heartbeat intervals are typically contained.

$$x\left(t\right)=0.8x_{t-1}-0.25x_{t-2}+\mathrm{\epsilon }(1)$$

$$x\left(t\right)=1.59x_{t-1}-0.96x_{t-2}+\mathrm{\epsilon }(2)$$

Heart Rate Variability Time Series

Study Protocol

Electrocardiography (ECG) recordings were obtained following the protocol described by Calderon-Juarez et al. (2020). These recordings were obtained during an active orthostatic test from healthy subjects and ESRD patients, as described below. Continuous one-channel ECG recordings were collected during 10 min in supine position followed by subsequent recordings during further 10 min of active standing. The final 5 min of each recording were selected as representative data segments of supine position and active standing, respectively. Patients maintained spontaneous breathing during all procedures. ECG recordings were obtained at 250 samples per second and the identification of R waves was achieved by a second derivative algorithm. The periods of consecutive heart cycles are commonly known as NN or RR intervals, which in turn form the HRV time series. Finally, a correction of artifacts was visually supervised in these series and any outliers by the existence of ectopic beats were replaced using linearly interpolated intervals.

Participants

Forty recordings were obtained in healthy subjects, age 32 years (27–37, CI 95%), body mass index (BMI) 22.06 kg/m² (20–24, CI 95%), proportion of males 34.5%. Twenty-nine ESRD patients were included, age 26 years (24–30, CI 95%, p = 0.084 vs. healthy), BMI 23.3 kg/m² (22–25, CI 95%, p = 0.053 vs. healthy), proportion of males 51.2% (p = 0.295 vs. healthy group). End stage renal disease patients were studied before and after treatment with hemodialysis following the same active orthostatic test protocol. These patients were studied in a previous work (Calderon-Juarez et al., 2020). Hemodialysis sessions had a mean duration of 3.6 ± 0.5 h with total volume removal of 3.1 ± 1.1 L. HD vintage was 12.5 ± 10.2 months with a residual renal function of 0.9 ± 1.5 mL/min. Laboratory results within 1 month prior to the study (obtained from blood samples taken on any day when hemodialysis was not performed) showed creatinine = 8.7 ± 2.5 mg/dL, potassium = 4.9 ± 0.7 mEq/L, phosphorous = 5.1 ± 1.5 mEq/dL, calcium = 8.9 ± 1.1 mg/dL, hemoglobin = 8.3 ± 2.7 g/dL, albumin = 3.9 ± 0.5 g/dL, cholesterol = 165 ± 41 mg/dL, and triglycerides = 145 ± 96 mg/dL. The ESRD etiologies for these patients were systemic lupus erythematosus (n = 1), focal segmental glomerulosclerosis (n = 1), or unknown (n = 27). All procedures following the ethical standards of the 1964 Helsinki’s declaration in its later amendments. Our protocol was approved by the Research and Ethics Committee of the Instituto Nacional de Cardiología Ignacio Chávez (protocol number 21-1236). Informed consent was obtained from all participants.

Hemodialysis Prescription

Hemodialysis (HD) sessions were delivered with volumetric dialysis machines (4008 H, Fresenius Medical Care, Bad Homburg, Germany) using ultrapure dialysate (HCO${}_3^{-}$ = 35 mmol/L, Na⁺ = 138 mmol/L, K⁺ = 2 mmol/L, Ca^{2 +} = 3.5 mEq/L, Mg^{2 +} = 1.0 mEq/L) and polysulfone membranes (F-60 and F-80, Fresenius Medical Care, Walnut Creek, CA, United States). Hypertension was controlled by strict prescription of dry body weight without using antihypertensive drugs following an approach of extracellular volume control by convection. Patients were on a non-restrictive diet and did not use erythropoietin.

Heart Rate Variability Time and Frequency Domain Indices

HRV traditional indices for this study protocol have been reported previously (Gonzalez et al., 2013; Calderon-Juarez et al., 2020). To provide a broad characterization of HRV in the subjects and patients, time domain and frequency domain indices were also calculated here. The meanNN index is the mean value of all RR intervals contained in the time series and SDNN is the standard deviation. Power spectral indices were computed by the Fourier transform method, resampling at 3 Hz and applying a non-overlapped Hamming window of 300 data points with 50% overlap. The Low Frequency band (LF) corresponds to frequencies 0.04–0.15 Hz, and the High Frequency band (HF) corresponds to 0.15–0.4 Hz. High Frequency band is tightly related with parasympathetic activity, whereas LF corresponds to a combination of sympathetic and parasympathetic influence (No authors listed, 1996). We report the LF/HF ratio to express the autonomic modulation as a succinct expression (No authors listed, 1996).

Recurrence Quantitative Analysis

The RQA is based on the construction of recurrence plots, defined by Marwan et al. (2007) :

$$R_{i,j}=\mathrm{\Theta }({\mathrm{\epsilon }}_i-||\overrightarrow{x_i}-\overrightarrow{x_j}||),\overrightarrow{x_i}\in {ℝ}^m,i,j=1,,N,(3)$$

where N is the number of considered states x_i, ε_i is a threshold distance, ||⋅|| a norm and Θ (⋅) is the Heaviside function.

As described thoroughly by Trauth et al. (2019), the representation of multidimensional systems from one-dimensional time series by the time delay embedding approach preserves the dynamic characteristics of the system (Packard et al., 1980). Given that the embedding dimension is sufficiently large, the reconstructed phase space does preserve the topological characteristics of the real phase space (Packard et al., 1980; Takens, 1981). The norm to establish the vicinity for the construction of the recurrence plot must be defined though. Probably, the most common one is the Euclidean norm (the neighborhood is a sphere), in which ε is the radius that contains a fixed number of states (Marwan et al., 2007; Trauth et al., 2019). To study the behavior of nonstationary signals, the fixed amount of nearest neighbors (FAN), in which the radius ε changes for each point, leads to an asymmetric recurrence plot in which all columns have the same recurrence density despite the nonstationary behavior, or trends, in the time series (Marwan, 2011). Therefore, to address this phenomenon, FAN norm has been recommended for analyzing HRV (Martín-González et al., 2018).

The embedding parameters of the time series in this work were calculated with the function of false nearest neighbors (embedding dimension – m) and correlation function (embedding delay – τ). The value of m and τ were selected at the point where the false nearest neighbors and the correlation function reached their first local minimum at zero, respectively (Calderon-Juarez et al., 2020). These parameters were calculated for each time series, and the same set of values were used for RQA of surrogate data. After applying the embedding method for reconstructing the attractor of each HRV time series into the phase space, RPs were obtained with an ε = 0.07, the FAN norm, a Theiler window = τ, window shift = 1 and minimal length of diagonal and vertical lines = 2.

We used the CRP toolbox for MATLAB provided by Marwan et al. (2007), available at (http://tocsy.agnld.uni-potsdam.de/crp.php). The following RQA indices were obtained: recurrence rate (RR), determinism (DET), averaged diagonal length (ADL), length of longest diagonal line (LLDL), entropy of diagonal length (ENT) (Marwan et al., 2007), laminarity (LAM), trapping time (TT) (Marwan and Kurths, 2002), length of longest vertical line (LLVL), recurrence time of the 1st type (T1), recurrence time of the 2nd type (T2) (Gao and Cai, 2000), recurrence period density entropy (RPDE), clustering coefficient (CC) (Marwan et al., 2009) and transitivity (TRANS) (Donner et al., 2010), see Appendix A for definition of RQA indices. For meanNN correlations with RQA indices in surrogate data, the mean values of the RQA indices from the 99 generated surrogates for every subject were obtained.

Stationarity Testing

The existence of restricted weak stationarity (i.e., steady mean and variance) was tested in the synthetic data and original HRV time series to assess the potential implications for surrogate testing of analyzing data with a nonstationary behavior. We followed the algorithm proposed by Porta et al. (2004). A Kolmogorov-Smirnov test goodness-of-fit was used to evaluate if a normal distribution was present in time series, otherwise a logarithmic transformation was applied. N-L+1 ordered sequences of length L were used to create a randomly selected M number of segments or subsets. The length N was set to 300 data points in accordance with the above-mentioned autoregressive processes. L was set to 50 data points to observe at last 5 cycles of LF (about 0.1 Hz); eight M subsets were taken at random to increase the possibility of selecting subsets covering the full extent of time series. After this selection, for the time series with a normal distribution, the stability of the mean and variance was checked by analysis of variance (ANOVA) and Bartlett tests, respectively. For time series with no normal distribution, the stability of the mean and variance was tested using Kruskal-Wallis and Levene tests, respectively. Statistical differences for all tests were considered at the confidence level of p < 0.05.

Surrogate Testing

The IAAFT described by Schreiber and Schmitz (1996) was used for the generation of stationary surrogates with MATLAB toolbox provided in Lancaster et al. (2018). PWIAAFT surrogates were generated with a threshold (ρ) of 0, 0.01, 0.03, and 0.3, which were the same explored in Keylock (2006, 2007). We followed the routine described in detail by Keylock et al. (2011) and used the MATLAB toolbox provided by this author available at (https://sites.google.com/site/chriskeylocknet/software/surrogate-generation-algorithms/pwiaaft). Ninety-nine surrogates were generated from each original time series, being either synthetic data [obtained from Eqs. (1) and (2)] or HRV data (obtained from participants), to achieve a two-sided α error of 0.01. Statistically significant differences of surrogate data testing were considered when the statistic of the original time series was p < 0.05.

Statistical Analysis

Categorical variables are reported as percentages and were compared between healthy subjects and patients by exact Fischer’s tests. For the comparison among the study groups (healthy, ESRD before HD and ESRD after HD), positions (supine and active standing) and surrogate technique (IAAFT and PWIAAFT) a post hoc correction was done by the Bonferroni method. In other words, we compared the proportion of nonlinear time series in IAAFT vs. PWIAAFT (same group and position), supine position vs active standing, healthy vs. ESRD before HD (same position), healthy vs ESRD after HD (same position) and ESRD before HD vs. ESRD after HD (same position). For continuous variables, normal distribution was assessed through Kolmogorov-Smirnov test, median (95% confidence interval) are expressed and were compared with Mann-Whitney U test. Bivariate correlations were tested by the Spearman correlation coefficient. The statistical analyses were performed with the Statistical Package for the Social Sciences (SPSS) version 26, and p-values <0.05 were considered as significant.

Results

Synthetic Data

Stationarity Testing

Figure 1 shows original data from the stationary second order autoregressive process (AR2s – panel A) and nonstationary second order autoregressive (AR2ns – panel D), which were appropriately identified as stationary and nonstationary by the restricted weak stationarity test, respectively (Section “Stationarity testing”). Illustrative examples of the stationarity testing as applied to IAAFT surrogates (middle column) and PWIAAFT (ρ = 0.01) surrogates (right column) are also shown in Figure 1. The IAAFT surrogates of both AR2s (panel B) and AR2ns (panel E) were regarded as stationary. In the PWIAAFT surrogate of AR2s (panel C), a restricted weak stationarity is detected, which was not identified in the PWIAAFT surrogate of AR2ns (panel F). These time series (original and surrogates) correspond to the time series shown in Figures 2, 3.

FIGURE 1

Figure 1. Full-size panels depict synthetic time series, original (left column), one Iterative Amplitude Adjusted Fourier Transform (IAAFT) surrogate time series (middle column) and one Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) (ρ = 0.01) surrogate time series (right column). AR2s Eq. (1) time series corresponds to top row, while the AR2ns one Eq. (2) is in the bottom row. (A) original AR2s, (B) IAAFT surrogate and (C) PWIAAFT (ρ = 0.01) surrogate. (D) AR2ns Eq. (2) and one example of (E) IAAFT surrogate and (F) PWIAAFT (ρ = 0.01) surrogate. Small-size panels show 8 randomly selected segments of 50 data points obtained from the whole time series. The dashed lines represent the means and the dotted lines indicate one standard deviation. AR2s original, IAAFT and PWIAAFT surrogate series were identified as stationary. Whereas AR2ns original and PWIAAFT surrogate series were regarded as nonstationary, the corresponding IAAFT surrogate was identified as stationary. Time series in all panels are shown as arbitrary units.

FIGURE 2

Figure 2. Time series and recurrence plots (m = 5, τ = 3) for synthetic data AR2s Eq. (1). Original time series (panel A), one surrogate obtained by Iterative Amplitude Adjusted Fourier Transform (IAAFT) (panel B) and one surrogate obtained by Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) (panel C). Time series in all panels are shown as arbitrary units.

FIGURE 3

Figure 3. Time series and recurrence plot (m = 5, τ = 3) for synthetic data AR2ns Eq. (2). Original time series (panel A), one surrogate obtained by Iterative Amplitude Adjusted Fourier Transform (IAAFT) (panel B) and one surrogate obtained by Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) (panel C). Time series in all panels are shown as arbitrary units.

Nonlinearity Testing of AR2—Stationary

The AR2s original time series RP is shown in Figure 2 (panel A); the IAAFT algorithm applied to this series generated a noisy pattern in RP (panel B). The following RQA indices obtained from the IAAFT surrogates falsely rejected the null hypothesis: ADL, DET, ENT, LAM, LLDL, LLVL, RR, T2, and TT. On the other hand, the null hypothesis is accepted by considering RPDE, T1, CC, and TRANS. PWIAAFT surrogates with ρ = 0.01 (panel C) assessed with all RQA indices were found consistent with the null hypothesis. We also explored more PWIAAFT surrogates generated by ρ equal to 0, 0.03, and 0.1. In all cases, the same results were obtained.

Nonlinearity Testing of AR2 – Nonstationary

Figure 3 shows the RPs of the AR2ns original data (panel A), IAAFT surrogate (panel B), and PWIAAFT surrogate with ρ = 0.01 (panel C). The results of all RQA indices as applied to IAAFT surrogates were not in accordance with the null hypothesis. These results obtained from the PWIAAFT algorithm were consistent with the null hypothesis for all RQA indices. Other PWIAAFT surrogates generated with the parameter ρ of 0, 0.03 and 0.1 reflected the same findings.

An example of the distribution of a tested RQA statistic (LAM) for ARs and ARns is presented in Figure 4. Regarding IAAFT technique, p = 0.01 for both stationary and nonstationary linear processes, conversely, PWIAAFT surrogates accept the null hypothesis for both linear time series (p > 0.05).

FIGURE 4

Figure 4. Histograms for the laminarity (LAM) values of recurrence plot from of 99 surrogates (orange) obtained with Iterative Amplitude Adjusted Fourier Transform (IAAFT) and Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) techniques. The LAM values measured from the original data are depicted in blue. AR2s, IAAFT (A), PWIAAFT (B); AR2ns, IAAFT (C), PWIAAFT (D).

Heart Rate Variability Data

Time Domain and Spectral Heart Rate Variability Indices

The meanNN index was larger (lower heart rate) in supine position compared with active standing in the healthy group and ESRD patients after HD (Table 1). LF/HF was smaller in supine position compared with active standing in the healthy group. A larger meanNN value was observed in the healthy group compared to ESRD before HD and after HD in supine position. Also, meanNN was larger compared to ESRD after HD during active standing. SDNN was larger in healthy individuals compared to ESRD patients before and after HD in both positions. LF/HF ratio was larger in healthy individuals when compared to ESRD patients before HD, but this difference was not found when compared to ESRD patients after HD. During active standing, LF/HF was different between ESRD patients before and after HD.

TABLE 1

Table 1. Time domain and spectral heart rate variability (HRV) indices shown as median values (95% confidence interval of the median).

Stationarity Testing

Figure 5 shows examples of stationary testing applied to HRV data in supine position (top row) and active standing (bottom row) from a healthy subject (left column), an ESRD patient before HD (middle column) and an ESRD patient after HD (right column). All the examples shown in Figure 5 were classified as nonstationary.

FIGURE 5

Figure 5. Examples of nonstationary heart rate variability (HRV) time series. Full-size panels show the whole HRV time series in supine position (top row) of (A) healthy subject, (B) end stage renal disease (ESRD) patient before hemodialysis (HD), and (C) ESRD patient after HD (same individual). HRV time series collected at active standing (bottom row) from (D) healthy subject, (E) ESRD patient before HD, and (F) ESRD patient after HD. HRV time series units in all panels are shown as seconds (s).

The original HRV time series were mostly classified as nonstationary; only 3 of the 196 analyzed time series were identified as stationary (about 1.5%). The 3 stationary HRV time series were obtained from a healthy subject (supine position), an ESRD patient after HD (supine position), and an ESRD patient after HD (active standing).

Surrogate Data Testing

Examples of HRV time series of healthy and ESRD subjects (before and after HD), RP and corresponding surrogates in supine position and active standing are displayed in Figures 6, 7, respectively (the same examples shown in Figure 5). While the recurrence points are dispersed over all the RP in the IAAFT surrogates (middle column), PWIAAFT surrogates (right column) provide a similar distribution of recurrence points compared to the original time series (left column). This is observed for healthy subjects and ESRD patients before and after HD in both supine position (Figure 6) and active standing (Figure 7). Recurrence quantitative analysis indices in almost all IAAFT surrogates lead to reject the null hypothesis (Table 2). However, in comparison the number of cases with null hypothesis rejections (the percentage of time series in which the surrogate data testing null hypothesis was rejected) was significantly lower using the PWIAAFT surrogates, with exception of CC and TRANS. Although the results of PWIAAFT surrogates shown in Figures 1–7 and Table 2 were generated with ρ = 0.01, we also explored the following values: 0.00, 0.03, and 0.10. We did not find statistically different proportions of rejection rates using these values (Supplementary Table 1).

FIGURE 6

Figure 6. (A–I) Examples of time series and recurrence plots for heart rate variability (HRV) data in supine position. Top row corresponds to a healthy subject (m = 4, τ = 1), (A) original data, (B) Iterative Amplitude Adjusted Fourier Transform (IAAFT) surrogate, and (C) Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) surrogate. Middle row, end stage renal disease (ESRD) patient before hemodialysis (HD) (m = 6, τ = 6), (D) original data, (E) IAAFT surrogate, and (F) PWIAAFT surrogate. Bottom row, ESRD patient after HD (m = 6, τ = 7), (G) original data, (H) IAAFT surrogate, and (I) PWIAAFT surrogate. HRV time series units in all panels are shown as seconds (s).

FIGURE 7

Figure 7. (A–I) Examples of time series and recurrence plots for heart rate variability (HRV) data in active standing. Top row corresponds to a healthy subject (m = 5, τ = 10), (A) original data, (B) Iterative Amplitude Adjusted Fourier Transform (IAAFT) surrogate, and (C) Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) surrogate. Middle row, end stage renal disease (ESRD) patient before hemodialysis (HD) (m = 6, τ = 10), (D) original data, (E) IAAFT surrogate, and (F) PWIAAFT surrogate. Bottom row, ESRD patient after HD (m = 8, τ = 6), (G) original data, (H) IAAFT surrogate, and (I) PWIAAFT surrogate. HRV time series units in all panels are shown as seconds (s).

TABLE 2

Table 2. Percentage (95% confidence Interval) of heart rate variability (HRV) time series in every group that reject the null hypothesis according to the results of different recurrence quantitative analysis (RQA) indices.

Figure 8 shows the percentage of nonlinear time series using LAM statistic. A trend toward lower rejection rates was found in ESRD patients before hemodialysis compared with healthy subjects; this trend can be observed in active standing compared to supine position. However, no statistically significant differences were found. The following rejection rates of LAM correspond to IAAFT surrogates: healthy group supine position 95% (90.72%–99.27%, CI 95%), active standing 100%. End stage renal disease group before HD at supine position and active standing 100% rejections; ESRD group after HD at supine position 96.6% (93.04%–99.9%, CI 95%) and active standing 100%. Rejections rates of LAM with PWIAAFT surrogates were: healthy group at supine position 47.5% (37.71%–57.28%, CI 95%) and active standing 35% (25.65%–44.34%, CI 95%). End stage renal disease group before HD at supine position 20.7% (12.75%–28.64%, CI 95%) and active standing 17.2% (9.8%–24.59%, CI 95%) rejections; ESRD group after HD at supine 37.9% (28.39%–47.4%, CI 95%) and standing 31% (21.93%–40.06%, CI 95%) rejections. The trend toward lower rejection rates in ESRD was also observed in other RQA indices (i.e., DET, ENT, LLVL, TT) (Table 2).

FIGURE 8

Figure 8. Percentage of heart rate variability (HRV) time series of every group that leads to reject the null hypothesis (nonlinearity demonstrated) using Iterative Amplitude Adjusted Fourier Transform (IAAFT) and Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) techniques (bars display the 95% confidence interval). § PWIAAFT vs IAAFT in the same group p < 0.001. There were no significant differences between groups (same position) nor within groups (supine vs active standing, same group).

Correlations With meanNN

The meanNN index is linearly correlated with embedding parameters in healthy subjects, as it is shown in Table 3. This correlation with m is lost in ESRD patients before hemodialysis; however, it is regained after hemodialysis treatment. The meanNN index is also correlated with many of RQA indices in original data (Table 4). These correlations are lost in most parameters (except for LAM, TT, and LLVL) in ESRD patients before hemodialysis. After treatment, meanNN is significantly correlated with all RQA indices.

TABLE 3

Table 3. Spearman correlation coefficient between meanNN and the embedding parameters, tau (τ) and dimension (m).

TABLE 4

Table 4. Spearman correlation coefficients between meanNN and recurrence quantitative analysis (RQA) indices for the original heart rate variability (HRV) time series.

In the same manner as above, we assessed the correlation with meanNN and mean values of RQA in both IAAFT (Table 5) and PWIAAFT (ρ = 0.01) surrogates (Table 6). Regarding IAAFT surrogates, the meanNN is correlated only with RR, T1, and RPDE in the healthy group. In ESRD patients before hemodialysis, meanNN was correlated only with LLDL, LLVL, T2, and RPDE. But after treatment, meanNN was correlated with almost all RQA indices, with the exception of ADL, TT, T1 and T2. Using the PWIAAFT surrogates, we found better preservation in comparison with IAAFT of the correlation between meanNN and RQA indices, as observed in Table 5.

TABLE 5

Table 5. Spearman correlation coefficients between meanNN and recurrence quantitative analysis (RQA) indices for Iterative Amplitude Adjusted Fourier Transform (IAAFT) surrogates of heart rate variability (HRV) time series.

TABLE 6

Table 6. Spearman correlation coefficients between meanNN and recurrence quantitative analysis (RQA) indices for Pinned Wavelet Iterative Amplitude Adjusted Fourier Transform (PWIAAFT) (ρ = 0.01) surrogate data of heart rate variability (HRV) time series.

Discussion

Contribution

We show the application of RQA indices as discriminative nonlinear statics in surrogate data testing and proved the presence of nonlinear structures in short-term HRV time series of healthy subjects and ESRD patients during an active standing test. Other contribution of this work is the implementation of PWIAAFT surrogates for the analysis of HRV data. This method facilitates nonlinear testing as the a priori demonstration of stationarity is not strictly needed. This condition is rarely identified in HRV data (Niccolai et al., 1995; Braun et al., 1998; Porta et al., 2004; Gao et al., 2013), particularly if these data are obtained from healthy subjects studied during daily or ambulatory conditions. Our findings show that even in controlled scenarios, most of healthy subjects and ESRD patients exhibit nonstationary behavior.

Recurrence quantitative analysis (RQA) has been widely used for assessing HRV data, its advantages for the analysis of short, noisy and nonstationary time series becomes a convenient feature for the study of cardiovascular physiology (Marwan et al., 2002). However, nonlinearity by itself, to our best knowledge has not been tested by means of RQA in short-term HRV recordings. Surrogate data testing is a well-known procedure to prove nonlinearity by contradiction. However, the presence of a nonstationary behavior may become a limitation to obtain either reliable HRV indices, such as those provided by the frequency domain analysis (Li et al., 2019), or even appropriate surrogates.

Synthetic Data

In this work we applied the IAAFT technique to linear synthetic stationary and real nonstationary data. It has been suggested that IAAFT surrogates lead to falsely accept the null hypothesis due to their small deviations of the applied statistic measure and rigid preservation of the linear properties in time series (Lancaster et al., 2018). However, some recurrence indices applied here lead to falsely reject the null hypothesis in stationary linear synthetic data. This finding may indicate that RQA is particularly sensible to the randomization of the data and rupture of their structure. In nonstationary synthetic time-series, all statistic measures falsely rejected the null hypothesis, even those that adequately lead to accept the linear hypothesis of stationary data. It is known that “stationarization” (the introduction of stationarity in the timeseries) is a property of surrogates obtained by the IAAFT technique and this may be a reason for higher false rejections (Borgnat and Flandrin, 2009; Lancaster et al., 2018). It is important to emphasize that this technique can lead to null hypothesis rejections because the original time series are either nonlinear or by contrast nonstationary. This phenomenon was previously observed in HRV time series using other discriminative statistics (Faes et al., 2009), finding that the actual rate of rejections decreases once that the technique for surrogate data generation considers nonstationarity. PWIAAFT takes into consideration this characteristic and preserves accurately the original linear structure of the data, as it is shown in this work, for both stationary and non-stationary data. This technique allowed the acceptance of the null hypothesis with all the RQA indices as applied to linear synthetic data.

Heart Rate Variability Data

Traditional HRV indices (Table 1) show the increased sympathetic predominance associated to active standing and ESRD (Gonzalez et al., 2013; Gonzalez-Gomez et al., 2018; Calderon-Juarez et al., 2020). Regarding nonlinear testing of HRV, in a previous study using data generated through Fourier transform-based surrogates (i.e., IAAFT) (Porta et al., 2007), a very low proportion of short-term HRV time series from healthy subjects was found to be nonlinear. But these series are not intuitively expected to be linear due to the nonlinear mechanisms modulating heart rate that are generally considered to be involved. It is possible therefore that such series in that study were too noisy or too short to clearly exhibit nonlinear dynamics. In addition, the activation of the sympathetic branch of the autonomic nervous system decreases the proportion of nonlinear time series, this has been corroborated by pharmacological stimulation and the gradual head-up tilt test (Porta et al., 2007). It has also been suggested that cardiorespiratory coupling confers nonlinear behavior to HRV, because the controlled respiration at a slow rate introduce nonlinear dynamics to HRV (Porta et al., 2000).

In this work, when the IAAFT surrogates were obtained from HRV data, a high rate of null hypothesis was confirmed in relation to RQA indices. Nonetheless, the results for synthetic data demonstrate that these findings can be misleading. Furthermore, only approximately 1.5% of all the HRV time series analyzed in this work were regarded as stationary. As explained above, this is a potential source leading to false nonlinearity detections. It is remarkable that PWIAAFT surrogates show an important decrease in the rate of rejection, similarly to the results shown by time-varying autoregressive surrogate series (Faes et al., 2009), which also involve nonstationary behavior. Added to the well-known PWIAAFT conservation of nonstationarity (Keylock, 2007, 2019; Keylock et al., 2011, 2015) and the ubiquitous presence of nonstationarity in the analyzed HRV time series, the dramatic drop of nonlinearity detection shown by PWIAAFT in comparison to IAAFT is thus likely related to the elimination of the instability of mean and variance in the IAAFT surrogates.

Depending on the RQA index, the percentage of short-term HRV recordings that are found to contain nonlinear properties can be as high as 60% in healthy subjects when DET is used as the statistic measure. For the ESRD patients, the rejection rate decreases to 31% before HD treatment and 34.5% after HD. Furthermore, this rejection rate tends to even lower values in active standing compared with supine position for healthy and ESRD patients, but there were not statistically significant differences regarding this position. These findings suggest that RQA is a suitable tool to detect nonlinearity in short-term series, even when these series manifest nonstationarity. Other pathophysiological conditions, such as acute myocardial infarction have been addressed (Faes et al., 2019) with the surrogate data approach. Patients with this condition tend to show lower proportions of nonlinear HRV times series, which is similar to ESRD patients studied in this work. All these findings suggest that some pathologies suppress nonlinearity from HRV dynamics.

It was proposed by a previous work (Calderon-Juarez et al., 2020) that the meanNN parameter as obtained from HRV data is linearly correlated with some RQA indices in healthy subjects. Notwithstanding that the underlying physiological mechanism of these correlations is not clearly known, an intricate multilayer of physiological interactions could be involved (Kooman et al., 2018). As previously identified (Calderon-Juarez et al., 2020), these correlations are known to be lost in ESRD patients and partially retrieved after hemodialysis. The correlations between meanNN and RQA indices are no longer present in IAAFT surrogates probably owing to the poor conservation of the original time series structure. Yet most of these correlations are preserved with the PWIAAFT surrogates, suggesting that these correlations are partially given by linear statistical and spectral parameters. Some authors have proposed to normalize HRV linear indices by dividing them with the mean heart rate to correct, by this approach, the influence of heart rate on HRV (Hayano et al., 1991). Monfredi et al. (2014) have also shown a robust correlation of mean heart rate and standard deviation of NN intervals; however, they claim that such normalization is insufficient to adequately correct the nonlinear influence of heart rate on HRV (Monfredi et al., 2014). Our work shows that the surrogates HRV time series, in which any nonlinearity structure is destroyed, such correlation of the mean heart rate with RQA is preserved. Notwithstanding that other factors such as age and sex also modify HRV, the meanNN is a determinant characteristic of these time series because it explains a significant dispersion of the RQA indices, thus these indices could also be subjected to normalization by the meanNN.

Limitations and Perspectives

The study of several types of nonlinear behaviors and other types of nonstationarities is beyond the scope of this work. Further research may be conducted to identify which RQA indices are suitable for testing different nonlinear structures. As proposed by Borgnat and Flandrin (2009), nonstationarity can be in fact tested by the generation of stationary surrogate data, which may be considered for future studies of HRV data. Longer HRV time series, which contain enough information to address slower fluctuations and therefore pose different physiological mechanisms of regulation (Lerma et al., 2017), were not explored in this work and these series should be assessed in future projects as well. We collected a small number of ESRD and active standing recordings, thus any potential lower rate of null hypothesis rejections for these data was not possible to be addressed. The respiratory cycle is another physiological factor that influences the HRV time series, its effect remains to be assessed with the combination of techniques presented here. Future studies are required to assess the nonlinear behavior with other HRV indices that are assumed to reflect nonlinearity and to compare them with the present findings.

Conclusion

Recurrence quantitative analysis (RQA) is a suitable framework for the analysis of short, noisy, nonstationary time series and here we also endorse that it is sensitive to capture nonlinear features despite the drawbacks in physiological data analysis that can be introduced by ubiquitous conditions such as the nonstationary behavior. We found that an important proportion of HRV time series from healthy subjects and ESRD patients do contain nonlinear information and hence may be studied from a nonlinear scope point of view to achieve a broader understanding of cardiovascular physiology.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available upon request to the corresponding author, provided pertinent legal requirements are met.

Ethics Statement

The studies involving human participants were reviewed and approved by Research and Ethics Committee of the Instituto Nacional de Cardiología Ignacio Chávez (protocol number 21-1236). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

MC-J, GG, JE, and CL: conceptualization, methodology, and writing—original draft preparation. MC-J and CL: software. GG, HP-G, and CL: resources. MC-J, GG, EQ, and CL: data gathering. MC-J, GG, JE, HP-G, and CL: writing—review and editing. GG: funding acquisition. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Universidad Nacional Autónoma de México (the National Autonomous University of Mexico), grant number DGAPA PAPIIT IN216318. MC-J was supported by the CONACyT and PECEM with a Ph.D. scholarship. The Instituto Nacional de Cardiología Ignacio Chávez and the Universidad Nacional Autónoma de México supported the open access funding.

Conflict of Interest

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2022.807250/full#supplementary-material

References

Bokhari, S. R. A., Inayat, F., Jawa, A., Virk, H. U. H., Awais, M., Hussain, N., et al. (2018). Cardiovascular autonomic neuropathy and its association with cardiovascular and all-cause mortality in patients with end-stage renal disease. Cureus 10:e3243. doi: 10.7759/cureus.3243

PubMed Abstract | CrossRef Full Text | Google Scholar

Borgnat, P., and Flandrin, P. (2009). Stationarization via surrogates. J. Stat. Mechanics 2009:01001.

Google Scholar

Braun, C., Kowallik, P., Freking, A., Hadeler, D., Kniffki, K. D., and Meesmann, M. (1998). Demonstration of nonlinear components in heart rate variability of healthy persons. Am. J. Physiol. 275, H1577–H1584. doi: 10.1152/ajpheart.1998.275.5.H1577

PubMed Abstract | CrossRef Full Text | Google Scholar

Calderon-Juarez, M., Gonzalez-Gomez, G. H., Echeverria, J. C., Perez-Grovas, H., and Lerma, C. (2020). Association between mean heart rate and recurrence quantification analysis of heart rate variability in end-stage renal disease. Entropy (Basel) 22:114. doi: 10.3390/e22010114

PubMed Abstract | CrossRef Full Text | Google Scholar

Carnethon, M. R., Liao, D., Evans, G. W., Cascio, W. E., Chambless, L. E., and Heiss, G. (2002). Correlates of the shift in heart rate variability with an active postural change in a healthy population sample: the atherosclerosis risk in communities study. Am. Heart J. 143, 808–813. doi: 10.1067/mhj.2002.121928

PubMed Abstract | CrossRef Full Text | Google Scholar

Donner, R. V., Zou, Y., Donges, J. F., Marwan, N., and Kurths, J. (2010). Recurrence networks-a novel paradigm for nonlinear time series analysis. New J. Phys. 12:033025.

Google Scholar

Echeverria, J. C., Infante, O., Perez-Grovas, H., Gonzalez, H., Jose, M. V., and Lerma, C. (2017). Effects of orthostatism and hemodialysis on mean heart period and fractal heart rate properties of chronic renal failure patients. Artif. Organs. 41, 1026–1034. doi: 10.1111/aor.12887

PubMed Abstract | CrossRef Full Text | Google Scholar

Faes, L., Gómez-Extremera, M., Pernice, R., Carpena, P., Nollo, G., Porta, A., et al. (2019). Comparison of methods for the assessment of nonlinearity in short-term heart rate variability under different physiopathological states. Chaos 29:123114. doi: 10.1063/1.5115506

PubMed Abstract | CrossRef Full Text | Google Scholar

Faes, L., Pinna, G. D., Porta, A., Maestri, R., and Nollo, G. (2004). Surrogate data analysis for assessing the significance of the coherence function. IEEE Trans. Biomed. Eng. 51, 1156–1166. doi: 10.1109/TBME.2004.827271

PubMed Abstract | CrossRef Full Text | Google Scholar

Faes, L., Zhao, H., Chon, K. H., and Nollo, G. (2009). Time-varying surrogate data to assess nonlinearity in nonstationary time series: application to heart rate variability. IEEE Trans. Biomed. Eng. 56, 685–695. doi: 10.1109/TBME.2008.2009358

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, J., and Cai, H. (2000). On the structures and quantification of recurrence plots. Phys. Lett. A 270, 75–87.

Google Scholar

Gao, J., Gurbaxani, B., Hu, J., Heilman, K., Emauele, V., Lewis, G., et al. (2013). Multiscale analysis of heart rate variability in non-stationary environments. Front. Physiol. 4:119. doi: 10.3389/fphys.2013.00119

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonzalez, H., Infante, O., Perez-Grovas, H., Jose, M. V., and Lerma, C. (2013). Nonlinear dynamics of heart rate variability in response to orthostatism and hemodialysis in chronic renal failure patients: recurrence analysis approach. Med. Eng. Phys. 35, 178–187. doi: 10.1016/j.medengphy.2012.04.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonzalez-Gomez, G. H., Infante, O., Martinez-Garcia, P., and Lerma, C. (2018). Analysis of diagonals in cross recurrence plots between heart rate and systolic blood pressure during supine position and active standing in healthy adults. Chaos 28:085704.

Google Scholar

Hayano, J., Sakakibara, Y., Yamada, A., Yamada, M., Mukai, S., Fujinami, T., et al. (1991). Accuracy of assessment of cardiac vagal tone by heart rate variability in normal subjects. Am. J. Cardiol. 67, 199–204. doi: 10.1016/0002-9149(91)90445-q

PubMed Abstract | CrossRef Full Text | Google Scholar

Heunis, T., Aldrich, C., Peters, J. M., Jeste, S. S., Sahin, M., Scheffer, C., et al. (2018). Recurrence quantification analysis of resting state EEG signals in autism spectrum disorder – a systematic methodological exploration of technical and demographic confounders in the search for biomarkers. BMC Med. 16:101. doi: 10.1186/s12916-018-1086-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Keylock, C. J. (2006). Constrained surrogate time series with preservation of the mean and variance structure. Phys. Rev. E 73:036707. doi: 10.1103/PhysRevE.73.036707

PubMed Abstract | CrossRef Full Text | Google Scholar

Keylock, C. J. (2007). A wavelet-based method for surrogate data generation. Physica D 225, 219–228.

Google Scholar

Keylock, C. J. (2019). Hypothesis testing for nonlinear phenomena in the geosciences using synthetic. Surrogate data. Earth Space Sci. 6, 41–58.

Google Scholar

Keylock, C. J., Stresing, R., and Peinke, J. (2015). Gradual wavelet reconstruction of the velocity increments for turbulent wakes. Phys. Fluids. 27:025104. doi: 10.1063/1.4907740

CrossRef Full Text | Google Scholar

Keylock, C. J., Tokyay, T. E., and Constantinescu, G. (2011). A method for characterising the sensitivity of turbulent flow fields to the structure of inlet turbulence. J. Turbul. 12:N45.

Google Scholar

Kooman, J. P., Katzarski, K., van der Sande, F. M., Leunissen, K. M., and Kotanko, P. (2018). Hemodialysis: a model for extreme physiology in a vulnerable patient population. Semin. Dial. 31, 500–506. doi: 10.1111/sdi.12704

PubMed Abstract | CrossRef Full Text | Google Scholar

Lancaster, G., Iatsenko, D., Pidde, A., Ticcinelli, V., and Stefanovska, A. (2018). Surrogate data for hypothesis testing of physical systems. Phys. Rep. 748, 1–60. doi: 10.1016/j.physrep.2018.06.001

CrossRef Full Text | Google Scholar

Lerma, C., Echeverria, J. C., Infante, O., Perez-Grovas, H., and Gonzalez-Gomez, H. (2017). Sign and magnitude scaling properties of heart rate variability in patients with end-stage renal failure: are these properties useful to identify pathophysiological adaptations? Chaos 27:093906. doi: 10.1063/1.4999470

PubMed Abstract | CrossRef Full Text | Google Scholar

Lerma, C., Gonzalez, H., Perez-Grovas, H., Jose, M. V., and Infante, O. (2015). Preserved autonomic heart rate modulation in chronic renal failure patients in response to hemodialysis and orthostatism. Clin. Exp. Nephrol. 19, 309–318. doi: 10.1007/s10157-014-0990-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Lerma, C., Infante, O., Pérez-Grovas, H., and José, M. V. (2003). Poincaré plot indexes of heart rate variability capture dynamic adaptations after haemodialysis in chronic renal failure patients. Clin. Physiol. Funct. Imaging 23, 72–80. doi: 10.1046/j.1475-097x.2003.00466.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, K., Rüdiger, H., and Ziemssen, T. (2019). Spectral analysis of heart rate variability: time window matters. Front. Neurol. 10:545. doi: 10.3389/fneur.2019.00545

PubMed Abstract | CrossRef Full Text | Google Scholar

Little, M. A., McSharry, P. E., Roberts, S. J., Costello, D. A. E., and Moroz, I. M. (2007). Exploiting nonlinear recurrence and fractal scaling properties for voice disorder detection. BioMed. Eng. 6:23. doi: 10.1186/1475-925X-6-23

PubMed Abstract | CrossRef Full Text | Google Scholar

Lucio, J. H., Valdés, R., and Rodríguez, L. R. (2012). Improvements to surrogate data methods for nonstationary time series. Phys. Rev. E 85:056202. doi: 10.1103/PhysRevE.85.056202

PubMed Abstract | CrossRef Full Text | Google Scholar

Martín-González, S., Navarro-Mesa, J. L., Juliá-Serdá, G., Ramírez-Ávila, G. M., and Ravelo-García, A. G. (2018). Improving the understanding of sleep apnea characterization using recurrence quantification analysis by defining overall acceptable values for the dimensionality of the system, the delay, and the distance threshold. PLoS One 13:e0194462. doi: 10.1371/journal.pone.0194462

PubMed Abstract | CrossRef Full Text | Google Scholar

Marwan, N. (2011). How to avoid potential pitfalls in recurrence plot based data analysis. Int. J. Bifurcat. Chaos 21, 1003–1017. doi: 10.1142/s0218127411029008

CrossRef Full Text | Google Scholar

Marwan, N., and Kurths, J. (2002). Nonlinear analysis of bivariate data with cross recurrence plots. Phys. Lett. A 302, 299–307. doi: 10.1016/s0375-9601(02)01170-2

CrossRef Full Text | Google Scholar

Marwan, N., Carmen Romano, M., Thiel, M., and Kurths, J. (2007). Recurrence plots for the analysis of complex systems. Phys. Rep. 438, 237–329. doi: 10.1016/j.physrep.2006.11.001

CrossRef Full Text | Google Scholar

Marwan, N., Donges, J. F., Zou, Y., Donner, R. V., and Kurths, J. (2009). Complex network approach for recurrence analysis of time series. Phys. Lett. A 373, 4246–4254.

Google Scholar

Marwan, N., Wessel, N., Meyerfeldt, U., Schirdewan, A., and Kurths, J. (2002). Recurrence-plot-based measures of complexity and their application to heart-rate-variability data. Phys. Rev. E 66:026702. doi: 10.1103/PhysRevE.66.026702

PubMed Abstract | CrossRef Full Text | Google Scholar

Monfredi, O., Lyashkov, A. E., Johnsen, A. B., 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

Naschitz, J. E., Itzhak, R., Shaviv, N., Khorshidi, I., Sundick, S., Isseroff, H., et al. (2003). Assessment of cardiovascular reactivity by fractal and recurrence quantification analysis of heart rate and pulse transit time. J. Hum. Hypertens. 17, 111–118. doi: 10.1038/sj.jhh.1001517

PubMed Abstract | CrossRef Full Text | Google Scholar

Niccolai, M., Varanini, M., Macerata, A., Pola, S., Emdin, M., Cipriani, M., et al. (1995). Analysis of non-stationary heart rate series by evolutionary periodogram. Comput. Cardiol. 1995, 10–13.

Google Scholar

No authors listed (1996). Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Circulation 93, 1043–1065.

Google Scholar

Ouyang, G., Li, X., Dang, C., and Richards, D. A. (2008). Using recurrence plot for determinism analysis of EEG recordings in genetic absence epilepsy rats. Clin. Neurophysiol. 119, 1747–1755. doi: 10.1016/j.clinph.2008.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Packard, N. H., Crutchfield, J. P., Farmer, J. D., and Shaw, R. S. (1980). Geometry from a time series. Phys. Rev. Lett. 45, 712–716.

Google Scholar

Pitsik, E., Frolov, N., Hauke Kraemer, K., Grubov, V., Maksimenko, V., Kurths, J., et al. (2020). Motor execution reduces EEG signals complexity: recurrence quantification analysis study. Chaos 30:023111. doi: 10.1063/1.5136246

PubMed Abstract | CrossRef Full Text | Google Scholar

Porta, A., Addio, G. D., Guzzetti, S., Lucini, D., and Pagani, M. (eds) (2004). Testing the presence of non stationarities in short heart rate variability series. Comput. Cardiol. 2004, 19–22.

Google Scholar

Porta, A., Bari, V., Marchi, A., De Maria, B., Cysarz, D., Van Leeuwen, P., et al. (2015). Complexity analyses show two distinct types of nonlinear dynamics in short heart period variability recordings. Front. Physiol. 6:71. doi: 10.3389/fphys.2015.00071

PubMed Abstract | CrossRef Full Text | Google Scholar

Porta, A., Baselli, G., Guzzetti, S., Pagani, M., Malliani, A., and Cerutti, S. (2000). Prediction of short cardiovascular variability signals based on conditional distribution. IEEE Trans. Biomed. Eng. 47, 1555–1564. doi: 10.1109/10.887936

PubMed Abstract | CrossRef Full Text | Google Scholar

Porta, A., Guzzetti, S., Furlan, R., Gnecchi-Ruscone, T., Montano, N., and Malliani, A. (2007). Complexity and nonlinearity in short-term heart period variability: comparison of methods based on local nonlinear prediction. IEEE Trans. Biomed. Eng. 54, 94–106. doi: 10.1109/TBME.2006.883789

PubMed Abstract | CrossRef Full Text | Google Scholar

Robles-Cabrera, A., Torres-Arellano, J. M., Fossion, R., and Lerma, C. (2021). Dependence of heart rate variability indices on the mean heart rate in women with well-controlled type 2 diabetes. J. Clin. Med. 10:4386. doi: 10.3390/jcm10194386

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

Schreiber, T., and Schmitz, A. (1996). Improved surrogate data for nonlinearity tests. Phys. Rev. Lett. 77, 635–638. doi: 10.1103/PhysRevLett.77.635

PubMed Abstract | CrossRef Full Text | Google Scholar

Schreiber, T., and Schmitz, A. (2000). Surrogate time series. Physica D 142, 346–382.

Google Scholar

Takens, F. (1981). “Detecting strange attractors in turbulence,” in Dynamical Systems and Turbulence, Warwick 1980 Lecture Notes in Mathematics, Vol. 898, eds D. Rand and L. S. Young (Berlin: Springer). doi: 10.1007/BF02368233

PubMed Abstract | CrossRef Full Text | Google Scholar

Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., and Doyne Farmer, J. (1992). Testing for nonlinearity in time series: the method of surrogate data. Physica D 58, 77–94. doi: 10.1016/j.jbiomech.2005.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Trauth, M. H., Asrat, A., Duesing, W., Foerster, V., Kraemer, K. H., Marwan, N., et al. (2019). Classifying past climate change in the Chew Bahir basin, southern Ethiopia, using recurrence quantification analysis. Clim. Dynam. 53, 2557–2572. doi: 10.1007/s00382-019-04641-3

CrossRef Full Text | Google Scholar

Yamamoto, Y., Hughson, R. L., Sutton, J. R., Houston, C. S., Cymerman, A., Fallen, E. L., et al. (1993). Operation everest II: an indication of deterministic chaos in human heart rate variability at simulated extreme altitude. Biol. Cybern. 69, 205–212. doi: 10.1007/BF00198960

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix A

Recurrence quantitative analysis indices computed by Cross Recurrence Plot Toolbox for MATLAB.

Recurrence rate (RR) (Marwan et al., 2007)

$$RR=\frac{1}{N^2}\sum _{i,j=1}^N{R}_{i,j}$$

Determinism (DET) (Marwan et al., 2007)

$$DET=\frac{{\sum }_{l=l_{min}}^{N}l{P}^{\in }\left(l\right)}{{\sum }_{i,j}^N{R}_{i,j}}$$

(where P^∈ (l) = {l_i; i = 1…N_i} is the frequency distribution of the lengths l of diagonal structures and N_l is the absolute number of diagonal lines).

Averaged diagonal length (ADL) (Marwan et al., 2007)

$$ADL=\frac{{\sum }_{l_{i=1}}^{l_N}{l}_i}{N}$$

Length of longest diagonal line (LLDL) (Marwan et al., 2007)

$$LLDL=max\left(\{l_i;i=1\mathrm{\dots }{N}_i\}\right)$$

Entropy of diagonal line, Shannon’s entropy (ENT) (Marwan et al., 2007)

$$ENT=-\sum _{l=l_{min}}^{N}p\left(l\right)lnp\left(l\right)$$

Laminarity (LAM) (Marwan et al., 2002)

$$LAM=\frac{{\sum }_{v=v_{min}}^{N}v{P}^{\mathrm{\epsilon }}\left(v\right)}{{\sum }_{v=1}^{N}v{P}^{\mathrm{\epsilon }}\left(v\right)}$$

(where P^ε(v) = {v_i; i = 1…N_v} denotes the frequency distribution of the l lengths of vertical structures).

Trapping time (TT) (Marwan et al., 2002)

$$TT=\frac{{\sum }_{v=v_{min}}^{N}v{P}^{\mathrm{\epsilon }}\left(v\right)}{{\sum }_{v=v_{min}}^N{P}^{\mathrm{\epsilon }}\left(v\right)}$$

Length of longest vertical line (LLVL) (Marwan et al., 2002)

$$LLVL=max\left(\{v_i;i=1\mathrm{\dots }{N}_i\}\right)$$

Recurrence time of the 1^st type (T1) (Gao and Cai, 2000)

$$T_j^1=\left|\left\{i,j:{\stackrel{\rightharpoonup }{x}}_i,{\stackrel{\rightharpoonup }{x}}_j\in R_i\right\}\right|$$

Recurrence time of the 2^nd type (T2) (Gao and Cai, 2000)

$$T_j^2=\left|\left\{i,j:{\stackrel{\rightharpoonup }{x}}_i,{\stackrel{\rightharpoonup }{x}}_j\in R_i;{\stackrel{\rightharpoonup }{x}}_{j-1}\notin R_i\right\}\right|$$

(where R_i are the recurrence points which belong to the state ${\stackrel{\rightharpoonup }{x}}_i$).

Recurrence period density entropy (RPDE) (Little et al., 2007)

$$H_{norm}=\frac{-{\sum }_{i=1}^{T_{max}}P(i)lnP(i)}{ln{T}_{max}}$$

(where P(i) is the recurrence period density, T_max is the maximum recurrence time found in the embedded state space).

Clustering coefficient (CC) (Marwan et al., 2009)

$$CC=\sum _{i=1}^N\frac{{\sum }_{i,j,k=1}^N{R}_{i,j}{R}_{j,k}{R}_{k,i}}{R{R}_i}$$

Transitivity (TT) (Donner et al., 2010)

$$TRANS=\sum _{i=1}^N\frac{{\sum }_{i,j,k=1}^N{R}_{i,j}{R}_{j,k}{R}_{k,i}}{{\sum }_{i,j,k=1}^N{R}_{i,j}{R}_{k,i}}$$