Why Solver Settings Belong in a Mechanistic Code Review
A QSP model specifies how biological states change, but a computer produces simulated trajectories by approximating those changes numerically. Two implementations of the same equations can therefore produce different outputs if they use different solvers, tolerances, step-size limits, or event conventions. Some differences may be negligible; others can change a predicted threshold crossing, peak concentration, or treatment response.
Reviewing these choices is not the same as deciding whether a model is biologically valid. The aim is to determine whether the numerical method and event logic are appropriate for the model, sufficiently documented, and capable of resolving the features that matter to the simulation protocol.
A solver review should begin by identifying the model's time scales and irregularities: rapid turnover, slow accumulation, dosing times, switches, and threshold-triggered changes. Then inspect whether the numerical settings can resolve those features. This requires the equations and event definitions to be unambiguous before judging whether a computed trajectory is plausible.
What to Inspect in Solver Configuration
Common choices include explicit or implicit methods, fixed or adaptive step sizes, and tolerances for estimating numerical error. Adaptive methods adjust the internal step size in response to estimated error; fixed-step methods use a specified increment. Explicit methods calculate the next state directly from current information. Implicit methods require solving an equation at each step, but can be more stable for some rapidly changing or stiff systems.
A typical adaptive solver controls estimated local error using absolute and relative tolerances. For a state component \(y_i\), a commonly used scale is $$ \mathrm{atol}_i+\mathrm{rtol}\,|y_i|. $$ Here \(\mathrm{atol}_i\) sets a scale for small or near-zero values, while \(\mathrm{rtol}\) sets a scale relative to the state magnitude. Implementations may use different error norms and details, so the exact solver and configuration should be recorded. Tolerances control an estimate of local numerical error; they do not guarantee that the total error in every reported quantity is below that value. Nor do they address uncertainty in model structure or parameters.
Review absolute tolerances with state units in mind. A single absolute tolerance may be inappropriate when one state is measured in nanomolar concentration and another in cell counts. Also inspect minimum and maximum step sizes, the output-time grid, solver termination status, and whether the calculation reached a maximum-step or convergence limit. A dense output grid does not necessarily mean that the solver used equally small internal steps.
Worked Example: Checking Step-Size Sensitivity for Drug Elimination
Consider a hypothetical concentration \(C(t)\) that follows first-order elimination, \(dC/dt=-kC\), with \(C(0)=2.00\) mg/L and \(k=0.500\) per hour. The analytic solution is \(C(t)=2.00e^{-0.500t}\) mg/L, so at \(t=2.00\) hours it is \(2.00e^{-1}=0.735759\) mg/L. Use forward Euler calculations to see why step-size checks matter.
Forward Euler updates the state by \(C_{j+1}=C_j+h(-kC_j)=C_j(1-kh)\), where \(h\) is the time step in hours. With \(h=0.500\) hours there are four steps to 2.00 hours:
This is \(0.102946\) mg/L below the analytic value. With \(h=0.100\) hours, there are 20 steps:
The second approximation is closer, with an absolute difference of \(0.018787\) mg/L. This comparison illustrates a diagnostic, not a universal error guarantee: for this smooth decay problem, reducing the step size improves these Euler results, but a production solver's error behavior depends on its method and error control. Review whether the chosen settings resolve the model's fastest relevant dynamics and whether important outputs remain stable under a tighter configuration.
Discontinuities and Scheduled Events
A scheduled bolus, a parameter switch, or a change in equations at a specified time makes the simulation piecewise-defined. The solution should be advanced to the event time, the event action applied according to a stated convention, and integration then continued using the post-event state or equations. A solver step should not silently average across a jump as if the dynamics were smooth.
For a bolus that adds dose \(D\) to concentration, a common convention is \(C(t_e^+)=C(t_e^-)+D\), where \(t_e\) is the event time and the superscripts denote the values immediately after and before the event. The event schedule should specify units, whether an event occurs at the initial time, and how simultaneous events are ordered if their order affects the state. Output at the event time also needs a clear convention: it may represent the left limit, right limit, or both.
Worked Example: Reviewing a Scheduled Bolus
Suppose a hypothetical drug concentration starts at \(2.00\) mg/L and follows \(dC/dt=-0.500C\) per hour. At \(t=2.00\) hours, a bolus instantaneously adds \(1.00\) mg/L. Find the concentration at \(t=3.00\) hours under the convention that the bolus is applied at exactly 2.00 hours.
Before the event, \(C(2.00^-)=2.00e^{-0.500(2.00)}=0.735759\) mg/L. The state immediately after the bolus is therefore \(C(2.00^+)=0.735759+1.00=1.735759\) mg/L. Elimination then continues for one hour:
If a numerical implementation applies the bolus late, at 2.20 hours, the concentration at 3.00 hours would instead be \(2.00e^{-1.5}+1.00e^{-0.500(3.00-2.20)}=0.446260+0.670320=1.116580\) mg/L. That is \(0.063789\) mg/L higher than the correctly timed result because the dose has had less time to eliminate. A code review should check that the solver reaches the scheduled time, applies the jump once, and restarts or continues with the correct post-event state.
State-Triggered Events and Event Detection
Some events occur when a state reaches a threshold, rather than at a prespecified time. For example, a treatment rule may switch when a biomarker falls below a limit. The implementation must specify the event function, crossing direction when relevant, and the action taken at the crossing. A solver generally locates a crossing between internal steps by root finding or another event-location procedure. Merely checking whether an output value is above or below a threshold can delay detection until the next output time.
Event logic also needs safeguards against unintended repeated triggering. If an event condition remains true after a reset, it may fire repeatedly unless it is made one-shot, deactivated, or defined with an appropriate reset or hysteresis rule. Whether a crossing is upward, downward, or either direction matters when the state can cross the same threshold more than once.
Worked Example: A Threshold Crossing Between Output Times
Suppose a hypothetical biomarker follows \(B(t)=12e^{-0.200t}\) units, with time in hours. A downstream process starts accumulating at a rate of 1.00 unit per hour when \(B(t)\) first falls to 8.00 units. Find the event time and compare event localization with checking only hourly output values.
Solving the threshold equation gives $$ 12e^{-0.200t_e}=8 \quad\Longrightarrow\quad t_e=\frac{\ln(12/8)}{0.200}=2.02733\ \mathrm{hours}. $$ At hourly output times, \(B(2)=12e^{-0.4}=8.04384\) units and \(B(3)=12e^{-0.6}=6.58574\) units. The crossing lies between those outputs. If the program only checks at output times and switches at the first value below 8.00, it starts the downstream process at 3.00 hours, a delay of \(3.00-\ln(12/8)/0.200\approx0.972674\) hours.
At \(t=3.00\) hours, the accumulated downstream amount should be \(3.00-t_e\approx0.972674\) units if the crossing is localized correctly and accumulation begins immediately. With the delayed output-grid trigger, it is zero at that time. This discrepancy is caused by event detection, not by a different biomarker equation. A review should verify crossing direction, event localization, and whether the event action takes effect at the crossing or at a later step.
Deriving a Practical Review and Diagnosing Failure
The reason to separate smooth integration from event handling follows directly from the model definition. On an interval with no events, the state evolves under one function \(f\). At an event time, a reset map or changed dynamics may apply. Integrating a smooth approximation across the event instead uses one update for an interval in which the governing rule is not constant. In general, that update need not equal the result of integrating to the event, applying the event action, and integrating onward.
where \(R\) is the specified reset map; for a bolus, it adds the dose to the concentration state. This sequence justifies treating event times as boundaries for the numerical integration or using a solver's correctly configured event mechanism.
Useful review diagnostics include checking solver warnings and status, plotting state trajectories around each event at a fine time resolution, comparing event counts and event times, and repeating key simulations with tighter tolerances or smaller maximum steps. Sudden jumps should match specified resets; smooth states should not show unexplained spikes or oscillations. Nonnegativity, conservation relationships, and expected qualitative trends can reveal problems, but agreement with those checks alone does not establish numerical accuracy.
When two runs differ, compare the model equations, initial conditions, parameter values, solver algorithm, tolerances, step-size controls, event definitions, and output conventions separately. A trajectory difference can arise from a coding change, a numerical approximation, or a different interpretation of event timing. Do not assume that tighter tolerances solve every discrepancy: a missed threshold crossing or incorrect event order is a logic error, not simply a tolerance problem.
Choosing Among Approaches and Recognizing Limits
| Choice | Useful when | Review consideration |
|---|---|---|
| Fixed-step integration | A transparent, controlled step grid is useful. | Confirm the step resolves fast dynamics and event times; small steps can be computationally expensive. |
| Adaptive-step integration | Dynamics vary in speed and error control is needed. | Record tolerances and step limits; verify that event handling and output interpolation are appropriate. |
| Explicit methods | The system is not too stiff and computational simplicity is useful. | Fast decay modes can require very small steps for stability. |
| Implicit methods | Stiff dynamics make explicit methods impractical. | Check solver convergence and computational settings; stability does not ensure adequate accuracy. |
| Output-grid threshold checks | Only when the intended event is explicitly defined on that grid. | Otherwise a crossing can be delayed or missed between output times. |
| Localized root-finding events | The event time is defined by a continuous state threshold. | Specify direction and repeat-trigger behavior, and verify the located event time. |
No solver setting makes a model numerically exact in all circumstances. Tolerances do not measure biological plausibility, and agreement across two solvers does not prove that both resolve every important feature. In addition, a highly accurate numerical solution can still be misleading if event timing, reset logic, units, or model equations are wrong. The review should therefore connect numerical settings to the model's time scales and the quantities used to interpret its outputs.
Common Pitfalls
- Treating a tolerance as a bound on total trajectory error. Tolerances usually control estimated local error; they do not guarantee a corresponding global error bound or address model uncertainty. Check sensitivity of important outputs under tighter settings.
- Assuming that a dense output grid means small integration steps. Solvers may interpolate between internal steps to report requested times. Inspect internal step controls and event-location behavior.
- Applying a jump without defining its timing convention. State whether the pre-event or post-event value is reported at the event time, and how simultaneous events are ordered.
- Detecting a threshold only at reporting times. A crossing can occur between outputs. Use event localization when the model defines a continuous-time crossing, or state explicitly that the rule is grid-based.
- Changing solver settings without recording them. Solver name, version, tolerances, and step limits are part of the computational specification needed to interpret a reproduced trajectory.
Check Your Understanding
Answer each question, then compare your reasoning with the brief solutions.
- Why can an absolute tolerance that is suitable for a concentration state be unsuitable for a cell-count state?
- A scheduled dose is applied at 2.00 hours, but the solver takes a step from 1.80 to 2.20 hours and does not process the event at 2.00 hours. What should a review investigate?
- For \(B(t)=10e^{-0.100t}\), at what time does \(B(t)\) first reach 7.00 units?
- Does a stable trajectory after tightening tolerances prove that a model is biologically valid? Explain.
Solutions. (1) Absolute tolerance sets an error scale in the units and magnitude of each state; the same numerical value can be negligible for one state and large for another. (2) Check whether the solver is instructed to stop at and process the scheduled event time, apply the dose once, and continue from the post-event state; a step must not silently smear or delay the jump. (3) Solve \(10e^{-0.100t}=7\), giving \(t=\ln(10/7)/0.100=3.56675\) hours. (4) No. Numerical stability under tighter tolerances supports a conclusion about sensitivity to those numerical settings, but not about biological validity, model structure, or parameter uncertainty.