HRV is a well recognized method to explore autonomic nervous system activity for laboratory and clinical uses (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996; Rajendra Acharya et al., 2006; Billman, 2011; Billman et al., 2015). Indeed, except for nonsinus rhythms induced by atrial fibrillation or extrasystoles, the main variations in HR are due to the permanent beat-to-beat adaptations of parasympathetic and sympathetic activities on the heart (Hon and Lee, 1963; Sayers, 1973; Billman, 2011). The quantification of these variations then allows to evaluate autonomic activity (Akselrod et al., 1981). This is achieved by various methods: time and frequency domain, nonlinear, and time-frequency (Rajendra Acharya et al., 2006; Pichot et al., 2016). Each calculated index provides different information, but some are redundant (Kautzner and Hnatkova, 1995). Globally, the rapid variations concerning several beats are linked to parasympathetic activity (Akselrod et al., 1981; Pomeranz et al., 1985); slower variations with a period 7–25 s are due to both sympathetic and parasympathetic activities (Saul et al., 1990; Pagani et al., 1997); and very slow variations with a period 25–300 s are usually assimilated to parasympathetic activity, renin angiotensine system and thermoregulation (Taylor et al., 1998). Ultraslow variations have been shown but required 24-h recordings to be calculated (Bigger et al., 1992). Then, for example, from the time and frequency domain methods, both SDNN and Ptot give information on the total beat-to-beat variations, i.e., an approximation of the global regulatory capacities of the autonomic nervous system; rMSSD, pNN50 and HF are linked to short-term beat-to-beat variations and are related to parasympathetic activity; the LF/HF ratio can be considered as an approximation of the sympatho-vagal balance (Pagani et al., 1984), although this index should be considered with caution (Billman, 2013).

Precautions must also be taken with regard to respiratory rate. Indeed, the physiological interpretations of the frequency domain indices described above are only valid if the subject’s respiratory frequency is in the high-frequency bandwidth, i.e., between 9 and 24 cycles per minute for HF bandwidth set between 0.15 and 0.40 Hz. For slower breathing rates, the peak of the LF will be contaminated by the effects of breathing and the HF will be artificially reduced. For respiration that is too fast, its effects will not be taken into account in the HF. It is also important to take into account the state in which the subject is: in the upright position or during hand-grip exercise, for example, the LF will be more influenced by sympathetic activity than in the supine position at rest, where parasympathetic activity will predominate.

A method that is often used to complement time domain analysis, the Poincaré plot is a graphical representation of the RR(n) as a function of RR(n + 1), modelled in the form of an elipse (Kamen et al., 1996). Then, the standard deviation of the RRs along the x = y-axis (index SD1) represents short-term variations in RR variability, while the standard deviation on the perpendicular axis (index SD2) represents longer term variability. To measure the relationship between these two indices, the value of SD1/SD2 is often used.

Also, unlike time and frequency domain methods, which identify known regulation mechanisms that respond to specific signal variations (the effect of breathing at high HF frequencies of RR and BP, or Meyer waves located at 0.1 Hz in the BP spectrum, for example,), nonlinear analyses permit to scrutinize the capacity of more complex regulatory mechanisms by the autonomic nervous system (Voss et al., 1995).

Among these, the analysis of fractality consists of quantifying the repetition of identical shapes but at different scales in the RR series (Peng et al., 1995; Bigger et al., 1996; Rajendra Acharya et al., 2006). Several approaches have been introduced to estimate the degree of fractality of HRV: the detrended fluctuation analysis (DFA) (Peng et al., 1995), the Hurst exponent, or the Higuchi (Higuchi, 1998) and Katz algorithms (Katz, 1988).

Entropy measures the complexity and regularity of patterns of different durations in the signal. Thus entropy will increase when a wide variety of patterns are identically distributed and decreases when these shapes are always the same. A number of algorithms have been developed to calculate these entropy indices: sample entropy, approximate entropy, Shannon entropy and the indices derived from it, conditional entropy, corrected conditional entropy, normalized corrected conditional entropy (Porta et al., 1998; Ferrario et al., 2004).

