HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Free All Versions of this Article: 100/3/1049    most recent 00712.2005v1 Submit a response Alert me when this article is cited Alert me when eLetters are posted Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Bernard, O. Articles by Francaux, M. Search for Related Content PubMed PubMed Citation Articles by Bernard, O. Articles by Francaux, M.

INNOVATIVE METHODOLOGY

On the modeling of breath-by-breath oxygen uptake kinetics at the onset of high-intensity exercises: simulated annealing vs. GRG₂ method

¹Laboratoire des Adaptations Physiologiques aux Activités Physiques, Université de Poitiers, Faculté des Sciences du Sport, Poitiers; ²Laboratoire Signal Image Communications, Université de Poitiers, Unité de Formation et de Recherche Sciences Fondamentales et Appliquées, Futuroscope Chasseneuil, France; and ³Institut d'Education Physique et de Réadaptation, Faculté de Médecine, Université Catholique de Louvain, Louvain-la-Neuve, Belgium

Submitted 16 June 2005 ; accepted in final form 19 October 2005

	ABSTRACT

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
Modeling in the time domain, the non-steady-state O₂ 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₂ algorithm (a gradient-descent method known for its robustness). Forty simulated O₂ 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₂ significantly biased two estimated kinetic parameters of the first exponential (the time delay td₁ and the time constant ₁) and impaired the precision (i.e., standard deviation) of the baseline A₀, td₁, and ₁ compared with SA. SA significantly improved the precision of the three parameters of the second exponential (the asymptotic increment A₂, the time delay td₂, and the time constant ₂). Nevertheless, td₂ 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₂ uptake on a cycle ergometer and O₂ 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₂ for A₁ although better with SA for A₂ and ₂. 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₂ uptake kinetics; parametric modeling; parameter estimation

In the context of empirical models applied to O₂ 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 O₂ (O_{2 max})], the O₂ 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% O_{2 max}), the estimation of the kinetic parameters refers to a more complex mathematical problem because the O₂ 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 O₂ 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 O₂ 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₂)] 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

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
Empirical Model to Describe the O₂ Kinetics

To compare SA and GRG₂, 40 simulated O₂ responses and 52 experimental O₂ responses are fitted by both algorithms, with a second-order exponential function of the form

(1)

where

is the unit step function. O₂ is the rate of O₂ uptake at time t. A₁ and A₂ are the asymptotic increments above the baseline A₀, ₁ and ₂ are the time constants, and td₁ and td₂ are the time delays for each exponential term after the start of the exercise. This function is commonly used to describe O₂ on-responses from light- to high-intensity (>60% 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 O₂ 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 O₂ increase (18). In the present study, the O₂ 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 O₂ data measured during the exercise transition are included in the curve fitting.

Estimation Algorithms

The mathematical concepts of both estimation methods (SA and GRG₂) and their implementations are extensively documented in the APPENDIX. In the present study, they are used to fit O₂ on-responses acquired breath by breath. Owing to breathing irregularities, these O₂ 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₂ 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₀, A₁, and A₂ have real positive values, 2) td₁ td₂, and 3) 0 < _{1 2}.

Furthermore, the initialization of GRG₂ 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 O₂ 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.

