Effect of Mathematical Modeling and Fitting Procedures on the Assessment of Critical Speed and Its Relationship With Aerobic Fitness Parameters

An accurate estimation of critical speed (CS) is important to accurately define the boundary between heavy and severe intensity domains when prescribing exercise. Hence, our aim was to compare CS estimates obtained by statistically appropriate fitting procedures, i.e., regression analyses that correctly consider the dependent variables of the underlying models. A second aim was to determine the correlations between estimated CS and aerobic fitness parameters, i.e., ventilatory threshold, respiratory compensation point, and maximal rate of oxygen uptake. Sixteen male runners performed a maximal incremental aerobic test and four exhaustive runs at 90, 100, 110, and 120% of the peak speed of the incremental test on a treadmill. Then, two mathematically equivalent formulations (time as function of running speed and distance as function of running speed) of three different mathematical models (two-parameter, three-parameter, and three-parameter exponential) were employed to estimate CS, the distance that can be run above CS (d′), and if applicable, the maximal instantaneous running speed (smax). A significant effect of the mathematical model was observed when estimating CS, d′, and smax (P < 0.001), but there was no effect of the fitting procedure (P > 0.77). The three-parameter model had the best fit quality (smallest Akaike information criterion) of the CS estimates but the highest 90% confidence intervals and combined standard error of estimates (%SEE). The 90% CI and %SEE were similar when comparing the two fitting procedures for a given model. High and very high correlations were obtained between CS and aerobic fitness parameters for the three different models (r ≥ 0.77) as well as reasonably small SEE (SEE ≤ 6.8%). However, our results showed no further support for selecting the best mathematical model to estimate critical speed. Nonetheless, we suggest coaches choosing a mathematical model beforehand to define intensity domains and maintaining it over the running seasons.


