In Statistical Inference in a Stochastic Epidemic SEIR Model with Control Intervention, a more complex model than the one we’ve seen yesterday was considered (and is called the SEIR model). Consider a population of size \(N\), and assume that \(S\) is the number of susceptible, \(E\) the number of exposed, \(I\) the number of infectious, and \(R\) for the number recovered (or immune) individuals, \(\displaystyle{\begin{aligned}{\frac {dS}{dt}}&=-\beta {\frac {I}{N}}S\\[8pt]{\frac {dE}{dt}}&=\beta {\frac {I}{N}}S-aE\\[8pt]{\frac {dI}{dt}}&=aE-b I\\[8pt]{\frac {dR}{dt}}&=b I\end{aligned}}\)Between \(S\) and \(I\), the transition rate is \(\beta I\), where \(\beta\) is the average number of contacts per person per time, multiplied by the probability of disease transmission in a contact between a susceptible and an infectious subject. Between \(I\) and \(R\), the transition rate is \(b\) (simply the rate of recovered or dead, that is, number of recovered or dead during a period of time divided by the total number of infected on that same period of time). And finally, the incubation period is a random variable with exponential distribution with parameter \(a\), so that the average incubation period is \(a^{-1}\).

Probably more interesting, Understanding the dynamics of ebola epidemics suggested a more complex model, with susceptible people \(S\), exposed \(E\), Infectious, but either in community \(I\), or in hospitals \(H\), some people who died \(F\) and finally those who either recover or are buried and therefore are no longer susceptible \(R\).

Thus, the following dynamic model is considered\(\displaystyle{\begin{aligned}{\frac {dS}{dt}}&=-(\beta_II+\beta_HH+\beta_FF)\frac{S}{N}\\[8pt]\frac {dE}{dt}&=(\beta_II+\beta_HH+\beta_FF)\frac{S}{N}-\alpha E\\[8pt]\frac {dI}{dt}&=\alpha E+\theta\gamma_H I-(1-\theta)(1-\delta)\gamma_RI-(1-\theta)\delta\gamma_FI\\[8pt]\frac {dH}{dt}&=\theta\gamma_HI-\delta\lambda_FH-(1-\delta)\lambda_RH\\[8pt]\frac {dF}{dt}&=(1-\theta)(1-\delta)\gamma_RI+\delta\lambda_FH-\nu F\\[8pt]\frac {dR}{dt}&=(1-\theta)(1-\delta)\gamma_RI+(1-\delta)\lambda_FH+\nu F\end{aligned}}\)In that model, parameters are \(\alpha^{-1}\) is the (average) incubation period (7 days), \(\gamma_H^{-1}\) the onset to hospitalization (5 days), \(\gamma_F^{-1}\) the onset to death (9 days), \(\gamma_R^{-1}\) the onset to “recovery” (10 days), \(\lambda_F^{-1}\) the hospitalisation to death (4 days) while \(\lambda_R^{-1}\) is the hospitalisation to recovery (5 days), \(\eta^{-1}\) is the death to burial (2 days). Here, numbers are from Understanding the dynamics of ebola epidemics (in the context of ebola). The other parameters are \(\beta_I\) the transmission rate in community (0.588), \(\beta_H\) the transmission rate in hospital (0.794) and \(\beta_F\) the transmission rate at funeral (7.653). Thus

epsilon = 0.001 
Z = c(S = 1-epsilon, E = epsilon, I=0,H=0,F=0,R=0)
p=c(alpha=1/7*7, theta=0.81, delta=0.81, betai=0.588,
    betah=0.794, blambdaf=7.653,N=1, gammah=1/5*7,
    gammaf=1/9.6*7, gammar=1/10*7, lambdaf=1/4.6*7,
    lambdar=1/5*7, nu=1/2*7)

If \(\boldsymbol{Z}=(S,E,I,H,F,R)\), if we write \(\frac{\partial \boldsymbol{Z}}{\partial t} = SEIHFR(\boldsymbol{Z})\)where \(SEIHFR\) is

SEIHFR = function(t,Z,p){
  S=Z[1]; E=Z[2]; I=Z[3]; H=Z[4]; F=Z[5]; R=Z[6]
  alpha=p["alpha"]; theta=p["theta"]; delta=p["delta"]
  betai=p["betai"]; betah=p["betah"]; gammah=p["gammah"]
  gammaf=p["gammaf"]; gammar=p["gammar"]; lambdaf=p["lambdaf"]
  lambdar=p["lambdar"]; nu=p["nu"]; blambdaf=p["blambdaf"]
  N=S+E+I+H+F+R
  dS=-(betai*I+betah*H+blambdaf*F)*S/N
  dE=(betai*I+betah*H+blambdaf*F)*S/N-alpha*E
  dI=alpha*E-theta*gammah*I-(1-theta)*(1-delta)*gammar*I-(1-theta)*delta*gammaf*I
  dH=theta*gammah*I-delta*lambdaf*H-(1-delta)*lambdaf*H
  dF=(1-theta)*(1-delta)*gammar*I+delta*lambdaf*H-nu*F
  dR=(1-theta)*(1-delta)*gammar*I+(1-delta)*lambdar*H+nu*F
  dZ=c(dS,dE,dI,dH,dF,dR)
  list(dZ)}

We can solve it, or at least study the dynamics from some starting values

library(deSolve)
times = seq(0, 50, by = .1)
resol = ode(y=Z, times=times, func=SEIHFR, parms=p)

For instance, the proportion of people infected is the following

plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="",col="red")
lines(resol[,"time"],resol[,"H"],col="blue")