Tutorials › Population Pharmacokinetics › Forecasting Under Nonlinear Pharmacokinetics

Therapeutic Drug Monitoring and Individual Predictions · Tutorial 559 of 599

Forecasting Under Nonlinear Pharmacokinetics

Use a concentration-dependent elimination model to calculate patient-specific forecasts and assess how nonlinear likelihoods and sparse observations shape forecast uncertainty.

Advanced 12 min read

What You'll Learn

  • Derive the concentration-time relationship for Michaelis–Menten elimination after a dose.
  • Calculate a future concentration by numerically solving the nonlinear model equation.
  • Construct a patient-specific posterior using a nonlinear concentration likelihood.
  • Translate posterior parameter uncertainty into uncertainty in a future concentration.
  • Explain why sparse observations may not distinguish maximum elimination rate from the Michaelis constant.
  • Diagnose weak identification and compare nonlinear forecasts with linearized approximations.

Why Nonlinear Elimination Changes a Forecast

With concentration-independent clearance, the rate of concentration decline is proportional to concentration, and the resulting concentration-time curve has an exponential form. Under saturable elimination, the rate depends on concentration through a nonlinear function. The same drug can therefore decline relatively slowly at high concentrations, where elimination is near its maximum rate, and more rapidly in proportional terms at low concentrations. A patient forecast must respect that changing relationship rather than extend a single exponential decline.

This tutorial uses a one-compartment Michaelis–Menten elimination model. The parameters are patient-specific: the maximum elimination rate and the concentration at which the elimination rate is half of that maximum. When those parameters are uncertain, a patient’s observations update their distribution through a nonlinear likelihood, and the forecast must carry the resulting uncertainty forward.

Definition: Let \(C(t)\) be concentration at time \(t\), \(V_{\max}>0\) the maximum elimination rate, and \(K_m>0\) the Michaelis constant. With no input during an elimination interval, the model is $$ \frac{dC(t)}{dt}=-\frac{V_{\max}C(t)}{K_m+C(t)}. $$ If \(C\) is measured in mg/L and \(t\) in hours, then \(V_{\max}\) has units mg/L/h and \(K_m\) has units mg/L. At \(C=K_m\), the elimination rate has magnitude \(V_{\max}/2\).

The forecast below assumes a well-mixed one-compartment system, a known concentration at the start of the interval, and constant \(V_{\max}\) and \(K_m\) throughout it. There is no input during the interval being forecast; a known dose can instead be represented as an instantaneous concentration jump, followed by a new elimination interval. The model does not itself account for ongoing absorption, distribution between compartments, changing parameters, or concentration measurement error. Those must be represented separately if they are relevant.

From the Elimination Equation to a Forecast

Separate variables and integrate from \(C_0\), the concentration at the start of the interval, to \(C\), the concentration after \(t\) time units. Since \(dC/dt=-V_{\max}C/(K_m+C)\), rearrangement gives \(dt=-(K_m+C)dC/(V_{\max}C)\). Integration yields the implicit concentration-time equation:

$$ V_{\max}t=C_0-C+K_m\ln\left(\frac{C_0}{C}\right). $$

For positive \(C_0\), \(V_{\max}\), and \(K_m\), the right-hand side is a continuous, strictly decreasing function of \(C\) on \(0<C\leq C_0\): its derivative with respect to \(C\) is \(-1-K_m/C<0\). It is zero at \(C=C_0\) and grows without bound as \(C\) approaches zero. Therefore, for any \(t>0\), the equation has exactly one solution between zero and \(C_0\). In general that solution is found numerically, not by treating the model as an exponential.

Worked Example: Forecasting With Known Nonlinear Parameters

Hypothetical data. Suppose a patient has concentration \(C_0=12.0\) mg/L at the start of an elimination interval. Use \(V_{\max}=4.00\) mg/L/h and \(K_m=2.00\) mg/L to forecast concentration after 2.00 h. The assumptions above are taken to hold for this interval.

Substitution into the integrated equation gives

$$ 8.00=12.0-C+2.00\ln\left(\frac{12.0}{C}\right). $$

