Time varying input for ODE solving (aka (?) forced ODE)

Modeling
1

Hi !
I’m currently working on a quite simple model for stan. It consists of time-series data and an ODE + noise model.

The equation for the ode relies on an input (the amount of “inducer”). In the current code (given below), the input is constant and given as data. What I would like to do is to give a time-varying input (also as data), for example two arrays, each one for the amount of one of the inducers at each time point, and solve the ODEs with it. Note that the amount of inducers is a step function After searching, i found several issues/post related to either using inference functions or the floor function to get inducers values from an array, but none actually gave a solution. How would you do this ?

The code

[details=Click to show the code] functions {
real hill_neg(real x,real theta,real eta) {
real aux=theta^eta;
return aux/(aux+x^eta);
}
real[] ode(real t,
real[] y,
real[] param,
real[] x,
int[] x_i) {
real dydt[2];
dydt[1] = hill_neg(y[2](1-x[1]),param[1],param[2])-y[1];
dydt[2] = hill_neg(y[1](1-x[2]),param[3],param[4])-y[2];
return dydt;
}
}
data {
int<lower=1> T;
real y[T,2];
real t0;
real ts[T];
real inducer[2];
}
transformed data {
int x_i[0];
}
parameters {
real<lower=0.0,upper=1> theta[2];
real<lower=1,upper=4> eta[2];
real<lower=0.00,upper=1> sigma;
}
transformed parameters {
real params[4] = {theta[1],eta[1],theta[2],eta [2]};
real vrnce = sigma^2;
}
model {
real y_hat[T,2];
//print("target (before) = ", target());
y_hat = integrate_ode_rk45(ode, {0.0,0.0}, t0, ts, params, inducer, x_i);
//print("y_hat = ",y_hat);
for (t in 1:T) {
y[t] ~ normal(y_hat[t], vrnce);
}
//print("target (after) = ",target());
}[/details]

Thanks a lot for your help

2

I think you’re looking for this post by @wds15: Forced ODEs, a start for a case study?

I wouldn’t say any time-series model is simple haha. Always room for surprises. Good luck!

3

Hi !

Thanks a lot for pointing me to this thread ! After reading it, I have some questions:

  1. I consider stan as a kind of magicl black box, the internals of which I don’t really understand, but could forcing ODE with a discontinuous function (step in my case) cause issues with the way stan uses the derivatives to do its computation ?
  2. Given that in my case I would know that the input can change only at each tick of clock (for example a each integer time t only), I would like to improve on the provided example’s complexity (which has to go through the input array each time) by doing something like x[floor(t)] to get the value of input x at time t, but I have read that there is no such function ? Is that correct ? Can it be by passed in any way ?

Thx !

4

If your inputs are step functions, any chance it makes sense to parameterize them as sigmoids or something and then you don’t have to pull your forcing function from data?

Just Googling around, looks like the answer in Matlab-land is yes, you should take care with step functions: https://www.mathworks.com/matlabcentral/answers/53190-ode-solver-with-heavaside-function . The tldr; in that thread is that the Matlab person recommends integrating up to the step, stopping the ODE solver, and then restarting it on top of the step. You could do this in Stan. It’s probably worth researching a little more than this though to understand what’s happening though and what your options are. Looks like a way to dodge the non-existent derivatives.

Yeah, there’s no floor function. In terms of interpolator functions for this use case in the ODEs, it’s an open issue: https://github.com/stan-dev/math/issues/58

5

Actually, I know for a fact @wds15 has ODEs with steps like this. Here’s some code he put together for some PKPD systems: https://github.com/stan-dev/example-models/tree/feature/issue-72-stan-pkpdlib/misc/pkpd

In that model patients are given doses at various timepoints. From the text it looks like he does the ODE restarts.

6

In short to deal with steps in ODES:

  • try to avoid if possible

  • if you have to, stop the integrator and restart at the step

  • or make the ode integrator work a bit more and place a few observation time points close by the step. Then the integrator will notice the step and just work a bit more.

However, as said, maybe you can avoid the step, possibly with a trick like Andrew posted recently on how to write down models when the linear slope changes at some point. I forgot the exact title, but what he has done there is a nice continuous approximation to steps.

Sebastian

7

Hi !

Thanks for taking the time to answer, There are still some points that remain unclear:

  1. Concerning the stop&start methods, I guess that it incurs a cost that can quickly get non negligible, do you have any opinion on that ?

  2. When you say “observation points” you mean the values given in the time vector to the ODE solver right ? If so, my understanding of this was that the solver would choose its own time points and only interpolates the given ones ? Is that wrong ? If not, would the method you are proposing still work ?

  3. To avoid steps, I guess I would have to use some kind of combination of sigmoids to continuously approximate my actual input. Intuitively, this kind of interpolation are at least linear to compute in the number of changes of the input, do you have any tips on how this cost would compare to the start & stop one ?

Thx a lot for your help :)

8

Glad to help!

1 & 2 are related. Stop & start becomes inefficient because what you say in 2 is correct (the ODE solver runs through space choosing its own time points, and just interpolates to spit out values at yours). If it stops and starts, it has to actually stop at your point, and when it starts again it has to figure out its time step and all from scratch.

3 if you have multiple steps, then probably best to just do the stop-restart thing.

Give it a go and see how it works :D. Might just all work fine with the stop-start!

9

Have a look at the wiki from sundials

http://sundials.wikidot.com/integrating-over-discontinuities

There its described how to deal with discontinuities, but really just get started!

10

I am still balancing the learning effort and the new capabilities of STAN. The build-in integrator and the possibility of estimating parameters from ODE is certainly interesting. I have playing with the predator-prey model in

https://mc-stan.org/users/documentation/case-studies/lotka-volterra-predator-prey.html

This example is very clear. Congratulations.

Toward the end, at “Exercises and Extensions #8”, it is proposed to evaluate the effects of external “covariates”, as Temperature, which typically change at smaller scales that the process of interest, are densely sampled and remains unaffected by the process. Thus, the example in

seems not to be the same case (…and it is hard to understand for me…).

Can someone provide a detailed (=beginner’s level) example of how to deal with this type of external variables?

For those using R and the deSolve library, I am looking for something that emulates aproxfunc.

Thanks
miquel

11

Unfortunately right now interpolation in Stan is difficult. We want to change that, but until then this tutorial is probably about as clear as it gets.

12

Thank you Ben