 Research Article
 Open Access
 Published:
Modeling the oxygen uptake kinetics during exercise testing of patients with chronic obstructive pulmonary diseases using nonlinear mixed models
BMC Medical Research Methodology volume 16, Article number: 66 (2016)
Abstract
Background
The sixminute walk test (6MWT) is commonly used to quantify exercise capacity in patients with several cardiopulmonary diseases. Oxygen uptake (\(\dot {\mathrm {V}}\)O_{2}) kinetics during 6MWT typically follow 3 distinct phases (rest, exercise, recovery) that can be modeled by nonlinear regression. Simultaneous modeling of multiple kinetics requires nonlinear mixed models methodology. To the best of our knowledge, no such curvefitting approach has been used to analyze multiple \(\dot {\mathrm {V}}\)O_{2} kinetics in both research and clinical practice so far.
Methods
In the present study, we describe functionality of the R package medrc that extends the framework of the commonly used packages drc and nlme and allows fitting nonlinear mixed effects models for automated nonlinear regression modeling. The methodology was applied to a data set including 6MWT \(\dot {\mathrm {V}}\)O_{2} kinetics from 61 patients with chronic obstructive pulmonary disease (disease severity stage II to IV). The mixed effects approach was compared to a traditional curvebycurve approach.
Results
A sixparameter nonlinear regression model was jointly fitted to the set of \(\dot {\mathrm {V}}\)O_{2} kinetics. Significant differences between disease stages were found regarding steady state \(\dot {\mathrm {V}}\)O_{2} during exercise, \(\dot {\mathrm {V}}\)O_{2} level after recovery and \(\dot {\mathrm {V}}\)O_{2} inflection point in the recovery phase. Estimates obtained by the mixed effects approach showed standard errors that were consistently lower as compared to the curvebycurve approach.
Conclusions
Hereby we demonstrate the novelty and usefulness of this methodology in the context of physiological exercise testing.
Background
The 6Minute Walk Test (6MWT) is routinely used to quantify submaximal exercise capacity in patients with chronic cardiopulmonary diseases [1]. It is considered as the test of choice to quantify functional capacity and patient’s daily life activities [2]. Portable wireless cardiopulmonary exercise testing devices enable to measure breathbybreath oxygen exchange kinetics during exercise performance. From the original acquired breathbybreath oxygen exchange data, curves can be generated and parameters can be estimated by curve model fitting [3, 4]. These parameters provide information about the interaction of the cardiovascular, cardiac autonomic, pulmonary, and metabolic system.
In the field of pulmonary medicine, there is a growing interest in estimating exercise parameters capable of objectively evaluating the functional capacity of patients with chronic obstructive pulmonary disease (COPD). Time course measurements of oxygen consumption (\(\dot {\mathrm {V}}\)O_{2}) provide key physiological determinants of exercise capacity which are hardly influenced by conditions other than the underlying disease. Modeling \(\dot {\mathrm {V}}\)O_{2} kinetics during exercise testing in patients with COPD is therefore a revelant topic.
Traditionally, individual \(\dot {\mathrm {V}}\)O_{2} uptake curves have been analyzed visually and parameters representing for example the time needed for a 50 % increase in oxygen uptake (T _{1/2}) were quantified in an inconsistent manner [5]. Nowadays, \(\dot {\mathrm {V}}\)O_{2} kinetics can be accurately modeled by fitting nonlinear regression models [6–8]. A modeling approach insures that parameters are estimated in an objective and systematic manner, fully exploiting the information available in the kinetics data.
Several difficulties arise when dealing with nonlinear regression models. Unlike linear regression where mathematical solutions exist for the estimation of the best fitted parameters, nonlinear regression often requires additional information from users who should specify among others, the equation of the model, a set of starting parameters for the minimization procedure, and pay a particular attention to the check of the models assumptions [8].
Routinely acquired oxygen uptake kinetics generate batches of curves whose parameters are generally estimated using curvebycurve model fits. The oxygen kinetics are sometimes originated from individuals which can be grouped into categories (e.g., distinct diagnostic group, disease severity, etc.) and one might be interested in testing differences between parameter estimates with regard to these predefined categories. Hypothesis testing is classically done in two integral steps. First, curve fitting is used and as a result a set of parameter estimates can be produced. In a second step, inference based on the estimated parameters is carried out. However, in this type of data, withingroup correlation typically arise and specific regression techniques must be used [9]. Mixed models are wellsuited for the analysis of grouped data. They allow to explicitely incorporate intragroup correlations by means of random effects and get joint estimates of the regression model parameters [10]. Fitting nonlinear mixed models comes with the same type of challenges as fitting nonlinear regression models.
The whole analytical workflow requires various statistical expertise. The aim of the current work is to demonstrate how recent methodological developments enable to carry out these series of successive procedures — namely the simultaneaous fit of nonlinear models within the mixed effects framework, followed by inferences based on parameter estimates — in an automated and allintegrated manner, by using recently developed packages from the R statistical software [11].
The next sections are organized as follows. The novel mixed effects methodology is described in the “Material and methods”. A detailed example of physiological exercise testing from patients with COPD is provided in the “Results”. The mixed models approach is compared with a traditional curvebycurve approach. The statistical and physiological relevance of the findings is exposed in the “Discussion”.
Material and methods
COPD dataset
Patients with COPD referred to the Department of Pulmonary Medicine of the University Hospital of Basel (Switzerland) for a 6MWT gave their informed consent to participate to the study. The data analyzed in the present study were fully anonymized and no individual clinical data are presented. The study was approved from the local institutional review board (Ethikkommission beider Basel). The study was conducted in accordance with the principles enunciated in the Declaration of Helsinki and the guidelines of Good Clinical Practice. Further details on the study can be found on previous publications [7, 12, 13]. The supporting data set is provided in Additional file 1.
Oxygen monitoring during 6MWT
We used the Oxycon Mobile ^{Ⓡ} (Viasys Healthcare, USA) portable, wireless cardiopulmonary exercise testing device to measure breathbybreath \(\dot {\mathrm {V}}\)O_{2} kinetics. Pulse rate was determined by using an ECGtriggered belt (Polar ^{Ⓡ} Electro OY T61). Blood oxgen saturation level (SpO_{2}) was measured by using a finger clip. \(\dot {\mathrm {V}}\)O_{2} and carbon dioxide output (\(\dot {\mathrm {V}}\)CO_{2}), tidal volumes and breathing frequency were assessed by using a facemask (dead space <70 mL) with a flow sensor and a gas analyzer. The patient carried data storage and transfer units by using a dedicated harness. Wireless transfer of breathbybreath data to a laptop computer allowed realtime monitoring. The additional weight (950 g) of the equipment has no effect on walking distance [12]. The exact 6MWT procedure with mobile telemetry has been previously described [12].
Original breathbybreath data were imported from the mobile telemetry device. Raw data were preprocessed by averaging the breathbybreath measurements over consecutive periods of 20 s.
Oxygen uptake kinetics during exercise testing
The American Thoracic Society (ATS) published detailed instructions regarding how to conduct the 6MWT as much standardized as possible. Accordingly, the 6MWT consists of 3 phases; five minutes of rest, six minutes of walking and five minutes of recovery. The oxygen uptake kinetics (\(\dot {\mathrm {V}}\)O_{2}) of each phase can be modeled by nonlinear regression as demonstrated in Fig. 1. A constant O_{2} consumption is expected in the resting phase, followed by a 6minute monotonic O_{2} increase, and a progressive recovery phase where O_{2} decreases back to its initial level.
Established parameters for the assessment of cardiopulmonary exercise capacity in patients with COPD during the 6MWT are the oxygen uptake at steady state (\(\dot {\mathrm {V}}\mathrm {O}_{2}ss\)) and the sixminute walking distance (6MWD) [14].
The incline of \(\dot {\mathrm {V}}\)O_{2} during the initial phase of lowintensity exercise (\(\dot {\mathrm {V}}\)O_{2} onkinetics) provides important information about oxygen delivery and muscle metabolism, and are found to be noticeably delayed in patients with several chronic cardiac and pulmonary diseases [15–19]. \(\dot {\mathrm {V}}\)O_{2} onkinetics can be quantified by the time (mean response time: MRT) required for \(\dot {\mathrm {V}}\)O_{2} to achieve 63 % of the \(\dot {\mathrm {V}}\mathrm {O}_{2}ss\) in response to exercise [20, 21].
Patients with COPD are limited during exercise by dyspnea and fatigue and also have difficulties to recover normal breathing after exercise. As the limited cardiopulmonary reserve in patients with COPD appears to affect exercise responses it may be postulated that it also affects the recovery phase. Recovery kinetics of oxygen uptake (\(\dot {\mathrm {V}}\)O_{2} offtransient kinetics) reflect the ability to recover from exercise that is indicative of daily life.
\(\dot {\mathrm {V}}\)O_{2} offtransient kinetics can be quantified by both the steepness of the curve as well as the time necessary for \(\dot {\mathrm {V}}\)O_{2} to recover by 50 % from its peak effort value (\(\mathrm {T}_{1/2}\dot {\mathrm {V}}\mathrm {O}_{2}\)).
An interpretable nonlinear regression model
A sixparameter nonlinear regression model describing the 3 main phases of the oxygen kinetics (before, during and after 6MWT) is defined as follows:
with \(\dot {\mathrm {V}}\mathrm {O}_{2}rest\), \(\dot {\mathrm {V}}\mathrm {O}_{2}ss\) and \(\dot {\mathrm {V}}\mathrm {O}_{2}recovery\) the oxygen level at rest, steady state during exercise and recovery, respectively; τ _{1} the growth rate of the monoexponential \(\dot {\mathrm {V}}\)O_{2} function during 6MWT; τ _{2} the steepness of the exponential decay during the recovery phase and \(\mathrm {T}_{1/2}\dot {\mathrm {V}}\mathrm {O}_{2}\) the time for half decrease of the \(\dot {\mathrm {V}}\)O_{2} level in the recovery phase. All the 6 previously mentioned parameters (\(\dot {\mathrm {V}}\mathrm {O}_{2}rest\), \(\dot {\mathrm {V}}\mathrm {O}_{2}ss\), τ _{1}, \(\mathrm {T}_{1/2}\dot {\mathrm {V}}\mathrm {O}_{2}\), τ _{2}, \(\dot {\mathrm {V}}\mathrm {O}_{2}recovery\)) have to be estimated by nonlinear regression procedure. λ is the length of the resting period. It is controlled by the experimenter who decides when the resting period ends (usually after 5 min of rest) by initiating the 6MWT. The experimenter reports manually the exact duration of the resting period (λ) which is therefore not estimated during the fitting procedure. λ _{ max } corresponds to the length of the longest resting period in the set of kinetics. This value is known a priori hence not estimated during the fitting procedure. It is simply used to “align” multiple kinetics by removing differences among the duration of individual resting phases.
This nonlinear model takes into account the basic experimental and physiological specificities of a 6MWT, including a resting phase whose duration is controlled by the experimenter, followed by an immediate monoexponential raise of oxygen during 6min exercise, followed by a progressive decline of \(\dot {\mathrm {V}}\)O_{2} during the recovery phase. It is worth noting that this sixparameter model including a single monoexponential function that describes the oxygen increase during exercise is a rather common choice [6] which suits well the purpose of our simple clinical application. More complex models describing the increase of oxygen consumption by means of two or threeterm exponential function have been proposed. However, they are only useful for the modeling of physiological processes occurring in specific situations such as heavyintensity exercise [22, 23].
Mixed effects modeling
Following the notation of Davidian and Giltinan [10], a nonlinear regression model for hierarchical data (e.g., patients and 6MWT measurements nested within patients) can be defined in two stages: First modeling the variability within the i ^{th} patient, and hereby incorporating betweenpatient variation.
Stage 1: For i=1,⋯,m individuals, the following models can be assumed:
where y _{ ij } are the response vectors of length j=1,⋯,n _{ i } with the corresponding vectors individual times t _{ ij }. The nonlinear function such as the above sixparameter model (Eq. 1) evaluated at time t _{ ij } is denoted by \(\dot {\mathrm {V}}\mathrm {O}_{2}(t_{ij}, \beta _{i})\) with a pdimensional individualspecific parameter β _{ i }. The residual vectors \(\epsilon _{ij} \sim \mathcal {N}(0, \sigma ^{2} \Lambda _{i})\) are assumed to be normally distributed with a correlation structure defined by the elements of matrices Λ _{ i }; for the COPD data we will assume that Λ _{ i } is the identity matrix.
The curve is described by the functions \(\dot {\mathrm {V}}\mathrm {O}_{2}(t_{ij}, \beta _{i})\) with an individualspecific (p×1) vector of parameters β _{ i }.
Stage 2: Betweenpatient effects are described by modeling the β _{ i }. These effects are separated into fixed and random effects as in an ordinary linear mixed model (except there is no residual error as it was already introduced in Stage 1):
where β is the vector of fixedeffects parameters, which may differ between patients according to recorded patient characteristics encoded in the design matrix A _{ i } (e.g., age or sex). Differences between patients, which are not captured by the recorded patient characteristics, are described by the patientspecific random effects vector b _{ i }; these random effects may possibly also be modified through explanatory variables encoded in the corresponding design matrix B _{ i } (e.g., time).
As the random effects are intended for capturing inexplicable effects, which may balance out on average, they may be assumed to follow a meanzero, possibly multivariate normal distribution: \(b_{i} \sim \mathcal {N}(0, G)\) where G denotes the betweenpatient variancecovariance matrix, which we assumed to be entirely unstructured such all entries (covariance and variance parameters) were estimated from the data.
In our example, the overall mixed model was parametrized as follows: each individual \(\dot {\mathrm {V}}\)O_{2} kinetics defines one cluster for which a different fixed effect curve is assumed depending on the disease severity stage (all 6 parameters: \(\dot {\mathrm {V}}\mathrm {O}_{2}rest\), \(\dot {\mathrm {V}}\mathrm {O}_{2}ss\), τ _{1}, \(\mathrm {T}_{1/2}\dot {\mathrm {V}}\mathrm {O}_{2}\), τ _{2}, \(\dot {\mathrm {V}}\mathrm {O}_{2}recovery\)); random effects were specified for the five parameters \(\dot {\mathrm {V}}\mathrm {O}_{2}rest\), \(\dot {\mathrm {V}}\mathrm {O}_{2}ss\), τ _{1}, \(\mathrm {T}_{1/2}\dot {\mathrm {V}}\mathrm {O}_{2}\), τ _{2}. In this particular case, A _{ i } is defined as the dummy coded design matrix specifying disease stage specific fixedeffect parameters, and B _{ i } is defined as the random effect design matrix using dummy coded patient identifiers for 5 of the nonlinear model parameters.
We used Akaike’s information criterion (AIC) for identifying the appropriate random effects structure [24, 25].
Sensitivity analysis
The curvebycurve approach provides parameter estimates from fitting a nonlinear regression model (Eq. 1) to each individual patients kinetics data. Subsequently, for each of the six model parameters the corresponding parameter estimates of all patients may be analyzed by means of analysis of variance models to obtain estimates for each disease stage. This twostep curvebycurve approach only relies on being able to fit the patient curves one by one. This approach resulted in estimated coefficients that could be interpreted as population level effects whereas nonlinear mixed effects models allow patientspecific interpretations.
Implementation
The whole analytical workflow relies on the R packages drc [26] and nlme [27] for the nonlinear regression with mixed effects. The newly developed package medrc [28] elegantly combines functionalities of the packages drc and nlme. Finally the package multcomp was used for the statistical inference [29]. The main features of these packages are summarized below.

