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

Modeling infectious diseases using R - age structured population

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 age structure population model

Contact pattern in age structure population

A special case of a SIR model with interacting subpopulations is an age-time dependent SIR model or age-structured SIR model. In such a model the population is divided into a finite number of age groups interacting with each other. Note that compared with the model discussed in the previous section, the flow of individuals into the susceptible class is possible only in the first age class. A second difference is that in the current model, individuals can pass to the same compartment in successive age classes, this means that a susceptible individual in age group \(i\) can remain susceptible and transfer (age passes) to the susceptible compartment in the successive age group, i.e., this individual passes from \(S_{i}\) to \(S_{i+1}\).

The contact matrix

The transmission model

For a population with \(K\) age groups, the system of ordinary differential equations for the first age group (\(i=1\)) in an age-structured SIR model is given by

\[ \left \{ \begin{array}{l} \frac{dS_{i}(t)}{dt} = - \left( \Sigma_{j=1}^{K}\beta_{ij}I_{j} \right )S_{i} +N\mu_{i}-\mu_{i}S_{i} -\eta_{i}S_{i},\\ \frac{dI_{i}(t)}{dt} = \left( \Sigma_{j=1}^{K}\beta_{ij}I_{j} \right )S_{i} - \nu_{i} I_{i} -\eta_{i}I_{i},\\ \frac{dR_{i}(t)}{dt} = \nu_{i}I_{i}-\mu_{i}R_{i}-\eta_{i}R_{i}. \end{array} \right. \]

Here, \(\eta_{i}\) is the rate at which individuals pass from \(S_{i},I_{i},R_{i}\) to \(S_{i+1},I_{i+1},R_{i+1}\). From the second age group onwards, the system of differential equations is similar to the system in one population but without births into the susceptible class (i.e. without \(N \mu_{i}\)) and with the additional flows \(\eta_{i-1}S_{i-1}\), \(\eta_{i-1}I_{i-1}\) and \(\eta_{i-1}R_{i-1}\) corresponding to the susceptible, infected and recovered compartment from the previous age group, respectively.

\[ \left \{ \begin{array}{l} \frac{dS_{i}(t)}{dt} = - \left( \Sigma_{j=1}^{K}\beta_{ij}I_{j} \right )S_{i} +\eta_{i-1}S_{i-1}-\mu_{i}S_{i} -\eta_{i}S_{i},\\ \frac{dI_{i}(t)}{dt} = \left( \Sigma_{j=1}^{K}\beta_{ij}I_{j} \right )S_{i} +\eta_{i-1}I_{i-1}- (\nu_{i}+\mu) I_{i}-\eta_{i}I_{i},\\ \frac{dR_{i}(t)}{dt} = \nu_{i}I_{i}-\mu_{i}R_{i}+\eta_{i-1}R_{i-1}. \end{array} \right. \]

Contact matrix age-structured SIR model with two age groups

The age-structured SIR model introduces the challenge of estimating the mixing matrix which, up to this point, was assumed to be known. In practice the mixing matrix is unknown and should be estimated. In order to estimate the WAIFW matrix, Anderson and May (1985) introduced a framework in which the mixing matrix itself was assumed unknown but its structure, the mixing pattern, was assumed to be known. Let us consider an age-structured SIR model with two age groups, \([0,a_{1})\) and \([a_{1},L)\), as described above. For each age-group there is an age-specific constant force of infection, \(\lambda_{i}, \; i=1,2\). Using the mixing matrix above it follows that \[ \left ( \begin{array}{c} \lambda_{1} \\ \lambda_{2} \end{array} \right ) = \left( \begin{array}{cc} \beta_{11} & \beta_{12}\\ \beta_{21} & \beta_{22} \end{array} \right ) \left ( \begin{array}{c} I_{1} \\ I_{2} \end{array} \right ). \]

The force of infection

Let us assume that both \(\lambda_{i}\) and \(I_{i}\) were estimated from pre vaccination cross-sectional serological data. In that case we can plug in the estimates \(\hat{I}_{1}\) and \(\hat{I}_{2}\) and it follows that

\[ \begin{array}{c} \hat{\lambda}_{1} = \beta_{11} \hat{I}_{1} + \beta_{12} \hat{I}_{2},\\ \hat{\lambda}_{2} = \beta_{21}\hat{I}_{1} + \beta_{22} \hat{I}_{2}. \end{array} \]

2. Example I: a model with two age groups

An age structure model in R

Definition of the model’s papameters

We consider the following transmission rates for a population with two age group \(\beta_11=\beta_22=0.0001\) and \(\beta_12=\beta_{21}=0.0000075\). The contact matrix is given by

\[ C= \left( \begin{array}{cc} 0.0001 & 0.0000075\\ 0.0000075 & 0.0001 \end{array} \right). \]

We assume two age groups, 0-20 and 20-75, life expectancy of 75 years (birth rate=death rate (mu =\(\frac{1}{75}\)) and a population size of 1000000 individuals. In R, the parameter vector is given by

parameters <- c(beta11=0.0001,beta12=0.0000075,beta21=0.0000075,beta22=0.0001,
                v1=4,v2=4,mu=1/75,mu2=1/20,N=1000000)
parameters
##       beta11       beta12       beta21       beta22           v1           v2 
## 1.000000e-04 7.500000e-06 7.500000e-06 1.000000e-04 4.000000e+00 4.000000e+00 
##           mu          mu2            N 
## 1.333333e-02 5.000000e-02 1.000000e+06

Note that the object mu2 is the rate in which individuals move from \(S_{i-1}\) to \(S_{i}\), \(I_{i-1}\) to \(I_{i}\) and \(R_{i-1}\) to \(R_{i}\) (\(\eta\)) and in our example it is equal to \(\frac{1}{20}\) since the first age group is 0-20.

The state vector

In the first age group, \(S(0)=266665, I(0)=1, R(0)=0\) and in the second age group \(S(0)=733334, I(0)=0, R(0)=0\). The state vector in R is

state <- c(Y1=266665,Y2=1,Y3=0,Y4=733334,Y5=0.0,Y6=0)
state
##     Y1     Y2     Y3     Y4     Y5     Y6 
## 266665      1      0 733334      0      0
times<-seq(0,60,by=0.01)

The transmission model

The SIR model with two age groups is implemented using the R function SIRtwo.

SIRtwo<-function(t,state,parameters)
{
with(as.list(c(state, parameters)),
{
dY1 <- -(beta11*Y2+beta12*Y5)*Y1+N*mu-mu*Y1-mu2*Y1
dY2 <- (beta11*Y2+beta12*Y5)*Y1-v1*Y2-mu*Y2-mu2*Y2
dY3 <- v1*Y2 - mu*Y3-mu2*Y3
dY4 <- -(beta21*Y2+beta22*Y5)*Y4-mu*Y4+mu2*Y1
dY5 <-  (beta21*Y2+beta22*Y5)*Y4-v2*Y5-mu*Y5+mu2*Y2
dY6 <- v2*Y5-mu*Y6+mu2*Y3
list(c(dY1,dY2,dY3,dY4,dY5,dY6))
}) 
}

Running the model

require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIRtwo,parms=parameters))
out1<-out #Example 1
head(out1)
##   time       Y1       Y2         Y3       Y4         Y5          Y6
## 1 0.00 266665.0 1.000000 0.00000000 733334.0 0.00000000 0.000000000
## 2 0.01 266629.2 1.254411 0.04487737 733369.5 0.08870454 0.001535006
## 3 0.02 266593.3 1.576325 0.10119126 733404.8 0.28878949 0.008602739
## 4 0.03 266557.3 1.986445 0.17202439 733439.8 0.71777656 0.027740930
## 5 0.04 266521.2 2.514482 0.26145085 733474.4 1.61272323 0.072303628
## 6 0.05 266484.9 3.205256 0.37500523 733507.9 3.45075564 0.169375119

Graphical output

Figure 1 and 2 show the number of susceptible and infected over time, respectively.

par(mfrow=c(1,1))
plot(times,out1$Y1+out1$Y4,ylim=c(0,500000),type="l",main=" ",xlab="time",ylab="S")
lines(times,out1$Y1,lty=2)
lines(times,out1$Y4,lty=3)
legend(10,500000,c("total","age group I","age group II"),lty=c(1,2,3))
SIR model with two age groups: susceptible.

Figure 1: SIR model with two age groups: susceptible. 

par(mfrow=c(1,1))

plot(times,out1$Y2+out1$Y5,ylim=c(0,50000),type="l",main=" ",xlab="time",ylab="I")
lines(times,out1$Y2,lty=2)
lines(times,out1$Y5,lty=3)
legend(10,50000,c("total","age group I","age group II"),lty=c(1,2,3))
SIR model with two age groups: infected.

Figure 2: SIR model with two age groups: infected.

Example 2: a change the parameter setting

A new contact matrix

Let us assume that the transmission rate within the first age group reduces from \(\beta{11}=0.0001\) to \(\beta_{11}=0.000001\) while the other parameters remain constant. The new contact matrix is given by

\[ C= \left( \begin{array}{cc} 0.000001 & 0.0000075\\ 0.0000075 & 0.0001 \end{array} \right). \]

in R, we updated the parameters vector

parameters <- c(beta11=0.000001,beta12=0.0000075,beta21=0.0000075,beta22=0.0001,
                v1=4,v2=4,mu=1/75,mu2=1/20,N=1000000)
parameters
##       beta11       beta12       beta21       beta22           v1           v2 
## 1.000000e-06 7.500000e-06 7.500000e-06 1.000000e-04 4.000000e+00 4.000000e+00 
##           mu          mu2            N 
## 1.333333e-02 5.000000e-02 1.000000e+06

Running the model

require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIRtwo,parms=parameters))
out05<-out #Example 2
head(out05)
##   time       Y1        Y2         Y3       Y4         Y5          Y6
## 1 0.00 266665.0 1.0000000 0.00000000 733334.0 0.00000000 0.000000000
## 2 0.01 266629.5 0.9634372 0.03924640 733369.5 0.07876116 0.001411269
## 3 0.02 266593.9 0.9304458 0.07707275 733404.8 0.23349479 0.007336380
## 4 0.03 266558.4 0.9030501 0.11365843 733440.0 0.54055317 0.022164270
## 5 0.04 266522.9 0.8853659 0.14930072 733474.8 1.15291406 0.054693785
## 6 0.05 266487.4 0.8856963 0.18453536 733509.0 2.37708411 0.122558953

