Using dede in deSolve package in RStudio to solve time-delayed ODE

Question

I'm trying to solve this SEIRVB ODE with time delay using dede in deSolve package of R in RStudio

I'm currently have this code

library(deSolve)

tau=0


# Parameters

Lambda = 0.5      
beta1 = 1.13921549    
beta2 = 2.68228986      
beta3 = 2.9627036      
beta4 = 1.19686956     
beta5 = 0.5    
beta6 = 1.34496108     
beta7 = 3.40332936   
beta8 = 1.48249240    
beta9 = 0.04681161     
beta10 = 2.32555864    
alpha = 0.5    
kappa = 0.02933240    
d = 0.0047876         
r = 0.09871        
tau=7 #time delay


# Function representing the SEIRVB model with time delay
seirvb_model <- function(t, y, parms) {
  if (t < tau)
    Ilag <- initial_conditions["I"]
  else
    Ilag <- lagvalue(t-tau,3)
  
  

  dS <- Lambda - beta1 * y[1] * y[2] - (beta3 + alpha) * y[1]
  dE <- beta4 * y[5] * y[2] + beta6 * y[6] * y[2] + beta1 * y[1] * y[2] - beta2 * y[2] * y[3] - (kappa + alpha) * y[2]
  dI <- beta5 * y[5] * Ilag + beta7 * y[6] * y[3] + beta2 * y[2] * y[3] - (d + alpha + r) * y[3]
  dR <- beta9 * y[5] + beta8 * y[6] + r * y[3] + kappa * y[2] - alpha * y[4]
  dV <- beta3 * y[1] - (beta9 + alpha + beta10) * y[5] - beta4 * y[5] * y[2] - beta5 * y[5] * Ilag
  dB <- beta10 * y[5] - beta6 * y[6] * y[2] - beta7 * y[6] * y[3] - (beta8 + alpha) * y[6]
  list(c(dS, dE, dI, dR, dV, dB))
}
  




# Initial conditions
initial_conditions <- c(
  S=0.3,
  E=0.1,
  I=0.006,
  R=0,
  V=0,
  B=0
)

# Time vector
times <- seq(0, 120, by = 1)

# Solve the SEIRVB model
ode_output <- dede(y = initial_conditions, times = times, func = seirvb_model, parms = NULL)

# Plotting the results
plot(ode_output, xlab = "Time", ylab = "Population", main = "SEIRVB Model for COVID-19")

Even though it works, the time delay was not considered in the solution as I got the same answer even though I changed the time delay already. Can someone help me how to fix this? Thank you!

Answer 1

It appears there is a difference:

tau <- 0 
ode_output2 <- dede(y = initial_conditions, times = times, func = seirvb_model, parms = NULL)
tau <- 40
ode_output3 <- dede(y = initial_conditions, times = times, func = seirvb_model, parms = NULL)

We can use all.equal() to see that the outputs are not the same.

Visualizing the differences between the baseline tau=7 and two alternatives (tau=0, tau=40);

png("dede_diff.png")
par(las = 1, bty = "l")
matplot(ode_output[,"I"] - cbind(ode_output3[,"I"],ode_output2[,"I"]),
        lty = 1:2, col = 1:2, ylab = "diff", type = "l")
legend("topright",
       legend = c("tau = 40", "tau=0"),
       lty = 1:2, col = 1:2)
dev.off()

You may wonder why the delay makes so little difference to the output, but it does make some difference.