Looping through parameters to get equilibrium with deSolve

Question

I struggle with loops intuitively. I have a simple consumer-resource model, and I want to loop through values of resource growth rate g to get final state values to then plot equilibrium as a function of the parameter values. This is what I have so far:

param.values = seq(from = 1, to = 10, by = 1)
variable = rep(0,length(param.values))
for (i in 1:length(param.values)){ 
  state <- c(r = 1, n = 1)
  parameters = c(g = variable[i],# resource growth rate
                 d = 0.5, # n mortality rate
                 k = 5, # r carrying capacity
                 c = 1, # consumption rate of n on r
                 e = 1, # conversion efficiency for n on r
                 h = 1 # handling time n on r
  )
  function1 <- function(times, state, parameters) {
    with(as.list(c(state, parameters)),{
      # rate of change
      dr = variable[i]*r*(1 - (r/k)) - (c*n*r/(1+(h*c*r))) 
      dn = (e*c*n*r/(1+(h*c*r)))- n*d
      
      # return the rate of change
      list(c(dr, dn))
    }) 
  }
  times <- seq(0, 100, by = 1)
  
  out <- ode(y = state, times = times, func = function1, parms = parameters)
  
  sol <- out[101, 2:3] # trying to get last equilibrium value to plot against param values...
  print(sol[i])
}

plot(sol[,1] ~ param.values)
plot(sol[,2] ~ param.values)

I think I have thinks right up until the ode function - where should I be indexing i after this? I hope this makes sense.

Answer 1

Your approach had several issues, so I tried to re-organize it so that it runs through. But, as your model shows a stable cycle, it does not reach an equilibrium.

Here a few hints

• The loop should only contain things that change during the simulation. Fixed code segments should come before the loop. This is easier to maintain and faster.
• First, run the model without the loop, to see whether it works.
• Then define a data structure (matrix or data frame) to store the results.

Here one approach how it can be implemented:

library("deSolve")

## define as much as possible outside the loop
function1 <- function(times, state, parameters) {
  with(as.list(c(state, parameters)),{
    # rate of change
    dr = g*r*(1 - (r/k)) - (c*n*r/(1+(h*c*r)))
    dn = (e*c*n*r/(1+(h*c*r)))- n*d

    # return the rate of change
    list(c(dr, dn))
  })
}

state <- c(r = 1, n = 1)

parameters = c(g = 1,   # resource growth rate
               d = 0.5, # n mortality rate
               k = 5,   # r carrying capacity
               c = 1,   # consumption rate of n on r
               e = 1,   # conversion efficiency for n on r
               h = 1    # handling time n on r
)
times <- seq(0, 100, by = 1)

## first test single run of model
out <- ode(y = state, times = times, func = function1, parms = parameters)
plot(out)
## It runs and we see a cycling model. I suspect it has no equilibrium!

param.values = seq(from = 1, to = 10, by = 1)

## define a matrix where results can be stored
sol <- matrix(0, nrow=length(param.values), ncol=2)

for (i in 1:length(param.values)){

  ## replace single parameter g with new value
  parameters["g"] <- param.values[i]
  out <- ode(y = state, times = times, func = function1, parms = parameters)

  ## store result of last value in row of matrix.
  ## Note that it may not be an equilibrium
  sol[i, ] <- out[101, 2:3] # trying to get last equilibrium value to plot against param values...
  print(sol[i, ])
}

plot(sol[,1] ~ param.values, type="l")
plot(sol[,2] ~ param.values, type="l")
## We see that the model has no equilibrium.

There are several other ways and, as said, the model has no equilibrium. Here another model example, a so-called chemostat with equilibrium.