Tag Archives: covid-19

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

Modeling pandemics (2)

When introducing the SIR model, in our initial post, we got an ordinary differential equation, but we did not really discuss stability, and periodicity. It has to do with the Jacobian matrix of the system. But first of all, we had three equations for three function, but actually\displaystyle{{\frac{dS}{dt}}+{\frac {dI}{dt}}+{\frac {dR}{dt}}=0}so it means that our problem is here simply in dimension 2. Hence\displaystyle {\begin{aligned}&X={\frac {dS}{dt}}=\mu(N-S)-{\frac {\beta IS}{N}},\\[6pt]&Y={\frac {dI}{dt}}={\frac {\beta IS}{N}}-(\mu+\gamma)I\end{aligned}}and therefore, the Jacobian of the system is\begin{pmatrix}\displaystyle{\frac{\partial X}{\partial S}}&\displaystyle{\frac{\partial X}{\partial I}}\\[9pt]\displaystyle{\frac{\partial Y}{\partial S}}&\displaystyle{\frac{\partial Y}{\partial I}}\end{pmatrix}=\begin{pmatrix}\displaystyle{-\mu-\beta\frac{I}{N}}&\displaystyle{-\beta\frac{S}{N}}\\[9pt]\displaystyle{\beta\frac{I}{N}}&\displaystyle{\beta\frac{S}{N}-(\mu+\gamma)}\end{pmatrix}We should evaluate the Jacobian at the equilibrium, i.e. S^\star=\frac{\gamma+\mu}{\beta}=\frac{1}{R_0}andI^\star=\frac{\mu(R_0-1)}{\beta}We should then look at eigenvalues of the matrix.

Our very last example was

times = seq(0, 100, by=.1)
p = c(mu = 1/100, N = 1, beta = 50, gamma = 10)
start_SIR = c(S=0.19, I=0.01, R = 0.8)
resol = ode(y=start_SIR, t=times, func=SIR, p=p)
plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="")

We can compute values at the equilibrium

mu=p["mu"]; beta=p["beta"]; gamma=p["gamma"]
N=1
S = (gamma + mu)/beta
I = mu * (beta/(gamma + mu) - 1)/beta

and the Jacobian matrix

J=matrix(c(-(mu + beta * I/N),-(beta * S/N),
         beta * I/N,beta * S/N - (mu + gamma)),2,2,byrow = TRUE)

Now, if we look at the eigenvalues,

eigen(J)$values
[1] -0.024975+0.6318831i -0.024975-0.6318831i

or more precisely 2\pi/b where a\pm ib are the conjuguate eigenvalues

2 * pi/(Im(eigen(J)$values[1]))
[1] 9.943588

we have a damping period of 10 time lengths (10 days, or 10 weeks), which is more or less what we’ve seen above,

The graph above was obtained using

p = c(mu = 1/100, N = 1, beta = 50, gamma = 10)
start_SIR = c(S=0.19, I=0.01, R = 0.8)
resol = ode(y=start_SIR, t=times, func=SIR, p=p)
plot(resol[1:1e5,"time"],resol[1:1e5,"I"],type="l",xlab="time",ylab="",lwd=3,col="red")
yi=resol[,"I"]
dyi=diff(yi)
i=which((dyi[2:length(dyi)]*dyi[1:(length(dyi)-1)])<0)
t=resol[i,"time"]
arrows(t[2],.008,t[4],.008,length=.1,code=3)

If we look carefully. at the begining, the duration is (much) longer than 10 (about 13)… but it does converge towards 9.94

