HRV analysis methods

In this section, HRV analysis methods available in Kubios HRV software products are introduced. The presented methods are mainly based on the guidelines given in ^{Task Force 1996}. The presentation of the methods is divided into three main categories, i.e. time-domain, frequency-domain, and nonlinear HRV analysis methods.

Time-domain HRV analysis methods

Mean RR, SDNN, RMSSD and pNN50

The time-domain HRV analysis methods are derived from the beat-to-beat RR interval values in time domain. Let the RR interval time series include $N$ successive beat intervals, i.e. ${\rm RR}=({\rm RR}_1, {\rm RR}_2, \ldots, {\rm RR}_N)$ . The mean RR interval ( $\overline{\rm RR}$ ) and the mean heart rate ( $\overline{\rm HR}$ ) are then defined as

$\begin{equation*} \overline{\rm RR} = \frac{1}{N}\sum_{n=1}^N {\rm RR}_n,\quad \overline{\rm HR} &=& \frac{60}{\overline{\rm RR}} \end{equation*}$

where ${\rm RR}_n$ denotes the value of $n$ ‘th RR interval. Several heart rate variability parameters that measure the variability within the RR time intervals in time-domain exist. The standard deviation of RR intervals (SDNN) is defined as

$\begin{equation*} {\rm SDNN} = \sqrt{\frac{1}{N-1}\sum_{n=1}^N ({\rm RR}_n-\overline{\rm RR})^2} \end{equation*}$

SDNN reflects the overall (both short-term and long-term) variation within the RR interval time series, whereas the standard deviation of successive RR interval differences (SDSD) given by

$\begin{equation*} {\rm SDSD} = \sqrt{E\{\Delta{\rm RR}_n^2\} - E\{\Delta{\rm RR}_n\}^2} \end{equation*}$

is a measure of short-term (beat-by-beat) variability. For stationary RR series $E\{\Delta{\rm RR}_n\} = E\{{\rm RR}_{n+1}\}-E\{{\rm RR}_n\} = 0$ and SDSD equals the root mean square of successive differences (RMSSD) given by

$\begin{equation*} {\rm RMSSD} = \sqrt{\frac{1}{N-1}\sum_{n=1}^{N-1} ({\rm RR}_{n+1} - {\rm RR}_n)^2}. \end{equation*}$

Another measure calculated from successive RR interval differences is the NN50 which is the number of successive intervals differing more than 50 ms or the corresponding relative amount

$\begin{equation*} {\rm pNN50} = \frac{\rm NN50}{N-1} \times 100\% . \end{equation*}$

Geometric measures

In addition to the above statistical measures, there are some geometric HRV analysis methods that are calculated from the RR interval histogram. The HRV triangular index is obtained as the integral of the histogram (i.e. total number of RR intervals) divided by the height of the histogram which depends on the selected bin width. In order to obtain comparable results, a bin width of 1/128 seconds is recommended ^{Task Force 1996}. Another geometric measure is the TINN which is the baseline width of the RR histogram evaluated through triangular interpolation, see Fig. 2A.

The Baevsky’s stress index (SI) is computed according to the formula ^{Baevsky 2009}

$\begin{equation*} {\rm SI} = \frac{{\rm AMo}\times 100\%}{2{\rm Mo\times MxDMn}} \end{equation*}$

where AMo is the so-called mode amplitude presented in percent, Mo is the mode (the most frequent RR interval) and MxDMn is the variation scope reflecting degree of RR interval variability (see Fig. 2B). The mode Mo is simply taken as the median of the RR intervals. The AMo is obtained as the height of the normalised RR interval histogram (bin width 50 msec) and MxDMn as the difference between longest and shortest RR interval values. In order to make SI less sensitive to slow changes in mean heart rate (which would increase the MxDMn and lower AMo), the very low frequency trend is removed from the RR interval time series by using the smoothness priors method ^{Tarvainen et al 2002}. In addition, the square root of SI is taken to transform the tailed distribution of SI values towards normal distribution.

Geometric heart rate variability (HRV) analysis methods

Figure 2: Computation of geometric measures of HRV: A) triangular index (HRVi) and triangular interpolation of RR intervals (TINN), and B) Baevsky’s stress index.

HR deceleration and acceleration capacity

Deceleration capacity (DC) of heart rate is a measure of cardiac parasympathetic modulation as it captures the lengthening of RR interval within 2-4 successive beats as explained below. Acceleration capacity (AC) of heart rate captures the opposite, i.e. shortening of RR interval within few successive beats. These parameters were originally proposed in ^{Bauer et al. 2006} and refined in ^{Nasario-Junior et al. 2014}.

