Recurrence Quantification Analysis Based Methodology in Automatic Aerobic Threshold Detection: Applicability and Accuracy across Age Groups, Exercise Protocols and Health Conditions

Abstract

A new method based on the Recurrence Quantification Analysis (RQA) of the heart rate (HR) offers an objective, efficient alternative to traditional methods for Aerobic Threshold (AerT) identification that have practical limitations due to the complexity of equipment and interpretation. This study aims to validate the RQA-based method’s applicability across varied demographics, exercise protocols, and health status. Data from 123 cardiopulmonary exercise tests were analyzed, and participants were categorized into four groups: athletes, young athletes, obese individuals, and cardiac patients. Each participant’s AerT was assessed using both traditional ventilatory equivalent methods and the automatic RQA-based method. Ordinary Least Products (OLP) regression analysis revealed strong correlations (r > 0.77) between the RQA-based and traditional methods in both oxygen consumption (VO₂) and HR at the AerT. Mean percentage differences in HR were below 2.5%, and the Technical Error for HR at AerT was under 8%. The study validates the RQA-based method, directly applied to HR time series, as a reliable tool for the automatic detection of the AerT, demonstrating its accuracy across diverse age groups and fitness levels. These findings suggest a versatile, cost-effective, non-invasive, and objective tool for personalized exercise prescription and health risk stratification, thereby fulfilling the study’s goal of broadening the method’s applicability.

1. Introduction

Identifying the aerobic threshold (AerT) is vital for crafting individualized exercise programs and optimizing both health and athletic performance [1,2,3]. The AerT is defined as the exercise intensity at which there is a shift from predominantly aerobic to anaerobic energy production [2]. Traditionally, the AerT has been determined through visually detected and manually uncomfortable methods during physical activities, including physiological variables like gas exchange, heart rate, blood lactate concentration, and rate of perceived exertion, respectively [4,5]. However, these approaches are often time-consuming, and their result is subject to operator-dependent variability during interpretation.

To address these limitations, we recently introduced an objective, efficient, nonlinear method based on Recurrence Quantification Analysis (RQA) that employs heart rate variability (HRV) for automatic AerT detection [6,7,8]. The autonomic nervous system, which controls the body’s physiological reactions to exercise, is reflected in HRV. Consequently, a more thorough and dynamic examination of the physiological reactions to different exercise intensities is provided by utilizing HRV to measure AerT. This connection is particularly significant in clinical settings, where monitoring and management of training, particularly in patients with heart problems or on medications that affect HRV, depend heavily on a thorough understanding of physiological responses to workload. Although HRV-based techniques might not always be more effective, they do have some clear benefits in terms of simplicity of usage and reduced bulky. We hope to simplify the evaluation procedure while preserving the association with recognized physiological measures by employing HRV. This approach allows for a more practical and user-friendly method of AerT evaluation, which could potentially enhance its applicability in various settings. Our previous study tested the HRV-based algorithm on a homogenous sample in terms of age (15.26 ± 1.87 years) and health status (healthy athletes) undergoing the cardiovascular test while performing a single type of activity (cycling) [6].

Indeed, the method’s initial focus was primarily on optimizing athletic performance, which underscores the significance of accurately identifying the AerT not only in sports but also in managing health conditions. Our first study demonstrated a strong correlation (r > 0.77) between the RQA-based method and traditional techniques when applied to a homogenous group of athletes, with a technical error of less than 5% for HR detection at AerT [6]. The precise detection of AerT is critical for assessing the risk of cardiovascular and metabolic disorders in obese individuals [9,10]. For cardiac patients, particularly those on beta-blockers, accurately determining AerT is vital because it is significantly influenced by medications that affect HRV, a crucial indicator of autonomic nervous system function [11,12,13]. Given the foundational importance in both athletic and health contexts, this study aims to expand the utility of our automatic AerT detection method. We are broadening its application to include various physical activities and a more diverse demographic encompassing different ages, genders, fitness levels, and, most importantly, specific health conditions such as obesity and cardiovascular issues treated with beta-blockers. Specifically, we focus on four case studies: athletes, young athletes ages 10–17, obese individuals, and cardiac patients. These groups were selected because they are paradigmatic for evaluating metabolic thresholds.

One of the main advantages of this algorithm is its ability to be integrated into wearable systems for real-time detection, which is becoming increasingly feasible [14]. Flexible bioelectronics, like wearable and epidermal electronics that can be attached to the skin, clothing, or even the body itself, would be ideal for an objective evaluation of physical activity [15]. Furthermore, the rapidly evolving framework of future smart gadgets for tracking physical activity may be built based on small multifunctional sensors [16]. With the help of these new technologies, this automatic algorithm could be integrated into a wearable device, advancing real-time detection. These improvements would make fitness levels easily adjustable and accessible, being beneficial for diseases like obesity [7] and heart problems [8].

2. Materials and Methods

2.1. Participants

This work expands and extends the one presented in [6] by encompassing groups with a broader range of ages, varied training states, and different health conditions and employing a graduated test protocol. We then analyzed 123 cardiopulmonary tests (CPET) from four distinct groups: 27 athletes (Group A [6]), 20 obese individuals (Group O, [7]), 19 cardiac patients (Group C, [8]), and 57 young athletes (Group Y, [17]). The participants (25.35 ± 16.62 years) were between 11 and 70 years old and were divided based on age, fitness level, and health conditions. The athletes and the young athletes’ participants were free from serious health conditions like neuropathy, autonomic dysfunction, and cardiovascular disorders, while the obese individuals and cardiac patients met specific health criteria. The cardiac patients underwent testing before and after rehabilitation sessions. The study was conducted in accordance with the Declaration of Helsinki and received approval from various ethical committees: the CAR-IRB—University of Rome ‘Foro Italico’ Committee (Approval No. CAR 37/2020), the Sapienza University of Rome Ethical Committee (Approval No. 70/11, 2011) for groups A and O, the Ethical Committee of Tours (France) for group C, and the Institutional Review Board at Appalachian State University (IRB no. 18-0147) for group Y.

Participants provided informed consent, and clinical examinations were conducted to ensure safety and eligibility for physical activity. The combined analysis offers a comprehensive view of the physiological responses across a wide age range and varying health and fitness levels, employing time series data collected during incremental exercise tests. The methodologies and protocols were carefully designed to ensure data reliability and ethical compliance, providing valuable insights into cardiopulmonary responses under different physical conditions.

Groups underwent different cardiopulmonary tests (CPET), which are described in Appendix A. During these tests, time-series data related to gas exchange (VO₂, VCO₂, oxygen saturation SpO₂, etc.) and heart rate (HR) were collected. In this study, we have focused exclusively on analyzing the HR and VO₂ time series, comparing them against the gold standard for validation purposes.

2.2. Exclusion Criteria

HR time series were excluded from the analysis under the following criteria:

• (a). Insufficient consecutive HR data points (fewer than 150, corresponding to a test duration of 5 min)

• (b). Excessive difference (>10 beats per minute) between consecutive HR data points

• (c). HR values lower than the average of the preceding 10 data points within the time series

• (d). The presence of outliers in the HR time series is defined as values exceeding 4 standard deviations from the mean of that series.

2.3. Preprocessing