The symbolic dynamics allow to study short RR or BP variability pattern behavior (Porta et al., 2001). The principle of this method is to divide the signal into parts containing a fixed number of points (typically L = 3) and to digitize it (typically E = 6 levels). For each epoch obtained, the presence of predefined patterns is quantified (0V: no variations in the signal; 1V: two successive equal points; 2LV: two variations in the same direction; 2UV: two variations in opposite directions) from which their rates of occurrence in the complete signal are deduced.

Finally, we should also mention the largest Lyapunov exponent, which is used to measure how chaotic is the analyzed signal (Wolf et al., 1985). Its values will tend towards zero for physiological signals showing very small and slow variations, whereas they will increase when the variations are quickely larger and larger.

In clinical practice, BP variability usually represents from hour-to-hour to days, weeks, months and years variations (Nardin et al., 2019). In this context, an increase of variability is correlated with an increased risk of cardiovascular and all-cause mortality due to the alteration in cardiovascular regulatory mechanisms (Parati et al., 2013). However, very short-term, beat-to-beat BP variations have been analyzed and can be performed to assess cardiovascular regulation (Parati et al., 2020).

BP regulation is realized by different neurohormonal systems including the baroreflex system, the renin-angiotensin system, the shear stress-induced release nitric oxide from the endothelium, and the myogenic vascular response (Stauss, 2007). In addition to these regulatory mechanisms, independent and spontaneous variation of BP sympathetically mediated so called Mayer waves, which arise around 0.1 Hz, has also been demonstrated (Julien, 2020). These physiological systems induce variations in BP, which can be identified in a specific bandwidth of spectral analysis (Parati et al., 1990). In humans, myogenic vascular function affects both the very low frequencies (0.02–0.07 Hz) and the low frequencies (0.075–0.15 Hz); and sympathetic modulation can be seen only in the low frequencies and the LF/HF ratio (Stauss, 2007). The usual interpretation of high frequencies (0.15–0.40 Hz) is that they are mainly due to the mechanical effects of breathing on blood pressure (Pagani et al., 1986) but some authors have also shown a link with the endothelial derived nitric oxide (Stauss, 2007). The renin-angiotensin system and shear stress-induced release nitric oxide have been shown to affect very low frequencies in rats but need further investigation to be confirmed in humans (Stauss, 2007). The main indices of very short-term BP variability arise from spectral analysis of beat-to-beat systolic, diastolic or mean pressure; but there is a lack of studies using nonlinear methods such as fractal or entropy.

Cardiac coherence takes place in the general field of biofeedback which considers that physiological, cognitive and emotional systems are intimately interrelated through ongoing reciprocal communication (McCraty et al., 2006). Many methods have been developed, but globally they all try to measure the synchrony between psychophysiological interconnected systems. In this context, cardiac coherence evaluates the stability of heart rhythm in the low frequencies. Thus the cardiac coherence index will be represented by the presence of a peak in these regions visible on the spectral analysis of RR variability (Shaffer et al., 2014). From the power spectral density of the RR, the cardiac coherence ratio is defined as Peak power/[Total power−Peak power]; where Peak power is the integral in a window 0.03 Hz wide, centered on the maximum peak identified in the 0.04–0.26 Hz range; and Total power is the integral between 0.0033 and 0.40 Hz (McCraty et al., 2006).

Subjects can enhance their cardiac coherence capacity during training sessions composed of paced respiration around 0.1 Hz. In the autonomic nervous system, such training induces a resetting of the baroreceptor sensitivity, an increase in vagal afferent traffic and a reduction in sympathetic outflow (McCraty et al., 2006).

