deSolve delay differential equation model of malaria not producing expected equilibrium

Question

I am trying to reproduce the following delay differential equation model of malaria in R: dde model

However, while I would expect the state values to reach an equilibrium, they just seem to immediately plummet... My code is the following:

model1 <- function(t, y, parms) {

  lagIh <- ifelse(t < taum, y[1], lagvalue(t - taum, 1))
  lagEm <- ifelse(t < taum, y[2], lagvalue(t - taum, 2))
  lagIm <- ifelse(t < taum, y[3], lagvalue(t - taum, 3))
      
  dIh <- (a*b*m*y[3])*(1 - y[1])-(r*y[1])
  dEm <- (a*c*y[1])*(1 - y[2] - y[3])-(a*c*lagIh)*(1 - lagEm - lagIm)*exp(-mu2*taum)-(mu2*y[2])
  dIm <- (a*c*lagIh)*(1 - lagEm - lagIm)*(exp(-mu2*taum))-(mu2*y[3])
    
    list(c(dIh, dEm, dIm))
  }

times <- seq(-15, 10000, by = 1)

#states
y <- c(Ih = 0.1,
       Em = 0.1,
       Im = 0.1)

#parameters
parms <- c(a = 0.5,
           b = 0.5,
           c = 0.5,
           m = 20,
           r = 0.05,
           mu2 = 0.05,
           taum = 5)

library(deSolve)
#solve
output <- dede(y = y, times = times, func = model1, parms = parms)

plot(output, xlab = "time", ylab = "-")

I have tried tinkering with the parameters and initial state values but nothing seems to change it. Any help would be greatly appreciated

Answer 1

The model missed a way to access its parameters in parms. After fixing this with a with()-construction, it runs through and reaches an equilibrium.

With a shorter simulation time (times times <- seq(-15, 200, by = 1)), it is easier to see how the equilibrium is reached.

library(deSolve)

model1 <- function(t, y, parms) {
  
  ## add with()-construction to access the parms by their name
  with(as.list(parms), {
    lagIh <- ifelse(t < taum, y[1], lagvalue(t - taum, 1))
    lagEm <- ifelse(t < taum, y[2], lagvalue(t - taum, 2))
    lagIm <- ifelse(t < taum, y[3], lagvalue(t - taum, 3))
    
    dIh <- (a*b*m*y[3])*(1 - y[1])-(r*y[1])
    dEm <- (a*c*y[1])*(1 - y[2] - y[3])-(a*c*lagIh) * 
             (1 - lagEm - lagIm)*exp(-mu2*taum)-(mu2*y[2])
    dIm <- (a*c*lagIh)*(1 - lagEm - lagIm)*(exp(-mu2*taum))-(mu2*y[3])
    
    list(c(dIh, dEm, dIm))
  }) # end with
}

#times <- seq(-15, 10000, by = 1)

## set it to shorter end time, to see the equilibrium better
times <- seq(-15, 150, by = 1)

#states
y <- c(Ih = 0.1,
       Em = 0.1,
       Im = 0.1)

#parameters
parms <- c(a = 0.5,
           b = 0.5,
           c = 0.5,
           m = 20,
           r = 0.05,
           mu2 = 0.05,
           taum = 5)


#solve
output <- dede(y = y, times = times, func = model1, parms = parms)

plot(output, xlab = "time", ylab = "-", mfrow=c(1, 3))