drc allows simultaneous fit of several nonlinear regression models [26]. It provides automated fit of a list of nonlinear models by directly specifying initial parameter values (selfstarters) for the estimation of the nonlinear model parameters. The main function in drc is drm(), which requires as arguments the name of the data set, the selfstarter model, the dependent and indepedent variables.

nlme provides tools for the fit of Gaussian nonlinear mixed models [27]. It supports specifying correlation structure for residuals and is particulary adapted for repeated measures designs. This package will be used indirectly through the package medrc described next.

medrc combines the automated nonlinear regression modeling framework of the package drc with the nonlinear mixed estimation framework of the package nlme. The function medrm() in medrc allows to fit nonlinear mixed models, by providing the following arguments: form the formula with the response (\(\dot {\mathrm {V}}\)O_{2}) on the left side as a function of the independent variable (time) on the right side; curveid the name of the categorical variable that divides the dataset into several clusters; data the name of the data set; fct the definition of the model; random the definition of the random effects. medrc is available at the following github repository: https://github.com/danielgerhard/medrc

multcomp [29] is used for parameter inference, by providing functionalities for multiple tests for the fixed effects of the mixed models.
A sample R code is provided in Additional file 2.
Results
Application to sixminute walk test oxygen kinetics in patients with COPD
\(\dot {\mathrm {V}}\)O_{2} kinetics were measured in 61 patients with COPD who were traditionally classified into 3 disease severity stages (GOLD II, III and IV). A summary of the patient characteristics is presented in Table 1.
The sixparameter nonlinear regression model was jointly fitted to the set of 61 oxygen kinetics using the function medrm() in the package medrc. The initial model selection of the random effects structure resulted in the minimum AIC for a model with random effects assigned to five out of the six model parameters; results from this model are reported below. The detailed results of the model selection process are provided in the Additional file 3: Table S1. It is worth noting that the model with random effects assigned to all six model parameters did not converge. Graphical check of the residuals is provided in the Additional file 4: Figure S1.
Figure 2 displays the 61 raw oxygen kinetics data together with the fitted curves summarized within each COPD disease severity stage. The five estimated variance parameters for the random effects were: \(SD_{\tau _{1}} = 44.09\); \(SD_{\dot {\mathrm {V}}\mathrm {O}_{2}rest} = 58.26\); \(SD_{\dot {\mathrm {V}}\mathrm {O}_{2}ss} = 250.51\); \(SD_{\tau _{2}} = 1.05\); \(SD_{\mathrm {T}_{1/2}\dot {\mathrm {V}}\mathrm {O}_{2}} = 44.31\).
Significant differences between disease stages were found regarding steady state oxygen uptake during exercise testing (\(\dot {\mathrm {V}}\mathrm {O}_{2}ss\) II vs. IV, adj. p=0.038), oxygen level after recovery (\(\dot {\mathrm {V}}\mathrm {O}_{2}recovery\) II vs. IV, adj. p=0.0013) and inflection point in the recovery phase (\(\mathrm {T}_{1/2}\dot {\mathrm {V}}\mathrm {O}_{2}\) II vs. IV, adj. p=0.088). There were significant differences when comparing the peak oxygen level reached during exercise (\(\dot {\mathrm {V}}\mathrm {O}_{2}ss\)) between moderate and very severe COPD disease stages. On average, patients with COPD stage IV reached a \(\dot {\mathrm {V}}\mathrm {O}_{2}ss\) 292.5 mL (s e=96.8) lower than patients with stage II during the 6MWT. These cardiopulmonary limitations are also reflected by the differences found in the recovery phase. Patients with stage IV needed on average 49.7 (18.1) additional seconds to reduce half of their oxygen consumption in the recovery phase as compared to patients with stage II.
Figure 3 shows in parallel the fitted curves obtained from the joint nonlinear mixed model and the curves obtained from the curvebycurve approach. The overall mean curve obtained by pooling all data from all patients is represented by a thick gray curve. The estimates obtained from both approaches are given in Table 2, which shows that the standard errors are consistently (with one exception) smaller when using the mixed model (up to three times smaller).
Discussion
A set of oxygen kinetics data that originated from a 6MWT in patients with COPD could be successfully analyzed both using a joint nonlinear mixed model and using a twostep curvebycurve approach. Joint modeling is challenging when it comes to selecting the appropriate random effects structure, which is an essential part of the model as it allows to separate out patient variation, as models may not converge. We found that the more random effects were included the better the model as judged by AIC. However, for mixed models there exist alternative information criteria such as the conditional AIC [30] which perhaps may strike a better balance between model complexity and computational feasability. In contrast, the curvebycurve approach is simpler, more robust, and operational. By averaging the patient specific curvebycurve estimates, the possibly varying standard errors of the individual estimates are not taken into account. Specifically, comparison of results of the joint mixed model approach with the classical curvebycurve approach revealed differences between the two approaches. The fitted regression curves obtained by the curvebycurve approach were shrunk towards the overall population mean curve, demonstrating the wellknown effect of attenuation on the marginal estimates as compared to conditional estimates [10].
However, the key difference between the two approaches was the reduction in estimated standard errors when using a mixed model approach, which successfully managed to remove a substantial part of the betweenpatient variation in the fixedeffects estimates. Our study indicates that the joint nonlinear mixed model is to be preferred over the twostep curvebycurve approach. In general joint modeling should be preferred if the betweenpatient variation is large as we would expect it to be in many biological and medical studies [10]. Another advantage of joint modeling is that curves with few time points may still be fitted as estimation in the mixed model borrows strength from curves with many time points to curves with few time points. The curvebycurve approach may fail to utilize such scarce data.
There is some evidence that misspecification of the distributions of the random effects (as normal distributions) may not severely impact inference on fixed effects parameters [31]. However, it is crucial to ensure approximately normally distributed residual errors. One approach is the transformbothsides approach, which was originally proposed for ordinary nonlinear regression models, but it may be extended to nonlinear mixedeffect regression models: it was used with the BoxCox family of power transformations, which include the logarithm [32]. A related approach is to model the residual variance in terms of a covariate [33]. An entirely different approach would be to consider models with other distributional assumptions than normality. To our knowledge such models are only readily available within a Bayesian framework (e.g., the MONOLIX software [34]) and still they involve considerable manual programming in specialized software.
Conclusion
Recent statistical developments provide experimenters with a variety of tools for fitting nonlinear mixed models. Within the statistical environment R the package medrc provides a comprehensive and flexible framework for the parametrisation and inference of hierarchical nonlinear mixedeffects regression models with various biological and medical applications, as exemplified by an application from pulmonary medicine.
References
 1
