ODE solver: initial or solution times that depend on parameters

[Bill Gillespie]

I'm hitting a road block to implementing PK models that include lag times or infusions that end at times that are model parameters. The lag times and times at which input rates change result in discontinuities in the ODEs. This results in a situation where I want to stop and restart the ODE solver at those times to get the solver to correctly deal with the discontinuity. Because those times depend on parameters I cannot use the current solvers because they require that the initial and solution times be data.

I am usually not interested in the solution at the parameter-dependent lag or rate change times. As a result there is no reason to calculate gradients w.r.t. those times. They are just nuisance "parameters" needed to make the ODE solver more correctly calculate the solutions at times relevant to observations. 

Is there any way to coerce those parameter-dependent time values to be treated as data other than hacking the C++ code to force coercion from var to double?

More generally we may want to reconsider the current limitation on the initial or solution times in the ODE solvers.

Thanks,

Bill

[Sebastian Weber]

Hi Bill!

lag times are a pain!

As you write, lag times influence the parameter values of the PK system quite a bit. Hence, the gradients wrt to the lag times are needed to my understanding; casting them away in C++ would be dangerous, I think.

So what you can do (this is in the archives of stan users), is to pass the tlag into the rhs of the ODE as a parameter. The rhs function then needs to handle the dosing event itself. Restarting the integrator at an unknown time is not possible, but you can just let the integrator handle the dosing discontinuity by himself. I know this is not ideal for the integrator, but the only solution I am aware of.

What I would attempt to do is to first fit a transit cmt Model first on some data which is informative about the tlag and then simply fix the lag and rerun with the full data. This is a hack, yes, but it should work. The transit cmt model can be written in Stan fully analytic, since you have all those nice gamma functions available.

I hope that helps. In case others have better ideas, I am very interested to hear these.

Best,
Sebastian

[Bill Gillespie]

Hi Sebastian,

Thanks for the ideas. 

To be clear I'm not proposing to neglect the gradient w.r.t. the lag times. Lag times would still be passed to the solver as part of the parameter array. I just want to incorporate additional initial and solution times that depend on such lag times, so the solver doesn't need to struggle with the discontinuities. In the vast majority of cases, that will not affect the true solution of the ODEs at times relevant to the likelihood. As a result the gradient with respect to the initial and solution times passed to integrate_ode would not be required even when they depend on model parameters.

Thanks,

Bill

[Bob Carpenter]

Did anyone ever answer this? I'm not an expert in this area, but
if there's a way to do it, we should put it in the manual or in an
example case study.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Sebastian Weber]

Hi!

There is no convenient solution right now to my knowledge. The best way to solve it at the moment is to use what is called a transit-cmt modeling approach. However, this is tedious if you have a multiple-dosing instead of a single dosing situation. Although, I have a plan in my head of how to make multiple-dosing with a transit-cmt work in Stan.

Another solution is to pass the lag time as parameter into the ODE RHS, but then the ODE solver has to deal with a discontinuity which is not a good idea to do (but should work).

A convenient solution for users would be an event handling mechanism for ODEs in Stan which is capable to allow events which happen at time-points which depend on data and parameters. I do not know if this is possible. The problem is that these dosing events cause discontinuous shifts in the state. However, often dosing occurs in cmt=1 while we observe another cmt such that for this case things maybe OK. Maybe Michael has a clear view on this.

Sebastian