Why Test Numerical Convergence?
A PBPK model specifies differential equations, but a computer returns an approximation to their solution. A plausible concentration-time profile is not, by itself, evidence that the numerical calculation is sufficiently accurate. Integration error can affect predictions at individual times, peak concentration, area under the curve, and other summaries. A convergence assessment asks whether these results remain acceptably stable when numerical settings are made more stringent or the integration method is changed.
This is a check on the numerical solution of a specified model. It is not evidence that the model structure or parameter values describe a real biological system correctly. The comparison is meaningful only when the runs use the same equations, parameters, initial conditions, doses, event times, units, and requested outputs. For repeated dosing, the same dosing schedule and the same handling of every dose event must be used in each run.
Let the model state be \(x(t)\), governed by \(dx/dt=f(t,x)\), and let \(C(t)\) denote a concentration derived from the state. A solver approximates the trajectory at requested times and may also integrate extra states for summaries such as AUC. Comparisons should use common times and matching definitions of outputs. A useful pointwise measure is the absolute difference, \(|C(t)-C_{\mathrm{ref}}(t)|\). A relative difference divides that quantity by a reference magnitude, but becomes unstable when the reference concentration is near zero. In that region, use an absolute tolerance or a scaled measure with a stated concentration floor.
What Is Being Refined?
Adaptive solvers estimate local integration error and adjust step sizes to meet their tolerances. In common implementations, the allowable error for state \(i\) is scaled using a quantity of the form \(\mathrm{atol}_i+\mathrm{rtol}|x_i|\), where \(\mathrm{atol}_i\) is an absolute tolerance and \(\mathrm{rtol}\) is a relative tolerance. Exact error-control details differ across software. The key practical point is that tolerances must be suitable for the scale and units of each state. A small absolute tolerance in one state may be negligible for a large state but important for a small one.
Three complementary checks are useful. First, tighten tolerances while keeping the method and other settings fixed. Second, compare a different appropriate solver method, using stringent settings for both. Third, refine the maximum step size while keeping the tolerances fixed. Each probes numerical sensitivity in a different way. None should change the model or the intended event schedule.
Fix the model, parameters, initial state, dose and event schedule, units, and output times. State which concentration and exposure metrics will be compared.
Use an analytic benchmark if available; otherwise use a stringent numerical run and seek confirmation from an independent method or additional refinement.
Tighten tolerances, change solver method, or reduce the maximum step size. Keeping the other settings fixed makes the source of any change easier to interpret.
Check concentrations at shared times and relevant summaries such as AUC and maximum concentration. Compare differences with prespecified accuracy requirements.
A Refinement Argument
For a smooth problem and a fixed numerical method in its asymptotic range, the global error often behaves approximately as \(e(h)=Kh^p\), where \(h\) is a representative step size, \(p\) is the method's convergence order, and \(K\) is a problem-dependent constant. If the step size is halved, this approximation gives \(e(h/2)\approx e(h)/2^p\). Thus, refinement should reduce error when the assumptions hold.
The difference between two approximations can also estimate the remaining error. Write \(y_h=y^*+Kh^p\) and \(y_{h/2}=y^*+K(h/2)^p\), where \(y^*\) is the exact result. Subtracting gives \(y_h-y_{h/2}\approx Kh^p(1-2^{-p})\), while \(y_{h/2}-y^*\approx Kh^p2^{-p}\). Therefore:
This is a justification for using refinement differences as an error indicator, not a universal error bound. The estimate relies on the same method, a smooth solution, a suitable range of step sizes, and approximately asymptotic error behavior. Adaptive step sequences and dose events can complicate that behavior; near discontinuities, steps should respect event times and the simple fixed-step argument may not apply.
Worked Example: Tightening Tolerances Against an Analytic Reference
Consider a hypothetical one-compartment elimination benchmark with concentration \(C(t)=10e^{-0.5t}\) mg/L, where \(t\) is in hours. This analytic profile is the reference. Its AUC from 0 to 8 hours is \(20(1-e^{-4})=19.6337\) mg·h/L. Two numerical runs use the same model and requested times; only the solver tolerances are changed. The reported AUC is obtained by integrating concentration as an additional state.
| Time (h) | Reference (mg/L) | Looser tolerances (mg/L) | Tighter tolerances (mg/L) |
|---|---|---|---|
| 0 | 10.0000 | 10.000 | 10.000 |
| 1 | 6.0653 | 6.100 | 6.066 |
| 2 | 3.6788 | 3.700 | 3.679 |
| 4 | 1.3534 | 1.360 | 1.3535 |
| 8 | 0.1832 | 0.184 | 0.1832 |
At 2 hours, the looser run differs from the reference by \(|3.700-3.6788|=0.0212\) mg/L, or \(0.0212/3.6788\times100=0.576\%\). The tighter run differs by \(0.0002\) mg/L, or \(0.00544\%\). The looser run's AUC is 19.780 mg·h/L, so its relative AUC difference is \(|19.780-19.6337|/19.6337\times100=0.745\%\). The tighter run's AUC is 19.635 mg·h/L, giving a difference of \(0.0013\) mg·h/L, or \(0.00662\%\). In this example, tightening tolerances substantially improves agreement for both a point concentration and the integrated exposure summary.
Comparing Methods and Exposure Summaries
A second solver method is a useful cross-check because different numerical methods can have different error patterns. Agreement is most informative when both methods are suitable for the problem, both are run with appropriately stringent settings, and event handling is identical. The method with the smaller difference from the reference is not necessarily more accurate at every time or for every output. Compare the outputs that matter for the intended analysis rather than relying on a single solver statistic.
For each chosen metric, define an acceptance requirement before inspecting results. For example, a project might require relative differences below 0.1% for AUC and peak concentration, alongside an absolute concentration difference below a specified limit near zero. Such thresholds are application-specific: there is no universal convergence percentage that is suitable for every PBPK use.
Worked Example: Comparing Two Stiff Solver Methods
Suppose a hypothetical repeated-dose PBPK simulation has a carefully refined reference with AUC over one dosing interval of 24.810 mg·h/L, maximum concentration \(C_{\max}=5.000\) mg/L, and concentration at 24 hours of 0.3421 mg/L. Two appropriate stiff methods are run on the same model, with identical events and requested outputs. For this example, the prespecified requirement is a relative difference no greater than 0.1% for each metric.
| Method | AUC (mg·h/L) | \(C_{\max}\) (mg/L) | \(C(24)\) (mg/L) |
|---|---|---|---|
| Method A | 24.790 | 4.998 | 0.3424 |
| Method B | 24.812 | 5.000 | 0.3420 |
For Method A, the relative AUC difference is \(|24.790-24.810|/24.810\times100=0.0806\%\). The relative \(C_{\max}\) difference is \(|4.998-5.000|/5.000\times100=0.0400\%\), and the relative difference at 24 hours is \(|0.3424-0.3421|/0.3421\times100=0.0877\%\). All three meet the stated 0.1% requirement.
For Method B, the corresponding differences are \(0.002/24.810\times100=0.00806\%\), \(0.000/5.000\times100=0.000\%\) at the displayed precision, and \(0.0001/0.3421\times100=0.0292\%\). These also meet the requirement and are smaller at the displayed precision. This supports numerical stability for these metrics under the tested settings; it does not establish that every point on the concentration-time curve is equally accurate.
Refining the Maximum Step and Checking the Profile
Tolerance control alone does not guarantee that every important feature is resolved. A maximum step limit can help ensure that the solver does not take steps too large to represent a rapid change, even if its local error estimate appears acceptable. This can matter around rapid distribution changes, sharp absorption phases, or changes in input. Step refinement should supplement, not replace, suitable tolerance settings.
Compare profiles at common times, including times near expected peaks and event boundaries. If concentrations are evaluated by interpolation between solver steps, ensure that the comparison uses consistent output times and interpolation practices. A useful scaled profile error is \(E_{\max}=\max_j |C(t_j)-C_{\mathrm{ref}}(t_j)|/\max(C_{\mathrm{ref}}(t_j),C_{\mathrm{floor}})\), where \(C_{\mathrm{floor}}\) is a positive, stated concentration scale. The floor prevents tiny reference values from making relative differences arbitrarily large. Report the absolute difference too, so the scale is transparent.
Worked Example: Refining a Maximum Step Size
A hypothetical reference profile and two simulations use the same PBPK equations, dose events, and tolerances. Only the maximum step size differs. Concentrations are compared at shared times, and the chosen floor for the scaled error is \(C_{\mathrm{floor}}=0.05\) mg/L.
| Time (h) | Reference (mg/L) | Maximum step 2 h (mg/L) | Maximum step 0.5 h (mg/L) |
|---|---|---|---|
| 0.5 | 4.80 | 4.72 | 4.799 |
| 1 | 3.20 | 3.11 | 3.198 |
| 2 | 1.75 | 1.68 | 1.749 |
| 4 | 0.62 | 0.60 | 0.620 |
| 8 | 0.08 | 0.079 | 0.080 |
At 2 hours, the 2-hour maximum-step run has an absolute difference of \(|1.68-1.75|=0.07\) mg/L and a scaled difference of \(0.07/1.75=0.0400\), or 4.00%. The 0.5-hour run has an absolute difference of \(0.001\) mg/L and scaled difference \(0.001/1.75=0.000571\), or 0.0571%. Across the displayed times, the largest scaled difference for the 2-hour run is 4.00%, while for the 0.5-hour run it is \(0.002/3.20\times100=0.0625\%\) at 1 hour. The finer maximum step therefore gives much closer agreement on this grid. This conclusion is limited to the sampled profile and the stated reference and settings.
Diagnostics, Limitations, and a Common Pitfall
A robust assessment checks whether refinement changes the results less and whether independent methods agree on relevant outputs. Also inspect basic numerical plausibility: concentrations should not become materially negative without a model-based reason; mass balance or other conservation relationships should be consistent with the model; and event times should be represented correctly. A failure of any such check can reveal implementation or resolution problems even if two reported summaries happen to be close.
A numerical reference can itself be inaccurate. Strengthen it by tightening its settings further, using an independent method, or comparing against an analytic benchmark where possible. Check important intermediate states as well as final summaries: similar AUCs can conceal compensating concentration errors. Likewise, agreement at a sparse set of times does not ensure agreement between them, especially near a narrow peak or a discontinuity. Report the reference construction, solver settings, event handling, comparison times, metrics, and observed differences so the evidence can be evaluated.
Common mistake: “The solver reported success, so the predictions are accurate.” A successful return usually means the solver completed according to its own error-control rules; it does not establish that settings were suitable for every state or output. Instead, compare refined runs and an appropriate reference, examine both concentrations and exposure summaries, and judge the differences against requirements relevant to the application. Conversely, convergence is not model validation: numerical agreement says nothing by itself about whether the PBPK assumptions or parameters are biologically adequate.
Check Your Understanding
Answer each question using the definitions and comparison principles above.
- Why should a solver comparison use the same dose-event schedule and requested output times?
- A reference concentration is 2.00 mg/L and a refined result is 2.02 mg/L. What is the relative difference?
- Why can relative error be misleading when the reference concentration is close to zero?
- Two runs have nearly identical AUCs but visibly different peaks. What does this imply about relying on AUC alone?
- Does agreement between two refined solvers show that a PBPK model is biologically valid?
Solutions. (1) Otherwise, differences can reflect changed inputs, event handling, or evaluation times rather than numerical integration. (2) The difference is \(0.02\) mg/L, so the relative difference is \(0.02/2.00\times100=1.00\%\). (3) Dividing by a very small reference can make a small absolute difference appear arbitrarily large; use an absolute tolerance or a stated positive scale floor in that region. (4) AUC alone can hide pointwise or peak discrepancies, so compare the profile and peak-related outputs too. (5) No. Agreement supports numerical stability for the tested calculations, not biological adequacy or external predictive performance.