Due to its simplicity to deploy and interpret, HRV became a gold standard to evaluate autonomic nervous system activity. Thus it is commonly used for research, clinical as well as leisure purposes. Indeed, virtually all research protocols concerning the autonomic nervous system use HRV. As some studies have demonstrated its ability to predict fatal and nonfatal events, the method is also used in clinics, especially for cardiac and diabetic populations to follow at-risk patients. Also, amateur to professional sportspersons are now accustomed to employ HRV indices to manage their training loads and to prevent overtraining, mainly through the use of HR monitor watches.

Very short-term BP variability is less popular, mainly because of the cost of the material needed to record continuous BP waveform. However, these measures have been applied successfully to cardiovascular disease, such as hypertension, heart failure and stroke; particularly to estimate sympathetic modulation of vascular tone and impaired cerebrovascular myogenic function that are predictors of stroke (Smeda et al., 1999; Kiviniemi et al., 2014; Grilletti et al., 2018; Javorka et al., 2021). As a matter of fact in these populations, very short-term BP variability also constitutes a morbidity and mortality predictor (Lucreziotti et al., 2000).

Cardiac coherence is used as an alternative treatment for certain chronic diseases, stress or pain management (Shaffer et al., 2014; Fournié et al., 2021).

The sequence method aims to detect spontaneous sequences of cardiovascular changes (Parati et al., 1988; Parlow et al., 1995) i.e., consecutive increase of BP resulting in an increase in RR (positive sequence) or a consecutive decrease in BP, followed by a decrease in RR (negative sequence) (Reyes Del Paso et al., 2017). For each detected sequence, a regression line between SBP and RR is plotted and permits to calculate the slope of the curve. Then BRS is calculated as the average of all individual slopes expressed in ms/mmHg. The positive and negative BRSs can be expressed separately, and the effectiveness of the baroreflex can be calculated as the proportion in which the baroreflex is able to induce changes in response to the oscillation in SBP (Di Rienzo et al., 2001).

Some parameters can be set for the detection of sequences and the calculation of slopes: the minimum and maximum number of beat in a sequence (generally 3–6 beats), the lag between SBP and RR response (usually set to 0 to pick up the fast vagal arm of the cardiac baroreflex), the correlation coefficient of the regression line to validate a sequence (usually >0.85), and the threshold values for SBP and R-R interval changes in a sequence (typically, 1 mmHg and 1 ms, respectively) and the maximum successive R-R variations (commonly 20%).

BRS can be calculated by estimating the transfer function between beat-to-beat RR and SBP series (Robbe et al., 1987; Maestri et al., 1998; Pinna and Maestri, 2001; Pinna et al., 2002). The modulus gives the gain between SBP and RR (thus corresponding to BRS in ms/mmHg), and the argument provides the phase shift between the two signals (in radian but more generally expressed as a delay in seconds). In practice, these indices are calculated in the LF and the HF frequencies as their average in the two bandwidths (Clayton et al., 1995; Pinna and Maestri, 2002) and validated when the coherence function between RR and SBP is superior to 0.50 in the corresponding LF or HF bandwidths. Also, the delays between SBP and RR are considered as valid only if they suit with physiological timing of a working baroreflex, i.e., between 0.24 and 4 s (Keyl et al., 2001; Milan-Mattos et al., 2018).

The phase-rectified signal averaging method (PRSA) was introduced to overcome some weakness in the robustness of the above methods against noise and nonstationnarities present in both RR and SBP signals (Bauer et al., 2010).

The method consists of finding each SBP value that is higher than the previous one and to mark them as anchors. Then all the RR sequences centered at these anchors are aligned and averaged, thus constituting the overall HR response to SBP rises. The averaged signal allows to consider only the coupled RR–SBP signals and to eliminate noncoupled oscillations such as noise and artefacts. Finally, BRS is expressed as the magnitude of RR variation solely (in ms) (Bauer et al., 2010) or normalized by the average SBP increase that triggered the RR changes (in ms/mmHg) (Muller et al., 2012).

