### Participants

Children with cardiopulmonary disease attending outpatient clinics at IWK Health Centre in Halifax, Canada were recruited to the study. Healthy control children were recruited from friends and relatives of hospital personnel, or siblings of these patients. Data were collected on 100 pediatric subjects with 4 to 12 ratings per subject. For reliable parameter estimation the data were limited to 79 children with at least 6 data points. The study was approved by the Research Ethics Board of the IWK Health Centre. Assent was obtained from all participants. Mature minors, or parents of younger children, signed informed consent.

### Procedure

Subjects performed continuous, graded, maximal, cycle ergometer (WE Collins) test exercise employing step increments of either, 50, 100, or 150 kpm per minute depending on size and age. Increments were chosen to achieve test duration of 6‐10 minutes ‐ until voluntary, symptom‐limited, exhaustion occurred as previously described [6]. Borg scale ratings were obtained during each 1‐minute increment. Work was normalized across subjects by expressing it as fraction of individual maximum work capacity (Wmax).

Both dyspnea and perceived leg exertion were measured, but we report only Borg ratings for perceived leg exertion since our principal aim was modeling the stimulus‐response function. We employed the Dalhousie pictorial scales [11] and the Borg scale [12], chosen because it had been used in previous investigations in adults with similar aims [13, 14]. Subjects were first given an explanation of the scale by the research assistant.

The Borg scale of perceived exertion ranges from 0 to 10 with verbal descriptors 0= nothing at all, 0.5=very,very slight, 1=very slight, 2=slight, 3=moderate, 4=somewhat severe, 5=severe, 7=very severe, 9=very, very severe, 10=maximal [12, 15]. The scale was mounted on a large clipboard in front of the participant, and had a sliding cursor located on the left margin. The cursor was moved manually from the top downward by the research assistant until it pointed to the rating that best described the subject’s degree of leg fatigue, at which point he or she activated a bell mounted on the handlebars.

### Models

Perceived leg exertion of individual children were modeled in terms of fraction of maximum work capacity. Two models have been proposed in the literature, the power model (P) which has been used for fitting individual curves in adults [10], and the power delay model (PD) [9] which has not been used for fitting study data to the best of our knowledge. Since an increase in leg fatigue was observed after some delay, a lag or delay model (D) was developed for this study. A “lag” or delay was defined as % Wmax at which point a clear increase in ratings of leg fatigue occurred. While the power model accounted for curvature of the model fits, but not for a delay, and the delay model accounted for a delay, but not a curvature, we also considered a quadratic‐delay model (QD) and a power model with delay (PD). The quadratic delay model includes an intercept, a linear term, and a quadratic term. This allows for flexibility that the observed pattern may be linear without a curvature. The parameters in each model were estimated using a quasi‐Newton method with box constraints, where a variable can be given a lower bound [16]. For the purpose of estimating the parameters constraints were introduced, namely the delay is constrained to be larger than the initial observed % maximum work capacity, the exponent must be positive. Since the ratings increase, and a minimum proportion of work capacity must be observed, these are reasonable assumptions. The four models were defined as follows. For each child (*i*=1,…,*n*) with *j*=1,…,*n*
_{
i
} observations, leg exertion (*Y*=(*y*
_{
ij
})) is modeled as a function of % maximum work capacity (*X*=(*x*
_{
ij
})).

\begin{array}{ll}\phantom{\rule{6.8em}{0ex}}\text{Power model (P)}:{y}_{\mathit{\text{ij}}}\hfill & ={a}_{i}+{b}_{2i}{x}_{\mathit{\text{ij}}}^{{d}_{i}}+{\epsilon}_{\mathit{\text{ij}}}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{7em}{0ex}}\text{Delay model (D)}:{y}_{\mathit{\text{ij}}}\hfill & ={a}_{i}+{b}_{1i}({x}_{\mathit{\text{ij}}}-{c}_{i}\vee 0)+{\epsilon}_{\mathit{\text{ij}}}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{3em}{0ex}}\text{Power\u2010delay model (PD)}:{y}_{\mathit{\text{ij}}}\hfill & ={a}_{i}+{b}_{2i}{({x}_{\mathit{\text{ij}}}\phantom{\rule{0.3em}{0ex}}-\phantom{\rule{0.3em}{0ex}}{c}_{i}\vee 0\phantom{\rule{1em}{0ex}})}^{{d}_{i}}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}{\epsilon}_{\mathit{\text{ij}}}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{0.8em}{0ex}}\text{Quadratic\u2010delay model (QD)}:{y}_{\mathit{\text{ij}}}\hfill & ={a}_{i}+{b}_{1i}({x}_{\mathit{\text{ij}}}-{c}_{i}\vee \phantom{\rule{0.3em}{0ex}}0)\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{2.5em}{0ex}}+{b}_{2i}{({x}_{\mathit{\text{ij}}}-{c}_{i}\vee \phantom{\rule{0.3em}{0ex}}0)}^{2}+{\epsilon}_{\mathit{\text{ij}}}\phantom{\rule{2em}{0ex}}\end{array}

where {a}_{i},{b}_{1i},{b}_{2i},{d}_{i}\in \mathbb{R},{d}_{i}\ge 0, and *c*
_{
i
}> min*j*(*x*
_{
ij
}). For all models it is assumed that {\epsilon}_{\mathit{\text{ij}}}\sim \mathcal{N}(0,{\sigma}_{i}^{2}).

#### Model assessment