Computation of DC (AC) is performed as follows. First, one RR interval value from every deceleration (acceleration) phase of RR interval time series is identified as an anchor interval. The anchor intervals are those that show steepest change compared to previous RR interval value, i.e. strongest deceleration (acceleration) within every deceleration-acceleration phase (see Fig. 3).

heart rate (HR) deceleration and acceleration

Figure 3: Illustration of heart rate deceleration capacity (DC) and acceleration capacity (AC) computations. The anchor RR intervals (●) and the surrounding RR intervals (○).

Next, all deceleration (acceleration) phases found from the RR interval data are aligned at the anchor points and an ensemble average is computed, i.e. phase rectified signal averaging (PRSA) is applied. Fig. 2 shows the ensemble average (black lines) and the standard deviation intervals (light-blue and light-grey areas) for both deceleration and acceleration phases. DC and AC are calculated from the ensemble averaged deceleration and acceleration phases as follows

$\begin{eqnarray*} {\rm DC} &=& \left[ {\rm RR}(0) + {\rm RR}(1) - {\rm RR}(-1) - {\rm RR}(-2)\right] / 2 \\ {\rm AC} &=& \left[ {\rm RR}(0) + {\rm RR}(1) - {\rm RR}(-1) - {\rm RR}(-2)\right] / 2 \end{eqnarray*}$

That is, DC and AC are computed as a four point difference from the deceleration and acceleration PRSA signals, respectively. To better capture the instantaneous deceleration and acceleration, a modified DC and AC parameters are also computed according to^{Nasario-Junior et al. 2014} as follows

$\begin{eqnarray*} {\rm DC}_{mod} &=& \left[ {\rm RR}(0) - {\rm RR}(-1) \right] \\ {\rm AC}_{mod} &=& \left[ {\rm RR}(0) - {\rm RR}(-1) \right]. \end{eqnarray*}$

Please note that DC and AC sums above are divided by 2 (instead of 4, as originally proposed in ^{Bauer et al. 2006}) in order to reflect the absolute change in RR interval (in msec) similar to the modified indices.

**Table 1:** Summary of time-domain HRV analysis parameters calculated by Kubios HRV software.
PARAMETER	UNITS	DESCRIPTION
EE	[kcal/min]	Energy expenditure (EE) divided into Basal Metabolic Rate (BMR) and activity related energy expenditure. BMR is estimated with Mifflin-St Jeor formula and energy expenditure (EE) using Keytel’s model^{Mifflin et al. 1990, Keytel et al. 2005}
TRIMP	[1/min]	Training impulse (TRIMP) is estimated using the exponential Banister’s model^{Morton et al. 1990}
RR	[ms]	The mean of RR intervals
STD RR (SDNN)	[ms]	Standard deviation of RR intervals
HR	[beats/min]	The mean heart rate
Min & Max HR	[beats/min]	Minimum and maximum HR computed using N beat moving average (default value: N=5)
RMSSD	[ms]	Square root of the mean squared differences between successive RR intervals
NNxx	[beats]	Number of successive RR interval pairs that differ more than xx ms (default value: xx=50)
pNNxx	[%]	NNxx divided by the total number of RR intervals
HRV triangular index	–	The integral of the RR interval histogram divided by the height of the histogram
TINN	[ms]	Baseline width of the RR interval histogram^{Task Force 1996}
Stress index	–	Square root of Baevsky’s stress index^{Baevsky 2009}
DC, AC	[ms]	Deceleration capacity (DC) and acceleration capacity (AC) of heart rate computed as a four-point difference^{Bauer et al. 2006, Nasario-Junior et al. 2014}
DC_mod, AC_mod	[ms]	Modified DC and AC computed as a two-point difference^{Bauer et al. 2006, Nasario-Junior et al. 2014}
SDANN	[ms]	Standard deviation of the averages of RR intervals in 5-min segments^{Task Force 1996}
SDNNI	[ms]	Mean of the standard deviations of RR intervals in 5-min segments^{Task Force 1996}

Frequency-domain HRV analysis methods

