Why Numerical Error Needs Its Own Check
A PBPK prediction generated by Monte Carlo has at least two distinct computational approximations. First, a finite number of parameter or patient draws represents the uncertainty distribution; this creates Monte Carlo sampling error. Second, each draw's differential equations and derived outputs are computed numerically, using a solver, tolerances, and often a time grid. This creates numerical error. A stable-looking prediction interval does not by itself show that both sources are small.
To assess numerical error in a propagated prediction, preserve the Monte Carlo draws and rerun selected draws with tighter solver tolerances or a refined time step. Because the inputs are paired, output differences for the same draw directly diagnose sensitivity to numerical settings. Then examine whether those changes could alter the output quantiles of interest. This is not a replacement for assessing Monte Carlo error: it is a separate verification step.
For a sample of \(n\) outputs, define the empirical \(p\)-quantile as \(\hat Q(p)=\inf\{y:\hat F_n(y)\ge p\}\), where \(\hat F_n\) is the empirical cumulative distribution function. With the nearest-rank convention, this is the \(k\)th sorted output, where \(k=\lceil np\rceil\). State the quantile convention when reporting a numerical comparison; software may use other interpolation conventions.
A Paired Refinement Check
The baseline and refined calculations should differ only in the numerical settings being investigated. Use the same parameter draws, model equations, units, dosing schedule, simulation horizon, event definitions, and output-extraction rules. Save the draw identifiers and both sets of outputs. If a draw fails under either configuration, record and investigate the failure rather than silently dropping it from one calculation.
Select draws deliberately. For a target quantile, include draws whose baseline outputs lie near its rank, with draws on both sides of the target. A small selection far from the quantile may reveal general solver instability but is weak evidence about that quantile. Additional draws can be selected from regions with rapid concentration changes, stiff dynamics, or event-sensitive outputs. Selection based on baseline outputs is useful for diagnosis, but it does not guarantee that unrerun draws have negligible numerical error.
Save the exact input draws, random seeds or draw records, and baseline outputs. Do not resample between the two solver runs.
For selected draws, tighten tolerances and, where relevant, reduce the maximum time step. Keep the scientific model and output definition unchanged.
Calculate \(\Delta_i\), inspect its size and pattern, and replace rerun outputs in the full baseline ensemble before sorting if reporting a hybrid quantile diagnostic.
Compare numerical changes with a prespecified, decision-relevant tolerance and investigate instability; do not treat a tighter setting alone as proof of accuracy.
Why Paired Differences Separate the Sources
Write \(Y_i^\ast\) for the output that would be obtained for draw \(i\) from an ideal numerical solution of the specified model, and let \(Q_p\) be the target \(p\)-quantile of the model's output distribution. Conceptually, the total difference between a computed sample quantile and the target can be decomposed as
The first bracket is numerical discrepancy for the finite set of draws; the second is Monte Carlo sampling error for that set. The ideal solution is not available in practice, so a refined run serves as an approximation to it. Paired differences help because each comparison holds \(\theta_i\) fixed: variation from drawing different inputs is not mixed into \(\Delta_i\). They do not, by themselves, estimate the Monte Carlo error of the quantile.
To see why, write the sorted baseline and refined outputs as \(Y_{(1)},\ldots,Y_{(n)}\) and \(Z_{(1)},\ldots,Z_{(n)}\). The coordinatewise bound implies \(Y_i\le Z_i+\epsilon\) for every draw. Sorting preserves this order bound, so \(Y_{(k)}\le Z_{(k)}+\epsilon\). Reversing the roles gives \(Z_{(k)}\le Y_{(k)}+\epsilon\). Therefore \(|Y_{(k)}-Z_{(k)}|\le\epsilon\), including the rank \(k=\lceil np\rceil\) used for the nearest-rank quantile. This guarantee requires a bound for every draw; checking only a subset cannot establish it for the entire ensemble.
Worked Example: A Fully Specified Hybrid Quantile
Suppose a hypothetical ensemble has \(n=10\) baseline predictions for a concentration output, in mg/L. The sorted values are 0.80, 0.91, 1.02, 1.10, 1.18, 1.26, 1.34, 1.42, 1.55, and 1.80. The target is the nearest-rank 80th percentile.
The rank is \(k=\lceil np\rceil=\lceil10(0.80)\rceil=8\), so the baseline estimate is \(1.42\) mg/L. Suppose the draws at ranks 7, 8, and 9 are rerun, with refined outputs 1.33, 1.46, and 1.54 mg/L, respectively. All other seven outputs are specified to remain at their baseline values for this hybrid diagnostic.
Replacing those three outputs and sorting the full hybrid ensemble gives 0.80, 0.91, 1.02, 1.10, 1.18, 1.26, 1.33, 1.46, 1.54, and 1.80 mg/L. Its rank-8 value is \(1.46\) mg/L. Thus the hybrid rank-8 estimate changes by \(1.46-1.42=0.04\) mg/L. This calculation is valid because the complete ensemble and rerun status are specified. The change is a numerical sensitivity diagnostic, not a Monte Carlo confidence interval or proof that the unrerun draws have zero numerical error.
Worked Checks at the Draw and Quantile Levels
Worked Example: Refining a Draw Near a Prediction Quantile
In a hypothetical 12-draw ensemble, the nearest-rank 75th percentile is at rank \(\lceil12(0.75)\rceil=9\). The complete sorted baseline outputs, in mg/L, are 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, and 2.00. The baseline quantile is therefore \(1.60\) mg/L.
Suppose the outputs at ranks 8, 9, and 10 are rerun at tighter tolerances and become 1.52, 1.67, and 1.68 mg/L. The other nine baseline outputs are known and retained. After replacement, the sorted hybrid outputs are 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.52, 1.67, 1.68, 1.80, and 2.00 mg/L. The hybrid rank-9 estimate is \(1.67\) mg/L, a change of \(1.67-1.60=0.07\) mg/L.
This is evidence that the selected draws can affect the target estimate under refinement. It is not an exact fully refined ensemble quantile: the nine unrerun outputs may also change. The calculation is useful precisely as a bounded-scope diagnostic, provided it is labelled as a hybrid result and not mistaken for a verified final quantile.
Worked Example: Checking a Derived Output Across Time Steps
For one hypothetical PBPK draw, a baseline run gives an AUC of \(42.60\) mg·h/L and a maximum concentration, \(C_{\max}\), of \(5.20\) mg/L. Runs with progressively smaller maximum time steps give AUC values of \(42.40\) and \(42.38\) mg·h/L, and \(C_{\max}\) values of \(5.08\) and \(5.07\) mg/L. Assume the equations, dosing events, time horizon, and output-extraction method are identical across runs.
The change from baseline to the finest run is \(|42.38-42.60|=0.22\) mg·h/L for AUC, or \(100(0.22/42.60)=0.5164\%\). For \(C_{\max}\), it is \(|5.07-5.20|=0.13\) mg/L, or \(100(0.13/5.20)=2.500\%\). AUC changes become smaller from the intermediate to the finest run: \(|42.38-42.40|=0.02\) mg·h/L. The corresponding \(C_{\max}\) difference is \(|5.07-5.08|=0.01\) mg/L.
These results suggest better stabilization of AUC than of \(C_{\max}\) for this draw. They do not prove a formal error bound: solver tolerances are not necessarily bounds on errors in every reported output, and a single draw cannot establish quantile-level stability. If a hypothetical acceptance criterion were a change below 1% for each output, the finest-run AUC change would meet it, while the \(C_{\max}\) change would not.
Diagnostics, Alternatives, and Limitations
Look for output differences that persist as settings are refined, grow near dosing or discontinuity events, or vary systematically across draws. Check that event times, interpolation rules, units, and derived-output calculations are identical. In particular, \(C_{\max}\) can be sensitive to a coarse grid or peak-extraction rule even when a smooth integrated quantity such as AUC appears stable. Solver tolerances control aspects of local numerical error; they should not be presented as guaranteed global bounds on a PBPK output.
A useful comparison is to rerun the entire ensemble at refined settings. This directly permits a fully refined empirical quantile comparison and, if every paired output difference is bounded, the order-statistic result above applies. Its cost may be substantial. Targeted reruns reduce cost and help locate sensitivity near important ranks, but leave numerical error in unrerun draws unresolved. Repeating the Monte Carlo simulation with new random draws addresses sampling variation instead; it does not isolate numerical error because both the inputs and outputs change. These approaches answer different questions and can be combined when resources allow.
No numerical refinement test proves that the implemented equations are scientifically correct. Nor does agreement across two settings rule out shared implementation errors, an unsuitable time horizon, incorrect event handling, or a mistaken output definition. Compare against analytic or independently verified solutions where available, and use several progressively refined settings when feasible. Prespecify an acceptable output change in context: a numerically small difference may still matter if it changes a decision or a reported interval endpoint.
- Do not infer an ensemble quantile from a few rerun values. Their positions among all outputs and the values around the target rank determine the order statistic.
- Do not call paired reruns Monte Carlo replicates. The input draws are held fixed; the comparison diagnoses numerical sensitivity.
- Do not assume tighter tolerances guarantee accuracy. Check output convergence, time-step effects, and numerical events.
- Do not confuse a stable mean with a stable quantile. The output ranks near a quantile may move even when an overall average barely changes.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- Why should a numerical refinement check reuse the same parameter draws?
- For \(n=20\) and a nearest-rank 90th percentile, which order-statistic rank is used?
- If every output in an ensemble changes by no more than 0.03 mg/L, what bound follows for the change in its nearest-rank quantile?
- Why can rerunning only three outputs near a baseline quantile fail to determine the refined ensemble quantile?
- A new Monte Carlo run produces a different quantile. Does that alone demonstrate solver error? Explain.
Solutions. 1. Reusing the same draws pairs outputs, so differences isolate changes associated with numerical settings rather than changes in sampled inputs. 2. The rank is \(\lceil20(0.90)\rceil=18\). 3. The quantile changes by at most 0.03 mg/L, provided the bound applies to every paired output and the same quantile convention is used. 4. The selected outputs may not occupy the relevant ranks after replacement, and unrerun outputs may also move; the full hybrid ensemble must be sorted to calculate its rank. 5. No. New draws change Monte Carlo sampling variation as well as potentially numerical calculations; a paired comparison holding draws fixed is needed to diagnose numerical sensitivity.