## Abstract

Modeling in the time domain, the non-steady-state O_{2} uptake on-kinetics of high-intensity exercises with empirical models is commonly performed with gradient-descent-based methods. However, these procedures may impair the confidence of the parameter estimation when the modeling functions are not continuously differentiable and when the estimation corresponds to an ill-posed problem. To cope with these problems, an implementation of simulated annealing (SA) methods was compared with the GRG_{2} algorithm (a gradient-descent method known for its robustness). Forty simulated V̇o_{2} on-responses were generated to mimic the real time course for transitions from light- to high-intensity exercises, with a signal-to-noise ratio equal to 20 dB. They were modeled twice with a discontinuous double-exponential function using both estimation methods. GRG_{2} significantly biased two estimated kinetic parameters of the first exponential (the time delay td_{1} and the time constant τ_{1}) and impaired the precision (i.e., standard deviation) of the baseline *A*_{0}, td_{1}, and τ_{1} compared with SA. SA significantly improved the precision of the three parameters of the second exponential (the asymptotic increment *A*_{2}, the time delay td_{2}, and the time constant τ_{2}). Nevertheless, td_{2} was significantly biased by both procedures, and the large confidence intervals of the whole second component parameters limit their interpretation. To compare both algorithms on experimental data, 26 subjects each performed two transitions from 80 W to 80% maximal O_{2} uptake on a cycle ergometer and O_{2} uptake was measured breath by breath. More than 88% of the kinetic parameter estimations done with the SA algorithm produced the lowest residual sum of squares between the experimental data points and the model. Repeatability coefficients were better with GRG_{2} for *A*_{1} although better with SA for *A*_{2} and τ_{2}. Our results demonstrate that the implementation of SA improves significantly the estimation of most of these kinetic parameters, but a large inaccuracy remains in estimating the parameter values of the second exponential.

- O
_{2}uptake kinetics - parametric modeling
- parameter estimation

for several decades, the study of O_{2} uptake (V̇o_{2}) kinetics at the onset of step exercises has attracted the attention of physiologists. With the introduction of automatic systems for computing V̇o_{2} on a breath-by-breath basis (1, 4), there has been an increasing interest in characterizing the kinetics of V̇o_{2} in the time domain. As the breath-by-breath analysis allows a high time resolution, the non-steady-state V̇o_{2} kinetics have been currently modeled with different mathematical approaches (32). Two different types of mathematical models can be derived to characterize V̇o_{2} kinetics: structural and empirical models. “A structural model consists of a series of mathematical relationships that have identifiable physiological correlates and describe the intermediate components of the physiological response to an input stress. Integration of these intermediate responses leads to a prediction of the output response of interest, which then can be compared with observed responses.” Inversely, in empirical models, “the best description of the output response is chosen by some (usually statistical) criteria from a set of competing mathematical models. Insight into the nature of the underlying intermediate physiological processes is sought from analysis of the model components and parameters” (2). Then, the strategy for fitting empirical models to experimental data differs from fitting structural models in that the simplest model that describes the data satisfactorily is sought. This implies that no physiological signification can a priori be assigned to the parameters of these empirical models, albeit physiological correspondence may be a posteriori hypothesized.

In the context of empirical models applied to V̇o_{2} kinetics to step increases in power, the estimation of the kinetic parameters has been commonly performed by nonlinear regression algorithms. These methods typically refer to gradient-descent-based (GD) algorithms (6, 7, 11, 17, 35), validated on problems with continuously differentiable functions and in the absence of numerical difficulties such as degeneracy or ill conditioning (i.e., slight fluctuations in the observed data can result in very large fluctuations in the estimated parameters) (14, 20, 21). For exercise transitions from light to moderate work rates [<60% maximal V̇o_{2} (V̇o_{2 max})], the V̇o_{2} time course is assumed to be well described by a first-order exponential function (2, 19, 26, 29, 35). Then, GD methods seem to be appropriate, because the model is well conditioned and continuously differentiable over the period of work. For exercise transitions from light to higher work rates (>60% V̇o_{2 max}), the estimation of the kinetic parameters refers to a more complex mathematical problem because the V̇o_{2} kinetics are generally best described with a second-order exponential model incorporating unit step functions (2, 3, 6, 7, 25). The second exponential is introduced in the model to describe a drift (rather than a steady state) in V̇o_{2} that is observed beyond the third minute of exercise. This model is not continuously differentiable owing to the introduction of step functions. And the problem of fitting a sum of exponentials to numerical data is known to be an ill-conditioned optimization problem (24). In this context, GD methods do not seem perfectly adapted to the mathematical properties of this double exponential function and may impair the estimation of the kinetic parameters. An additional drawback is that GD methods need to be initialized by the investigator (6, 7, 11, 17): the complex properties of the double exponential model may compromise or influence the convergence of the estimation when the investigator introduces starting values far from the optimal solution. An alternative approach to solve this fitting problem consists in using stochastic methods of optimization such like simulated annealing (SA), which belongs to Monte Carlo Markov chain algorithms (15, 24): they are designed for locating the global optimum of the function to be optimized regardless of its mathematical properties. They are also unsupervised in the sense that there is no need for intervention of the investigator because the starting values are randomly chosen in a parameter space.

