Implementing dosing regimens in diffeq

Specific Domains Biology, Health, and Medicine
1

Hey all,

I’ve been working to implement a complex dosing regimen consisting of both IV bolus doses and continuous IV infusions.

An example would be:
t = -80, dose = 1000, rate = 0
t = -2, dose = 2500, rate = 100
t = 23, dose = 500, rate = 0
t = 23, dose = 4500, rate = 50
etc…
bolus doses are represented by rate = 0, while infusions do have rate > 0.

I decided on creating a function to return the rate over time. This would mean that for any timepoint t, the function returns the total rate for all the different doses given.
In the above example, if t = 23, this would mean that the rate is 100 + 500 + 50 = 650. If we step an Inf small time step further we would get a rate of only 50 (since both the first infusion and the bolus dose have ended).

This is put into a differential equation like so:

function dCdt(C, p, t)
    CL, V = p
    rate = get_rate_over_t(t, dosage_info)
    I = rate / V
    I - (CL / V) * C

I then define the ODEProblem normally and solve with explicitly telling the solver to step to the dosing time steps like so:

solution = solve(prob, tstops=[-80, -2, 23, etc..])

When I plot the solution however I can see that the dosing events from the bolus doses are not registered. The infusions are displayed as expected.
When I print t, rate, C, and I - (CL / V) * C (dC) from dCdt during solving I see that the solver does indeed step through the timepoints correctly and also finds the rates correctly. The solver however steps through bolus timepoints multiple times, first calculating a correct jump in dC and C, in the end settling on a high dC, but a moslty unaffected C.

I figured this could be because the solver thought this was an abberation (since it is so specifically linked to a absolute timepoint) and just decreased the effect of dC. A way to have the bolus dose show up was to set a very high rate for the bolus doses, so that the duration of the dose was not an exact point in time. Unfortunately, this makes the rate per bolus dose a hyperparameter that significantly affects the resulting C for the solution.

Is there a way to program these sudden shifts in DifferentialEquations.jl? It would be best if I could keep them linked to a specific point in time.

Anyone have suggestings or experience with this?

2

If you’re looking to do nonlinear mixed effects modeling, I would highly recommend https://pumas.ai/ since it has a whole system for dosing (Dosage Regimens, Subjects, and Populations · Pumas) in nonlinear mixed effects models. Note licenses are free for academics.

Otherwise you can hack together what you need from the callback system: Callbacks and Events

3

Thanks!

Figured out how to do it using a TimePresetCallback.
I’ve also been looking into Pumas, but looking to hold the reins and learn by doing for now.
Pumas looks promising compared to nonmem though!

1 Like