plot(diff(t[seq(2,40,by=2)]),type="b")
abline(h=2 * pi/(Im(eigen(J)$values[1]))

So here, theoretically, every 10 weeks (assuming that our time length is a week), we should observe an outbreak, smaller than the previous one. In practice, initially it is every 13 or 12 weeks, but the time to wait between outbreaks decreases (until it reaches 10 weeks).

Modeling pandemics (1)

The most popular model to model epidemics is the so-called SIR model – or Kermack-McKendrick. Consider a population of size N, and assume that S is the number of susceptible, I the number of infectious, and R for the number recovered (or immune) individuals, \displaystyle {\begin{aligned}&{\frac {dS}{dt}}=-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-\gamma I,\\[6pt]&{\frac {dR}{dt}}=\gamma I,\end{aligned}}so that \displaystyle{{\frac{dS}{dt}}+{\frac {dI}{dt}}+{\frac {dR}{dt}}=0}which implies that S+I+R=N. In order to be more realistic, consider some (constant) birth rate \mu, so that the model becomes\displaystyle {\begin{aligned}&{\frac {dS}{dt}}=\mu(N-S)-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-(\gamma+\mu) I,\\[6pt]&{\frac {dR}{dt}}=\gamma I-\mu R,\end{aligned}}Note, in this model, that people get sick (infected) but they do not die, they recover. So here, we can model chickenpox, for instance, not SARS.

The dynamics of the infectious class depends on the following ratio:\displaystyle{R_{0}={\frac {\beta }{\gamma +\mu}}} which is the so-called basic reproduction number (or reproductive ratio). The effective reproductive ratio is R_0S/N, and the turnover of the epidemic happens exactly when R_0S/N=1, or when the fraction of remaining susceptibles is R_0^{-1}. As shown in Directly transmitted infectious diseases:Control by vaccination, if S/N<R_0^{-1} the disease (the number of people infected) will start to decrease.

Want to see it  ? Start with

mu = 0
beta = 2
gamma = 1/2

for the parameters. Here,  R_0=4. We also need starting values

epsilon = .001
N = 1
S = 1-epsilon
I = epsilon
R = 0

Then use the ordinary differential equation solver, in R. The idea is to say that \boldsymbol{Z}=(S,I,R) and we have the gradient \frac{\partial \boldsymbol{Z}}{\partial t} = SIR(\boldsymbol{Z})where SIR is function of the various parameters. Hence, set

p = c(mu = 0, N = 1, beta = 2, gamma = 1/2)
start_SIR = c(S = 1-epsilon, I = epsilon, R = 0)

The we must define the time, and the function that returns the gradient,

times = seq(0, 10, by = .1)
SIR = function(t,Z,p){
S=Z[1]; I=Z[2]; R=Z[3]; N=S+I+R
mu=p["mu"]; beta=p["beta"]; gamma=p["gamma"]
dS=mu*(N-S)-beta*S*I/N
dI=beta*S*I/N-(mu+gamma)*I
dR=gamma*I-mu*R
dZ=c(dS,dI,dR)
return(list(dZ))}

To solve this problem use

library(deSolve)
resol = ode(y=start_SIR, times=times, func=SIR, parms=p)

We can visualize the dynamics below

par(mfrow=c(1,2))
t=resol[,"time"]
plot(t,resol[,"S"],type="l",xlab="time",ylab="")
lines(t,resol[,"I"],col="red")
lines(t,resol[,"R"],col="blue")
plot(t,t*0+1,type="l",xlab="time",ylab="",ylim=0:1)
polygon(c(t,rev(t)),c(resol[,"R"],rep(0,nrow(resol))),col="blue")
polygon(c(t,rev(t)),c(resol[,"R"]+resol[,"I"],rev(resol[,"R"])),col="red")

We can actually also visualize the effective reproductive number is R_0S/N, where

R0=p["beta"]/(p["gamma"]+p["mu"])

The effective reproductive number is on the left, and as we mentioned above, when we reach 1, we actually reach the maximum of the infected,

plot(t,resol[,"S"]*R0,type="l",xlab="time",ylab="")
abline(h=1,lty=2,col="red")
abline(v=max(t[resol[,"S"]*R0&gt;=1]),col="darkgreen")
points(max(t[resol[,"S"]*R0&gt;=1]),1,pch=19)
plot(t,resol[,"S"],type="l",xlab="time",ylab="",col="grey")
lines(t,resol[,"I"],col="red",lwd=3)
lines(t,resol[,"R"],col="light blue")
abline(v=max(t[resol[,"S"]*R0&gt;=1]),col="darkgreen")
points(max(t[resol[,"S"]*R0&gt;=1]),max(resol[,"I"]),pch=19)

And when adding a \mu parameter, we can obtain some interesting dynamics on the number of infected,

times = seq(0, 100, by=.1)
p = c(mu = 1/100, N = 1, beta = 50, gamma = 10)
start_SIR = c(S=0.19, I=0.01, R = 0.8)
resol = ode(y=start_SIR, t=times, func=SIR, p=p)
plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="")