The purpose of the present paper is to propose an implementation of SA specifically dedicated to fit second-order exponential models with step functions to breath-by-breath V̇o_{2} data, in exercise transitions from light- to high-intensity work rates. Applications on simulated and experimental data are performed and compared with the results obtained with a supervised GD algorithm [the generalized reduced gradient 2 algorithm (GRG_{2})] already used in previous studies (6, 7). The choice for this last algorithm has been dictated by the fact that it has a reputation of robustness on problems involving numerical difficulties, compared with other nonlinear optimization methods (14).

## METHODS

### Empirical Model to Describe the V̇o_{2} Kinetics

To compare SA and GRG_{2}, 40 simulated V̇o_{2} responses and 52 experimental V̇o_{2} responses are fitted by both algorithms, with a second-order exponential function of the form (1) where is the unit step function. V̇o_{2} is the rate of O_{2} uptake at time *t*. *A*_{1} and *A*_{2} are the asymptotic increments above the baseline *A*_{0}, τ_{1} and τ_{2} are the time constants, and td_{1} and td_{2} are the time delays for each exponential term after the start of the exercise. This function is commonly used to describe V̇o_{2} on-responses from light- to high-intensity (>60% V̇o_{2 max}) work rates (2, 3, 6, 7, 25).

When fitting the experimental data with this model, previous studies usually excluded the first 30 s after the start of exercise (19, 35). The abrupt increase in V̇o_{2} occurring during this first period reflects cardiodynamic gas exchange. The rationale for its exclusion is that it is functionally distinct from the subsequent exponential response (33, 35). Its magnitude is also recognized to decrease dramatically when the step transition is from a light-to-moderate exercise intensity compared with a rest-to-exercise transition (35). And the higher the baseline of the transition, the less distinguishable the cardiodynamic component will be from the subsequent exponential V̇o_{2} increase (18). In the present study, the V̇o_{2} baseline is fixed at ≈1.5 l/min (i.e., a light-intensity exercise). Then, the cardiodynamic component is considered to be negligible, and all V̇o_{2} data measured during the exercise transition are included in the curve fitting.

### Estimation Algorithms

The mathematical concepts of both estimation methods (SA and GRG_{2}) and their implementations are extensively documented in the appendix. In the present study, they are used to fit V̇o_{2} on-responses acquired breath by breath. Owing to breathing irregularities, these V̇o_{2} data fluctuate from breath to breath and so can be viewed as the sum of an underlying physiological information and “noise.” Because this noise is mainly white and Gaussian (19, 27, 30), the criterion to be optimized by SA and GRG_{2} is the residual sum of squares (RSS) between the raw data and the model (RSS has to be minimized). In a physiological sense, the configuration of both algorithms imposes some constraints on the estimated parameters: *1*) *A*_{0}, *A*_{1}, and *A*_{2} have real positive values, *2*) td_{1} ≤ td_{2}, and *3*) 0 < τ_{1} ≤ τ_{2}.

Furthermore, the initialization of GRG_{2} is made with the same starting values of kinetic parameters, whatever the response to fit. However, to minimize the influence of initialization on the “optimum” solution, each V̇o_{2} response is fitted three times with three different sets of starting values (Table 1). Over the three “optimum” solutions obtained for a given response, the solution retained is the one producing the lowest RSS.

### Comparison of Both Methods on Simulated Data

#### Simulated data.

Simulated V̇o_{2} responses are generated on the basis of real time course and temporal responses. First, experimental V̇o_{2} data have been measured during a period of 2 min of a light exercise followed by 8 min of constant-load exercise (80% V̇o_{2 max}) on two different subjects (see *Experimental data*). Then, the second-order exponential model (*Eq. 1*) has been used to fit the V̇o_{2} responses of both subjects to generate two different sets of kinetic parameters for the transition from light to high intensity (Table 2). At high intensity (>60% V̇o_{2 max}), it is recognized that the V̇o_{2} time course described by the second component of *Eq. 1* (and called the “slow V̇o_{2 max} component”) either delays the attainment of a steady state or leads V̇o_{2} to increase inexorably until the end of the exercise (or until V̇o_{2} is reached) (2, 3, 6, 7). Therefore, both subjects have been carefully chosen to take into account both shapes of time course. *Set 1* of kinetic parameters produces a second exponential component during which V̇o_{2} reaches the asymptotic value of the component before the end of the exercise. *Set 2* produces a second exponential component during which V̇o_{2} increases almost linearly until the exercise stops. Therefore, both reference V̇o_{2} kinetics have been noised by generating 20 Gaussian white noise vectors with a variance equal to 52,000 for *set 1* and 20 Gaussian white noise vectors with a variance equal to 125,000 for *set 2*. As a whole, this produced 40 different simulated V̇o_{2} on-responses with an approximately constant SNR (SNR_{dB} in decibels) equal to 20 dB. Figure 1 shows examples of simulated V̇o_{2} on-responses obtained with *set 1* (*A*) and *set 2* (*B*) of kinetic parameters (the dark lines indicate the underlying second-order exponential curves).