ATS. ATS/ACCP Statement on cardiopulmonary exercise testing. Am J Respir Crit Care Med. 2003; 167(2):211–77.
 2
ATS. ATS statement: guidelines for the sixminute walk test. Am J Respir Crit Care Med. 2002; 166(1):111–7.
 3
Weber KT, Kinasewitz GT, Janicki JS, Fishman AP. Oxygen utilization and ventilation during exercise in patients with chronic cardiac failure. Circulation. 1982; 65(6):1213–1223.
 4
Solal AC, Chabernaud JM, Gourgon R. Comparison of oxygen uptake during bicycle exercise in patients with chronic heart failure and in normal subjects. J Am Coll Cardiol. 1990; 16(1):80–5.
 5
Whipp BJ, Ward SA, Lamarra N, Davis JA, Wasserman K. Parameters of ventilatory and gas exchange dynamics during exercise. J Appl Physiol Respir Environ Exerc Physiol. 1982; 52(6):1506–1513.
 6
Koga S, Shiojiri T, Kondo N. Measuring vo2 kinetics: the practicalities In: Jones A. M, Poole D. C, editors. Oxygen Uptake Kinetics in Sport, Exercise and Medicine. London and New York: Routledge: 2005. p. 39–61.
 7
Kern L, Condrau S, Baty F, Wiegand J, van Gestel AJ, Azzola A, Tamm M, Brutsche M. Oxygen kinetics during 6minute walk tests in patients with cardiovascular and pulmonary disease. BMC Pulm Med. 2014; 14:167.
 8
