Modeling Pandemics (3)

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"]

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

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


3 thoughts on “Modeling Pandemics (3)”

  1. Dear Arthur,
    Thank you for sharing your work. One small code adjustment maybe needed, in the equation for dH/dt:
    I think it should be:
    according with the model equations.

  2. Dear Arthur,

    Your blog is amazing. Thank you for all your work.
    I think I have found a minor mistake in some expressions though.

    In the equation for dI/dt the term theta*gammah*I is with positive sign. I think it should be negative (the R code is ok).

    In the equation for dF/dt the term there is the term (1-theta)*(1-delta)*gammar*I. I think it should be (1-theta)*delta*gammar*I instead. The R code may be also wrong.

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.