In the frequency-domain HRV analysis methods, a power spectrum density (PSD) estimate is calculated for the RR interval series. The regular spectrum estimators implicitly assume equidistant sampling and, thus, the RR interval series is converted to equidistantly sampled series by interpolation methods prior to spectrum estimation. In Kubios HRV software, a cubic spline interpolation method is used. In HRV analysis, the spectrum is generally estimated using either a Fast Fourier transformation (FFT) based methods or parametric autoregressive (AR) modeling based methods. For details on these methods see, e.g., ^{Marple 1987}. The advantage of FFT based methods is the simplicity of implementation, while the AR spectrum yields improved resolution especially for short samples. Another property of AR spectrum that has made it popular in HRV analysis is that it can be factorized into separate spectral components.

In Kubios HRV software, the HRV spectrum is calculated with FFT based Welch’s periodogram method and with the AR method (see Fig. 4). Spectrum factorization in AR method is optional. In the Welch’s periodogram method the HRV sample is divided into overlapping segments. The spectrum is then obtained by averaging the spectra of these segments, which decreases the variance of the FFT spectrum.

Kubios HRV software includes also the Lomb-Scargle periodogram ^{van Dongen et al. 1999}, which differs from the Welch’s periodogram in the sense that it does not assume equidistant sampling and is thus computed directly from the non-interpolated RR interval time series. The variance of the Lomb-Scargle periodogram is decreased by smoothing the periodogram using MA-filering (the window width of the MA-filter can be adjusted in software preferences).

Figure 4: HRV spectrum estimates using FFT based Welch’s periodogram method (left) and autoregressive (AR) modelling based spectrum estimation method.

The generalized frequency bands in case of short-term HRV recordings are the very low frequency (VLF, 0–0.04 Hz), low frequency (LF, 0.04–0.15 Hz), and high frequency (HF, 0.15–0.4 Hz). The frequency-domain measures extracted from a spectrum estimate for each frequency band include absolute and relative powers of VLF, LF and HF bands; LF and HF band powers in normalized units; the LF/HF power ratio; and peak frequencies for each band. In the case of FFT spectrum, absolute power values for each frequency band are obtained by simply integrating the spectrum over the band limits. In the case of AR spectrum, on the other hand, if factorization is enabled distinct spectral components emerge for each frequency band with a proper selection of the model order and the absolute power values are obtained directly as the powers of these components. If factorization is disabled the AR spectrum powers are calculated as for the FFT spectrum. The band powers in relative and normalized units are obtained from the absolute values as described in Table 2.

**Table 2:** Summary of frequency-domain HRV analysis parameters calculated by Kubios HRV software.
PARAMETER	UNITS	DESCRIPTION
Peak frequency	[Hz]	VLF, LF, and HF band peak frequencies
Absolute power	[ms²]	Absolute powers of VLF, LF, and HF bands
Absolute power	[log]	Natural logarithm transformed values of absolute powers of VLF, LF, and HF bands
Relative power	[%]	Relative powers of VLF, LF, and HF bands: VLF [%] = VLF [ms²] / total power [ms²] x 100% LF [%] = LF [ms²] / total power [ms²] x 100% HF [%] = HF [ms²] / total power [ms²] x 100%
Normalized power	[n.u.]	Powers of LF and HF bands in normalised units: LF [n.u.] = LF [ms²] / (total power [ms²] – VLF [ms²]) x 100% HF [n.u.] = HF [ms²] / (total power [ms²] – VLF [ms²]) x 100%
LF/HF	–	Ratio between LF and HF band powers
EDR	[Hz]	ECG derived respiration (available only if ECG data is used for HRV analysis)

Nonlinear HRV analysis methods

Considering the complex control systems of the heart it is reasonable to assume that nonlinear mechanisms are involved in heart rate regulation. Nonlinear HRV analysis methods include Poincaré plot ^{Brennan et al. 2001, Carrasco et al. 2001}, approximate and sample entropy ^{Richman & Moorman 2000, Fusheng et al. 2001}, detrended fluctuation analysis ^{Peng et al. 1995, Penzel et al. 2003}, correlation dimension ^{Guzzetti et al 1996, Henry et al. 2001}, and recurrence plots ^{Webber et al. 1994, Trulla et al. 1996, Zbilut et al. 2002}. During the last years, the number of studies utilizing such methods have increased substantially. The downside of these methods is still, however, the difficulty of physiological interpretation of the results.

Poincaré plot

One commonly used nonlinear HRV analysis method that is simple to interpret is the so-called Poincaré plot. It is a graphical representation of the correlation between successive RR intervals, i.e. plot of ${\rm RR}_{n+1}$ as a function of ${\rm RR}_n$ . The shape of the plot is essential and a common approach to parameterize the shape is to fit an ellipse to the plot as shown in the figure below.