Baty F, Ritz C, Charles S, Brutsche M, Flandrois JP, DelignetteMuller ML. A toolbox for nonlinear regression in R: The package nlstools. J Stat Soft. 2015; 66(5):1–21.
 9
Pinheiro JC. Conditional versus Marginal Covariance Representation for Linear and Nonlinear Models. Austrian J Stat. 2006; 35(1):31–44.
 10
Davidian M, Giltinan DM. Nonlinear Models for Repeated Measurement Data. New York: Chapman & Hall/CRC; 1995.
 11
R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2015. R Foundation for Statistical Computing. http://www.Rproject.org/. Accessed 24 May 2016.
 12
Tueller C, Kern L, Azzola A, Baty F, Condrau S, Wiegand J, Tamm M, Brutsche M. Sixminute walk test enhanced by mobile telemetric cardiopulmonary monitoring. Respiration. 2010; 80(5):410–8.
 13
van Gestel A. J, Baty F, RauschOsthof AK, Brutsche MH. Cardiopulmonary and gasexchange responses during the sixminute walk test in patients with chronic obstructive pulmonary disease. Respiration. 2014; 88(4):307–14.
 14
Troosters T, Vilaro J, Rabinovich R, Casas A, Barbera JA, RodriguezRoisin R, Roca J. Physiological responses to the 6min walk test in patients with chronic obstructive pulmonary disease. Eur Respir J. 2002; 20(3):564–9.
 15
