Detecting Metabolic Thresholds from Nonlinear Analysis of Heart Rate Time Series: A Review

t rate time series are widely used to characterize physiological states and athletic performance. Among the main indicators of metabolic and physiological states, the detection of metabolic thresholds is an important tool in establishing training protocols in both sport and clinical elds. This paper

of Audiology & ENT, University of Ferrara, 44121 Ferrara, Italy *Correspondence: giovanna.zimatore@uniecampus.it (G.Z.); giuseppe.maulucci@unicatt.it (G.M.) Abstract: Heart rate time series are widely used to characterize physiological states and athletic performance. Among the main indicators of metabolic and physiological states, the detection of metabolic thresholds is an important tool in establishing training protocols in both sport and clinical elds. This paper reviews the most common methods, applied to heart rate (HR) time series, aiming to detect metabolic thresholds. These methodologies have been largely used to assess energy metabolism and to identify the appropriate intensity of physical exercise which can reduce body weight and improve physical tness. Speci cally, we focused on the main nonlinear signal evaluation methods using HR to identify metabolic thresholds with the purpose of identifying a method which can represent a useful tool for the real-time settings of wearable devices in sport activities. While the advantages and disadvantages of each method, and the possible applications, are presented, this review con rms that the nonlinear analysis of HR time series represents a solid, robust and noninvasive approach to assess metabolic thresholds. Keywords: metabolic threshold; heart rate variability; sport; recurrence quanti cation analysis; nonlinear dynamic; Poincar²plot; wearable devices 1. Introduction During incremental intensity exercises, metabolic changes are observed and therefore the identification of metabolic thresholds allows personalized workload prescription milestones to be assessed. According to the Skinner diagram [1], during increasing intensity exercise, three phases of the body's energy usage and two threshold points might be established. With increasing exercise intensity, the aerobic threshold (AerT) [1] describes the transi- tion from aerobic metabolism (production of energy using oxygen) to anaerobic metabolism (production of energy without oxygen); the anaerobic threshold (AnT) [2] describes the transition in which an over-proportional increase in ventilation occurs as a result of in- creased carbon dioxide production. During incremental exercise, these thresholds suggest two different workloads. The rst threshold (AerT) de nes the upper limit of an almost exclusive aerobic exercise with a low blood lactate concentration (about 2 mM/L). The Int. J. Environ. Res. Public Health2022,19, 12719.

ventilation occurs as a result of in- creased carbon dioxide production. During incremental exercise, these thresholds suggest two different workloads. The rst threshold (AerT) de nes the upper limit of an almost exclusive aerobic exercise with a low blood lactate concentration (about 2 mM/L). The Int. J. Environ. Res. Public Health2022,19, 12719.

Int. J. Environ. Res. Public Health2022,19, 12719 2 of 24 threshold AnT is related to the highest workload of maximal lactate production rate in a steady state, i.e., when the production and the removal of blood lactate concentration are in equilibrium (about 4 mM/L). From a practical point of view, the AerT detection is widely used both in medicine and in the sports elds, where the threshold detection is used for investigating sedentary subjects and the active recovery processes. On the contrary, the AnT detection concerns almost exclusively amateur and professional athletes, and its role is to create targeted training programs [3,4]. These thresholds are usually detected by ventilatory (gas exchange) or metabolic (blood lactate) parameters [2]. The gas exchange testing (cardiopulmonary exercise test, CPET) and the blood lactate sampling are the gold standards in threshold detection [5,6]. Threshold estimation occurs by assessing heart and lungs performance at rest and during exercise [7], or evaluating blood lactate concentration [4]. Despite their ef ciency, the widespread use of these methods is severely impaired due to them being highly invasive (as they require blood draws) and impractical (as they require a laboratory equipped with speci c devices and trained operators) [8]. Several papers have demonstrated the validity of analysis of heart rate temporal series during sports activities [9–22] (see Table). This is an important advancement in the eld allowing many of the limitations already reported to be over- come. Indeed, although HR is commonly detected by using a conventional ECG device that requires the use of on-body electrodes wired to a recorder [23], nowadays, a num- ber of wearable smart devices (wrist, watch [24], phone [25], nger ring [26], etc.) are available, allowing the measurement of HR quickly, continuously and in all situations of everyday life including tness activities. These devices optically detect HR trough measurements of volumetric changes of blood in a vein, usually located at the level of the wrist (photoplethysmography—PPG) [27]. The detection accuracy of some of these devices [28,29], that can also contextually detect position, velocity and acceleration, has been certi ed by the Food and Drug Administration