INTRODUCTION
The prescription of exercise intensity, one of the most important criteria to induce specific adaptations to training (Maclnnis and Gibala, 2017), is often based on the percentage of the maximal rate of oxygen uptake (VO 2 max) or maximal heart rate (American College of Sports Medicine, 2000;Burgomaster et al., 2007;Roy et al., 2018). However, among individuals, the lactate threshold, the respiratory compensation point (RCP), and critical power (CP)/speed (CS) were located at different percentages oḟ VO 2 max (Fontana et al., 2015), leading to substantial differences between participants in terms of characteristics of metabolic responses and duration of exercise at a common percentage of the maximum. Therefore, using an exercise prescription based on percentages of maximum values does not guarantee control of exercise intensity (DiMenna and Jones, 2009;Lansley et al., 2011). Instead, a model considering exercise intensity domains for exercise prescription has been recommended (Iannetta et al., 2020). Parameters such as oxygen uptake kinetics (Whipp and Mahler, 1980), ventilatory threshold (VT) (Wasserman et al., 1973), maximum lactate steady-state (Iannetta et al., 2018), and CP/CS (Vanhatalo et al., 2007;Constantini et al., 2014;Jones et al., 2019) can be used to define these various intensity domains.
CP/CS allows defining of the boundary between heavy and severe intensity domains (Jones et al., 2019;Galán-Rioja et al., 2020). Therefore, having an accurate estimation of CP/CS is important. This is obtained by fitting the experimental data to a mathematical model, chosen among several possibilities that differ with respect to their mathematical forms and number of parameters (Monod and Scherrer, 1965;Wilkie, 1980;Moritani et al., 1981;Whipp et al., 1982;Morton, 1986Morton, , 1990Morton, , 1996Peronnet and Thibault, 1989). The original linear model formulation was proposed by Monod and Scherrer (1965). This model was applied to cycle ergometry and relates the work performed during an exhaustive bout and its duration through two parameters (two-parameter model): CP (Monod and Scherrer, 1965) or threshold of fatigue (Bigland-Ritchie and Woods, 1984) and the sustainable work of exercise above that metabolic rate (W ) (Monod and Scherrer, 1965). Power has been related to time by dividing the original formulation by the exercise duration (Poole et al., 1986;Gaesser and Wilson, 1988;Housh et al., 1989) while Gaesser et al. (1990) proposed expressing this exercise duration as function of power, which led to the well-known hyperbolic formulation (Morton and Hodgson, 1996). Another model variant, proposed by Morton (2006), expresses the work performed as function of power, since this work (power multiplied by time to exhaustion) is also a dependent variable. However, this model has, to our knowledge, never been used so far.
Major shortcomings of the two-parameter model are the assumptions 1) of infinite running speed as time to exhaustion approaches zero, and 2) that at the point of fatigue, d has been completely covered (Gaesser et al., 1995;Morton, 1996). To overcome these limitations, Morton (1996) proposed a threeparameter model including an additional parameter, the maximal instantaneous running speed (s max ), and a d that can be only partly covered for a running speed between CS and s max . Alternatively, Hopkins et al. (1989) proposed a three-parameter exponential model based on CS and s max , but where d was replaced by an undefined time constant (τ). The authors reported that their three-parameter exponential model gave better fits than the two-parameter model for inclined treadmill running of short duration (<3 min) (Hopkins et al., 1989). These two-or three-parameter models can be formulated as either distance as function of time, time as function of distance, running speed as function of time, time as function of running speed, distance as function of running speed, and running speed as function of distance, which are mathematically equivalent.
To obtain a statistically appropriate estimation of the model parameters, the correct choice of model formulation and regression analysis should be chosen (Patoz et al., 2021). Such choice is based on the data provided by the experiment and the knowledge of the independent and dependent variables. For the treadmill CS test, running speed is the independent variable while time to exhaustion and distance (implicitly, because it is given by running speed multiplied by time to exhaustion) are the dependent variables. To minimize the error of a model formulation expressing the dependent and independent variables on the vertical and horizontal axes, respectively, the least squares (LS) loss function can be used and requires that the dependent variable be observed with additive error while the independent one would have no additive error (Morton and Hodgson, 1996). Statistical theory has shown that errors in the independent variable are of minor importance, making error minimization in the dependent variable sufficient (Morton and Hodgson, 1996). However, due to heteroscedasticity of the dependent variable (McLellan and Skinner, 1985;Poole et al., 1988;Faude et al., 2017), Morton and Hodgson (1996) suggested to use weighted LS (WLS).
Several researchers have compared the estimation of the parameters provided by the three different models (twoparameter, three-parameter, and three-parameter exponential) and some of their different formulations for cycle ergometry (Gaesser et al., 1995;Bull et al., 2000;Bergstrom et al., 2014) and running on a treadmill (Housh et al., 2001). Significant differences were obtained between the different formulations of the two-parameter model (Gaesser et al., 1995;Bull et al., 2000;Housh et al., 2001;Bergstrom et al., 2014). The three models also differed significantly from one another and the threeparameter model gave the lowest estimation of CP (Gaesser et al., 1995;Bull et al., 2000;Bergstrom et al., 2014) and CS (Housh et al., 2001). However, these studies (Gaesser et al., 1995;Bull et al., 2000;Housh et al., 2001;Bergstrom et al., 2014) did not consider time to exhaustion as the dependent variable, as honestly highlighted by Gaesser et al. (1995). Moreover, these previous studies (Gaesser et al., 1995;Bull et al., 2000;Housh et al., 2001;Bergstrom et al., 2014) are not methodologically exhaustive. Indeed, none of these studies acknowledged heteroscedasticity of the dependent variable.
Hence, the purpose of this study was twofold. First, we compared the estimations of the model parameters obtained by statistically appropriate fitting procedures (the combination of model formulation and regression analysis) applied to the three different models (two-parameter, three-parameter, and threeparameter exponential). We hypothesized that the estimations of CS, d , and s max would be significantly different between the mathematical models employed, but not between the fitting procedures. We also hypothesized that the three-parameter model would give the lowest estimation of CS, as already observed by Housh et al. (2001) for statistically inappropriate fitting procedures. Second, we determined the correlations between estimated CS and aerobic fitness parameters, i.e., VT, RCP, andVO 2 max, as well as the standard error of estimate (SEE) of these relations. We hypothesized that lower quality of the fit [determined by Akaike information criterion (AIC)] would be associated with lower correlations between CS and aerobic fitness parameters and higher SEE.

Participant Characteristics
Sixteen male runners gave written informed consent to participate in the present experiment (age: 25.6 ± 3.9 years old; height: 1.79 ± 0.05 m; body mass: 69.2 ± 5.3 kg). For study inclusion, participants were required to be in good selfreported general health with no symptoms of cardiovascular disease or major coronary risk factors, no current or recent lowerextremity injury that could prevent them from giving 100% of their capacity during the test or from meeting a certain level of running performance. More specifically, runners were required to have a speed associated withVO 2 max (sV O 2 max ) greater or equal to 4.44 m/s (16 km/h). The study protocol was approved by the Ethics Committee (CER-VD 2018-01814) and adhered to the latest Declaration of Helsinki of the World Medical Association.