Define \(g(C)=12.0-C+2.00\ln(12.0/C)-8.00\). For example, \(g(5.54)\) is approximately \(0.006\), while \(g(5.55)\) is approximately \(-0.008\). The unique root lies between these values. Numerical solution gives \(C\approx5.544\) mg/L. Thus the forecast is about 5.54 mg/L after 2.00 h, conditional on the stated parameter values and model assumptions.

The nonlinear term matters: the elimination rate is not constant over the interval, nor is the fractional decline constant. This is a point forecast; it does not include uncertainty about the patient’s parameters or future measurement error.

Patient-Specific Nonlinear Likelihoods

Let \(\theta=(V_{\max},K_m)\) denote the patient’s nonlinear parameter vector, and let \(y_j\) be a measured concentration at time \(t_j\). Given the patient’s dose history, let \(C(t_j;\theta)\) be the model-predicted concentration. A common observation model assumes independent additive normal errors with known standard deviation \(\sigma\):

$$ y_j\mid\theta\sim N\left(C(t_j;\theta),\sigma^2\right), \qquad L(\theta)=\prod_{j=1}^{m} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{[y_j-C(t_j;\theta)]^2}{2\sigma^2}\right\}. $$

Here \(L(\theta)\) is the likelihood for the patient’s \(m\) observations. This formulation assumes the observation errors are independent, have the stated constant variance, and are additive on the concentration scale. If the assay error is better represented by a different model, that likelihood should replace the normal one; the nonlinear concentration predictions remain inside the likelihood either way.

With a population-informed prior density \(p(\theta)\), the patient-specific posterior is proportional to \(L(\theta)p(\theta)\). The posterior is generally not normal: the model prediction is nonlinear in both parameters, and a sparse data set can leave a curved region of plausible parameter combinations. For a future time \(t_f\), calculate \(C(t_f;\theta)\) over posterior parameter values. The posterior mean of these conditional forecasts summarizes their centre; their spread measures parameter-related uncertainty. A full posterior predictive distribution also includes future observation error if the target is a future measured concentration rather than the underlying model concentration.

Worked Example: Updating a Sparse Patient Forecast

Hypothetical data. Suppose a population prior places equal probability on two plausible parameter combinations, labelled A and B. Numerical evaluation of the nonlinear model gives predicted concentration 8.000 mg/L at the observed time under A and 9.1176 mg/L under B. The observed concentration is \(y=8.000\) mg/L, and the additive normal observation error has known \(\sigma=0.500\) mg/L. For a future time, the same parameter combinations give forecasts of 6.87411 mg/L and 5.67143 mg/L, respectively. These are illustrative posterior candidates, not estimates from a real study.

For equal prior probabilities, the posterior odds are the likelihood ratio. The normalizing constants cancel:

$$ \frac{L(A)}{L(B)} = \exp\left( \frac{(y-9.1176)^2-(y-8.000)^2}{2(0.500)^2} \right) \approx 12.16. $$

Thus the posterior probability of A is \(12.16/(1+12.16)\approx0.924\); that of B is approximately 0.076. The posterior mean of the future model concentration is

$$ (0.924)(6.87411)+(0.076)(5.67143) =6.782706\approx6.783\text{ mg/L}. $$

For two possible values \(a\) and \(b\) with probabilities \(w\) and \(1-w\), the variance is \(w(1-w)(a-b)^2\). Applying this identity gives

$$ \sqrt{(0.924)(0.076)(6.87411-5.67143)^2} \approx0.3187\text{ mg/L}. $$

The observation shifts probability toward A, reducing forecast uncertainty compared with the equal-prior spread, but it does not eliminate uncertainty. The numbers summarize uncertainty across these two candidate parameter combinations only; a continuous posterior would ordinarily be summarized using many posterior draws.

Why Sparse Observations May Not Identify Both Parameters

The two parameters can trade off. At concentrations much smaller than \(K_m\), \(K_m+C\) is approximately \(K_m\), so the model’s rate is approximately proportional to \(V_{\max}/K_m\). Data restricted to that range may constrain the ratio more strongly than either parameter separately. At concentrations much larger than \(K_m\), the rate approaches \(V_{\max}\), providing more direct information about the maximum rate. Observations that span different concentration ranges can therefore help distinguish the parameters, although the actual information also depends on timing and measurement precision.

This approximation follows directly by dividing numerator and denominator by \(K_m\):

