从手稿复制 ODE 食物网模型

Question

我正在尝试复制此处发布的湖泊食物网络模型。该模型代表两条食物链（沿海与远洋），由顶级捕食者（鱼类）连接。我已经对模型进行了编码，但是当我在 2-3 个时间步长后运行它时，模型会生成NaN. 我已经多次检查我的代码，寻找括号等问题，但找不到问题。

如果我将fish初始丰度设置为 0，模型就会运行，所以我认为问题一定出在模型的鱼组件上。

以下是方程式：

Ap = 中上层资源，Z = 中上层浮游动物，Pp = 中上层捕食者，F = 鱼类，Al = 沿岸资源，I = 无脊椎动物，Pl = 沿岸捕食者。

这是我对模型进行编码的尝试：

library(deSolve)

# define the model
vade_2005_model <- function(Time, State, Pars){
  
  with(as.list(c(State, Pars)), {

# pelagic components -----------------------------------------------

# resource
pel_res_dt <- (rPel * AP * (1 - (AP/(KT * q)))) - (aZP * ((Z * AP)/(AP + bZP)))

# zooplankton
pel_zoo_dt <-  (ef * aZP * ((Z * AP)/(AP + bZP))) - (dZP * Z) - (aPP * ((PP * Z)/(Z + bPP)))

# pelagic predator (PP)
pel_PP_dt <- (ef * aPP * ((PP * Z)/(Z + bPP))) - (dPP * PP) - (aFP * del * ((fish * PP)/(PP + bFA)))


# top predator - fish ---------------------------------------------
fish_dt <- (ef * ((del * aFP * ((PP * fish)/(PP * bFA))) + ((1 - del) * aFL * ((PL * fish)/(PL * bFA))))) - (dFA * fish)

# Littoral component -----------------------------------------------

# resource
lit_res_dt <- (rLit * AL * (1 - (AL/(KT * (1 - q))))) - (aIL * ((I * AL)/(AL + bIL)))

# littoral invert
lit_inv_dt <-  (ef * aIL * ((I * AL)/(AL + bIL))) - (dIL * I) - (aPL * ((PL * I)/(I + bPL)))

# littoral predator (PL)
lit_PL_dt <- (ef * aPL * ((PL * I)/(I + bPL))) - (dPL * PL) - (aFL * (1 - del) * ((fish * PL)/(PL + bFA)))

list(c(pel_res_dt, pel_zoo_dt, pel_PP_dt,
       fish_dt,
       lit_res_dt, lit_inv_dt, lit_PL_dt))
 })
}
  

# model parameters (taken from the manuscript)  
pars = c(rPel = 1.0,   # per capital growth rate of pelagic resource
         rLit = 0.8,   # per capital growth rate of littoral resource
         aZP = 1.55,    # Attack rate of zooplankton (in pelagic)
         aPP = 1.35,    # Attack rate of PP (in pelagic)
         aFP = 1.05,    # Attack rate of fish (in pelagic)
         aFL = 1.0,    # Attack rate of fish (in littoral)
         aIL = 1.45,    # Attack rate of invert (in littoral)
         aPL = 1.25,    # Attack rate of PL (in littoral)
         bZP = 0.2,    # Half saturation rate of zooplankton (in pelagic)
         bPP = 0.2,    # Half saturation rate of PP (in pelagic)
         bFA = 0.2,  # Half saturation rate of fish (in all)
         bIL = 0.2,    # Half saturation rate of invert (in littoral)
         bPL = 0.2,    # Half saturation rate of PL (in littoral)
         dZP = 0.6,    # biomass loss rate of zooplankton (in pelagic)
         dPP = 0.15,    # biomass loss rate of PP (in pelagic)
         dPL =  0.15,   # biomass loss rate of PL (in littoral)
         dIL = 0.6,    # biomass loss rate of invert (in littoral)
         dFA = 0.1,    # biomass loss rate of fish (all habitat)
         ef = 0.8,     # conversion efficiency of resource biomass into consumer biomass
         del = 0.5,    # fish preference (1 = PP, 0 = PL)
         KT = 1,    # lake carrying capacity
         q = 0.5    # prop. productivity in pelagic food chain
)      

# initial densities (assumed as not given in the manuscript)
yini <- c(AP = 0.5,
          Z = 0.5,
          PP = 0.5,
          fish = 0.5,
          AL = 0.5,
          I = 0.5,
          PL = 0.5
)

# time steps
times <- seq(0, 1000, by = 1)
  
# run model
out <- ode(yini, times, vade_2005_model, pars, method = "ode45")
out

如果有人能看到我哪里出错了，我将不胜感激！

Answer 1

tl;dr你的“鱼”方程中有一个拼写错误（你乘以而不是添加分母）。（我错过了您在问题中已经说过您已将问题定位到这个方程！不过，也许下面的调试过程会有所帮助......）

我在更小的时间范围内以更小的 delta-t 运行模型，试图找出哪些状态变量存在问题。

out <- ode(yini, seq(0, 4, length.out = 101), vade_2005_model, pars, method = "ode45")
matplot(out[,1], out[,-1], type = "l", log = "y", ylim = c(1e-6, 1e6),
        lty = 1:7, col = 1:7)
legend("topright",
       legend = colnames(out)[-1],
       col = 1:7,
       lty = 1:7)

看起来鱼类种群正在爆炸式增长，而 PP/PL 种群正在崩溃（这是鱼类种群爆炸式自然发生的结果；原则上，如果方程成立，这不会导致数学问题，但这并不奇怪这会导致数值问题）。

回去查看 dF/dt 方程，果然发现了错字。

使用更正后的分母从 0 重新运行到 1000：