Explicit and Implicit Solvers for PBPK Systems

Why Stability Matters in a PBPK Integration

A PBPK model is represented as a system of differential equations for its state variables, such as amounts or concentrations in compartments. A numerical solver advances those states in time using finite steps. The solver must do more than produce a plausible value at the end of a step: its numerical updates must not turn a genuinely decaying process into a growing one.

This issue is especially clear when a system contains processes on different time scales. For example, rapid exchange between blood and a tissue can coexist with much slower elimination. Even when the overall PBPK trajectory changes gradually, a numerical method may need a small step to handle the rapid mode. The linear test equation provides a simple way to see why.

For a PBPK system, the test equation is not a replacement for the full model. It describes what can happen to a mode of the system, such as a fast exchange mode, when the equations are linearized. The eigenvalues of the linearized dynamics play the role of \(\lambda\). Stable decaying modes have negative real parts; a mode with a large negative eigenvalue is particularly restrictive for many explicit methods.

Explicit and Implicit Euler: Two Stability Conditions

Explicit Euler evaluates the derivative at the start of a step. Applied to \(y'=\lambda y\), it gives \(y_{n+1}=y_n+h\lambda y_n\), so the amplification factor is \(R(z)=1+z\). Its absolute-stability condition is the disk \(|1+z|\leq 1\), centered at \(-1\) with radius 1. For a negative real \(\lambda\), this reduces to a direct upper bound on \(h\).

Implicit Euler instead evaluates the derivative at the end of the step. Its update is \(y_{n+1}=y_n+h\lambda y_{n+1}\), so \(R(z)=1/(1-z)\). The stability condition is \(|1-z|\geq 1\). Every \(z\) with negative real part satisfies this inequality, so implicit Euler is stable for any positive step on a decaying real mode. It is also stable throughout the left half of the complex plane, a property called A-stability.

Runge–Kutta and Backward Differentiation Methods

Runge–Kutta methods take several derivative evaluations within a step and combine them to advance the state. For the classical explicit fourth-order Runge–Kutta method, the amplification factor on the test equation is

On the negative real axis, this method is stable approximately for \(-2.785\leq z\leq 0\). That interval is wider than explicit Euler's \([-2,0]\), but it is still bounded. A sufficiently fast decaying mode therefore forces an explicit Runge–Kutta solver to use small steps, even if the slower states change little over those steps.

Backward differentiation formulas (BDFs) use the current derivative and several past state values. They are implicit: at each step, the new state must be found from an equation involving that state. BDF1 is implicit Euler. BDF2, for a constant step \(h\), uses

For a BDF method, stability is assessed by checking the roots of its characteristic equation, rather than by a single amplification factor. For BDF2 on the test equation, those roots satisfy \((3-2z)r^2-4r+1=0\). BDF1 and BDF2 are A-stable: they can stably integrate modes with negative real eigenvalues without the explicit methods' finite negative-real-axis step restriction. Higher-order BDF methods do not all retain A-stability; their stability regions cover sectors of the left half-plane rather than the whole half-plane.

Implicit stability does not remove the need to choose a suitable step. A method can remain stable while representing a rapidly changing transient poorly. Nor does stability alone ensure the PBPK solution is physiologically plausible or accurate. The distinction is central: stability concerns unwanted numerical amplification, while accuracy concerns how closely the numerical trajectory approximates the differential-equation solution.

A PBPK System with Fast and Slow Modes

Consider a hypothetical PBPK system whose linearized dynamics can be expressed in modal coordinates as two independent decaying components:

Here \(q_f\) and \(q_s\) are modal amplitudes in mg, and time is in hours. They represent a fast and a slow mode of a coupled system, not necessarily the amount in any one physical compartment. Let \(q_f(0)=q_s(0)=1\ \mathrm{mg}\). This simplified example isolates the numerical effect of widely separated rates.

The example shows the practical trade-off. An explicit method is straightforward to apply, but its stability may impose many small steps to carry a slow process across a long time interval. An implicit BDF method can often take larger stable steps in the presence of rapidly decaying modes, but each step requires solving an implicit equation. For nonlinear PBPK equations, this generally means solving a system of nonlinear equations, which has computational cost and can itself require careful implementation.

Choosing and Checking a Solver

Explicit Runge–Kutta methods are a natural choice when the relevant dynamics allow steps large enough for efficient integration and explicit evaluations are inexpensive. Implicit BDF methods are useful when decaying fast modes constrain explicit steps and a solver can efficiently handle the required implicit solves. Implicit Runge–Kutta methods are another option, but their implicit stage calculations can be more involved. No method is universally best: the equations, required accuracy, and computational cost matter.

Several numerical warning signs merit investigation. Unexpected oscillations or negative compartment amounts can indicate an overly large step or a method poorly suited to the dynamics, though they can also arise from model definitions or inputs. Repeatedly rejected steps or a solver that must take extremely small steps can signal that an explicit method is constrained by a fast mode. A stable-looking trajectory is not proof of accuracy: for example, implicit Euler gave a stable but poor fast-mode value in the worked example. Comparing results after reducing the step, or comparing an explicit and an implicit method at adequate accuracy, can help reveal numerical sensitivity.

These checks do not substitute for verifying the PBPK model itself. A numerical method integrates the equations supplied to it; it cannot correct an incorrect flow coupling, a missing process, or an implausible parameter. Likewise, changing to an implicit solver does not guarantee the solution will preserve every desirable property, such as nonnegative amounts, at every step. Stability regions describe behavior for the linear test equation, and provide guidance—not a universal guarantee for a nonlinear, time-varying PBPK system.

Deriving the Explicit Euler Restriction

The step restriction follows directly by requiring the numerical multiplier to have magnitude no greater than 1. For explicit Euler with a real decaying rate \(\lambda<0\), the update is \(y_{n+1}=(1+h\lambda)y_n\). Absolute stability therefore requires

Subtracting 1 gives \(-2\leq h\lambda\leq 0\). Since \(\lambda<0\) and \(h\geq 0\), dividing by \(\lambda\) reverses the inequality that gives the upper bound, yielding \(0\leq h\leq -2/\lambda=2/|\lambda|\). This is the explicit Euler stability condition. The restriction is more severe for modes with larger \(|\lambda|\), which explains the challenge of explicit integration when a PBPK system has a very fast process alongside a slow one.