Running head: HEART RATE VARIABILITY FOR PSYCHOLOGISTS 1 Heart Rate Variability for Psychologists: A Review of HRV Indices and an Analysis Tutorial

The use of heart rate variability (HRV) in research has been greatly popularized over the past decades due to the ease and affordability of HRV collection, coupled with its clinical relevance and significant relationships with psychophysiological constructs and psychopathological disorders. Despite the wide use of electrocardiogram (ECG) in research and advancement in sensors technology, the analytical approach and steps applied to obtain HRV measures can be seen as complex. Thus, posing a challenge to users who may not have the adequate background knowledge to obtain the HRV indices reliably. To maximize the impact of HRV-related research and its reproducibility, parallel advances in users’ understanding of the indices and the standardization of analysis pipelines in its utility will be crucial. This paper addresses this gap and aims to provide an overview of the most up-to-date and commonly used HRV indices, as well as common research areas that these indices have shown to be very useful, particularly in psychology. In addition, we also provide a step-by-step guide on how to perform HRV analysis using an integrative neurophysiological toolkit, NeuroKit2.

The physiological mechanism underlying HRV exemplifies the tight, yet complex coupling between the brain and the body. Heartbeats are initiated by the impulse generated from the sinoatrial (SA) node -a bundle of specialized cells in the heart. The SA node, dictating the rhythm of heart contractions, is in turn controlled by the complex interactions between different regulatory systems and mechanical activities such as respiration. As the interdependent sources of HR operate on different time scales, HRV is itself unevenly distributed in different frequency bands, with high and low frequency changes corresponding to different neurophysiological mechanisms. Rapid changes of HR reflect the cardiac regulatory influences of the autonomic nervous system (ANS), in tandem with its dynamic interaction with cardiovascular and respiratory activities Shaffer & Ginsberg (2017). On the other hand, slow HR changes reflect bodily fluctuations of lower frequency, such as circadian rhythms, sleep cycles, metabolism, and changes in body temperature and hormonal systems. The autonomic nervous system (ANS), which influences HR on short-time intervals, is subdivided into two distinct components, namely the sympathetic (SNS) and the parasympathetic nervous system (PNS). The SNS, known as the "quick response" system, predominates during elevated activity and stressful states and its stimulation results in an increase in HR. The PNS, known as the "relaxed response" system, predominates in the quiet and relaxing states, and is related to a decrease in HR. In healthy individuals, the two systems work together (promptly alternating between each other) to maintain the regulatory balance in physiological autonomic function.
In addition to the neural contributions of the two ANS systems, the short-term periodic changes of HR are also influenced by the cyclic activation of the baroreceptors -a mechanical sensor that works to maintain a homeostatic level of blood pressure (Berntson et al., 1997). When there is an increase in blood pressure, baroreceptors fire to decrease HR and increase the diameter of blood vessels, which subsequently brings the blood pressure back to normal. The opposite effect is observed when the blood pressure decreases. This ability of the body to regulate blood pressure is referred to as baroreflex. The stronger the baroreflex, the more frequently our HR is adjusted to maintain a stable level of blood pressure. Thus, a stronger baroreflex is associated with greater HRV.
The baroreflex is also the autonomic reflex that underlies the physiologic interactions between the cardiovascular and respiratory systems. Specifically, the arterial baroreceptors are also affected by breathing (arterial pressure decreases during inhalation and increases during exhalation) and results in variation in HR during each breathing cycle. This specific influence of respiration on HRV is referred to as the respiratory sinus arrhythmia (RSA). RSA can be seen as one of the aspects of HRV and encapsulates the coupling between respiratory and cardiac

HRV Indices
When HRV was first recognized as a potential indicator of ANS abnormalities (Hon, 1965), the main method used to quantify it was to measure the amount of variance in inter-beat intervals (IBI) -the time intervals between successive heartbeats. In 1981, Akselrod et al. (1981) introduced the power spectral analysis of HR to assess the power (or energy) found within different frequency bands. As different regulatory systems influence HRV at different frequencies, frequency-domain measures allow the assessment of their unique contributions to changes of HR. A few years later, L. Cohen (1989) introduced the time-frequency domain measures to capture both the time and frequency variabilities simultaneously. In the same decade, Goldberger and West (1987) argued that the HR signal was not simply a product of a regular periodic oscillator; instead, it displayed characteristics of a deterministic chaos -a seemingly random sequence of events with an underlying structure (Bigger Jr et al., 1996;Denton et al., 1990). Therefore, in order to capture the complex non-linear dynamic behaviours of the "chaotic" fluctuations of HR, new analytical approaches based on entropy and fractal dimensions were proposed (Bassingthwaighte et al., 2013;Bigger Jr et al., 1996;Huikuri et al., 2003;Lombardi, Sandrone, et al., 1996). Since then, the evolution of HRV research has been associated with the emergence of new and supposedly better indices or approaches to quantify and reliably assess HRV, making it a rapidly evolving field heavily connected to mathematical and computational developments.

