Numerical Integration of Feedback Pharmacokinetic Models

Why Numerical Integration Matters

In a feedback pharmacokinetic (PK) model, drug concentration can change enzyme abundance, and enzyme abundance can in turn change drug clearance. The resulting equations are coupled: concentration at one time affects the rate of enzyme change, while the changing enzyme level affects the later concentration profile. Except in special cases, the state values must be approximated by numerical integration of the differential equations.

Numerical integration is not merely a software detail. If the solver takes steps that are too large, it may distort peak and trough concentrations, delay or exaggerate induction, or return predictions that depend noticeably on the solver settings. Conversely, unnecessarily strict settings can increase computation without materially improving the predictions. The goal is to make numerical error small relative to the scientific and statistical precision required for the analysis.

A Coupled PK-Induction System

Consider a hypothetical one-compartment model. Let \(A(t)\) be the drug amount in the central compartment, in mg, and let \(E(t)\) be enzyme abundance relative to its baseline value, so \(E(t)\) is dimensionless and baseline abundance is \(E=1\). Concentration is \(C(t)=A(t)/V\), where \(V\) is the apparent volume in L. Let \(CL_0\) be baseline clearance in L/h, \(k_{\mathrm{out}}\) the enzyme turnover rate in h\(^{-1}\), \(E_{\max}\) the maximum fractional induction, and \(EC_{50}\) the concentration producing half-maximal induction, in mg/L.

Here \(R_{\mathrm{in}}(t)\) is the input rate in mg/h. With no drug present, the enzyme equation returns \(E\) toward 1; with drug present, the concentration-dependent term raises its production target. The clearance at time \(t\) is \(CL_0 E(t)\). This system assumes the stated one-compartment disposition, a well-mixed concentration, the specified induction function, and nonnegative concentration. It is an illustrative model, not a claim that every induction mechanism follows this form.

For an NLME analysis, the same numerical problem is solved for each subject’s parameter values and dosing history. A numerical solution supplies predicted concentrations at observation times; the likelihood or estimation procedure then uses those predictions. Solver error can therefore affect model predictions even when the statistical estimation algorithm itself converges.

What the Solver Tolerances Control

An adaptive solver estimates the error associated with a proposed step and adjusts its step size. A common conceptual acceptance scale for state component \(y_j\) is

The relative tolerance, \(\mathrm{rtol}\), controls error relative to the magnitude of a state. The absolute tolerance, \(\mathrm{atol}_j\), provides a scale near zero and has the units of state \(y_j\). With multiple states, absolute tolerances should reflect their units and scientifically relevant scales: an amount in mg and a dimensionless enzyme state should not automatically be assigned the same absolute tolerance. Implementations differ in how they combine componentwise error estimates, so this expression is a useful interpretation, not a guarantee of a particular solver’s exact acceptance rule.

Tolerances usually control estimated local integration error, not a direct bound on the total error at every later time. Accumulated error can depend on the dynamics, the solver method, and discontinuities. A reliable practice is therefore to repeat the calculation with stricter tolerances or a smaller maximum step and compare the predictions that matter, such as concentrations at observation times and the timing or magnitude of induction.

Step Size, Accuracy, and a Simple Derivation

A solver’s step size is the time interval advanced in one numerical step. To see why it matters, consider the scalar elimination equation \(dy/dt=-ky\), with \(k>0\), and initial amount \(y(0)=y_0\). Its exact solution is \(y(t)=y_0e^{-kt}\). The forward Euler method advances by

where \(h\) is the step size. The exact one-step multiplier is \(e^{-kh}\), while Euler’s multiplier is \(1-kh\). Taylor expansion gives \(e^{-kh}=1-kh+(kh)^2/2+O(h^3)\); thus, for small \(h\), the one-step discrepancy is of order \(h^2\). Over a fixed time interval, the accumulated error of this first-order method is generally of order \(h\), provided the solution is sufficiently smooth and the method is operating in its stable range.

