Recurrence quantification analysis of heart rate variability to detect both ventilatory thresholds

Aims of this study were: to verify if Recurrence Quantification Analysis (RQA) of Heart Rate Variability (HRV) time series could determine both ventilatory thresholds in individuals with different fitness levels, and to assess the validity of RQA method compared to gas-exchange method (GE). The two thresholds were estimated in thirty young individuals during incremental exercise on cycle-ergometer: Heart rate (HR), Oxygen consumption (VO2) and Workload were measured by the two methods (RQA and GE). Repeated measures ANOVA was used to assess main effects of methods and methods-by-groups interaction effects for HR, VO2 and Workload at aerobic (AerT) and anaerobic (AnT) thresholds. Validity of RQA at both thresholds was assessed for HR, VO2 and Workload by Ordinary Least Products (OLP) regression, Typical Percentage Error (TE), Intraclass Correlation Coefficients (ICC) and the Bland Altman plots. No methods-by-groups interaction effects were detected for HR, VO2 and Workload at AerT and AnT. The OLP analysis showed that at both thresholds RQA and GE methods had very strong correlations (r >0.8) in all variables (HR, VO2 and Workload). Slope and intercept values always included the 1 and the 0, respectively. At AerT the TE ranged from 4.02% (5.48 bpm) to 10.47% (8.53 Watts) (HR and Workload, respectively) and in all variables ICC values were excellent (≥0.85). At AnT the TE ranged from 2.53% (3.98 bpm) to 6.64% (7.81 Watts) (HR and Workload, respectively) and in all variables ICC values were excellent (≥0.90). Therefore, RQA of HRV time series is a new valid approach to determine both ventilatory thresholds in individuals with different physical fitness levels, it can be used when gas analysis is not possible or not convenient.


Introduction
In exercise physiology the "threshold "concept is frequently used to assess both aerobic fitness and performance. The gas exchange (GE) method allows detecting ventilatory and gas exchange response and it is highly utilized for threshold identification [1,2]. Although the GE method represents the gold standard, more sustainable and inexpensive approaches can be easily applied to determine the thresholds when gas analysis is non convenient or not possible.
In a previous work we proposed a new non-linear method based on Recurrence Quantification Analysis (RQA) of Heart Rate Variability (HRV) time series to estimate the aerobic threshold (AerT) in obese subjects [3]. In this special population the AerT (or first lactate threshold according to Binder and colleagues [2]), represents a useful parameter to identify the most appropriate physical exercise intensity in order to reduce body weight (BW) and to improve physical fitness [4]. With increasing exercise intensity above the AerT, the lactate production rate becomes higher than the metabolizing capacity in the working muscle. This leads to an increased production of carbon dioxide and to a steeper increase of ventilation, while the increase of oxygen consumption continues to rise linearly with increasing workload. This onset of exercise-induced hyperventilation (i.e.an overproportioned increase ventilation in relation to carbon dioxide production) represents the Anaerobic threshold (AnT) or second lactate threshold [1,2]. Above the AnT the muscular lactate production rate exceeds the systemic lactate elimination rate. Therefore, this phase is represented by a more pronounced ventilation and a further increase in carbon dioxide production in an attempt to compensate for the marked rise in lactate [1,2,5].
The AerT and AnT are not always clearly distinguishable in subjects with low level of fitness such as obese individuals. This population cannot perform an incremental physical exercise for enough long time and consequently usually the AerT is usually the only ventilatory parameter that can be determined [6]. On the other side, the AerT and AnT can be more easily observed in a healthy population with a higher fitness level as athletes. Therefore, this population should be considered for developing a new method to estimate both ventilatory thresholds. Previous studies in athletes proposed several methods to provide information about the relationship between heart rate and the AerT and AnT [7][8][9][10][11][12]. RQA can be defined as a graphical, statistical and analytical tool used by several disciplines from physiology [13][14][15], to earth science [16] and economics [17,18] to detect phase transitions. The method based on RQA of HRV time series to distinguish the two ventilatory thresholds still needs to be established. Nowadays it is extremely easy to monitor the heart rate using low-cost, non-invasive, and mobile systems [19][20][21]. Since the use of these devices does not necessarily require the presence of the subject in the laboratory, from a logistic point of view, the RQA of HRV time series could be an extremely useful method for assessing sport performance and for planning training intensity in athletes in situations where gas analysis is not permitted.
The purpose of this study was to verify if the RQA of HRV time series proposed for obese individuals can be applied also to healthy young subjects. In particular, the first aim was to assess if the new non-linear method based on RQA could detect both ventilatory thresholds in individuals with different physical fitness levels. The second aim of the study was to investigate the validity of the RQA method compared to the GE method in thresholds detection.