#### Statistical analysis.

Each of the 40 simulated responses is fitted with the second-order exponential function (*Eq. 1*) by using GRG_{2} and the SA algorithms with the configurations presented above. It is unlikely that an estimation method produces estimates equal to the reference value. But the most “optimum” estimation method produces estimates that are the closest to the reference value for a given set of noise statistics. This is evaluated on the basis of two criteria: *1*) The bias (*b*) is an index of accuracy of the estimation. It is the mean of the estimation errors, i.e., the mean difference between the reference value and the estimates. *2*) The standard deviation of the estimation errors (σ) is an index of the precision of the estimation.

For each kinetic parameter, results produced by each estimation procedure are expressed as the bias (*b*), its 95% confidence interval, and the standard deviation of the estimation error (σ). A Student's *t*-test is used to analyze whether each bias is significant (i.e., whether *b* is significantly different from 0). Comparison of the precision of an estimated parameter is done with an *F*-test on the variances (σ^{2}) of this parameter produced by both estimation procedures. Comparison of biases is done with a paired Student's *t*-test. The critical level accepted for statistical significance is set at *P* = 0.05.

### Comparison of Both Methods on Experimental Data

#### Overview.

The purpose of this section is to verify that the results obtained on simulated data are in accordance with those obtained on experimental data. As presented above, SA and the GRG_{2} algorithms determine the kinetic parameters by minimizing RSS. The method that provides the lowest RSS could be considered as the better one (see the appendix). Therefore, the first criterion used to compare both algorithms on experimental data is RSS obtained on the same signals.

In parallel, it is interesting to compare the repeatability of each estimation procedure. Several repetitions of the same exercise by humans offer an intraindividual variability of the response due, for one part, to the variability of the underlying time course of V̇o_{2} and, for another part, to the variability of the noise. We make the assumption that, in healthy subjects, the intraindividual variability of the underlying physiological response is low. Therefore, the most pertinent algorithm should be the one that is the less sensible to the intraindividual variability of the noise; this is the estimation method that produces the lowest difference of estimation between *test 1* and *test 2*. The second criterion of comparison is therefore the coefficient of repeatability adopted by the British Standards Institution (8).

#### Subjects.

Twenty-six healthy young men, students in physical education, volunteered to participate in this study after being informed of the purpose, the requirements, and their rights as subjects. They were (means ± SD) 22 ± 1 yr of age, their height was 182 ± 8 cm, and their body mass was 73.0 ± 7.9 kg. The experimental protocol was approved by the Ethics Committee of the Faculty of Medicine of the Université Catholique de Louvain, and all subjects gave informed, written consent before participating in the study. The subjects were fully familiar with laboratory exercise.

#### Experimental protocol.

The subjects were instructed to arrive at the laboratory in a rested and normally hydrated state, at least 2 h postprandial, and to avoid strenuous exercise in the 24 h preceding a test session. For each subject, tests took place at the same time of day. All exercises were performed on a friction-loaded cycle ergometer (Monark type 818E, Stockholm, Sweden). For all tests, seat height was held constant within a subject and the feet were strapped to the pedals. In a preliminary test session, the subjects performed a ramp exercise test (20 W/min; pedal rate 80 rpm) to volitional exhaustion, preceded by a 5-min unloaded cycling work. The breath-by-breath data were collected and displayed at 30-s intervals to determine peak V̇o_{2}. During each subsequent session, subjects performed a square-wave exercise transition from a 4-min baseline at 80 W to an 8-min heavy exercise at 80% of peak V̇o_{2} (pedal rate 80 rpm). Pulmonary gas exchange was measured breath by breath during the entire duration of the square-wave exercise. This test was realized twice in separate sessions to test the repeatability of the kinetic parameters.

#### Pulmonary gas-exchange analysis.

The subjects breathed through a face mask during the entire duration of the square-wave exercise, and pulmonary gas exchanges were monitored breath by breath by use of an integrated computerized system (Quark Cosmed, Roma, Italy). Precision-analyzed gas mixtures and a 3-liter Rudolph syringe were used to calibrate the rapid gas analyzers and the pneumotachograph. Then, for each test within a subject, the breath-by-breath responses of V̇o_{2} were time aligned to the start of the heavy work.

#### Statistical analysis.

