Ordinary differential equations (ODEs) - Is there any way to prevent negative values?

Question

I am trying to apply a system of ordinary differential equations (ODEs) at each spatial grid cell. Thus, each landscape cell has an associated ODE model. The number of susceptible mosquitoes (Sv), exposed mosquitoes (Se) and infected mosquitoes (St) are updated at each time step and the ODE model is coupled with an discrete agent-based model of animal movement. Here is an example to run the ODE model:

library(deSolve)

mod1 <- function(out_tab, time_step, var){
 Sv <- out_tab[time_step,c("Sv")]
  Ev <- out_tab[time_step,c("Ev")]
  Iv <- out_tab[time_step,c("Iv")]

  Nh <- out_tab[time_step,c("Nh")]
  Ih <- out_tab[time_step,c("Ih")]
  bv <- 100
  dv <- 0.07
  betav <- 0.33
  av <- 0.5
  muv <- 0.1

  mod <- function(times, states, parameters) {
    with(as.list(c(states, parameters)), {

      dSv <- bv - dv*Sv - betav*av*(Ih/Nh)*Sv
      dEv <- betav*av*(Ih/Nh)*Sv - dv*Ev - muv*Ev
      dIv <- muv*Ev - dv*Iv

      return(list(c(dSv, dEv, dIv)))

    })
  }

  states <- c(Sv = Sv, Ev = Ev, Iv = Iv)
  input_parameters <- c(bv = bv, dv = dv, betav = betav, av = av, muv = muv)

  ## Solve ODEs
  out <- ode(func=mod, y=states, times=seq(time_step, time_step + 1, by=1), parms=input_parameters, method="iteration")
  out <- as.data.frame(out)

  out_tab[time_step + 1, c("Sv", "Ev", "Iv")] <- out[dim(out)[1], c("Sv", "Ev", "Iv")]
  # out_tab[time_step + 1, c("Nh")] <- var
  # out_tab[time_step + 1, c("Ih")] <- var

  return(out_tab)

}
out_tab <- data.frame(Sv = 10, Ev = 3, Iv = 2, Nh = 2, Ih = 1)
Nh <- c(5, 1, 8, 0, 5)
Ih <- c(2, 0, 4, 0, 1)

for(time in 1:2){
  out_tab <- mod1(out_tab = out_tab, time_step = time)
  ## print(out_tab)
  out_tab[time + 1, c("Nh")] <- Nh[time + 1]
  out_tab[time + 1, c("Ih")] <- Ih[time + 1]
}

At time = 2 (in days), Ih (number of animals in a cell) is equal to 0 so that the value for Ev is negative. Is there any way to prevent negative values from ODEs? I used the Method "iteration" because the ODE model is incorporated in a discrete agent-based model.

Answer 1

There are a few major issues here; to cut a long story short you have defined your model system incorrectly to use with method = "iteration", so no amount of light tinkering will get you sensible results. I'll get to this in the second part of my answer, but first, I'll answer your original question.

Using Events to Obtain Non-Negative State Values

You can force non-negative population sizes in deSolve using events. I suggest you read the deSolve documentation further, as there are many ways to trigger events, but we'll integrate one into your code that is triggered during each time step. Because it is time-based, it'll need some sort of maximum time to refer to, so I've created a maxtime value that you can also use in your for-loop. You should define this value somewhere accessible to your other functions; I simply put it before your mod1 function declaration.

# Here we declare the maximum time to which your system will evaluate
maxtime <- 2

# This is where we define your event function
# Add this directly above your call to ode()
  posfun <- function(t, y, parms){
    with(as.list(y), {
      y[which(y<0)] <- 0  
      return(y)
    })
  }

# Here's your original call to ode(), with a small addition
# Notice that we added events; iteration is missing, more on that later
  out <- ode(func=mod, 
             y=states, 
             times=seq(time_step, time_step + 1, by=1), 
             parms=input_parameters, 
             events=list(func = posfun, time = c(0:maxtime)))
  
  out <- as.data.frame(out)

# Don't forget to add maxtime to your for-loop
for(time in 1:maxtime){
  out_tab <- mod1(out_tab = out_tab, time_step = time)
  ## print(out_tab)
  out_tab[time + 1, c("Nh")] <- Nh[time + 1]
  out_tab[time + 1, c("Ih")] <- Ih[time + 1]

Looking at posfun, we see that it simply checks each of our state variables at each time step, and sets any negative values to zero. If we check our output we see it gives us non-negative population densities:

 out_tab
        Sv       Ev       Iv Nh Ih
1  10.0000  3.00000 2.000000  2  1
3 179.7493 16.67784 3.244288  1  0
4 416.2958 10.01499 6.133576  8  4

Great, right? Well, not quite. Unfortunately, there is currently no way to use events when method = iteration. Given what you've shared about your model so far, there are certainly arguments for modelling it in continuous time. Just because dispersal events are discretized does not necessarily mean that birth, death, and infection events must be discrete as well. It is important to conceptualize the different timescales at which these phenomena are occurring in real life, and ask yourself whether or not you are capturing them accurately with your model. However, this is already getting beyond the scope of StackOverflow, so on to the second part:

Coding discrete-time models into R with deSolve

Looking at the documentation for the iteration method in deSolve, we see that: "Method "iteration" is special in that here the function func should return the new value of the state variables rather than the rate of change."

You have clearly coded in a continuous-time model that returns the value of derivatives, and not a discrete-time model that returns the value of your state variables. Your birth component in dSv is a bare constant, so you are using an intrinsic birth rate; in a discrete-time model your birth component would be some constant number of offspring (almost always an integer) multiplied by the number of "parents" from the previous timestep.

Looking at your system of equations, we can see how this creates a problem: Instead of having each individual of Sv give birth to 100 individuals at each time step, you are setting Sv to 100. It is inevitable that Sv will plummet down past zero as you quickly accumulate a net loss.

An example discrete-time model would look something like:

# discrete-time host-parasite model
Parasite <- function(t, y, ks) {
P <- y[1]
H <- y[2]
f <- A * P / (ks + H)
Pnew <- H * (1 - exp(-f))
Hnew <- H * exp(rH * (1 - H) - f)
list (c(Pnew, Hnew))
}

Notice how we constantly refer to the "present" values of P and H in order to calculate their population dynamics for the next time step. I hope this has been helpful to you, and I wish you luck with your modelling adventures!