Experimental Procedure
Each participant completed five experimental sessions interspersed by at least 2 days in the laboratory. All participants were advised to avoid strenuous exercise the day before a test but to maintain their usual training program otherwise. During the first session, participants completed a maximal incremental aerobic test on a treadmill (Arsalis T150-FMT-MED, Louvainla-Neuve, Belgium). This test consisted of a 10-min warm-up at 2.78 m/s followed by an incremental increase in the running speed of 0.28 m/s every 2 min until exhaustion. Throughout the test, participants breathed into a mask connected to a gas analyzer (Quark, Cosmed, Italy). Pulmonary gas exchange variables [expired minute ventilation (VE), oxygen uptake (VO 2 ), and carbon dioxide output (VCO 2 )] were measured breath-by-breath and subsequently averaged over 10-s intervals throughout the test. Before each test, the O 2 and CO 2 analyzers were calibrated using room air and known concentrations of calibration gas (16.00% O 2 , 5.02% CO 2 , and the remainder N 2 ), and the turbine was calibrated using a 3-L syringe (Hans Rudolph, Germany).
This test was used, first, to determine the peak speed (PS) of the incremental test of each participant. PS is defined as the running speed of the last fully completed increment (sV O 2 max ) plus the fraction of time spent in the following uncompleted increment (α) multiplied by the running speed increment ( s = 0.28 m/s) (Kuipers et al., 2003): PS = sV O 2 max + α s. Second, theVO 2 max was defined as the highest measuredVO 2 value corresponding to (1) a plateau ofVO 2 with increased running speed ( V O 2 between the last two increments smaller than 50% of the average V O 2 during the submaximal phase of the test) and/or (2) an heart rate greater than 90% of the theoretical maximum heart rate given by 220-age associated with a respiratory quotient greater than 1.1 and a rate of perceived exertion greater than 17. Third, VT and RCP were determined based on gas exchange data and using the method proposed by Wasserman et al. (1973).
The other four tests were performed in a randomized order and consisted of exhaustive runs at a given percentage of the participant's PS (90, 100, 110, or 120%). These tests were as follows: after a 10-min warm-up at 2.78 m/s and a 5-min rest period, the running speed was increased to a given percentage of PS, and the participant had to maintain the pace until exhaustion. The time to exhaustion was collected for each of the four sessions. No information about the timings or running speed was given to any participant, who were strongly encouraged, during any of the five experimental sessions. All participants were familiar with running on a treadmill.

Mathematical Modeling
The estimations of CS, d , and s max were obtained from two different but mathematically equivalent formulations for the three different models. Gaesser et al. (1990) proposed the two-parameter model formulation given by Eq. 1 (non-linear, time-running speed) while Eq. 2 (non-linear, distance-running speed) represents the formulation proposed by Morton (2006). The three-parameter model formulation proposed by Morton (1996) and the inverse of the three-parameter exponential model formulation proposed by Hopkins et al. (1989) are given by Eqs. 3 and 5 (non-linear, time-running speed), respectively, while Eqs. 4 and 6 (non-linear, distance-running speed) represent their distances as a function of running speed formulations.
t, s, and d stand for time, running speed, and distance, respectively. Of note, distance as a function of running speed was simply given by multiplying time as function of running speed by running speed, i.e., d(s) = s t(s).
The three-parameter exponential model does not provide a direct estimation of d because the distance that can be run above CS is time-dependent in such a model. Indeed, rearranging the two-parameter model formulation proposed by Whipp et al. (1982) and given by Eq. 7 (i.e., the inverse of Eq. 1) Then, applying this result to the three-parameter exponential model gives an equation where the left-hand side is time-dependent, i.e., of this equation appears where its first derivative is equal to zero, which is at t = τ and is given by d max = τ(s max − CS)e −1 . This parameter (d max ) was used as an estimate of d for the three-parameter exponential model when comparing the d provided by the different models.

