SIR/SEIR モデルの初期の指数増加

expgrowth_sir.R

#reference :
#Junling Ma.(2020).Estimating epidemic exponential growth rate and basic reproduction number,Infectious Disease Modelling,5,129-141.
#https://doi.org/10.1016/j.idm.2019.12.009.

library(deSolve)
library(ggplot2)
library(tidyr)
library(dplyr)
library(patchwork)

SIRmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dS <- - beta*S*I
    dI <- beta*S*I - gamma*I
    dR <- gamma*I 
    return(list(c(dS, dI, dR)))
  })
}

pars_sir  <- c(beta=1, gamma=0.25)
ini_sir = c(S=0.999, I=0.001, R=0)
times <- seq(0, 50, by = 0.1)
ode_out <- ode(y=ini_sir, times=times, func=SIRmod, parms=pars_sir)

sir_out <- data.frame(ode_out) %>% 
  pivot_longer(S:R) %>% 
  mutate(I = name=="I")

lambda_sir <- (pars_sir["beta"]-pars_sir["gamma"])

p1 <- ggplot(sir_out, aes(x=time, y=value))+
  geom_line(aes(group=name,colour=I), show.legend = FALSE)+
  scale_color_manual(values = c("lightgrey","black"))+
  stat_function(fun = function(x)ini_sir[2]*exp(lambda_sir*x),
                linetype=2, colour="royalblue", n=1001, linewidth=1)+
  theme_bw()+labs(title="SIR") +
  ylim(c(0,1.05))

####

SEIRmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dS <- - beta*S*I
    dE <- beta*S*I - sigma*E
    dI <- sigma*E - gamma*I
    dR <- gamma*I 
    return(list(c(dS, dE, dI, dR)))
  })
}

pars_seir  <- c(beta=1, sigma=1, gamma=0.25, N=100)
ini_seir = c(S=0.999, E=0, I=0.001, R=0)
times <- seq(0, 50, by = 0.1)
ode_out <- ode(y=ini, times=times, func=SEIRmod, parms=pars)

seir_out <- data.frame(ode_out) %>% 
  pivot_longer(S:R) %>% 
  mutate(I = name=="I")

lambda <- 0.5*(-(pars_seir["sigma"]+pars_seir["gamma"]) +
  sqrt((pars_seir["sigma"]-pars_seir["gamma"])^2+4*(pars_seir["sigma"]*pars_seir["beta"])))

p2 <- ggplot(seir_out, aes(x=time, y=value))+
  geom_line(aes(group=name,colour=I), show.legend = FALSE)+
  scale_color_manual(values = c("lightgrey","black"))+
  stat_function(fun = function(x)ini_seir[3]*exp(lambda*x),
                linetype=2, colour="royalblue", n=1001, linewidth=1)+
  theme_bw()+labs(title="SEIR") +
  ylim(c(0,1.05))

p1+p2
ggsave("exp_growth.png",width = 7,height = 5)