situations of everyday life including tness activities. These devices optically detect HR trough measurements of volumetric changes of blood in a vein, usually located at the level of the wrist (photoplethysmography—PPG) [27]. The detection accuracy of some of these devices [28,29], that can also contextually detect position, velocity and acceleration, has been certi ed by the Food and Drug Administration (FDA). Despite the undoubted advantage in using HR time series, we should note that, typically, data analysis techniques are based on the assumption of signal stationarity. However, the inherently nonstationary nature of HR [30], which physiologically changes continuously to adapt to external stimuli, may generate unreliable conclusions. Although several signal preprocessing methodologies have been proposed to overcome suchissues [9,11,30–34], approaches based on nonlinear analysis appear to be the more common approach. This review is structured as follows: in Section, the main characteristics of HR time series analysis are described. In Section, a comparison is presented among the most representative nonlinear methods, such as the Poincar²geometry [21], the Detrended Fluctuation Analysis (DFA) [16,35], the Entropy [36–38] and the Recurrence Quanti cation Analysis (RQA) [30,33]. In Section, the applications of these techniques in the detection of metabolic thresholds are reported. In Section, the advantages and limitations of these methods are presented. Finally, in the Appendix procedures used in each section. 2. Heart Rate Time Series Analysis Heart physiology is one of the rst areas in Biology where the ideas of complexity and deterministic chaos (familiar to studies in Physics) produced practical applications [31]. As previously mentioned, nonlinear methods of HR classi cation have been developed to highlight and emphasize HR nonlinear uctuations. The oscillation patterns of a healthy heart are complex and constantly changing, allowing the cardiovascular system to rapidly adjust the sudden physical and psychological challenges to homeostasis [30]. The autonomic nervous system (ANS) controls HR through the balancing action of the sympathetic (S) and the parasympathetic (PS) nervous system. Increased S or diminished PS activity results in cardiac acceleration. Conversely, a low S activity or a high PS activity causes cardiac deceleration [23,39]. An additional ne HR regulation is maintained

sudden physical and psychological challenges to homeostasis [30]. The autonomic nervous system (ANS) controls HR through the balancing action of the sympathetic (S) and the parasympathetic (PS) nervous system. Increased S or diminished PS activity results in cardiac acceleration. Conversely, a low S activity or a high PS activity causes cardiac deceleration [23,39]. An additional ne HR regulation is maintained by the respiratory control centers which modulate the vagal