Materials and methods
In this paper the two evaluation methods were called Gas Exchange method (GE) and Recurrence Quantification Analysis method (RQA). In addition, the two thresholds, Aerobic (AerT) and Anaerobic (AnT), were called Aerobic Gas Exchange Threshold (AerT GE ), Anaerobic Gas Exchange Threshold (AnT GE ), Aerobic Recurrence Quantification Analysis Threshold (Aer-T RQA ), and Anaerobic Recurrence Quantification Analysis Threshold (AnT RQA ).

Participants
A priori power analysis indicated that a total sample size of 30 subjects were required to detect a medium effect size (f = 0.25) given a coefficient of correlation p = 0.60 with 80% power and alfa = 0.05, using ANOVA repeated measure, within factors. Thirty subjects (2 females; 28 males) (age = 15.7 ± 2.7 years) were recruited for this study. Participants were divided in three groups: competitive rowers (group A, n = 8), recreational rowers (group B, n = 8) and other recreational sports (group C, n = 14). The A group was made up of competitive rowers, i.e. athletes who performed a minimum of 5 workouts per week and have participated in a regional and/or national competitions in the previous year. The subjects of group B attended rowing activities twice a week. Both the A and B groups were recruited from the "Circolo Canottieri Tirrenia Todaro" of Rome. The C group was made up of subjects who practice recreational activities different to rowing (twice a week).
All participants underwent clinical examinations to exclude any side effects to physical activity. However, a medical certificate for competitive sport activities or noncompetitive sport activities was requested. Moreover, subjects of all groups had the same main anthropometric characteristics except height (see Table 1). The exclusion criteria were neuropathy, autonomic dysfunction, cardiovascular diseases. All subjects (or their parents if <18 years) provided a written informed consent before the beginning of the study. This study was conducted in accordance with the Declaration of Helsinki and approved by the CAR-IRB-University of Rome "Foro Italico" Committee (Approval N˚CAR 37/2020).

Procedures
After clinical examination and anthropometric measurements, all subjects performed a graded incremental exercise test on cycle-ergometer. All subjects were tested in the morning (between 9:00 and 12:00 a.m.), and under similar environmental conditions (temperature 21-22˚C; humidity 50-60%). Subjects had their usual breakfast at least 90 minutes before the test session. All tests were performed at the Department of Movement, Human and Health Sciences at the University of Rome "Foro Italico".

Anthropometric measurements
The following participants' anthropometric measurements were assessed: body-weight, height, Body Mass Index (BMI) and percent of fat mass (FM%). Weight and height were measured using a scale and a stadiometer to the nearest 0.1 kg and 0.1 cm, respectively. BMI was calculated as the ratio between weight in kg and the square of height in meters (kg/m 2 ). FM% was

Incremental exercise test
The incremental exercise test on the cycle-ergometer was performed with the following protocol: participants started with a 1 minute rest period sitting on the bike, followed by 1 minute of unloaded pedaling (0 Watt). The workload was then increased by 20 Watts/minute for the A group (protocol 1) and by 15 Watts/minute for the B and C group (protocol 2), in order to maintain the total exercise time within about 15 minutes. Participants were asked to keep a cadence of 60-70 revolutions per minute (rpm). During the test, perception of physical exertion was assessed using OMNI Scale of Perceived Exertion (0-10 scale, [22]) 15 seconds prior to the end of every stage. The participants were asked to rate on a scale from 0 (extremely easy) to 10 (extremely hard) their subjective intensity of effort. The test ended when one of the following conditions was reached: a value of 10 on OMNI Scale of Perceived Exertion, the 90% of the subject's predicted HRmax (beats/min) or a respiratory exchange ratio equal to 1.1. The HR was recorded by a chest belt (HRM-Dual™, Garmin1) contemporarily to Oxygen consumption (VO 2 , mL/min), carbon dioxide production (VCO 2 , mL/min), and pulmonary ventilation (VE, mL/min) that were measured by an automatic gas analyzer (Quark RMR-CPET Cosmed™, Rome, Italy) [23] (See S2 Table).