The first and the last steps of the graded exercise were discarded. In the context of a graded exercise test, whether to discard the first and last steps of the exercise for the analysis of the aerobic threshold depends on various factors, including the specific protocol of the test, the purpose of the analysis, and the characteristics of the individual being tested. Graded exercise typically involves a progressive increase in exercise intensity, and it is used to assess cardiovascular fitness and determine aerobic and anaerobic thresholds. Here are some considerations regarding the first and last steps.

The initial stage of the graded exercise is often a warm-up phase where the individual is just beginning to exercise, and the body is transitioning from a state of rest to one of physical activity. Moreover, physiological responses (like HR and oxygen consumption) might not have stabilized yet. This could lead to less reliable data for determining the aerobic threshold.

The final stage of the graded exercise usually represents the point of maximal or near-maximal effort. Here, the individual is often in a state of anaerobic metabolism. For aerobic threshold analysis, this phase may not be as relevant because the focus is on the point where the body transitions from aerobic to anaerobic metabolism, which typically occurs at a lower intensity level. Data from the last step can be more relevant for determining the anaerobic threshold or maximal aerobic capacity, which is one of our study’s possible developments. In this study, which aimed to generalize the algorithm in [6] between heterogeneous groups and protocols, we decided to discard the first and last steps for every group and protocol (see Figure 1 in [6]).

2.4. Data Cleaning

To cleanse the time series data of Oxygen Consumption (VO₂) and Heart Rate from outliers resulting from malfunctions in the CPET (Cardiopulmonary Exercise Testing) equipment and the Garmin heart rate monitor while minimally altering the raw data, we employed a smoothing technique based on a moving average, which is particularly useful in smoothing out short-term fluctuations and highlighting longer-term trends or cycles (see in details Appendix B). An example of a moving average applied to an HR time series is reported in Figure 1.

2.5. Gas Exchange Method (Gex Method) for AerT Detection

The gold standard method for the identification of critical physiological thresholds, such as the AerT, is the Gas Exchange Method (Gex method), which requires a Cardiopulmonary Exercise Test (CPET), an advanced and comprehensive test that measures pulmonary and cardiovascular response during exercise [18,19]. During a CPET, individuals typically perform graded exercises on devices such as those manufactured by Cosmed^®. These devices are engineered to meticulously monitor and record various physiological parameters. The experimental setup diagram and a detailed description of the test system are reported in Figure 2.

The Gex method specifically examines the relationship between the ventilatory equivalent for oxygen (VE/VO₂) and oxygen consumption (VO₂). The ventilatory equivalent for oxygen is a dimensionless number that represents the ratio of the volume of air ventilated to the amount of oxygen consumed. In practice, the method involves graphing VE/ VO₂ against VO₂ values to identify the AerT. The Aerobic Threshold is discerned (by visual inspection) at the point where the VE/VO₂ ratio reaches its minimum value before it begins to rise again during a progressively strenuous exercise protocol.

The Gex procedure for determining AerT was validated using a second method known as V-slope [2]. This method involves plotting a graph with VCO₂ on the y-axis and VO₂ on the x-axis. Two regression lines are then fitted to represent the upper and lower segments of this relationship. The intersection of these lines corresponds to the Aerobic Threshold determined by the Gex procedure.

To accurately determine the HR and VO₂ at the Aerobic Threshold, a precise method is employed. The average values of HR and VO₂ are calculated over the final 30 s of the exercise stage just before the AerT is reached. This approach is taken because these last 30 s are seen as reflective of the individual’s physiological state at the threshold. Capturing the average values during this critical period is thought to provide a stable and reliable set of data points. These data points are crucial for subsequent statistical analysis, offering a solid foundation for comparing various conditions or interventions that might affect the Aerobic Threshold.

2.6. RQA-Based Method for AerT Automatic Detection

RQA is a powerful method used to analyze nonlinear dynamical systems. It involves reconstructing time series data into a phase space using delay embedding techniques, where vectors are created based on an embedding dimension m and time delay τ. A matrix R is then formed by comparing points in this phase space using a predefined threshold distance ϵ. The resulting matrix is graphically represented as a Recurrence Plot (RP), showing patterns of recurrence, such as diagonal lines indicating repeating behaviors. Key parameters like ϵ, τ, and m are chosen carefully to suit the characteristics of the data and study goals. Measures derived from the recurrence plot, such as Determinism (DET) and Laminarity (LAM), quantify the predictability and structure of the system based on the lengths of diagonal and vertical (or horizontal) lines in the plot, respectively. For a more detailed overview of this nonlinear analysis, see Appendix C.

RQA has been used to provide a straightforward framework for describing the coupling of cardio-respiratory physiology under increasing physical exertion in terms of both complexity and dynamical transitions. This RQA-based methodology has been validated using DET and LAM to identify ventilatory thresholds, as outlined in references [6,7].

In particular, the assessment of the AerT is based on identifying a minimum in DET (DET_min), which is corroborated by a corresponding minimum in LAM (LAM_min).

Automatic Detection Algorithm

The RQA-based method has its limitations. Indeed, these thresholds have been identified so far by visually inspecting where DET/LAM underwent a sudden change or reached saturation. The need now is to determine these thresholds without the requirement of external user intervention but rather in an automated manner. To achieve this result, we developed a novel algorithm written in Python (version 5.3.3) that can automatically detect the most convex minima of a time series. This algorithm has been already presented in [6]. For more details, see Appendix D.

Here we report the main steps:

1. Data cleaning and preprocessing on HR and VO₂ time series.

2. The HR time series is divided into partially overlapping epochs of width = 100 points, shifted by 1 point, and for every epoch, RQA is performed, and DET values are calculated.

3. The new DET time series is then analyzed to identify the most convex minima. After calculating the first and second derivatives, a point is considered a minimum if its second derivative is greater than or equal to the mean of the second derivative plus two times its standard deviation:

(1)$f″\left(t\right)>mean\left(f″\left(t\right)\right)+2\times SD\left(f″\left(t\right)\right)$

The value $mean\left(f″\left(t\right)\right)+2\times SD\left(f″\left(t\right)\right)$ is named cut-off or radius.

This condition selects points where the curvature is significantly convex. Finally, the function returns a list of indices (minima) corresponding to these minima points.

• 4.. The AerT is identified by locating the most convex minimum of DET using the formula:

(2) ${\mathrm{DET}}_{\min }=\max \left({\displaystyle \frac{d^2\mathrm{DET}\left(t\right)}{d{t}^2}}\right)$

This selection is among all relative minima of DET. The point at which AerT occurs is termed DETmin, identified at the time point tVT1.

2.7. Statistics

In this context, variables are named following a specific convention to denote the method used for their calculation. Variables appended with “_RQA” are those calculated using the RQA-based method. This suffix indicates that the analysis and values of these variables are derived from RQA techniques. On the other hand, variables with the suffix “_Gex” are calculated using the Gas Exchange (Gex) method. This naming convention helps in differentiating the source and methodological background of each variable, ensuring clarity and precision in our analysis and discussions.

The Kolmogorov–Smirnov Test of Normality [20] confirmed the normal distribution of VO₂ and HR values at the AerT. The outcomes were as follows: VO₂_RQA (D = 0.09, p = 0.26) and HR_RQA (D = 0.07, p = 0.64).

Repeated measures analysis of variance (RM ANOVA) was utilized to assess the main effects of methods (RQA vs. Gex) and the interaction effects.

Agreement between the RQA-based and Gex methods was assessed using an Ordinary Least Products (OLP) regression analysis [21].

