Validity of dynamical analysis to characterize heart rate and oxygen consumption during effort tests

Performance is usually assessed by simple indices stemming from cardiac and respiratory data measured during graded exercise test. The goal of this study is to characterize the indices produced by a dynamical analysis of HR and VO2 for different effort test protocols, and to estimate the construct validity of these new dynamical indices by testing their links with their standard counterparts. Therefore, two groups of 32 and 14 athletes from two different cohorts performed two different graded exercise testing before and after a period of training or deconditioning. Heart rate (HR) and oxygen consumption (VO2) were measured. The new dynamical indices were the value without effort, the characteristic time and the amplitude (gain) of the HR and VO2 response to the effort. The gain of HR was moderately to strongly associated with other performance indices, while the gain for VO2 increased with training and decreased with deconditioning with an effect size slightly higher than VO2 max. Dynamical analysis performed on the first 2/3 of the effort tests showed similar patterns than the analysis of the entire effort tests, which could be useful to assess individuals who cannot perform full effort tests. In conclusion, the dynamical analysis of HR and VO2 obtained during effort test, especially through the estimation of the gain, provides a good characterization of physical performance, robust to less stringent effort test conditions.

Scientific RepoRtS | (2020) 10:12420 | https://doi.org/10.1038/s41598-020-69218-1 www.nature.com/scientificreports/ Regarding standard analysis of respiratory parameters, the main indicators of athlete's performing capacities are the maximal VO 2 reached during the exercise, the maximal aerobic power or the maximal speed reached, and the values of power or speed at the two Ventilatory Thresholds (VTs), corresponding to the lactic apparition (VT1) and the accumulation (VT2) threshold 7 . Although these VO 2 parameters are currently considered among the best indices of aerobic fitness evaluation 8 , several drawbacks exist. First, determining them requires most of the time a visual analysis of the data. Second, they make use of only a part of the gas consumption dynamics, discarding the majority of the information contained in the entire effort test.
The second approach, based on dynamical system modeling, could allow to more accurately characterize the HR or VO 2 response during effort. Dynamical analysis based on differential equations is an active subject of research in the behavioral field since the seminal work of Boker 9 and has led to numerous studies in the field of psychology and to several methodological advances 10 . A first order differential equation approach has the potential ability to adjust HR measurement 11,12 and VO 2 dynamics during variable effort loads 13 . We propose here to use a simple first order differential equation coupled with a mixed effect regression to quantify the link between the exercise load during effort tests and the resulting HR or VO 2 dynamics. Because dynamical models use all the information measured during the effort test, it may allow to accurately assess performance using non-maximal effort tests.
The aim of this study is to characterize the indices produced by the dynamical analysis of HR and VO 2 for different effort test protocols. The construct validity of these new dynamical indices will be provided by testing the link with their standard counterpart. Their ability to detect performance change over two different contexts of training load will determine their predictive validity and sensitivity to change. We will therefore analyze longitudinal data measured for two groups of young athletes with two different protocols. One group should show a performance increase following a three months training period, and the second group should have a performance decrease after an off-season of 6 weeks. The possibility to apply the proposed dynamical analysis to submaximal effort tests will be studied by comparing the result of the analysis performed on the full tests with the one performed on only the first part of the test.