Detection of the AerT and AnT (GE method)
Gas exchange method (GE) is the gold standard to detect both AerT and AnT as reported by Meyer et al. (2005) [1]. The AerT GE (that occurred at time named T 1 ), was determined offline for each subject by plotting the ventilatory equivalent of oxygen (VE/VO 2 ) as a function of VO 2 to identify the point where the VE/VO 2 reached its lowest value during the exercise test [24][25][26]. This procedure was confirmed by a graph with VCO 2 on the y axis and VO 2 on the x axis: two regression lines are fitted for the upper and the lower part of the relation and their intersection represents the AerT GE (V-slope method, [1]).
Similarly, the AnT GE (that occurred at time named T 2 ) was determined by plotting the ventilatory equivalent of carbon dioxide (VE/VCO 2 ) as a function of VO 2 to identify the point where the VE/VCO 2 reached its lowest value during the incremental test. This procedure was confirmed by a graph with VE on the y axis and VCO 2 on the x axis: two regression lines are fitted for the upper and the lower part of the relation and their intersection represents the AnT GE [1].
HR e VO 2 values corresponding at the two thresholds were reckoned as mean values of the last 30 seconds of the workload.

Recurrence Quantification Analysis method (RQA)
RQA-based method is widely explained in a previous paper [3]: it was observed that by recurrence analysis, in epoch-by-epoch (with sliding windows) mode (RQE), the rapid shifts from high to low (and vice versa) of percent of determinism (DET, the percentage of recurrence points which form diagonal lines) is usually an indicator of regime changes and phase transitions [27,28]. In fact, DET (High DET means high autocorrelation) is a universal marker of crisis in fields ranging from physiology to finance [29]. Laminarity (LAM, the percentage of recurrence points which form vertical lines), analogous to DET, measures the number of recurrence points which form vertical lines and indicates the amount of laminar phases (intermittency) in the system studied. In the present work we use both DET and LAM to study the increase in correlation and to identify the ventilatory thresholds; the present study was focused on chaos-order transitions and physical fatigue was considered as an order parameter acting at physiological level. More details can be found at http://www.recurrence-plot.tk/ and in [28,[30][31][32].

Data preprocessing
As previously explained by   [3], since the heart rate (HR) was continuously recorded and collected breath-by-breath, the time series of RR interval (temporal variation between the sequences of consecutive heartbeats) was reckoned by the ratio 60000/HR. Moreover, since HR physiologically increases in relation to the workload, this physiological trend was removed in order to analyze the cardiac regime changes (detrended RR interval) for each subject (the theoretical value of best linear fit on original time series were subtracted to every point belong to original time series). The optimization procedure of input parameters was discussed in a previous work [3], where the following same values were used: • delay (lag) is set to 1, • embedding dimension is set to 7, • cut-off distance (radius) to 50% of mean distance between all pairs of points in time, • line (minimal number of consecutive recurrences to score a determinism line) set to 4.
RQA is computed across consecutive distance matrices corresponding to consecutive and overlapping sliding windows (epochs) along the series, and this mode of analysis is called RQE (Recurrence Quantification by Epochs). The occurrence of an abrupt change of DET and LAM in adjacent windows corresponds to a transition in the dynamical regime of the signal. RQE analysis was carried out by adopting the same parameters setting as for the global mode (in the present study, when an epoch coincides with all points of the time series, the result is called RQA), plus the definition of windows having a length of 100 points (1 point correspond to 1 breath) and shifting of 1 point between consecutive windows, respectively.
RQA software, RQE.exe and RQC.exe version 8.1, are those included in RQA 14.1 [30]. RQE on sliding windows partially overlapped are carried on a detrended RR interval and then RQA measures are obtained. The software used for the RQA analysis, that generates Recurrence plots (RP) and other RQA utilities, is available from http://cwebber.sites.luc.edu.