The regression’s coefficients of determination (R²), regression parameters (slope and intercept), and 95% confidence intervals (CI) were calculated to identify fixed and proportional biases. The following criteria were adopted to interpret the magnitude of the correlation r between the estimates: <0.1 trivial, 0.1–0.3 small, 0.3–0.5 moderate, 0.5–0.7 large, 0.7–0.9 very large, and 0.9–1.0 almost perfect [20]. The hypothesis of proportional and fixed bias was rejected when the 95% CI for the slope included 1 and the intercept included 0.

Percentage differences between the RQA-based and Gex methods at AerT for HR and VO₂ were calculated using Typical Percentage Error (TE(%)) and Mean Percentage Difference (Meandiff (%)).

Intraclass Correlation Coefficients (ICC) were employed as measures of criterion validity for the RQA-based method compared to the Gex method at AerT [22].

The magnitude of differences between the RQA-based and Gex methods was assessed using Cohen’s d effect size statistics. Effect sizes were classified as follows: d < 0.2 was considered trivial, 0.2 < d < 0.6 as small, 0.6 < d < 1.2 as moderate, 1.2 < d < 2.0 as large, and d > 2.0 as very large. Additionally, the interpretation of ICC ranges was defined as excellent (ICC > 0.9), moderate (0.75 < ICC < 0.90), fair (0.50 < ICC < 0.75), and poor (ICC < 0.5) [23].

Bland–Altman plot analyses were employed to visualize the agreement between RQA-based and Gex methods [24].

Statistical significance was set at p < 0.05 in both studies, and all statistical analyses were conducted using SPSS version 24.0 software (SPSS Inc., Chicago, IL, USA).

3. Results

3.1. Variability of Heart Rate and Oxygen Consumption at Aerobic Threshold across Participant Groups

The final number of time series considered in this study was 111, since 12 time series did not meet the inclusion criteria introduced in the Materials and Methods section. Table 1 serves as a comprehensive overview of participant profiles and study parameters, allowing for a detailed comparison and analysis across different demographic and health-related categories within the study population.

Figure 3 presents HR and VO₂ values at the AerT for the study participants. It showcases a wide spectrum of HR values, ranging from 60 to 180 beats per minute (bpm), and VO₂ values, spanning from 600 to 2500 mL/min at AerT. Notably, young athletes typically display lower VO₂ values, which can be attributed to their smaller body sizes. Conversely, cardiac patients often exhibit both lower VO₂ and HR values, reflecting their distinct physiological condition.

Figure 4 illustrates two representative HR time series and outlines the key steps of the automatic detection process. These time series were analyzed using the Python script with our automatic detection algorithm. The series on the right pertains to a young athlete from group Y, while the series on the left is from a cardiac patient in group C. The blue graded lines vividly depict the incremental test’s workload. As the intensity of physical exercise increases, so does the HR, as observed in panels (a) and (d). Panel (d) demonstrates how beta blockers effectively regulate the increase in HR. Through the RQA of HR, the determinism for each epoch and its second derivative were computed, as shown in panels (b) and (e). The colored points mark relative minima identified by our automatic algorithm. In panels (c) and (f), the most convex minimum, which aligns with the AerT value, is highlighted by a red vertical line.

3.2. Comparative Analysis of Traditional and RQA-Based Methods in Assessing Aerobic Threshold

The concordance between the traditional Gex method and the RQA-based automated technique for determining the AerT has been examined for both HR and VO₂. An analysis of variance with repeated measures (ANOVA RM) showed no significant differences between the two methods for HR and VO₂ at AerT (p > 0.05), as outlined in Table 2.

Additionally, the ANOVA RM revealed no significant differences in HR and VO₂ across various groups when considering factors such as weight status (72 Normal or 39 Overweight), health status (76 Healthy from groups Y and A; 35 Diseased from groups C and O), and gender (73 Male or 38 Female).

The substantial difference in participant numbers between the treadmill (19 time series) and cycle-ergometer (92 time series) protocols could influence the statistical outcomes and significance observed in the comparison methodXprotocol for HR values (p = 0.03). A larger sample size in the cycle-ergometer group likely enhances statistical power, increasing the ability to detect significant differences between the automated detection methods. Moreover, larger samples generally yield more precise estimates and may better represent the target population. However, it is important to acknowledge potential biases and variations in participant characteristics across these protocols, which could impact the generalizability and interpretation of results.

Moreover, the comparison methodXgroup for HR values (p = 0.04) is influenced by the difference between Group Y (coupled data t-test, p < 0.05) and the other groups (coupled data t test p > 0.1 for A, C, and O group). Young athletes typically exhibit distinct physiological characteristics compared to adults, including higher heart rates at rest and during exercise due to factors such as increased cardiac output, lower resting heart rates, and greater cardiovascular efficiency. These differences can influence HR measurements and may contribute to variations observed in the comparison of HR metrics between groups.

Table 3 reports various parameters and statistics related to the comparison between two methods (RQA and Gex) for measuring HR and VO₂ in the total sample of 111 time series. The R-squared (R2) values (0.68 for HR and 0.6 for VO₂) suggest a moderate to good level of agreement between the two measurement methods. The slope indicates how measurements from one method relate to measurements from the other method. For both HR and VO₂ the slope value is around 1, which indicates perfect agreement. The intercept provides insight into any systematic bias or offset between the two methods. The values for HR and VO₂ are within the accepted range, suggesting that there is no bias. Meandiff(%) and TE(%) reflect the magnitude and variability of differences between measurements obtained by the two methods.

The mean percentage difference between the two assessment methods was below 2.5% for HR and under 8% for VO₂, with these discrepancies not proving statistically significant (p > 0.05) at the Aerobic Threshold. TE was considered acceptable, with values below 8% for HR and around 18% for VO₂, indicating a high level of precision in the measurements [24,25].

The Intraclass Correlation Coefficients (ICC) demonstrated excellent reliability, with scores of 0.90 for HR and 0.87 for VO₂, indicating strong agreement between methods [22]. The effect size, calculated using Cohen’s d, was found to be trivial for both HR and VO₂, indicating minimal differences between the methods, as detailed in Table 4.

The Ordinary Least Products (OLP) regression analysis plots for HR and VO₂ at the AerT between both methods are visually depicted in Figure 5. Here, the OLP line of best fit is drawn as a solid line, while the line of unity is depicted as a dotted line, facilitating a clear comparison between the traditional and automated techniques.

Furthermore, Bland-Altman plots, illustrated in Figure 6, showed no systematic bias for either HR or VO₂. The plots indicated that the dispersion of differences was more contained for HR than VO₂, highlighting VO_2′s greater variability. Only a few outliers were present, but they did not significantly affect the overall concordance between the methods, thereby reinforcing the robustness of the findings.

4. Discussion

The most recent published study introduced an automated method for identifying the aerobic threshold using Heart Rate (HR) data gathered during the Cardiopulmonary Exercise Test (CPET), specifically in adult athletes, including rowers and amateurs [6]. This current research aimed to expand the application of this innovative method by validating its reliability across a broader sample. Specifically, we tested the method’s effectiveness across different demographic groups, including variations in age, gender, and fitness levels, as well as among individuals with diverse health conditions, such as obesity and cardiovascular diseases, and those undergoing treatment with beta-blockers. By examining these variations, we aimed to establish the method’s robustness and adaptability in detecting the AerT under a wide range of physiological conditions, ensuring its relevance for both general health and specialized athletic training. Our study extended the automated detection method to 111 time series from individuals between 11 and 70 years old and including various health conditions (healthy children, obese individuals, and cardiac patients). The tests were administered on both treadmills and cycle-ergometers, necessitating different incremental protocols to achieve a consistent total exercise time and thereby yielding comparable HR time series [6,7,8].

