具有 deSolve 的 R 中具有动态增长率的消费者资源模型

Question

我有以下 3 段主要代码，首先绘制降雨图，然后绘制降雨对猎物（资源）增长率的影响，然后使用恒定增长率绘制消费者资源（食草动物植物）动态图。我的目标是将动态增长率（通过等式 dg 实施到消费者资源模型中。我的最终目标是消费者资源人口动态随时间变化的图表。

第一段代码（绘制降雨量）：

### 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

第二段代码（绘制降雨对资源增长率的影响）：

### 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")

第三段代码（绘制消费者资源动态 - 这是我尝试实现上面创建的动态增长率 dg，而不是静态增长率 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

因此，基本上，在每个时间步，我希望根据该时间步的增长率更新消费者资源动态。这可能吗？如果是这样，如何？预先感谢您的友好回复。

Answer 1

您可以直接在 model 函数中包含降雨方程。这里 t 始终是实际时间步长。如果我们将它们作为 return 值的第 3 和第 4 个元素添加到导数向量之后，它们也会作为“辅助变量”出现在图中。

如果dg是无量纲值，可以乘以增长率，如果是增长率，将g换成dg。请检查这个！

来自降雨子模型的附加参数可以保持全局，也可以在参数向量中移动。我还重命名了一些标识符（以 1 结尾的标识符）并将求解器更改为 lsoda 而不是 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

具有 deSolve 的 R 中具有动态增长率的消费者资源模型

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

r

ode

population

desolve