This argument explains refinement, but does not mean that every solver uses Euler’s method. Practical solvers often use higher-order methods and estimate errors by comparing approximations of different orders. Nor does reducing \(h\) guarantee a useful improvement if the error is dominated by incorrectly specified parameters, events, or model structure.

Stiffness and Numerical Instability

A system is called stiff when it contains substantially different time scales, often with a rapidly decaying mode alongside a much slower mode. In the induction model, drug disposition may change over minutes or hours while enzyme turnover changes over days. A solver may then need small steps to represent the fast dynamics or maintain numerical stability, even during periods when the overall trajectory appears to change slowly.

For the scalar decay equation \(dy/dt=-ky\), forward Euler is stable only when its amplification factor satisfies \(|1-kh|<1\), which gives \(0<kh<2\). If \(k=10\) h\(^{-1}\) and \(h=0.25\) h, the factor is \(1-10(0.25)=-1.5\): repeated steps alternate sign and grow in magnitude, although the exact solution decays. A real PK state cannot become negative under simple elimination, so such behavior is a clear numerical warning. In a coupled system, local eigenvalues of the Jacobian characterize the modes; the scalar condition is an illustration, not a complete stability test for every solver and system.

Explicit Runge–Kutta methods can be efficient for nonstiff systems. When stiffness forces extremely small explicit steps, an implicit method, such as a BDF or implicit Runge–Kutta method, may be more efficient because it can remain stable at larger steps. Implicit methods require more work per step and are not automatically more accurate in every application. Compare methods using convergence checks and relevant predictions, rather than assuming that one solver class is always preferable.

Dosing Events Must Be Represented at the Right Time

Bolus doses create jumps in an amount state, whereas the continuous ODE describes its evolution between doses. If a bolus of size \(D\) is given at time \(t_d\), the amount satisfies \(A(t_d^+)=A(t_d^-)+D\); the enzyme state remains continuous unless the model explicitly specifies an instantaneous enzyme change. The solver should integrate to \(t_d\), apply the jump, and continue from the post-dose state. Infusion starts and stops do not generally create jumps in amount, but they change \(R_{\mathrm{in}}(t)\), so the integrator must account for the change in the right-hand side at the scheduled time.

Diagnostics and Practical Choices

A well-behaved numerical solution should be reproducible to the required precision under reasonable solver refinement. Inspect more than whether the program returned a successful status. Useful checks include concentration predictions at observation times, peak and trough values, enzyme trajectories, and values immediately before and after each event. Unexpected negative amounts, oscillations, abrupt changes in a continuous state, or predictions that shift materially with tighter tolerances are warning signs.

A maximum step size can be useful when a short-lived feature or rapidly changing input needs resolution, even if the adaptive error estimate would otherwise permit larger steps. Requested output times alone do not necessarily force a solver to take an internal step at each time. Event times must be supplied or otherwise handled explicitly. Solver refinement is also not a remedy for structural model inadequacy, poor identifiability, or inappropriate residual assumptions.

Common Pitfalls and Limitations

• Using one absolute tolerance for states with different units. A numerical value suitable for an amount in mg may be too loose or too strict for a dimensionless enzyme state. Set scales by state.
• Assuming a successful solver return proves accuracy. Successful completion means the algorithm met its own numerical criteria, not that the global prediction error is negligible. Check refinement.
• Treating a dose as a smooth input when it is a bolus. A bolus changes the amount discontinuously. Apply the jump at the prescribed time rather than smearing it over a solver step.
• Interpreting solver tolerances as measurement precision. Tolerances concern numerical approximation of the ODE solution; they are not the residual error model or confidence in a parameter.
• Reducing the step without checking stability or cost. Very small steps may be inefficient, and explicit solvers can struggle with stiff systems. A stiff solver may help, but its results still require validation.

Numerical integration has finite cost and finite accuracy. Tight tolerances can increase the time needed for fitting many subjects, while loose settings can perturb the objective function used by estimation. The appropriate settings depend on the model scales, time scales, events, and the required precision. No single tolerance or solver is universally correct; demonstrate that the resulting predictions are stable enough for the intended analysis.