Why Simulation Times Need to Match the Dosing History
A simulation can solve the intended model and still produce misleading summaries if its output times do not match the events in that model. A dose may cause an instantaneous jump in amount or concentration; an infusion start or stop changes the rate of input; and a time-varying covariate may change the model’s dynamics. If these times are omitted or treated ambiguously, the recorded profile can conceal an important change even when the numerical solver has integrated the equations successfully.
The key distinction is between integration and reporting. A numerical solver may use many adaptive internal steps to approximate a solution, while returning values only at times requested by the user. Conversely, requesting a value at an event time does not by itself ensure that the solver handled the event with the correct pre-event and post-event rules. A reliable workflow specifies both the model’s event handling and the times at which outputs are required.
Constructing an Event-Aware Grid
For one simulated individual, let \(T_{\mathrm{req}}\) be the requested prediction times, \(T_{\mathrm{dose}}\) the dose times, and \(T_{\mathrm{event}}\) the times of other changes that affect the model. These can include infusion starts and stops, changes in a time-varying covariate, or a change in an input or model rule. A basic event-aware grid is the sorted union of these sets and any additional times needed for the desired output resolution.
The union is a useful starting point, not a complete event-handling rule. If two sets contain the same time, the output need not contain duplicate rows just because the time arose from two sources. However, a state-changing event can require two records at the same numerical time: one immediately before the event and one immediately after. In that case, record the event side explicitly, for example as \(t^-\) and \(t^+\), rather than treating the two values as interchangeable.
For a bolus dose, a simple state update is \(A(t^+) = A(t^-) + D\), where \(A\) is amount and \(D\) is the dose. If volume \(V\) is unchanged, concentration jumps by \(D/V\). The pre-dose value is the state just before applying the update; the post-dose value is the state just after it. By contrast, when an infusion starts or stops, the amount is usually continuous at that instant but the derivative of amount changes because the input rate changes. The solver and output process must follow the event rule appropriate to the event type.
An event-aware grid should not be confused with a demand to force every solver step onto every output time. Many solvers can evaluate their interpolated numerical solution at requested times without using those times as internal integration steps. What matters is that event times are communicated to the solver as events or integration boundaries when required by the implementation, and that outputs are requested on the correct side of any discontinuity.
Gather prediction times, doses, and other changes for each individual. Specify what the model does at each event and whether a pre-event value, a post-event value, or both are needed.
Sort the unique numerical times, retain event labels, and preserve separate pre-event and post-event records where the state jumps. Check that every required prediction time is present.
Check that doses and changes were applied at the intended times and in the intended order. Refine output times where profile features need resolution, and separately check numerical solver accuracy.
Worked Example: Combining Prediction Times With Bolus Doses
Consider a hypothetical one-compartment model with volume \(V=25\) L and elimination rate constant \(k=0.15\) h\(^{-1}\). A 500 mg bolus is given at time 0 and again at 8 h. The requested prediction times are 1, 3, 9, and 12 h, and a regular output schedule includes 0, 4, 8, and 12 h. Construct the output times and illustrate the event-side issue.
The bolus increment in concentration is
The sorted numerical times are 0, 1, 3, 4, 8, 9, and 12 h. At 0 and 8 h, however, a single concentration value cannot represent both sides of the bolus. The event-aware records are \(0^+\), 1, 3, 4, \(8^-\), \(8^+\), 9, and 12 h; include \(0^-\) as well if a pre-first-dose output is needed. The post-dose value at time 0 is 20 mg/L, assuming the initial concentration before dosing is zero.
Immediately before the second dose, the first dose has decayed for 8 h:
The second dose adds 20 mg/L, so \(C(8^+) \approx 6.024+20=26.024\) mg/L. One hour later, the concentration is \(C(9)=26.024e^{-0.15}\approx22.399\) mg/L. This calculation depends on applying the second dose after evaluating the pre-dose state. If the output stores only one unlabelled value at 8 h, it is unclear whether that value is about 6.024 mg/L or 26.024 mg/L.
Why Alignment Matters: A Short Justification
Suppose the state follows a smooth differential equation between events, but an event at \(t_e\) changes it according to \(x(t_e^+)=G(x(t_e^-))\). Integrating across \(t_e\) as though the same smooth equation applied throughout can fail to apply \(G\) at all, or can apply it at the wrong point in the numerical sequence. The sound approach is to integrate to \(t_e\), evaluate the pre-event state, apply the event update, and continue from the post-event state. The output record must preserve these distinctions when the state is discontinuous.
This also explains why merely adding \(t_e\) to a requested output list is not always enough. Output evaluation specifies where to record results; it does not necessarily specify the state update or force the integrator to stop and restart. The event must be represented in the model’s event logic, with the solver informed in the way its implementation requires.
Worked Example: Infusion Start and Stop Times
Suppose a hypothetical patient starts a constant infusion at 2 h and stops it at 6 h. The infusion rate is 20 mg/h, volume is 10 L, and elimination is first-order with \(k=0.10\) h\(^{-1}\). Assume the initial amount is zero and no other doses occur. Requested predictions are at 3, 5, and 8 h; a regular output schedule is 0, 4, and 8 h. Construct the grid and calculate selected concentrations.
The sorted grid is 0, 2, 3, 4, 5, 6, and 8 h. The infusion starts and stops at 2 and 6 h, so both times must be included as events. Unlike a bolus, the infusion start does not cause an instantaneous amount jump: the amount is continuous, but its rate of change becomes \(20-0.10A\) mg/h. At the stop time, the amount remains continuous and the rate changes to \(-0.10A\) mg/h.
Before 2 h, \(A=0\). During the infusion, for \(2\leq t\leq6\), the amount is
At 3 h, \(A(3)=200(1-e^{-0.1})\approx19.03\) mg, giving \(C(3)=A(3)/10\approx1.903\) mg/L. At 5 h, \(A(5)=200(1-e^{-0.3})\approx51.84\) mg, so \(C(5)\approx5.184\) mg/L. At 6 h, the infusion ends with \(A(6)=200(1-e^{-0.4})\approx65.94\) mg, or \(C(6)\approx6.594\) mg/L.
After stopping, there is no input and the amount decays from its value at 6 h. Thus \(C(8)=6.594e^{-0.10(2)}\approx5.399\) mg/L. The concentration is continuous at 6 h, but its slope changes there. If the event time were omitted, a solver or post-processing step could incorrectly treat the infusion as continuing or stopping at a nearby time.
Interpolation, Peaks, and Troughs
An output grid is often sparse because only a limited number of prediction times are needed. Interpolating between those outputs can be useful for visualization or approximate evaluation, but interpolation does not recover features that the grid failed to resolve. This is especially consequential near abrupt dose changes and when a peak or trough occurs between requested times.
For a bolus followed by first-order elimination, concentration between doses is \(C(t)=C_0e^{-kt}\), with \(k>0\). Its second derivative is \(C''(t)=k^2C(t)>0\), so the curve is convex. A straight line connecting two points on this curve lies above the curve between them. Therefore, linear interpolation between two decay-profile values overestimates the true concentration at an interior time. At a bolus time, interpolation across the jump is not meaningful at all: the pre-dose and post-dose values belong to different sides of the event.
A reported peak is also grid-dependent. If the dose creates an instantaneous concentration increase, the post-dose value is a candidate peak and should be explicitly represented. For an infusion, the peak may occur at infusion stop or later, depending on the model and regimen; a coarse grid can miss it. Similarly, the trough for repeated bolus dosing is often the value immediately before the next dose. Recording only the post-dose value at that time can hide the trough entirely.
Worked Example: A Sparse Grid Distorts a Trough and an Interpolated Value
Consider a hypothetical bolus profile with concentration 10 mg/L immediately after a dose at 0 h, followed by elimination with \(k=0.20\) h\(^{-1}\). The dosing interval is 12 h. Suppose output is requested at 0, 4, 8, and 12 h, with the 12 h record stored only after the next bolus.
Before the next dose, the true trough is
At 8 h, \(C(8)=10e^{-1.6}\approx2.019\) mg/L. The post-dose value at 12 h is 10 mg/L, not 0.9072 mg/L. From the stored post-dose values at 0, 4, 8, and 12 h, the smallest is approximately 2.019 mg/L. That sampled minimum misses the actual pre-dose trough of 0.9072 mg/L.
Interpolation across 8 h and 12 h is worse: those endpoints lie on opposite sides of the 12 h bolus, so the segment does not represent the decay curve alone. Even within the decay interval, linear interpolation is only approximate. For example, between 8 h and \(12^-\), linear interpolation gives at 10 h
The approximately 0.110 mg/L difference is the consequence of replacing a convex exponential curve with a straight chord. Including \(12^-\) captures the trough in this example; using a sufficiently refined grid or evaluating the model at the desired time avoids relying on interpolation for the interior value.
Diagnostics, Alternatives, and Limitations
A useful check is to compare the expected event schedule with the actual output records for each individual. Verify that every dose and event appears at its intended time, that the event order is correct when times coincide, and that requested prediction times have not been dropped during sorting or deduplication. For state-changing events, check that the pre-event and post-event values satisfy the event rule. For example, a bolus should produce the specified amount increment, while an infusion stop should not create an instantaneous concentration jump in a model where amount is continuous.
Time representation can cause subtle problems. Floating-point values intended to represent the same time may differ slightly, creating duplicate-looking records or a missed match between an event and a requested prediction time. Use a documented time convention and tolerance for matching values, and preserve distinct events that are genuinely close together. Do not merge nearby times automatically if doing so would change dose order or erase a short infusion.
There are several alternatives to a single sparse output grid. A fixed, dense grid is simple and useful for profile plots, but can generate many unnecessary outputs and still fail to represent a discontinuity correctly if event sides are not distinguished. Event-focused output is more economical and ensures that important boundaries are recorded, but it may not resolve smooth peaks between events. Adaptive numerical solvers control integration error according to their algorithms and tolerances; this does not guarantee that a sparse set of reported outputs captures a peak or trough. A practical workflow combines event-focused times with requested predictions, adds local output points where needed, and checks solver accuracy separately.
Refining the output grid and tightening solver tolerances address different sources of error. Grid refinement improves the temporal resolution of reported values and reduces reliance on interpolation. Solver-tolerance checks assess whether the computed numerical solution is stable to integration settings. Neither check substitutes for correct event logic. For a feature that matters clinically or scientifically, compare results under finer local output resolution and stricter solver settings, and confirm that the event-side convention matches the intended definition of the feature.
Common Mistake
Wrong reasoning: “The solver uses adaptive steps, so a sparse output grid cannot miss a peak or trough.” Adaptive integration controls numerical approximation of the solution; it does not ensure that the returned output times include the extremum or the pre-dose side of a discontinuity. Another error is to treat the single value at a bolus time as self-evidently pre-dose or post-dose.
Correction: Define event behavior explicitly, request the appropriate event-side values, and add output points where the profile needs temporal resolution. Treat interpolation as an approximation with limitations, not as a substitute for event-aware simulation output.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A bolus is given at 5 h, and a prediction is requested at 5 h. Why may one unlabelled output at 5 h be insufficient?
- An infusion starts at 1 h and stops at 4 h. Requested predictions are at 2 and 6 h, and regular outputs are at 0 and 6 h. What numerical output times are needed before adding any extra resolution points?
- For first-order elimination with \(k>0\), is the linear interpolation between two points on a decay segment above or below the true concentration curve? State why.
- A trough is defined as the concentration immediately before a scheduled bolus at 12 h. Which event-side value should be recorded?
- Does requesting an output at an event time necessarily force the solver to apply the event correctly? Explain.
Solutions. (1) The concentration may differ immediately before and after dosing; the pre-dose and post-dose values must be distinguished. (2) The union is 0, 1, 2, 4, and 6 h. The 1 h infusion start and 4 h stop must be included as event times; add both event sides only if the state or output convention requires them. (3) It is above the curve on the interior because \(C''(t)=k^2C(t)>0\), so the decay curve is convex and its chord lies above it. (4) Record \(12^-\), the value just before applying the dose. (5) No. The event rule must be represented and communicated to the solver; requesting a report time alone does not necessarily apply the state update or ensure correct event handling.