Why Numerical Implementation Matters
An effect-compartment model links plasma concentration to effect-site concentration through a differential equation. Although closed-form solutions are available for simple inputs, real computations often use plasma concentrations that vary over time and are supplied at discrete times. A numerical solver approximates the effect-site trajectory from that input. Its predictions can depend on the solver, the integration step size, and the error tolerances—not only on the pharmacometric parameters.
The aim is not to select the most sophisticated solver automatically. It is to implement the stated model faithfully, make the numerical error small relative to the question being asked, and check the result against cases with known answers. The earlier material on interpolating plasma concentrations addresses how to define the input between observation times; here, assume that a continuous input \(C_p(t)\) has already been specified.
A numerical method advances an approximate value of \(C_e\) from one time to the next. At each step, it evaluates the derivative using the current or intermediate values of \(C_p\) and \(C_e\). The input must be evaluated consistently with its specified interpolation, and the time grid must use units compatible with \(k_{e0}\). A unit mismatch can produce a plausible-looking trajectory that is wrong by orders of magnitude.
Solvers and a Useful Analytic Benchmark
For a step of length \(h\), forward Euler uses the slope at the start of the step:
The classical fourth-order Runge–Kutta method (RK4) evaluates the derivative at four points within each step and combines those slopes. It is usually more accurate than Euler at the same step size for a smooth input, but it still has numerical error. Adaptive solvers estimate error by comparing approximations and adjust step size to meet specified tolerances.
A simple exact solution provides a benchmark for all these implementations. If \(C_p(t)=C_p\) is constant from \(t_0\), multiply the differential equation by the integrating factor \(e^{k_{e0}t}\). The product rule gives
Integrating from \(t_0\) to \(t\), then dividing by \(e^{k_{e0}t}\), gives
This derivation applies whether or not the initial effect-site concentration equals the plasma concentration. In particular, if \(C_e(t_0)=C_p\), the formula gives the constant solution \(C_e(t)=C_p\); no division by \(C_e-C_p\) or logarithm of zero is needed.
For a linearly changing input, let \(C_p(t_0+u)=C_{p,0}+mu\), where \(m\) is the plasma-concentration slope and \(0\leq u\leq h\). Substitution into the ODE shows that \(C_p(t_0+u)-m/k_{e0}\) is a particular solution: its derivative is \(m\), and \(k_{e0}[C_p-(C_p-m/k_{e0})]=m\). The difference between any solution and that particular solution therefore decays exponentially. Consequently,
This second benchmark checks not only the solver but also how a continuously changing input is passed to it.
Worked Example: Checking Fixed-Step Euler Against the Exact Solution
Suppose a hypothetical model has constant \(C_p=10\) mg/L, \(C_e(0)=0\) mg/L, and \(k_{e0}=0.5\) h\(^{-1}\). Compare Euler estimates at 2 hours for step sizes of 1 hour and 0.5 hours with the exact answer. The input is constant, the stated initial condition is used, and all step sizes are in hours, so the constant-input benchmark applies.
With \(h=1\) h, \(hk_{e0}=0.5\). The Euler updates are \(C_e(1)=0+0.5(10-0)=5\) mg/L and \(C_e(2)=5+0.5(10-5)=7.5\) mg/L.
With \(h=0.5\) h, each update multiplies the remaining gap to 10 mg/L by \(1-0.25=0.75\). After four steps, \(C_e(2)=10[1-(0.75)^4]=6.8359\) mg/L.
The exact value is \(C_e(2)=10[1-e^{-0.5(2)}]=10(1-e^{-1})=6.3212\) mg/L. The absolute errors are \(7.5-6.3212=1.1788\) mg/L for the 1-hour step and \(6.8359-6.3212=0.5147\) mg/L for the 0.5-hour step. Reducing the step improves this example, but neither Euler result should be mistaken for the exact trajectory.
Step Size, Stability, and Tolerances
For constant plasma concentration, Euler updates the deviation from equilibrium according to \(C_{e,n+1}-C_p=(1-hk_{e0})(C_{e,n}-C_p)\). Thus the numerical deviation shrinks in magnitude only when \(|1-hk_{e0}|<1\), or \(0<hk_{e0}<2\). If \(1<hk_{e0}<2\), it alternates around equilibrium while shrinking. For a non-oscillating update in this constant-input case, require \(0<hk_{e0}\leq1\). These are stability and shape checks for Euler, not a guarantee of accuracy; a stable step can still be too coarse.
An adaptive solver typically controls an estimated local error using a scale such as \(\mathrm{atol}+\mathrm{rtol}|C_e|\), where atol is an absolute tolerance in concentration units and rtol is dimensionless. The solver accepts or rejects steps and changes their size to keep its estimated local error within that scale. The tolerances do not directly guarantee that the total error over an entire simulation is below the same threshold. Check convergence by tightening tolerances and confirming that predictions of interest change negligibly.
Worked Example: Using Step Doubling to Assess Euler Error
Consider a hypothetical constant input \(C_p=8\) mg/L, with \(C_e(0)=2\) mg/L and \(k_{e0}=0.4\) h\(^{-1}\), over 1 hour. The inputs satisfy the constant-input assumptions, so the exact result is \(8+(2-8)e^{-0.4}=8-6(0.670320)=3.9781\) mg/L.
A single Euler step of 1 hour gives \(2+1(0.4)(8-2)=4.4\) mg/L. Instead take two half-hour steps. The first gives \(2+0.5(0.4)(6)=3.2\) mg/L; the second gives \(3.2+0.5(0.4)(8-3.2)=4.16\) mg/L. The difference between the one-step and two-step estimates is \(4.4-4.16=0.24\) mg/L.
For a first-order method such as Euler, this step-doubling difference estimates the error scale of the refined result, though it is not an exact error bound. If a target error scale were 0.01 mg/L absolute plus 0.001 times the concentration, its value near 4.16 mg/L would be \(0.01+0.001(4.16)=0.01416\) mg/L. The estimated difference, 0.24 mg/L, is much larger; the step should be reduced and the comparison repeated. In either case, the discrepancy from the exact 3.9781 mg/L confirms that these coarse Euler calculations are not yet sufficiently accurate for that target.
Changing Inputs and Numerical Verification
When plasma concentration is piecewise linear, use the prescribed interpolation within each interval rather than holding the value at the left endpoint for the whole step. If the input is supplied at knots, it is prudent to end integration steps at those knots so the solver does not treat a change in slope as though it were smooth across a step. If the input has a jump, split integration at the jump and apply the specified pre- and post-jump input on the appropriate sides. The effect-site state remains governed by its ODE; a plasma jump does not itself imply an instantaneous jump in \(C_e\).
Worked Example: Verifying RK4 for a Linear Plasma Input
For a hypothetical interval of 2 hours, let \(C_p\) rise linearly from 2 to 6 mg/L, so \(m=(6-2)/2=2\) mg/L/h. Let \(C_e(0)=1\) mg/L and \(k_{e0}=0.5\) h\(^{-1}\). The input is continuous and linear over the interval, with consistent units, so the linear-input solution applies.
The exact result is \(2+2(2-1/0.5)+[1-2+2/0.5]e^{-0.5(2)} =2+3e^{-1}=3.1036\) mg/L.
For RK4, subtract the particular solution \(C_p-m/k_{e0}\) from \(C_e\). The resulting deviation follows \(x'=-0.5x\) and starts at \(x(0)=1-[2-2/0.5]=3\) mg/L. For a step with \(z=k_{e0}h\), RK4 multiplies this deviation by \(R(-z)=1-z+z^2/2-z^3/6+z^4/24\). With \(h=1\) h, \(z=0.5\), so \(R(-0.5)=0.6067708\). After two steps, \(C_e(2)=2+3(0.6067708)^2=3.1045\) mg/L. With \(h=0.5\) h, \(z=0.25\), \(R(-0.25)=0.7788086\); after four steps, \(C_e(2)=2+3(0.7788086)^4\approx3.1037\) mg/L. Both approach the analytic value, and the smaller-step result is closer in this calculation. The benchmark also verifies that the RK4 stages use the changing plasma input rather than a single held value throughout each interval.
Practical Checks, Alternatives, and Limitations
A reliable implementation benefits from several complementary checks. First, verify units, the sign of the derivative, and the initial condition: when \(C_e<C_p\), the derivative must be positive, and when \(C_e>C_p\), it must be negative. Second, reproduce the constant-input analytic solution and then a changing-input benchmark. Third, reduce fixed step sizes or tighten adaptive tolerances and compare predicted concentrations, peak times, or other quantities used in the analysis. Material changes indicate that the computation has not yet demonstrated numerical convergence.
| Approach | Useful feature | Important limitation |
|---|---|---|
| Forward Euler | Simple to implement and inspect. | Low order; accuracy and stability can require small steps. |
| Fixed-step RK4 | Good accuracy for smooth inputs at moderate step sizes. | Does not adapt its step or automatically certify total error. |
| Adaptive ODE solver | Adjusts step size using estimated local error. | Tolerances are not a substitute for convergence checks or correct input handling. |
| Closed-form solution | Exact for inputs with an available analytic solution. | Not generally convenient for arbitrary inputs or more complex coupled systems. |
Common warning signs include oscillating or implausibly negative Euler estimates, a trajectory that changes when the output grid changes, and failure to reproduce a constant-input benchmark. A solver can report small estimated error while still solving the wrong problem—for example, if concentrations and time use inconsistent units, if an input knot is skipped, or if the wrong interpolation is supplied. Numerical integration also addresses only the ODE calculation: it does not resolve uncertainty in the plasma input, parameter values, or the adequacy of the effect-compartment model.
Common mistake: A user may conclude that a simulation is accurate because an adaptive solver completed without warning or because its output is smooth. Solver completion and smoothness do not establish accuracy. Compare against an analytic case, inspect the input handling, and verify that predictions are stable under tighter tolerances or smaller steps.
Check Your Understanding
Answer each question, then check the solutions below.
- If \(C_p\) is constant and \(C_e(t_0)=C_p\), what does the exact solution predict, and why is this case important when deriving the formula?
- For Euler with \(k_{e0}=0.8\) h\(^{-1}\) and \(h=1\) h, calculate the multiplier on the deviation from constant \(C_p\). Is the update stable, and does it preserve a non-oscillating approach?
- For \(C_p=5\) mg/L, \(C_e(0)=1\) mg/L, and \(k_{e0}=1\) h\(^{-1}\), calculate the exact effect-site concentration at 1 hour.
- Why should an integration step end at a knot where a piecewise-linear plasma input changes slope?
Solutions. (1) The solution is \(C_e(t)=C_p\) for all \(t\). The integrating-factor result covers this case directly; separation followed by division by \(C_e-C_p\) would not. (2) The multiplier is \(1-hk_{e0}=1-0.8=0.2\). Since its absolute value is less than 1, the update is stable; since it is positive, the deviation does not alternate in sign. (3) \(C_e(1)=5+(1-5)e^{-1}=5-4(0.367879)=3.5285\) mg/L. (4) The input’s slope changes at the knot, so a step spanning it would not represent one smooth linear segment; ending there ensures the solver evaluates the specified input on each segment correctly.