The post-hoc power of the F-test increased from 0.88 for a 27-subject sample [6] to 0.999 for our expanded sample of 111 subjects. In the present research, the HR values derived from the automated RQA-based method were not significantly different from those obtained using the gold standard ventilatory equivalents-based method (Gex), as evidenced by a p-value greater than 0.05. Furthermore, a strong correlation (r > 82%) and an acceptable Technical Error (TE) of 8% were observed. Importantly, despite the potential effects of beta-blockers in cardiac patients altering the HR variability during exercise, the RQA-based method proved to be reliable for this group. However, it’s worth noting that in individuals with heart conditions, the use of heart rate as an exercise threshold marker is compromised by medication effects. For a more accurate exercise prescription in these cases, reliance on the Rating of Perceived Exertion (RPE) is advised [4]. General physiological considerations could be confirmed through intergroup comparisons. For instance, regular exercise enhances cardiovascular efficiency, reflected in higher HR and maximum oxygen consumption (VO_2max) values among athletes and young athletes when compared to obese or cardiac patients [1,3,8]. Age-related factors like lung volume also exhibited expected trends.

The RQA-based method shows promising future applications in the detection of AerT, a fundamental and complex challenge of exercise physiology [26]. By analyzing the recurrence patterns in physiological signals, RQA could provide deeper insights into the complex interplay of physiological processes and the critical thresholds demarcating different metabolic states. This enables personalized threshold identification, capturing unique patterns specific to each individual test. For the young athletes’ group, it’s crucial to consider the physiological uniqueness of adolescence, including growth spurts and changes in aerobic capacity. Tailored reference values and assessment methods are needed for this population [27]. Additional applications include identifying nonlinear markers related to fatigue and performance degradation, aiding in training optimization and performance prediction. As wearable technologies advance, the integration of RQA-based methods for real-time metabolic monitoring is increasingly feasible [28]. Machine learning algorithms could further refine threshold detection, particularly when combined with wearable device data or during tele-exercise sessions. Such methods would make fitness levels readily accessible and customizable, which are also useful for conditions like obesity [7] and cardiac issues [8]. Advanced machine learning techniques can efficiently explore functional network variations based on extensive data [28]. Recently, Kaufmann et al. 2023 compared established lactate (LT) and ventilatory thresholds (VT) with HR-derived thresholds for exercise intensity prescription in endurance athletes [29]. Their systematic review concluded that HR-derived thresholds offer a promising alternative, for example, can be found in the results of Schaffarczyk et al. 2023; they have shown good agreement between the heart rate and oxygen uptake at first and second ventilatory threshold by a nonlinear analysis named detrended fluctuation analysis (DFA) in a cohort of randomly selected women [30].

Advantages and Limitations

The main advantages of the RQA method are simplicity, non-invasiveness, and automation, which reduces subjectivity. The method is less expensive and more practical compared to traditional methods that require complex equipment and human intervention for interpretation. The issue of time consumption is likewise addressed by the automatic method. After 30 to 40 min (including the preparation, exercise, and recovery phases) of the CPET test, the time needed to determine metabolic thresholds drops from 10s of minutes for an operator to about 3 s. This is the amount of time needed for the computational time of the Python program, which includes all stages (preprocessing, data cleaning, and threshold calculation). This proves that it is more than possible to include this automatic system for real-time detection in wearable technology.

Limitations of the study include the number of subjects: a larger and more diverse sample could provide a more comprehensive understanding of the methodology’s applicability and accuracy. Group C is limited to a specific cardiovascular pathology; further investigations are needed to verify the applicability of the methodology to individuals with different health challenges. The within-subject variability could be studied if multiple tests on the same subject were conducted. The sample consisted of individuals from a predominantly European background. We acknowledge the importance of ethnic diversity and plan to conduct further studies on more heterogeneous populations to investigate any potential differences and further validate our method.

5. Conclusions

In summary, nonlinear heart rate variability analysis (RQA) during graded exercises holds tremendous promise for future research and practical applications in the identification of metabolic thresholds based on HR. This method can facilitate personalized training plans, performance monitoring, and real-time adjustments through wearable technology integration. Further research is likely to deepen our understanding of exercise physiology, thereby advancing the fields of sports science, exercise physiology, and healthcare.

Author Contributions

Conceptualization, G.Z., C.S., M.C.G. and L.G.; methodology, M.M., D.F., G.Z. and C.S.; software, G.Z. and C.S.; validation, G.M. and M.D.S.; data curation, D.F., L.M. and M.M.; writing—original draft preparation, G.Z., C.S. and M.C.G.; writing—review and editing, G.M., M.D.S. and C.B.; visualization, D.F.; supervision, C.B. and L.G. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

The study was conducted in accordance with the Declaration of Helsinki and received approval from various ethical committees: the CAR-IRB—University of Rome ‘Foro Italico’ Committee (Approval no. CAR 37/2020), the Sapienza University of Rome Ethical Committee (Approval No. 70/11, 2011) for groups A and O, the Ethical Committee of Tours (France) for group C, and the Institutional Review Board at Appalachian State University (IRB no. 18-0147) for group Y.

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Data Availability Statement

The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. An exemplificative graph depicting the original heart rate time series for a graded exercise over 500 points, along with the smoothed version using a 5-point moving average. The red line represents the smoothed heart rate.

Figure 2. Standard CPET equipment: (a) mask, volume sensor, and gas analyzer tubing; (b) blood pressure monitor; (c) smartwatch for HR data collection; (d) ergometer (cycle or treadmill); (e) gas analyzer; and (f) display.

Figure 3. HR (bpm) vs. VO2 (mL/min) at AerT from 111 CPET tests. The participants are divided into athletes (group A, blue cross), obese (group O, red circle), cardiac patients (group C, light-blue diamond), and young athletes (group Y, yellow circle), respectively.

Figure 4. Representative HR time series, panels (a,d), belong to groups C and Y, respectively. Panels (b,e) show DET(%), and the second derivative where the dashed red horizontal line corresponds to the cut-off; in panels (c,f), DET (%) (in black), and workload (in blue) are shown point by point, respectively. The red vertical line corresponds to AerT (the most convex minimum of DET (%), as explained in the Materials and Methods section).

Figure 5. OLP regression line on (a) HR and (b) VO2 at AerT (VT1), all values correspond to the minima DETmin. The participants are divided into athletes (group A, blue cross), obese (group O, red circle), cardiac patients (group C, light-blue diamond), and young athletes (group Y, yellow circle).

Figure 6. Bland Altman plot for (a) HR (bpm) and (b) VO2 (mL/min) at AerT. The participants are divided into athletes (group A, blue cross), obese (group O, red circle), cardiac patients (group C, light-blue diamond), and young athletes (group Y, yellow circle), respectively. The horizontal black line corresponds to bias (mean difference), and the dashed horizontal black line to the lower and upper limits of agreement (LOA); LOAs are calculated as the mean difference ±1.96 standard deviations (SD).

Table 1

