I am using the deSolve
package to solve a differential equation describing predator-prey dynamics. As an example, below is a simple L-V predator-prey model. I would like some of the parameters in the model to vary through time. I can vary state variable (e.g. prey density) no problem using the event
argument in the ode
function.
But I cannot use the event
argument to alter parameters.
Here is the simple L-V model with no events added (works fine)
# Lotka-Volterra Predator-Prey model from ?deSolve::ode
# define model
LVmod <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
Ingestion <- rIng * Prey * Predator
GrowthPrey <- rGrow * Prey * (1 - Prey/K)
MortPredator <- rMort * Predator
dPrey <- GrowthPrey - Ingestion
dPredator <- Ingestion * assEff - MortPredator
return(list(c(dPrey, dPredator)))
})
}
# parameters
pars <- c(rIng = 0.2, # rate of ingestion
rGrow = 1.0, # growth rate of prey
rMort = 0.2 , # mortality rate of predator
assEff = 0.5, # assimilation efficiency
K = 10) # carrying capacity
# initial densities (state variables)
yini <- c(Prey = 1, Predator = 2)
# time steps
times <- seq(0, 200, by = 1)
# run model
out <- ode(yini, times, LVmod, pars)
## plot
plot(out)
Here is the L-V model with state variable Prey
multiplied by some rnorm()
every 5 timesteps (works fine).
# add prey every 5 timesteps using events
add_prey <- function(t, var, parms){
with(as.list(var),{
Prey <- Prey * rnorm(1, 1, 0.05)
return(c(Prey, Predator))
})
}
# run ode - works fine
out <- ode(y = yini,
times = times,
func = LVmod,
parms = pars,
method = "ode45",
events = list(func = add_prey, time = seq(0, 200, by = 5)))
plot(out)
Here is my attempt to increase K every 5 timesteps (does not work)
# vary K through time
add_k <- function(t, var, parms){
with(as.list(var),{
K <- K + 2
return(c(Prey, Predator))
})
}
# run ode
out <- ode(y = yini,
times = times,
func = LVmod,
parms = pars,
method = "ode45",
events = list(func = add_k, time = seq(0, 200, by = 5)))
Which produces this error:
Error in eval(substitute(expr), data, enclos = parent.frame()) :
object 'K' not found
Based on the error K is not being passed to add_k
, in add_k
the line with(as.list(var)
is obviously only accessing the variables Prey
and Predator
. In the ode
and event
helpfiles I can only find information regarding altering state variables (Prey
and Predator
in this case), and no information about altering other parameters. I am new to ODEs, so maybe I am missing something obvious. Any advice would be much appreciated.
Events are used to modify state variables. This means that at each event time, the ode solver is stopped, then states are modified and then the solver is restarted. In case of time-dependent *parameters, this can be done much easier without events. We call this a "forcing function" (or forcing data).
A modified version of the original code is shown below. Here a few explanations:
approxfun
is a function that returns a function. We name itK_t
that interpolates between the data timestepst
and the parameter valuesK
. The argumentrule = 2
is important. It describes how interpolation is to take place outside the range oft
, where2
means that the closest value is returned, because some solvers overshoot the end of the simulation and then interpolate back. Themethod
can beconstant
orlinear
, whatever is more appropriate.K_t
can be added to the function as an optional argument at the end. It is also possible to define it globally or to include it inparms
if it is defined as a list, not a vector.LVmod
function checks if this additional parameter exists and if yes overwrites the default ofK
.LVmod
function we return the actual value ofK
as an optional (auxiliary) variable, so that it is included in the model output.