Conversely to HRV which provides a global status of the autonomic nervous system of a subject, BRS indices focus on the short-term neuronal regulation system of BP. Due to this specificity, such indices are preferably used on cardiac and hypertensive populations (Maestri et al., 1998; Pinna et al., 2015), in which it has been demonstrated to be a predictor of cardiovascular and all-cause morbidity and mortality (La Rovere et al., 1998; La Rovere et al., 2009; Bauer et al., 2010).

Assessing causality from the Granger approach is the most popular, including the ‘MVGC multivariate Granger causality’ toolbox (Barnett and Seth, 2014). Granger causality or classical Granger causality (cGC) is a popular tool for the user for assessing the presence of directional interactions between two time series of a multivariate data set (Wiener, 1956; Granger, 1969). However, cGC only includes the time-lagged effects between processes. In respiratory and cardiovascular physiology, significant instantaneous effects are present. Subsequently, cGC may lead to an incomplete description of the real phenomenon between processes. As a possible solution, the utilization of an extended model accounting for both instantaneous and lagged effects has been proposed. This modelling is named extended Granger causality (eGC). In a formal point of view, the cGC( i → j ) is the logarithmic measure of the ratio between the variance of the residuals U ∼ j of the restricted regression and the variance of the residuals U j of the unrestricted regression.

Accordingly, the eGC is the logarithmic measure of the ratio between the variance of the residuals W ∼ j of the restricted regression including the zero-lag effects matrix B 0 related to instantaneous effects (Faes, 2014; Schiatti et al., 2015) and the variance of the residuals W j of the unrestricted regression also including B 0 .

Notions of causality are commonly formalized in the context of a multivariate autoregressive (MVAR) representation of time series in order to allow time and frequency domain pictures (Faes et al., 2013a). Several frequency domain measures of causality have been introduced. Actually, measures to quantify causality in the frequency domain have been proposed from strictly causal MVAR representation: directed transfer function (DTF) (Kamiński and Blinowska, 1991), directed coherence (DC) (Baccala et al., 1998) and partial directed coherence (PDC) (Baccala and Sameshima, 2001). Faes and Nollo have extended DC and PDC measures to extended causal MVAR representation to provide extended directed coherence (eDC) and extended partial directed coherence (ePDC) (Faes and Nollo, 2010). A synthesis of all measures of causality and coupling has been proposed by Faes and Nollo (Faes and Nollo, 2011). Table 5 shows the frequency domain measures of causality presented in this paper.

Pharmacological and basic physiological manoeuver studies have given some information on the underlying mechanisms supported by each RR, SBP and RE causal relationships (Faes et al., 2011; Faes et al., 2013b; Faes and Nollo, 2011; Porta et al., 2013b; Porta et al., 2014; Porta et al., 2015). Then it has been demonstrated that the causal link between SBP and RR is due to baroreflex involvement, as the opposite link between RR and SBP has been hypothesized to be due to Starling law and arterial Windkessel effect (Porta et al., 2011). The coupling between respiratory centers and the cardiac vagal center in brainstem, the activation of cardiopulmonary reflexes and direct mechanical pulmonary stimulation of the sinus node tissue through inspiratory inflations, constitute the physiological mechanisms of the relationship between RR and RE (Saul et al., 1991; Eckberg, 2003). The causal relationship from RE to SBP is due to variations of stroke volume resulting from mechanical effects on intrathoracic pressure and venous return (Toska and Eriksen, 1993; Caiani et al., 2000). Finally, no physiological mechanisms were show for the causality link between SBP to RE (Porta et al., 2013b).

The main studies using Granger causality analyses concern the physiological validation of the methods and some mathematical improvements, such as the extended version of the Granger causality which is an adaptation of the method to physiological constraints (Faes and Nollo, 2010; Schiatti et al., 2015). Some authors tested successfully the causal method on anesthesia (Bassani et al., 2012; Porta et al., 2013a), syncope (Faes et al., 2013b; Faes et al., 2015; Schiatti et al., 2015), pre-eclampsia (Riedl et al., 2010), or fibromyalgia patients (Zamunér et al., 2017) but there remains a lack of major clinical studies using these indices as prognostic values.