Time-domain Analysis
Time-domain measures reflect the total variability of HR and are relatively indiscriminate when it comes to precisely quantifying the respective contributions of different underlying regulatory mechanisms (Acharya et al., 2006). However, this "general" sensitivity can be seen as a positive feature (e.g., in exploratory studies or when specific underlying neurophysiological mechanisms are not the focus). Moreover, as they are easy to compute and interpret, time-domain measures are still among the most commonly reported HRV indices. SDNN, measured in the unit of milliseconds (ms), is the standard deviation of the normal beats intervals (denoted NN intervals), i.e., excluding beats related to technical artefacts and physiological artefacts (e.g., like ectopic beats -beats that do not arise from SA node depolarizations -and atrial fibrillation). Mathematically, SDNN is the square root of the total variance in the entire recording and therefore, it encompasses both short-term high frequency variation (mostly due to parasympathetically-mediated RSA; Shaffer and Ginsberg (2017)) and long-term low frequency components of the HR signals. Some proponents argue that the accuracy of SDNN is higher for longer recording periods, especially over 24-hour periods where values of SDNN have been used for cardiac risk stratification in medical settings (Shaffer & Ginsberg, 2017).
Stemming from the concept of SDNN, there are other statistical variables that similarly capture the variance of NN but in a more fine-grained manner. Firstly, the standard deviation of average NN is extracted for each 5-minute segment of the long-term HR time series (e.g., 24hour recordings) (SDANN). It is highly correlated with SDNN and has been said to provide little additional information beyond what SDNN is able to (Shaffer & Ginsberg, 2017). On the other hand, the SDNN index (SDNNI) is the mean of the standard deviations of NN calculated per 5minute segment (Shaffer & Ginsberg, 2017). While SDANN is an estimate of long-term components contributing to HRV, SDNNI assesses the variability within the short intervals and hence is an estimate of short-term components that modulate HR signals.

Deviation-based Indices.
While these traditional indices have been recommended by official guidelines (Force, 1996) and have been extensively used in past research, new metrics are being explored in parallel with the development of the mathematics and statistics fields. Some of them involve certain forms of normalization, such as the Coefficient of Variation of NN (CV; i.e., the ratio of the standard deviation to the mean), which could be more stable when comparing individuals or conditions with different HR baselines (Abdi, 2010). Others involve "robust" alternatives to traditional dispersion indices, such as the Median Absolute Deviation (MAD; a median-based index of variability) or the interquartile range (IQR), that are less influenced by outliers and thus, less sensitive to measurement errors (Almazaydeh et al., 2012;Yilmaz et al., 2010).
While the aforementioned indices are calculated directly from the NN sequence, other approaches are derived from the difference between successive NN intervals. These measures include the standard deviation of the successive difference (SDSD), the root mean square of successive difference (RMSSD), the proportion of successive NN intervals that are larger than a given threshold (e.g., 20 ms or 50 ms -pNN20 and pNN50 respectively). As these measures are derived from beat-to-beat differences, they are mostly indicative of short-term variations in HR.
HRV can also be quantified from geometric (i.e., graphical) representations of NN, in the form of histograms and scatter plots. Such indices include, among others, the HRV Triangular Index (HTI), the triangular interpolation of NN histogram (TTIN) and other geometrical properties of the NN scatter plots (e.g. Lorenz plots). While geometric methods are less sensitive to measurement errors as compared to basic dispersion-based indices (M. Malik et al., 1993), a large number of NN intervals is required for reliable assessment. Therefore, even though 5minutes recordings have successfully been used in past research (Jovic & Bogunovic, 2011), recordings of at least 20 minutes are recommended for accurate computation (Force, 1996).