Figure 5: Poincaré plot analysis with the ellipse fitting procedure. SD1 and SD2 are the standard deviations perpendicular to and along the the line-of-identity ${\rm RR}_n = {\rm RR}_{n+1}$ , respectively.

The ellipse is oriented according to the line-of-identity ( ${\rm RR}_n = {\rm RR}_{n+1}$ ). The standard deviation of the points perpendicular to the line-of-identity denoted by SD1 describes short-term variability which is mainly caused by RSA. It can be shown that SD1 is related to the time-domain measure SDSD according to^{Brennan et al. 2001}

$\begin{equation*} \rm SD1^2 = \frac{1}{2} SDSD^2 \end{equation*}$

where SDSD is defined as ${\rm SDSD}=\sqrt{E[\Delta {\rm RR}_n^2] - \overline{\Delta {\rm RR}_n}^2}$ , which is equal to RMSSD for stationary time series. The standard deviation along the line-of-identity denoted by SD2, on the other hand, describes long-term variability and has been shown to be related to time-domain measures SDNN and SDSD by ^{Brennan et al. 2001}

$\begin{equation*} \rm SD2^2 = 2\, SDNN^2 - \frac{1}{2} SDSD^2. \end{equation*}$

The standard Poincaré plot can be considered to be of the first order. The second order plot would be a three dimensional plot of values ( ${\rm RR}_n,\RR_{n+1},{\rm RR}_{n+2}$ ). In addition, the lag can be bigger than 1, e.g., the plot ( ${\rm RR}_n,{\rm RR}_{n+2}$ ).

Approximate entropy

Approximate entropy (ApEn) measures the complexity or irregularity of the signal ^{Richman & Moorman 2000, Fusheng et al. 2001}. Large values of ApEn indicate high irregularity and smaller values of ApEn more regular signal. The ApEn is computed as follows.

First, a set of length $m$ vectors $u_j$ is formed

$\begin{equation*} u_j = ({\rm RR}_j,{\rm RR}_{j+1},\ldots,{\rm RR}_{j+m-1}),\quad j= 1,2,\ldots N-m+1 \end{equation*}$

where $m$ is called the embedding dimension and $N$ is the number of measured RR intervals. The distance between these vectors is defined as the maximum absolute difference between the corresponding elements, i.e.

$\begin{equation*} d(u_j,u_k) = \max\left\{\vert {\rm RR}_{j+l} - {\rm RR}_{k+l}\vert \,\bigl\vert\, l = 0,\ldots,m-1\right\}. \end{equation*}$

Next, for each $u_j$ the relative number of vectors $u_k$ for which $d(u_j,u_k)\le r$ is calculated. This index is denoted with $C_j^m(r)$ and can be written in the form

$\begin{equation*} C_j^m(r) = \frac{{\rm nbr\, of\,}\left\{ u_k\,\bigl\vert\, d(u_j,u_k)\le r \right\}}{N-m+1}\quad\forall\, k. \end{equation*}$

Due to the normalization, the value of $C_j^m(r)$ is always smaller or equal to 1. Note that the value is, however, at least $1/(N-m+1)$ since $u_j$ is also included in the count. Then, take the natural logarithm of each $C_j^m(r)$ and average over $j$ to yield

$\begin{equation*} \Phi^m(r) = \frac{1}{N-m+1}\sum_{j=1}^{N-m+1} \ln{C_j^m(r)}. \end{equation*}$

Finally, the approximate entropy is obtained as

$\begin{equation*} {\rm ApEn}(m,r,N) = \Phi^m(r)-\Phi^{m+1}(r). \end{equation*}$

Thus, the value of the estimate ApEn depends on three parameters, the length $m$ of the vectors $u_j$ , the tolerance value $r$ , and the data length $N$ . In Kubios HRV software the default value of $m$ is set to be $m=2$ . The length $N$ of the data also affects ApEn. When $N$ is increased the ApEn approaches its asymptotic value. The tolerance $r$ has a strong effect on ApEn and it should be selected as a fraction of the standard deviation of the RR interval data (SDNN). This selection enables the comparison of RR data from different subjects. A common selection for $r$ is $r = 0.2 {\rm SDNN}$ , which is also the default value in Kubios HRV software.

Sample entropy

Sample entropy (SampEn) is similar to ApEn, but there are two important differences in its calculation ^{Richman & Moorman 2000, Lake et al. 2002}. For ApEn, in the calculation of the number of vectors $u_k$ for which $d(u_j,u_k)\le r$ also the vector $u_j$ itself is included. This ensures that $C_j^m(r)$ is always larger than 0 and the logarithm can be applied, but at the same time it causes bias to ApEn. In sample entropy the self-comparison of $u_j$ is eliminated by calculating $C_j^m(r)$ as