Data Analysis
Two different fitting procedures were used on the data set obtained for each participant to estimate CS, d , and s max . More specifically, t(s) and d(s) using WLS were evaluated. These two fitting procedures are statistically appropriate, i.e., they minimize the error along the axes corresponding to the dependent variables (Vinetti et al., 2020) and should overcome heteroscedasticity Morton and Hodgson (1996). Weights were applied to the corresponding dependent variables, i.e., time to exhaustion in t(s), and distance in d(s). Following Morton and Hodgson (1996), weights were set proportional to the inverse of the variance of the dependent variable, where the variance was itself set proportional to the dependent variable. Noteworthy, the model variants d(t) and t(d) have not been used. The reason being that in these cases, distance and time to exhaustion should be considered as dependent variables. However, the errors of both variables are correlated, i.e., the error of distance is given by the product of speed and the error of time to exhaustion variable, since speed does not carry any error. This is known as endogeneity and, to the best of our knowledge, there exists no regression method that can handle such case (Antonakis et al., 2014). Error minimization was performed iteratively using the Levenberg-Marquardt algorithm (Levenberg, 1944;Marquardt, 1963). After inspecting residual plots, deviations from homoscedasticity were present for the two fitting procedures applied to the three different models, the three-parameter model with d(s) and the two-parameter model with t(s) showing the least and the most heteroscedasticity, respectively (Supplementary Figure 1).
To obtain theVO 2 values at the CS estimates for each participant, first the CS estimates were converted to the times at which these running speeds occurred during the maximal incremental aerobic test assuming a linear relation between running speed and time [i.e., s = 2.78 + 0.14t, leading to t = (s − 2.78)/0.14, where t and s stand for time and running speed, respectively]. Then, theVO 2 values at the CS estimates were simply given by placing these corresponding times into the computed linear regression ofVO 2 as a function of time recorded during the maximal incremental aerobic test. Data analysis was performed using Python (version 3.7.4, Python Software Foundation 1 ).

Statistical Analysis
Descriptive statistics were expressed as the mean ± standard deviation. The 90% confidence intervals (CI) of CS, d and if applicable, s max , the combined standard error of the estimate (%SEE), i.e., the sum of SEE preliminary transformed to percent units of CS, d and if applicable, s max , and the AIC of the fitting procedure were computed to assess the quality of the fit. For the linear regression ofVO 2 as a function of time, its coefficient of determination (R 2 ) was calculated to examine its accuracy.
After inspecting residual plots, no obvious deviations from homoscedasticity and normality were present. Linear mixed models fitted by restricted maximum likelihood were used to compare CS, d , and s max obtained from the three mathematical models (two-parameter if applicable, three-parameter, and threeparameter exponential) and two fitting procedures [t(s) and d(s)]. The fixed effects included the mathematical models, fitting procedures, and their interaction. The within-subject nature was controlled for by including random effects for participants. The variance explained by the fixed effects over the total expected variance was given by R 2 marginal while R 2 conditional represented the variance explained by the fixed and random effects together over the total variance (Johnson, 2014). Intraclass correlation coefficients (ICC) of the random effects were computed as the ratios of the variance of the random coefficient divided by the sum of itself and the residual variance. On the basis of commonly used thresholds, poor, moderate, good, and excellent ICCs are given by ICC values <0.5, 0.5-0.75, 0.75-0.90, and ≥0.90, respectively (Koo and Li, 2016). Pairwise post hoc comparisons of any significant fixed effects were performed using Holm corrections.
Correlations, 90% CI, SEE (in %), and systematic differences of predicted value ( , in %) were computed among the three mathematical models and two fitting procedures with regard to CS, d , and s max and similarly between CS and aerobic fitness parameters. Data were log transformed as suggested by Hopkins et al. (2009). Correlations were computed using Pearson's correlation coefficients (r). Very high, high, moderate, low, and negligible correlations were given by r values of 0.90-1.00, 0.70-0.90, 0.50-0.70, 0.30-0.50, and 0.00-0.30, respectively (Hinkle et al., 2003). Statistical analysis was performed using Python, Jamovi (version 1.0.8, [Computer Software] 2 ), and R 3.5.0 (The R Foundation for Statistical Computing, Vienna, Austria) with a level of significance set at P ≤ 0.05.
The regression analysis for one representative participant and for each of the three mathematical models as well as the two fitting procedures [t(s) and d(s)] is presented in Figure 1. Table 1 depicts the time to exhaustion corresponding to the four exhaustive runs performed at 90, 100, 110, and 120% of the participant's PS. Table 2 depicts CS, d , and s max , together with their corresponding 90% CI, %SEE, and AIC obtained from the three mathematical models and two fitting procedures.
The linear mixed model with random effects explained almost all variance in the data for CS while a large part of variance in the data was still unexplained for d and s max even with random effects ( Table 3). These results were reinforced by the ICC of the random effects, which was excellent for CS but poor and moderate for d and s max , respectively ( Table 3). A significant mathematical model effect was obtained for CS, d , and s max (P < 0.001; Table 3). CS was significantly faster for the three-parameter exponential model compared with CS determined by two-(P < 0.001) and three-parameter (P < 0.001) models and it was significantly faster for the twothan for the three-parameter model (P < 0.001; Table 2). d was significantly lower for the two-and three-parameter exponential model than for the three-parameter model (P < 0.001; Table 2). The three-parameter exponential model had a significant slower estimation of s max than the three-parameter model (P < 0.001; Table 2).
No significant fitting procedure effect or significant mathematical model x fitting procedure interaction effect were reported for CS, d , and s max (P ≥ 0.77; Table 3).
On a group level, the average AIC was lower for the threeparameter model for both fitting procedures; however, it was very close to the average AIC for the three-parameter exponential model ( Table 2). Note that, because the units of the residual sum of squares error (RSS) depend on the fitting procedure itself, the AICs can be compared between models within a given fitting procedure but not between the two fitting procedures. On an individual level, t(s) and d(s) fitting procedures gave the lowest AIC when using the three-parameter model for 12 participants while 4 participants obtained the lowest AIC when using the three-parameter exponential model.
The three-parameter model reported the highest 90% CI as well as the highest %SEE (Table 2). However, %SEE can in general not be compared between the two-and three-parameter models because they do not have the same number of parameters to estimate. Nevertheless, the 90% CI of CS and d in the threeparameter and three-parameter exponential models were higher than in the two-parameter model, even if expressed in percent units. Therefore, the two models with three parameters carried more error on their estimates than the two-parameter model. The 90% CI and %SEE were similar when comparing the two fitting procedures for a given model ( Table 2).
TheVO 2 at the CS estimates expressed as a percentage ofVO 2 max as well as the CS expressed as a percentage of sV O 2 max for the three mathematical models and two fitting procedures are given in Table 4. TheVO 2 corresponding to the CS estimates were based on linear regression, therefore, the significant differences betweenVO 2 values were the same as those for the CV estimates ( Table 2; Housh et al., 2001).
Correlations, 90% CI, SEE, and between CS and aerobic fitness parameters are given in Table 5. Correlations were high and very high, and all statistically significant (P ≤ 0.001).

