2.6. Recurrence Plots and Recurrence Quantification Analysis

Recurrence plots (RP) allow visualising system states in the phase space. In complex dynamical systems, recurrence is related to the temporal evolution of dynamical systems trajectories in the phase space.

Generally, to compute an RP, an embedding dimension (m) and a time-delay (τ) are necessary. Delay, τ, is the minimum time lag to minimise the autocorrelation of a time series. Then, m represents the number of independent variables needed to characterise the dynamics of the system. Finally, RP is a square matrix, with time on both axes, of pairwise Euclidean distances between the reconstructed system states to which a distance threshold (ε) is applied [43]. Mathematically RP is defined as:

where N is the number of measured states x→i, Θ is the Heaviside step function (i.e., Θ(x) = 1, if ‖x→i−x→j‖≤ε, and Θ(x) = 0 otherwise), ‖·‖ is a norm, and ε is a threshold previously defined based on the time-series properties. In this study, the phase space trajectories are based on the Euclidean distance between x→i and x→j of the series. If Rij=1 at a time (i, j), is marked as a black dot in the position (i, j). Otherwise, if Rij=0, the corresponding states will be represented as white dots.

The same principle is maintained in the cross-recurrence plot (CRP) methodology. However, in CRP, two different time series are analysed simultaneously, and black dots represent the co-occurrence of similar states between two time series. Mathematically, a CRP (x1, x2,⋯xi,⋯xN) and (y1, y2,⋯yj,⋯yN), is calculated by:

Several measures of complexity have been proposed to be quantified by the RQA, though, in this work, we focused on Determinism (DET), Average length of structures (LT), Shannon’s Entropy (ENT), Laminarity (LAM) and trapping time (TT), the extended formulas are added in the Appendix A. Furthermore, RQA was extended by computing the diagonal-wise recurrence quantification profile [55]. The recurrence rate around the line of coincidence (LOC) and the surroundings time lags was calculated to measure the two time-series coupling as a lag function. The maximum number of lags to be analysed was six, the same as the cross-correlation method.

CRQA R package [55,56] was used to construct RP, obtain RQA measures and compute diagonal-wise recurrence profile. First, the VIs series were normalised using a z-score normalisation; then, the distance matrix was rescaled based on the maximum value following the recommendations of Webber and Zbilut [57]. Optimizeparam function is then computed to find the three parameters’ optimal values (τ, m, and ε). The delay (τ) was found by obtaining the local minimum where mutual information drops to both series [58]. The embedding dimension (m) was calculated by the false nearest neighbours’ algorithm [59]. The threshold ε was estimated by an iterative process based on the time-series’ standard deviation (SD). In this work, ε was limited to 5% of the recurrence rate (RR) in all the cases. When multiples values of m, τ, were obtained by the optimisation, the maximum of them was selected as the optimal value to apply in the construction of CRPs.

The quantification of RP and CRPs structures was calculated with the Crqa function using the values obtained from the optimisation function. Then, the drpfromts function was computed to plot diagonal-wise recurrence profiles in the RPs and CRPs.