循环遍历参数以使用 deSolve 达到平衡

Question

我直觉地与循环作斗争。我有一个简单的消费者资源模型，我想遍历资源增长率 g 的值以获得最终状态值，然后将均衡绘制为参数值的函数。这是我目前所拥有的：

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)

我想在 ode 函数之前我一直在思考 - 在这之后我应该在哪里索引 i？我希望这是有道理的。

Answer 1

您的方法有几个问题，因此我尝试重新组织它以使其运行得以通过。但是，由于您的模型显示了一个稳定的循环，它没有达到平衡。

这里有一些提示

循环应该只包含在模拟过程中发生变化的东西。固定代码段应该在循环之前。这样更容易维护，速度也更快。
首先，运行没有循环的模型，看是否可行。
然后定义一个数据结构（矩阵或数据框）来存储结果。

这是一种实现方法：

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.

还有其他几种方法，如前所述，该模型没有平衡。这里还有一个 model example，所谓的恒化器 with 平衡。

循环遍历参数以使用 deSolve 达到平衡

Looping through parameters to get equilibrium with deSolve

r

loops

ode

desolve