我试图复制一个湖泊食物网模型here。这个模型代表了两个食物链（沿岸的和远洋），由一个顶级捕食者（鱼）连接。我已经编写了这个模型，但是当我运行它2-3个时间步后，这个模型产生NaN。我已经通过我的代码多次寻找括号等的问题，但是不能找到问题。
如果我将fish的初始丰度设置为0，模型就会运行，所以我假设问题一定出在模型的fish组件上。
以下是公式：
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

如果有人能看到我哪里出了问题，那就太感谢了！