Detection of the AerT and AnT (RQA method)
The method of threshold detection based on RQA of RR intervals (HRV time series) can be described as follows: specifically, AerT RQA was determined by the time point when the statistically relevant minimum of percent of determinism and laminarity occurred. AnT RQA , instead, occurred when percent of determinism is saturated (the saturation occur when DET was above 90% and the maximum value did not change for at least 20 seconds); both threshold time points correspond to evident changes in RP texture (see Fig 1). HR e VO 2 values corresponding to the two thresholds were reckoned as mean values on the last 10 values of the original time series before the value corresponding to the two threshold times (T 1 and T 2 ).

Statistical analysis
Repeated measures analysis of variance (RM ANOVA) was used to assess main effects of methods (RQA vs GE) and methods-by-groups (A, B and C) interaction effects for HR, VO 2 and Workload at the AerT and the AnT. Agreements between RQA and GE methods at the AerT and the AnT were assessed for HR, VO 2 and Workload by Ordinary Least Products (OLP) regression analysis [33]. Coefficients of determination (R 2 ) and regression parameters (slope and intercept) with the 95% of confidence intervals (95% CI) were calculated for the OLP regression equation to determine fixed and proportional biases. The hypothesis of proportional and fixed bias was rejected when the 95% CI contained the value 1 for the slope and the 0 for the intercept, respectively. Percentage differences between the RQA and GE methods at the AerT and the AnT (reported as mean and range values) were calculated for HR, VO 2 and Workload. RQA method validity was also assessed by comparing HR, VO 2 and Workload at the AerT and the AnT vs the same variables measured by the GE method with a paired samples t-test. Measurement error was expressed in Typical percentage Error (TE) which was calculated by dividing the standard deviation of the difference percentage by ffi ffi ffi 2 p . Intraclass correlation coefficients (ICC) was used as a parameter for criterion validity of the RQA method The red line enlightens the statistically relevant minimum that corresponds to the AerT and the green line corresponds to the AnT where the percent of determinism reaches saturation. Subject #1 is a competitive rower (group A, protocol 1) and subject #3 is a recreational rower (group B, protocol 2). The test starts at 00:00; AerT occurred at time T 1 , 8:38 and 5:49 for #1 and #3, respectively; AnT occurred at time T 2 , 11:24 and 11:01 for #1 and #3, respectively. A square pattern in the middle is easily distinguishable, the vertical red line is manually placed nearby before the first ventilatory transition, the green one nearby before the second one. The colors of RP reflect different Euclidean distances between trajectories and are analogous to geographical relief maps going from blue to red at increasing Euclidean distance. Distances greater than cut-off correspond to black regions.
https://doi.org/10.1371/journal.pone.0249504.g001 compared to the GE method at the AerT and the AnT for all variables (HR, VO 2 and Workload).
Bland-Altman plots were applied to determine the 95% limits of agreement (LoA) between the RQA and GE methods at the AerT and the AnT for the HR, VO2 and Workload [36].
Statistical significance was defined as p �0.05. All statistical analysis was performed by SPSS version 24.0 software (SPSS Inc., Chicago, IL).

