# SEIR model in continuous and discrete mode # (C) Minato Nakazawa # SEIR is given as matrix, days as rows, compartment as cols. # dS/dt = -beta SI + gamma R (gamma is loss of immunity rate) # dE/dt = beta SI - incident E (incident is inverse of latent period) # dI/dt = incident E - recover I (recover is inverse of the disease duration) # dR/dt = recover I - gamma R beta <- 0.2; incident <- 1/7; recover <- 1/14; gamma <- 1/30 days <- 365 pop <- 100 init.E <- 1 SEIRd <- matrix(c(pop-init.E, init.E, rep(0, days*4-2)), days, 4, byrow=TRUE) SEIRc <- SEIRd/pop head(SEIRd) # Showing initial population for (i in 2:days) { # continuous model nSc <- gamma * SEIRc[i-1, 4] nEc <- beta * SEIRc[i-1,1] * SEIRc[i-1, 3] nIc <- incident * SEIRc[i-1, 2] nRc <- recover * SEIRc[i-1, 3] SEIRc[i, 1] <- SEIRc[i-1, 1] - nEc + nSc SEIRc[i, 2] <- SEIRc[i-1, 2] + nEc - nIc SEIRc[i, 3] <- SEIRc[i-1, 3] + nIc - nRc SEIRc[i, 4] <- SEIRc[i-1, 4] + nRc - nSc # discrete model nSd <- rbinom(1, SEIRd[i-1, 4], gamma) nEd <- rbinom(1, SEIRd[i-1, 1], beta*(SEIRd[i-1, 3]/pop)) nId <- rbinom(1, SEIRd[i-1, 2], incident) nRd <- rbinom(1, SEIRd[i-1, 3], recover) SEIRd[i, 1] <- SEIRd[i-1, 1] - nEd + nSd SEIRd[i, 2] <- SEIRd[i-1, 2] + nEd - nId SEIRd[i, 3] <- SEIRd[i-1, 3] + nId - nRd SEIRd[i, 4] <- SEIRd[i-1, 4] + nRd - nSd } layout(1:2) matplot(1:days, SEIRc, type="l", col=1:4, lty=1:4, lwd=1, main="Continuous SEIR") matplot(1:days, SEIRd, type="l", col=1:4, lty=1:4, lwd=1, main="Discrete SEIR") legend("topright", lty=1:4, col=1:4, lwd=1, legend=c("S","E","I","R"))