Frequency-domain analysis
As the different regulatory systems modulate HR at distinct frequencies, frequencydomain indices that reflect the distribution of power across different frequencies bands are particularly suited to assess the specific components of HRV. According to Force (1996), the HRV power spectrum can be meaningfully divided into four bands: ultra-low frequency (ULF; ≤0.003 Hz), very low frequency (VLF; 0.0033-0.04 Hz), low frequency (LF; 0.04-0.15 Hz) and high frequency (HF; 0.15-0.4 Hz).
Generally, due to the inherent differences in their signaling mechanisms, the parasympathetic (vagal) modulation of HR, including RSA, is faster than the sympathetic activities. Thus, research suggests that HF and LF are prominent reflections of parasympathetic and sympathetic activity, respectively (Acharya et al., 2006;Seyd et al., 2008). The ratio of LF to HF power (LF/HF) has also been commonly used as an index of the sympathovagal balance -a measure of the relative contributions of SNS to PNS activity.
HF component is relatively established as an index of the parasympathetic modulation of HR and under controlled conditions, its derivative, LnHF (the natural logarithm of HF) as an estimator of vagal tone (Acharya et al., 2006;Shaffer & Ginsberg, 2017). However, the concept of tagging a frequency component to a particular division of ANS activities has been controversial for other components. Particularly, the interpretation of LF power as an index of sympathetic activities has been repeatedly challenged. Research has shown that LF power is modulated by both ANS branches as well as baroreceptor activities (Acharya et al., 2006;Force, 1996;Goldstein et al., 2011). There are also mixed findings regarding the exact physiological mechanisms underlying ULF and VLF components (Kleiger et al., 2005;Shaffer & Ginsberg, 2017) with the former often being associated with very low-frequency biological processes like circadian rhythms or metabolisms (Barrett et al., 2001;Takabatake et al., 2001) and the latter with body temperature regulation and vasomotor activity (Bouzida et al., 2009;Coenen et al., 1977;Machado et al., 2015).
Another limitation of the frequency-domain approach is pertaining to the inherent differences in the outputs of different spectral analysis algorithms, which makes it difficult to compare the absolute values of frequency components across studies (Force, 1996; K. K. . For instance, a non-parametric method such as the fast Fourier transform (FFT) periodogram has been argued to overestimate the LF/HF as compared to the Lomb-Scargle (LS) periodogram (K. K. Kim et al., 2009). Even with the same algorithm, the values of the frequency components can greatly vary depending on the choices of the parameters such as the length of the estimation window or the model order (Marek Malik et al., 1996).
To allow for more reliable comparison of frequency components across studies, the spectral indices are often normalized by dividing the individual component by the total power. The normalization process allows the normalized components (normalized LF, LFn and normalized HF, HFn) to be less varied by the methodological differences and less sensitive to the changes in total power (Burr, 2007;Force, 1996). Depending on the length of the recording, the normalizing denominator (total power) can be with or without the ULF component (sum of power includes ULF, VLF, LF and HF for 24 recordings, and includes VLF, LF and HF for shortterm recordings; Shaffer and Ginsberg (2017)). Even though the normalized spectral HRV indices can be interpreted in the same manner as the absolute values, the general guideline recommends reporting the low-frequency and high-frequency components in both absolute (e.g. LF, HF) and normalized forms (e.g. LFn, HFn) to present a complete picture of the power distribution (Force, 1996).

Time-frequency methods of analysis
The spectral analysis algorithms used in frequency-domain approaches can only compute the power spectrum density (PSD) over the entire recording, failing to account for the temporal variation of the frequency components. To address this limitation, tools to simultaneously capture both time and frequency changes (i.e., a time-frequency approach, TF), have been proposed (L. Cohen, 1989Cohen, , 1995. By evaluating the frequency spectrum over short moving windows, TF methods quantify the instantaneous PSD at any given time (Force, 1996). The TF methods include, among others, linear approaches such as the short-time Fourier Transform (STFT), the wavelet transform, and quadratic approaches such as the Wigner Ville Distribution (WVD).
As it allows for the integration over both the time and frequency domain, the TF HRV indices usually include the maximum and minimum energy of the time windows in HR signal, the standard deviation between energy of time windows, the total energy of signal in different frequency bands (e.g. energy of VLF, LF or HF), and the corresponding average energy of each band (derived by dividing total energy by length of band). While the TF approaches are believed to be more informative since they provide a means to quantitatively assess the changes in frequency components across time, they are not as widely used as other domains. One current limiting factor is that some of them are computationally expensive, while others, being easier to compute, can hardly provide high accuracy in both time and frequency domain (Elsenbruch et al., 2000;Mainardi, 2009). The appeal of TF approaches could be greatly enhanced by addressing these technical challenges, together with future works that investigate the interpretations of relevant TF features extracted from the time-varying spectrum.

Non-linear dynamics
Recent research has emphasized non-linear interactions of physiological systems underlying cardiovascular regulation (Voss et al., 2009). Hence, HRV indices derived from nonlinear methods were proposed to adequately characterize the dynamical properties of HR that other methods are unable to (Force, 1996).
The Poincaré plot, whose principles were derived from nonlinear dynamics theory, is one of the most common methods to measure non-linear HRV (T. Henriques et al., 2020). It corresponds to a scatterplot of each NN interval plotted against its corresponding preceding interval, which approximates the cardiac system's evolution (Brennan et al., 2001).
The points are dispersed around the identity line and converge into an ellipsoid configuration.
Points above the line represent HR decelerations (NN intervals that are longer than preceding ones) and points below the line of identity indicate HR accelerations (Brennan et al., 2002).
Qualitatively, the Poincaré plot's shape provides a visual summary of the heart's behaviour (T. Henriques et al., 2020). From the plot, outliers such as premature heartbeats and Poincaré plot.
other technical artefacts can be easily removed, something which cannot be done with spectral and most time-domain analyses. On a quantitative level, indices extracted from the plot correspond to time-domain indices. SD1 is the standard deviation of points perpendicular to the identity line and is equivalent to RMSSD (Ciccone et al., 2017), describing short-term NN variability. SD2 is the standard deviation of points parallel to the line of identity and is equivalent to SDNN, representing long-term NN variability. The former reflects short-term vagal activity whereas the latter represents sympathetic modulation. A third index, SD1/SD2, describes the ratio of short to long term variations in NN interval fluctuations and reflects sympathovagal balance (Hsu et al., 2012). SD1/SD2 has been shown to additionally capture complexity in HR patterns (Claudia et al., 2003;Hoshi et al., 2013; e.g. healthy heart displays a "comet" shape Poincaré plot while cardiac abnormalities display atypical "fan" or "complex" shapes, Woo et al., 1992).
However, SD1/SD2 seems to be sensitive to the time lag used, with some recent research recommending a lag of 5 or 6 heartbeats (Claudia et al., 2003;Koichubekov et al., 2017).
Nevertheless, there is no consensus on the optimal time lag yet, which itself could be subject to inter-individual differences.
Entropy-based methods find their roots in information theory in which they are conceptualized as measures of complexity (or "unpredictability") of a signal. While they originated from fields of thermodynamics (Delgado-Bonal & Marshak, 2019), their applications have expanded towards the quantification of complexity in physiological systems and have been used widely as a diagnostic tool in biomedicine. The underlying concept of entropy is that it quantifies repetitions of patterns in a signal, with greater entropy indicating higher randomness and unpredictability (and thus, complexity) and lower entropy values implying that the cardiac system is predictable, (i.e., periodic and does not produce much new information) (Faes et al., 2019). Common entropy-based methods include Approximate Entropy (ApEn) (Pincus, 1991),

Entropy.
Sample Entropy (SampEn) (Richman & Moorman, 2000), and Multiscale Entropy (MSE) (Costa et al., 2002).SampEn was built on ApEn (Richman & Moorman, 2000) and addresses the latter's tendency to overestimate the amount of regularity in the signal (Al-Angari & Sahakian, 2007;Delgado-Bonal & Marshak, 2019;Xie et al., 2008), being in turn more consistent in measuring complexity across different time series lengths (T. Henriques et al., 2020). However, another important difference is that ApEn computes regularity based on short time sequences whereas However, as higher irregularity is not always equivalent to higher complexity (Costa et al., 2003), these traditional measures have been shown to falsely quantify random HRV changes in certain pathological conditions as being complex (Costa et al., 2002(Costa et al., , 2005. In order to estimate complexity from irregularity statistics more accurately, Costa et al. (2005) suggested to instead examine irregularity in multiple spatio-temporal scales where SampEn is estimated on more than one scale, giving rise to a new entropy-based method called MSE. Compared to single-scale entropy estimates, MSE has been able to accurately account for lower HRV observed in different cardiac abnormalities, such as atrial fibrillation where there are random outputs in HRV (Costa et al., 2005). As MSE formulation does not allow for reliable computation of entropy in short time windows, new methods such as refined MSE (RMSE) (Wu et al., 2013) and refined composite MSE (RCMSE) (Wu et al., 2014) have been introduced, though few studies have tested them on HRV (Blons et al., 2019;Deschodt-Arsac et al., 2020).
The entropy analysis of HRV is an actively evolving field, with more novel approaches being proposed, each claiming to confer certain advantages over more traditional approaches. For example, new measures like Fuzzy Entropy (FuzzEn) (Chen et al., 2007) and Fuzzy Measure Entropy (FuzzyMEn) (Liu et al., 2013) that were developed to overcome the poor statistical stability of ApEn and SampEn (Castiglioni & Di Rienzo, 2008;Humeau-Heurtier, 2015), have been gaining traction. While a comparison of clinical validity across these entropy measures has been demonstrated in a few studies (e.g., Graff et al., 2012Graff et al., , 2015Liu et al., 2013), the choice of entropy type and the appropriate parameters for entropy calculation is not a straightforward process and investigation into how these factors influence HRV interpretations is still ongoing (Mayer et al., 2014). where a higher CD reflects greater complexity in HRV. Other algorithms for computing fractal dimension also exist, such as the Box-counting dimension, and the algorithm by Katz (Katz, 1988) and Higuchi (Higuchi, 1990) but CD is the most popular by far (T. Henriques et al., 2020).
In recent years, the most extensively used method in measuring HRV complexity is Detrended Fluctuation Analysis (DFA) (Peng et al., 1995). DFA is a measure of the correlations between successive NN intervals at different time scales, and is used to analyze recordings spanning several hours. Its computation disregards nonstationarities in the cardiac time series and Fractal Measures.
hence any spurious correlations due to artifacts of nonstationarity behaviour (e.g., trends due to external stimuli) is avoided, an advantage conferred over other fractal methods (Echeverria et al., 2003;Penzel et al., 2003). DFA eventually extracts both long-term correlations and short-term correlations, with anywhere between 4 and 16 heartbeats for the former, and 16 and 64 heartbeats for the latter (Penzel et al., 2003). The short-term correlations correspond to the baroreceptor reflex whereas the long-term correlations measure regulatory mechanisms in the cardiac system (Shaffer & Ginsberg, 2017). See Table 1 for a summary of all indices. RSA refers to the fluctuations of HR at the respiratory frequency, or in other words, the effect of respiratory-gated PNS activity on the heart. As RSA commonly occurs at highfrequency heart rate fluctuations, most studies quantify RSA using the HF component of the spectral measures (Grossman & Taylor, 2007). Other methods to quantify RSA include the Porges-Bohrer method (PB) (Porges & Bohrer, 1990) and the Peak-to-trough algorithm (P2T) (Grossman et al., 1990;Katona & Jih, 1975). PB extracts RSA by removing variances in the HR time series that are not attributable to spontaneous breathing frequencies, and has high accuracy only when there is low signal-to-noise ratio in the data. On the other hand, P2T derives RSA by finding the difference between the shortest heart period during inhalation and the longest heart period during exhalation per single breath.
Even though the three different RSA indices are highly correlated (Grossman et al., 1990), PB has been proposed to be the most sensitive to vagal activity.

HRV in Psychological Research
As mentioned above, the different HRV indices are believed to reflect different physiological mechanisms. It is therefore unsurprising that they are associated with distinct psychological constructs, such as cognitive abilities, neural processes, personality traits and psychopathological disorders.

HRV and Self-Control
Based on the neurovisceral integration model (Thayer & Lane, 2000, the neural circuitry implicated in the effective regulation of psychological processes (i.e., emotional and cognitive self-regulation) are also related to the regulation of cardiac autonomic activity. The model has gained ample support over the past decades as mounting evidence emerges, showing a significant relationship between HRV and prefrontal cortex activity (Friedman, 2007;Thayer & Brosschot, 2005;Thayer & Friedman, 2002;Thayer & Lane, 2009). This is in line with studies investigating the relationship between HRV and executive functioning, showing that participants with higher resting HRV (particularly short-term changes) performed better at different executive tasks, such as sustained attention (Porges & Raskin, 1969;Richards & Casey, 1991), working memory (Hansen et al., 2003(Hansen et al., , 2009 highlighted the executive-specific advantage of higher vagally-mediated resting HRV which was not observed in non-executive tasks based on simple reaction time. Autonomic flexibility (i.e., the ability of ANS to modify bodily states), indexed by HRV, is crucial for adjusting between high and low arousal states, and has thus been related to emotional regulation. In studies that investigated emotional responses to stress, individuals with higher resting RSA had lower perceived (negative) emotional arousal and made greater use of adaptive coping strategies (Fabes et al., 1993;Fabes & Eisenberg, 1997). Similarly, individuals with higher HF experience more positive moods, in particular cheerfulness and calmness, and this relationship was mediated by emotion regulation abilities (e.g., ability to engage in situational planning, refocusing, and reappraisal) (Geisler et al., 2010). This parallels findings reporting that anxious participants experienced greater decrease in RCMSE during stress exposure (Blons et al., 2019). This is in line with neuroimaging findings highlighting the relationship between HRV and the activity of brain regions involved in threat perception and saliency detection such as the amygdala, media prefrontal cortex . Although the evidence in favour of a link between cognitive and affective dysregulation and reduced HRV is robust, the exact physiopsychological mechanisms underlying this association remain to be uncovered (Appelhans & Luecken, 2006).
While most of the literature focuses on the role of resting HRV in cognitive control and emotional regulation, another level of analysis of HRV, the phasic HRV, has been less examined but is also of great interest for health and well-being. Resting or tonic HRV is measured at the state of resting, whereas phasic HRV is measured during the state of stress where changes in HRV, relative to baseline, represent the reactivity of the cardiac system to different stimuli.
While the advantages of higher tonic HRV have been repeatedly shown in the literature , the desirability of phasic HRV changes (increase versus decrease phasic HRV) is context-dependent (G. Park et al., 2014). Researchers have highlighted that higher phasic HRV suppression (T. P. Beauchaine et al., 2007;G. Park et al., 2014;Schwerdtfeger & Derakshan, 2010; autonomic stress response to prepare for fight-or-flight situations, Thayer et al., 1996) act as a protective factor during stress exposure as it has been linked to higher levels of tonic HRV and phasic HRV post-stress, indicating faster recovery (Benvenuti et al., 2015;El-Sheikh et al., 2011;Rottenberg et al., 2005;Weber et al., 2010). However, when the stressor requires the engagement of executive functioning Segerstrom & Nes (2007), higher phasic HRV enhancement is preferred as it indicates greater attentional resources and self-regulatory effort (Butler et al., 2006;Laborde et al., 2014Laborde et al., , 2015Segerstrom & Nes, 2007). Interestingly, by looking at participants' performance under different levels of perceptual loads, G. Park et al. (2014) observed that individuals with higher tonic HRV can exert regulatory effort (presence of phasic HRV enhancement) and avoid making autonomic stress response (absence of phasic HRV suppression) even under high loads where the task could be a potential stressor. In contrast, individuals with lower tonic HRV experienced autonomic stress response (phasic HRV enhancement) even in low loads condition, suggesting a tendency to appraise even trivial challenges as stressors. Given the significant interaction between tonic and phasic HRV and its potential to provide insight on individual differences in regulatory processes, future studies should investigate not only HRV at rest but also HRV reactivity (changes between rest and event) and HRV recovery (changes between event and post-event) (Laborde et al., 2017).

HRV and Psychopathology
Research focusing on the relationship between HRV and mood disorders suggests that, in general, a lower vagal tone (and thus, a lower HRV) is associated with negative outcomes (T. Beauchaine, 2001). Short-term HRV indices such has RMSSD, SDANN, HF and non-linear measures like MSE were found to be significantly reduced in patients with depression (Kemp et al., 2010;Leistedt et al., 2011;Schulz et al., 2010) and negatively correlated with symptom severity (Koenig et al., 2016). Short-term HRV indices are also lower (and correlated with symptoms severity, Piccirillo et al., 1997) in anxiety disorders, including Generalized Anxiety Disorder (GAD), Social Anxiety Disorder, Panic Disorder and Post-traumatic Stress Disorder (Chalmers et al., 2014;Dalack & Roose, 1990;Dimitriev et al., 2016;Hartmann et al., 2019;Kemp et al., 2010;K. Kim et al., 2016;Klein et al., 1995;Lakusic et al., 2007;Prasko et al., 2011;Stein et al., 2000;Thayer et al., 1996;Tulen et al., 1996). Additionally, depression with comorbid anxiety has been shown to be associated with greater reduction in HRV as compared to depression alone (Kemp et al., 2012). In this context, lower HRV has been associated with changes in the sympathovagal balance in favour of a sympathetic predominance, which may subsequently increase cardiovascular risk in these pathologies (Dimitriev et al., 2016). Bipolar Disorder is also associated with a decreased HRV across a variety of HRV indices, including RMSSD, pNN50, SDNN, HF, LF, and SampEn (H.-A. Chang et al., 2014;H. Cohen et al., 2003;Henry et al., 2010;Lee et al., 2012). A similar pattern would be expected with Borderline Personality Disorder, given the shared physiopsychological mechanisms with bipolar disorder (Carr et al., 2018). However, the few studies currently report mixed findings Weinberg et al. (2009). Research also suggests a lower HRV with reduced RMSSD, HF, ApEn and SampEn in patients suffering from schizophrenia (K. Bar et al., 2005; K.-J. Bar et al., 2007;J. S. Chang et al., 2009;Henry et al., 2010;Jindal et al., 2005), and a correlation with symptoms such as apathy, withdrawal and delusions (Valkonen-Korhonen et al., 2003), suggesting a relationship between psychotic symptoms and cardiac autonomic function (Henry et al., 2010;Toichi et al., 1999).
HRV has been presented as a useful clinical marker of psychopathology, in which it represents the combined effect of parasympathetic inhibition as well as excessive sympathetic activation, reflecting a sub-optimal adaptive response to external and internal changes. Although the specific pathophysiological mechanisms underlying the associations between psychiatric conditions and HRV remain poorly understood, the majority of studies observed significant relationships between the different pathologies and parasympathetic parameters (i.e., such as HF and RMSSD that are particularly sensitive to short-term HR changes). This has led researchers to propose that a restrained vagal activity could be a transdiagnostic biomarker across psychiatric disorders (K. Bar et al., 2005;Chalmers et al., 2014;Einvik et al., 2011;Pavlov & Tracey, 2012).
This hypothesis is further supported by the observation of worsening parasympathetic tone during the chronic course of psychiatric conditions (K. Bar et al., 2005;Moon et al., 2013), and parallel improvements in parasympathetic parameters as symptom severity reduces over the course of treatment (Hage et al., 2017;Hartmann et al., 2019;Kemp et al., 2010;S. Park et al., 2018).

HRV as a Therapeutic Target
With the progressive establishment of HRV as a convenient peripheral (and correlational) index of self-regulatory abilities (Porges, 2007;Thayer & Lane, 2000), its usage in applied settings grew. Recently, the psychophysiological cardiac coherence model has been introduced to promote training to enhance HRV. This model emphasizes the importance of having reciprocal interactions and coherence between the regulatory systems to maximize the adaptability to environmental demands (McCraty et al., 2009;McCraty & Zayas, 2014). Cardiac coherence refers to a stable state of increased synchrony between the rhythms of different oscillatory systems, typically HR, respiration and blood pressure. In practice, this coherent state is reflected by a more ordered and sine wave-like HRV waveform, oscillating at a frequency approximating 0.1 Hz, the frequency at which the phase synchrony of HR and respiration occurs (P. Lehrer et al., 2013;E. G. Vaschillo et al., 2006).
According to the model, cardiac coherence can be induced by breathing at a steady pace approximating 0.1 Hz (i.e., 6 breaths per minute). Such deliberate slow-paced breathing (also referred to as "deep breathing," P. Lehrer et al., 2013) maximizes both baroreflex sensitivity and RSA on HR, giving rise to a high-amplitude oscillation of HR (P. M. Lehrer et al., 2003;E. Vaschillo et al., 2002). According to P. M. Lehrer et al. (2003), regular breathing exercising techniques help to continuously stimulate and exercise autonomic reflexes, especially the baroreflex. In fact, deep breathing has been shown to produce significant improvements in investigations are needed to understand, clarify and validate its potential clinical utility.

HRV in Practice: A Tutorial using Python
Over the years, HRV has become a convenient and reliable marker to detect ANS and cardiovascular abnormalities. However, due to the complex mathematical models and theories of signal processing involved, there is a practical barrier for users without the prior technical knowledge, such as many psychology students and researchers, to effectively understand, interpret, produce and reproduce HRV indices. In the following, we will present a step-by-step guide to obtain HRV indices using NeuroKit2 (Makowski et al., 2021), an open-source Python package that is designed to be friendly to both novice and advanced programmers. Instructions on how to get Python and NeuroKit2 up and running can be found at https://neurokit2.readthedocs.io/en/latest/installation.html.

Data Preprocessing
Cardiac activity data stored on the disk can be loaded using functions appropriate to the format (e.g., read_bitalino()or read_acqknowledge()). For the purposes of this tutorial, we will download and use an example dataset available in NeuroKit2, which consists of recordings of 4 participants during 8 minutes of resting state. Their cardiac activity (ECG) and respiration (RSP) signals were recorded and downsampled to 200 Hz for the purpose of this demonstration (note that higher sampling rates are recommended to maximize the precision of NN intervals).
The data from the first participant is then extracted, stored into an object, and two of its variables (the ECG and RSP signals) are passed to a preprocessing function that will clean the signals and, importantly, find the location of R peaks and respiration cycles. Details about the options and algorithms used are available in the package's documentation. ## array ([165, 345, 506, 672, 838]) The bio_process() function allows simultaneous processing of different signals using signal-specific pipelines. For instance, for ECG signals, the steps involve minimizing noise in the data, removing artifacts due to environmental or biological sources, as well as extracting all potential peaks (e.g. R-peaks) in the signal before removing peaks that are not of interest (i.e., ectopic peaks).
NeuroKit2 offers flexibility in choosing different processing methods. For this tutorial, we will use the bio_process() function which allows for less customization of the pipeline and relies on sensible default parameters for performing signal-specific operations (i.e., cleaning, peak extraction, rate computation). Alternatively, users can customize their own processing pipeline following the tutorial available on the documentation of NeuroKit2.

HRV Analysis
HRV indices of different domains can then be obtained using specific functions such as hrv_time(), hrv_frequency(), hrv_nonlinear() and hrv_rsa(), to which the list of R-peaks location, the sampling rate (200 Hz in this example), and respiration phases (specifically for hrv_rsa()) have to be provided. However, all the available indices can also be computed at once using hrv(), alongside a figure combining useful diagnostic plots (such as the Power Spectral Density and the Poincaré plot) when the show argument is set to True.

Multiple Participants Analysis
As the function(s) to compute HRV indices for a participant return a single row dataframe, the process can easily scale up for studies that contain multiple participants by, essentially, binding all the rows together into one dataframe. In the example below, we will demonstrate how to loop through participants, compute HRV indices at each turn and append the data. For more complex or tailored analyses, it is within such loops that additional participant-specific operations can be inserted (e.g., apply specific data preprocessing steps or extract other types of information  Acharya et al., 2006;Saboul et al., 2014). Additionally, although the use of LF/HF has also been quite popular amongst researchers, definitive conclusions about sympathovagal balance based solely on this parameter are ill-advised due to the unclear sources of LF.
Finally, while both time-and frequency-domain indices are correlated with non-linear indices, the latter are more reflective of the general complexity of HR signal rather than any specific underlying physiological mechanism (Sassi et al., 2015). However, theoretical and methodological concerns related to their usage have been raised. Firstly, these indices are very sensitive to the parameters chosen (e.g., window length), with unclear guidelines, rules, or standards for their (optimal) selection. Secondly, while some studies attest to the clinical utility of non-linear HRV dynamics when applied longitudinally (i.e., ability to differentiate progression of cardiovascular diseases), it is not yet clear why certain non-linear methods are better at distinguishing between pathologies. Apart from methods covered in this paper, other methods such as largest Lyapunov exponent (LLE), Hurst exponent (HE) and Symbolic Dynamics have also gained wide acceptance for assessing various complex systems, but their applications in HRV are substantially less popular than the aforementioned non-linear indices (T. Henriques et al., 2020). All in all, it is recommended to treat these indices as complementary measures to traditional linear indices (Laborde et al., 2017;Sassi et al., 2015).
Finally, other considerations influence the choice and interpretation of HRV indices, such as recording length, quality of the recorded signal, breathing methods and participant variables (e.g., age, gender, cardio-related medication, …) (Catai et al., 2020;Force, 1996). Related to this, researchers have provided a comprehensive list of practical recommendations (Sassi et al., 2015) to improve the use of standardized protocols for HRV research. All in all, the focus of the current paper is aimed at providing an overview of HRV indices, their usage in psychological research, and a practical demonstration on how to obtain these indices using NeuroKit2.

Author Contributions
DM conceived and TP coordinated the study. TP and ZL participated in the manuscript drafting. DM and AC performed a critical review of the manuscript. All authors read and approved the final manuscript.

Funding
The authors did not receive support from any organization for the submitted work.

Conflict of Interest Statement
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.