Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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; E=Z; I=Z; H=Z; F=Z; R=Z
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") 