Time-domain analysis is applied to quantify HRV using indices based on means or standard deviations, generally calculated over long-term recordings, typically 24 h (Kleiger et al., 1995; Balocchi et al., 2006; Rajendra Acharya et al., 2006). Because these indices represent short- to long-term variations in RR intervals, some indices are dependent on the recording length. It is therefore recommended to perform the calculations on recordings of similar duration to allow comparisons between measures or subjects.

The standard deviation of normal-to-normal intervals (SDNN) represents the variability over the entire recording period, giving the overall autonomic modulation regardless of sympathetic or parasympathetic arm. Some indices have been developed to quantify high frequency variations arising from parasympathetic activity such as the pNN50, calculated as the percent difference between adjacent normal RR intervals greater than 50 ms, and the rMSSD, calculated as the square root of the mean of the sum of the squared differences between adjacent normal RR intervals. The SDANN index, calculated as the standard deviation of the mean of all normal RR intervals for 5-min segments, quantifies the changes in heart rate due to cycles longer than 5 min. The SDNNIDX index, calculated as the mean of the standard deviation of all normal RR intervals for all 5-min segments, quantifies the changes in heart rate due to cycles shorter than 5 min (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996).

The geometrical indices are calculated on the sample density distribution of the RR intervals, which corresponds to the assignment of the number of equally long RR intervals to each value of their length (Malik, 1995; Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996; Rajendra Acharya et al., 2006). The main studies that have validated this method have used bins of 7.8125 ms width that were adapted for RR series obtained from EKG at a sampling rate of 128 Hz. To allow comparisons, this bin width must be the same, regardless of the EKG sampling rate. Thus, for sampling rates other than 128 Hz, the sample density distribution graph presents discontinuities, and is then reconstructed by smoothing the curve using a moving average function. The HRV triangular index is calculated as the integral of the density distribution divided by the maximum of the density distribution (Y). The histogram can be interpolated as a triangle, using the minimum square difference. The triangular interpolation of the RR interval histogram (TINN) is the baseline width of this triangle.

These measures quantify overall HRV, and are influenced mainly by slow—but not rapid—oscillations of RR intervals (Malik et al., 1989). This geometrical method provides a good-quality analysis but requires a reasonable number of beats for effective application.

The power spectral density (PSD) defined using Fourier transform is used to assess the different frequencies in the RR series. The spectrum is calculated using Welch's periodogram algorithm with a Hamming window of 256 points, an overlap of 50%, and at a precision of 256 points/Hz. The analyzed series can also be resampled from 1 to 8 Hz. The frequency indices are calculated as the integral of the PSD in specific bands of interest, as standardized by the Task Force in 1996 (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). After that time, frequency-domain analyses of HRV were recommended to be performed on short-term 5-min recording segments and on long-term segments of typically 24 h. The two durations involve different bands of interest. The default bandwidths used to calculate the total power (Ptot), ultra-low frequencies (ULF), very low frequencies (VLF), low frequencies (LF), and high frequencies (HF) are summarized in Table 1. For short-term analyses, LF and HF are also expressed as normalized values LFnu = 100^*LF/(Ptot-VLF) and HFnu = 100^*HF/(Ptot-VLF), and the LF/HF ratio is calculated (Table 1).

Total spectral power (Ptot) indicates overall HRV, and is used to assess overall autonomic cardiac modulation as well as the SDNN parameter. HF power represents short-term HR variation. Studies have shown that injected atropine completely eliminated HF power (Akselrod et al., 1981; Pomeranz et al., 1985). Thus, HF power is modulated by parasympathetic activity only, corresponding to peak respiratory rate. Pharmacological studies have shown that muscarinic cholinergic blocker and beta-adrenergic blocker lowered LF power, enhanced by dual blockade (Akselrod et al., 1981; Pomeranz et al., 1985). Both parasympathetic and sympathetic cardiac activity would therefore be associated with HR power in the LF band. Saul et al. (1990) and others (Pagani et al., 1997) showed a concomitant increase in LF power and muscle sympathetic nerve activity measured by microneurography. Furthermore, under atropine, LF power increased during orthostatic testing (Taylor et al., 1998). Although, these studies showed sympathetic cardiac modulation in LF power, changes in LF power can be interpreted only in relation to HF power. Accordingly, normalized indexes such as the LF/HF ratio, LFnu, and HFnu are used to examine this relationship.

