Handling increasing complexities for Infectious Disease Modeling using R
Dr Gurpreet Singh, Dr Biju Soman
2020-11-20
- 1 Introduction.
- 2 Revision: Framework for mathematical modeling
- 3 Increasing complexity
- 4 Increasing complexities: Addition of interventions.
- 5 Increasing complexities: Modeling age groups in SIR models.
- 6 Conclusion.
1 Introduction.
This document is an attempt to demonstrate the handling of increasing complexities in IDM using R. It is further to the last session wherein we learned about using deSolve package for development and analysis of Infected Cohort models.
2 Revision: Framework for mathematical modeling
It is considered as a good practice to draw the model by hand/ by using any computer software for conceptualising the model. Following steps can be considered as minimum framework/ steps required for modelling:-
- Evidence synthesis
- Create model input using parameters
- Create a model using input values
- Return a dataframe as output using model equations
- Exploratory Data Analysis of the output
3 Increasing complexity
To demonstrate increasing complexity, let us now develop a model for a disease wherein three compartments viz susceptible, infected, and recovered are present. There are no births (inflows) or deaths (outflows) in the model. Further, all those who are infected are also infectious, and once recovered, immunity is long lasting. The basic assumptions of compartment models such as homogeneity and random mixing are also considered.
3.1 Visualisation of the model
Figure 3.1: SIR model
3.2 Adding changes in the IR model to develop SIR model
We shall develop an SIR model for a susceptible population of 10000000 wherein one infectious case with infectious period of seven days and R0 of 2.28 has led to an outbreak of the disease. The individuals in the population are susceptible at outset, gets infected followed by recovery. We shall follow the population dynamics at daily interval for 300 days.
## 1. Evidence synthesis
initial_number_susceptible <- 10000000 # Adding susceptibles to IR model
initial_number_infected <- 1 # Considering one infected at start time
initial_number_recovered <- 0
Duration_of_infectiousness <- 7
recovery_rate <- 1/7
Reproduction_number <- 2.28 #R0== beta/gamma,
infection_rate <- Reproduction_number*(1/Duration_of_infectiousness) #beta
follow_up_duration <- 300
## 2. Create model input using parameters
initial_state_values <- c(S = initial_number_susceptible, # Adding susceptibles
I = initial_number_infected,
R = initial_number_recovered)
parameters <- c(beta = infection_rate, # Adding beta parameter
gamma = recovery_rate)
times <- seq(from = 0, to = follow_up_duration, by = 1)
## 3. Create a model using input values
cohort_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S+I+R # Total population
lambda <- beta * I/N #Calculating lambda as function of beta
dS <- -lambda * S
dI <- lambda * S - gamma * I
dR <- gamma * I
return(list(c(dS, dI, dR)))
})
}
## 4. Return a dataframe as output using model equations
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = cohort_model,
parms = parameters))
## 5. Exploratory Data Analysis (Visualisation) of the output
dat <- output %>% pivot_longer(cols = c(S,I,R),
names_to = "state",
values_to = "value")
dat %>% ggplot(aes(x = time,
y = value,
color = state,
group = state))+
geom_line()+
xlab("Time(Days)")+
ylab("Number of persons")+
labs(title = "Changes in Suceptible, Infected and Recovered")
4 Increasing complexities: Addition of interventions.
Let us now assume that we have a vaccine with 100% effectiveness against acquiring infection.
4.1 Will the epidemic be averted by vaccinating 10% population before the infectious individual enters the population?
## 1. Evidence synthesis
prop_vaccinated <- 0.1
initial_number_susceptible <- (1-prop_vaccinated)*10000000
initial_number_infected <- 1
initial_number_recovered <- prop_vaccinated*10000000
Duration_of_infectiousness <- 7
recovery_rate <- 1/7
Reproduction_number <- 2.28 #R0== beta/gamma,
infection_rate <- Reproduction_number*(1/Duration_of_infectiousness) #beta
follow_up_duration <- 300
## 2. Create model input using parameters
initial_state_values <- c(S = initial_number_susceptible,
I = initial_number_infected,
R = initial_number_recovered)
parameters <- c(beta = infection_rate,
gamma = recovery_rate)
times <- seq(from = 0, to = follow_up_duration, by = 1)
## 3. Create a model using input values
cohort_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S+I+R # Total population
lambda <- beta * I/N #Calculating lambda as function of beta
dS <- -lambda * S
dI <- lambda * S - gamma * I
dR <- gamma * I
return(list(c(dS, dI, dR)))
})
}
## 4. Return a dataframe as output using model equations
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = cohort_model,
parms = parameters))
## 5. Exploratory Data Analysis (Visualisation) of the output
dat <- output %>% pivot_longer(cols = c(S,I,R),
names_to = "state",
values_to = "value")
dat %>% ggplot(aes(x = time,
y = value,
color = state,
group = state))+
geom_line()+
xlab("Time(Days)")+
ylab("Number of persons")+
labs(title = "Changes in Suceptible, Infected and Recovered at 10% vaccine coverage")
4.2 Will the epidemic be averted by increasing vaccination coverage? What is the critical vaccination threshold?
Critical Vaccination Threshold is defined as minimum percentage of population which needs to be vaccinated to prevent occurence of epidemic. Going manually by changing the vaccination coverage is a start point in understanding the concept.
4.2.1 Modelling for 20% coverage:
## 1. Evidence synthesis
prop_vaccinated <- 0.2
initial_number_susceptible <- (1-prop_vaccinated)*10000000
initial_number_infected <- 1
initial_number_recovered <- prop_vaccinated*10000000
Duration_of_infectiousness <- 7
recovery_rate <- 1/7
Reproduction_number <- 2.28 #R0== beta/gamma,
infection_rate <- Reproduction_number*(1/Duration_of_infectiousness) #beta
follow_up_duration <- 300
## 2. Create model input using parameters
initial_state_values <- c(S = initial_number_susceptible, # Adding susceptibles
I = initial_number_infected,
R = initial_number_recovered)
parameters <- c(beta = infection_rate, # Adding beta parameter
gamma = recovery_rate)
times <- seq(from = 0, to = follow_up_duration, by = 1)
## 3. Create a model using input values
cohort_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S+I+R # Total population
lambda <- beta * I/N #Calculating lambda as function of beta
dS <- -lambda * S
dI <- lambda * S - gamma * I
dR <- gamma * I
return(list(c(dS, dI, dR)))
})
}
## 4. Return a dataframe as output using model equations
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = cohort_model,
parms = parameters))
## 5. Exploratory Data Analysis (Visualisation) of the output
dat <- output %>% pivot_longer(cols = c(S,I,R),
names_to = "state",
values_to = "value")
dat %>% ggplot(aes(x = time,
y = value,
color = state,
group = state))+
geom_line()+
xlab("Time(Days)")+
ylab("Number of persons")+
labs(title = "Changes in Suceptible, Infected and Recovered at 20% vaccine coverage")
4.2.2 Modelling for 30% coverage:
4.2.3 Modelling for 40% coverage:
4.2.4 Modelling for 50% coverage:
5 Increasing complexities: Modeling age groups in SIR models.
Now, we shall learn on mechanisms to structure SIR model into two age groups and how to define the force of infection in a population with age-specific mixing.
Rationale: From understanding of disease epidemiology, we came to know that age-specific mixing patterns are different among age groups and thus different ages experience different forces of infection. In the model, susceptible children shall experience a force of infection \(\lambda 1\) and susceptible adults become infected at a rate \(\lambda 2\) .
Assumptions: There is no ageing in the population. Hence, there are no arrows going from compartments 1 (children) to compartments 2 (adults). Additionally, we assume that children and adults all recover from the infection at the same rate \(\gamma\).
Remember, in previous models we had summarized the infection rate with the \(\beta\) parameter. However, to model age-specific mixing patterns in a population, it is helpful to break \(\beta\) down into its 2 key components:
\(\beta\) : Probability of infection per contact which is same in all age groups and equals 5%.
\(c_{ij}\) : Number of contacts that a susceptible person in age group i makes with people in age group j per unit time.
Additional Information. We assume that, children make 13 contacts per day on average, of which 7 are with other children, and that adults make 11 contacts per day on average, of which 10 are with other adults. This gives the following contact parameters:
\(c_{11}\)=7
\(c_{12}\)=6
\(c_{21}\)=1
\(c_{22}\)=10
People remain infected on average for 5 days.
The total population size is 1 million, and children consitute 20% of total population.
5.1 Visualisation of the model
Figure 5.1: SIR model
## 1. Evidence synthesis
S1 = 200000 # 20% of the population are children - all susceptible
I1 = 1 # the outbreak starts with 1 infected person (can be either child or adult)
R1 = 0
S2 = 800000 # 100%-20% of the population are adults - all susceptible
I2 = 0
R2 = 0
b = 0.05 # the probability of infection per contact is 5%
c_11 = 7 # daily number of contacts that children make with each other
c_12 = 6 # daily number of contacts that children make with adults
c_21 = 1 # daily number of contacts that adults make with children
c_22 = 10 # daily number of contacts that adults make with each other
gamma = 1/5
## 2. Create model input using parameters
initial_state_values <- c(S1 = 200000, # 20% of the population are children - all susceptible
I1 = 1, # the outbreak starts with 1 infected person (can be either child or adult)
R1 = 0,
S2 = 800000, # 100%-20% of the population are adults - all susceptible
I2 = 0,
R2 = 0)
parameters <- c(b = 0.05, # the probability of infection per contact is 5%
c_11 = 7, # daily number of contacts that children make with each other
c_12 = 6, # daily number of contacts that children make with adults
c_21 = 1, # daily number of contacts that adults make with children
c_22 = 10, # daily number of contacts that adults make with each other
gamma = 1/5) # the rate of recovery is 1/5 per day
times <- seq(from = 0, to = 90, by = 0.1)
## 3. Create a model using input values
sir_age_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N1 <- S1+I1+R1 # the total number of children in the population
N2 <- S2+I2+R2 # the total number of adultsin the population
# Defining the force of infection
# Force of infection acting on susceptible children:
lambda_1 <- b * c_11 * I1/N1 + b * c_12 * I2/N2
# Force of infection acting on susceptible adults:
lambda_2 <- b * c_21 * I1/N1 + b * c_22 * I2/N2
# The differential equations
# Rate of change in children:
dS1 <- -lambda_1 * S1
dI1 <- lambda_1 * S1 - gamma * I1
dR1 <- gamma * I1
# Rate of change in adults:
dS2 <- -lambda_2 * S2
dI2 <- lambda_2 * S2 - gamma * I2
dR2 <- gamma * I2
# Output
return(list(c(dS1, dI1, dR1, dS2, dI2, dR2)))
})
}
## 4. Return a dataframe as output using model equations
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_age_model,
parms = parameters))
## 5. Exploratory Data Analysis (Visualisation) of the output
# Turn output into long format
output_long <- output %>% pivot_longer(cols = c(S1,I1,R1, S2, I2,R2),
names_to = "state",
values_to = "value")
# Plot number of people in all compartments over time
ggplot(data = output_long,
aes(x = time, y = value, colour = state, group = state)) +
geom_line() +
xlab("Time (days)")+
ylab("Number of people") +
labs(colour = "Compartment") 
6 Conclusion.
To conclude, the complexity of the models will increase further as more parameters and compartments are added. Incorporating age structure will introduce compartments and processes for each age group. Contact rates between age groups and gender may vary calling for further compartmentalisation.Vector borne diseases shall include separate but interacting compartments between and within host and vectors. However, the basic priciples/ framework shall remain the same. Also, from the existing data, caliberation of parameters such as reproduction number can be carried out using Maximum Log-Likelihood and other methods. Do we never have outbreaks when reproduction number is less than one (think stochastically)..and the journey continues.. Enjoy modelling! Thank You!!