Time and frequency domain, geometrical, empirical mode decomposition, moments and Lyapunov exponent, entropy, Poincaré plot, symbolic dynamics, and fractal indices of variability are displayed in Figure 4. The user can select the signal to analyze (RR, SBP, MBP, and DBP) using a popup menu, after which the results are automatically updated. Some graphs of the most recognized analyses are plotted. Additionally, the user can select the cardiac coherence analysis from the same popup menu (Figure 5).

The results provided by the three methods described above (sequences, transfer function and phase-rectified signal averaging) are displayed in Figure 6. For the sequence method, the user can exclude the outlier slopes manually or automatically (exclusion of the sequences out of mean ± 3*SD). Positive, negative and excluded sequences are plotted in blue, black and red, respectively. Also, the localization of each detected sequence can be displayed on an additional window by selecting the “Show seq.” pushbutton.

HRV, BPV, and baroreflex indices are calculated using wavelet transform, as described in the Signal processing section. Three graphs allow plotting changes in RR, HR, SBP, MBP, DBP and the LF, HF, LFnu, HFnu, and the LF/HF ratio for RR and BP signals along time (Figure 7). The displayed signals are selected from three popup menus on the left side of each graph. If entered, the events and selected areas are plotted and the user can zoom in on the data using a synchronized scale for the three graphs. On the bottom of this Figure, a popup allows the user to select the number of points of the sliding average performed on the signals.

The upper left part of Figures 8, 9 displays the original beat-to-beat signals and their modelizations (red color), while the upper right part gives the lag and the fit model values for each of them, allowing the user to evaluate the adequacy of the model for the calculation of Granger causality indices.

On the lower left part, two popup menus permit to select the type of Granger causality analysis to display (time or frequency domain and classical or extended) and the number of signals to consider (bivariate or trivariate analysis). Results are displayed in three manners: 1) a symbolic representation of the causality indices using arrows between each signal, with the thicker the arrow, the greater the causality and vice versa; 2) the numerical values of each causality; and 3) the causality indices represented as bar graphs (time domain) or spectrums (frequency domain) for easier comparisons.

This analysis displays the beat-to-beat evolution before and after an event of RR, HR, SBP, MBP and the variables Ptot, VLF, LF, HF, LFnu, HFnu, LF/HF alphaLF and alphaHF calculated by the wavelet method. If several events with the same name have been entered, the program calculates the average of the changes in the chosen variable synchronized around the event. In this window, the user selects the variable to be analyzed, the event around which the analysis is to be performed, and the durations to be considered before and after the event using the appropriate popup menu (Figure 10). When the user selects the “Plot data” pushbutton, the program searches all the sequences that meet the entered criteria and plots the results. The variables are resampled at 1 Hz to allow a regular 1-s step plotting. Comparisons are possible by adding multiple analyses using different HRV, BPV and baroreflex indices and/or different events on the same graph. Results can be plotted with all sequences or as the mean ± SD, and an outlier filtering can be applied.

All analyses described above can by computed in batch mode to avoid repetitive operations and to enable the analysis of a wide number of files. Before using this tool, users have to prepare all files as follows (all steps are described above): 1) import all data to create the .cvr files; 2) make all corrections and exclusions of the beat-to-beat data; and 3) set all the areas of interest and/or events with which the batch analyses will be perform. Then after having selected the File list analysis in the ‘Analysis’ menu or by pressing the corresponding toolbar, the batch analysis window opens. First, the user selects the folder containing the .cvr files to analyze. For each analysis, the user selects the options in the windows; for the analyses, the software will use the parameters set in the preferences window. After pressing the pushbutton corresponding to the desired analysis, the software will analyze all areas and save the results in a.txt file directly importable and ready for use in statistical software environments.