Int. J. Environ. Res. Public Health2022,19, 12719 3 of 24 out ow in the brainstem [10,31,40]. In this context, HR and heart rate variability (HRV) are strictly regulated by ANS and can provide relevant details on its dynamics and control mechanisms. To analyze HR signals, feature-based time-domain methods, frequency- domain methods and non-linear time-domain methods are typically used. 2.1. Feature-Based Time-Domain Methods The methods based on the analysis of the time evolution of an ECG signal or a beat-by-beat HR record (as a discrete time series, FigureInt. J. Environ. Res. Public Health 2022, 19, 12719 3 of 24 Conversely, a low S activity or a high PS activity causes cardiac deceleration [23,39]. An additional fine HR regulation is maintained by the respiratory control centers which modulate the vagal outflow in the brainstem [10,31,40]. In this context, HR and heart rate variability (HRV) are strictly regulated by ANS and can provide relevant details on its dynamics and control mechanisms. To analyze HR signals, feature-based time-domain methods, frequency-domain methods and non-linear time-domain methods are typically used. 2.1. Feature-Based Time-Domain Methods The methods based on the analysis of the time evolution of an ECG signal or a beat- by-beat HR record (as a discrete time series, Figure 1) are termed as time-domain analysis methods. Figure 1. Heart rate variation (in beats per minute—bpm) of a healthy female subject during an incremental exercise (for more details see Appendix A.4). HR record can be expressed in bpm (beats/minute), but the main descriptors of the time-domain analysis are expressed in the literature in terms of the RR interval (milliseconds, ms). RR is the distance between two successive R waves of the QRS signal on the electrocardiogram. Time-domain features are denominated as [23,40]: • RMSSD, the root mean square of difference between adjacent RR intervals • SDNN, the standard deviation of all RR intervals • SDSD or SDDS, the standard deviation of differences between adjacent RR intervals (or RR series) • MSD, the mean successive difference between adjacent RR intervals • MASD, the mean absolute successive difference between adjacent RR intervals 2.2. Frequency-Domain Methods The power spectrum

mean square of difference between adjacent RR intervals • SDNN, the standard deviation of all RR intervals • SDSD or SDDS, the standard deviation of differences between adjacent RR intervals (or RR series) • MSD, the mean successive difference between adjacent RR intervals • MASD, the mean absolute successive difference between adjacent RR intervals 2.2. Frequency-Domain Methods The power spectrum analysis techniques transform time series data into frequency- domain data. The fast Fourier transform (FFT) and the short time Fourier transform (STFT) are the most common power-spectrum analysis techniques [41,42]. These spectral methods require a mathematical manipulation of the data and large time series (i.e., large amounts of data points). Conventional analysis methods based on the Fourier transform technique are not very suitable for the analysis of non-stationary signals, since these signals change over time. The main principle of STFT consists in choosing a small enough window for analysis in which the signal can be considered stationary. Moreover, when slowly and rapidly changing transient events occur, the STFT introduces significant errors. As the analysis window of the time-domain data gets narrower, the time resolution and the assumption of stationarity is improved, but the frequency resolution is reduced. Frequency-domain analysis can distinguish between high frequency (HF) and low frequency (LF) components in the data and can identify peaks related to autonomic control [43] (see Figure 2, the peaks in green and blue). Since the HF and the LF correspond Figure 1. Heart rate variation (in beats per minute—bpm) of a healthy female subject during an incremental exercise (for more details see Appendix). HR record can be expressed in bpm (beats/minute), but the main descriptors of the time-domain analysis are expressed in the literature in terms of the RR interval (millisec- onds, ms). RR is the distance between two successive R waves of the QRS signal on the electrocardiogram. Time-domain features are denominated as [23,40]: RMSSD, the root mean square of difference between adjacent RR intervals SDNN, the standard deviation of all RR intervals SDSD or SDDS, the standard deviation of differences between adjacent RR intervals (or RR series) MSD, the mean

ms). RR is the distance between two successive R waves of the QRS signal on the electrocardiogram. Time-domain features are denominated as [23,40]: RMSSD, the root mean square of difference between adjacent RR intervals SDNN, the standard deviation of all RR intervals SDSD or SDDS, the standard deviation of differences between adjacent RR intervals (or RR series) MSD, the mean successive difference between adjacent RR intervals MASD, the mean absolute successive difference between adjacent RR intervals 2.2. Frequency-Domain Methods The power spectrum analysis techniques transform time series data into frequency- domain data. The fast Fourier transform (FFT) and the short time Fourier transform (STFT) are the most common power-spectrum analysis techniques [41,42]. These spectral methods require a mathematical manipulation of the data and large time series (i.e., large amounts of data points). Conventional analysis methods based on the Fourier transform technique are not very suitable for the analysis of non-stationary signals, since these signals change over time. The main principle of STFT consists in choosing a small enough window for analysis in which the signal can be considered stationary. Moreover, when slowly and rapidly changing transient events occur, the STFT introduces signi cant errors. As the analysis window of the time-domain data gets narrower, the time resolution and the assumption of stationarity is improved, but the frequency resolution is reduced. Frequency-domain analysis can distinguish between high frequency (HF) and low frequency (LF) components in the data and can identify peaks related to autonomic con- trol [43] (see Figure, the peaks in green and blue). Since the HF and the LF correspond to short and long-latency components in the HR time series, they can be, respectively, considered as qualitative indicators of sympathetic and parasympathetic activity [39].

Int. J. Environ. Res. Public Health2022,19, 12719 4 of 24Int. J. Environ. Res. Public Health 2022, 19, 12719 4 of 24 to short and long-latency components in the HR time series, they can be, respectively, considered as qualitative indicators of sympathetic and parasympathetic activity [39]. During heavy exercise (i.e., work intensity above AerT), there is a prevalence of HF components in contrast to LF components. The opposite effect has been observed during moderate exercise, i.e., prevalence of LF components compared to HF [18]. Unfortunately, during long-term recordings, there is a lower HR modulation stability which introduces significant difficulties in the interpretation of the frequency-domain analysis data. For this reason, time-domain methods are primarily used for data analysis. Figure 2. Spectral analysis of heart rate (HR) of a healthy female subject (the same as shown in Figure 1, for more details see Appendix A.4). To overcome the lower HR modulation stability, Cottin et al. (2004) [18] proposed a method where the mean high frequency peak (HF peak, Hz) and the area of HF power (HFp, ms 2 ) could be computed automatically from the spectrogram. Then, the HF product [18,44] can be estimated by multiplying the HF peak (Hz) by HFp (ms 2 ). STFT and a time–frequency domain method were applied in the detection of the met- abolic thresholds, and their applications are discussed in Section 4. 3. Nonlinear Analysis for HR Time Series 3.1. Poincarè Geometry The Poincaré plot is a geometrical and nonlinear analysis method with numerous applications in different scientific domains [45], and it is extensively used for a qualitative visualization of the physiological signals involved [12–14]. It represents one of the most used methods in the assessment of heart rate dynamics in recent literature. In the review by Henriques et al. (2020) [34], it was reported that over 300 publications involving this analysis method on HR time series have been published, resulting in over 18,000 citations. In the Poincaré plot, all values of a time series are plotted against previous values, leading to an ellipsoidal point cloud (Figure 3). Plot shapes can be categorized into func-

