EasieRR: An open-source software for non-invasive heart rate variability assessment

1. The assessment of heart rate (HR) and heart rate variability (HRV) based on electrocardiograms (ECG) is considered a good proxy for stress in a wide range of animal species. However, problems can occur e.g., when measuring ECG in ambulatory


| INTRODUC TI ON
Measures of heart rate (HR) and heart rate variability (HRV) have a long history as indicators of cardiac health, and stress in humans (Carney, Freedland, & Veith, 2005;Delaney & Brodie, 2000;Fleisher, 1996;Hall et al., 2004), and have gained considerable interest in animal behavioural studies during the last 30 years (Hopster, Werf, & Blokhuis, 1998;Kovács et al., 2014;Mohr, Langbein, & Nürnberg, 2002;von Borell, 2000).HRV is based on the analysis of normal fluctuations in the time intervals of consecutive heartbeats = interbeat intervals (IBIs) or RR-intervals.It reflects the interplay between sympathetic (SNS) and parasympathetic nervous system (PNS) and thus also provides information about cardiac vagal tone (Porges, 1995;von Borell et al., 2007).Analysing HRV is regarded as a suitable approach to determine the activity of the autonomic nervous system (ANS) in the context of stress, affect and emotion (Appelhans & Leucken, 2006;Boissy et al., 2007;Thayer, Åhs, Fredrikson, Sollers, & Wager, 2012), and is progressively emerging as a suitable indicator of welfare states in farm animal research.
Furthermore, animals of the same species cope differently with environmental challenges.As the ANS has a major impact on the regulation of fundamental physiological functions related to coping and stress resilience, many of these phenotypic differences in stress response are mediated by different activation of the ANS.ANS activity in different autonomic phenotypes of free-living streaked shearwater Calonectris leucomelas was linked to stress level (Muller, Vyssotski, Yamamoto, & Yoda, 2018), consistent across time and context.Evans et al. (2016) used HR and HRV to investigate the impact of ecological mechanisms on the hibernation process in free-ranging brown bears.
These studies have speculated ANS may be a key mechanism driving phenotypic variation in animal populations, and is therefore a potentially important mediator in the evolution of life history (Evans et al., 2016;Muller et al., 2018).
The process of artefact correction of the electrocardiograms (ECG) signal has been shown to be essential for the appropriate analysis of HRV data (Berntson & Stowell, 1998;Shaffer & Combatalade, 2013).This makes HRV analysis especially challenging for data derived from unrestrained animals where, due to technical or physiological interference, a high occurrence of artefacts is likely.However, even when recording ECG during resting states, visual inspection is important to detect artefacts like atrioventricular blocks (AV blocks type 1-3), due to disturbances in impulse conduction at the heart.
These type of artefacts can be relatively commonly found in horses, but also in other animals and if not excluded, can result in very high HRV.In veterinary and behavioural animal research, human HR monitors, such as POLAR ® devices (POLAR Electro Oy) has been widely used in unrestrained animals.The major disadvantage of many of these devices has been that only RR-intervals were extracted but TA B L E 1 Common software available for ECG processing and HRV analysis compared to EasieRR

Analysis reports and export of results Availability
Kubios HRV (Premium)  1).However, these programs require experience in the field and can be time-consuming to learn to use correctly.
In this article, we describe an open-source software with an intuitive graphical user interface (GUI) that imports ECG from any csv file or txt files exported from AcqKnowledge software, detects the prominent peaks, makes correction of spurious detections easy and transparent, and calculates a variety of HRV parameters in the time-and nonlinear domain.

| E asieRR : COMPUTATIONAL BACKG ROUND AND THEORY
EasieRR is an open-source software developed to assist researchers in the use of HR parameters and their processing and analysis.
Special emphasis was put on an intuitive GUI to ease detection and manual correction of artefacts.
EasieRR has been programmed using MATLAB2018b (Mathworks) and compiled as a stand-alone application with the MATLAB compiler 7.0 for Microsoft Windows operating systems (Version 7 and upwards).EasieRR is hosted at https://figsh are.com/proje cts/Easie RR/68831 and distributed under the terms of the GNU General Public License.