Subjects.
To test the reliability of the dynamical analysis model, data were acquired in two different populations (Guadeloupe and Spanish athletes) subjected to two different profiles of exercise (step-by-step cycling and continuous intensity running increase) and physiological conditions (training and deconditioning), presented in Table 1.
Group 1 consists of 32 young athletes (19 males and 13 females; 15.1 ± 1.5 year-old) of the Regional Physical and Sports Education Centre (CREPS) of French West Indies (Guadeloupe, France), belonging to a national division of fencing, or a regional division of sprint kayak and triathlon. GET was performed at the end of the off-competition season, and after 3 months of intense training (3-7 sessions/week). All athletes completed a medical screening questionnaire, and a written informed consent from the participants and the legal guardians was obtained prior to the study. The study was approved by the CREPS Committee of Guadeloupe (Ministry of Youth and Sports) and the CREPS Ethics Committee and performed according to the Declaration of Helsinki.
Group 2 consists of 14 young males, (15.4 ± 0.8 year-old) amateur soccer players from Malaga (Spain), performing three weekly training sessions and one weekly competition. A first GET was performed at the end of the soccer season and a second 6 weeks after. All participants were warned to avoid any training activity during this time. The measurements have been used in a previous publication 14 , they were approved by the Research Ethics Committee of the University of Málaga, Spain (EMEFYDE UMA: 2012-2015 report) and were carried out according to the principles of the Declaration of Helsinki. Participation in the study was voluntary, and prior to its initiation, written informed consent was obtained from the participants and the legal guardians of those under 18 years of age.
Effort test measurement. Group 1 performed an incremental testing on an SRM Indoor Trainer electronic cycloergometer (Schoberer Rad Meßtechnik, Jülich, Germany) associated to a Metalyzer 3B gas analyzer system (CORTEX Biophysik GmbH, Leipzig, Germany). The SRM cycloergometer is directly computer supervised to automatically maintain a constant mechanic workload by adjustment of the brake in accordance to the number of revolutions per minute. Cardiorespiratory parameters were recorded cycle-to-cycle during all the test to obtain HR and VO 2 all along the test session. The effort protocol used consisted of a 3 min rest phase, followed by a 3 min cycling period at 50 watts, followed by an incremental power testing of + 15 Watts by minute until exhaustion. At the end of the test, measurements were prolonged during a 3 min period to record the physiological recovery of athletes. www.nature.com/scientificreports/ Group 2 performed GET on a PowerJog J series treadmill connected to a CPX MedGraphics gas analyzer system (Medical Graphics, St Paul, MN, USA) with cycle-to-cycle measurements of respiratory parameters -including VO 2 , and HR-with a 12 lead ECG (Mortara). The stress test consisted of an 8-10 min warm up period of 5 km h −1 followed by continuous 1 km h −1 by minute speed increase until the maximum effort was reached. Power developed during the effort test was calculated using the formula described by the American College of Sport Medicine (ACSM). The latter determines an approximate VO 2 of runners 15 associated to the Hawley and Noakes equation that links oxygen consumption to mechanical power 16 . Truncated effort tests. In order to test the robustness of the dynamical analysis, truncated effort tests were generated from the maximal effort tests for both groups. It consisted in removing the measurements of the test for power (or speed) above 2/3 of the maximum power (or maximum speed) value, so that the maximum power (or speed) achieved during the truncated test lies between the two ventilatory thresholds. The recovery period was set as the recovery measurements of the full effort test with values below the maximum value reached during the truncated exercise. An example of truncated effort is presented in Fig. 1, for a VO 2 measurement during an effort test of group 1.