Each of the 52 experimental responses (26 for the first test and 26 for the second test) is fitted with the second-order exponential function (*Eq. 1*) by using GRG_{2} and the SA algorithms with the configurations presented above. For RSS analysis, a linear regression analysis is computed to compare the RSS values produced by both algorithms on the 52 experimental responses. For the analysis of the repeatability of each kinetic parameter, each algorithm produces the mean of the 26 intraindividual differences of repeated measurements and the corresponding coefficient of repeatability. This coefficient of repeatability is defined as twice the standard deviation of the 26 intraindividual differences of repeated measurements (8). The prerequired condition for repeatability analysis is that the mean of the intraindividual differences be not significantly different from zero. For each estimation procedure, the mean difference for each kinetic parameter is compared with zero by a Student's *t*-test. Comparison of the repeatability of both estimation procedures is done with an *F*-test used on the variances of the intraindividual differences. The critical level accepted for statistical significance is set at *P* = 0.05.

## RESULTS

### Simulated Data

For each parameter of the second-order exponential function, Fig. 2 and Table 3 present the bias, the limits for its 95% confidence interval, and the standard deviation of the estimation error produced either by GRG_{2} or by the SA algorithm on the whole 40 simulated responses.

Concerning the baseline and the first exponential component of the model (called phase II), the results demonstrate that the SA algorithm produces better estimates than GRG_{2}. For three of four kinetic parameters (*A*_{0}, td_{1}, τ_{1}), the precision of the estimations is significantly higher with the SA than with the GRG_{2} algorithm (*P* < 0.05), and two of four parameters are significantly biased by the GRG_{2} algorithm (td_{1} and τ_{1}; *P* < 0.05). Moreover, the bias for td_{1} and τ_{1} is significantly lower with the SA than with the GRG_{2} algorithm (*P* < 0.05). The bias for *A*_{1} is significantly higher with the SA than with the GRG_{2} algorithm (*P* < 0.05), but its value remains small compared with the mean magnitude of phase II (<4%).

Concerning the second exponential component of the model (called the slow component), the precision is significantly higher with the SA algorithm for all three parameters (*A*_{2}, td_{2}, and τ_{2}; *P* < 0.05). Nevertheless, the high values for the confidence intervals of the biases and for the standard deviations of the estimation errors of the three kinetic parameters indicate that all estimations are impaired, even when produced by the SA algorithm. Finally, td_{2} is significantly biased by both algorithms (*P* < 0.05).

### Experimental Data

Figure 3 presents the linear relationship between RSS values (*r*^{2} = 0.9995; *P* < 0.05), with the values produced by the SA algorithm being the independent variables and the values produced by GRG_{2} being the dependent variables. The *y*-intercept (b) is significantly higher than 0 (b = 278,244; *P* < 0.05), whereas the slope of the relationship (*a*) is not significantly different from the unity (*a* = 0.9999; *P* > 0.05). In details, 46 of the 52 exercise tests show a lower RSS value when produced by the SA algorithm. On these 46 tests, the SA values of RSS are on average 4% lower than with GRG_{2}. Conversely, in the six other situations, the average value produced by the GRG_{2} algorithm is only 0.2% lower than the value produced by SA*.*

For each kinetic parameter, Table 4 presents the mean values (SD) obtained on the first and second repetition of the exercise, the mean intraindividual difference, and the repeatability coefficient for each estimation procedure. It must be emphasized that the estimated average SNR is ≈21 dB for each test (*test 1* and *test 2*), which is comparable to the SNR used to generate simulated data. Except for *A*_{0}, the repeated measurements (*test 1* vs. *test 2*) are not significantly different from each other, whatever the estimation procedure (*P* > 0.05). The repeatability coefficients are not significantly different between estimation procedures for td_{1}, τ_{1}, and td_{2} (*P* > 0.05). GRG_{2} gives a significantly better coefficient of repeatability for *A*_{1} (*P* < 0.05). The repeatability coefficients are significantly better with SA for *A*_{2} and τ_{2} (*P* < 0.05). As for simulated data, the large repeatability coefficients for the kinetic parameters of the second-exponential component show that both estimation procedures may produce large differences of parameter values from the same (or comparable) V̇o_{2} on-response.

## DISCUSSION

The main result of the present study is that the SA algorithm, compared with the GRG_{2} algorithm, improves the estimation of the kinetic parameters of the second-order exponential function describing the V̇o_{2} time course in a transition from light- to high-intensity exercises. This is shown by results on simulated data and also by the global better estimation of parameters on experimental data (SA returns lower RSS than GRG_{2}). To our knowledge, this study is the first to focus specifically on the improvement of the estimation procedure for these second-order kinetic parameters.