Results
In Table 2 all the main parameters at the two thresholds and at the peak exercise for each group, are reported.
Workload at the AerT and HR, VO 2 , and Workload at the AnT were not statistically different between RQA and GE methods. Significant differences were found between the two methods for HR at the AerT (HR at the AerT GE = 137.77 ± 15.25 beats/min; HR at the AerT RQA = 141.44 ± 15.35 beats/min; p = 0.026 d = 0.24, small) and for VO 2 at the AerT (VO 2 at the Table 2. Main parameters for group (Mean value ± SD).

HR, VO 2 and workload at AerT
The agreements between HR, VO 2 and Workload at the AerT GE and the AerT RQA are presented in Table 4. In all variables (HR, VO 2 and Workload) OLP regression analysis showed that the AerT GE and the AerT RQA had very strong correlations (r >0.8). Slope and intercept values always included the 1 and the 0, respectively. Mean percentage differences were < 6%. HR and VO 2 at the AerT resulted statistically different between the two evaluation methods (p <0.05) while no differences were assessed for Workload at the AerT. The TE for HR and VO 2 at the AerT was considered acceptable (< 10%), moreover the TE for the Workload at AerT was slightly higher than the acceptable 10% limit [35]. ICC values were excellent (� 0.85) for all variables at the AerT [36,37]. The OLP regression analysis plots of HR, VO 2 and Workload values at the AerT GE and at the AerT RQA are graphically shown in Fig 2A-2C, respectively. In each graph it is described the OLP regression plot with the linear regression (solid line), the identity (dashed line), the equation, the correlation coefficient (r) and the absolute mean differences. For the HR, the VO 2 and the workload the Bland-Altman plots with the mean difference in the solid line and the 95% LoA in the dashed lines are shown in Fig 3A-3C, respectively.

HR, VO 2 and workload at AnT
From Table 2 the value of RER at peak of exercise was higher than the RER at AnT in all three groups and it confirm that exercising to the individual peak capacity allows to measure the AnT in all subjects. The agreement between HR, VO 2 and Workload at the AnT GE and the AnT RQA are presented in Table 4. In all variables (HR, VO 2 and Workload) OLP regression analysis showed that the AnT GE and the AnT RQA had very strong correlations (r >0.8) [39]. Slope and intercept values always included the 1 and the 0, respectively. Mean percentage differences were <3% and all variables corresponding to the AnT did not result statistically different between the two evaluation methods. For HR, VO 2 and Workload at the AnT the TE was considered acceptable (<10%) and ICC values were excellent (�0.90) [36,38]. The OLP regression analysis plots of HR, VO 2 and Workload at the AnT and at the AnT are graphically shown in Fig 2D-2F, respectively. For the HR, the VO 2 and the workload the Bland-Altman plots with the mean difference in the solid line and the 95% LoA are shown in Fig 3D-3F, respectively.

Discussion
The first aim of the study was to verify if in individuals with different physical fitness levels the new non-linear method based on RQA could detect both ventilatory thresholds. The second aim was to verify the validity of the RQA method compared to the GE method in thresholds detection. To our knowledge, this is the first study assessing the accuracy of the RQA method during the determination of both the AerT and the AnT. In particular, in the present investigation we examined the RQA epoch-by-epoch procedure in order to analyze HRV time series during phase transitions.
Even though for HR and VO 2 at the AerT significant differences were found between the two methodologies (RQA and GE methods), no differences were detected for Workload at the AerT and for HR, VO 2 and Workload at the AnT. Differences between RQA and GE methods assessed using effect size statistics (as Cohen's d) were trivial (<0.2) for all parameters except for HR at AerT (small effect size = 0.24). Furthermore, the absence of methods-by-group interactions in all variables at both the thresholds indicated that the RQA method is not influenced by the level of physical fitness of the subjects. Therefore, this method can be applied in individuals practicing different sports and with different physical fitness levels in situations where gas analysis is not available. These results permitted us to move to the second aim of the study which consisted in verifying the validity of the RQA method against the GE method during thresholds identification through different validation parameters (OLP regression, TE, ICC, and the Bland Altman plots) individually for HR, VO 2 and Workload variables corresponding at the AerT and the AnT.
Regarding the two methods of AerT detection (AerT GE and AerT RQA ), the OLP regression in all variables (HR, VO 2 and Workload) showed very strong correlations [39]. Neither fixed or proportional biases were present in all parameters and mean percentage differences were less than 5% for HR and Workload and less than 6% for VO 2 . As observed in the RM ANOVA results, HR and VO 2 at the AerT showed significant differences between the two methods, even if the differences were respectively small and trivial. Indeed, both these parameters showed an acceptable TE (HR 4.02%; VO 2 8.40%) and excellent ICC values (HR 0.85; VO 2 0.90) [36][37][38][39]. Therefore, for all these reasons the statistical differences detected between the two methods can be considered not relevant for practical applications. In particular, the low percentage of TE in the heart rate indicated that this variable at the AerT RQA represents the most accurate parameter to be monitored. To conclude, the Workload corresponding to the AerT did not result statistically different between the two evaluation methods and showed in the ICC an excellent value of 0.95. However, this parameter showed a TE slightly higher than the acceptable limit (10.47%) [39].
These first results about the validity of the RQA method for the detection of the AerT could be summarized stating that the AerT RQA is a valid method to use as alternative of gas analysis. In particular, the heart rate at AerT is considered the most valid parameter to monitor using the RQA method, and the determination of the AerT with this approach loses accuracy with VO 2 and Workload variables. Between the different exercise variables, the HR seems to be not only the most accurate but also the most sustainable, indeed this parameter can be easily recorded during training using non-invasive, not expensive, time-efficient devices which can be applied routinely and simultaneously in a large number of athletes [19]. In the field of exercise assessment, the heart rate variability has been largely used to identify the AerT in several specific populations as diabetic and cardiac patients [40,41], healthy individuals [42][43][44], competitive cyclists and triathletes [45], and professional soccer players [46]. However, little is known about the accuracy of the alternative approach of the RQA of HRV time series. Indeed, this method so far has been used only for the identification of the AerT in obese subjects [3]. The present work enlarges the boundaries of this methodology confirming the efficacy of RQA analysis also in individuals with a higher physical fitness level. Therefore, the utilization of the AerT RQA may provide an important training intensity guidance for the identification of the exercise zone which can be covered predominantly by the aerobic metabolism [1,5,25].
Even regarding the two methods of AnT detection (AnT GE and AnT RQA ), the OLP regression in all parameters (HR, VO2 and Workload) showed very strong correlations. [39] and neither fixed or proportional biases were detected. In addition, mean percentage differences were less than 1% for HR and VO 2 and less than 3% for Workload and no significant differences were found in all variables. Moreover, the acceptable TE values (HR 2.53%; VO 2 4.67%, Workload 6.64%) and the excellent ICC values (HR 0.90; VO 2 0.97; Workload 0.97) confirmed the validity of all parameters corresponding at the AnT RQA [35][36][37][38]. As it was observed in the AerT, the lowest TE found in HR demonstrated that this parameter is the most accurate to use for the detection of the AnT with the RQA method.
Based on the aforementioned findings we believe that the RQA approach is a valid method to assess AnT. In addition, comparing the accuracy of the different variables, the heart rate at the AnT seems to be the most valid and sustainable parameter to use with the RQA method. The HRV has been utilized to estimate the AnT in different exercise evaluation fields [9-11, 41, 45, 46]. However, to our knowledge, this is the first study which attempted to directly examine the non-linear approach of the RQA of HRV time series in the estimation of the AnT. Considering the high importance that the AnT detection has for sport performance [1,7,8], these results may provide a useful contribution in training prescription and evaluation.
Comparing the validity of the AerT RQA and AnT RQA (with OLP regression, TE, ICC and the Bland Altman plots) it appeared that the AnT RQA has a better accuracy than AerT RQA and other studies showed similar results [9][10][11]45]. This lower accuracy in the AerT RQA may be explained by methodological problems in the determination of this threshold using gas analyzers. Indeed, it seems that the accuracy, the intra-and inter-observer reliability, the repeatability and the percentage of indeterminate thresholds of the AerT GE had been questioned [1,47]. Therefore, the lower correspondence between AerT RQA and AerT GE may be caused by methodological problems on the GE method in some particular individuals. In view of all these aspects, the use of the RQA method could be used as a supportive approach to gas analysis for all the cases in which the detection of the AerT is ambiguous.
Some limitations may affect our data analysis and have to be mentioned: i. the number of participants of each group is relatively small, although the total sample size was determined by an "a priori" power analysis; ii. the RQA is a non-linear method in which the correct choice of the parameters is not always trivial and can influence the results [32,48]; iii. to synchronize hearth rates and amounts of exchanged gas (to be sure that HR was recorded at the same change of workload), we detected HR values breath-by-breath. Even if, in principle, this procedure could have had an influence on our data since breathing frequency and depth are connected to HRV, we could detect both thresholds similarly to other studies that utilized beat-by-beat HR detection [11,45,49]; iv. our approach was validated using a step protocol, with a workload increase every 60 sec.
Therefore, further investigations including a larger number of subjects, a different sampling of HR (i.e. beat-by-beat) together with a test and re-test procedure, allowing to determine results reliability, are needed. To conclude, future physiological exercise studies should assess if the findings of the present work could be adopted in field tests which do not necessarily require the use of a cycle-ergometer and in other populations different than young, healthy and active individuals.

Conclusion
The recurrence quantification analysis of heart rate variability time series is a new approach for the determination of ventilatory thresholds in individuals with different physical fitness levels, therefore, it can be used as a valid method for threshold detection. Thus, during an incremental exercise test it is possible to assess sport performance and delineate the intensity of training zones using non-invasive and low-cost devices. Based on the present results, coaches could obtain information about the training trend of their athletes even when gas analysis is non convenient or not possible.