For model assessment we calculated root‐mean‐square error RMSE{}_{i}=\sqrt{\sum _{j=1}^{{n}_{i}}{({y}_{\mathit{\text{ij}}}-{\u0177}_{\mathit{\text{ij}}})}^{2}/{n}_{i}}. In place of the Akaike information criterion (AIC) we used the corrected AIC (AICc), which imposes a penalty on the number of parameters for finite or sparse samples. The Bayesian information criterion (BIC) was also calculated. These quantities are defined as

\begin{array}{ll}\phantom{\rule{1em}{0ex}}{\text{AIC}}_{i}& =2k-log\left({L}_{i}\right)\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{1em}{0ex}}{\text{AICc}}_{i}& ={\text{AIC}}_{i}+\frac{2k(k+1)}{{n}_{i}-k-1}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{1em}{0ex}}{\text{BIC}}_{i}& =klog{n}_{i}-{n}_{i}log\left({L}_{i}\right)\phantom{\rule{2em}{0ex}}\end{array}

where *L*
_{
i
} is the maximum of the individual likelihood function. The number of coefficients are *k*=3,3,4 and 4 for the four models. For the purpose of RMSE comparisons (Figure 1) one subject with ratings at only two levels was removed (cluster 9 in Figure 2).

#### Mixed effects models

Nonlinear mixed effects models can be used to study average parameters for delay and growth. The models are similar to the PD and QD models but with additional terms to model overall variation and individual variation (MPD or MQD). This assumes that the delay is constant across all individuals. More general mixed effects models include the delay as a random effect, a mixed effects quadratic‐delay model with varying delays for each subject (MPDV or MQDV).

\begin{array}{ll}\phantom{\rule{0.8em}{0ex}}\text{MQD}:{y}_{\mathit{\text{ij}}}& =a+{b}_{1}({x}_{\mathit{\text{ij}}}\phantom{\rule{0.3em}{0ex}}-\phantom{\rule{0.3em}{0ex}}c\vee \phantom{\rule{0.3em}{0ex}}0)+{b}_{2}{({x}_{\mathit{\text{ij}}}-c\vee \phantom{\rule{0.3em}{0ex}}0)}^{2}+{u}_{i}+{\epsilon}_{\mathit{\text{ij}}},\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{0.8em}{0ex}}\text{MPD}:{y}_{\mathit{\text{ij}}}& =a+{b}_{2}{({x}_{\mathit{\text{ij}}}\phantom{\rule{0.3em}{0ex}}-\phantom{\rule{0.3em}{0ex}}c\vee \phantom{\rule{0.3em}{0ex}}0)}^{d}+{u}_{i}+{\epsilon}_{\mathit{\text{ij}}}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{0.8em}{0ex}}\text{MQDV}:{y}_{\mathit{\text{ij}}}& =a+{b}_{1}({x}_{\mathit{\text{ij}}}\phantom{\rule{0.3em}{0ex}}-\phantom{\rule{0.3em}{0ex}}{c}_{i}\vee \phantom{\rule{0.3em}{0ex}}0)+{b}_{2}{({x}_{\mathit{\text{ij}}}-{c}_{i}\vee \phantom{\rule{0.3em}{0ex}}0)}^{2}\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}{u}_{i}+{\epsilon}_{\mathit{\text{ij}}}\phantom{\rule{2em}{0ex}}\\ \phantom{\rule{0.8em}{0ex}}\text{MPDV}:{y}_{\mathit{\text{ij}}}& =a+{b}_{2}{({x}_{\mathit{\text{ij}}}\phantom{\rule{0.3em}{0ex}}-\phantom{\rule{0.3em}{0ex}}{c}_{i}\vee \phantom{\rule{0.3em}{0ex}}0)}^{d}+{u}_{i}+{\epsilon}_{\mathit{\text{ij}}}\phantom{\rule{2em}{0ex}}\end{array}

where a,{b}_{1},{b}_{2},d\in \mathbb{R}, and *c*∈[ min*i*,*j*(*x*
_{
ij
}), max*i*,*j*(*x*
_{
ij
})] for MQD and MPD, or *c*
_{
i
}=*c*+*γ*
_{
i
}, *c*
_{
i
}∈[ min*j*(*x*
_{
ij
}), max*j*(*x*
_{
ij
})] for MQDV and MPDV.

The overall variation is modeled by \u03f5\sim \mathcal{N}(0,{\sigma}^{2}), and the random intercept by {u}_{i}\sim \mathcal{N}(0,{\tau}^{2}). For MPDV and MQDV the random delay is assumed to follow a normal distribution {\gamma}_{i}\sim \mathcal{N}(0,{\nu}^{2}). All random components are assumed to be independent. Parameters are estimated using a Bayesian approach by choosing conjugate priors for *a*,*b*
_{1},*b*
_{2},*σ*
^{2},*τ*
^{2}, and non‐informative uniform priors for the delay *c* and the exponent *d* (see Additional file 1). A dispersed prior density was chosen for the distribution of the parameters to allow for more data‐driven estimators.

These four competing models were implemented, the mixed effects quadratic‐delay model, the mixed effects power‐delay model, each with common delay, and both models with delay as random effects. Under each implementation three Monte Carlo Markov Chains were run. Convergence was well achieved after 10,000 iterations, namely when Potential Scale Reduction Factor \sqrt{R}<1.2 for all parameters [17]. Further 1,000 iterations were sampled every 10 iterations to obtain a total of 300 posterior samples for all parameters and random effects. For model comparisons the Deviance Information Criterion (DIC) for mixed‐effects models based on the complete likelihood was used [18]. A smaller DIC indicates a better predictive power of the model.

### Computing environment

For estimating the model parameters of models P, D, PD, QD the function *optim* was used in the statistical software R 2.15.3 [19]. The nonlinear mixed effects models were fitted using MATLAB R2013a (Mathworks Inc, Natick Massachusetts).