Why an Integration Fails
A PBPK integration can stop with an error, complete while producing implausible values, or appear successful while mishandling a dosing event. These outcomes do not have one common cause. A solver warning is evidence about what happened during computation, not by itself a diagnosis of the model or the numerical method. Troubleshooting requires connecting three sources of information: the error message, the internal step-size history, and the trajectories of the states.
For a continuous portion of a PBPK model, write the state equations as \(dy/dt=f(t,y)\), where \(y\) contains amounts or other model states. A numerical solver approximates their solution over time, subject to specified error tolerances. Earlier material on solver tolerances, stiffness, and discrete dosing events provides the relevant background; here the emphasis is on identifying why an integration misbehaves and choosing a targeted check.
The most useful first question is not simply “Which solver should I try?” but “What pattern would each possible cause produce?” A rapid sequence of shrinking steps near a fast transient suggests a different problem from a jump in a concentration at an infusion stop, or from a state becoming negative while the solver continues with otherwise ordinary step sizes.
Read Messages, Steps, and Trajectories Together
Messages such as “minimum step size reached,” “maximum number of steps exceeded,” or “Newton iteration failed” describe a computational difficulty. A minimum-step warning means the solver kept reducing its step but still could not satisfy its requirements. Exceeding a maximum-step limit may mean that the requested interval is long or that many small steps were necessary; raising the limit can allow the run to continue, but does not explain why those steps were needed. Newton failures in an implicit solver may occur in a stiff region, but can also result from invalid or poorly scaled equations.
Inspect the step history around the first problematic time. Repeated step reduction in a smooth region may indicate stiffness, excessive tolerances, or a difficult nonlinear solve. A sudden change in step size at a known infusion start or stop may be appropriate: the derivative changes there. If small steps occur only after a state becomes negative or non-finite, first investigate the state and the equations that use it. A trajectory plot should include every state relevant to the failure, not only the plasma concentration.
Record the first warning, non-finite value, implausible state, or event time. Later symptoms may be consequences of an earlier problem.
Determine whether steps shrink gradually, collapse near a scheduled input change, or remain ordinary while a state becomes invalid.
Verify units, parameter signs, initial conditions, derivative calculations, and the intended input value on each side of an event.
Compare a suitable solver or tighter tolerances, or segment the integration at a known event. Confirm that the change resolves the suspected cause rather than merely hiding the message.
Numerical Instability and Step Size
A simple linear equation shows why a step can be too large for a rapidly changing state. For \(dy/dt=-ky\), with \(k>0\), the exact solution is \(y(t)=y(0)e^{-kt}\). Forward Euler, a simple explicit method, updates the state by multiplying it by \(1-kh\), where \(h\) is the step size. Its amplification factor must have magnitude less than one for successive errors in this test equation to decay.
Thus the Euler stability condition is \(0<kh<2\). If \(kh>2\), the computed solution grows in magnitude even though the true solution decays. This condition is a diagnostic illustration, not a universal step-size rule for a full PBPK system: coupled equations, solver methods, and tolerances affect actual behavior.
Worked Example: A Fast State Produces an Unstable Explicit Step
Hypothetical model. Suppose a rapidly exchanging PBPK state follows \(dy/dt=-1000y\), with \(y(0)=1.00\) amount unit and \(k=1000\ \mathrm{h}^{-1}\). Consider an explicit Euler step of \(h=0.00300\) h. The model and calculation are chosen to illustrate a numerical diagnostic, not to represent a particular tissue.
Here \(kh=(1000\ \mathrm{h}^{-1})(0.00300\ \mathrm{h})=3.00\), so the Euler multiplier is \(1-kh=-2.00\). The first updates are \(y_1=-2.00\), \(y_2=4.00\), and \(y_3=-8.00\) amount units. The numerical magnitude grows, alternates in sign, and is inconsistent with the exact solution.
After one step, the exact value is \(y(0.00300)=1.00e^{-3.00}=0.04979\) amount units. The discrepancy and oscillating trajectory point to a step that is too large for this explicit method on the fast timescale. A diagnostic response is to inspect whether the solver is taking or attempting steps of this scale and, if necessary, use smaller steps or a suitable implicit stiff solver. Merely increasing a step limit would not address the stability problem.
The calculation also shows why a negative state needs context. With \(kh=1.20\), Euler's multiplier is \(-0.20\): the magnitude shrinks, but the computed state changes sign. A solver may therefore produce a nonphysical negative amount before it shows outright growth. For an amount known to remain nonnegative under the model equations, either pattern warrants scrutiny.
Discontinuous Inputs and Scheduled Events
An infusion start or stop changes the continuous input and can make \(f(t,y)\) discontinuous at a known time. If a solver takes a step across that time as though the input had one value throughout, its numerical approximation may be inaccurate. The appropriate response is generally to end a continuous integration segment at the event, update the input, and continue from the event time. A bolus is different: it is represented as a discrete state jump at its scheduled time, followed by integration from the updated state. It is not, by itself, a change in the definition of the continuous derivative function.
Worked Example: Checking an Infusion-Stop Event
Hypothetical model. A one-compartment amount \(A\), in mg, starts at zero. An infusion enters at \(R=12.0\ \mathrm{mg/h}\) from \(t=0\) to \(t=1.00\) h, and elimination follows \(k=0.500\ \mathrm{h}^{-1}\). During infusion, \(dA/dt=R-kA\); after the infusion stops, \(dA/dt=-kA\). The infusion is assumed constant over its stated interval.
While infusion is active, the solution starting from \(A(0)=0\) is \(A(t)=(R/k)(1-e^{-kt})\). At the stop time, \(A(1.00)=\{12.0\ \mathrm{mg/h}/0.500\ \mathrm{h}^{-1}\}(1-e^{-(0.500\ \mathrm{h}^{-1})(1.00\ \mathrm{h})})=24.0(1-e^{-0.500})=9.443\ \mathrm{mg}\).
After the stop, the amount decays without input. At \(t=1.20\) h, \(A(1.20)=9.443e^{-(0.500)(0.200)}=8.545\ \mathrm{mg}\), with the exponential argument dimensionless. If the input were incorrectly treated as continuing through \(1.20\) h, the calculation would instead be \(24.0(1-e^{-(0.500)(1.20)})=24.0(1-e^{-0.600})=10.83\ \mathrm{mg}\). The discrepancy is \(10.83-8.545=2.284\ \mathrm{mg}\), using unrounded values for the subtraction.
A trajectory that follows the latter value suggests the stop event was not applied when intended, rather than ordinary instability. Integrating to \(1.00\) h, switching off the input, and then continuing provides a direct check. If the scheduled event is a bolus instead, apply the specified jump to the relevant state at that time and resume integration from the post-jump state.
Nonphysical States and Model-Related Problems
A negative amount or concentration is often a warning, but it does not identify its own cause. Possible causes include a step that is too large, an error in the derivative function, an incorrect initial condition, a sign or unit error in a flux, or a model equation that does not enforce the expected physical constraint. Inspect whether the negative value is a brief numerical excursion, whether it first appears at an event, and whether it persists when the same equations are integrated with a more suitable method and tighter tolerances.
Mass-balance checks can help localize errors. For example, if two compartments exchange material but have no input or elimination, their exchange terms should cancel in the derivative of their total amount. If \(dA_1/dt=-J\) and \(dA_2/dt=J\), then \(d(A_1+A_2)/dt=0\). A nonzero computed change may indicate a coding or implementation problem, though numerical integration error must also be considered. Do not repair a negative result by silently truncating the state to zero: that changes the equations and can conceal the cause.
Worked Example: Distinguishing a Numerical Undershoot from a Persistent Invalid State
Hypothetical model. An isolated amount obeys \(dA/dt=-kA\), with \(A(0)=0.0200\) mg and \(k=4.00\ \mathrm{h}^{-1}\). For a single forward Euler step of \(h=0.300\) h, the update is \(A_1=A_0(1-kh)=0.0200(1-4.00\times0.300)=-0.00400\) mg.
The exact solution at the same time is \(A(0.300)=0.0200e^{-(4.00)(0.300)}=0.006024\) mg, which is positive. In this case \(kh=1.20\), so the Euler step is within the magnitude stability limit \(kh<2\), but its negative multiplier shows that magnitude stability alone does not guarantee a physically acceptable state. The single-step calculation therefore points to an overly coarse step for preserving positivity in this example.
A useful diagnostic is to repeat the integration with smaller steps or a solver that better handles the system, then compare the trajectory and the amount at the same time. If the negative value disappears and the result approaches the exact positive value, numerical undershoot is supported. If negative values persist across appropriate methods and tolerances, check the implementation, parameter units, flux signs, and initial states instead; numerical refinement cannot make an incorrectly specified model physically valid.
Choose the Remedy That Matches the Evidence
| Observed pattern | Likely issue to investigate | Targeted check |
|---|---|---|
| Oscillation or growth during a fast transient; steps repeatedly shrink | Step-size instability or stiffness | Compare a suitable implicit solver, inspect the fast timescale, and check tolerance scaling |
| Mismatch begins at an infusion start or stop | Continuous input changed inside a step or at the wrong time | Segment at the event and verify input values before and after it |
| Negative or non-finite state appears first, then solver fails | Numerical excursion or invalid model calculation | Trace the first affected state and check equations, units, fluxes, and step sizes |
| Many small steps occur without a clear state anomaly | Stiffness, strict tolerances, or difficult nonlinear equations | Review time scales, tolerance choices, Jacobian information, and solver suitability |
The main alternatives are not interchangeable fixes. A stiff implicit solver can help with stability restrictions from fast dynamics, but it cannot correct a missed infusion stop or a wrong elimination sign. Smaller steps or event segmentation address resolution around rapid changes and known input events, but may be inefficient over a long simulation. Tightening tolerances can test sensitivity to numerical error, yet may increase computation and will not fix incorrect equations. An explicit solver may be efficient for nonstiff systems, while an implicit method is often preferable when stability limits dominate; neither method removes the need to check the trajectory and event logic.
A convergence comparison is informative when it is interpreted carefully. If a quantity of interest changes substantially as tolerances are tightened or the solver is changed, its numerical reliability is not established. Similar results across settings provide evidence of numerical robustness for that calculation, not proof that the model structure, parameters, or biological assumptions are correct. The earlier material on testing solver convergence and numerical error addresses that comparison in more detail.
Check Your Understanding
Use the diagnostics above to identify what each pattern supports and what should be checked next.
- For \(dy/dt=-50y\), an Euler step uses \(h=0.050\) h. Calculate the amplification factor and decide whether its magnitude is below one.
- A concentration profile agrees with expectation until the exact infusion-stop time, then remains too high. What implementation check is most direct?
- A solver reports “maximum number of steps exceeded.” Does this message alone establish that the system is stiff? Explain.
- A two-compartment exchange-only model has \(dA_1/dt=-J\) and \(dA_2/dt=J\). What should happen to \(A_1+A_2\), and what does a substantial change suggest?
- A state becomes slightly negative with one coarse explicit integration but remains positive with refined steps and approaches a stable trajectory. Which explanation is supported, and what should not be done automatically?
Solutions. (1) The factor is \(1-kh=1-(50\ \mathrm{h}^{-1})(0.050\ \mathrm{h})=-1.50\); its magnitude is \(1.50>1\), so this Euler step is unstable for the linear test equation. (2) Check that the integration is segmented at the stop and that the infusion input is set to zero afterward; a bolus, if present, is handled as a state jump instead. (3) No. It indicates that the configured step limit was reached; inspect the step history, states, events, tolerances, and model calculations before attributing the issue to stiffness. (4) The total is conserved because \(d(A_1+A_2)/dt=-J+J=0\). A substantial change suggests an implementation or numerical problem to investigate. (5) The pattern supports a numerical undershoot from coarse stepping. Do not automatically clip the state to zero, because clipping alters the model and can conceal the source of the error.