Standard indices.
The HRR calculated is the standard HRR60, which is the difference between the HR at the onset of the recovery and the HR 60 s later. The ventilatory thresholds 1 (VT1) and 2 (VT2) are calculated using the Wasserman method using the minute ventilation VE/VO 2 for determining VT1 and VE/VCO 2 for VT2 17 . The rHRI is derived by performing a sigmoidal regression of HR before and during the first 3 min effort step (only in group 1) and calculating the maximum derivative from the estimated parameters, as described in 18 . Maximum aerobic power (MAP) is the maximum power spent during the maximal effort test. HRmax and VO 2 max are the maximum values of the rolling mean of HR and VO 2 over 5 points.
New indices using dynamical analysis. A first order differential equation describes a relation between a time dependent variable, its change in time and a possible time dependent excitation mechanism. For a variable Y (HR or VO 2 ), it reads: where Ẏ (t) is the time derivative of Y (i.e. its instantaneous change over time), Y 0 its equilibrium value (i.e. its value in the absence of any exterior perturbation) and P(t) the excitation variable, that is the time dependent variable accounting for the exogenous input setting the system out of equilibrium. Equation 1 describes the dynamics of a self-regulated system that has a typical exponential response of characteristic time τ and an equilibrium value Y 0 in the absence of excitation (i.e. when P(t) = 0 ). For a constant excitation (i.e. a constant P(t) = P ), the system stabilizes at a value KP after several τ . This value depends on both the system and the excitation amplitude (see Fig. 2 left panel).
HR and VO 2 are two self-regulated features of our body: they respond to an effort with a certain characteristic time to reach a value corresponding to the energy demand 19 . Equation 1, as already demonstrated in 13 for VO 2 , can reproduce the dynamics of these two measures when considering that P(t) is the power developed by the body during effort. Assuming that HR or VO 2 follow Eq. 1, only three time-independent parameters are needed to characterize and to predict their dynamics for any time dependent effort:   www.nature.com/scientificreports/ at 10 km/h and who would have a total increase of HR of 60 beats/min for that effort, the decay time would be the time needed to increase his heartbeat by 38 beats/min (60 beats/min × 63%). • K , the gain, is the proportionality coefficient between a given effort increase and the corresponding total HR or VO 2 increase ( HR and VO 2 ). An illustration is provided in Fig. 2 left panel: a HR gain of K HR = 1 beat/ min/W leads to a HR increase of 100 beats/minute for a 100 W effort increase, and to HR = 200 beats/min for a 200 W effort increase.
An example of the dynamics for HR following Eq. 1 is given in Fig. 2 considering HR 0 = 50 beats min −1 , τ HR = 30 s, K HR = 1 and two efforts types. These three coefficients tightly characterize the dynamics of HR and allow us to generate the response to any effort.
The estimation of the three parameters characterizing the dynamics according to Eq. 1 is done in a two-step procedure, consisting in first estimating the first derivative of the variable studied over a given number of points with a Functional Data Analysis (FDA) regression spline method 10,20 . It consists in generating a B-spline function that fits the outcome to be studied and then estimating the derivative of that function. In order for the generated B-spline function to be differentiable, it needs to be smooth. This is achieved through a penalty function controlled by a smoothing parameter. This parameter was chosen to maximize the R 2 , which is the goodness of fit of the model to the data.
Once the derivative is estimated, a multilevel regression is performed to estimate the linear relation between the derivative, the variable and the excitation (summarized by the three parameters presented before).
This two-step estimation procedure has been extensively tested and described in a recent simulation study 12 . It can be applied to data with non-constant time sampling if it contains more than 5 points per typical decay time (in our case, at least one point every 20 s) and has a measurement noise below 50% of the signal amplitude.
Once the three dynamical parameters are estimated, an estimated curve can be reconstructed by performing a numerical integration of Eq. 1 (using the deSolve package in R 21 ).
This procedure (parameter estimation and estimated curve reconstruction) has been embedded and described in the open-source package doremi 22 available in the open source software R. Example code reproducing the analysis presented in this article can be found in the package vignettes.
Statistical analysis. HR measurements with a rate of change higher than 20 beat min −1 from one measurement to the next one were first removed as they were considered spurious results from the sensors.
Indices difference within each group between the first and the second measurement was assessed using paired t tests, and effect sizes were estimated by Cohen's d index. Associations between standard physical performance indices and the results of our dynamical analysis were assessed using Spearman rank correlation coefficients for continuous variables and logistic regression for dichotomous variables. Training was operationalized as a binary variable set to 0 for measurements before training for group 1 and after deconditioning for group 2 (untrained situation), and to 1 for measurements after training for group 1 and before deconditioning for group 2 (trained situation).
All analyses were performed using R version 3.4.2 23 , the package doremi 12,22 for the dynamical analysis and the packages data.table, Hmisc and ggplot2 for the data management and statistical indicators.