Comparison of both algorithms (SA vs. GRG_{2}) on simulated data is linked to the postulate that the intersample fluctuations of pulmonary gas exchanges are white and Gaussian. Previous studies have already characterized these fluctuations as being white and Gaussian in healthy and chronic obstructive pulmonary disease adults (19, 27, 30) but focused only on exercises of moderate intensity (for which the pulmonary gas exchanges are appropriately described by a first-order exponential function at the onset of a step exercise). In these studies, the restriction to this intensity domain facilitated the statistical analysis of the breath-to-breath fluctuations because it could be performed at V̇o_{2} steady state. In the present study (with a more complex time course of pulmonary gas exchanges than at moderate intensity), the same statistical analysis was difficult to perform because the attainment of a V̇o_{2} steady state was a priori not known. Therefore, we made the assumption that the statistical characteristics of the noise of both intensity domains (moderate and high) were comparable, and we generated simulated data using additive white Gaussian noise. This postulate could have then produced different results between simulated and experimental data if the statistical characteristics of both kinds of data had been dramatically different. However, the best estimations produced by the SA algorithm (compared with GRG_{2}) on simulated V̇o_{2} kinetics are confirmed on experimental data.

The superiority of SA may be explained by two main factors. The pertinence of the SA algorithm depends on the configuration of intrinsic parameters, such as the transition kernel (see the appendix). Our results suggest that the configuration chosen for these parameters is appropriate enough to improve the estimations made by the GRG_{2} algorithm. The second explaining factor is that the SA algorithm is a stochastic unsupervised method, whereas GRG_{2} needs to be initialized. Therefore, the estimated parameters found by GRG_{2} are supposed to be influenced by the initialization. It must be kept in mind that, in the present study, an attempt to minimize this drawback was made by testing several initializations. Despite this precaution, the SA algorithm produces better estimates than GRG_{2}.

Our results on simulated data also demonstrate that the confidence levels observed for the kinetic parameters of the first exponential component (td_{1} and τ_{1}) are in accordance with the results published elsewhere for the corresponding kinetic parameters of moderate-intensity step exercises, i.e., for kinetic parameters of phase II (19). With an average SNR equal to 20.3 dB, the SA algorithm produces a standard deviation of the estimation error equal to 11.4 and 9.0 s for td_{1} and τ_{1}, respectively. Lamarra et al. (19) have previously proposed a prediction equation of this standard deviation incorporating two factors: the SNR value of the fitted signal and a constant *L* that depends on the value of the reference time constant and on the amount of data available for fitting. This equation, validated for a first-order exponential model when fitted with a GD method, predicts standard deviations around 10.5 and 8.8 s for td_{1} and τ_{1}, respectively, when taking into account our average SNR_{dB} of 20.3 dB (Lamarra et al. did not directly express the SNR in dB but used an index of the noise-to-signal ratio defined as the ratio between the standard deviation of the fluctuations and the magnitude of the response; in the present case, this gives an average value of ≈21%). This is an interesting result because it is well known that the parameter confidences of a model degrade with increasing the model order. Because we focused on a second-order exponential function describing V̇o_{2} time course, this quality of estimation was not expected compared with the estimations made by Lamarra et al. for a first-order function. This is in evidence an improvement due to the SA algorithm because the quality of estimation for td_{1} and τ_{1} is comparatively less with GRG_{2} (the standard deviations of the estimation error are 30.8 and 24.9 s, respectively). Concerning the experimental data, the repeatability coefficients of phase II parameters seem a priori larger than those previously published. Puente-Maestu et al. (26) have tested the repeatability of the time constant (τ) and of the magnitude of the asymptotic response (*A*) for the phase II of the V̇o_{2} time course of moderate-intensity step exercises. The noise-to-signal ratio [expressed as defined by Lamarra et al. (19)] was 10%. Using the equation of Lamarra et al., a prediction of their repeatability coefficient for the same SNR_{dB} than in our study (SNR_{dB} = 20.3 dB) gives values equal to 21.2 s and 220 ml/min for τ and *A*, respectively. Our values of repeatability coefficient for the corresponding kinetic parameters of phase II are found to be larger (at least equal to 34 s and 370 ml/min for τ_{1} and *A*_{1}, respectively). A possible explanation may be linked to the time constant value. Lamarra et al. have showed that the higher the time constant value, the higher the confidence level of a kinetic parameter for a given SNR_{dB}. Because the population studied by Puente-Maestu et al. exhibits time constant values close to 80 s (more than twice the values of the present population), this may partly explain the present worst result of repeatability.