DISCUSSION
Conventional statistical approaches demonstrated a significant effect of the mathematical model when estimating CS, d , and s max , but no significant effect of the fitting procedure. These results validated our first hypothesis that the estimates of CS, d , and s max would be significantly different between mathematical models employed, but not between fitting procedures. Moreover, the three-parameter model gave the lowest estimation of CS, in accordance with our first hypothesis. Lower SEE and higher correlations between aerobic fitness parameters and CS estimated using a given mathematical model and fitting procedure were not necessarily associated with a lower AIC for these models and procedures, which refuted our second hypothesis. The linear mixed model showed interindividual differences in CS, d , and s max , as depicted by the larger R 2 conditional than R 2 marginal ( Table 3), but with a higher impact for CS than for d and s max , as depicted by the excellent ICC of the random effects for CS but poor and moderate ICCs for d and s max , respectively ( Table 3). In addition, a large part of the variance was still unexplained for d  The variance explained by the fixed effects over the total expected variance was given by R 2 marginal while R 2 conditional represented the variance explained by the fixed and random effects together over the total variance. Statistical significances (P ≤ 0.05) are indicated in bold.  [VO 2 ; expressed as a percentage of maximal rate of oxygen uptake (VO 2 max)] at the critical speed (CS) estimates as well as CS [expressed as a percentage of speed associated withVO 2 max (sV O 2 max )] for the three mathematical models [two-parameter (2-param), three-parameter (3-param), and three-parameter exponential (3-param exp)] and the two fitting procedures [t(s) and d(s) using weighted least squares]. and s max (R 2 conditional ≤ 72.0%; Table 3). Therefore, CS could be well estimated by using the mathematical model and fitting procedure, but this is not the case for d and s max . Furthermore, the high and very high betweenmodel correlations (r ≥ 0.93) obtained for CS suggest that the estimation of CS provided by each model qualitatively represents the same, as already pointed out by Gaesser et al. (1995). By contrast, some between-model correlations were high and moderate for d and s max , respectively (r ≥ 0.67), suggesting less link between estimations of d and s max from the different regression analyses. Overall, the regression analyses provided more robust estimates of CS than of d and s max .
The three-parameter model gave the lowest AIC on a group level as well as for 75% of the participants for both t(s) and d(s) fitting procedures. Nevertheless, the AICs of both threeparameter models were very close to one another ( Table 2). The AICs for the two-parameter model were 34% and 21% higher than the three-parameter model ones for t(s) and d(s), respectively. Therefore, the two-parameter model gave the lowest quality of the fit, while the three-parameter model seemed to be the most accurate one for both fitting procedures even though it was only slightly better than the three-parameter exponential model. These observations contradict previous findings that obtained similar R 2 values between different mathematical models (Gaesser et al., 1995;Bull et al., 2000;Housh et al., 2001;Bergstrom et al., 2014) [except for a two-parameter linear model expressing power as function of 1/time (Gaesser et al., 1995;Bull et al., 2000)]; this might be explained several ways. First, comparing the accuracy of regression analyses for models based on a different number of parameters (e.g., two vs. three parameters) requires an adjusted R 2 to normalize with respect to the number of parameters within the model. However, these studies (Gaesser et al., 1995;Bull et al., 2000;Housh et al., 2001;Bergstrom et al., 2014) did not mention such usage. Second, R 2 was shown to be an unfavorable measure to describe the validity of a non-linear regression (e.g., both model formulations of the three-parameter model) (Spiess and Neumeyer, 2010) and when using weights in the regression analysis (Willet and Singer, 1988). Therefore, one remaining possibility to compare the quality of the fit of different mathematical models is to use RSS or a parameter that depends on it such as AIC. However, the units of RSS (and thus AIC) being dependent on the fitting procedure (i.e., on the choice of model formulation and axes on which the errors are minimized), the AICs of the various fitting procedures cannot be compared, i.e., AIC of t(s) cannot be compared to the one of d(s).
Another option to compare the quality of the fit of different mathematical models is to use %SEE (Triska et al., 2021). However, as already mentioned, %SEE depends on the number of parameters to estimate and is therefore not optimal to compare two-and three-parameter models. Nevertheless, in our case, the 90% CI of CS and d in the three-parameter and three-parameter exponential models were higher than in the two-parameter model, even if expressed in percent units. Therefore, the two-parameter model gave the lowest %SEE (9.7%) and SEE for CS (0.7%) and d (9%), while the threeparameter gave the highest %SEE (25%; CS: 2.2%, d :17%, and s max : 5.8%), but only 1% higher than the three-parameter exponential model (24%; CS: 0.8%, d :20.8%, and s max :2.4%). Based on %SEE, the three-parameter model seemed to be the least accurate model, which is in contradiction with the results based on AIC.
CS was thought to reflect an inherent characteristic of the aerobic energy supply system (Hughson et al., 1984;Gaesser and Wilson, 1988;Poole et al., 1988). Such a characteristic is supported by the small SEE and high and very high correlations obtained between CS and aerobic fitness parameters such as VT, RCP, andVO 2 max (SEE ≤ 6.84; r ≥ 0.77; Table 5). These results additionally confirm previous observations that showed that CS correlated withVO 2 max (Hughson et al., 1984;Gaesser and Wilson, 1988;Poole et al., 1988) and RCP (Moritani et al., 1981). However, the three-parameter model reported the highest SEE [if we do not consider SEE for the two-parameter model and d(s)] and smallest correlations, which were associated with the largest 90% CI (4.43 ≤ SEE ≤ 5.82; 0.77 ≤ r ≤ 0.85; Table 5). This is in line with the fact that the three-parameter model reported the highest %SEE (Table 2). Nonetheless, SEE and correlations were still small and high, respectively, for this model. The linear mixed model provided a significant effect of the mathematical model when estimating CS, d , and s max ( Table 3). These results accord with those of previous observations that depicted considerable differences in the estimation of parameters among different models (Gaesser et al., 1995;Bull et al., 2000;Housh et al., 2001;Bergstrom et al., 2014). The three-parameter model provided the lowest estimation of CS on a group level ( Table 2) as well as on an individual level. CS estimated using the three-parameter model were 6% and 9% smaller than when using the two-parameter and three-parameter exponential models, respectively. The two-parameter model was shown to produce overestimated CS (Pepper et al., 1992). The authors observed that the time to exhaustion at a running speed set at CS estimated by the two-parameter model was much smaller than expected. Indeed, participants were able to run only 16 min instead of a theoretically indefinite time. Because CS predicted by the three-parameter exponential model was faster than CS predicted by the two-parameter model (+3% ; Table 2), we could conclude that the three-parameter exponential model also produced overestimated CS. The observed between model differences for the CS estimates (up to 0.44 m/s, Table 2) are not negligible and would certainly have an impact when prescribing a training session based on exercise intensity. Therefore, we encourage coaches prescribing exercise based on critical intensity to choose a mathematical model beforehand to estimate CS and maintain it over the running seasons, so that CS is always estimated in the same way. Moreover, even though the estimated CS should be a very good approximation of the critical intensity but not the critical intensity per se, we suggest to physiologically verify that the estimated CS represents the upper boundary of sustainable exercise. In addition, coaches should not hesitate to make small adjustments based on the observed performance. Moreover, given the day-to-day variation of human performance and the CI of the estimated CS, i.e., about 5% of its value ( Table 2), it would be justified to prescribe exercise intensity outside these confidence limits to avoid being in the phase transition between the heavy and severe intensity domains (Anderson et al., 2019). Jones and Vanhatalo (2017) found that CS occurred at 70-90% ofVO 2 max, depending on training status (the higher the training status, the higher the CS in %VO 2 max). In the present study, theVO 2 at the CS estimates for the threeparameter model were close to the middle of the range defined by Jones and Vanhatalo (2017) (83%; Table 4), while thė VO 2 at the CS estimates for the two-parameter and threeparameter exponential models were in the higher end of the range (≥88.2 %VO 2 max; Table 4). HigherVO 2 at the CS estimates were already reported by Housh et al. (2001) for the two-parameter and three-parameter exponential models (≥94 %VO 2 max) than for the three-parameter model (89 %VO 2 max). These authors even reportedVO 2 at the CS estimates that exceededVO 2 max for the exponential model (105 %VO 2 max). In the present study, CS corresponded to 81, 87, and 90 %sV O 2 max for the three-parameter, two-parameter, and threeparameter exponential models, respectively. Billat et al. (1995) observed that CS corresponded to 86% of sV O 2 max for runners having 75 ml/min/kg ofVO 2 max and 6.22 m/s of sV O 2 max . These values were higher than those of the participants of this study (+16 and +19%, respectively). Therefore, we could speculate that CS estimated by the three-parameter model (81 %sV O 2 max ) is closer to reality than CS estimated by the other two models (≥87 %sV O 2 max ). Both arguments reinforce the idea that both two-parameter and three-parameter exponential models overestimate CS. In any case, a future study involving exhaustive runs below, at, and above CS whilst assessing oxygen uptake responses to exercise would be needed to quantitatively validate this suggestion.
The estimation of d using the three-parameter model were roughly 2.5 times larger than those from the other two models. These findings are consistent with those of previous studies (Gaesser et al., 1995;Morton, 1996;Bull et al., 2000;Housh et al., 2001;Bergstrom et al., 2014). Morton (1996) suggested that such a model overcomes physiological assumptions of the two-parameter model such as an infinite power when time approaches zero and that at d , the muscular energy reserve is empty. Assuming an sV O 2 max of 6 m/s and a time to exhaustion of ∼5 min at 100% of sV O 2 max (Billat et al., 1995), the corresponding total distance covered is 1,800 m. The anaerobic contribution was shown to represent approximately 10% of the total distance covered, i.e., approximately 200 m (Billat, 2001). Therefore, because Morton (1996) suggested that the three-parameter model allows d to be only partly covered for a running speed between CS and s max , this statement causes the estimate of d that is larger in the three-parameter model than in the other two models to not be unrealistically high. This idea is reinforced by an explanation based on anaerobic energy calculation by Gaesser et al. (1995). Buchheit and Laursen (2013) found that athletes with similar sV O 2 max to those of the present study had a s max ranging from 161 to 183 %sV O 2 max . Higher level athletes (sV O 2 max = 6.36 m/s) were shown to have a lower relative s max (149 %sV O 2 max ) (Sandford and Stellingwerff, 2019). Therefore, the estimation of s max using the three-parameter exponential model seemed to be unrealistically too small (∼136 %sV O 2 max ) whereas the one obtained using the three-parameter model seemed closer to reality (∼155 %sV O 2 max ). Nonetheless, this has to be nuanced by the fact that most of the running trials at 120 %PS gave a time to exhaustion below 2 min, which is below the usual recommendation (Jones and Vanhatalo, 2017) and could have influenced the estimation of the parameters present in the mathematical models. In addition, participants were long distance runners, meaning that they are not accustomed to running at high speeds (i.e., >100 %sV O 2 max ) and that they actually did not have a high s max . This assumption is supported by the observations of Sandford and Stellingwerff (2019), who showed that a 400-m elite runner (sV O 2 max = 6.23 m/s) had a s max of 158 %sV O 2 max while a 1,500-m elite runner (sV O 2 max = 6.45 m/s) had a s max of 141 %sV O 2 max . No significant fitting procedure or mathematical model x fitting procedure interaction effects were reported for the estimations of CS, d , and s max ( Table 3). Gaesser et al. (1995) proposed that differences between the estimation of parameters among models could come from the designation of the dependent and independent variables, the number of parameters in each model, and the choice of model (e.g., two-parameter, threeparameter, or three-parameter exponential). Moreover, two mathematically equivalent model formulations requiring linear vs. non-linear regressions were shown to provide different estimations of their underlying parameters (Colquhoun, 1971). In this study, we observed that using different but statistically appropriate fitting procedures, i.e., that correctly attribute the dependent and independent variables, applied to a given model did not have an impact on the estimations of CS, d , and s max , as long as all the model formulations are non-linear or linear.
Heteroscedasticity of the dependent variable was explicitly depicted by Hinckson and Hopkins (2005) when using usual LS fitting procedure. Indeed, these authors demonstrated systematic and non-uniform deviation from their models by showing the residuals as function of predicted values. In this study, the suggestion made by Morton and Hodgson (1996) to overcome heteroscedasticity, i.e., weights proportional to the inverse of the values of the dependent variable, were applied. However, the absolute weighted residuals as function of predicted values for the two fitting procedures applied to the three different models depicted clear deviations from homoscedasticity (Supplementary Figure 1). Therefore, considering weights in the fitting procedure did not overcome the heteroscedasticity problem. Nonetheless, a future study considering different weighting schemes should be performed in order to observe if a specific weighting scheme, different from the one proposed by Morton and Hodgson (1996), could overcome heteroscedasticity of the dependent variable.
Some limitations to the present study exist and need to be addressed. On the one hand, the participant should complete five experimental sessions interspersed by at least 2 days, which could be slightly unpractical. On the other hand, performing a regression analysis with only four measurement points is already quite few, especially when dealing with heteroscedasticity. Nonetheless, the estimation of CS based on four points is considered as the best practice (Poole et al., 2021). Moreover, there is a well-known large variability in the time to exhaustion during treadmill running at CS (Pepper et al., 1992). Furthermore, due to the proximity between CS and RCP in terms of %VO 2 max (CS: 87.6 %VO 2 max; RCP: 89.3 %VO 2 max) and the high and very high correlations between them (r ≥ 0.85), one could wonder the relevance of CS. However, the recent meta-analysis of Galán-Rioja et al. (2020) showed that CS and RCP are not synonymous. Besides, CS can be estimated using personal best times, which does not require the participant to go to the laboratory (Jones et al., 2019). Finally, a recent study demonstrated that using estimations of CS from raw training data can be sufficient to successfully predict marathon performance and provide useful pacing information (Smyth and Muniz-Pumares, 2020).
To conclude, this study demonstrated that CS, d , and s max estimated from three different mathematical models (twoparameter, three-parameter, and three-parameter exponential model) differed significantly, but that no difference in the estimation of CS, d , and s max was reported between different statistically appropriate fitting procedures applied to a given model. Weights did not help overcoming heteroscedasticity of the dependent variable. CS estimates from the three different models were correlated with aerobic fitness parameters, i.e., VT, RCP, andVO 2 max. Moreover, small SEE was obtained. The three-parameter model gave the lowest AIC on a group level and the smallest CS estimates. However, the three-parameter model reported the highest %SEE and 90% CI. Therefore, our results showed no further support for selecting the best mathematical model to estimate critical speed. Nevertheless, our results showed that statistically appropriate fitting procedures gave the same estimates for a given model. For these reasons, we suggest coaches choosing a mathematical model with appropriate fitting procedure beforehand to define CS and intensity domains and maintaining it over the running seasons. Moreover, our findings suggest that each CS estimation during season should be physiologically verified and training prescription should be done around CS (±10%) for taking into account CI of its estimation and the day-to-day variation of human performance.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.