Why Repeated Doses Need an Explicit Schedule
A repeated-dose PBPK simulation is not just one long integration of a smooth differential equation. A bolus changes one or more state variables instantly, and an infusion start or stop changes the input to the system abruptly. The solver must not carry a pre-event derivative or an unapplied dose across such a boundary as though nothing changed. A reliable implementation therefore treats event times as boundaries between integration segments and records exactly what happens at each boundary.
The schedule is the union of all relevant bolus times and infusion start and stop times. It is sorted chronologically, with events at the same time grouped for processing. Between consecutive event times, the active inputs are constant or otherwise specified smoothly, and the ODE is integrated on that segment. At an event time, the code applies the scheduled state jumps and updates the active infusion inputs before integrating onward.
Constructing the Event Schedule
For each prescribed bolus, record a unique event identifier, its time, the dose, and the state component or components receiving it. For each infusion, record its start time, stop time, and input rate. The event table is then formed from the union of the bolus times and all infusion start and stop times. Sort the distinct times, but retain all event records at any shared time; combining records must not discard their identities.
A useful convention is to integrate from one event time up to, but not through, the next event time. At the boundary, first retain the solver's left-hand state, then process every event scheduled there, update the active infusion rates, and begin the next integration segment using the right-hand input. Boluses at the initial simulation time are processed before integration begins. Events at the final simulation time are still processed if the requested simulation includes them.
Worked Example: Sorting a Mixed Repeated-Dose Schedule
Consider a hypothetical simulation from 0 to 18 h. A 100 mg bolus is scheduled at 0 h, followed by 100 mg boluses every 8 h. An infusion runs from 3 to 5 h and again from 11 to 13 h. Construct the event times and verify the event inventory.
The bolus times are 0, 8, and 16 h. The infusion contributes start and stop times 3, 5, 11, and 13 h. Sorting their union gives 0, 3, 5, 8, 11, 13, and 16 h. The segments to integrate are therefore [0,3), [3,5), [5,8), [8,11), [11,13), [13,16), and [16,18].
| Time (h) | Event | Bolus records processed | Infusion status after event |
|---|---|---|---|
| 0 | Bolus | 1 | Off |
| 3 | Infusion start | 0 | On |
| 5 | Infusion stop | 0 | Off |
| 8 | Bolus | 1 | Off |
| 11 | Infusion start | 0 | On |
| 13 | Infusion stop | 0 | Off |
| 16 | Bolus | 1 | Off |
The audit totals are three bolus records, two starts, and two stops. Each start has a matching stop, and each bolus identifier appears once. This schedule specifies when integration must be interrupted; it does not require a particular numerical solver or imply that the state changes at an infusion boundary.
Integrating and Processing Each Boundary
On a segment with no events in its interior, integrate the ODE using the input active on that segment. At its right endpoint \(t_j\), save the left-hand state \(x(t_j^-)\). Apply the bolus jump map or maps to obtain \(x(t_j^+)\), and set the infusion input to its value immediately after \(t_j\). Restart the solver from \(x(t_j^+)\) with that right-hand input. When the boundary is an infusion start or stop only, the restart state equals the left-hand state.
For a bolus that adds dose \(D\) to an amount state with bioavailability \(F\), the jump is \(A(t_j^+)=A(t_j^-)+FD\). In a PBPK state vector, the corresponding amount may be added to a particular compartment or distributed according to the model's specified dose-deposition rule. For simultaneous boluses with additive effects, apply each jump once or apply their sum once; these give the same resulting state. Infusion starts and stops update the input rate, not the stored amount.
Assign each bolus, infusion start, and infusion stop a unique identifier, time, and required dose or rate information.
Sort event times and group records at exactly the same scheduled time, retaining every identifier.
Use the currently active inputs and stop the segment at the next event time rather than stepping through it with the old input.
Apply each scheduled jump, update infusion rates, and mark each event identifier as processed.
Restart from the right-hand state and input; at the end, verify that every required event was processed once and only once.
The reason for segmenting can be stated through the ODE solution map. Let \(\Phi_{b,a}\) denote evolution from time \(a\) to \(b\) under the inputs valid between those times, and let \(J_j\) denote the state jump at event \(t_j\). For a bolus followed by a later event \(t_{j+1}\), the correct update is
If the event is an infusion boundary, \(J_j\) is the identity on the state, but the map \(\Phi\) after \(t_j\) must use the new infusion input. Thus the schedule identifies where either the state map or the differential equation changes. Restarting ensures the solver's subsequent evaluations use the correct post-event state and right-hand side.
Worked Amount Calculation Across Infusion Boundaries
Worked Example: Repeated Infusions and a Bolus at a Shared Time
For a hypothetical one-compartment amount model, let elimination be first order with \(k=0.10\ \mathrm{h}^{-1}\). An infusion at \(25\ \mathrm{mg/h}\) runs from 0 to 4 h, is off from 4 to 8 h, and runs again from 8 to 12 h. A 50 mg bolus is scheduled at 12 h, when the second infusion stops. Assume the initial amount is 0 mg and the bolus enters this amount state directly. The event times are 0, 4, 8, and 12 h.
For a constant infusion rate \(R\) over a duration \(\Delta t\), the amount at the end of the interval is \(A_{\mathrm{start}}e^{-k\Delta t}+(R/k)(1-e^{-k\Delta t})\). At 4 h, using \(R/k=25/0.10=250\ \mathrm{mg}\),
The infusion stop at 4 h does not jump the state, so \(A(4^+)=82.420\ \mathrm{mg}\). During the off interval, only elimination operates. Retaining full precision in the intermediate calculations,
At 8 h the infusion restarts; the amount remains \(55.248\ \mathrm{mg}\), while the input rate changes to \(25\ \mathrm{mg/h}\). Integrating to 12 h gives
At 12 h, process the stop and the bolus at the same boundary: set the infusion rate to zero and add 50 mg once. The infusion stop does not change the amount, so \(A(12^+)=119.454+50\approx169.454\ \mathrm{mg}\). This calculation checks both the event ordering and the distinction between a rate change and a bolus jump.
Simultaneous Events and Exactly-Once Processing
A single time can contain multiple records. For example, one infusion can stop as another starts while one or more boluses are administered. Treat the shared time as one integration boundary, retain each event record, and process the complete group before restarting. For additive boluses, summing the doses is valid if all are assigned to the same state component and have the same deposition rule. If their destinations or rules differ, apply each specified jump separately.
Worked Example: Auditing Simultaneous Events
Suppose a hypothetical amount state is \(A(6^-)=40\ \mathrm{mg}\). At 6 h, infusion 1 stops at \(12\ \mathrm{mg/h}\), infusion 2 starts at \(8\ \mathrm{mg/h}\), and two separately recorded boluses of 30 mg and 20 mg enter the same state. The elimination rate constant is \(0.05\ \mathrm{h}^{-1}\). Determine the post-event amount, the active infusion rate, and the initial post-event rate of change.
The two bolus identifiers are both processed at 6 h, giving \(A(6^+)=40+30+20=90\ \mathrm{mg}\). The old infusion is off and the new one is on, so the post-event infusion rate is \(8\ \mathrm{mg/h}\), not \(12+8\ \mathrm{mg/h}\). The post-event derivative is therefore
The event ledger should show one processed stop, one processed start, and two processed bolus records at 6 h. The amount increases by exactly 50 mg across the event; changing the infusion rate contributes no instantaneous amount. Counting the records and checking the jump separately would expose either a duplicated bolus or an incorrect addition of both infusion rates after the boundary.
Diagnostics, Alternatives, and Limitations
An event schedule should be checked independently of the solver's success flag. Useful checks include confirming that times are finite and within the simulation interval, that each infusion stop follows its start, and that every start has the intended stop. After integration, compare the set of processed identifiers with the scheduled set: every required event should have a processed count of one. For an amount state, also check the jump balance at each bolus boundary, \(A(t_j^+)-A(t_j^-)=FD\) for a single bolus under the stated deposition rule.
Common warning signs include a dose-sized discontinuity at an unexpected time, an amount jump twice the prescribed dose, an infusion continuing after its stop, or a missing restart after a boundary. A small numerical tolerance may be needed when comparing times generated by arithmetic, but it should not be used to merge distinct scheduled events merely because they are close. Prefer a canonical event-time table or explicit event identifiers; otherwise a broad tolerance can collapse genuinely separate doses.
The main alternative is to use a solver with native event handling that can locate or stop at specified times and apply user-defined updates. This can be convenient, but the implementation still needs explicit semantics for left- and right-hand inputs, jump maps, simultaneous events, and exactly-once processing. Another approach is to approximate a bolus by a very short infusion or to use tiny fixed steps around a boundary. Those choices change the modeled input or add unnecessary computation; they are not substitutes for applying a discrete jump when the intended dose is a bolus.
Restarting at every discontinuity is a robust general strategy, especially when a solver's stages would otherwise evaluate the derivative on both sides of an abrupt change. A specialized solver may manage these boundaries internally, but its event conventions must be verified. In either approach, event scheduling does not guarantee that the pharmacological model is correct, and a clean event ledger does not by itself establish adequate integration accuracy. It verifies the timing and accounting layer, not every aspect of the PBPK simulation.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- A 40 mg bolus is given at 2 h, an infusion starts at 2 h, and it stops at 7 h. What are the integration boundaries, and what changes at 2 h?
- An infusion stops at 5 h while a bolus is given at 5 h. Does the infusion stop change the state amount instantaneously?
- At one shared time, two boluses of 10 mg and 15 mg enter the same amount state. What is the total state jump?
- Why is sorting event times insufficient unless event identifiers are retained?
- What should be checked if the code reports that all scheduled times were visited but one bolus is missing from the final event ledger?
Solutions. (1) The boundaries are 2 and 7 h, together with the simulation start and end. At 2 h, apply the 40 mg jump and set the infusion active before integrating onward. (2) No; the amount is continuous across an infusion stop, although the right-hand input rate changes. (3) The total jump is \(10+15=25\ \mathrm{mg}\), applied once. (4) Grouping or deduplicating times without retaining records can erase simultaneous events, while processing a time without tracking its records can duplicate them. (5) Compare scheduled and processed identifiers, not only times; then inspect the event loop's boundary conditions and processed-state bookkeeping to locate the skipped record.