Overview of participant profiles and study parameters.

	Group A	Group O	Group C	Group Y	All
Name	Athletes	Obese	Cardiac Patients	Young athletes	Athlete, obese
					Young athletes
					cardiac patients
n	27	19	16	49	111
Age(years)	15.30 ± 1.87	42.90 ± 12.67	52.80 ± 10.09	15.20 ± 2.02	25.35 ±16.62
BMI (kg/m2)	22.06 ± 3.28	39.98 ± 4.95	25.92 ± 4.29	22.65 ± 4.97	25.78 ± 7.60
Range (years)	14–20	19–63	35–70	11–17	11–70
Gender	25/2	6/13	16/0	26/23	73/38
(M/F)
Wstatus (NW/OW)	27/0	0/19	13/3	32/17	72 NW/39 OW
Hstatus	Healthy (H)	Disease (D)	Disease (D)	Healthy (H)	76 H/35 D
Incremental	20 or	1 km/h every	10 (Watts/min)	10 (Watts/min)	10, 15, 20
protocol	15(Watts/min)	3 min			(Watts/min)
					1 km/h /3 min
Type	Cycle-ergometer (C)	Treadmill (T)	Cycle-ergometer (C)	Cycle-ergometer (C)	90 C/19 T

Wstatus is the weight status (NW = normal weight, OW = overweight); Hstatus is health status.

Table 2

Main Parameters at AerT by two methods (mean ± SD) in athletes (A), cardiac patients (C), obese (O), and young athletes (Y).

Group	HR_RQA (bpm)	HR_Gex (bpm)	VO2_RQA (mL/Min)	VO2_Gex (mL/Min)
A (27)	136.81 ± 15.46	137.34 ± 15.26	1776.56 ± 445.95	1767.57 ± 335.89
O (19)	110.22 ± 14.06	112.37 ± 14.09	1379.28 ± 335.19	1359.12 ± 316.4
C (16)	87.31 ± 12.53	84.69 ± 13.01	1051.09 ± 462.49	902.79 ± 242.67
Y (49)	126.36 ± 14.13	121.16 ± 14.18	1134.54 ± 343.67	1056.57 ± 260.37
All (111)	120.5 ± 21.36	118.62 ± 21.44	1320.57 ± 472.98	1263.70 ± 427.71
p method *	0.50		0.10
p method X group	0.04		0.38
p method X protocol	0.03		0.33
p method X Wstatus *	0.97		0.62
p method X Hstatus *	0.87		0.86
p method X gender *	0.98		0.61

Two factors (Gex vs. RQA) repeated measures for HR and VO₂ values across all participants, groups (A, C, O, Y), protocol (Cycle-ergometer or treadmill), Wstatus (obese and normal weight), Hstatus (healthy and disease) and gender (male, female). * Gex method vs. RQA-based method all p > 0.05; post-hoc power = 0.998, n = 111 subjects.

Table 3

Agreement with the gold standard (RQA-based vs. Gex method) for HR (bpm) and VO₂ (mL/min), Pearson correlation (R²), slope and intercept values, mean difference (Mean_diff(%)), t-test (p) and Typical Percentage Error (TE(%)) of the whole sample.

Parameters	R2	Slope	Intercept	Meandiff (%)	p	TE (%)
HR	0.68	1.0 (0.57 to 1.75)	2.31 (−86.77 to 53.11)	2.23	0.121	7.71
VO2	0.60	1.11 (0.72 to 1.69)	−76.86 (−815.94 to 406.55)	6.51	0.053	17.54

Table 4

Agreement with Gex method for HR (bpm) and VO₂ (mL/min) value.

Parameters	r Pearson	ICC	Effect Size (d)
HR (bpm)	0.82 **	0.90 (0.86 to 0.93)	0.10 *
VO2 (mL/min)	0.77 **	0.87 (0.81 to 0.91)	0.14 *

Pearson correlation coefficient (r Pearson), Intraclass Correlation Coefficient (ICC), and Cohen’s Effect Size (d); ** all large; * d < 0.2 trivial; (0.2 < d < 0.6 moderate); ICC > 0.9 excellent, 0.75 < ICC < 0.90 moderate.

Appendix A

Cardiopulmonary Tests

Groups underwent different cardiopulmonary tests (CPET), which are summarized below:

Group A (Athletes): In the morning, under temperate conditions (21–22 °C, 50–60% humidity), rowers with varied training undergo a cycling exercise test. Following a standard breakfast and clinical assessments, they start cycling at 60–70 rpm. The test begins with a minute of rest, then a minute at 0 W, escalating the workload by 20 W/min for one group and 15 W/min for the others. Their exertion, respiratory exchange, heart rate, and gas exchange are closely monitored using the Quark RMR-CPET Cosmed™ (Rome, Italy) and a Garmin^® (Olathe, KS, USA) heart rate monitor, ensuring a detailed evaluation of their physical response [6].

Group O (Obese): During similar morning conditions, obese participants engage in a treadmill exercise following their usual breakfast and anthropometric measurements. The exercise begins at a leisurely 3 km/h, gradually increasing in speed and slope. Their cardiovascular and respiratory responses are meticulously recorded, utilizing a Woodway PRO (Woodway, Waukesha, WI, USA) treadmill and Quark RMR-CPET Cosmed™ equipment [7].

Group C (Cardiac patients): These individuals underwent their exercise regimen in the morning. Using a cycle-ergometer, the exercise starts at a low intensity of 20 watts, increasing by 10 watts every minute until the participant reaches exhaustion, indicated by a plateau in oxygen uptake. This protocol is especially cautious, considering their heart condition, and includes comprehensive monitoring of heart rate and gas exchange using specialized equipment such as the ERG 900 and GE Medical System’s CASE Exercise Testing System (GE HealthCare Technologies, Inc., Chicago, IL, USA) [31,32].

Group Y (Young Athletes): Conducted either late in the morning or late afternoon, young athletes undergo their tests in controlled conditions (20–22 °C, 40–50% humidity, 1015 mb above sea level). After ensuring no physical activity or stimulant intake prior to the test, they engage in a cycling exercise that starts with a resting phase, followed by unloaded pedaling, and then an incremental increase of 15 watts per minute until exhaustion. Their physical responses, including heart rate and respiratory parameters, are closely monitored using advanced equipment like the LODE Corival cycle-ergometer (Groningen, The Netherlands) and COSMED K5 Portable Metabolic Analysis Technology System (Rome, Italy), providing vital data on their endurance and physical capabilities [17].

During each test, the perception of physical exertion was assessed using the OMNI Scale of Perceived Exertion (ranging from 0 to 10) [4], measured 15 s prior to the end of each stage. Participants were asked to rate their subjective intensity of effort on a scale from 0 (extremely easy) to 10 (extremely hard). The test was terminated when any of the following conditions were met: a rating of 10 on the OMNI Scale, 90% of the subject’s predicted Maximum Heart Rate (HRmax, calculated using the well-known Fox formula HRmax = 220-age, beats/min [18]), or a respiratory exchange ratio (the ratio between the volume of CO₂ being produced by the body and the amount of O₂ being consumed) equal to 1.1.

Appendix B

Data Cleaning