Concerning the second exponential component of the model, the SA algorithm improves the confidence intervals of the biases and the standard deviations of the estimation errors of the three parameters (*A*_{2}, td_{2}, and τ_{2}) and also the repeatability of *A*_{2} and τ_{2} (*P* < 0.05 compared with the GRG_{2} algorithm). Unfortunately, these values remain high whatever the algorithm used. Therefore, their interpretation is experimentally limited. One possible explanation for this poor confidence level is linked to the shape of the underlying V̇o_{2} kinetics. On simulated data, *set 2* of reference kinetic parameters induces a linear increase of V̇o_{2} (rather than an exponential one), over the duration corresponding to the second exponential. This is due to both high values for the magnitude and the time constant of this exponential component (Table 2). In this context, the poor confidence level may be explained by the fact that the *A*_{2}-to-τ_{2} ratio is a good approximation of the slope of a linear function describing this increase in V̇o_{2}; thus infinity of values for both parameters may produce this ratio (2, 3). A second possible explanation is linked to the low SNR_{dB} value of the data. We simulated breath-by-breath V̇o_{2} time courses with a large breath-by-breath variability of the data (corresponding to a low SNR_{dB}) to mimic the typical breath-by-breath variability of V̇o_{2} data obtained with commercial gas-exchange analysis systems. Although this does not prevent doing a comparison of both algorithms (which is the purpose of the present study), it is well recognized that the lower SNR_{dB}, the lower the confidence level of the estimation of a model parameter. To improve further this confidence level, reducing the magnitude of the noise of V̇o_{2} data has therefore to be envisaged before fitting the kinetics. Concerning the repeatability of the estimated parameters, a previous study has already tested it using a GD method (25). But it neither expressed the repeatability in the same way than in the present study nor provided SNR_{dB} values; thus only qualitative comparisons are possible. Their results are in accordance with the present ones: these authors stated that the repeatability of *A*_{1} and τ_{1} is well higher than for *A*_{2} and τ_{2} and that the interpretation of these two last indexes is limited. An interesting finding of the present study is that the SA algorithm improves the repeatability of *A*_{2} and τ_{2}, compared with GRG_{2}. Concomitantly, the better precision for SA in estimating *A*_{2} and τ_{2} on simulated data (compared with GRG_{2}) means that SA is less sensible than GRG_{2} to the intersignal differences of noise. This may partly explain the better repeatability coefficients produced by SA if we make the assumption that the within-subject variability of the V̇o_{2} response is due rather to the variability of the noise than to the variability of the underlying physiological response (19). Thus the less sensibility of SA to the intersignal differences of noise on simulated data enables a lesser sensibility to the within-subject variability between repeated experimental measurements.

Another result of the present study is the significant overestimation of td_{2} on simulated data, whatever the estimation algorithm (see Table 3). This may be due to the characteristics of the simulated noise rather than an intrinsic incapacity of both algorithms to estimate td_{2} accurately. Readers must be aware of the fact that a realization of a Gaussian white noise over a long time period may contain short series of values that look like “local perturbations” (i.e., the mean residual substantially deviates from zero on very short local periods). Therefore, the beginning of the second exponential (td_{2}) may be attributed to the time of appearance of a local perturbation when it looks like a “step” in the V̇o_{2} time course. In the present study, the probability that these local perturbations appear beyond the reference td_{2} (td_{2ref}) seems high because the exercise duration beyond td_{2} is long compared with the time spent from the start of exercise to td_{2}; this may explain the overestimation of this time delay by both algorithms. An example of this phenomenon is observed in Fig. 4. Figure 4*A* shows the reference V̇o_{2} kinetics produced by *set 1* of parameters (dark line) and the associated simulated data. Figure 4*B* shows the residual differences between these kinetics and data. Note a local perturbation well apparent between 200 and 220 s (i.e., beyond the time of td_{2ref}). Consequently, a value of 215 s is attributed to td_{2} for the fitted V̇o_{2} kinetics (dotted line in *A*). In real contexts, these local perturbations may happen from perturbations of the breathing pattern such as coughing, sneezing, or swallowing, or from a temporary incapacity for some algorithms to compute accurate V̇o_{2} values (V̇o_{2} is locally overestimated or underestimated) (5). Thus the probability for td_{2} to be biased in real contexts may be high.

