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