To summarize, whereas HF power is modulated by parasympathetic modulation, LF power is controlled by both sympathetic and parasympathetic activity, and normalized indexes or the LF/HF ratio are used to estimate sympathetic modulation and autonomic equilibrium (Pagani et al., 1984) although this last index remains criticized (Billman, 2013).

None of the above-described methods allow temporal localization of sudden changes in RR signal behavior. To overcome these limitations, we applied wavelet transform (Pichot et al., 1999), which obtains a temporally localized sliding analysis of the signal, providing access to the heart rate variability at any time, for example, when the autonomous nervous system equilibrium is suddenly disrupted by an acute clinical situation or a physiological laboratory intervention. This method has been used successfully to assess autonomic reactivity to tilt tests (Jasson et al., 1997), exercise (Tiinanen et al., 2009), sleep apneas (Chouchou et al., 2014), experimental pain (Chouchou et al., 2011), generalized interictal EEG discharge (Sforza et al., 2014), pharmacological blockades (Pichot et al., 1999) and, anesthesia (Pichot et al., 2001).

For the analysis, a temporal sliding window of different weights (corresponding to different analysis levels or frequency ranges) containing the wavelet function is shifted along the signal. The weight characterizes a family member with a particular dilatation factor. The calculation gives a serial list of wavelet coefficients, which represent changes in the correlation between the RR signal and a given wavelet at different levels of analysis (or different frequency ranges) along the signal. The correspondence between wavelet coefficients and frequencies allows calculating the standard, ULF, VLF, LF, HF, LFnu, HFnu, and the LF/HF ratio along time. A complete description of this method, which is implemented in the software, is available in two previously published articles (Pichot et al., 1999; Wiklund et al., 2002).

An alternative, nonlinear approach was proposed to examine cardiac autonomic control (Voss et al., 1995). In recent years, interest has grown in the nonlinear dynamics that characterize autonomic cardiovascular control (Rajendra Acharya et al., 2006; Maestri et al., 2007b). Investigations of the complex feedback loops that impact the cardiac function have led to new indices designed to reflect signal complexity. We therefore incorporate certain of these indices that have proven useful for HRV interpretation and for health and mortality prediction (Maestri et al., 2007a; Huikuri et al., 2009).

Cardiac activity is controlled by the sympathetic and parasympathetic systems as well as hormonal system (Guyenet, 2006). These systems induce heart rate oscillations at different rhythms and HRV methods are used to study these rhythms and consequently autonomic cardiac modulations.

Slow and fast oscillations do not carry the same accuracy to represent sympathetic and parasympathetic activities. HRV analysis has been demonstrated reliable in assessing parasympathetic activity (Sayers, 1973; Akselrod et al., 1981; Pomeranz et al., 1985) while this method remains less accurate to assess sympathetic activity as heart rate slow changes are related partly to sympathetic and parasympathetic activity (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). A better approach of sympathetic activity relies on examining relative changes in the LF/HF ratio and LF%, or SD1/SD2 (Pagani et al., 1984; Tulppo et al., 1996; Malliani et al., 1998). Several studies have shown that relative changes in fast- and slow-oscillations allow approaching relative sympathetic activity in response to tilt testing under atropine (Taylor et al., 1998) or experimental pain (Burton et al., 2009; Chouchou et al., 2011).

Furthermore, some indices, such as temporal indices for example, depend on the length of the selected RR series to analyse while other are not, such as PSD. Ideally, except for time-frequency analyses, the data should be stationary (Chouchou et al., 2009), not less than 5-min duration or 250 beats, recorded over comparable time period, and in similar situations as sitting, upright, lying for example (Bahjaoui-Bouhaddi et al., 2000).

Also, interpretation of HRV indices remains especially debated in long-term recordings, because major determinants of HRV such as environmental factors, physical activity, and sleep duration vary significantly over time (Lombardi and Stein, 2011). Thus, HRV indices can be calculated on the whole RR sample or as the mean of successive epochs. Users might be aware of that and can refer to Table 2 for these specifications.

