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.
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 ...)