$$ \frac{V_{\max}C}{K_m+C} = \frac{V_{\max}}{K_m}C \frac{1}{1+C/K_m} \approx \frac{V_{\max}}{K_m}C \quad\text{when } C/K_m\text{ is small}. $$

Worked Example: Seeing the Effect of Sampling Concentration

Hypothetical parameter comparison. Consider two plausible parameter pairs: A has \(V_{\max}=4.00\) mg/L/h and \(K_m=2.00\) mg/L; B has \(V_{\max}=8.00\) mg/L/h and \(K_m=4.00\) mg/L. Both have \(V_{\max}/K_m=2.00\) h\(^{-1}\). Compare their instantaneous elimination rates at low concentration \(C=0.100\) mg/L and high concentration \(C=20.0\) mg/L.

At \(C=0.100\) mg/L, the rate magnitudes are \(4.00(0.100)/(2.00+0.100)=0.1905\) mg/L/h for A and \(8.00(0.100)/(4.00+0.100)=0.1951\) mg/L/h for B. These are close, as expected when both parameter pairs share the same ratio and concentration is low relative to \(K_m\).

At \(C=20.0\) mg/L, the corresponding magnitudes are \(4.00(20.0)/(2.00+20.0)=3.636\) mg/L/h and \(8.00(20.0)/(4.00+20.0)=6.667\) mg/L/h. The difference is now substantial. Measurements informative about the high-concentration decline could separate these parameter pairs more effectively than low-concentration observations alone. These instantaneous rates illustrate the model’s sensitivity; actual parameter inference must use the predicted concentrations at the observation times and their likelihood.

Diagnostics, Alternatives, and Limitations

A nonlinear forecast should be checked against both the data and the parameter uncertainty. Useful warning signs include a posterior that remains broad or strongly correlated in \(V_{\max}\) and \(K_m\), forecasts that vary greatly across plausible parameter values, and systematic differences between observed concentrations and predictions at particular times or concentration ranges. Numerical optimization or sampling problems can also produce apparent uncertainty that is an artifact of failed computation; numerical stability is a separate concern. When a posterior is weakly informed by the patient’s data, its dependence on the population prior should be made explicit rather than hidden by reporting a single precise-looking estimate.

A linearized approximation can be convenient when parameter uncertainty is small: locally approximate the nonlinear forecast by a first-order expansion around a parameter estimate, then propagate uncertainty through the local derivatives. Its limitation is that it represents curvature poorly when uncertainty is broad, the posterior is skewed, or plausible parameter values lie along a curved ridge. Direct posterior calculation, using numerical integration or posterior draws, retains the nonlinear mapping from parameters to concentrations and is generally preferable when that curvature matters. Neither approach repairs a misspecified structural model or an inappropriate observation-error model.

A common mistake is to fit a straight line to concentration against time, or to extrapolate an exponential decline, simply because a few observed concentrations appear nearly linear on a chosen scale. That reasoning assumes a concentration-independent fractional decline; Michaelis–Menten elimination does not generally have that property. Use the specified nonlinear elimination model for each parameter value, inspect whether the observation range informs both parameters, and communicate uncertainty when it does not.

Key takeaway: A patient-specific forecast under nonlinear elimination is a distribution of model-based future concentrations across plausible patient parameters, not merely a single curve through the observed concentrations. Sparse data may update that distribution without identifying \(V_{\max}\) and \(K_m\) separately.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

  1. Under the stated model, what is the elimination rate magnitude at \(C=K_m\)?
  2. Why can observations only at concentrations far below \(K_m\) provide weak information about \(V_{\max}\) and \(K_m\) separately?
  3. For the two-candidate example, what is the posterior probability of B after the observation?
  4. Name one assumption of the additive normal likelihood that may need to be reconsidered for a different assay-error pattern.

Solutions. (1) Substituting \(C=K_m\) gives \(V_{\max}K_m/(K_m+K_m)=V_{\max}/2\). (2) At low concentration the rate is approximately \((V_{\max}/K_m)C\), so data may mainly identify the ratio rather than the two parameters separately. (3) The probabilities sum to one, so B has posterior probability \(1-0.924=0.076\). (4) The likelihood assumes independent additive normal errors with constant standard deviation; for example, an error pattern whose variability changes with concentration would not satisfy the constant-variance assumption.

Clinical Biostats