This is a note for SIR model with vaccination\(^1\). Suppose p proportion of population was vaccinated with 100% vaccine efficacy, we can write the dynamic system as the following:

\[ \left\{\begin{array}{l} \frac{ds(t)}{dt}=\mu(1-p)-(\lambda(t)+\mu)s(t) & (1)\\ \frac{d\lambda(t)}{dt}=\nu\lambda(t)(R_0s(t)-1) &(2)\\ \end{array} \right. \] Here, \(p\) is the proportion of vaccinated, \(\lambda(t)\) is the infection force function, we will not treat it as a constant here. The equation \((2)\) can be derived from the following equations

\[\begin{align*} &\frac{di(t)}{dt}=\lambda(t)s(t)-\nu i(t)\Leftrightarrow\\ &\frac{d\lambda(t)}{\tilde{\beta}dt}=\lambda(t)s(t)-\nu \frac{\lambda(t)}{\tilde{\beta}} \text{ This is because } (\lambda(t)=\tilde{\beta}i(t)\Rightarrow i(t)=\frac{\lambda(t)}{\tilde{\beta}})\\ &\Leftrightarrow \frac{d\lambda(t)}{dt}=\tilde{\beta}\lambda(t)s(t)-\nu \lambda(t)\\ &\Leftrightarrow \frac{d\lambda(t)}{dt}=R_0\nu\lambda(t)s(t)-\nu \lambda(t) \text{ this is because } R_0=\frac{\tilde{\beta}}{\nu}\\ &\Leftrightarrow \frac{d\lambda(t)}{dt}=\nu\lambda(t)(R_0s(t)-1) \text{ This is quation} (2) \end{align*}\] \(\square\)

Note \(\tilde{\beta}=\beta N\) is so called frequency dependent or effective contact rate, \(R_0\) is the basic reproduce number where \(R_0=\frac{\tilde{\beta}}{\nu}\) when ignore natural mortality, \(1/\nu\) is the mean infectious period.

R code provided by the book\(^1\) using RK4 to solve the above system.

fx1<- function(mu,lambda.t,x.t,probii)
  {
         fx.t <- mu*(1-probii) -(lambda.t+mu)*x.t
         return(fx.t)
  } 

gl1<- function(v,mu,r0,lambda.t,x.t)
  {
         gx.t <- (v)*lambda.t*(r0*x.t-1)
         return(gx.t)
  } 

rk4.1 <- function(step,nstep,L,lambda0,v,r0,x0,pp)
{
mu      <-1/L
x       <-c(x0,c(1:nstep)*0)
pi      <-c(0,c(1:nstep)*0)+pp
lambda  <-c(lambda0,c(1:nstep)*0)
timei <- c(1:(nstep+1))-1
x.star <- 1/r0
for(j in 1:nstep){
 mx1      <- fx1(mu,lambda[j],x[j],pi[j])
 ml1      <- gl1(v,mu,r0,lambda[j],x[j])
 mx2      <- fx1(mu,lambda[j]+(0.5*step)*ml1,x[j]+(0.5*step)*mx1,pi[j])
 ml2      <- gl1(v,mu,r0,lambda[j]+(0.5*step)*ml1,x[j]+(0.5*step)*mx1)
 mx3      <- fx1(mu,lambda[j]+(0.5*step)*ml2,x[j]+(0.5*step)*mx2,pi[j])
 ml3      <- gl1(v,mu,r0,lambda[j]+(0.5*step)*ml2,x[j]+(0.5*step)*mx2)
 mx4      <- fx1(mu,lambda[j]+(step)*ml3,x[j]+(step)*mx3,pi[j])
 ml4      <- gl1(v,mu,r0,lambda[j]+(step)*ml3,x[j]+(step)*mx3)
 x[j+1]     <- x[j]+(step/6)*(mx1+2*(mx2+mx3)+mx4)
 lambda[j+1]<- lambda[j]+(step/6)*(ml1+2*(ml2+ml3)+ml4)
}
#par(mar=c(5.1,4.1,4.1,2.1))
plot(timei*step,x,type="l",xlab="time",ylab="proportion susceptible")
lines(c(0,nstep)*step,c(x.star,x.star),lty=2)
mtext("time",side=1,line=3)
#par(mar=c(7.1,4.1,2.1,2.1))
plot(timei*step,lambda,type="l",xlab="time",ylab="lambda(t)")
lines(c(0,nstep)*step,c(mu*(r0-1),mu*(r0-1)),lty=2)
mtext("time",side=1,line=3)
}

### FIGURE 3.5
par(mfrow=c(2,1),lwd=3,las=1,cex.axis=1.1,cex.lab=1.1,mgp=c(4.8, 0.5, 0),mar=c(5.1,5.5,4.1,2.1))

rk4.1(0.01,20*(1/0.01),70,0.2,25,15,0.06666,0.0)

### FIGURE 3.6
par(mfrow=c(2,1),lwd=3,las=1,cex.axis=1.1,cex.lab=1.1,mgp=c(3, 0.5, 0))
rk4.1(0.01,20*(1/0.01),70,0.2,25,15,0.06666,0.4)

### FIGURE 3.7
par(mfrow=c(2,1),lwd=3,las=1,cex.axis=1.1,cex.lab=1.1,mgp=c(3, 0.5, 0))
rk4.1(0.01,20*(1/0.01),70,0.2,25,15,0.06666,0.8)

Next, to use deSolve package to solve the system

#######################################################
### Dynamic aspects of vaccination in the SIR Model ###
#######################################################

### define parameters
parameters <- c(mu=1/75,v=25,r0=15,probii=0.8)
state <- c(y1=0.06666,y2=0.2)
times<-seq(0,20,by=0.01)

### model 
AM146<-function(t,state,parameters)
{
    with(as.list(c(state, parameters)),
    {
    dy1 = mu*(1-probii) -(y2+mu)*y1
    dy2 =v*y2*(r0*y1-1) 
    list(c(dy1,dy2))
    }) 
}

library(deSolve)
out <- as.data.frame(ode(y=state,time=times,func=AM146,parms=parameters))
head(out)
##   time         y1        y2
## 1 0.00 0.06666000 0.2000000
## 2 0.01 0.06654460 0.1999517
## 3 0.02 0.06642950 0.1998171
## 4 0.03 0.06631477 0.1995966
## 5 0.04 0.06620045 0.1992906
## 6 0.05 0.06608660 0.1989000

This is for vaccination of 80% percent

plot(out$time,out$y1,type="l",xlab="time",ylab="proportion susceptible")

plot(out$time,out$y2,type="l",xlab="time",ylab="lambda(t)")

The followings are the SIR model with vaccination and a constant effective contact rate \(\beta\).

\[ \left\{\begin{array}{l} \frac{dS(t)}{dt}=N\mu(1-p)-(\beta I(t)+\mu)S(t) \\ \frac{dI(t)}{dt}=\beta S(t)I(t)-\nu I-\mu I\\ \frac{dR}{dt}=\nu I(t)-\mu R(t)+N \mu p \end{array} \right. \] Set the effective contact rate \(\beta=0.0005\), recover rate \(\nu=1\) year\(^{-1}\), i.e the infection period is one year, the life expectancy is set to 75 years old, the the natural mortality is \(\mu=\frac{1}{75}\). The following R code showed the disease of dynamic when vaccination coverage was 0%, 20% and 40% and the process reached to the equilibrium state.

SIR<-function(t,state,parameters)
{
    with(as.list(c(state, parameters)),
    {
    dX <- N*mu*(1-p)-beta*Y*X - mu*X
    dY <- beta*Y*X - v*Y - mu*Y
    dZ <- v*Y -mu*Z+N*mu*p
    list(c(dX, dY, dZ))
    }) 
}

### define parameters
parameters <- c(mu=1/75,beta=0.0005,v=1)
state <- c(X=4999,Y=1,Z=0)
times<-seq(0,800,by=0.01)
N<-5000

require(deSolve)
p<-0
outp0 <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
p<-0.2
outp02 <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
p<-0.4
outp04 <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))

### FIGURE 3.8 - CHECK
par(mfrow=c(1,2),lwd=2,las=1,cex.axis=0.8,cex.lab=0.8,mgp=c(2, 0.45, 0),
mar=c(3.2, 2.8, 3.5, 1.5))

plot(times,outp0$X,type="l",main=" ", xlab="time", ylab="Number of susceptibles")
lines(times,outp02$X,lty=2,col='red')
lines(times,outp04$X,lty=3,col='blue')
#legend(100,5000,c("p=0","p=0.2","p=0.4"),lty=c(1,4,2))

plot(outp0$X,outp0$Y,main=" ",type="l",xlab="Number of susceptibles",ylab="Number of infected")
lines(outp02$X,outp02$Y,lty=2,col='red')
lines(outp04$X,outp04$Y,lty=3,col='blue')

References

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.