$\begin{equation*} C_j^m(r) = \frac{{\rm nbr\, of\,}\left\{ u_k\,\bigl\vert\, d(u_j,u_k)\le r \right\}}{N-m}\quad\forall\, k\neq j. \end{equation*}$

Now the value of $C_j^m(r)$ will be between 0 and 1. Next, the values of $C_j^m(r)$ are averaged to yield

$\begin{equation*} C^m(r) = \frac{1}{N-m+1}\sum_{j=1}^{N-m+1} C_j^m(r) \end{equation*}$

and the sample entropy is obtained as

$\begin{equation*} {\rm SampEn}(m,r,N) = \ln{(C^m(r)/C^{m+1}(r))}. \end{equation*}$

The default values set for the embedding dimension $m$ and for the tolerance parameter $r$ in Kubios HRV software are the same as those for the ApEn computation. Both ApEn and SampEn are estimates for the negative natural logarithm of the conditional probability that a data of length $N$ , having repeated itself within a tolerance $r$ for $m$ points, will also repeat itself for $m+1$ points. SampEn was designed to reduce the bias of ApEn and has a closer agreement with the theory for data with known probabilistic content ^{Lake et al. 2002}.

Multiscale entropy (MSE)

Multiscale entropy (MSE) is an extension of SampEn in the sense that it incorporates two procedures ^{Costa et al. 2005}

A coarse-graining process is applied to the RR interval time series. Multiple coarse-grained time series are constructed for the time series by averaging the data points within non-overlapping windows of increasing length $\tau$ , where $\tau$ represents the scale factor and is selected to range between $\tau=1,2,\ldots,20$ . The length of each coarse-grained time series is $N/\tau$ , where $N$ is the number of RR intervals in the data. For scale $\tau=1$ , the coarse-grained time series is simply the original beat-to-beat RR interval time series.
SampEn is calculated for each coarse-grained time series. SampEn as a function of the scale factor produces the MSE. MSE for scale factor $\tau=1$ returns standard SampEn (computed from the original data points).

Detrended fluctuation analysis

Detrended fluctuation analysis (DFA) measures the correlation within the signal. The correlation is extracted for different time scales as follows ^{Peng et al. 1995}. First, the RR interval time series is integrated

$\begin{equation*} y(k) = \sum_{j=1}^k ({\rm RR}_j - \overline{\rm RR}),\quad k = 1,\ldots,N \end{equation*}$

where $\overline{\rm RR}$ is the average RR interval. Next, the integrated series is divided into segments of equal length $n$ . Within each segment, a least squares line is fitted into the data. Let $y_n(k)$ denote these regression lines. Next the integrated series $y(k)$ is detrended by subtracting the local trend within each segment and the root-mean-square fluctuation of this integrated and detrended time series is calculated by

$\begin{equation*} F(n) = \sqrt{\frac{1}{N}\sum_{k=1}^N (y(k)-y_n(k))^2}. \end{equation*}$

