24-09-2026                       >eR-BioStat and STAEM-3D

Modeling infectious diseases using R - vaccination in the SIR model

Ziv Shkedy and Rashider Aloni based on Chapter 2-4 in the book Modeling infectious diseases parameters based on serological and social contact data (https://www.simid.be/modeling-infectious-disease-parameters-based-on-serological-and-social-contact-data/) 

## load/install libraries
.libPaths(c("./Rpackages",.libPaths()))
library(knitr)
library(tidyverse)
library(deSolve)
library(minpack.lm)
library(ggpubr)
library(readxl)
library(gamlss)
library(data.table)
library(grid)
library(png)
library(nlme)
library(gridExtra)
library(mvtnorm)
library(e1071)
library(lattice)
library(ggplot2)
library(dslabs)
library(NHANES)
library(plyr)
library(dplyr)
library(nasaweather)
library(ggplot2)
library(gganimate)
library(av)
library(gifski)
library(foreach)
library("DAAG")
library(DT)
library(TeachingDemos)
library(gridExtra)

1. The SIR model with vaccination in open population

Incorporating a vaccination program in the SIR model can be done by adding a direct path from the susceptible compartment to the immune compartment. We illustrate this using a SIR model that is age-independent and we look at its dynamic behavior over time. Let us assume that a proportion \(p\) of each cohort of newborns is vaccinated at birth, the SIR model with demography can then be written as

\[ \left\{ \begin{array}{l} \frac{dS}{dt} = N\mu(1-p) - (\lambda(t)+\mu)S,\\ \frac{dI}{dt} =(\lambda(t)+\mu)S - (\nu+ \alpha+\mu)I,\\ \frac{dR}{dt} =N \mu p + \nu I - \mu R. \end{array} \right. \]

The SIR model formulated above implies that, for each cohort, a proportion p enter the population directly in the immune class and therefore do not take a part in the transmission process.

2. The SIR model with vaccination in R: post vaccination equilibrium

Recall that the basic reproductive number, \(R_{0}\), is defined as

\[ R_{0}=\frac{N \times \beta}{\nu+\mu}. \]

If \(R_{0}\) is above the threshold of one , the disease will maintain itself while for \(R_{0} < 1\) the disease will die out.

The Critical Proportion of Vaccination

The direct effect of a vaccination program implies that at birth a proportion \(p\) of each cohort of new born is transferred directly to the immune class. This means that in the new equilibrium after vaccination the proportion of the population in the susceptible class is at most \(1-p\). The basic reproductive rate at equilibrium is given by \[ R_{0}=\frac{1}{\tilde{S}(\infty)}, \] or \[ R_{0}\tilde{S}(\infty)=1. \]

Hence, at equilibrium in a vaccinated population \[ R_{0}\tilde{S}(\infty)=R_{0}(1-p). \] For \(R_{0}(1-p) < 1\) the infection will not be able to maintain itself. Therefore, the critical proportion of individuals that has to be vaccinated in order to eliminate the disease, \(p_{c}\), is the one satisfied \(R_{0}(1-p_{c}) < 1\) or \[ p_{c}=1-\frac{1}{R_{0}}. \]

Post vaccination parameters

The force of infection and the basic reproductive number

For a vaccination program (\(p\) vaccinated at birth) the post vaccination force of infection at the new equilibrium decreases linearly with \(p\) and given by

\[ \lambda'=\mu R_{0}(p_{c}-p), \]

and the post vaccination basic reproductive rate is equal to \[ R_{0}=\frac{\lambda'+ \mu}{(1-p) \mu}. \]

As \(p \rightarrow p_{c}\) then \(\lambda' \rightarrow 0\) and the disease will be eliminated. Note that, \(1-p\) is the proportion of individuals who are not vaccinated but will not transfer to the infected compartment since the force of infection is converging to 0.

The Average Age at Infection

The new average age at infection, i.e., the average time spent in the susceptible compartment is given by

\[ A`=\frac{A}{1-p}. \]

3. Example 1: a model without vaccination (\(p=0\))

library(deSolve)

Our starting point is a scenario of a population without vaccination. In this case, the SIR model is given by

\[ \left\{ \begin{array}{l} \frac{dS}{dt} = N\mu - (\lambda(t)+\mu)S,\\ \frac{dI}{dt} =(\lambda(t)+\mu)S - (\nu+ \alpha+\mu)I,\\ \frac{dR}{dt} = \nu I - \mu R. \end{array} \right. \]

This implies that individuals are entering to the immune compartment from the infected compartment.

Defintion of the model’s parameters

We assume that life expectancy is 75 years, \(\mu=\frac{1}{75}, \beta=0.0005\) and \(\nu=1\). For a population without vaccination we assume \(p=0\). Note that in the object parameters, the vaccination parameter is p=0, i.e., vaccination coverage of \(0\%\).

parameters <- c(mu=1/75,beta=0.001/2, v=1, P=0.0)
parameters
##         mu       beta          v          P 
## 0.01333333 0.00050000 1.00000000 0.00000000

We assume population size is 5000 and that at \(t=0\), \((S(0)=4999,I(0)=1, R(0)=0)\).

state <- c(X=4999,Y=1,Z=0)
state
##    X    Y    Z 
## 4999    1    0
times<-seq(0,1600,by=0.01)
lastt<-length(times)

For this parameter setting, the basic reproductive number \(R_{0}\) is equal to \[ R_{0}=\frac{N \times \beta}{\nu+\mu}=\frac{5000 \times 0.0005}{1+\frac{1}{75}} \]

(5000*0.0005)/(1+1/75)
## [1] 2.467105

Definition of the model in R

We define a function, SIR which is used to implement the model formulated above in R.

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

Running the model

We use the ode to solve the equation system of the SIR model.

require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
head(out)
##   time        X        Y          Z
## 1 0.00 4999.000 1.000000 0.00000000
## 2 0.01 4998.975 1.014973 0.01007419
## 3 0.02 4998.950 1.030170 0.02029788
## 4 0.03 4998.924 1.045594 0.03067330
## 5 0.04 4998.898 1.061249 0.04120244
## 6 0.05 4998.871 1.077138 0.05188798

Model’s output

Equilibrium values

Without vaccination, \(S_{\infty}=2026.667\), \(I_{\infty}=39.12281\) and \(R_{\infty}=2934.211\) (for more details see Section 5 below). Note that the R object lastt is the last time point for the numerical integration of the equation system of the model. These values can be clearly seen in Figure 1 and 2.

out[lastt,]
##        time        X        Y        Z
## 160001 1600 2026.667 39.12281 2934.211

Graphical output

Susceptible and infected individuals

par(mfrow=c(2,1))
plot (times,out$X ,type="l",main="S", xlab="time", ylab="-")
plot (times,out$Y ,type="l",main="I", xlab="time", ylab="-")
SIR model without vaccination (p=0).

Figure 1: SIR model without vaccination (p=0).

The equilibrium plot

par(mfrow=c(1,1))
plot(out$X,out$Y,type="l",xlab="S",ylab="I")
Equilibrium plot for a SIR model without vaccination (p=0).

Figure 2: Equilibrium plot for a SIR model without vaccination (p=0).

4. Example 2: model with \(40\%\) vaccination at birth

Definition of the model parameters

In this example we assume that \(40\%\) of the population are vaccinated at birth. The state variables remain the same as before

state <- c(X=4999,Y=1,Z=0)
state
##    X    Y    Z 
## 4999    1    0

but we need to specify the proportion of vaccination. For a population size of 5000, \(\beta = 0.0005\), \(\nu = 1\), life expectancy of 75 years and vaccination vaccination rate of \(40\%\) at birth, \(p=0.4\) we have

parameters <- c(mu=1/75,beta=0.001/2, v=1, P=0.40)
parameters
##         mu       beta          v          P 
## 0.01333333 0.00050000 1.00000000 0.40000000

Model solution

Figure 3 shows the number of susceptible over time for two scenarios with \(p = 0\) and \(p=0.4\) and reveals the post vaccination equilibrium patterns. The inter epidemic period is longer when the proportion of vaccinated individuals increases. This implies that individuals stay, on average, longer time in the susceptible class and as a result the average age in which individuals are infected (the average age at infection) increases. In the long run, after vaccination, the population will research a new endemic equilibrium. Notice that although the pre-vaccination force of infection \(\lambda(0)\) is equal under the two scenarios, in the new equilibrium the force of infection decreases as the proportion of vaccinated individuals increases. Since $(t) decreases as \(p\) increases the rate in which individuals leave the susceptible class decreases as \(p\) increases and as a result, as shown in Figure 4 the inter epidemic period is longer as \(p\) increases. This implies that a vaccination program has two effects on the population:

  • The direct effect is the transfer of a proportion p from the susceptible directly to the immune class.

  • The indirect effect is related to the decline of the magnitude of the force of infection and its affect both vaccinated and un vaccinated individuals.

require(deSolve)
outp40 <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
head(outp40)
##   time        X        Y         Z
## 1 0.00 4999.000 1.000000 0.0000000
## 2 0.01 4998.708 1.014972 0.2767231
## 3 0.02 4998.416 1.030167 0.5535601
## 4 0.03 4998.124 1.045588 0.8305133
## 5 0.04 4997.831 1.061237 1.1075846
## 6 0.05 4997.538 1.077120 1.3847767
par(mfrow=c(1,1))
plot (times,outp40$X ,type="l",main="S", xlab="time", ylab="-")
lines(times,out$X,col=2)
legend(250,4500,c("P=0","P=0.4"),lty=c(1,1),col=c(2,1))
SIR model with vaccination (p=0.4, black line): Susceptibale.

Figure 3: SIR model with vaccination (p=0.4, black line): Susceptibale. 

par(mfrow=c(1,1))
plot (times,outp40$Y ,type="l",main="I", xlab="time", ylab="-")
lines(times,out$Y,col=2)
legend(250,900,c("P=0","P=0.4"),lty=c(1,1),col=c(2,1))
SIR model with vaccination (p=0.4, black line): Infected.

Figure 4: SIR model with vaccination (p=0.4, black line): Infected.

Post vaccination equilibrium values (model’s output)

As can be seen in Figure 5, the equilibrium the number of susceptible for \(p = 0\) and \(p = 0.4\) will be the same (2026.667). However, the number of infectious individuals (lower panel) for the case with \(p = 0.4\;(12.807)\) will be smaller than the number of infectious individuals for \(p = 0\;(39.122)\).

out[lastt,]
##        time        X        Y        Z
## 160001 1600 2026.667 39.12281 2934.211
outp40[lastt,]
##        time        X        Y        Z
## 160001 1600 2026.667 12.80705 2960.526
plot(outp40$X,outp40$Y,type="l",ylim=c(0,1250))
lines(out$X,out$Y,col=2)
legend(4000,900,c("P=0","P=0.4"),lty=c(1,1),col=c(2,1))
Post vaccination equilibrium with p=0,0.4.

Figure 5: Post vaccination equilibrium with p=0,0.4.

5. Example 3: The critical proportion of vaccination

Vaccination rate of \(50\%\)

For a population with \(\mu=\frac{1}{75}\), \(\beta=0.0005\), and \(\nu=1\), the basic reproductive number is equal to 2.47. In this case, the critical proportion of vaccination is equal to

\[ P_{c}=1-\frac{1}{2.46}=0.593 \]

First, we implement a vaccination program in which \(50\%\) of the population are vaccinated at birth. Figure 6 shows the susceptible at each time under the scenarios with \(40\%\) and \(50\%\) vaccination coverage.

require(deSolve)
parameters <- c(mu=1/75,beta=0.001/2, v=1, P=0.5)
parameters
##         mu       beta          v          P 
## 0.01333333 0.00050000 1.00000000 0.50000000
outp50 <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
head(outp50)
##   time        X        Y         Z
## 1 0.00 4999.000 1.000000 0.0000000
## 2 0.01 4998.642 1.014972 0.3433853
## 3 0.02 4998.283 1.030166 0.6868756
## 4 0.03 4997.924 1.045586 1.0304732
## 5 0.04 4997.565 1.061235 1.3741801
## 6 0.05 4997.205 1.077115 1.7179988
par(mfrow=c(1,1))
plot (times,outp50$X ,type="l",main="S", xlab="time", ylab="-")
lines(times,outp40$X,col=2)
legend(250,4500,c("P=0.4","P=0.5"),lty=c(1,1),col=c(2,1))
SIR model with vaccination (p=0.5, black line): Susceptibale.

Figure 6: SIR model with vaccination (p=0.5, black line): Susceptibale.

Figure 7 shows the number of infected individuals under the three scenarios with \(p=0,0.4,0.5\). Notice how the inter-epidemic period increase as \(p\) increases.

require(deSolve)
par(mfrow=c(1,1))
plot (times,outp50$Y ,type="l",main="I", xlab="time", ylab="-")
lines(times,outp40$Y,col=2)
lines(times,out$Y,col=3)
legend(250,1000,c("P=0","P=0.4","P=0.5"),lty=c(1,1),col=c(3,2,1))
SIR model with vaccination (p=0.5, black line): Susceptibale.

Figure 7: SIR model with vaccination (p=0.5, black line): Susceptibale.

Equilibrium values

As shown in Figure 8, for all vaccination coverage, \(S_{\infty}=2026.67\) but \(I_{\infty}\) decrease as \(p\) increases. As a results, the inter-epidemic period that was observed in Figure 7 increases with \(p\).

out[lastt,]
##        time        X        Y        Z
## 160001 1600 2026.667 39.12281 2934.211
outp40[lastt,]
##        time        X        Y        Z
## 160001 1600 2026.667 12.80705 2960.526
outp50[lastt,]
##        time       X        Y        Z
## 160001 1600 2026.66 6.227879 2967.112
plot(outp50$X,outp50$Y,type="l",ylim=c(0,1250))
lines(outp40$X,outp40$Y,col=2)
lines(out$X,out$Y,col=3)
legend(4000,1200,c("P=0","P=0.4","P=0.5"),lty=c(1,1),col=c(3,2,1))
Post vaccination equilibrium with p=0.4 and p=0.5.

Figure 8: Post vaccination equilibrium with p=0.4 and p=0.5.

Vaccination rate of \(60\%\)

In this example, \(60\%\) of the population are vaccinated at birth, i.e. \(p > p_{c}\). Under this scenario we expect that the disease will die out. Figure 9 presents the number of susceptible over time under the three scenarios and shows that for \(p=0.6\), after an initial outbreak, the susceptible class cannot build up to exceed the threshold for a second outbreak and the number of susceptible will be equal to \(p \times N\) (2000) in our example a can be seen in Figure 10 that presents a zoom-in results between 1000 to 1400 times units.

require(deSolve)
parameters <- c(mu=1/75,beta=0.001/2, v=1, P=0.6)
parameters
##         mu       beta          v          P 
## 0.01333333 0.00050000 1.00000000 0.60000000
outp60 <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
head(outp60)
##   time        X        Y         Z
## 1 0.00 4999.000 1.000000 0.0000000
## 2 0.01 4998.575 1.014972 0.4100475
## 3 0.02 4998.150 1.030166 0.8201912
## 4 0.03 4997.724 1.045585 1.2304332
## 5 0.04 4997.298 1.061232 1.6407756
## 6 0.05 4996.872 1.077111 2.0512210
par(mfrow=c(1,1))
plot (times,outp60$X ,type="l",main="S", xlab="time", ylab="-")
lines(times,outp50$X,col=2)
lines(times,outp40$X,col=3)
legend(250,4500,c("P=0.4","P=0.5","P=0.6"),lty=c(1,1),col=c(3,2,1))
SIR model with vaccination (p=0.6, black line): Susceptibale.

Figure 9: SIR model with vaccination (p=0.6, black line): Susceptibale. 

tlo<-100000
tup<-140000
par(mfrow=c(1,1))
plot (times[tlo:tup],outp60$X[tlo:tup],ylim=c(1900,2100) ,type="l",main="S", xlab="time", ylab="-")
lines(times[tlo:tup],outp50$X[tlo:tup],col=2)
lines(times[tlo:tup],outp40$X[tlo:tup],col=3)
legend(1100,1975,c("P=0.4","P=0.5","P=0.6"),lty=c(1,1),col=c(3,2,1))
SIR model with vaccination (p=0.6, black line): Susceptibale, 1000 <t < 1400.

Figure 10: SIR model with vaccination (p=0.6, black line): Susceptibale, 1000 <t < 1400.

Figure 11 shows the number of immune. The decrease in the number of immune individuals is due to the mortality rate in the population. Note that, as shown in Figure 12, the number of immune at equilibrium for \(p=0.6\) is equal to 3000.

tlo<-1
tup<-40000
par(mfrow=c(1,1))
plot (times[tlo:tup],outp60$Z[tlo:tup] ,type="l",main="R", xlab="time", ylab="-")
lines(times[tlo:tup],outp50$Z[tlo:tup],col=2)
lines(times[tlo:tup],outp40$Z[tlo:tup],col=3)
legend(50,2000,c("P=0.4","P=0.5","P=0.6"),lty=c(1,1),col=c(3,2,1))
SIR model with vaccination (p=0.6, black line): Immune, 0 <t < 400.

Figure 11: SIR model with vaccination (p=0.6, black line): Immune, 0 <t < 400. 

tlo<-40000
tup<-60000
par(mfrow=c(1,1))
plot (times[tlo:tup],outp60$Z[tlo:tup],ylim=c(2500,3250) ,type="l",main="R", xlab="time", ylab="-")
lines(times[tlo:tup],outp50$Z[tlo:tup],col=2)
lines(times[tlo:tup],outp40$Z[tlo:tup],col=3)
legend(450,2800,c("P=0","P=0.4","P=0.5"),lty=c(1,1),col=c(3,2,1))
SIR model with vaccination (p=0.6, black line): Immune, 400 < t < 600.

Figure 12: SIR model with vaccination (p=0.6, black line): Immune, 400 < t < 600.

Figure 13 shows that the for \(p=0.6\) there is only one initial outbreak.

tup<-40000
par(mfrow=c(1,1))
plot (times[1:tup],outp60$Y[1:tup] ,type="l",main="I", xlab="time", ylab="-")
lines(times[1:tup],outp50$Y[1:tup],col=2)
lines(times[1:tup],outp40$Y[1:tup],col=3)
legend(50,800,c("P=0","P=0.4","P=0.5"),lty=c(1,1),col=c(3,2,1))
SIR model with vaccination (p=0.6, black line): Infected.

Figure 13: SIR model with vaccination (p=0.6, black line): Infected.

Equilibrium values

For a population without vaccination

\[ S_{\infty}=\frac{(\nu+\mu)}{\beta} \]

((1+1/75)/0.0005)
## [1] 2026.667

and

\[ I_{\infty}=\frac{\mu}{\beta}(R_{0}-1) \]

((1/75)/0.0005)*(2.467105-1)
## [1] 39.1228

This can be seen in the models’ output below as well. Note that the object lastt is the last time point that was used for integration.

out[lastt,]
##        time        X        Y        Z
## 160001 1600 2026.667 39.12281 2934.211

For a vaccination program with \(p \ge p_{c}\) the equilibrium values are equal to \(S_{\infty}=(1-P) \times N\), \(R_{\infty}=p \times N\) and \(I_{\infty}=0\). In our example \(P=0.6\) and \(N=5000\). As shown in the model output below, \(S_{\infty}=0.4 \times 5000=2000\), \(R_{\infty}=0.6 \times 5000=3000\) (see also Figure 10, 12 and 13).

outp60[lastt,]
##        time    X            Y    Z
## 160001 1600 2000 3.772673e-30 3000