Poole DC, Hirai DM, Copp SW, Musch TI. Muscle oxygen transport and utilization in heart failure: implications for exercise (in)tolerance. Am J Physiol Heart Circ Physiol. 2012; 302(5):1050–1063.
 16
Laveneziana P, Valli G, Onorati P, Paoletti P, Ferrazza AM, Palange P. Effect of heliox on heart rate kinetics and dynamic hyperinflation during highintensity exercise in COPD. Eur J Appl Physiol. 2011; 111(2):225–34.
 17
PuenteMaestu L, Sanz ML, Sanz P, Nunez A, Gonzalez F, Whipp BJ. Reproducibility of the parameters of the ontransient cardiopulmonary responses during moderate exercise in patients with chronic obstructive pulmonary disease. Eur J Appl Physiol. 2001; 85(5):434–41.
 18
Hughson RL. Oxygen uptake kinetics: historical perspective and future directions. Appl Physiol Nutr Metab. 2009; 34(5):840–50.
 19
Chiappa GR, BorghiSilva A, Ferreira LF, Carrascosa C, Oliveira CC, Maia J, Gimenes AC, Queiroga F, Berton D, Ferreira EM, Nery LE, Neder JA. Kinetics of muscle deoxygenation are accelerated at the onset of heavyintensity exercise in patients with COPD: relationship to central cardiovascular dynamics. J Appl Physiol. 2008; 104(5):1341–1350.
 20