| Peak detection
EasieRR applies by default a band pass filter (4.Order 0 degree Butterworth filter with low cut-off at 1 Hz and high cut-off at 20 Hz), hence removing high frequency peaks in the signal as well as removing any voltage offset.These two cut-off frequencies have been shown to be applicable for HRV analysis in goats, but can be modified by the user for optimal use on other species.The extraction of time intervals between heart beats is traditionally done by detecting R-peaks (RR-intervals), because they are usually the most distinct peaks (Lippmann, Stein, & Lerman, 1994).However, as long as variability within the HR cycles is calculated, other peaks and thus intervals i.e.S-S intervals can be used in case they are more prominent.It is essential that the type of peak used for detection is consistent within and between all individuals of an experiment (G.G.Berntson, pers. comm.).There are various algorithms for R-peak detection applied within different programs, i.e. global/local threshold detection (Berntson, Quigley, Jang, & Boysen, 1990), Pan-Tompkins QRS detector (Pan & Tompkins, 1985) or template matching (Friesen et al., 1990).In EasieRR, the peak detection used for heart cycle interval determination is based on the peak prominence.This allows for robust R-peak detection when combined with the band pass filtering and a predetermined minimum timespan between RR-intervals (i.e. heart refractory time of the species used).

| Artefact detection and processing
Currently available software uses different approaches to correct artefacts such as deletion of spurious RR-intervals, interpolation of missing or extra beats, i.e. cubic spline interpolation, linear interpolation, but also manual correction (see Table 1).To maximize transparency and reliability, EasieRR does not detect artefacts automatically.Instead, the GUI is specifically designed for easy visual examination and detection of artefacts.For this purpose, EasieRR is In contrast to most existing software which use interpolation algorithms for outlier correction (i.e.Kubios), EasieRR is based on a deletion algorithm.This method has been shown to be best suited for time-domain measures (Rincon Soler, Silva, Fazan, & Murta, 2017) and analysis of shorter sequences of ECG recordings (Lippman et al., 1994).The latter is often the case when data is obtained in ambulatory settings.Besides deleting outliers, EasieRR also allows selection of missed peaks and movement of misplaced marks.and the major axis illustrates SD2.Thus the major advantages of using Poincaré plots is their suitability as a quantitative visual tool allowing immediate recognition of artefacts (Myers, Workman, Birkett, Ferguson, & Kienzle, 1992) as well as estimation of the activity of the PNS (Kamen, Krum, & Tonkin, 1996;Woo, Stevenson, Moser, & Middlekauff, 1994).

| Analysis in the time and nonlinear domain
To validate EasieRR, we analysed different HRV parameters (Table 2) using identical ECG data lasting 20 s for the programs described in Table 1.While all programs calculated almost identical values for HR and mean-RR, there are small deviations in the time-domain and nonlinear parameters.This is probably due to the different underlying algorithms for artefact correction used in the different programs.

| PROG R AM DE SCRIP TI ON
The GUI is split into three separate interactive windows (Figure 5):  Abbreviations: bpm, beats per minute; ECG, electrocardiograms; HRV, heart rate variability; mean RR, mean interbeat interval; RMSSD, root mean square of successive differences between interbeat intervals; SD1, standard deviation along minor axis of the ellipse; SD2, standard deviation along major axis of the ellipse; SDNN, standard deviation of interbeat intervals.

TA B L E 2
Comparison of basic HRV parameters calculated using EasieRR and other commonly used HRV software.
(ECG data was the same for all software; X = not available in the software analysis) buttons with the options of either manually mark (delete) artefacts, move marks of spuriously detected peaks or insert marks for missed peaks.Their corresponding 'Undo' buttons allow the user to reverse their last action.After each artefact correction, the 'Recalculate' button will create a new tachogram as well as Poincaré plot.After artefact correction is finished the data range can be saved using the 'save range' button and viewed with the 'Check data' button.When the analysis is finalized, pressing the 'save to file' button will save all analysed time ranges.A more detailed user manual can be found in the documentation hosted here: https://figsh are.com/proje cts/Easie RR/68831.

| Synchronization with behaviour
Heart rate variability data can also be synchronized with observed behaviour via import it into common software for behavioural video analysis (e.g.The Observer ® ).More information is provided in the user manual available online.

