Replicating an ODE food-web model from a manuscript

Question

I am trying to replicate a lake food web model published here. The model represents two food chains (littoral vs pelagic), connected by a top predator (fish). I have coded the model but when I run it after 2-3 time steps the model produces NaN. I have been through my code many times looking for issues with parentheses etc. but can't find the problem.

If I set the fish initial abundance to 0 the model runs, so I assume the issue must be with the fish component of the model.

Here are the equations:

Ap = pelagic resource, Z = pelagic zooplankton, Pp = pelagic predator, F = fish, Al = littoral resource, I = invertebrate, Pl = littoral predator.

Here is my attempt at coding the model:

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

If anyone can see where I have gone wrong it would be much appreciated!

Not an expert in the field, but at first look i cant' see a derivative depending on `Time`. is it ok? Are all derivatives constant with respect to time? — Ric, Dec 01 '22 at 20:33
I am no expert either. But as far as I understand it is correct regarding time.. I have set up and run similar, but more simple models (such as basic Lotka-Voltera) in this is exact way and they have run as expected. — flee, Dec 01 '22 at 20:48

Ben Bolker · Accepted Answer · 2022-12-01T21:29:18.433

5

tl;dr there's a typo in your 'fish' equation (you multiplied instead of adding in the denominator). (I missed that you already said in your question that you had localized the problem to this equation! Nevertheless, maybe the debugging procedure below will be helpful ...)

I ran the model with a much smaller delta-t over a smaller time horizon to try to see which state variables were problematic.

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)

It looks like the fish population is exploding and the PP/PL populations are crashing (which would follow naturally from an exploding fish population; in principle if the equations are well-posed this wouldn't lead to a mathematical problem, but it's not surprising that this is causing numerical problems).

Went back and looked at the dF/dt equation, and sure enough, found the typo.

Re-running from 0 to 1000 with the corrected denominators:

edited Dec 01 '22 at 21:29

answered Dec 01 '22 at 21:23

Ben Bolker

211,554
25
370
453

Great thanks @Ben Bolker for answering. One small hint: instead of `matplotout([,1], out[,-1] ....)` one can also use deSolve's `matplot.0D(out, ...)` with the whole output object. – tpetzoldt Dec 01 '22 at 22:14
Nice. (Something I wanted to comment on somewhere, at some point, is that using `all_vals <- as.list(c(State, Pars)); attach(all_vals); on.exit(detach(all_vals))` can make debugging within the gradient function (if it needs to be done) *much* easier than the standard `with(as.list(...), { ... })` paradigm ... – Ben Bolker Dec 01 '22 at 22:20
@BenBolker Thanks for taking the time to spot my mistake. Also for the `attach` tip... the debugging for ODEs is a bit different to the normal loop-based models I have used previously. – flee Dec 01 '22 at 23:31
I agree that `with` -blocks are difficult to debug., but I would vote against the `attach`-method. In fact, I do wonder why `attach` was not yet banned and removed from R. The version with `on.exit` should be safe, but experience shows that it is often used without the necessary care and understanding and will then lead to errors that are even more difficult to trace. Not sure what to use as an alternative: use globals? Assign everything to local variables? Use closures? Ideas are welcome. – tpetzoldt Dec 02 '22 at 06:41
2

I agree that `attach()` is **usually** a bad idea (indeed, that's why it's not allowed in package code on CRAN), but I think this is a case where it works well - it's the most transparent and convenient solution. I wrote an `unpack()` function that has similar effects: https://github.com/mac-theobio/McMasterPandemic/blob/2c512006225f69628cbba813dd0f050b37d3aeeb/R/utils.R#L38-L41 – Ben Bolker Dec 02 '22 at 14:01

Replicating an ODE food-web model from a manuscript

1 Answers1