Cross AM, Higginbotham MB. Oxygen deficit during exercise testing in heart failure. Relation to submaximal exercise tolerance. Chest. 1995; 107(4):904–8.
 21
Sietsema KE, BenDov I, Zhang YY, Sullivan C, Wasserman K. Dynamics of oxygen uptake for submaximal exercise and recovery in patients with chronic heart failure. Chest. 1994; 105(6):1693–1700.
 22
Macdonald M, Pedersen PK, Hughson RL. Acceleration of VO2 kinetics in heavy submaximal exercise by hyperoxia and prior highintensity exercise. J Appl Physiol. 1997; 83(4):1318–1325.
 23
Koga S, Shiojiri T, Shibasaki M, Kondo N, Fukuba Y, Barstow TJ. Kinetics of oxygen uptake during supine and upright heavy exercise. J Appl Physiol. 1999; 87(1):253–60.
 24
Pinheiro J, Bates D. MixedEffects Models in S and SPLUS. Statistics and Computing: Springer; 2000.
 25
Li Y, Jiang L. Fitting logistic growth curve with nonlinear mixedeffects models. Adv J Food Scie Tech. 2013; 5:392–7.
 26
Ritz C, Streibig JC. Bioassay analysis using r. J Stat Soft. 2005; 12(5):1–22.
 27
