Consumer-resource model with dynamic growth rate in R with deSolve

I have the following 3 main pieces of code that first plots rainfall, then plots rainfall's effect on prey (resource) growth rate, then plots consumer-resource (herbivore-plant) dynamics using a constant growth rate. My goal is to implement the dynamic growth rate (via the equation dg into the consumer-resource model. My end goal is a graph of consumer-resource population dynamics over time with the dynamic growth rate.

First piece of code (plotting rainfall):

### Rainfall plot ###
t = seq(0, 50, 1) # time period
avgrain = 117.4 # average rainfall 
A = 100 
w = 0.6 
phi = 0.1 

rain = avgrain + (A*sin((w*t)+phi)) # rainfall function
plot(t, rain, type="l", xlab="time", ylab="Rainfall") # rainfall plot

Second piece of code (plotting rainfall's effect on resource growth rate):

### Rainfall's effect on growth rate, g ###
ropt1 = 117.4 # optimal rainfall for resource growth
s1 = 120 # standard deviation for resource growth rate as a function of rainfall

dg = exp(-(rain - ropt1)^2/(s1^2)) # rain's effect on plant growth rate
plot(t, dg, type="l", xlab="time", ylab="Plant growth rate as a function of rainfall")

Third piece of code (plotting consumer-resource dynamics - this is where I am trying to implement the dynamic growth rate created above, dg, instead of the static growth rate, g):

### Consumer-resource model as a function of time ###
library(deSolve)
states <- c(r=1, # resource (plant population) state variable
            c=1) # consumer (herbivore population) state variable

parameters <- c(g=1, # resource growth rate )
                K=25, # resource carrying capacity
                a=0.5, # consumer attack rate (between 0-1)
                h=1, # consumer handling time
                e=0.9, # consumer conversion efficiency
                m=0.5) # consumer mortality rate


function1 <- function(times1, states, parameters) {
  with(as.list(c(states, parameters)),{
   # rate of change of state variables
    dr = g*r*(1-(r/K))-((c*a*r)/(1+(a*h*r)))
    dc = ((c*e*a*r)/(1+(a*h*r)))- c*m
    
    # return rate of change
    list(c(dr, dc))
  }) 
}

times1 <- seq(0, 100, by = 1)

out1 <- ode(y = states, times = times1, func = function1, parms = parameters, method="ode45")

plot(out1) # plot state variable change across time

So, essentially, at each time step, I want the consumer-resource dynamics to be updated according to the growth rate at that time step. Is this possible? If so, how? Thank you in advance for your kind response.

Solution

You can directly include the rainfall equations in the model function. Here t is always the actual time step. They appear also in the plot as "auxiliary variables" if we add them as 3rd and 4th element of the return value after the vector of derivatives.

If dg is a dimensionless value, it can be multiplied with the growth rate, or in case it is a growth rate, replace g with dg. Please check this!

The additional parameters from the rainfall submodel may remain global or can be moved in the parameter vector. I also renamed some identifiers (the ones ending with 1) and changed the solver to lsoda instead of ode45.

library(deSolve)
states <- c(r=1, # resource (plant population) state variable
            c=1) # consumer (herbivore population) state variable

parameters <- c(g=1, # resource growth rate )
                K=25, # resource carrying capacity
                a=0.5, # consumer attack rate (between 0-1)
                h=1, # consumer handling time
                e=0.9, # consumer conversion efficiency
                m=0.5,  # consumer mortality rate
                s1=120,
                avgrain = 117.4, # average rainfall 
                A = 100, 
                w = 0.6, 
                phi = 0.1, 
                ropt1 = 117.4, # optimal rainfall for resource growth
                s1 = 120 # sd for resource growth rate ...
                )

# note that "t" is only one time point from times
model <- function(t, states, parameters) {
  with(as.list(c(states, parameters)), {
    # rainfall time series function
    rain <- avgrain + (A*sin((w*t)+phi)) # rainfall function
    # functional response equation
    dg <- exp(-(rain - ropt1)^2/(s1^2)) # rain's effect on plant growth rate

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

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

## lsoda is faster and more robust that ode45
out <- ode(y = states, times = times, func = model, parms = parameters, method="lsoda")

plot(out) # plot states and auxiliary variables across time