Why Discontinuities Need Their Own Checks
A pharmacometric model may evolve smoothly between events and change abruptly at a dose, reset, or input switch. These are different operations: the differential equation describes the evolution between events, while an event rule changes a state or the model input at a specified time. Treating both as ordinary smooth evolution can give inaccurate predictions even when a solver reports successful integration.
This tutorial focuses on diagnosing failures around those boundaries. The essential question is not only whether the ODE solver converged, but also whether it propagated the state to the correct event time, applied the correct update, and then continued from the correct post-event state. Because likelihoods are calculated from model predictions, an event-handling error can also appear as a sudden likelihood change.
Piecewise Integration and Event Updates
Between events, a state vector \(x(t)\) follows an ODE such as \(dx/dt=f(x,t,\theta,u(t))\). Here \(u(t)\) may represent an input such as an infusion rate. If the input switches at a known event time, the ODE is generally different on the two sides of that time. A bolus dose instead leaves the between-event ODE unchanged but applies a jump map to the state.
For a one-compartment amount \(A(t)\) with first-order elimination, the between-event equation is \(dA/dt=-kA\), where \(k\) is an elimination rate constant. An IV bolus of amount \(D\) at time \(t_d\) is represented by \(A(t_d^+)=A(t_d^-)+D\). The solver integrates to \(t_d\), applies this update once, and continues from the post-dose state. It should not approximate the bolus by integrating through a narrow, artificial pulse unless that approximation has been deliberately specified and validated.
The same principle applies to resets and switches. For an infusion that ends at \(t_s\), the infusion input is active before \(t_s\) and inactive afterward. The state itself is continuous at the stop time, but its derivative may change immediately. For a reset, the state may jump while the differential equation on either side remains the same.
This procedure assumes the between-event ODE has a well-defined solution over each interval—for example, its right-hand side is sufficiently regular in the state for the initial-value problem to have a unique local solution. It also assumes that event times, event types, and update rules are specified consistently. A solver cannot infer whether a recorded dose time is wrong or whether a reset should occur before a dose.
Worked Example: A Bolus Dose at the Wrong Time
Worked Example: A Bolus Dose at the Wrong Time
Hypothetical model. An amount \(A\), measured in mg, follows \(dA/dt=-0.100A\) per hour between boluses. A 100 mg bolus is given at time 0, and a 20 mg bolus is scheduled at time 4 h. Find the amount at time 6 h. Then assess the effect of mistakenly applying the second bolus at 4.5 h.
With the scheduled event, the first dose has been eliminated for 6 h and the second for 2 h. The solution to \(dA/dt=-kA\) over an interval of length \(\Delta t\) is \(A_{\mathrm{end}}=A_{\mathrm{start}}e^{-k\Delta t}\). Thus:
If the second dose is applied at 4.5 h, it has only 1.5 h to be eliminated by time 6 h:
The predictions differ by \(20(e^{-0.15}-e^{-0.2})=0.8395\) mg, calculated from the unrounded expressions before rounding the result. The delayed event raises the predicted amount because the second dose has less time to be eliminated.
To see the likelihood consequence, suppose a hypothetical observation at 6 h is 71.0 mg with a known residual standard deviation of 1.00 mg, and use the Gaussian negative log-likelihood contribution after omitting its common constant, \(q=\frac{1}{2}[(y-\hat y)/\sigma]^2\). For the scheduled event, \(y-\hat y=71.0-71.2557787=-0.2557787\) mg, so \(q=\frac{1}{2}(-0.2557787/1.00)^2=0.03271\). For the delayed event, the unrounded prediction is 72.0953231 mg, giving \(y-\hat y=-1.0953231\) mg and \(q=\frac{1}{2}(-1.0953231/1.00)^2=0.5999\). The larger contribution is a consequence of the changed prediction, not evidence by itself that the ODE solver failed.
Why Crossing an Event Can Fool a Solver
Adaptive ODE solvers estimate local integration error and adjust their step sizes. Their error estimates are designed for the differential equation being integrated; they do not automatically make a discontinuous right-hand side smooth. If a step spans an infusion stop, for example, the solver may use information from both sides of the switch as though one smooth local approximation described the whole interval. An error estimate can then be unreliable, and the numerical state can depend on step placement or interpolation across the boundary.
A tight tolerance is not a substitute for event handling. It controls estimated numerical error within the integration procedure, but it does not correct a wrongly specified event time, a missing dose, an incorrectly ordered pair of updates, or a step that is allowed to treat a known switch as smooth. The practical check is to ensure the solver stops at known boundaries and restarts with the correct state and input.
The justification for splitting is direct. On an interval with no event, the state is obtained by solving the relevant ODE from its initial state. At the event, the model definition supplies the boundary condition \(x(t^+)=J(x(t^-),\theta,e)\), or specifies a new input for the next interval. Applying those steps in order reproduces the model as defined. A single smooth propagation across the event instead solves a different problem unless the state and equation happen to be unaffected by that event.
Worked Example: An Infusion Stop Changes the Derivative
Worked Example: An Infusion Stop Changes the Derivative
Hypothetical model. An amount \(A\) starts at 0 mg. During an infusion at 10 mg/h, \(dA/dt=10-0.5A\) per hour; after the infusion stops at 2 h, \(dA/dt=-0.5A\) per hour. Calculate \(A(3)\), and quantify the error if the solver effectively leaves the infusion on until 2.1 h.
For a constant input \(R\), the solution over an interval of length \(t\) is \(A(t)=A(0)e^{-kt}+(R/k)(1-e^{-kt})\). At the correct stop time, \(A(2)=20(1-e^{-1})\) mg. The amount then decays for one hour:
If input continues for an extra 0.1 h, the extra amount present at 2.1 h relative to the correctly stopped infusion is the contribution of input from 2 to 2.1 h. Each increment of input is subsequently eliminated, so at 3 h the prediction error is:
Therefore the late-stop prediction is approximately \(7.6680+0.62195=8.28996\) mg, using unrounded values for the sum. The amount remains continuous at the stop time, but its derivative changes from \(10-0.5A\) to \(-0.5A\). A diagnostic should check both that the state is continuous there and that the infusion input used after 2 h is zero.
State Updates, Event Order, and Diagnostics
Event handling is easiest to audit when the model’s event schedule is explicit. For each event, inspect its time, type, amount or reset value, affected state, and ordering relative to other events at the same time. Then compare the state immediately before and after the update with the specified jump rule. For an additive bolus, the difference must be \(D\); for a reset to zero, the post-event state must be zero. These checks are more informative than a solver success flag alone.
Order matters when two events occur at the same time and their updates do not commute. For example, applying a dose and then resetting an amount to zero gives a different result from resetting first and then applying the dose. The model or software convention must define the intended order; do not assume that sorting events by time resolves ties correctly.
Worked Example: Dose and Reset at the Same Time
Hypothetical model. Immediately before time 5 h, an amount is 12 mg. At 5 h, a reset sets the amount to zero and a 30 mg bolus is due. The intended order is reset first, then bolus. The amount then follows \(dA/dt=-0.2A\) per hour. Find the amount at 6 h and compare it with the prediction from the reversed order.
With the specified order, the reset gives \(0\) mg, and the bolus gives \(A(5^+)=0+30=30\) mg. After one hour:
If the dose is applied first, the amount becomes 42 mg; the subsequent reset then makes it zero. The 6 h prediction is therefore 0 mg. With a hypothetical observation of 24.0 mg and residual standard deviation 2.00 mg, the Gaussian negative log-likelihood contribution, omitting its common constant, is \(q=\frac{1}{2}[(y-\hat y)/\sigma]^2\). The specified order gives \(q=\frac{1}{2}[(24.0-24.5619)/2.00]^2=0.03947\); the reversed order gives \(q=\frac{1}{2}[(24.0-0)/2.00]^2=72.00\). This large difference is explained by the event order, not by a small integration-tolerance change.
A useful diagnostic sequence is to localize the first interval where the predicted path becomes inconsistent with the model. Compare solver output immediately before and after each event; verify event times and units; check that every event is applied exactly once; and rerun with explicit integration boundaries at all known discontinuities. Then tighten tolerances or refine the integration settings and compare state predictions and the likelihood. If results change materially after enforcing event boundaries, the original calculation may have crossed a discontinuity improperly. If predictions remain stable but the state update is wrong, investigate the event specification rather than the integrator.
Record every dose, reset, and input switch with its exact time, units, and any same-time ordering rule.
Inspect the pre-event and post-event state and verify the specified jump map or input change directly.
Force propagation to stop at known events, apply each update, and restart from the updated state.
Compare predictions and likelihood contributions under tighter settings, while keeping the event schedule and model unchanged.
Alternatives, Failure Modes, and Limits
Event-aware ODE integration is the appropriate default when event times and update rules are known. An analytical solution on each interval can provide a useful independent check for simple models, as in the examples above. A general-purpose numerical solver is still needed for many models, but its integration intervals should respect known discontinuities.
Using very small steps without explicitly aligning them to an event may reduce some numerical error, but it does not guarantee correct event timing or update order. Smoothing a bolus into a short infusion changes the model and can change predictions; it is appropriate only if that smoothed input represents the intended model. Likewise, interpolating a state across a jump can produce values that are not states the model ever occupies.
The diagnosis has limits. A prediction or likelihood that changes with solver settings may indicate inadequate numerical resolution, but it can also expose a coding or event-schedule inconsistency. Conversely, stable output does not prove that the model’s event times or reset rules are scientifically correct. Numerical checks establish that the specified model is being evaluated consistently; they cannot validate the scientific appropriateness of that specification.
A common mistake is to reason, “The solver returned success and the likelihood is finite, so the event handling must be correct.” Solver status concerns the numerical procedure, not whether the right event was applied at the right time in the right order. The correction is to audit the event schedule and state updates explicitly, then test numerical stability with event-aligned integration. Another mistake is to infer solver failure from a large likelihood change alone. As the worked dose example shows, a timing error can substantially change a likelihood even when each piece of the ODE is solved accurately.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- An amount follows \(dA/dt=-0.1A\) per hour and is 50 mg immediately before a 10 mg bolus. What is the state immediately after the bolus, and what is the state one hour later?
- Why can an adaptive solver’s tolerance fail to protect a calculation whose integration step crosses a known infusion stop?
- A reset to zero and a dose of 5 mg occur at the same time. What additional model information is needed to determine the post-event state?
- What should be checked if predictions change substantially when integration settings are tightened around a dose?
Solutions. (1) The additive bolus gives \(A(t^+)=50+10=60\) mg. After one hour, \(A=60e^{-0.1}=54.29\) mg. (2) The local error estimate is based on a smooth propagation within a step; crossing the input switch violates that assumption, so the estimated error may not represent the error from using the wrong input on part of the step. (3) The event-ordering rule is needed: reset then dose gives 5 mg, whereas dose then reset gives 0 mg. (4) Verify the event time and units, ensure integration stops at the event, confirm the update occurs once and in the intended order, and compare pre- and post-event states before attributing the change to ordinary integration error.