The aim of the present paper was to demonstrate the superiority of SA (compared with GRG_{2}) for fitting a given V̇o_{2} signal, but not to discuss the validity of this raw signal for insight into the nature of the underlying physiological processes. However, the above discussion points out that not only the estimation procedure may influence the confidence level of the estimated parameters; the SNR level and the appearance of local perturbations in the raw signal are also of importance. Therefore, special attention must be paid to the quality of this raw signal (i.e., as little noise as possible) before fitting it. A frequently used procedure to reduce the magnitude of the noise and local perturbations is the superposition of several repeated measurements (2, 3, 29). Unfortunately, one of the disadvantages of this filtering procedure is that a great number of repeated measurements is sometimes needed to obtain a desired confidence level in an estimated kinetic parameter. For example, Lamarra et al. (19) proposed a prediction equation of the number of square-wave repetitions required to achieve a desired confidence level when the modeling function is a first-order exponential equation. Applied to our data, it could be predicted that ∼40 repeated measurements would have been needed to reduce the standard deviation of τ_{1} from 9.0 to 1.0 s. Rather than filtering the data a posteriori, attempts may be tried to improve the methodology for computing breath-by-breath V̇o_{2} data. At this time, two main approaches for calculating V̇o_{2} have been developed and used in publications. The first one calculates breath-by-breath O_{2} taken up at the mouth by considering the difference between inspired and expired O_{2} volumes (4, 22). The second one estimates the breath-by-breath alveolar-capillary transfer of O_{2} by accounting for the breath-to-breath changes in alveolar O_{2} stores (1, 5, 9, 16, 31, 34). Some of these algorithms (1, 5, 10, 12, 34) have been shown to largely reduce the magnitude of the breath-by-breath variability of the data (especially those that compute the alveolar gas exchanges). Unfortunately, their use for insight into the nature of the underlying physiological processes is still unsure because none of them induces the same V̇o_{2} time course during transient states (1, 5, 11, 12), and so none of them can be considered as a gold standard. For example, when comparing the different algorithms that compute V̇o_{2} at the alveolar level, the time constant of the phase II (τ) at the onset of a square wave exercise may vary from 34.3 to 46.8 s (11). These considerations lead us to suggest that future research should concomitantly deal with the improvement of both the procedures for assessing V̇o_{2} data and the procedures for estimating the kinetic parameters, when the final goal is to provide insight into the nature of the physiological processes underlying the V̇o_{2} kinetics.

In conclusion, the present investigation demonstrates that the proposed implementation of a global stochastic optimization method (SA) improves the estimations of the V̇o_{2} kinetic parameters of square-wave high-intensity exercises, compared with a GD method (GRG_{2}). However, despite tight limits of confidence for the baseline and the parameters of the first exponential component of the model, a large inaccuracy remains in estimating the parameter values of the second exponential component. This may be due to the shape of the underlying V̇o_{2} response showing, in some cases, linear rather than exponential increase in V̇o_{2} and also due to the large breath-to-breath irregularities of the V̇o_{2} response. Therefore, one of the perspectives of this work concerns an attempt to improve the SNR, by reenvisaging the algorithms used to compute V̇o_{2} data.

## APPENDIX

### Estimation Problem

The breath-by-breath values of V̇o_{2} are sampled with a variable sampling period as the respiratory frequency changes due to the breath-to-breath irregularities and also due to the need for increasing the ventilatory flow rate. V̇o_{2} values can then be viewed as a vector **vo**_{2} = [vo_{2,l}, … , vo_{2,N}], sampled over a time vector **t** = [*t*_{0}, … , *t*_{N}]^{T}. Each *t*_{k} − *t*_{k−1}, *k* = 1, … , *N*, corresponds to the laps of time of a breath, and the first value of **vo**_{2} is obtained at *t*_{1}. **vo**_{2} is supposed to incorporate an underlying physiological response (i.e., the O_{2} exchanges occurring at the alveolar site) plus an additive noise **e** = [*e*_{1}, … , *e*_{N}]^{T}. In the present context, *Eq. 1* is supposed to describe this underlying V̇o_{2} response. Thus it is possible to write **vo**_{2} following the general linear model (GLM) (24): (2) The vectors **g**_{m}, *m* = 0, 1, 2, have a length of *N* and can be written as follows: (3) (4) (5) and (6) (7) The goal of the estimation methods is to estimate the following parametric model: (8) The estimation of θ is done classically by procedures that minimize the RSS between the experimental data points and the model, i.e., the procedures are searching for the set of parameters that minimizes the Euclidean norm of **e** (see *Eq. 2*). In the present case, this is a constrained optimization problem (some constraints exist on the parameter values; see *Estimation algorithm* in methods) on a function that is not continuously differentiable. Moreover, the problem of fitting a sum of exponential to numerical data is known to be ill conditioned. In the following, we present the two algorithms we compared on this estimation problem.

### Simulated Annealing

SA is based on theory of Markov chains (MC) and belongs to the set of stochastic global methods of optimization. SA is an estimation procedure in the maximum likelihood (ML) sense that consists in optimizing the likelihood function *L*. In the case of V̇o_{2} kinetics, *L* can be expressed as follows: (9) where *L*(θ; **vo**_{2}) is the likelihood function and *p*(**vo**_{2} |θ) is theprobability density function of **vo**_{2} conditionally to θ. Optimizing *L* is known to be adapted to the most general contexts in optimization and provides consistent parameter estimation, whatever the probabilistic characteristics of the noise **e**. In the present case, **e** can be considered as a realization of zero mean independently and identically distributed temporal series of Gaussian random variables, i.e., **e** is a realization of a Gaussian white noise with variance σ_{e}^{2}(19, 27, 30). Then, each *e*_{k} has the same probability density function, which is of the form: (10) Because **e** is additive and independently and identically distributed, *L* can be expressed as: (11) Because **e** is Gaussian, using the log-likelihood function can do simplification: (12) Therefore, the ML optimization criterion becomes: (13) where θ̂ is the estimation of θ and σ̂ is the estimation of the variance of the noise.

