setting the lower limit of a derivative in R desolve

Question

I am creating a mathematical model using R's package deSolve. I have three derivatives: dx (the proportion infected hosts), dY (the proportion infected vectors), and dm (which is the ratio of vectors to hosts in the population). The point of my model is to show the effects of a certain insecticide treatment on the population (with the effects represented as parameter "z". To incorporate this time dependent covariate into the model, the approxfun function was used. The model is working correctly, however I would like to set a lower limit for dm (assuming that not all the vectors in the population would be effected). Without setting the lower limit, my code and graph look like this:

Initial vectors for days post treatment, % killed
x <- c(4, 30, 60, 90, 120, 210, 360)
z <- c(1.0, 0.99, 0.99, 0.79, 0.7, 0.02, 0) 

plot(z ~ x)

#=============================================
#  fit data with logistic curve
#       -extract fit values using equation y = Asym / (1 + exp((xmid - input) / scal))
#=============================================
fit2 <- nls(z ~ SSlogis(x, Asym, xmid, scal), data = data.frame(x, z))
summary(fit2)

lines(seq(0, 400, length.out = 400),
      predict(fit2, newdata = data.frame(x = seq(0.5, 400, length.out = 400))))

Asym<-summary(fit2)$parameters[1,1]
xmid<-summary(fit2)$parameters[2,1]
scal<-summary(fit2)$parameters[3,1]

times <- seq(0, 1000, by = 1)
signal <- data.frame(times = times, import = rep(0, length(times)))
signal$import=  Asym / (1 + exp((xmid - times) / scal))

#Force time dependent covariate into the model
input <- approxfun(signal, rule = 2)


RMTx2 <- function(times, stateTx2, parametersTx2)   
{
  with(
    as.list(c(stateTx2, parametersTx2)), 
    {
      z <- input(times)
      dX <- ((m*a*b*Y)+(p*k*(a*m*z*Y)))*(1-X)-r*X 
      dY <- a*c*X*(exp(-g*n)-Y)-((g*(1-m/K)*Y)+(m*a*z*Y)) 
      dm <- ((R*(1-m/K)*m )+(-m*a*z))
      return(list(c(dX, dY, dm)))
    }
  )
}


initTx2 <- c(X = 0.01, Y= 0, m=40) 
parametersTx2 <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09)
outTx2 <- as.data.frame(ode(y = initTx2, times = times, func = RMTx2, parms = parametersTx2))
RESULTS2<-data.frame(outTx2$X,outTx2$Y)
RESULTS2m <-data.frame(outTx2$m, outTx2$Y*outTx2$Y)
timesTx2 <- seq(0, 10000, by = 1)

And here is the plot. What I would like to do is limit the drop in dm over time so that it will not drop below a certain value during treatment

I am trying to set the parameter dm so that, for example, the value cannot drop below 15. I have tried a few codes to do this including:

RMTx2 <- function(times, stateTx2, parametersTx2)   
{
  with(
    as.list(c(stateTx2, parametersTx2)), 
    {
      z <- input(times)
      dX <- ((m*a*b*Y)+(p*k*(a*m*z*Y)))*(1-X)-r*X 
      dY <- a*c*X*(exp(-g*n)-Y)-((g*(1-m/K)*Y)+(m*a*z*Y))
      dm <- if (isTRUE (((R*(1-m/K)*m )+(-m*a*z)) > MM)) ((R*(1-m/K)*m )+(-m*a*z)) else 15 
      return(list(c(dX, dY, dm)))
    }
  )
}

initTx2 <- c(X = 0.01, Y= 0, m=40) 
parametersTx2 <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09, MM= 15)
outTx2 <- as.data.frame(ode(y = initTx2, times = times, func = RMTx2, parms = parametersTx2))
RESULTS2<-data.frame(outTx2$X,outTx2$Y)
RESULTS2m <-data.frame(outTx2$m, outTx2$Y*outTx2$Y)

Unfortunately, for some reason this is causing the population to just increase indefinitely:

Is there something fundamental about this approach that won't work? Or is this more a coding error? Thanks!

Answer 1

You may try a Monod-type safeguard for this. MM is your threshold and km the Monod-constant. Its value is quite uncritical, but it should be not too small to avoid numerical difficulties for the solver.

RMTx2 <- function(times, stateTx2, parametersTx2)   
{
  with(
    as.list(c(stateTx2, parametersTx2)), 
    {
      z <- input(times)
      dX <- ((m*a*b*Y)+(p*k*(a*m*z*Y)))*(1-X)-r*X 
      dY <- a*c*X*(exp(-g*n)-Y)-((g*(1-m/K)*Y)+(m*a*z*Y)) 

      #dm <- (R*(1-m/K)*m )+(-m*a*z)
      dm <- (R*(1-m/K)*m )+(-m*a*z) * (m - MM) / (km + (m - MM))  

      list(c(dX, dY, dm), dm=dm) ## return derivative dm as additional output
    }
  )
}


initTx2 <- c(X = 0.01, Y= 0, m=40) 

parameters_old <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09, MM = 0, km=0)
parameters_new <- c(a=1/14, b=0.00068, n=45, g= 0.005, c=0.28, k= 0.10, r= 1/(3*365), p=0, K=40, R= 0.09, MM = 15, km=0.1)

out_old <- ode(y = initTx2, times = times, func = RMTx2, parms = parameters_old)
out_new <- ode(y = initTx2, times = times, func = RMTx2, parms = parameters_new)

plot(out_old, out_new)