MSEIR
Enwu Liu
2023-05-05
This is the note for the MSEIR model (Models with maternal antibodies and latent periods) from the book “Modeling infectious disease parameters based on serological and social contact data: a modern statistical perspectiv” in page 49, the authors showed the ODE system and the solution of the ODE by the figure 3.18. The ODE system can be sovled by R deSolve pacakge, however, the authors sloved the system manually (showed here by the code they provided), the manual solution is quite tedious, here I show how to get the manual solution step by step.
Here is the code and results:
### MSEIR Model
epi11<- function(last.age,d,lambda,sigma,ni)
{
#d <- (1/0.5)
#lambda <- (1/5)
#sigma <- (34.76)
#ni <- (31.74)
N0 <- 1000
a <- seq(from=1,to=last.age,length=500)
la <- 1
N1a <- N0*la
ia <- exp(-d*a)
kk1 <- (d/(d-lambda))
kk2 <- (exp(-lambda*a)-exp(-d*a))
xa <- (d/(d-lambda))*(exp(-lambda*a)-exp(-d*a))
#browser()
ha <- (
(lambda*d)/(d-lambda)
)*
(
((exp(-sigma*a)-exp(-lambda*a))/(lambda-sigma))
-((exp(-sigma*a)-exp(-d*a))/(d-sigma))
)
ya <- (sigma*lambda*d)*
(
((exp(-ni*a)-exp(-sigma*a))/((lambda-sigma)*(d-sigma)*(sigma-ni)))
+((exp(-ni*a)-exp(-lambda*a))/((lambda-d)*(lambda-sigma)*(lambda-ni)))
+((exp(-ni*a)-exp(-d*a))/((d-lambda)*(d-sigma)*(d-ni)))
)
a <- c(0,a)
ia <- c(1,ia)
xa <- c(0,xa)
ha <- c(0,ha)
ya <- c(0,ya)
za <- 1 - ia - xa - ha -ya
plot(a,ia,type="l",xlab="Age",ylab="M(a)",pch=0.5,xlim=c(0,last.age))
#title("a:Proportion of host with maternal antibodies",adj=0,cex=0.35)
plot(a,xa,type="l",xlab="Age",ylab="S(a)",pch=0.5,xlim=c(0,last.age))
#title("b:Proportion of susceptibles",adj=0,cex=0.35)
plot(a,ha,type="l",xlab="Age",ylab="E(a)",pch=0.5,xlim=c(0,last.age))
#title("c:Proportion of host in the latent class",adj=0,cex=0.35)
plot(a,ya,type="l",xlab="Age",ylab="I(a)",pch=0.5,xlim=c(0,last.age))
#title("d:Proportion of infected",adj=0,cex=0.35)
plot(a,za,type="l",xlab="Age",ylab="R(a)",pch=0.5,xlim=c(0,last.age))
#title("e:Proportion host in the immune class",adj=0,cex=0.35)
plot(a,(ia+za+ya),type="l",xlab="Age",ylab="seroprevalence",pch=0.5,xlim=c(0,last.age))
#title("f:Proportion of sero-positive",adj=0,cex=0.35)
}
### FIGURE 3.18
par(mfrow=c(3,2),lwd=2,las=1,cex.axis=1.1,cex.lab=1.1,mgp=c(3, 0.5, 0))
epi11(40,1/0.5,0.2,26.07,36.5)
Note the above R code did not solve the ODE system of the Equation 3.12 in the book. In fact,they solved a simpler ODE system without considering the nature mortality rate, i.e the following system
\[ \left\{ \begin{array}{l} \frac{dM(t)}{dt}=-\gamma M(t) & (1) \\ \frac{dS(t)}{dt}=\gamma M(t)-\lambda S(t) &(2)\\ \frac{dE(t)}{dt}=\lambda S(t)-\sigma E(t) &(3)\\ \frac{dI(t)}{dt}=\sigma E(t)-\nu I(t) &(4)\\ \frac{dR(t)}{dt}=\nu I(t) &(5) \end{array} \right. \] Note, I use \(t\) to replace \(a\), \(t\) stands for age variable here. Now let us solve the system manually.
For equation \((1)\), we can solve it by separation of variables:
\[\frac{dM(t)}{M(t)}\frac{1}{dt}=-\gamma\] Integrate both sides in terms of \(dt\) we get:
\[\int\frac{dM(t)}{M(t)}\frac{1}{dt} dt=-\int\gamma dt \] \[\Rightarrow ln M(t)=-\gamma t+c\] During the the time period from 0 to \(a\) we have \[\Rightarrow M(a)=e^{-\gamma a} \tag{6}\]
a<- seq(from=1,to=40,length=500)
gamma<-2 #this is d in the R code
Ma<-exp(-gamma*a)
a <- c(0,a)
Ma <- c(1,Ma)
plot(a,Ma,type="l",xlab="Age",ylab="M(a)",pch=0.5,xlim=c(0,40))
Next, let us solve the second equation,
\[\begin{align*} &\frac{S(t)}{dt}=\gamma M(t)-\lambda S(t)\\ &\Leftrightarrow \frac{S(t)}{dt}+\lambda S(t)=\gamma M(t)\\ &\Leftrightarrow \frac{S(t)}{dt}+\lambda S(t)=\gamma e^{-\gamma t} \text{ replace M(t) with the result from (1)} \end{align*}\]
This is just a first order linear ODE with the following form
\[\frac{dy}{dt} + p(t)y = g(t)\] Therefore, \[\begin{align*} S(a)&=\frac{\int_0^a e^{\lambda t}\gamma e^{-\gamma t}dt}{e^{\int_0^a\lambda dt}}\\ &=\frac{\gamma\int_0^a e^{(\lambda-\gamma)t}dt}{e^{\lambda a}}\\ &=\frac{\frac{\gamma}{\lambda-\gamma}\left [ e^{(\lambda-\gamma)a}-1\right ]}{e^{\lambda a}}\\ &=\frac{\gamma}{\gamma-\lambda}(e^{-\lambda a}-e^{-\gamma a}) \tag{7} \end{align*}\]
This is the solution for the equation \((2)\), we can use the same method to solve equation \((3),(4),(5)\), next, let us plot the solution for the equation \((2)\).
a<- seq(from=1,to=40,length=500)
gamma<-2
lambda <- (1/5)
Sa<-(gamma/(gamma-lambda))*(exp(-lambda*a)-exp(-gamma*a))
#add star values
a<- c(0,a)
Sa <- c(0,Sa)
#plot
plot(a,Sa,type="l",xlab="Age",ylab="S(a)",pch=0.5,xlim=c(0,40))
To solve equation \((3)\), we replace \(S(t)\) with the solution from the equation \((2)\)
\[\frac{dE(t)}{dt}=\lambda \frac{\gamma}{\gamma-\lambda}(e^{-\lambda t}-e^{-\gamma t})-\sigma E(t)\] It is still a first oder linear differential equation, therefore the solution can be written as
\[\begin{align*}E(a)&=\frac{\int_0^a e^{\sigma t} \lambda \frac{\gamma}{\gamma-\lambda}(e^{-\lambda t}-e^{-\gamma t})dt }{e^{\int_0^a\sigma dt}}\\ &=\frac{\frac{\lambda \gamma}{\gamma-\lambda} \int_0^a \left[e^{(\sigma -\lambda)t}-e^{(\sigma-\gamma)t}\right]dt }{e^{\int_0^a\sigma dt}}\\ &=\frac{\frac{\lambda \gamma}{\gamma-\lambda}\left[\frac{1}{\sigma-\lambda}(e^{(\sigma -\lambda)a}-1)-(\frac{1}{\sigma-\gamma}e^{(\sigma-\gamma)a}-1)\right]}{e^{\sigma a}}\\ &=\frac{\lambda \gamma}{\gamma-\lambda}\left[\frac{1}{\sigma-\lambda}(e^{ -\lambda a}-e^{-\sigma a})-\frac{1}{\sigma-\gamma}(e^{-\gamma a}-e^{-\sigma a})\right]\\ &=\frac{\lambda \gamma}{\gamma-\lambda}\left[\frac{1}{\lambda-\sigma}(e^{-\sigma a}-e^{ -\lambda a})-\frac{1}{\gamma-\sigma}(e^{-\sigma a}-e^{-\gamma a})\right] \end{align*}\]
Now we can plot the solution for the equation \((3)\)
a<- seq(from=1,to=40,length=500)
gamma<-2
lambda <-1/5
sigma<-26.07
Ea<-((lambda*gamma)/(gamma-lambda))*(
((exp(-sigma*a)-exp(-lambda*a))/(lambda-sigma))
-((exp(-sigma*a)-exp(-gamma*a))/(gamma-sigma))
)
#add star values
a<- c(0,a)
Ea <- c(0,Ea)
#plot
plot(a,Ea,type="l",xlab="Age",ylab="E(a)",pch=0.5,xlim=c(0,40))
We use the results from equation \((3)\) to solve equation \((4)\)
\[\frac{dI(t)}{dt}=\sigma E(t)-\nu I(t)\\ \Leftrightarrow \frac{dI(t)}{dt}=\frac{\sigma \lambda \gamma}{\gamma-\lambda}\left[\frac{1}{\lambda-\sigma}(e^{-\sigma t}-e^{ -\lambda t})-\frac{1}{\gamma-\sigma}(e^{-\sigma t}-e^{-\gamma t})\right]-\nu I(t)\]
Still, this is a first order linear differential equation, the solution is: \[\begin{align*} I(a)&=\frac{\int_0^a e^{\nu t}\frac{\sigma \lambda \gamma}{\gamma-\lambda}\left[\frac{1}{\lambda-\sigma}(e^{-\sigma t}-e^{ -\lambda t})-\frac{1}{\gamma-\sigma}(e^{-\sigma t}-e^{-\gamma t})\right]dt}{e^{\int_0^a \nu dt}}\\ &=\frac{\frac{\sigma \lambda \gamma}{\gamma-\lambda}\int_0^a\left[\frac{1}{\lambda-\sigma}(e^{(\nu-\sigma)t}-e^{ (\nu-\lambda) t})-\frac{1}{\gamma-\sigma}(e^{(\nu-\sigma) t}-e^{(\nu-\gamma)t})\right]dt}{e^{\int_0^a \nu dt}}\\ &=\frac{\frac{\sigma \lambda \gamma}{\gamma-\lambda}\left[\frac{1}{(\lambda-\sigma)(\nu-\sigma)}(e^{(\nu-\sigma)a}-1)-\frac{1}{(\lambda-\sigma)(\nu-\lambda)}(e^{(\nu-\lambda) a}-1)-\frac{1}{(\gamma-\sigma)(\nu-\sigma)}(e^{(\nu-\sigma) a}-1)+\frac{1}{(\gamma-\sigma)(\nu-\gamma)}(e^{(\nu-\gamma)a}-1)\right]}{e^{\nu a}}\\ &=\frac{\sigma \lambda \gamma}{\gamma-\lambda}\left[\frac{1}{(\lambda-\sigma)(\nu-\sigma)}(e^{-\sigma a}-e^{-\nu a})-\frac{1}{(\lambda-\sigma)(\nu-\lambda)}(e^{-\lambda a}-e^{-\nu a})-\frac{1}{(\gamma-\sigma)(\nu-\sigma)}(e^{-\sigma a}-e^{-\nu a})+\frac{1}{(\gamma-\sigma)(\nu-\gamma)}(e^{-\gamma a}-e^{-\nu a})\right]\\ &=\frac{\sigma \lambda \gamma}{\gamma-\lambda}\left[\frac{\gamma-\sigma-\lambda+\sigma}{(\lambda-\sigma)(\nu-\sigma)(\gamma-\sigma)}(e^{-\sigma a}-e^{-\nu a})-\frac{1}{(\lambda-\sigma)(\nu-\lambda)}(e^{-\lambda a}-e^{-\nu a})+\frac{1}{(\gamma-\sigma)(\nu-\gamma)}(e^{-\gamma a}-e^{-\nu a})\right]\\ &=\frac{\sigma \lambda \gamma}{\gamma-\lambda}\left[\frac{\gamma-\lambda}{(\lambda-\sigma)(\nu-\sigma)(\gamma-\sigma)}(e^{-\sigma a}-e^{-\nu a})-\frac{1}{(\lambda-\sigma)(\nu-\lambda)}(e^{-\lambda a}-e^{-\nu a})+\frac{1}{(\gamma-\sigma)(\nu-\gamma)}(e^{-\gamma a}-e^{-\nu a})\right]\\ &=\frac{\sigma \lambda \gamma}{(\lambda-\sigma)(\nu-\sigma)(\gamma-\sigma)}(e^{-\sigma a}-e^{-\nu a})-\frac{\sigma \lambda \gamma}{(\gamma-\lambda)(\lambda-\sigma)(\nu-\lambda)}(e^{-\lambda a}-e^{-\nu a})+\frac{\sigma \lambda \gamma}{(\gamma-\lambda)(\gamma-\sigma)(\nu-\gamma)}(e^{-\gamma a}-e^{-\nu a})\\ &=\frac{\sigma \lambda \gamma}{(\lambda-\sigma)(\gamma-\sigma)(\sigma-\nu)}(e^{-\nu a}-e^{-\sigma a})+ \frac{\sigma \lambda \gamma}{(\lambda-\gamma)(\lambda-\sigma)(\lambda-\nu)}(e^{-\nu a}-e^{-\lambda a})+\frac{\sigma \lambda \gamma}{(\gamma-\lambda)(\gamma-\sigma)(\gamma-\nu)}(e^{-\nu a}-e^{-\gamma a}) \end{align*}\]
This is the solution for the equation \((4)\), next we plot the solution.
a<- seq(from=1,to=40,length=500)
gamma<-2
lambda <-1/5
sigma<-26.07
nu<-36.5
Ia<- (sigma*lambda*gamma)*
(
((exp(-nu*a)-exp(-sigma*a))/((lambda-sigma)*(gamma-sigma)*(sigma-nu)))
+((exp(-nu*a)-exp(-lambda*a))/((lambda-gamma)*(lambda-sigma)*(lambda-nu)))
+((exp(-nu*a)-exp(-gamma*a))/((gamma-lambda)*(gamma-sigma)*(gamma-nu)))
)
#add star values
a<- c(0,a)
Ia <- c(0,Ia)
#plot
plot(a,Ia,type="l",xlab="Age",ylab="I(a)",pch=0.5,xlim=c(0,40))
For the last equation, we can calculate \[R(a)=1-M(a)-S(a)-E(a)-I(a)\], here we can see we are modelling the proportions, not the actual number of cases, usually when we modelling proportions we use lower case letters in equations, however, I used capital letters here.
a<- seq(from=1,to=40,length=500)
Ra<-1-Ma-Sa-Ea-Ia
#add star values
a<- c(0,a)
#plot
plot(a,Ra,type="l",xlab="Age",ylab="R(a)",pch=0.5,xlim=c(0,40))
The last plot
a<- seq(from=1,to=40,length=500)
Sp<-Ma+Ra+Ia #sero positive
#add star values
a<- c(0,a)
#plot
plot(a,Sp,type="l",xlab="Age",ylab="seroprevalence",pch=0.5,xlim=c(0,40))
In fact, the above ODE system can be solved by the deSolve package
simultaneously. Next, I will show how to solve the MSEIR using the R
package and I will add a natural mortality rate into the system, i.e we
will solve the system below. And we assume everyone was born with
passive immunity and the birth rate is equal to natural mortality rate.
We assume the natural mortality rate is 1/75, i.e the inverse of the
life expectacny which we suppose to be 75 years old. 
\[ \left\{ \begin{array}{l} \frac{dM(t)}{dt}=N\mu-(\gamma+\mu) M(t) & (1) \\ \frac{S(t)}{dt}=\gamma M(t)-(\lambda +\mu)S(t) &(2)\\ \frac{dE(t)}{dt}=\lambda S(t)-(\sigma+\mu) E(t) &(3)\\ \frac{dI(t)}{dt}=\sigma E(t)-(\nu +\mu) I(t) &(4)\\ \frac{dR(t)}{dt}=\nu I(t)-\mu R(t) &(5) \end{array} \right. \]
SIR<-function(t,state,parameters)
{
with(as.list(c(state, parameters)),
{
N=1000
dM <-N*mu-gamma*M-mu*M
dS <- gamma*M-lambda*S-mu*S
dE <- lambda*S-sigma*E-mu*E
dI <- sigma*E - nu*I-mu*I
dR <- nu*I-mu*R
list(c(dM,dS, dE, dI,dR))
})
}
### define parameters
parameters <- c(gamma=2,lambda = 0.2, mu=1/75,nu=36.5,sigma=26.07)
state <- c(M=1000,S=0,E=1,I=0, R=0)
times<-seq(0,40,by=0.01)
require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
out$SP<-out$M+out$I+out$R
plot(out$time,out$M,type="l",xlab="Age",ylab="M(a)",xlim=c(0,40),ylim=c(0,1000))
plot(out$time,out$S,type="l",xlab="Age",ylab="S(a)",xlim=c(0,40),ylim=c(0,1000))
plot(out$time,out$E,type="l",xlab="Age",ylab="E(a)",xlim=c(0,40))
plot(out$time,out$I,type="l",xlab="Age",ylab="I(a)",xlim=c(0,40))
plot(out$time,out$R,type="l",xlab="Age",ylab="R(a)",xlim=c(0,40))
plot(out$time,out$SP,type="l",xlab="Age",ylab="Seroprevalence",xlim=c(0,40))
Reference
1.Hens, Niel, et al. Modeling infectious disease parameters based on serological and social contact data: a modern statistical perspective. Springer Science & Business Media, 2012.