[R-sig-dyn-mod] Altering params for Simecol model at time step

Hi, all. 

I am working on a simecol model for which I need to be able to change the value of a parameter at a particular time step. 

Here?s the model: 

sir <- new("odeModel",
                 main = function(time, init, parms, ...){
                   with(as.list(c(init,parms)),{

                     dS <- -beta*S*I
                     dI <- beta*S*I-I*gamma
                     dR <- I*gamma
                     list(c(dS,dI,dR))
                   })},

                 parms = c(beta=0.18,gamma=0.167),
                 times = c(from=0,to=300,by=0.01),
                 init = c(S=10000000,I=25,R=0),
                 solver = "lsoda"
)

I?d like to be able to lower beta by 0.1 at, for example, time=60. 

It looks like the events method in deSolve would accomplish this, but I?m not sure how to integrate the deSolve method into the simecol model.

Thanks for any suggestions.

Hi,

please excuse the delay, it is weekend ;-)

Events are used to manipulate state variables, but you need to
change parameters. We call this a "forcing function" or an "input".

Inputs are, in principle, explained in Section 5.3 of the simecol
package vignette. However, instead of approx or approxTime I would
recommend to use approxfun, that is now much faster:



sir <- new("odeModel",
   main = function(time, init, parms, inputs, ...){
     with(as.list(c(init, parms)), {
       beta <- inputs(time)
       dS <- -beta*S*I
       dI <- beta*S*I-I*gamma
       dR <- I*gamma
       list(c(dS, dI, dR), beta=beta)
     })
   },
   parms = c(gamma=0.167),
   times = c(from=0, to=300, by=1),
   init = c(S=10000000, I=25, R=0),
   inputs = approxfun(c(0, 60, 300),
     c(0.18, 0.08, 0.08), method="constant", rule=2),
   solver = "lsoda"
)


sir <- sim(sir)
plot(sir)


If you need more than one input, make inputs a list:

inputs = list(
   beta=approxfun(c(0, 60, 300), c(0.18, 0.08, 0.08),
     method="constant", rule=2),
   gamma = approxfun ......
),


then call it like this:

beta <- inputs$beta(time)


Another hint: "by" should not be too small for performance reasons, but
on the other hand small enough so that solver does not jump over the
inputs. In your case, it is safe to increase "by" to 1.0

Does this help you?

Thomas

Dr. Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

E-Mail: thomas.petzoldt at tu-dresden.de
http://tu-dresden.de/Members/thomas.petzoldt

Thomas, 

Thanks for this, it?s really clear and does exactly what I was hoping to accomplish. 

Thanks also for anticipating my next question (multiple inputs). 

Cheers,
GD

sir <- new("odeModel",
 main = function(time, init, parms, inputs, ...){
   with(as.list(c(init, parms)), {
     beta <- inputs(time)
     dS <- -beta*S*I
     dI <- beta*S*I-I*gamma
     dR <- I*gamma
     list(c(dS, dI, dR), beta=beta)
   })
 },
 parms = c(gamma=0.167),
 times = c(from=0, to=300, by=1),
 init = c(S=10000000, I=25, R=0),
 inputs = approxfun(c(0, 60, 300),
   c(0.18, 0.08, 0.08), method="constant", rule=2),
 solver = "lsoda"
)


sir <- sim(sir)
plot(sir)