Pinheiro J, Bates D, DebRoy S, Sarkar D. R Core Team nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1–120. 2015. http://CRAN.Rproject.org/package=nlme.
 28
Gerhard D, Ritz C. 2015. Medrc: Mixed Effect Doseresponse Curves. R package version 0.0–73.
 29
Hothorn T, Bretz F, Westfall P. Simultaneous inference in general parametric models. Biometrical J. 2008; 50(3):346–63.
 30
Vaida F, Blanchard S. Conditional akaike information for mixedeffects models. Biometrika. 2005; 92:351–70.
 31
Verbeke G, Lesaffre E. The effect of misspecifying the randomeffects distribution in linear mixed models for longitudinal data. Comput Stat Data Anal. 1997; 23(4):541–56.
 32
Nielsen OK, Ritz C, Streibig JC. Nonlinear mixedmodel regression to analyze herbicide doseresponse relationships 1. Weed Technol. 2004; 18(1):30–7.
 33
Nielsen OK, Chikoye D, Streibig JC. Efficacy and costs of handheld sprayers in the subhumid savanna for cogongrass control 1. Weed technology. 2005; 19(3):568–74.
 34
Lavielle M, Mentré F. Estimation of population pharmacokinetic parameters of saquinavir in hiv patients with the monolix software. J Pharmacokinet Pharmacodyn. 2007; 34(2):229–49.
Acknowledgements
The study was supported by an unconditional research grant by the Lungenliga St. Gallen and an institutional grant by the Kantonsspital St. Gallen.
Availability of data and materials
The data set supporting the conclusions of this article is provided in Additional file 1.
Authors’ contributions
FB wrote the manuscript, carried out the analysis and was involved in the implementation of the methodology. CR and DG developed the methodology. AVG provided expertise on the physiological aspects of the work. MB designed the experiments, and supervised all research experiments. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Ethics approval and consent to participate
The study was conducted in accordance with the principles enunciated in the Declaration of Helsinki and the guidelines of Good Clinical Practice.
Author information
Additional files
Additional file 1
Supporting oxygen kinetics data set. (RDA 151 kb)
Additional file 2
Sample R code. (R 388 kb)
Additional file 3
Model selection for the random effects structure. (PDF 661 kb)
Additional file 4
Graphical check of the residuals. (EPS 777 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Baty, F., Ritz, C., van Gestel, A. et al. Modeling the oxygen uptake kinetics during exercise testing of patients with chronic obstructive pulmonary diseases using nonlinear mixed models. BMC Med Res Methodol 16, 66 (2016). https://doi.org/10.1186/s1287401601738
Received:
Accepted:
Published:
Keywords
 Nonlinear mixed effects
 Modeling
 Chronic obstructive pulmonary disease
 Exercise testing
 Oxygen kinetics