Simulating ODE model for different initial conditions in R

Question

I have a model, and I want to generate random initial conditions, run the model, and save the output so that each simulation is a replicate. But I have a hard time interpreting and implementing loops (and I also know they are not always the best to use in R), so I'm struggling.

My ultimate goal is to iterate the simulation across 10 different random initial conditions, and save the output of the ODE including a column for simulation number.

First I have my random initial conditions:

library(deSolve)
states <- c(r=runif(1, min=0.1, max=25), # resource state variable
            c=runif(1, min=0.1, max=10)) # consumer state variable

Then I have my parameters and model:

parameters <- c(g=5, # resource growth rate )
                K=25, # resource carrying capacity
                a=1, # consumer attack rate
                h=1, # consumer handling time
                e=0.9, # consumer conversion efficiency
                m=0.5,  # consumer mortality rate
                avgrain = 1500, # average rainfall 
                A = 1000, 
                w = 0.6, 
                phi = 8.5, 
                ropt1 = 1500, # optimal rainfall for resource growth
                s1 = 1000, # standard deviation for plant growth rate as a function of rainfall
                ropt2 = 1000, # optimal rainfall for herbivore attack (feeding) rate
                s2 = 500, # standard deviation for herbivore attack rate as a function of rainfall
                avgtemp = 20, # average temperature 
                A_temp = 7, 
                w_temp = 0.5, 
                phi_temp = 0.5, 
                topt1 = 13, # optimal temperature for resource growth
                ts1 = 10 # standard deviation for plant growth rate as a function of temperature
                )

model <- function(t, states, parameters) {
  with(as.list(c(states, parameters)), {
    # rainfall time series function
    rain <- avgrain + (A*sin((w*t)+phi)) # rainfall function
    # temperature time series function
    temp = avgtemp + (A_temp*sin((w_temp*t)+phi_temp))
    
    # dynamic g and a equations
    dg_both <- (exp(-(rain - ropt1)^2/(s1^2))) + (exp(-(temp - topt1)^2/(ts1^2)))
    da = exp(-(rain - ropt2)^2/(s2^2))

    # rate of change of state variables
    dr <- dg_both*r*(1-(r/K)) - ((c*da*r)/(1+(da*h*r)))
    dc <- ((c*e*da*r)/(1+(da*h*r)))- c*m
    # return rate of change
    list(c(dr, dc), rain=rain, temp=temp, dg_both=dg_both, da=da)
  }) 
}

times <- seq(0, 200, by = 1) 

out <- ode(y = states, times = times, func = model, parms = parameters, method="lsoda")

Would I do this with a for loop? Thank you in advance!

Answer 1

Here one of the other approaches, mentioned by @Ben Bolker. Here we use replicate instead of a loop. This has the advantage, that we don't need to create a list() for the results beforehand.

N <- 10
res <- replicate(N, ode(y = c(r = runif(1, min = 0.1, max = 25),
                              c = runif(1, min = 0.1, max = 10)),
                        times = times, func = model, 
                        parms = parameters, method="lsoda"), 
                 simplify = FALSE)

plot(out, res)

As an additional goody, we can also plot the results using deSolve's built-in plotting function. This works of course also with res in Ben's approach. The resulting data structure can then be simplified to something like a matrix or array, either with do.call(rbind, res) as in Ben's example, or with option simplify directly in replicate.