Graphical output

Figure 3 and 4 show the number of susceptible and infected individuals in the population under the new scenario.

par(mfrow=c(1,1))
plot(times,out05$Y1+out05$Y4,ylim=c(0,500000),type="l",main=" ",xlab="time",ylab="S")
lines(times,out05$Y1,lty=2)
lines(times,out05$Y4,lty=3)
legend(10,500000,c("total","age group I","age group II"),lty=c(1,2,3))
SIR model with two age groups: susceptible, beta_11=0.000001.

Figure 3: SIR model with two age groups: susceptible, beta_11=0.000001. 

plot(times,out05$Y2+out05$Y5,ylim=c(0,50000),type="l",main=" ",xlab="time",ylab="I")
lines(times,out05$Y2,lty=2)
lines(times,out05$Y5,lty=3)
legend(10,50000,c("total","age group I","age group II"),lty=c(1,2,3))
SIR model with two age groups: infected (beta_11=0.000001.

Figure 4: SIR model with two age groups: infected (beta_11=0.000001.

Comparison between the two examples

Figure 5 shows the total number of susceptible under the new scenarios. note that due to the reduction in the transmission rate in the first age group, the susceptible class is buldit slowly and it takes more time to observed an outbreak. In addition, the new equilibrium is reached with more individual in the susceptible class as can be seen in Figure 7 as well.

par(mfrow=c(1,1))
plot(times,out1$Y1+out1$Y4,ylim=c(0,500000),type="l",main=" ",xlab="time",ylab="S")
lines(times,out05$Y1+out05$Y4,lty=2)
legend(10,500000,c("Example 1 (beta_11=0.0001)","Example 2 (beta_11=0.000001)"),lty=c(1,2))
SIR model with two age groups: susceptible.

Figure 5: SIR model with two age groups: susceptible.

The longer inter-epidemic period under scenario 2 can be clearly seen in Figure 6.

par(mfrow=c(1,1))
plot(times,out1$Y2+out1$Y5,ylim=c(0,50000),type="l",main=" ",xlab="time",ylab="S")
lines(times,out05$Y2+out05$Y5,lty=2)
legend(10,50000,c("Example 1 (beta_11=0.0001)","Example 2 (beta_11=0.000001)"),lty=c(1,2))
SIR model with two age groups: infected.

Figure 6: SIR model with two age groups: infected. 

par(mfrow=c(1,2))
plot(out1$Y1,out1$Y2,type="l")
lines(out05$Y1,out05$Y2,lty=2)
title("age group 1")
plot(out1$Y4,out1$Y5,type="l")
lines(out05$Y4,out05$Y5,lty=2)
title("age group 2")
Equilibrium plot.

Figure 7: Equilibrium plot.