Simulated data.   Simulated O₂ responses are generated on the basis of real time course and temporal responses. First, experimental O₂ data have been measured during a period of 2 min of a light exercise followed by 8 min of constant-load exercise (80% O_{2 max}) on two different subjects (see Experimental data). Then, the second-order exponential model (Eq. 1) has been used to fit the O₂ 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% O_{2 max}), it is recognized that the O₂ time course described by the second component of Eq. 1 (and called the "slow O_{2 max} component") either delays the attainment of a steady state or leads O₂ to increase inexorably until the end of the exercise (or until O₂ 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 O₂ reaches the asymptotic value of the component before the end of the exercise. Set 2 produces a second exponential component during which O₂ increases almost linearly until the exercise stops. Therefore, both reference O₂ 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 O₂ on-responses with an approximately constant SNR (SNR_dB in decibels) equal to 20 dB. Figure 1 shows examples of simulated O₂ on-responses obtained with set 1 (A) and set 2 (B) of kinetic parameters (the dark lines indicate the underlying second-order exponential curves).

View this table: [in this window] [in a new window]   Table 2. Kinetic parameter sets used to generate simulated O2 on-responses

View larger version (32K): [in this window] [in a new window]   Fig. 1. Typical simulated breath-by-breath O2 uptake (O2) kinetics during a step exercise. Examples of simulated O2 kinetics obtained with set 1 and set 2 of kinetic parameters. The second-order exponential model is superimposed (dark line). A: set 1 delays the attainment of steady state. B: set 2 induces an increase in O2 until the end of the exercise.

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 (²) 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₂ 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 O₂ 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 O₂. 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 O₂ (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 O₂ 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₂ 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

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
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₂ or by the SA algorithm on the whole 40 simulated responses.

View larger version (15K): [in this window] [in a new window]   Fig. 2. Biases and their 95% confidence intervals for the estimation errors of the kinetic parameters of the simulated O2 responses. The bias is the mean difference between the value estimated with either the simulated annealing (SA) algorithm or the generalized reduced gradient 2 (GRG2) algorithm, less the reference value. The 95% confidence intervals are expressed as the inferior and superior limits of confidence. A0, A1, A2, baseline and asymptotic increments for each exponential term of the empirical model (in ml O2/min); 1 and 2, time constants; td1 and td2, time delays for the first and second exponential terms of the empirical model, respectively (in s).

View this table: [in this window] [in a new window]   Table 3. Bias (i.e., mean difference) and precision (i.e., standard deviation) for the estimation of the kinetic parameters of the 40 simulated O2 responses

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₂, td₂, and ₂; 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₂ is significantly biased by both algorithms (P < 0.05).

Experimental Data

Figure 3 presents the linear relationship between RSS values (r² = 0.9995; P < 0.05), with the values produced by the SA algorithm being the independent variables and the values produced by GRG₂ 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₂. Conversely, in the six other situations, the average value produced by the GRG₂ algorithm is only 0.2% lower than the value produced by SA.

View larger version (9K): [in this window] [in a new window]   Fig. 3. Linear relationship between residual sum of squares (RSS) values produced by the SA and GRG2 algorithms, on experimental O2 responses.

View this table: [in this window] [in a new window]   Table 4. Kinetic parameters estimated on experimental O2 responses

	DISCUSSION

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

Comparison of both algorithms (SA vs. GRG₂) 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 O₂ 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 O₂ 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₂) on simulated O₂ 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₂ algorithm. The second explaining factor is that the SA algorithm is a stochastic unsupervised method, whereas GRG₂ needs to be initialized. Therefore, the estimated parameters found by GRG₂ 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₂.

Our results on simulated data also demonstrate that the confidence levels observed for the kinetic parameters of the first exponential component (td₁ and ₁) 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₁ and ₁, 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₁ and ₁, 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 O₂ 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₁ and ₁ is comparatively less with GRG₂ (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 O₂ 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 ₁ and A₁, 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₂, td₂, and ₂) and also the repeatability of A₂ and ₂ (P < 0.05 compared with the GRG₂ 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 O₂ kinetics. On simulated data, set 2 of reference kinetic parameters induces a linear increase of O₂ (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₂-to-₂ ratio is a good approximation of the slope of a linear function describing this increase in O₂; 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 O₂ 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 O₂ 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 O₂ 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₁ and ₁ is well higher than for A₂ and ₂ 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₂ and ₂, compared with GRG₂. Concomitantly, the better precision for SA in estimating A₂ and ₂ on simulated data (compared with GRG₂) means that SA is less sensible than GRG₂ 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 O₂ 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₂ 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₂ 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₂) may be attributed to the time of appearance of a local perturbation when it looks like a "step" in the O₂ time course. In the present study, the probability that these local perturbations appear beyond the reference td₂ (td_2ref) seems high because the exercise duration beyond td₂ is long compared with the time spent from the start of exercise to td₂; this may explain the overestimation of this time delay by both algorithms. An example of this phenomenon is observed in Fig. 4. Figure 4A shows the reference O₂ kinetics produced by set 1 of parameters (dark line) and the associated simulated data. Figure 4B 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₂ for the fitted O₂ 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 O₂ values (O₂ is locally overestimated or underestimated) (5). Thus the probability for td₂ to be biased in real contexts may be high.

View larger version (39K): [in this window] [in a new window]   Fig. 4. A: comparison of the reference O2 kinetics obtained with set 1 of parameters (dark line) and the fitted O2 kinetics obtained with SA (dotted line) to 1 simulated O2 response (incorporating added noise). B: residual differences between the simulated data and the reference kinetics. Note that the reference time delay 2 (td2ref) equals 140 s and that the mean residual between 200 and 220 s deviates clearly from zero. SA estimates a td2 value (td2est) equal to 215 s.

In conclusion, the present investigation demonstrates that the proposed implementation of a global stochastic optimization method (SA) improves the estimations of the O₂ kinetic parameters of square-wave high-intensity exercises, compared with a GD method (GRG₂). 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 O₂ response showing, in some cases, linear rather than exponential increase in O₂ and also due to the large breath-to-breath irregularities of the O₂ response. Therefore, one of the perspectives of this work concerns an attempt to improve the SNR, by reenvisaging the algorithms used to compute O₂ data.

	APPENDIX

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
Estimation Problem

The breath-by-breath values of O₂ 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. O₂ values can then be viewed as a vector vo₂ = [vo_2,l, ... , vo_2,N], sampled over a time vector t = [t₀, ... , 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₂ is obtained at t₁. vo₂ is supposed to incorporate an underlying physiological response (i.e., the O₂ exchanges occurring at the alveolar site) plus an additive noise e = [e₁, ... , e_N]^T. In the present context, Eq. 1 is supposed to describe this underlying O₂ response. Thus it is possible to write vo₂ 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 O₂ kinetics, L can be expressed as follows:

(9)

where L(; vo₂) is the likelihood function and p(vo₂ |) is theprobability density function of vo₂ 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²(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²:

(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 O₂ 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₂) 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₂) p(_vo2|vo₂), 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₂) 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₂ Algorithm)

The principle of GD methods like the GRG₂ 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₂ 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₂ 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₂ 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 + _iis 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₂, 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₂ method is implemented in the EXCEL solver (14).

	GRANTS

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

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

	ACKNOWLEDGMENTS

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

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

	FOOTNOTES

 

Address for reprint requests and other correspondence: O. Bernard, Laboratoire des Adaptations Physiologiques aux Activités Physiques (LAPHAP EA 3813), Université de Poitiers, Faculté des Sciences du Sport, 4 allée Jean Monnet, 86000 Poitiers, France (e-mail: olivier.bernard{at}univ-poitiers.fr)

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.

	REFERENCES

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 

1. Auchincloss JH Jr, Gilbert R, and Baule GH. Effect of ventilation on oxygen transfer during early exercise. J Appl Physiol 21: 810–818, 1966.[Free Full Text]
2. Barstow TJ and Molé PA. Linear and nonlinear characteristics of oxygen uptake kinetics during heavy exercise. J Appl Physiol 71: 2099–2106, 1991.[Abstract/Free Full Text]
3. Barstow TJ, Casaburi R, and Wasserman K. O₂ uptake kinetics and the O₂ deficit as related to exercise intensity and blood lactate. J Appl Physiol 75: 755–762, 1993.[Abstract/Free Full Text]
4. Beaver WL, Wasserman K, and Whipp BJ. On-line computer analysis and breath-by-breath graphical display of exercise function tests. J Appl Physiol 34: 128–132, 1973.[Free Full Text]
5. Beaver WL, Lamarra N, and Wasserman K. Breath-by-breath measurement of true alveolar gas exchange. J Appl Physiol 51: 1662–1675, 1981.[Abstract/Free Full Text]
6. Bernard O, Maddio F, Ouattara S, Jimenez C, Charpenet A, Melin B, and Bittel J. Influence of the oxygen uptake slow component on the aerobic energy cost of high-intensity submaximal treadmill running in humans. Eur J Appl Physiol 78: 578–585, 1998.[CrossRef][Web of Science]
7. Bernard O, Ouattara S, Maddio F, Jimenez C, Charpenet A, Melin B, and Bittel J. Determination of the velocity associated with O₂. Med Sci Sports Exerc 33: 464–470, 2000.
8. British Standards Institution. Precision of Test Methods I: Guide for the Determination and Reproducibility for a Standard Test Method (NS 5497, Part I). London: BSI, 1979.
9. Busso T and Robbins PA. Evaluation of estimates of alveolar gas exchange by using a tidally ventilated nonhomogeneous lung model. J Appl Physiol 82: 1963–1971, 1997.[Abstract/Free Full Text]
10. Capelli C, Cautero M, and di Prampero PE. New perspectives in breath-by-breath determination of alveolar gas exchange in humans. Pflügers Arch 441: 566–577, 2001.[CrossRef][Web of Science][Medline]
11. Cautero M, Beltrami AP, di Prampero PE, and Capelli C. Breath-by-breath alveolar oxygen transfer at the onset of step exercise in humans: methodological implications. Eur J Appl Physiol 88: 203–213, 2002.[CrossRef][Web of Science][Medline]
12. Cautero M, di Prampero PE, and Capelli C. New acquisitions in the assessment of breath-by-breath alveolar gas transfer in humans. Eur J Appl Physiol 90: 231–241, 2003.[Medline]
13. Fitzgerald WJ. Markov Chain Monte Carlo methods with applications to signal processing. Signal Processing 81: 3–18, 2001.[CrossRef]
14. Fylstra D, Lasdon L, Watson J, and Waren A. Design and use of the Microsoft Excel Solver. INTERFACES 28: 29–55, 1998.[Abstract/Free Full Text]
15. Gilks WR, Richardson S, and Spieglehalter DJ. Markov Chain Monte Carlo in Practice. New York: Chapman and Hall, 1996.
16. Grønlund J. A new method for breath-by-breath determination of oxygen flux across the alveolar membrane. Eur J Appl Physiol 52: 167–172, 1984.[CrossRef]
17. Hughson RL, Sherrill DL, and Swanson GD. Kinetics of O₂ with impulse and step exercise in humans. J Appl Physiol 64: 451–459, 1988.[Abstract/Free Full Text]
18. Inman MD, Hughson RL, Weisiger KH, and Swanson GD. Estimate of mean tissue O₂ consumption at onset of exercise in males. J Appl Physiol 63: 1578–1585, 1987.[Abstract/Free Full Text]
19. Lamarra N, Whipp BJ, Ward SA, and Wasserman K. Effect of interbreath fluctuations on characterizing exercise gas exchange kinetics. J Appl Physiol 62: 2003–2012, 1987.[Abstract/Free Full Text]
20. Lasdon LS, Fox RL, and Ratner MW. Nonlinear optimization using the generalized reduced gradient method. RAIRO 3: 73–104, 1974.
21. Lasdon LS, Waren AD, Jain A, and Ratner M. Design and testing of a generalized reduced gradient code for nonlinear programming. ACM Transactions on Mathematical Software 4: 34–49, 1978.
22. Linnarsson D and Lindborg B. Breath-by-breath measurement of respiratory gas exchange using on-line analog computation. Scand J Clin Lab Invest 34: 219–224, 1974.[Medline]
23. Neal RM. Probabilistic Inference using Markov Chain Monte Carlo Methods. Technical Report CRG-TR-93-1. University of Toronto, Ontario, Canada, 1993.
24. O'Ruanaidh JK and Fitzgerald WJ. Numerical Bayesian Methods Applied to Signal Processing. New York: Springer-Verlag, 1996.
25. Özyener F, Rossiter HB, Ward SA, and Whipp BJ. Influence of exercise intensity on the on- and off-transient kinetics of pulmonary oxygen uptake in humans. J Physiol 533: 891–902, 2001.[Abstract/Free Full Text]
26. Puente-Maestu L, Sánz ML, Sánz P, Nuñez A, González F, and Whipp BJ. Reproducibility of the parameters of the on-transient cardiopulmonary responses during moderate exercise in patients with chronic obstructive pulmonary disease. Eur J Appl Physiol 85: 434–441, 2001.[CrossRef][Web of Science][Medline]
27. Puente-Maestu L, Abad MJB, Godoy R, Pérez-Parra JM, Cubillo JM, and Whipp BJ. Breath-by-breath fluctuations of pulmonary gas exchange and ventilation in COPD patients. Eur J Appl Physiol 87: 535–541, 2002.[Medline]
28. Robert CP and Casella G. Monte Carlo Statistical Method. New York: Springer Verlag, 1999.
29. Rossiter HB, Ward SA, Doyle VL, Howe FA, Griffiths JR, and Whipp BJ. Inferences from pulmonary O₂ uptake with respect to intramuscular [phosphocreatine] kinetics during moderate exercise in humans. J Physiol 518: 921–932, 1999.[Abstract/Free Full Text]
30. Rossiter HB, Howe FA, Ward SA, Kowalchuk JM, Griffiths JR, and Whipp BJ. Intersample fluctuations in phosphocreatine concentration determined by ³¹P-magnetic resonance spectroscopy and parameter estimation of metabolic responses to exercise in humans. J Physiol 528: 359–369, 2000.[Abstract/Free Full Text]
31. Swanson GD and Sherrill DL. A model of estimates of breath-by-breath alveolar gas exchange. J Appl Physiol 55: 1936–1941, 1983.[Abstract/Free Full Text]
32. Tschakovsky ME and Hughson RL. Interaction of factors determining oxygen uptake at the onset of exercise. J Appl Physiol 86: 1101–1113, 1999.[Abstract/Free Full Text]
33. Weissman ML, Jones PW, Oren A, Lamarra N, Whipp BJ, and Wasserman K. Cardiac output increase and gas exchange at the start of exercise. J Appl Physiol 52: 236–244, 1982.[Medline]
34. Wessel HU, Stout RL, Bastabier CK, and Paul MH. Breath-by-breath variation of FRC: effect on O₂ and O₂ measured at the mouth. J Appl Physiol 46: 1122–1126, 1979.[Abstract/Free Full Text]
35. Whipp BJ, Ward SA, Lamarra N, Davis JA, and Wasserman K. Parameters of ventilatory and gas exchange dynamics during exercise. J Appl Physiol 52: 1506–1513, 1982.[Abstract/Free Full Text]

This Article Abstract Full Text (PDF) Free All Versions of this Article: 100/3/1049    most recent 00712.2005v1 Submit a response Alert me when this article is cited Alert me when eLetters are posted Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Bernard, O. Articles by Francaux, M. Search for Related Content PubMed PubMed Citation Articles by Bernard, O. Articles by Francaux, M.