| CON CLUS I ON AND FUTURE DIREC TIONS
EasieRR is a free software for evaluating ECG in non-restrained animals and allows for calculation of recommended standard HRV parameters in both, the time-and the nonlinear domains.
The intuitive GUI facilitates the detection and correction of artefacts through the visualization of RR-intervals in tachogram and Poincaré plot.
displaying the raw (and band pass filtered) ECG signal, a tachogram as well as a Poincaré plot simultaneously.Significant deviations between the lengths of successive interbeat intervals are often caused by spurious peak detection and artefacts.These divergent intervals are easily recognizable as distinct peaks in the tachogram and can also be seen as outliers in the typically elliptical shaped Poincaré plot.Figures1-3display in detail the procedure of artefact detection and correction using the three correction options: move mark, insert mark, mark outlier.
Assessment of HRV has been commonly done using time-domain parameters which are based on the calculation of successive RR-intervals, also called normal-to-normal (NN) intervals.EasieRR calculates the time-domain measures SDNN (standard deviation of NN-intervals), RMSSD (root mean square of successive differences between normal heartbeats) and mean HR.While HR is a good indicator for overall arousal or activity, it does not allow to draw inferences on the activity of the two autonomic nervous branches, the PNS and the SNS.In contrast, RMSSD has been found to reflect PNS-mediated HRV and can quantify the instantaneous beat-to-beat variance in HR.SDNN rather reflects the long-term variability of beat-to-beat intervals and is usually interpreted as an indicator of the sympatho-vagal balance (Task Force of the European Society of Cardiology & the North American Society of Pacing & Electrophysiology, 1996).

F
Example of artefact correction using 'mark outlier' button (cycle detection method = negative peak: correction is displayed on red marks and lines only).(a) Before correction: Circles mark signs of spuriously detected peak.(b) After correction: Circles mark location of changes in ECG (deleted data point marked with 'x' and tachogram (coloured yellow for deleted data).Outliers have disappeared in Poincaré plot In recent years, the use of nonlinear methods has been gaining more interest.Nonlinear methods like the Poincaré plot make it possible to include non-stationary data in the analysis (Guzik et al., 2007).A Poincaré plot is a scatter plot where each RR-interval is plotted as a function of the previous RR-interval.This typically results in an elliptical-shaped scatter plot tilted 45 degrees counter clockwise relative F I G U R E 2 Example of artefact correction using 'insert mark' button (cycle detection method = negative peak: correction is displayed on red marks and lines only).(a) Before correction: Circles mark signs of undetected negative peak in ECG, tachogram and Poincaré plot.(b) After correction: Circles mark location of changes in ECG ('+' for new data point) and tachogram (dotted line for new data).Outliers have disappeared in Poincaré plot to the x-axis (see Figure 4).From these data points, it is possible to calculate the standard deviation in two dimensions: Along the minor axis of the ellipse (SD1) and the major axis (SD2).SD1 is describing the instantaneous variability of the RR-interval time series (see RMSSD) and reflecting parasympathetic efferent activity at the sinus node.SD2 is describing the long-term variability of the RR-interval time series.F I G U R E 3 Example of artefact correction using 'move mark' button (cycle detection method = negative peak: correction is displayed onred marks and lines only).(a) Before correction: Circles mark signs of spuriously detected negative peak.(b) After correction: Circles mark location of changes in ECG ('+' for moved data point, 'x' for deleted data point) and tachogram (dotted line for new data).Poincaré plot is a more defined ellipse Commonly these two standard deviations are visualized within the Poincaré plot via an ellipse where the minor axis illustrates the SD1 an upper window displaying the ECG signal (raw and filtered) which allows to visually inspect automatically detected QRS complexes and to correct artefacts.The lower left window shows the corresponding tachogram facilitating the efficient location of artefacts and a Poincare plot at the lower right.The EasieRR GUI is operated via a range of buttons.The 'Choose data range' button allows the user to select a sequence to be analysed.Time-domain as well as nonlinear parameters are then automatically calculated and the corresponding Poincaré plot is generated in the lower right window.For artefact correction, the GUI offers three F I G U R E 4 Poincaré plot of RR-data series with typical elliptical distribution and SD1 (red) and SD2 (green)