Replacing negative values in a model (system of ODEs) with zero

Question

I'm currently working on solving a system of ordinary differential equations using deSolve, and was wondering if there's any way of preventing differential variable values from going below zero. I've seen a few other posts about setting negative values to zero in a vector, data frame, etc., but since this is a biological model (and it doesn't make sense for a T cell count to go negative), I need to stop it from happening to begin with so these values don't skew the results, not just replace the negatives in the final output.

Answer 1

My standard approach is to transform the state variables to an unconstrained scale. The most obvious/standard way to do this for positive variables is to write down equations for the dynamics of log(x) rather than of x.

For example, with the Susceptible-Infected-Recovered (SIR) model for infectious disease epidemics, where the equations are dS/dt = -beta*S*I; dI/dt = beta*S*I-gamma*I; dR/dt = gamma*I we would naively write the gradient function as

gfun <- function(time, y, params) {
   g <- with(as.list(c(y,params)),
       c(-beta*S*I,
          beta*S*I-gamma*I,
          gamma*I)
       )
   return(list(g))
}

If we make log(I) rather than I be the state variable (in principle we could do this with S as well, but in practice S is much less likely to approach the boundary), then we have d(log(I))/dt = (dI/dt)/I = beta*S-gamma; the rest of the equations need to use exp(logI) to refer to I. So:

gfun_log <- function(time, y, params) {
   g <- with(as.list(c(y,params)),
       c(-beta*S*exp(logI),
          beta*S-gamma,
          gamma*exp(logI))
       )
   return(list(g))
}

(it would be slightly more efficient to compute exp(logI) once and store/re-use it rather than computing it twice ...)

Answer 2

If a value doesn’t become negative in reality but becomes negative in your model, you should change your model or, equivalently, modify your differential equations such that this is not possible. With other words: Do not try to constrain your dynamical variables but their derivatives. Everything else will only lead to problems with your solver, while it should not care about a change in the differential equation.

For a simple example, suppose that:

• you have a one-dimensional differential equation ẏ = f(y),
• y shall not become negative,
• your initial y is positive.

In this case, y can only become negative if f(0) < 0. Thus, all you have to do is to modify f such that f(0) ≥ 0 (and it is still smooth).

For a proof of principle, you can multiply f with an appropriately modified sigmoid function (which allows you to compose every logical operation with smooth functions). This way, nothing would change for most values of y, and you only change your differential equation if y is close to 0, i.e., when you were going to manipulate things anyway.

However, I would not really recommend using sigmoids without thinking about your model. If your model is totally wrong near y = 0, it will very likely already be useless for nearby values. If your simulations venture in this terrain and you want the results to be meaningful, you should fix this.