To cleanse the time series data of Oxygen Consumption (VO₂) and Heart Rate from outliers resulting from malfunctions in the CPET (Cardiopulmonary Exercise Testing) equipment and the Garmin heart rate monitor, while minimally altering the raw data, we employed a smoothing technique based on a moving average, which is particularly useful in smoothing out short-term fluctuations and highlighting longer-term trends or cycles. The moving average is calculated over a specified number of data points, known as the window size. For example, a 5-point rolling mean would calculate the average of data points 1 to 5, then 2 to 6, 3 to 7, and so on. In each step, the mean (average) of the values within the window is computed. This mean then represents the value of the middle data point of the window. The window is then moved one data point forward, and the process is repeated. This continues until the window has moved across the entire data set. An example of moving average applied to a heart rate time series is reported in the manuscript in Figure 1.

By averaging the data, the rolling mean smooths out short-term fluctuations and allows for the identification of longer-term trends. The choice of window size greatly affects the degree of smoothing: a larger window will produce a smoother line but can obscure more detailed fluctuations. For our time series we selected a window width of 5 data points.

Appendix C

Recurrence Quantification Analysis

RQA is a statistical, graphical, and analytical tool [33,34] that has been successfully used in several fields, such as finance [35], physiology [36,37,38,39,40], earthquakes, and geoscience [41,42]. It is useful for analysing nonlinear dynamic systems. It has been demonstrated that RQA is a reliable tool for analysing the stochastic and chaotic dynamics of physiological signals [43,44].

The steps for the RQA analysis are explained below:

• Phase Space Reconstruction: the time series data is often reconstructed into a phase space to reveal the underlying dynamics. This is done using delay embedding techniques. Let’s supposed to have a time series x(t), where $t=1,2,\dots ,N$. The phase space can be reconstructed by creating a set of vectors Y(i), where each vector is given by:

(A1)$$Y\left(i\right)=\left[x\left(i\right),x\left(i+\tau \right),x\left(i+2\tau \right),\dots ,x\left(i+\left(m-1\right)\tau \right)\right]$$

Here, m is the embedding dimension, and τ is the time delay, both of which need to be appropriately chosen.

• Constructing the Recurrence Matrix: A matrix R is then constructed by comparing each point in the phase space with every other point. The elements of the matrix R are defined by:

(A2)$$Ri,j=\Theta \left(\epsilon -\mid \mid Y\left(i\right)-Y\left(j\right)\mid \mid \right)$$ where Θ is the Heaviside step function, ε is a predefined threshold distance (or radius), and ∣∣⋅∣∣ denotes a norm (often the Euclidean norm). Ri,j has value 1 if the distance between Y(i) and Y(j) is less than ε; otherwise, its value is 0.

The recurrence matrix is then obtained:(A3)$$\left[\begin{array}{c}\begin{array}{ccc}R{R}_1& R{R}_3& \dots \\ R{R}_2& R{R}_4& \dots \\ \dots & \dots & \dots \\ R{R}_N& \dots & \dots \end{array}\end{array}\right]=\left[\begin{array}{cccc}{e}_{1,1}& e_{1,2}& \dots & e_{1,m}\\ e_{2,1}& e_{2,2}& \dots & e_{2,m}\\ \dots & \dots & \dots & \dots \\ e_{N-m+1,1}& e_{N-m+1,2}& \dots & e_{N-m+1,m}\end{array}\right]$$

• Visualizing the Recurrence Plot: The recurrence plot is a two-dimensional representation of the matrix R. Each point (i,j) is plotted as a small dot if Ri,j = 1. The result is a plot that reveals patterns of recurrence within the dynamical system. Diagonal lines in the plot, for example, indicate a repeating pattern or periodic behaviour (Figure A1).

Figure A1. Unthresholded recurrence plot of a HR time series recorded as a breath-by-breath, from a cardiopulmonary exercise test (CPET) device (Cosmed, Rome, Italy).

The choice of parameters like ε, τ, and m is crucial, as it can significantly affect the structure and interpretation of the RP. Typically, these parameters are chosen based on the specific characteristics of the time series being analyzed and the objectives of the study. For our study, the following RQA parameters were used and held constant: line = 4, embedding dimension = 7, radius or cut-off = 5. For further details about parameters optimization procedure, see [6].

RQA involves quantifying the patterns observed in the recurrence plot. Several measures are used to describe the patterns, including:

• Recurrence Rate (REC): This is the simplest RQA measure and is defined as the fraction of points in the recurrence plot that are recurrences.

• Determinism (DET): This measure reflects the predictability of the system. It is the fraction of recurrence points that form diagonal lines of at least a minimum length l_min in the RP.

• Laminarity (LAM): This measure reflects the proportion of recurrence points that form vertical (or horizontal, depending on the orientation) lines of at least a minimum length v_min in the RP.

The choice of the minimum line length is crucial and depends on the specific application or the characteristics of the data being analyzed.

These measures provide a quantitative way to describe the patterns observed in a recurrence plot, offering insights into the dynamics of the system under study. They are particularly useful in identifying and characterizing complex, nonlinear behaviours in time series data.

RQA is used in a wide range of fields including physics, biology, finance, and environmental sciences. It is particularly useful in analysing complex, nonlinear, or non-stationary time series data. One of the main advantages of RQA is its ability to analyse short and non-stationary data. It does not require large datasets to identify patterns and can be used to detect changes or transitions in dynamical systems. Interpretation of RQA requires expertise, as the patterns in a recurrence plot can be complex. The choice of parameters, like embedding dimension and delay, which are required to construct the recurrence plot, can significantly influence the results and must be chosen carefully.

Appendix D

Python Code for Automatic Detection Algorithm

The main core of the Automatic Detection algorithm lies in the “find_convex_minima” procedure. The function first calculates the first derivative (dev1) of the DET time series and then it calculates the second derivative (dev2). This is the rate of change of the first derivative, giving a sense of the curvature of the data. The function then computes the mean (mean_val) and standard deviation (devst) of the second derivative. The function then iterates through the data frame to identify points that are considered the most convex minima.

The code for the function written in Python code language is reported below.

Figure A2. Python code.

References