review by Henriques et al. (2020) [34], it was reported that over 300 publications involving this analysis method on HR time series have been published, resulting in over 18,000 citations. In the Poincaré plot, all values of a time series are plotted against previous values, leading to an ellipsoidal point cloud (Figure 3). Plot shapes can be categorized into func- tional classes to indicate cardiac functioning [13,15]. Moreover, this diagram may be ana- lyzed quantitatively by the standard deviations of the distances RR–RR(i) to the lines y = x and y = −x + 2(RR–RR m), where RR–RR m is the mean of all RR–RR(i). The standard devi- ations are referred to as SD1 (green line) and SD2 (blue line) (Figure 3). Starting with the relationship between successive beats and the behavior of the heart, the graphic gives summary relevant information: the transverse axis of the ellipse (SD1) is a measure of the short-term changes in the RR intervals and is considered as an indicator of the parasym- pathetic activity; the longitudinal axis (SD2) reflects the long-term changes in the RR in- tervals, and it is considered as an inverse indicator of the sympathetic activity. Figure 2. Spectral analysis of heart rate (HR) of a healthy female subject (the same as shown in Figure, for more details see Appendix). During heavy exercise (i.e., work intensity above AerT), there is a prevalence of HF components in contrast to LF components. The opposite effect has been observed during moderate exercise, i.e., prevalence of LF components compared to HF [18]. Unfortunately, during long-term recordings, there is a lower HR modulation stability which introduces signi cant dif culties in the interpretation of the frequency-domain analysis data. For this reason, time-domain methods are primarily used for data analysis. To overcome the lower HR modulation stability, Cottin et al. (2004) [18] proposed a method where the mean high frequency peak (HF peak, Hz) and the area of HF power (HFp, ms 2 ) could be computed automatically from the spectrogram. Then, the HF product [18,44] can be estimated by multiplying the HF peak(Hz) by HFp(ms 2

for data analysis. To overcome the lower HR modulation stability, Cottin et al. (2004) [18] proposed a method where the mean high frequency peak (HF peak, Hz) and the area of HF power (HFp, ms 2 ) could be computed automatically from the spectrogram. Then, the HF product [18,44] can be estimated by multiplying the HF peak(Hz) by HFp(ms 2 ). STFT and a time–frequency domain method were applied in the detection of the metabolic thresholds, and their applications are discussed in Section. 3. Nonlinear Analysis for HR Time Series 3.1. Poincar±Geometry The Poincar²plot is a geometrical and nonlinear analysis method with numerous applications in different scienti c domains [45], and it is extensively used for a qualitative visualization of the physiological signals involved [12–14]. It represents one of the most used methods in the assessment of heart rate dynamics in recent literature. In the review by Henriques et al. (2020) [34], it was reported that over 300 publications involving this analysis method on HR time series have been published, resulting in over 18,000 citations. In the Poincar²plot, all values of a time series are plotted against previous values, leading to an ellipsoidal point cloud (Figure). Plot shapes can be categorized into functional classes to indicate cardiac functioning [13,15]. Moreover, this diagram may be analyzed quantitatively by the standard deviations of the distancesRR–RR(i)to the lines y = x and y = x + 2(RR–RRm), whereRR–RRmis the mean of allRR–RR(i).The standard deviations are referred to as SD1 (green line) and SD2 (blue line) (Figure). Starting with the relationship between successive beats and the behavior of the heart, the graphic gives summary relevant information: the transverse axis of the ellipse (SD1) is a measure of the short-term changes in the RR intervals and is considered as an indicator of the parasympathetic activity; the longitudinal axis (SD2) re ects the long-term changes in the RR intervals, and it is considered as an inverse indicator of the sympathetic activity.

This review discusses methods for detecting metabolic thresholds using heart rate time series.