The software is optimized for analysing RR signal arizing from human recordings. R-peaks detection as well as HRV anayses can be performed on animal dataset by setting the adapted parameters in the preferences window. However, users must verify the accuracy of the methods for animals while some of them have been previously validated, such as Fourier analysis (Cerutti et al., 1991), others are not commonly use.

Finally, HRV analysis requires clean RR series while noisy recordings may lead to invalid physiological interpretations (Saul et al., 1988). An accurate preprocessing based on interpolation is usually utilize to correct for missing and extra beats.

Data can be imported from EKG or RR files. The available formats for EKG data are EDF, ISHNE, binary, text, Matlab. A generic import window allows to pre-visualize the EKG trace while the user sets the data file specifications. Preformatted RR files can be imported from text and binary files. Several preformated configurations are proposed (Suunto, Polar, Aria Delmar, and Vista Novacor Holters).

R peaks are detected from EKG signals using a laboratory-developed algorithm largely inspired by Chen et al.'s method (Chen et al., 2006). This moving average-based method, initially developed for real-time QRS detection, has demonstrated a 99.5% detection rate in the MIT-BIH Arrhythmia Database. First, the raw EKG signals are upsampled at 1000 Hz to ensure good time resolution for the RR intervals, even at a low initial EKG sampling rate. The raw EKG signal is first denoised using wavelet transform. The denoised signal is then passed through a moving averaging-based linear high-pass filter in order to highlight the QRS complex, which is processed by full-wave rectification and nonlinear amplification followed by sliding-window summation. The resulting pulse train (i.e., waveform) is used to calculate an adaptive threshold for QRS complex detection.

EKG recordings can contain ectopic beats—due to cardiac dysrhythmia—and missing or spurious beats—generally due to poor EKG signal quality. These nonsinusal beats must be corrected, because they induce errors in the calculation of HRV indices (Saul et al., 1988), which could result in invalid physiological interpretations.

In order to replace invalid beats, the software therefore provides automatic correction of the RR series, inspired by the algorithm developed by Kamath and Fallen (Kamath and Fallen, 1995). First, false beats are detected using Cheung's algorithm (Cheung, 1981): a high and low threshold are set for the relative variation in successive RR intervals (+32.5% and −24.5%, respectively). Second, for each detected error, the number of missing beats is estimated by comparing the total time duration within the error period with the duration of the immediately preceding beat. If the number of beats to recalculate is 3 or less, a cubic spline interpolation is done. This configuration generally originates from a single supraventricular or ventricular ectopic beat, or from an isolated R-peak missdetection. For 4 or more successive errors, the missing beats are interpolated by copying and inserting the same number of previous RRs between the first and last valid RR Although, these corrections generally results in a clean RR signal, this type of automatic algorithm can produce inconsistent values when the original EKG or RR signal is too corrupted, with lengthy portions of successive false beats. It is therefore recommended to visually inspect and review each series before performing HRV analyses. Additionally to RR correction, the software enables excluding parts of the signal from analysis, and the number of corrections is displayed in the results.

RR can also be manually corrected by two methods. First, the user can select an RR interval and change its value either by entering a new value or by using spline cubic interpolation to determine it. Second, when the EKG trace is available, R peaks can be manually inserted and removed by editing the EKG signal.

EKG recordings commonly contain excessive artifacted data that cannot be correctly interpolated. If the number of these data is low, HRV analysis can be performed on the remaining valid signals. HRVanalysis 1.0 allows excluding certain data from the analyses to avoid abnormal HRV values. This can be done by manually selecting the area to exclude or by using an automatic method based on successive beat-to-beat variations and RR labels (if present).

Calculations are performed according to the parameters set in the Preferences menu. For each type of analysis, a menu allows saving the results as a.txt file, saving the figure, and printing the figure. When saving the results, if the user selects an already existing.txt file, the data will be appended to the file. This allows combining analyses derived from different parts of the signal and/or different files and using the results directly in a statistical analysis program.

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.