Why Convergence Matters in Stiff QSP Models
A numerical solution is an approximation to the solution of a model’s differential equations. In QSP simulations, a pathway may contain both a rapidly changing binding or signaling state and a slowly changing biomarker. When those timescales differ greatly, the system is stiff: an explicit solver may have to take very small steps for numerical stability, even when the state of scientific interest changes slowly. An implicit solver can handle many stiff systems more efficiently, but neither solver type makes its result automatically accurate.
A convergence diagnostic tests whether the quantities used for interpretation remain stable as the numerical calculation is made more stringent. The central practical question is not simply whether two runs look similar. It is whether differences in important states and derived predictions are small relative to their scales and to the intended use, such as comparing a prediction with a decision threshold.
A useful comparison holds the model equations, parameters, initial conditions, units, and output times fixed, and changes only the solver’s numerical error controls. Relative tolerance, often written as rtol, controls error relative to state magnitude; absolute tolerance, atol, supplies an error scale when a state is near zero. Solver documentation defines precisely how these tolerances enter its step-size and error-estimation procedure. They are not direct guarantees of a particular global error in every state or prediction.
A Two-Timescale Pathway With an Analytical Reference
Consider a hypothetical pathway with a fast upstream activation state \(A(t)\) and a slower downstream state \(B(t)\). Time is measured in hours, \(A\) is dimensionless, and \(B\) is measured in normalized activity units. The initial states are zero.
Take \(k_f=200\ \text{h}^{-1}\) and \(k_s=0.02\ \text{h}^{-1}\). The characteristic times are \(1/k_f=0.005\) hours and \(1/k_s=50\) hours; their ratio is \(10{,}000\). This separation makes the system a useful stiff-system example. The equations are linear, the inputs and parameters are constant, and there are no discontinuities over the interval considered.
Solving the first equation gives \(A(t)=1-e^{-k_f t}\). Substituting this into the equation for \(B\) gives a linear equation, which can be solved with an integrating factor. For \(k_f\ne k_s\), its solution is:
This explicit solution provides a reference against which a numerical trajectory can be checked. Integrating it from zero to 24 hours gives the area under the downstream-state curve:
The first term is approximately \(2.469585\), and the subtracted term is approximately \(0.0009529\). Thus the area is about \(2.468632\), not a value near 11.81. Keeping the units and the integration limits attached to this calculation is essential: a numerical solver’s output should be compared with the reference for the same state, interval, and units.
Worked Example: Computing an Analytical Reference
For the pathway above, find \(B(1)\) and use the 0-to-24-hour area as a reference for a solver comparison. The stated assumptions apply: constant rate parameters, zero initial states, and smooth dynamics on the interval.
At one hour, the fast exponential \(e^{-200}\) is negligible at the displayed precision, while \(e^{-0.02}=0.980199\). Substitution gives:
The analytical checks for a solver are therefore \(B(1)\approx0.009852\) and \(\mathrm{AUC}_{0,24}\approx2.468632\) normalized activity-hours. A solver’s agreement at the final time alone would not establish that its intermediate trajectory or integrated prediction is accurate; these are separate quantities to compare.
Comparing Tolerances and Predictions
For a practical refinement check, run the same deterministic model at successively tighter tolerances. Compare trajectories at common output times and compare every derived quantity that matters, such as an AUC, peak value, or value at a specified time. One possible state discrepancy for state \(x_j\) at output time \(t_i\) is:
Here \(s_j\) is a stated, positive, meaningful scale for that state. For example, it could be a reference magnitude or a scientifically relevant change. A scale prevents units from being ignored; it also avoids treating a relative difference near zero as meaningful when the state itself is negligible. The chosen scale and acceptable discrepancy should be justified for the intended use. There is no universal convergence cutoff appropriate to every QSP model.
Worked Example: Refining Tolerances for an AUC
Suppose three hypothetical deterministic solver runs produce the following AUC values for the pathway, with the same model, parameters, initial conditions, and 0-to-24-hour integration interval. The displayed solver values are illustrative; the analytical reference is the value derived above.
| Run | Relative tolerance | Computed AUC (normalized activity-hours) |
|---|---|---|
| Loose | 1e-3 | 2.451000 |
| Intermediate | 1e-5 | 2.467900 |
| Tight | 1e-7 | 2.468610 |
Against the analytical reference \(2.468632\), the tight-run absolute discrepancy is:
The relative discrepancy is approximately \(0.000022/2.468632=0.00000891\), or \(0.000891\%\). The difference between the intermediate and tight runs is \(0.000710\) normalized activity-hours, about \(0.0288\%\) of the tight-run value. These results show refinement toward the analytical value. Whether this is sufficient depends on the prediction’s intended precision; it is not established by the tolerance labels alone.
A convergence check should examine trajectories as well as summaries. Two runs can have similar AUCs while differing at clinically important times, because positive and negative trajectory differences can partly cancel in an integral. Conversely, tiny changes in a state may have little consequence for a derived prediction. State-level and output-level comparisons answer different questions and should both be used when relevant.
Target Binding and Decision-Relevant Error
A fast target-binding state can create the same numerical challenge. Consider a hypothetical occupancy-like state \(\theta\) approaching a fixed equilibrium value of 0.8 with a rate of \(200\ \text{h}^{-1}\). A slower response \(Y\) follows it with a rate of \(0.02\ \text{h}^{-1}\). Both states begin at zero, and the response is driven by half the occupancy:
This simplified system is a numerical diagnostic example, not a full receptor-binding model. It assumes a constant effective equilibrium occupancy and linear response; it omits changing free ligand and other binding or turnover processes. Its purpose is to test whether a fast state and a slow output are both resolved adequately.
Worked Example: Checking a Binding-Driven Prediction Near a Threshold
Because the equations are the same cascade as before, scaled by 0.8, the analytical 0-to-24-hour area for \(Y\) is \(0.8\) times the earlier area:
Suppose a hypothetical solver comparison gives AUCs of 1.9810, 1.9747, and 1.97491 at successively tighter tolerances. The tight result differs from the analytical value by approximately \(0.000004\) fraction-hours; the loose result differs by approximately \(0.006094\) fraction-hours. For a hypothetical decision threshold of 1.980 fraction-hours, the loose result is above the threshold, whereas the analytical value and tight result are below it. This example illustrates why a numerically small error in general may still matter when a prediction is close to a decision boundary. The solver outputs are hypothetical, not observations from a study.
A suitable conclusion is that the tight run agrees closely with the analytical prediction and supports the below-threshold classification for this example. One should not infer that the model is biologically correct, that other outputs have converged, or that all parameter settings will behave similarly.
Diagnostics, Alternatives, and Limitations
Evidence of trouble includes trajectories that change materially under tolerance refinement, oscillations or implausible negative values in states expected to remain nonnegative, abrupt jumps in a smooth model, and derived predictions that change enough to alter a decision. Check which states and time regions contribute to a discrepancy. A fast transient may be missed if output times are too sparse, so use common, sufficiently informative evaluation times; plotting only a coarse grid can hide trajectory differences.
A solver’s reported success flag means that its own numerical procedure met its internal criteria, not that every global error is acceptably small for a particular application. Tolerance scales also matter: a single absolute tolerance may be inappropriate when states have very different units or magnitudes. Set state-specific absolute tolerances where supported, inspect near-zero states with absolute differences, and verify that numerical outputs respect known bounds and other model checks.
Tolerance refinement is one approach among several. An analytical reference is strongest when available for the same equations and conditions. Step-size refinement can be useful with a fixed-step method, but changing step size may not diagnose every error source. Comparing a suitable implicit stiff solver with another solver can provide additional evidence, although agreement between algorithms is not proof of correctness. Regression tests detect changes relative to saved outputs; they do not by themselves establish that either output is accurate. Each approach addresses a different risk.
These diagnostics assess numerical error for a specified implementation and scenario. They do not assess uncertainty in parameters, biological structural discrepancy, or agreement with clinical observations. Those are different questions from whether the chosen equations have been solved accurately. For models with dosing events or other discontinuous inputs, comparisons must also use consistent event handling and evaluation conventions; smooth-interval tolerance refinement alone cannot verify an event implementation.
A Practical Convergence Workflow
Hold the equations, parameters, initial conditions, units, time interval, and output times constant while changing numerical error controls.
Run at least two, preferably several, progressively tighter tolerance settings or step sizes. Record solver settings and warnings.
Calculate discrepancies at common times for important states, then separately compare AUCs, peaks, or other decision-relevant predictions.
State the accepted numerical discrepancy and check whether refinement could change the interpretation or decision. Report limitations if no adequate reference or precision criterion is available.
Check Your Understanding
Answer each question, then check the concise solutions below.
- For the two-timescale pathway, what are the fast and slow characteristic times, and why does their ratio suggest stiffness?
- A loose-tolerance run gives an AUC of 2.4510 and a tight run gives 2.46861 normalized activity-hours. What is their absolute difference?
- Why is comparing only the final state insufficient when an AUC is a prediction of interest?
- A solver reports successful completion, but a tighter run changes a threshold classification. What should be concluded?
Solutions. (1) The times are \(1/200=0.005\) hours and \(1/0.02=50\) hours; the 10,000-fold separation creates fast and slow scales that can constrain numerical integration. (2) The absolute difference is \(|2.46861-2.4510|=0.01761\) normalized activity-hours. (3) The final state does not determine the entire trajectory or its integral; trajectory differences can accumulate or cancel in an AUC. (4) The prediction is not yet numerically stable for that decision. Refine and investigate the calculation, and report the sensitivity; successful completion alone does not resolve the discrepancy.