Results
The associations between standard indices were high, especially between the maximum value of oxygen consumption (VO 2 max), the MAP achieved and the ventilatory threshold powers for VO 2 (correlations ranging from 0.73 to 0.93). There was also a significant negative correlation between rHRI and VO 2 max (correlation coefficient of − 0.42, p = 0.023), meaning that a higher maximum aerobic power reached during effort or a higher  The estimated curve deviates from the experimental data mainly at high effort intensity and at rest before the effort. The ensemble of the estimated values compared to the true observed ones are presented in Supplementary Fig. 1 (online).
The dynamical analysis estimation of resting values overestimated the measured values (HR measures averaged approximately 20 s before the first effort increase, see Supplementary Table 2 online), partly because the participants did not provide enough values before the start of the test. Thus, we will discard this index for the rest of the study. VO 2 max and K VO 2 , the gain of VO 2 (i.e. proportionality coefficient between an effort increase and the final VO 2 increase caused by this supplementary energy expenditure), increased significantly during the 3 months training period of group 1, and decreased significantly during the 6 weeks of deconditioning of group 2 ( Table 2). The effect size was slightly higher for K VO 2 than VO 2 max in the two groups and was higher for deconditioning than for training for both variables.
A small decrease of the power of the first ventilatory threshold (power VT1) is also observed in population 2. τ VO 2 , the response time τ of VO 2 to the effort is shorter than τ HR , the response time of HR, in both populations. The HR gain ( K HR ) is remarkably similar in both groups, and unaffected by training or detraining. The relative standard deviation of τ VO 2 (between 35 and 50%) is higher than the relative standard deviation of the associated gain ( K VO 2 ). None of the dynamical parameters (gain K or decay time τ ) displayed significant correlation with the length of the experimental data record.
In univariable analysis, τ HR was correlated with measures of HRmax and HRR (Table 3), and K HR (i.e., proportionality coefficient between effort increase and final HR increase caused by this supplementary energy expenditure) was negatively correlated with weight, maximal aerobic power, maximum O 2 consumption, and the two ventilatory thresholds. Only in group 1, K HR was also negatively correlated with age, height and rHRI, whereas a correlation with HRmax is found only in group 2. In other words, a decrease of K HR , (i.e. a decrease of HR for a given effort) was linked with an improvement of oxygen maximal consumption, maximal aerobic power and the power corresponding at the two transition thresholds. Overall, correlations with new indices were higher than the correlations found between standard HR indices and other performance variables (see Supplementary Table 1 online). In a multivariable analysis performed in each group including age, weight, height, VO 2 max and power at ventilatory thresholds, only weight remained significantly associated with K HR (see Supplementary Table 3   www.nature.com/scientificreports/ VO 2 decay time ( τ VO 2 ) was globally independent of physiological variables and standard indices (Table 3), whereas K VO 2 was strongly associated with VO 2 max. In a multivariable analysis performed on each group including age, weight, height, training, VO 2 max and power at ventilatory thresholds, VO 2 max and training remained significantly associated with K VO 2 (see Supplementary Table 3 online). In group 1, training increased the K VO 2 of 1.1 mL/min/W on average and an increment of 1L/min of VO 2 max increased K VO 2 by 2.7 mL/min/W on average. In group 2, the deconditioning decreased K VO 2 by 2.1 mL/min/W and the decrease of 1L of VO 2 max lowered the VO 2 gain by 1.8 mL/min/W.

Truncated effort test.
When performing the dynamical analysis on the truncated effort tests (see Fig. 1 Fig. 4. The gains estimated on the truncated effort were slightly higher than the ones estimated on the entire effort test. Correlation between the gain (for VO 2 and HR) and the other performance indices remained similar to the ones observed in Table 3. The VO 2 gain K V 0 2 estimated on the truncated effort test still significantly changed between the two time points for both groups: from 8.9 (1.6) to 10.2 (1.8) mL/min/W for group 1 (p < 0.01, Cohen's d = 0.754), and from 15.0 (2.3) to 12.2 (1.7) ML/min/W for group 2 (p < 0.01, Cohen's d = 1.38). In summary, the VO 2 gain presented higher values but still significantly increased with training and decreased with deconditioning.

Discussion
Main findings. Modeling the evolution of HR and VO 2 during effort tests with a first order differential equation driven by the power spent during the effort, produced an estimation able to reproduce in average 95% of the observed variance of HR or VO 2 . The model was successfully tested in two different populations (Guadeloupe and Spanish athletes) subjected to two different profiles of exercise (step-by-step cycling and continuous intensity running increase) and physiological conditions (training and deconditioning). The dynamical analysis provided three indices: the equilibrium value or resting state, the decay time, and the gain or proportionality between a given effort increase and the corresponding total increase in HR or VO 2 . HR gain was correlated to the main indices of athlete's performance (MAP, VO 2 max, VT1 and VT2), which was not the case of other Table 2. Comparison of the classical indices and the indices stemming from the dynamical analysis of VO 2 and HR: the gain K and the decay time τ. Effort test measurements have been performed before and after the 3-month training in group 1 and before and after the 6-week deconditioning in group 2. VO 2 : O 2 consumption; HR: Heart Rate; MAP: Maximal Aerobic Power; HRR: Heart Resting Rate; rHRI: rate of Heart Rate Increase; VT: Ventilatory Threshold. Group 2 does not have rHRI because of the protocol used: the linear increase of power does not allow proper calculation of rHRI. Effect size is given by the Cohen's d of the t test.   www.nature.com/scientificreports/ standard HR indices. Furthermore, VO 2 gain was sensible to training or physical deconditioning. Finally, the indices obtained when modeling truncated effort tests (using about the first 2/3 of the effort test data) had similar characteristics, showing the robustness and usefulness of such approach to analyze incomplete effort tests. Such incomplete tests could occur due to lack of time but also when assessing older or sick individuals.

Standard indices.
Using standard performance indices, it was possible to assess the relevance of the training/deconditioning conditions used for this study. Results were in line with those obtained by other studies 6,19 , thus confirming the quality of the effort test results in the two groups of athletes. In particular, the relationships between ventilatory thresholds (VT), maximal aerobic power (MAP) and maximum oxygen consumption (VO 2 max), as well as the change in VO 2 max after 3 months of training and after 6 weeks of deconditioning, were in accordance with expected results 24 . VO 2 max variation was also more pronounced in the deconditioning group than in the training one, as reported in previous observations 25,26 . Concerning rHRI, the negative correlations with VO 2 max and MAP was reported previously and is due to a parasympathetic withdrawal with sympathetic activation causing a relatively slower HR increase in response to intensity increase for well-trained athletes when compared to untrained 3,6 .

Dynamical analysis.
There was a moderate correlation between VO 2 gain ( K VO 2 ) and VO 2 max 27 . Under an assumption of linearity between mechanical workload and O 2 consumption, VO 2 max corresponds to the oxygen consumption for the MAP expenditure and is directly linked to K VO 2 : However, VO 2 max is estimated via a single experimental measurement, supposed to be the VO 2 at the maximum effort achieved by the athletes. The ability to reach maximum capacities during effort test is subject to several internal and external factors such as athlete's engagement, mood state, fatigue and many others. Furthermore, the linear relation between energy demand and O 2 consumption may not hold for high power expenditure 28 , and thus VO 2 max may not be representative of physical performance for intermediate efforts.
In contrast, K VO 2 is estimated from the entire VO 2 dynamics during the effort test, yielding a robust estimate of the VO 2 response to effort. As a consequence, the VO2 gain estimated on truncated effort tests was still sensible to training and deconditioning and seems a promising performance index for submaximal effort test, such as those employed for patients suffering chronic disease or for elderly patients. The typical response time (i.e. the decay time τ ) of VO 2 was shorter than the HR one, in agreement with previous results 29 . This temporal delay of HR compared to VO 2 kinetics is due to the fact that heart flow regulation is partially driven by the oxygen demand of the organism detected via chemoreceptors, causing the HR increase to be a consequence of the VO 2 increase 30 .
The high variability of the decay time estimated, compared to the one of the gains, may find its roots in the wide range of the athletes' sport profile in our study. Indeed the different energetic profiles of the athletes according to their sport discipline 31,32 or soccer field position 33 could modify the kinetic of the VO 2 curve, leading to the variability of the decay time observed. On the other hand, the variability of the gain is the result of the aerobic metabolic efficiency, which is constant according to the substrate 34 , and the cycling or running efficiency, which is globally similar in a homogenous population of athletes.
The negative link between the K HR and subject weight may be explained by the known association between fat-free mass weight and heart's left ventricular size and mass 35 . This association, reflecting a well-trained heart in heavier athletes, results in a lower HR for a given effort and so a lower K HR .

Strength and weakness.
The main strength of this study is the use of two different populations of athletes, with two different effort tests and two different training schemes, showing its potential generalizability. Nevertheless, further study will need to extend these results to older adults, young children, and people with strong sedentary habits. A second strength is related to the analyses used, which allowed the estimation of performance indices without a maximum effort test. These analyses pave the way to obtaining accurate performance indices and information on training or deconditioning among larger groups of the population, such as the elderly, or patients at risk of cardiovascular events. The availability of ready to use, open source, tools for such analysis should facilitate its use for researchers and sport coaches 22 .
As for limitations, the dynamic model used in this study made three assumptions that led to slightly suboptimal fits. First, the assumption that the equilibrium value is constant before and after the effort does not hold and led to the overestimation of these values. Indeed, HR and VO 2 are known to decrease back to their resting value on a longer time scale due to the reduction of blood volume (i.e. dehydration), the evacuation of the heat accumulated during the muscular contractions, or the over-activation of the sympathetic system during exercise 36 . The second assumption is that the entire dynamics has one unique characteristic exponential time, making the model unable to account for cardiac drift associated to prolonged effort or any long-term modification of the variable dynamics. The third assumption is that the gain of VO 2 and HR is constant along the effort, i.e. that an increase of exerted power leads to the same final increase of HR or VO 2 . However, it is known that the VO 2 dynamics saturates at high effort intensities 37 and that the HR response to effort diminishes after the second ventilatory threshold (the inflexion point of the heart rate performance curve 38 ). The simple model proposed in this article cannot account for such changes on the dynamics, and instead estimates an average gain over the entire effort test. This leads to the increased difference between the estimated curve and the experimental data at high effort intensity. It also explains the higher gains estimated for the truncated effort test, which do not include the part of the effort where the real gain is actually diminishing.
(2) VO 2max = VO 2resting + MAP × K VO 2 Scientific RepoRtS | (2020) 10:12420 | https://doi.org/10.1038/s41598-020-69218-1 www.nature.com/scientificreports/ Possibility to release the restrictions listed above is of high interest and is the subject of current research. However, despite the fact that the model can still be improved, it already provides indices with good sensibility to performance change and cardio-respiratory indices used to measure fitness.

conclusion
The dynamical analysis of heart rate (HR) and oxygen consumption (VO 2 ) during effort appears to be relevant to evaluate performing capacities of athletes and their evolution. It reproduced in average 95% of HR or VO 2 dynamics using only three estimated cardiovascular indices. It was more sensitive to training and deconditioning than classic indices. Furthermore, its ability to extrapolate VO 2 and HR indices from truncated effort tests using only the first steps of the exercise could place it as a valuable tool to evaluate functional capacity from participants unwilling or unable to do maximal exercise testing.