1. Bossi, A.H.; Cole, D.; Passfield, L.; Hopker, J. Conventional Methods to Prescribe Exercise Intensity Are Ineffective for Exhaustive Interval Training. Eur. J. Appl. Physiol.; 2023; 123, pp. 1655-1670. [DOI: https://dx.doi.org/10.1007/s00421-023-05176-6]

2. Meyer, T.; Lucía, A.; Earnest, C.P.; Kindermann, W. A Conceptual Framework for Performance Diagnosis and Training Prescription from Submaximal Gas Exchange Parameters—Theory and Application. Int. J. Sports Med.; 2005; 26, (Suppl. S1), pp. S38-S48. [DOI: https://dx.doi.org/10.1055/s-2004-830514] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15702455]

3. Cuomo, G.; Lorenzo, A.D.; Tramontano, A.; Iannone, F.P.; D’Angelo, A.; Pezzella, R.; Testa, C.; Parlato, A.; Merone, P.; Pacileo, M. et al. Exercise Training in Patients with Heart Failure: From Pathophysiology to Exercise Prescription. Rev. Cardiovasc. Med.; 2022; 23, 144. [DOI: https://dx.doi.org/10.31083/j.rcm2304144] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39076241]

4. Robertson, R.J.; Goss, F.L.; Dube, J.; Rutkowski, J.; Dupain, M.; Brennan, C.; Andreacci, J. Validation of the Adult OMNI Scale of Perceived Exertion for Cycle Ergometer Exercise. Med. Sci. Sports Exerc.; 2004; 36, pp. 102-108. [DOI: https://dx.doi.org/10.1249/01.MSS.0000106169.35222.8B] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/14707775]

5. Cauwenberghs, N.; Sente, J.; Sabovčik, F.; Ntalianis, E.; Hedman, K.; Claes, J.; Goetschalckx, K.; Cornelissen, V.; Kuznetsova, T. Cardiorespiratory Fitness Components in Relation to Clinical Characteristics, Disease State and Medication Intake: A Patient Registry Study. Clin. Physiol. Funct. Imaging; 2023; 43, pp. 441-452. [DOI: https://dx.doi.org/10.1111/cpf.12842]

6. Zimatore, G.; Serantoni, C.; Gallotta, M.C.; Guidetti, L.; Maulucci, G.; De Spirito, M. Automatic Detection of Aerobic Threshold through Recurrence Quantification Analysis of Heart Rate Time Series. Int. J. Environ. Res. Public Health; 2023; 20, 1998. [DOI: https://dx.doi.org/10.3390/ijerph20031998]

7. Zimatore, G.; Gallotta, M.C.; Innocenti, L.; Bonavolontà, V.; Ciasca, G.; De Spirito, M.; Guidetti, L.; Baldari, C. Recurrence Quantification Analysis of Heart Rate Variability during Continuous Incremental Exercise Test in Obese Subjects. Chaos; 2020; 30, 033135. [DOI: https://dx.doi.org/10.1063/1.5140455]

8. Mourot, L.; Tordi, N.; Bouhaddi, M.; Teffaha, D.; Monpere, C.; Regnard, J. Heart Rate Variability to Assess Ventilatory Thresholds: Reliable in Cardiac Disease?. Eur. J. Prev. Cardiol.; 2012; 19, pp. 1272-1280. [DOI: https://dx.doi.org/10.1177/1741826711423115]

9. Emerenziani, G.P.; Ferrari, D.; Vaccaro, M.G.; Gallotta, M.C.; Migliaccio, S.; Lenzi, A.; Baldari, C.; Guidetti, L. Prediction Equation to Estimate Heart Rate at Individual Ventilatory Threshold in Female and Male Obese Adults. PLoS ONE; 2018; 13, e0197255. [DOI: https://dx.doi.org/10.1371/journal.pone.0197255]

10. Deboeck, G.; Vicenzi, M.; Faoro, V.; Lamotte, M. Aerobic exercise capacity is normal in obesity with or without metabolic syndrome. Respir. Med.; 2023; 210, 107173. [DOI: https://dx.doi.org/10.1016/j.rmed.2023.107173]

11. Gouw, A.H.; Van Guilder, G.P.; Cullen, G.G.; Dalleck, L.C. Is the Tyme Wear Smart Shirt Reliable and Valid at Detecting Personalized Ventilatory Thresholds in Recreationally Active Individuals?. Int. J. Environ. Res. Public Health; 2022; 19, 1147. [DOI: https://dx.doi.org/10.3390/ijerph19031147] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35162172]

12. Magal, M.; Franklin, B.A.; Dwyer, G.B.; Riebe, D. Back to Basics: A Critical Review of the Methodology Commonly Used to Estimate Cardiorespiratory Fitness. ACSM’s Health Fit. J.; 2023; 27, pp. 12-19. [DOI: https://dx.doi.org/10.1249/FIT.0000000000000847]

13. Zimatore, G.; Gallotta, M.C.; Campanella, M.; Skarzynski, P.H.; Maulucci, G.; Serantoni, C.; De Spirito, M.; Curzi, D.; Guidetti, L.; Baldari, C. et al. Detecting Metabolic Thresholds from Nonlinear Analysis of Heart Rate Time Series: A Review. Int. J. Environ. Res. Public Health; 2022; 19, 12719. [DOI: https://dx.doi.org/10.3390/ijerph191912719] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36232025]

14. Zignoli, A.; Fornasiero, A.; Rota, P.; Muollo, V.; Peyré-Tartaruga, L.A.; Low, D.A.; Fontana, F.Y.; Besson, D.; Pühringer, M.; Ring-Dimitriou, S. et al. Oxynet: A collective intelligence that detects ventilatory thresholds in cardiopulmonary exercise tests. EJSS; 2022; 22, pp. 425-435. [DOI: https://dx.doi.org/10.1080/17461391.2020.1866081] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33331795]

15. Jiang, S.; Zhang, T.; Zhou, Y.; Lai, P.; Huang, Y. Wearable ultrasound bioelectronics for healthcare monitoring. Innov.; 2023; 4, 100447. [DOI: https://dx.doi.org/10.1016/j.xinn.2023.100447]

16. You, C.Y.; Hu, B.F.; Xu, B.R.; Zhang, Z.Y.; Wu, B.M.; Huang, G.S.; Mei, Y.F. Foldable-circuit-enabled miniaturized multifunctional sensor for smart digital dust. Chip; 2022; 1, 100034. [DOI: https://dx.doi.org/10.1016/j.chip.2022.100034]

17. Meucci, M.; Nandagiri, V.; Kavirayuni, V.S.; Whang, A.; Collier, S.R. Correlation Between Heart Rate at Maximal Fat Oxidation and Aerobic Threshold in Healthy Adolescent Boys and Girls. Pediatr. Exerc. Sci.; 2021; 33, pp. 139-143. [DOI: https://dx.doi.org/10.1123/pes.2020-0210]

18. Hagan, R.D.; Smith, M.G. Pulmonary Ventilation in Relation to Oxygen Uptake and Carbon Dioxide Production during Incremental Load Work. Int. J. Sports Med.; 1984; 5, pp. 193-197. [DOI: https://dx.doi.org/10.1055/s-2008-1025904]

19. Hollmann, W. 42 Years Ago—Development of the Concepts of Ventilatory and Lactate Threshold. Sports Med.; 2001; 31, pp. 315-320. [DOI: https://dx.doi.org/10.2165/00007256-200131050-00002]

20. Mishra, P.; Pandey, C.M.; Singh, U.; Gupta, A.; Sahu, C.; Keshri, A. Descriptive Statistics and Normality Tests for Statistical Data. Ann. Card. Anaesth.; 2019; 22, pp. 67-72. [DOI: https://dx.doi.org/10.4103/aca.ACA_157_18]

21. Ludbrook, J. Linear Regression Analysis for Comparing Two Measurers or Methods of Measurement: But Which Regression?. Clin. Exp. Pharmacol. Physiol.; 2010; 37, pp. 692-699. [DOI: https://dx.doi.org/10.1111/j.1440-1681.2010.05376.x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20337658]

22. Koo, T.K.; Li, M.Y. A Guideline of Selecting and Reporting Intraclass Correlation Coefficients for Reliability Research. J. Chiropr. Med.; 2016; 15, pp. 155-163. [DOI: https://dx.doi.org/10.1016/j.jcm.2016.02.012] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27330520]

23. Hopkins, W.G.; Marshall, S.W.; Batterham, A.M.; Hanin, J. Progressive Statistics for Studies in Sports Medicine and Exercise Science. MSSE; 2009; 41, pp. 3-13. [DOI: https://dx.doi.org/10.1249/MSS.0b013e31818cb278] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19092709]

24. Bland, J.M.; Altman, D.G. Statistical Methods for Assessing Agreement between Two Methods of Clinical Measurement. Lancet; 1986; 1, pp. 307-310. [DOI: https://dx.doi.org/10.1016/S0140-6736(86)90837-8]

25. Pühringer, M.; Ring-Dimitriou, S.; Iglseder, B.; Frey, V.; Trinka, E.; Paulweber, B. Sequencing Patterns of Ventilatory Indices in Less Trained Adults. Front. Sports Act Living; 2023; 4, 1066131. [DOI: https://dx.doi.org/10.3389/fspor.2022.1066131]

26. Zignoli, A. Machine Learning Models for the Automatic Detection of Exercise Thresholds in Cardiopulmonary Exercising Tests: From Regression to Generation to Explanation. Sensors; 2023; 23, 826. [DOI: https://dx.doi.org/10.3390/s23020826]

27. Ferretti, G.; Fagoni, N.; Taboni, A.; Vinetti, G.; di Prampero, P.E. A Century of Exercise Physiology: Key Concepts on Coupling Respiratory Oxygen Flow to Muscle Energy Demand during Exercise. Eur. J. Appl. Physiol.; 2022; 122, pp. 1317-1365. [DOI: https://dx.doi.org/10.1007/s00421-022-04901-x]

28. Ortega, F.B.; Leskošek, B.; Blagus, R.; Gil-Cosano, J.J.; Mäestu, J.; Tomkinson, G.R.; Ruiz, J.R.; Mäestu, E.; Starc, G.; Milanovic, I. et al. European Fitness Landscape for Children and Adolescents: Updated Reference Values, Fitness Maps and Country Rankings Based on Nearly 8 Million Test Results from 34 Countries Gathered by the FitBack Network. Br. J. Sports Med.; 2023; 57, pp. 299-310. [DOI: https://dx.doi.org/10.1136/bjsports-2022-106176]

29. Kaufmann, S.; Gronwald, T.; Herold, F.; Hoos, O. Heart Rate Variability-Derived Thresholds for Exercise Intensity Prescription in Endurance Sports: A Systematic Review of Interrelations and Agreement with Different Ventilatory and Blood Lactate Thresholds. Sports Med. Open; 2023; 9, 59. [DOI: https://dx.doi.org/10.1186/s40798-023-00607-2]

30. Schaffarczyk, M.; Rogers, B.; Reer, R.; Gronwald, T. Validation of a Non-Linear Index of Heart Rate Variability to Determine Aerobic and Anaerobic Thresholds during Incremental Cycling Exercise in Women. Eur. J. Appl. Physiol.; 2023; 123, pp. 299-309. [DOI: https://dx.doi.org/10.1007/s00421-022-05050-x]

31. Fletcher, G.F.; Balady, G.J.; Amsterdam, E.A.; Chaitman, B.; Eckel, R.; Fleg, J.; Froelicher, V.F.; Leon, A.S.; Piña, I.L.; Rodney, R. et al. Exercise standards for testing and training: A statement for healthcare professionals from the American Heart Association. Circulation; 2001; 104, pp. 1694-1740. [DOI: https://dx.doi.org/10.1161/hc3901.095960] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11581152]

32. Wannamethee, S.G.; Shaper, A.G. Physical Activity in the Prevention of Cardiovascular Disease: An Epidemiological Perspective. Sports Med.; 2001; 31, pp. 101-114. [DOI: https://dx.doi.org/10.2165/00007256-200131020-00003]

33. Giuliani, A.; Piccirillo, G.; Marigliano, V.; Colosimo, A. A Nonlinear Explanation of Aging-Induced Changes in Heartbeat Dynamics. Am. J. Physiol. Heart Circ. Physiol.; 1998; 275, pp. H1455-H1461. [DOI: https://dx.doi.org/10.1152/ajpheart.1998.275.4.H1455] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9746497]

34. Marwan, N.; Carmen Romano, M.; Thiel, M.; Kurths, J. Recurrence Plots for the Analysis of Complex Systems. Phys. Rep.; 2007; 438, pp. 237-329. [DOI: https://dx.doi.org/10.1016/j.physrep.2006.11.001]

35. Addo, P.M.; Billio, M.; Guégan, D. Nonlinear Dynamics and Recurrence Plots for Detecting Financial Crisis. N. Am. J. Econ. Financ.; 2013; 26, pp. 416-435. [DOI: https://dx.doi.org/10.1016/j.najef.2013.02.014]

36. Rajendra Acharya, U.; Paul Joseph, K.; Kannathal, N.; Lim, C.M.; Suri, J.S. Heart Rate Variability: A Review. Med. Biol. Eng. Comput.; 2006; 44, pp. 1031-1051. [DOI: https://dx.doi.org/10.1007/s11517-006-0119-0]

37. Censi, F.; Calcagnini, G.; Cerutti, S. Coupling Patterns between Spontaneous Rhythms and Respiration in Cardiovascular Variability Signals. Comput. Methods Programs Biomed.; 2002; 68, pp. 37-47. [DOI: https://dx.doi.org/10.1016/S0169-2607(01)00158-4]

38. Pitsik, E.; Frolov, N.; Hauke Kraemer, K.; Grubov, V.; Maksimenko, V.; Kurths, J.; Hramov, A. Motor Execution Reduces EEG Signals Complexity: Recurrence Quantification Analysis Study. Chaos; 2020; 30, 023111. [DOI: https://dx.doi.org/10.1063/1.5136246]

39. Webber, C.L., Jr.; Zbilut, J.P. Dynamical assessment of physiological systems and states using recurrence plot strategies. J. Appl. Phisiol.; 1994; 76, pp. 965-973. [DOI: https://dx.doi.org/10.1152/jappl.1994.76.2.965]

40. Wright, S.E.; Palmer, C. (03 May): Cardiac Dynamics in Auditory-Motor Synchronization: Roles of Short-Term Training and Rhythm Complexity. Audit. Percept. Cogn.; 2024; 7, pp. 86-109. [DOI: https://dx.doi.org/10.1080/25742442.2024.2345566]

41. Lin, M.; Zhao, G.; Wang, G. Recurrence Quantification Analysis for Detecting Dynamical Changes in Earthquake Magnitude Time Series. Int. J. Mod. Phys. C; 2015; 26, 1550077. [DOI: https://dx.doi.org/10.1142/S0129183115500771]

42. Chelidze, T.; Matcharashvili, T. Dynamical Patterns in Seismology. Recurrence Quantification Analysis: Theory and Best Practices; Webber, C.L., Jr.; Marwan, N. Springer International Publishing: Cham, Switzerland, 2015; pp. 291-334. ISBN 978-3-319-07155-8

43. Araújo, N.S.; Reyes-Garcia, S.Z.; Brogin, J.A.F.; Bueno, D.D.; Cavalheiro, E.A.; Scorza, C.A.; Faber, J. Chaotic and Stochastic Dynamics of Epileptiform-like Activities in Sclerotic Hippocampus Resected from Patients with Pharmacoresistant Epilepsy. PLoS Comput. Biol.; 2022; 18, e1010027. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1010027] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35417449]

44. Chatain, C.; Ramdani, S.; Vallier, J.M.; Gruet, M. Recurrence quantification analysis of force signals to assess neuromuscular fatigue in men and women. Biomed. Signal Process. Control; 2021; 68, 102593. [DOI: https://dx.doi.org/10.1016/j.bspc.2021.102593]