popPD model estimation with step function covariate

[wrholmes]

I have a question about improving efficiency of estimation (or really the ODE simulation) with a covariate that discontinuously flips from 0 to 1 in time.

Here is the backdrop. I’m simulating a rollover study where people that were initially in Placebo condition are switched to taking the drug after an initial time period. Due to the nature of the data I only know that they switched from placebo to drug. As such, I effectively have a covariate that switches discontinuously in time (placebo = 0, drug = 1 with a switch at some time).

In this context, I am fitting a population model with IIV on some parameters with a simple ODE model at the individual level. Standard popPD sort of thing. Some of the parameters that have IIV on them are associated with this discontinuous covariate.

I am finding that when I try to estimate this model, the estimation is very slow relative to what’s expected (order of magnitude slow). This is almost certainly due to the discontinuity in the covariate and its impact on the ODE solver. The ODE model itself is quiet simple and not at all stiff. The data itself is also fairly regular. Further, the number of estimation steps is reasonable. So I’m pretty sure it is the covariate discontinuity.

Any suggestions to address this? I’ve used the SciML toolkit outside of Pumas in the past and am fairly familiar. I assume you all have plugged in some stiff solvers. However I know that Pumas is making some choices in the backend that I don’t know about. Any suggestions on how to handle this sort of issue? Is it as simple as specifying a particular solver, and if so is there one recommended for this case where stiffness is local in time?

[ChrisRackauckas]

tl;dr, the solution

Probably the simplest way to handle this is to setup this as a dose time optimization. Setup a new ODE variable x where x’ = 0, and then setup this switch as just a regular dose with amt=1, then just use x in the other equations. Because dosing always introduces a discontinuity, the Pumas system is already equipped to handle that efficiently, and then this becomes the same as any other dosing time optimization.

A bit of an explanation of the behavior

But to give a flavor of what’s going on here, ODE solvers generally cannot keep high accuracy when stepping over a discontinuity. Here’s a snippet from Hairer’s Solving Ordinary Differential Equations v1:

What you see here is an ODE which has a discontinuity at t=0.62, and what the solver will do is suggest a dt that goes over the discontinuity, reject that step (the star), then use a smaller dt, reject that, then smaller dt, accept that, then do this over again, until it takes a really tiny step near the discontinuity to get over it, and recovers. This is because stepping over a discontinuity breaks the higher order convergence rules, so error adaptivity will handle discontinuities automatically but it will do so by trial and error.

But if you tell the system where the discontinuity is, it will step directly on it, and then update its state, and continue on. That’s a ton faster. Pumas does that internally for all of its dosing mechanisms, and so this is why it’ll be more efficient to just hook into that (and not to mention more robust in estimation and error control).

Some More Musing for Future Notes

Fully handling it on the user end should be making a callback that then switches a parameter during the run, so using a PresetTimeCallback that is a bit modified. What you’d do is, put a callback into the solve which has an initialization phase that does add_tstop!(integrator, integrator.p.switch_time), and then you’d modify integrator.p for the discontinuous switch. This slide in a recent workshop I gave shows how to define a callback init that adds a tstop:

A Deep Dive Into DifferentialEquations

If you do that, then it’ll set the solver to automatically handle the discontinuity time. But the difficulty you’ll run into there is passing the callback into the Pumas infrastructure.

But don’t go this route unless you’re looking to do something crazy. Which Pumas of course allows, but the 99% user shouldn’t do this.

[devmotion]

Hi!

If I understand you correctly, you use a time-varying covariate in your model to describe the switch from placebo to drug?

Pumas automatically ensures that ODE solvers step to the time points of the discontinuities in the covariates and that ODE solvers know about the discontinuous changes at these points (the need for this is nicely explained in the post by @ChrisRackauckas above). Hence, if you model this switch as a time-varying covariate you don’t need to make any of the additional adjustments suggested above.

@wrholmes

I am finding that when I try to estimate this model, the estimation is very slow relative to what’s expected (order of magnitude slow).

If you can provide an example of the model and the slow estimation, we can look into what’s causing the problem you’re observing.

@wrholmes

if so is there one recommended for this case where stiffness is local in time?

The hybrid ODE solvers used by Pumas by default switch between a stiff and non-stiff solver during the ODE simulation, depending on whether the problem appears to be stiff or non-stiff at the current integration time point. Thus this choice/switch is already local in time. You can read more about the (default) ODE solvers and their (default) settings in Differential Equations in Pumas and Analytical Solutions and Differential Equations · Pumas.