This computation is repeated over different segment lengths to yield the index $F(n)$ as a function of segment length $n$ . Typically $F(n)$ increases with segment length. A linear relationship on a double log graph indicates presence of fractal scaling and the fluctuations can be characterized by scaling exponent $\alpha$ (the slope of the regression line relating $\log{F(n)}$ to $\log{n}$ . Different values of $\alpha$ indicate the following

$\alpha = 1.5$ : Brown noise (integral of white noise)
$1 < \alpha < 1.5$ : Different kinds of noise
$\alpha = 1$ : $1/f$ noise
$0.5 < \alpha < 1$ :Large values are likely to be followed by large value and vice versa
$\alpha = 0.5$ : White noise
$0 < \alpha < 0.5$ : Large value is likely to be followed by small value and vice versa

Typically, in DFA the correlations are divided into short-term and long-term fluctuations. In Kubios HRV software, the short-term fluctuations are characterized by the slope $\alpha_1$ obtained from the ( $\log{n}$ , $\log{F(n)}$ ) graph within range $4\le n \le 12$ (default values). Correspondingly, the slope $\alpha_2$ obtained by default from the range $13\le n \le 64$ characterizes long-term fluctuations, see Fig. 6.

detrended fluctuation analysis (DFA) of HRV

Figure 6: Detrended fluctuation analysis. A double log plot of the index $F(n)$ as a function of segment length $n$ . $\alpha_1$ and $\alpha_2$ are the short term and long term fluctuation slopes, respectively.

Correlation dimension

Another method for measuring the complexity or strangeness of the time series is the correlation dimension which was proposed in ^{Grassberger et al. 1983}. The correlation dimension is expected to give information on the minimum number of dynamic variables needed to model the underlying system and it can be obtained as follows.

Similarly as in the calculation of approximate and sample entropies, form length $m$ vectors $u_j$

$\begin{equation*} u_j = ({\rm RR}_j,{\rm RR}_{j+1},\ldots,{\rm RR}_{j+m-1}),\quad j= 1,2,\ldots, N-m+1 \end{equation*}$

and calculate the number of vectors $u_k$ for which $d(u_j,u_k)\le r$ , that is

$\begin{equation*} C_j^m(r) = \frac{{\rm nbr\, of\,}\left\{ u_k\,\bigl\vert\, d(u_j,u_k)\le r \right\}}{N-m+1}\quad\forall\, k \end{equation*}$

where the distance function $d(u_j,u_k)$ is now defined as

$\begin{equation*} d(u_j,u_k) = \sqrt{\sum_{l=1}^m \left( u_j(l)-u_k(l)\right)^2}. \end{equation*}$

Next, an average of the term $C_j^m(r)$ is taken

$\begin{equation*} C^m(r) = \frac{1}{N-m+1}\sum_{j=1}^{N-m+1} C_j^m(r) \end{equation*}$

which is the so-called correlation integral. The correlation dimension $D_2$ is defined as the limit value

$\begin{equation*} D_2(m) = \lim_{r\rightarrow 0}\lim_{N\rightarrow \infty} \frac{\log{C^m(r)}}{\log{r}}. \end{equation*}$

In practice this limit value is approximated by the slope of the regression curve $(\log{r},\log{C^m(r)})$ ^{Henry et al. 2001}. The slope is calculated from the linear part of the log-log plot (see Fig. 7). The slope of the regression curves tend to saturate on the finite value of $D_2$ when $m$ is increased. In the software, a default value of $m=10$ was selected for the embedding.

Figure 7: Approximation of the correlation dimension $D_2$ from the ( $\log r$ , $\log C^m(r)$ ) plot.

Recurrence plot analysis

Yet another approach, included in the software, for analyzing the complexity of the time series is the so-called recurrence plot (RP) analysis^{Marwan et al. 2007}. In this approach, vectors

$\begin{equation*} u_j = ({\rm RR}_j,{\rm RR}_{j+\tau},\ldots,{\rm RR}_{j+(m-1)\tau}),\quad j = 1,2,\ldots,N-(m-1)\tau \end{equation*}$

where $m$ is the embedding dimension and $\tau$ the embedding lag. The vectors $u_j$ then represent the ${\rm RR}$ interval time series as a trajectory in $m$ dimensional space. A recurrence plot is a symmetrical $[N-(m-1)\tau]\times [N-(m-1)\tau]$ matrix of zeros and ones. The element in the $j$ ‘th row and $k$ ‘th column of the RP matrix, i.e. RP( $j$ , $k$ ), is 1 if the point $u_j$ on the trajectory is close to point $u_k$ . That is

$\begin{equation*} {\rm RP}(j,k) = \left\{\begin{array}{ll} 1, & d(u_j-u_k)\le r \\ 0, & {\rm otherwise} \end{array}\right. \end{equation*}$

where $d(u_j,u_k)$ is the Euclidean distance (see above) and $r$ is a fixed threshold. The structure of the RP matrix usually shows short line segments of ones parallel to the main diagonal. The lengths of these diagonal lines describe the duration of which the two points are close to each other. An example RP for heart rate variability time series is presented in figure below.

Figure 8: Recurrence plot matrix for HRV time series (black = 1 and white = 0).

In Kubios HRV software, the embedding dimension and lag were selected to be $m=10$ (default value) and $\tau=1$ (fixed), respectively. The threshold distance $r$ was selected to be $\sqrt{m}\, {\rm SD}$ (default value), where SD is the standard deviation of the RR time series. These selection are similar to those made in^{Dabire et al 1998}. Methods for quantifying recurrence plots were proposed in ^{Webber et al. 1994}. The methods included in Kubios HRV software are introduced below.

The first quantitative measure of RP is the recurrence rate (REC) which is simply the ratio of ones and zeros in the RP matrix. The number of elements in the RP matrix for $\tau=1$ is equal to $N-m+1$ and the recurrence rate is simply given as

$\begin{equation*} {\rm REC} = \frac{1}{(N-m+1)^2}\sum_{j,k = 1}^{N-m+1} {\rm RP}(j,k). \end{equation*}$

The recurrence rate can also be calculated separately for each diagonal parallel to the line-of-identity (main diagonal). The trend of REC as a function of the time distance between these diagonals and the line-of-identity describes the fading of the recurrences for points further away.

The rest of the RP measures consider the lengths of the diagonal lines. A threshold $l_{\rm min} = 2$ is used for excluding the diagonal lines formed by tangential motion of the trajectory. The maximum line length is denoted $l_{\rm max}$ and its inverse, the divergence,

$\begin{equation*} {\rm DIV} = \frac{1}{l_{\rm max}} \end{equation*}$

has been shown to correlate with the largest positive Lyapunov exponent^{Trulla et al. 1996}. The average diagonal line length, on the other hand, is obtained as

$\begin{equation*} l_{\rm mean} = \frac{\sum_{l=l_{\min}}^{l_{\rm max}} l N_{l}}{\sum_{l=l_{\rm min}}^{l_{\rm max}} N_{l}} \end{equation*}$

where $N_{l}$ is the number of length $l$ lines. The determinism of the time series is measured by the variable

$\begin{equation*} {\rm DET} = \frac{\sum_{l=l_{\rm min}}^{l_{\rm max}} lN_{l}}{\sum_{j,k=1}^{N-m+1} {\rm RP}(j,k)}. \end{equation*}$

Finally, the Shannon information entropy of the line length distribution is defined as

$\begin{equation*} {\rm ShanEn} = -\sum_{l=l_{\rm min}}^{l_{\rm max}} n_{l} \ln n_{l} \end{equation*}$

where $n_{l}$ is the number of length $l$ lines divided by the total number of lines, that is

$\begin{equation*} n_{l} = \frac{N_{l}}{\sum_{l'=l_{\rm min}}^{l_{\rm max}} N_{l'}}. \end{equation*}$

**Table 3:** Summary of nonlinear HRV analysis parameters calculated by Kubios HRV software.
PARAMETER	UNITS	DESCRIPTION
SD1	[ms]	In Poincaré plot, the standard deviation perpendicular to the line-of-identity ^{Brennan et al. 2001, Carrasco et al. 2001}
SD2	[ms]	In Poincaré plot, the standard deviation along the line-of-identity
SD2/SD1	–	Ratio between SD2 and SD1
ApEn	–	Approximate entropy^{Richman & Moorman 2000, Fusheng et al. 2001}
SampEn	–	Sample entropy^{Richman & Moorman 2000}
DFA, α₁	–	In detrended fluctuation analysis, short term fluctuation slope^{Peng et al. 1995, Penzel et al. 2003}
DFA, α₂	–	In detrended fluctuation analysis, long term fluctuation slope
D₂	–	Correlation dimension^{Guzzetti et al. 1996, Henry et al. 2001}
RPA		Recurrence plot analysis^{Webber et al. 1994, Dabire et al. 1998, Zbilut et al. 2002, Marwan et al. 2007}
Lmean	[beats]	Mean line length
Lmax	[beats]	Maximum line length
REC	[%]	Recurrence rate
DET	[%]	Determinism
ShanEn	–	Shannon entropy
MSE	–	Multiscale entropy for scale factor values τ=1,2,…,20 ^{Costa et al. 2005}

References

R. M. Baevsky. Methodical recommendations use kardivar system for determination of the stress level and estimation of the body adaptability standards of measurements and physiological interpretation. 2009.
A. Bauer, J.W. Kantelhardt, P. Barthel, R. Schneider, T. Mäkikallio, K. Ulm, K. Hnatkova, A. Schöming, H. Huikuri, A. Bundle, M. Malik, and G. Schmidt. Deceleration capacity of heart rate as a predictor of mortality after myocardial infarction: cohort study. Lancet, 367:1647–81, 2006.
M. Brennan, M. Palaniswami, and P. Kamen. Do existing measures of Poincaré plot geometry reflect nonlinear features of heart rate variability. IEEE Trans Biomed Eng, 48(11):1342–1347, 2001.
S. Carrasco, M.J. Caitán, R. González, and O. Yánez. Correlation among Poincaré plot indexes and time and frequency domain measures of heart rate variability. J Med Eng Technol, 25(6):240–248, November/December 2001.
M. Costa, A.L. Goldberger, and C.-K. Peng. Multiscale entropy analysis of biological signals. Physical Rev E, 71:021906, 2005.
H. Dabire, D. Mestivier, J. Jarnet, M.E. Safar, and N. Phong Chau. Quantification of sympathetic and parasympathetic tones by nonlinear indexes in normotensive rats. amj, 44:H1290–H1297, 1998.
Y. Fusheng, H. Bo, and T. Qingyu. Approximate entropy and its application in biosignal analysis. In M. Akay, editor, Nonlinear Biomedical Signal Processing: Dynamic Analysis and Modeling, volume II, chapter 3, pages 72–91. IEEE Press, New York, 2001.
P. Grassberger and I. Procaccia. Characterization of strange attractors. Phys Rev Lett, 50:346–349, 1983.
S. Guzzetti, M.G. Signorini, C. Cogliati, S. Mezzetti, A. Porta, S. Cerutti, and A. Malliani. Non-linear dynamics and chaotic indices in heart rate variability of normal subjects and heart-transplanted patients. Cardiovascular Research, 31:441–446, 1996.
B. Henry, N. Lovell, and F. Camacho. Nonlinear dynamics time series analysis. In M. Akay, editor, Nonlinear Biomedical Signal Processing: Dynamic Analysis and Modeling, volume II, chapter 1, pages 1–39. IEEE Press, New York, 2001.
L.R. Keytel, J.H. Goedecke, T.D. Noakes, H. Hiilloskorpi, R. Laukkanen, L. Van Der Merwe, and E.V. Lambert. Prediction of energy expenditure from heart rate monitoring during submaximal exercise. J Sports Sci, 23(3):289–297, 2005.
D.E. Lake, J.S. Richman, M.P. Griffin, and J.R. Moorman. Sample entropy analysis of neonatal heart rate variability. ajp, 283:R789–R797, September 2002.
S.L. Marple. Digital Spectral Analysis. Prentice-Hall International, 1987.
N. Marwan, M.C. Romano, M. Thiel, and J. Kurths. Recurrence plots for the analysis of complex systems. Physics Reports, 438:237-329, 2007.
M.D. Mifflin, S.T. St Jeor, L.A. Hill, B.J. Scott, S.A. Daugherty, and Y.O. Koh. A new predictive equation for resting energy expenditure in healthy individuals. Am J Clin Nutr, 51:241-7, 1990.
R.H. Morton, J.R. Fitz-Clarke and E.W. Banister. Modeling human performance in running. J Appl Physiol, 69(3):1171-7, 1990.
O. Nasario-Junior, P.R. Benchimol-Barbosa, and J. Nadal. Refining the deceleration capacity index in phase-rectified signal averaging to assess physical conditioning level. J Electrocardiol, 47(3):306–310, 2014.
C.-K. Peng, S. Havlin, H.E. Stanley, and A.L. Goldberger. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos, 5:82–87, 1995.
T. Penzel, J.W. Kantelhardt, L. Grote, J.-H. Peter, and A. Bunde. Comparison of detrended fluctuation analysis and spectral analysis for heart rate variability in sleep and sleep apnea. IEEE Trans Biomed Eng, 50(10):1143–1151, October 2003.
J.A. Richman and J.R. Moorman. Physiological time-series analysis using approximate entropy and sample entropy. Am J Physiol, 278:H2039–H2049, 2000.
M.P. Tarvainen, P.O. Ranta-aho, and P.A. Karjalainen. An advanced detrending method with application to HRV analysis. IEEE Trans Biomed Eng, 49(2):172–175, February 2002.
Task force of the European society of cardiology and the North American society of pacing and electrophysiology. Heart rate variability – standards of measurement, physiological interpretation, and clinical use. Circulation, 93(5):1043–1065, March 1996.
L.L. Trulla, A. Giuliani, J.P. Zbilut, and C.L. Webber Jr. Recurrence quantification analysis of the logistic equation with transients. Phys Lett A, 223(4):255–260, 1996.
H.P.A. Van Dongen, E. Olofsen, J.H. VanHartevelt, and E.W. Kruyt. Searching for biological rhythms: peak detection in the periodogram of unequally spaced data. J Bioloogical Rhythms, 14(6):617–620, 1999.
C.L. Webber Jr. and J.P. Zbilut. Dynamical assessment of physiological systems and states using recurrence plot strategies. J Appl Physiol, 76:965–973, 1994.
J.P. Zbilut, N. Thomasson, and C.L. Webber. Recurrence quantification analysis as a tool for the nonlinear exploration of nonstationary cardiac signals. Med Eng Phys, 24:53–60, 2002.