Why does my deSolve model in R stop integrating when I incorporate a conditional source of mortality in my population model?

Question

I have constructed a population model in R to identify the efficacy of potential control strategies for the cockroach population in my apartment. I have introduced a term that will add an additional source of mortality to the adult cockroach population when they encounter traps I have set. I introduce this term into the model with an if-else statement, only allowing the additional source of mortality to occur once the adult population (broken into a pregnant and non-pregnant population) has crossed a threshold of 100 individuals. I defined this as a conditional source of mortality because I believe my propensity to deploy traps will depend on my ability to detect the cockroaches in the first place, which will be a function of their adult population size. However, when I write this conditional source of mortality into my series of ODE's and apply the deSolve package, the model integration stops when the condition is met.

Below is my R code for the model structure, starting conditions, parameters, and output:

library(tidyr)
library(deSolve)

#The model:
cockpop.t <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dA <- ((lambda * N) + (theta * P) - (phi * A * (1 - (A+P)/K)) - (mua * A) - (if(A+P < det) {0}
                                                                               else(muta*A)))
    dP <- ((phi * A * (1 - (A+P)/K)) - (theta * P) - (if(A+P < det) {0}
                                                      else(muta*P)))
    dE <- (40*((theta * P) - (gamma * E) - (mue * E)))
    dN <- ((gamma * E) - (lambda * N) - (mun * N))
    return(list(c(dA, dP, dE, dN)))
  })
}

times <- seq(0, 1000, by = 1)
#Times

parameters.t <- c(lambda = 0.007142857, theta = 0.03571429, phi = 0.07142857, mua = 0.004761905, gamma = 0.04, mue = 0.004, mun = 0.0007, K = 1000, muta = 0.05, det=50)
#Model parameters

initial.A <- 9
initial.P <- 1
initial.E <- 0
initial.N <- 0
initial.values <- c(A = initial.A, P = initial.P, E = initial.E, N = initial.N) 
#Initial conditions

output.t<- as.data.frame(ode(func = cockpop.t, y = initial.values, parms = parameters.t, times = times)) %>%
  pivot_longer(!time, names_to="state",values_to="population")
#Model output 

The parameter 'det' is my detection threshold for the adult population (currently set at 100 individuals). Once I define the initial population sizes and parameters, the model runs until about the 200th timestep, at which point A + P == 100, and the model fails to integrate further. If I remove the conditional source of mortality from either A or P and leave it in the other, the model successfully integrates. What is causing this error and how might I modify my code to address it?

I tried re-writing the conditional statement as an ifelse() argument, assuming something might be written incorrectly in the way it was written above, but that did not fix the issue. I also tried defining the statement differently by constructing a separate population class, T, that kept track of the combined A and P population and using T in the conditional mortality statement instead of the (A+P) term. I thought that potentially incorporating the two class sizes in the conditional statement that contributes to the change in their respective class sizes may have caused issues, but that also did not solve the issue.

Thank you for your help!

Answer 1

The sharp on/off nature of your control function makes it hard to integrate — the way your system is set up makes it cycle on and off frequently (this is a conclusion from looking at the results below, not a careful analysis). I was able to get the integration to work, sort of, by switching to method = "rk4" which doesn't try to detect and handle stiffness, but I would not guarantee that the results are mathematically sound. (They may well be good enough for your purpose.) You could try a smoother control function ...

I wrapped the machinery in a function for more convenient experimentation.

f <- function(tmax, method = "lsoda") {
    output.t <- ode(func = cockpop.t, y = initial.values,
       parms = parameters.t, times = 0:tmax, method = method) |>
        as.data.frame() |>
        tidyr::pivot_longer(-time, names_to="state",values_to="population") |>
        ggplot(aes(time, population, colour = state)) + geom_line()
}

f(600) only makes it to to about t=285, as you mentioned.

print(gg1 <- f(600, method = "rk4"))

Zoom in:

gg1 + coord_cartesian(xlim = c(250, 350), ylim = c(15, 35))

Answer 2

Here another approach that works also with lsoda. Instead of a hard limit with an if-construction, one can use a smooth onset of the "conditional mortality". Different saturation functions can be used here, e.g. a Monod-term that goes from zero to one.

The provided solution is not only numerically stable and more precise, but also faster.

library(deSolve)

## original model with if()
cockpop.t <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dA <- ((lambda * N) + (theta * P) - (phi * A * (1 - (A + P)/K)) -
             (mua * A) - (if(A + P < det) {0} else(muta * A)))
    dP <- ((phi * A * (1 - (A+P)/K)) - (theta * P) -
             (if(A + P < det) {0} else(muta * P)))
    dE <- (40*((theta * P) - (gamma * E) - (mue * E)))
    dN <- ((gamma * E) - (lambda * N) - (mun * N))
    return(list(c(dA, dP, dE, dN)))
  })
}

## modified model with smooth transition
cockpop.t2 <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {

    ## use a Monod-term, smoothly increasing from zero to one
    x <- max(0, A + P - det)
    limiter <- x / (Km + x)

    dA <- lambda * N + theta * P - phi * A * (1 - (A + P) / K) -
             mua * A -  muta * A * limiter
    dP <- (phi * A * (1 - (A + P) / K)) - theta * P -
             muta * P * limiter

    dE <- 40 * (theta * P - gamma * E - mue * E)
    dN <- gamma * E - lambda * N - mun * N
    list(c(dA, dP, dE, dN))
  })
}


times <- seq(0, 1000, by = 1)
parameters.t <- c(lambda = 0.007142857, theta = 0.03571429,
                  phi = 0.07142857, mua = 0.004761905,
                  gamma = 0.04, mue = 0.004, mun = 0.0007,
                  K = 1000, muta = 0.05, det=50, Km = 0.1)
## The Monod constant Km is uncritical. It should be small enough
## to start early after A + P >= det, but not too small to avoid
## numerical instability.


initial.values <- c(A = 0, P = 1, E = 0, N = 0)

## Ben Bolker's solution with rk4
system.time(
o1 <- ode(func = cockpop.t, y = initial.values, parms = parameters.t,
          times = times, method = "rk4")
)
#   user  system elapsed 
#   0.09    0.00    0.09 

## tpetzoldt's approach with model reformulation
system.time(
o2 <- ode(func = cockpop.t2, y = initial.values, parms = parameters.t,
          times = times, method = "lsoda")
)
#   user  system elapsed 
#   0.06    0.00    0.06 

plot(o1, o2)
legend("topleft", c("if with rk4", "smooth transition"), 
       lty = 1:2, col = 1:2)

Instead of ggplot one can also use the built-in plot function of deSolve, that allows direct comparison of scenarios.