Let us remark that RSS can be directly obtained from σ̂_{e}^{2}: (14) Used in that context, SA is a stochastic iterative method of RSS minimization. SA is stochastic as it is based on the simulation of a MC. Indeed a new state of a MC is sampled each iteration by using a probabilistic transition kernel that models the transition between two states. When sufficient iterations have been done (and so sufficient states of the chain have been sampled), the probability law of the states converges if the transition kernel has some properties (13, 23, 28).

In our case, a state is defined by a given set of θ_{vo2}: from a random initial set θ_{vo2,0, a new set of parameters θvo2,i} is sampled each iteration *i* of the algorithm. To estimate θ_{vo2} with SA, the convergent probability law of the states above mentioned must be chosen uniform on the sets of parameters that provide global maxima of *Eq. 9* (or global minima of *Eq. 13*) and zero elsewhere. To this end, SA does not optimize directly the likelihood function but an energy function of the states, i.e., the model parameter values, which is proportional to the likelihood function. By definition, this energy function *En*(θ_{vo2}) has a Gibbs distribution (24): (15) where *Z* is a constant of normalization, *B*, the Boltsmann's constant, and *T* the temperature. SA consists in sampling V̇o_{2} according to *Eq. 15* with decreasing values of *T*.

We provide now with the expression of the energy function. From Bayes' theorem, the likelihood function can be expressed as follows: (16)

Because *p*(**vo**_{2}) is an unknown constant and *p*(θ_{vo2}) could also be a constant if an a priori uniform probability law is supposed for the set of parameters, *L*(θ_{vo2};**vo**_{2}) ≈ *p*(θ_{vo2}|**vo**_{2}), which can be expressed after few integration steps (24) (see *Eq. 6* for **G**): (17) The energy function *En* can now be defined as: (18) According to *Eq. 18*, the states that maximize *L*(θ_{vo2};**vo**_{2}) are the states of lowest energies. Then, minimizing *En*(θ_{vo2}) provides ML estimation of θ_{vo2} and so RSS minimization (see *Eqs. 13* and *14*). For the special case *T* = 0, the only possible states according to *Eq. 15* are the ones that produce the global minima of the energy. Therefore, the Gibbs distribution at *T* = 0 is the desired probability law above mentioned.

To iterate SA, two main approaches exist to sample a MC with a transition kernel that has appropriate properties: the Metropolis-Hastings algorithm and the Gibbs sampler. In this paper, we used the Metropolis-Hastings algorithm (for more details, see Ref. 24).

### Gradient-Descent-Based Method (the GRG_{2} Algorithm)

The principle of GD methods like the GRG_{2} algorithm is to find iteratively an optimum for θ_{vo2} by the use of gradient information. The GRG method (20, 21) is an evolution of the Wolfe's reduced gradient method. Although the validity of the reduced gradient method is limited to unconstrained nonlinear problem, with first-order information only, GRG_{2} is adapted to constrained nonlinear optimizations and incorporates second-order information based on reduced Hessian. Although it is guaranteed to find an optimum only on problems with continuously differentiable functions, and then only in the absence of numerical difficulties (such as degeneracy or ill conditioning), GRG_{2} has a reputation for robustness, compared with other nonlinear optimization methods, on difficult problems where these conditions are not fully satisfied. This is the case of estimating θ_{vo2} (see *Estimation Problem*).

GRG_{2} is based on Newton-Raphson iterations to find a root of the gradient of the function *f*(θ_{vo2}), ∇*f*(θ_{vo2}), where *f* is the cost function to optimize. Let suppose θ_{vo2,i + 1} = θ_{vo2,i} + δ_{i}, where δ_{i} is the increment between two consecutive sets of parameters, and consider a linear Taylor expansion to ∇*f*(θ) (24): (19) If θ_{vo2,i} + δ_{i}is the location of the minimum along the direction of θ_{vo2,i}, then (20) where [∇∇^{T} *f*(θ_{vo2},_{i})]^{−1} is the inverse of the Hessian matrix. By definition, a variable metric algorithm estimates iteratively this matrix. In GRG_{2}, the variable metric algorithm in use is the Broyden-Fletcher-Goldfarb-Shanno algorithm. The algorithm converges when ‖ δ_{i} ‖ undergoes a threshold fixed by the user. We fixed with the threshold at 10^{−12}. The GRG_{2} method is implemented in the EXCEL solver (14).

## GRANTS

This work was supported by the Association de Recherche en Physiologie des Activités Physiques et Sportives.

## Acknowledgments

We are grateful to Hervé Billette for technical assistance, as well as the subjects for patient participation.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2006 the American Physiological Society