DRAFT Notes: Building on the SIR Model
Infectious Disease Modelling Specialization, Imperial College London
analysisinsightdata
Friday, November 19, 2021
Reference: Variables, Greek Letters et al.
S = number of susceptible individuals in the population (compartment)
I = number of infected individuals in the population (compartment)
R = number of recovered/removed individuals in the population (compartment)
N = total number of individuals in the population
c = average number of contacts a susceptible makes per time unit
p (or b) = probability of a susceptible becoming infected on contact with infected person
beta β = infection rate; function of c * p (or b)
gamma γ = recovery rate
lambda λ = rate of movement from exposed to infectious (beta * I/N)
R0 R0 R-naught = basic reproduction number of infectious agents (beta / gamma)
Reff Reff Rt = Effective Reproduction Number (R0 * S/N)
pvacc (or p) = vaccination rate or coverage
pvacc_threshold pc = critical vaccination threshold; (1 - (1/R0)); or (1 - (1/R0))/(1 - cs); or (1 - (1/R0))/(1 - (cs * ci))
veff veff = vaccine efficacy; reduces susceptibility by (1 - cs); reduces infection by (1 - ci)
peff peff = effective vaccine coverage (in all-or-nothing = vaccination rate p * veff vaccine efficacy)
cs c_s cS = effectiveness of the vaccine in reducing susceptibility of vaccinated individuals; vaccination reduces susceptibility by (1 - cs)
ci c_i ci = effectiveness of the vaccine in reducing infectivity of vaccinated individuals; vaccination reduces infectivity by (1 - ci)
sigma σ = waning immunity rate (some references uses xi ξ)
mu μ = background mortality rate or infection-related excess mortality rate (depending on model)
birth b = immigration or birth rate (same as mu if stable population); b is sometimes also used for probability of infection per contact
M = number of deceased (compartment)
T = separate compartment for infected people on treatment; if used, I becomes the compartment for untreated infected people
h = rate of transition from untreated (I) infected people to treated (T) infected people
gammaT γT = recovery rate from treated (T) infected people
betaT βT = infection rate for treated (T) infected people
V = number of vaccinated people; can be used to track vaccine-induced immunity and model leaky scenarios (compartment)
Iv Iv = number of infected vaccinated individuals (compartment)
I0 I0 I-naught = initial prevalence at time t equals 0
E = separate compartment for exposed people in the latent period between being infected and becoming infectious (infected but not yet infectious)
delta δ = rate of progression to infectious from exposed (coursework uses alpha α, other references use sigma σ)
bite (sometimes b, sometimes a) = biting rate; average number of bites per mosquito per day
bv bv βV = probability of infection from an infected host to a susceptible vector; human-to-mosquito transmission probability
bh bh βH = probability of infection from an infected vector to a susceptible host; mosquito-to-human transmission probability per bite
muv μv = mortality rate of the vector
gammah γh (sometimes r) = recovery rate of the host
VC = Vectorial capacity (coursework uses V, which is also sometimes used to denote the vaccinated compartment)
p = the probability of the vector surviving through one day (p is also sometimes used for probability of a susceptible becoming infected on contact with infected person)
n = the average duration of the extrinsic incubation period in days
Stochastic simulations of a novel pathogen
So far, we have only been working with deterministic models, where we get the same output every time we run the model with a given input. You have now been introduced to the concept of stochasticity, or random events, in modelling. Random events can be particularly important in determining the course of an epidemic at the beginning of the outbreak, when only a few people are infected.
We are modelling introduction of a novel pathogen into a completely susceptible population. The “GillespieSSA” package allows you to quickly set up a Gillespie algorithm. The focus is on understanding the behaviour of stochastic models, and how this compares to that of deterministic models!
Part 1: Comparing deterministic and stochastic model output when R0 = 0.75
The following cell contains code for a simple deterministic SIR model, which you will be very familiar with. Familiarise yourself with the input: what are the initial conditions, and the parameter values? Simulate and plot the deterministic model output; the generated plot shows the prevalence of infection.
Confirm that the infection goes extinct without causing an epidemic, as you would expect for a deterministic model with R0 < 1.
nicesubtitle <- "SIR Model v1a11 SIR Basic Model, Deterministic Model c3w1"
print(nicesubtitle)## [1] "SIR Model v1a11 SIR Basic Model, Deterministic Model c3w1"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000001 # population size
duration <- 100 # total number of days
tsteps <- 1 # chunk into 1 day intervals
beta <- 0.3 # infection rate day^-1
gamma <- 0.4 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 0.30 0.40 0.75initial_state_values <- c(S = N-1,
I = 1,
R = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)),
cumul_incid = round(N-S,digits=5)
)
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff cumul_incid
## 1 0 1e+06 1 0 100 0 0 0.7499993 1print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff cumul_incid
## 1 0 1e+06 1 0 100 0 0 0.7499993 1print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R still_Su preval_Inf propor_Re Reff
## 101 100 999997 4.539671e-05 3.999788 100 0 0 0.749997
## cumul_incid
## 101 3.99983print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff, -S, -R) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","deeppink4"),
labels=c("Prevalent infections", "Cumulative incidence")) +
scale_shape_manual(values = c(4, 20)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") 
Calculating incidence in a deterministic model Calculating the cumulative incidence for a simple SIR model - capture the amount by which susceptibility has fallen in that same time period:
paste0("Number susceptible & infected at time 0: ", round(output$S[output$time == 0] + output$I[output$time == 0],3))## [1] "Number susceptible & infected at time 0: 1000001"paste0("Number susceptible & infected at time ",max(output$time),": ", round(output$S[output$time == max(output$time)] + output$I[output$time == max(output$time)],3))## [1] "Number susceptible & infected at time 100: 999997"paste0("Difference: ", round(max(output$cumul_incid),3))## [1] "Difference: 4"The following code block performs a single stochastic simulation with the same parameters as the deterministic model, and plots the resulting prevalence of infection. Every time you run the second cell, another iteration of the model will be simulated, and the output added to the plot in a different colour. Run the second code block at least 10 times in a row, then look at the plot of the prevalence of infection.
Question: What do you observe? How does this compare to the deterministic model output?
- The deterministic model provides consistent output given fixed parameters; the number of infected persons starts at 1 and then decreases, never turning into a epidemic.
- The stochastic output appears to provide slightly different results every time it is run; the time variable t is calculated differently and shows up as non-integer values, and the resulting time span is also randomly generated. We also observe records where the stochastic number infected I exceeds the maximum value from the corresponding deterministic model.
# Simulate stochastic algorithm once:
# Defining the model and input
a <- c("beta*S*I/1000000","gamma*I")
nu <- matrix(c(-1,+1,0,0,-1,+1),
nrow=3,
ncol=2,
byrow=FALSE)
tf <- 100
# Simulation
sir_out <- ssa(initial_state_values, a, nu, parameters, tf=tf, simName="SIR")
while(sir_out$stats$nSteps==1){
sir_out <- ssa(initial_state_values,a,nu, parameters,tf=tf,simName="SIR")
}
# Record number of simulations
n_sims <- 1
# Plot
stoch_plot <- ggplot(as.data.frame(sir_out$data)) +
geom_line(aes(x = t, y = I)) +
xlab("Time (days)")+
ylab("Prevalence of infection") +
labs(title = paste("Stochastic model output for R0 =",parameters["beta"]/parameters["gamma"]),
subtitle = paste(n_sims, "simulations")) +
ylim(c(0,150)) +
xlim(c(0,100))
stoch_plot
print(paste0("Stochastic records where (I)nfected value is greater than the maximum I value (",max(output$I),") from the deterministic model:" ))## [1] "Stochastic records where (I)nfected value is greater than the maximum I value (1) from the deterministic model:"as_tibble(sir_out$data) %>% filter(I > max(output$I)) %>% arrange(t)## # A tibble: 29 x 4
## t S I R
## <dbl> <dbl> <dbl> <dbl>
## 1 2.55 999999 2 0
## 2 2.92 999998 3 0
## 3 2.98 999997 4 0
## 4 3.02 999996 5 0
## 5 3.49 999995 6 0
## 6 3.55 999995 5 1
## 7 4.27 999995 4 2
## 8 4.90 999995 3 3
## 9 5.06 999995 2 4
## 10 5.16 999994 3 4
## # ... with 19 more rowsprint("Stochastic model maximum t time value: ")## [1] "Stochastic model maximum t time value: "max(as_tibble(sir_out$data)$t)## [1] 19.62135Run this following code block at least 10 times, and examine the plot each time:
sir_out <- ssa(initial_state_values,a,nu,parameters,tf=tf,simName="SIR")
while(sir_out$stats$nSteps==1){
sir_out <- ssa(initial_state_values,a,nu,parameters,tf=tf,simName="SIR")
}
n_sims <- n_sims+1
stoch_plot <- stoch_plot +
geom_line(data = as.data.frame(sir_out$data), aes(x = t, y = I), col = sample(rainbow(20),1)) +
labs(title = paste("Stochastic model output for R0 =",parameters["beta"]/parameters["gamma"]),
subtitle = paste(n_sims, "simulations"))
stoch_plot
print(paste0("Stochastic records where (I)nfected value is greater than the maximum I value (",max(output$I),") from the deterministic model:" ))## [1] "Stochastic records where (I)nfected value is greater than the maximum I value (1) from the deterministic model:"as_tibble(sir_out$data) %>% filter(I > max(output$I)) %>% arrange(t)## # A tibble: 7 x 4
## t S I R
## <dbl> <dbl> <dbl> <dbl>
## 1 1.51 999999 2 0
## 2 2.98 999998 3 0
## 3 3.70 999997 4 0
## 4 3.87 999997 3 1
## 5 7.03 999997 2 2
## 6 7.07 999996 3 2
## 7 7.50 999996 2 3print("Stochastic model maximum t time value: ")## [1] "Stochastic model maximum t time value: "max(as_tibble(sir_out$data)$t)## [1] 10.54396Part 2: Stochastic simulations with a range of R0
Now that you are familiar with what a stochastic model prevalence output looks like for R0 < 1 in comparison to that from a deterministic model, we will explore the stochastic model behaviour for different values of R0.
We are also interested in the cumulative incidence of infection - in each iteration, how many people have become infected in total over our simulation period of 100 days? So far throughout the course, to keep things simple we have only looked at prevalence output - the number or proportion of people that are infected at a given timepoint. In many situations, we also want to know what the incidence of infection is - how many people have become infected in a given time period. In our model output, the cumulative incidence is just the total number of infections that have occurred by a given point in time. For example, if 1 person becomes infected on day 1 and 2 people become infected on day 2, then the cumulative incidence at day 2 equals 3.
We have defined a function called simulate_stoch_model() that will simulate the stochastic model. In the function call, you can choose between plotting the prevalence of infection or the cumulative incidence, for given parameter values repeatedly, rather than you having to run it manually like in the example above. If you save the function output in a variable, it will also show you the cumulative incidence by 100 days (or by the time the epidemic goes extinct, if this happens before 100 days). We also need to define a new vector of initial conditions, to save the cumulative incidence as an output. Of course, at the beginning of the simulation, the cumulative incidence is 0.
# start - source("w9_function.R")
simulate_stoch_model <- function(beta, gamma, n_sims, plot = "prevalence", final_time = 100, initial_conditions = initial_state_values) {
cum_inc <- rep(0,n_sims)
parms <- c(beta = beta,
gamma = gamma)
a <- c("beta*S*I/(S+I+R)","gamma*I")
nu <- matrix(c(-1,+1,0,+1,0,-1,+1,0),
nrow=4,
ncol=2,
byrow=FALSE)
tf <- final_time
sir_out <- ssa(initial_conditions, a, nu, parms, tf=tf, simName="SIR")
cum_inc[1] <- sir_out$data[nrow(sir_out$data),"cum_inc"]
while(sir_out$stats$nSteps==1){
sir_out <- ssa(initial_conditions, a, nu, parms, tf=tf, simName="SIR")
}
if(plot == "prevalence") {
myplot <- ggplot(as.data.frame(sir_out$data)) +
geom_line(aes(x = t, y = I), na.rm=TRUE) +
xlab("Time (days)")+
ylab("Prevalence of infection") +
labs(title = paste("Stochastic model output for R0 =",parms["beta"]/parms["gamma"]),
subtitle = paste(n_sims, "simulations")) +
xlim(c(0,final_time))
} else if(plot == "cumulative_incidence") {
myplot <- ggplot(as.data.frame(sir_out$data)) +
geom_line(aes(x = t, y = cum_inc), na.rm=TRUE) +
xlab("Time (days)")+
ylab("Cumulative incidence of infection") +
labs(title = paste("Stochastic model output for R0 =",parms["beta"]/parms["gamma"]),
subtitle = paste(n_sims, "simulations")) +
xlim(c(0,final_time))
} else {
print("plot can be either prevalence or cumulative_incidence")
}
for (i in 2:n_sims) {
sir_out <- ssa(initial_conditions, a, nu, parms, tf=tf, simName="SIR")
while(sir_out$stats$nSteps==1){
sir_out <- ssa(initial_conditions, a, nu, parms, tf=tf, simName="SIR")
}
cum_inc[i] <- sir_out$data[nrow(sir_out$data),"cum_inc"]
if(plot == "prevalence") {
myplot <- myplot +
geom_line(data = as.data.frame(sir_out$data), aes(x = t, y = I), col = sample(rainbow(20),1),
na.rm=TRUE) +
labs(title = paste("Stochastic model output for R0 =",parms["beta"]/parms["gamma"]),
subtitle = paste(n_sims, "simulations"))
} else if(plot == "cumulative_incidence") {
myplot <- myplot +
geom_line(data = as.data.frame(sir_out$data), aes(x = t, y = cum_inc), col = sample(rainbow(20),1),
na.rm=TRUE) +
labs(title = paste("Stochastic model output for R0 =",parms["beta"]/parms["gamma"]),
subtitle = paste(n_sims, "simulations"))
} else {
print("plot can be either prevalence or cumulative_incidence")
}
i <- i+1
}
print(myplot)
return(cum_inc)
}
# end - source("w9_function.R")
initial_state_values <- c(S = 1000000,
I = 1,
R = 0,
cum_inc = 0)The input arguments for the function are \(\beta\), \(\gamma\), the number of simulations and the output we want to plot. To repeat part 1 of the activity you would run the following command:
stoch_output1 <- simulate_stoch_model(beta = 0.3, gamma = 0.4, n_sims = 10, plot = "prevalence")
Note that the limits of the y axis are not pre-defined this time, so they will be automatically chosen based on the highest peak prevalence from any iteration.
Plotting the cumulative incidence:
(stoch_output1 <- simulate_stoch_model(beta = 0.3, gamma = 0.4, n_sims = 10, plot = "cumulative_incidence"))
## [1] 2 7 5 1 4 1 1 1 15 1With R0 = 0.75, you may observe the occasional epidemic with a cumulative incidence of close to 100, but in most iterations the total number of people infected over the simulation period will not exceed 10-20 people.
Keeping \(\gamma\) fixed at 0.4 per day, simulate the stochastic model for the following scenarios:
R0 = 0.1
R0 = 0.9
R0 = 1.1
For each scenario, simulate 100 iterations (it may take a little time for the output to appear).
Question: For each of the scenarios, how often do you get an epidemic?
Defining events as epidemic when the cumulative incidence is greater than 1person infected. Results based on a prior stochastic simulation:
- gamma = 0.4 R0 = 0.1 beta = 0.04 results in 17 epidemics (83 have a cumulative incidence of 0 or 1)
- gamma = 0.4 R0 = 0.9 beta = 0.36 results in 67 epidemics (33 have a cumulative incidence of 0 or 1)
- gamma = 0.4 R0 = 1.1 beta = 0.44 results in 82 epidemics (18 have a cumulative incidence of 0 or 1)
R0 <- 0.1
gamma <- 0.4
beta <- R0 * gamma
print(paste0("c3w1 gamma = ", round(gamma, 3), " R0 = ", round(R0, 3), " beta = ", round(beta, 3)))## [1] "c3w1 gamma = 0.4 R0 = 0.1 beta = 0.04"stoch_output2 <- simulate_stoch_model(beta = beta, gamma = gamma, n_sims = 100, plot = "cumulative_incidence")
print(paste0("Number of epidemics (Cumulative Incidences > 1): ", length(stoch_output2[stoch_output2>1])))## [1] "Number of epidemics (Cumulative Incidences > 1): 17"print("Lowest 2 Cumulative Incidences: ")## [1] "Lowest 2 Cumulative Incidences: "head(table(stoch_output2), 2)## stoch_output2
## 1 2
## 83 14R0 <- 0.9
gamma <- 0.4
beta <- R0 * gamma
print(paste0("c3w1 gamma = ", round(gamma, 3), " R0 = ", round(R0, 3), " beta = ", round(beta, 3)))## [1] "c3w1 gamma = 0.4 R0 = 0.9 beta = 0.36"stoch_output2 <- simulate_stoch_model(beta = beta, gamma = gamma, n_sims = 100, plot = "cumulative_incidence")
print(paste0("Number of epidemics (Cumulative Incidences > 1): ", length(stoch_output2[stoch_output2>1])))## [1] "Number of epidemics (Cumulative Incidences > 1): 67"print("Lowest 2 Cumulative Incidences: ")## [1] "Lowest 2 Cumulative Incidences: "head(table(stoch_output2), 2)## stoch_output2
## 0 1
## 1 32R0 <- 1.1
gamma <- 0.4
beta <- R0 * gamma
print(paste0("c3w1 gamma = ", round(gamma, 3), " R0 = ", round(R0, 3), " beta = ", round(beta, 3)))## [1] "c3w1 gamma = 0.4 R0 = 1.1 beta = 0.44"stoch_output2 <- simulate_stoch_model(beta = beta, gamma = gamma, n_sims = 100, plot = "cumulative_incidence")
print(paste0("Number of epidemics (Cumulative Incidences > 1): ", length(stoch_output2[stoch_output2>1])))## [1] "Number of epidemics (Cumulative Incidences > 1): 82"print("Lowest 2 Cumulative Incidences: ")## [1] "Lowest 2 Cumulative Incidences: "head(table(stoch_output2), 2)## stoch_output2
## 1 2
## 18 6Developing an age-structured model
Structure a SIR model into two age groups and define the force of infection in a population with age-specific mixing. Put this into practice to simulate an outbreak.
The structure for a model stratified into children and adults:
You will already be familiar with some assumptions underlying this structure. There are no arrows going from compartments 1 (children) to compartments 2 (adults), which means we assume there is no ageing in this population. Of course, this is only realistic if we are modelling an infection over a relatively short timespan, where ageing from childhood into adulthood doesn’t matter. Additionally, we assume that children and adults all recover from the infection at the same rate \(\gamma\).
The reason why we are stratifying the population into different age groups is because we assume age-specific mixing patterns, meaning that different ages experience different forces of infection. In our model, susceptible children experience a force of infection \(\lambda_1\) and susceptible adults become infected at a rate \(\lambda_2\).
In previous models, we have summarised 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: - \(b\), the probability of infection per contact - \(c_{ij}\), the number of contacts that a susceptible person in age group i makes with people in age group j per unit time
We are modelling an outbreak in a population based in the UK. Based on the Mossong contact study, we assume that in Great Britain, 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\)
We also assume that the probability of infection per contact, \(b\), is the same for children and adults and equals 5% for the disease we are modelling, and that people remain infected on average for 5 days. The total population size is 1 million, and children make up 20% of that population.
Based on this information, simulate an infection in a totally naïve (unexposed) population over the course of 3 months and plot the full model output.
Questions: What was the cumulative incidence of infection during this epidemic? What proportion of those infections occurred in children? Which age group was most affected by the epidemic?
- The cumulative incidence of infection is 927,447 (92.7% of the total population 1,000,000).
- From a population count perspective, adults make up 80% of the total population with a cumulative incidence of 736,937 (92% of adults).
- Children make up 20% of the total population, but 95% of children (190,511) got infected.
- Children were most affected by percentage, while adults were most affect by raw count.
nicesubtitle <- "SIR Model age2c5 Age-structured Model c3w2"
print(nicesubtitle)## [1] "SIR Model age2c5 Age-structured Model c3w2"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000 # total population size
N1 <- round(N * .2,0) # children make up 20% population
N2 <- N - N1 # adults make up the rest
duration <- 90 # total number of days - 3 months
tsteps <- .1 # chunk into 1/10th day intervals
# children-1 and adult-2 contact parameters:
c11 <- 7 # children make 13 contacts per day on average of which 7 are with other children
c12 <- 6 # children make 13 - 7 = 6 adult contacts per day on average
c22 <- 10 # adults make 11 contacts per day on average of which 10 are with other adults
c21 <- 1 # adults make 11 - 1 = 10 contacts with children per day on average
b <- .05 # probability of infection per contact equals 5%, same for children and adults (note: prior exercises used the variable name p for for vaccination rate)
gamma <- 5 ** -1 # recovery rate - people remain infected on average for 5 days
# beta = infection rate day^-1 is a function of b * c)
# with 4 values of c, there will be 4 betas and thus 4 R0s
beta11 <- b * c11 # children with children
beta12 <- b * c12 # children with adults
beta21 <- b * c21 # adults with children
beta22 <- b * c22 # adults with adults
R011 <- b * c11 / gamma # children with children
R012 <- b * c12 / gamma # children with adults
R021 <- b * c21 / gamma # adults with children
R022 <- b * c22 / gamma # adults with adults
(parameters <- c(
gamma = gamma, # recovery rate
b = b,
c11 = c11, # infection rate - children w children
c12 = c12, # infection rate - children w adults
c21 = c21, # infection rate - adults w children
c22 = c22, # infection rate - adults w adults
R011 = R011, # basic repro number - children w children
R012 = R012, # basic repro number - children w adults
R021 = R021, # basic repro number - adults w children
R022 = R022 # basic repro number - adults w adults
)) ## gamma b c11 c12 c21 c22 R011 R012 R021 R022
## 0.20 0.05 7.00 6.00 1.00 10.00 1.75 1.50 0.25 2.50initial_state_values <- c(S1 = N1-1,
I1 = 1,
R1 = 0,
S2 = N2,
I2 = 0,
R2 = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_age2 <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N1 <- S1 + I1 + R1
N2 <- S2 + I2 + R2
lambda1 <- (b * c11 * I1/N1) + (b * c12 * I2/N2)
dS1 <- -(lambda1 * S1)
dI1 <- (lambda1 * S1) -(gamma * I1)
dR1 <- (gamma * I1)
lambda2 <- (b * c21 * I1/N1) + (b * c22 * I2/N2)
dS2 <- -(lambda2 * S2)
dI2 <- (lambda2 * S2) -(gamma * I2)
dR2 <- (gamma * I2)
return(list(c(dS1, dI1, dR1, dS2, dI2, dR2)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_age2,
parms = parameters)) %>%
mutate(
still_Su1 = round(S1/(S1+I1+R1),digits=5)*100,
preval_Inf1 = round(I1/(S1+I1+R1),digits=5)*100,
propor_Re1 = round(R1/(S1+I1+R1),digits=5)*100,
cumul_incid1 = round(N1-S1,digits=5)
, cumul_incid1_pct = round(cumul_incid1/N1,digits=5)*100
,
still_Su2 = round(S2/(S2+I2+R2),digits=5)*100,
preval_Inf2 = round(I2/(S2+I2+R2),digits=5)*100,
propor_Re2 = round(R2/(S2+I2+R2),digits=5)*100,
cumul_incid2 = round(N2-S2,digits=5)
, cumul_incid2_pct = round(cumul_incid2/N2,digits=5)*100
,
still_Su = round((S1+S2)/N,digits=5)*100,
preval_Inf = round((I1+I2)/N,digits=5)*100,
propor_Re = round((R1+R2)/N,digits=5)*100,
cumul_incid = round(N-(S1+S2),digits=5)
, cumul_incid_pct = round(cumul_incid/N,digits=5)*100
)
print("Initial record: ")## [1] "Initial record: "output %>%
arrange(time) %>%
head(1)## time S1 I1 R1 S2 I2 R2 still_Su1 preval_Inf1 propor_Re1 cumul_incid1
## 1 0 199999 1 0 8e+05 0 0 99.999 0 0 1
## cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 cumul_incid2
## 1 0 100 0 0 0
## cumul_incid2_pct still_Su preval_Inf propor_Re cumul_incid cumul_incid_pct
## 1 0 100 0 0 1 0print("peak infection day when child I1 is at its max: ")## [1] "peak infection day when child I1 is at its max: "output %>%
arrange(-I1, time) %>%
head(1)## time S1 I1 R1 S2 I2 R2 still_Su1
## 1 38.8 63933.86 59547.8 76518.34 331651.3 217187.2 251161.5 31.967
## preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2
## 1 29.774 38.259 136066.1 68.033 41.456 27.148
## propor_Re2 cumul_incid2 cumul_incid2_pct still_Su preval_Inf propor_Re
## 1 31.395 468348.7 58.544 39.559 27.673 32.768
## cumul_incid cumul_incid_pct
## 1 604414.8 60.441print("peak infection day when adult I2 is at its max: ")## [1] "peak infection day when adult I2 is at its max: "output %>%
arrange(-I2, time) %>%
head(1)## time S1 I1 R1 S2 I2 R2 still_Su1 preval_Inf1
## 1 39.7 54084.74 58729.24 87186.02 289307.9 220088 290604 27.042 29.365
## propor_Re1 cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2
## 1 43.593 145915.3 72.958 36.163 27.511 36.326
## cumul_incid2 cumul_incid2_pct still_Su preval_Inf propor_Re cumul_incid
## 1 510692.1 63.837 34.339 27.882 37.779 656607.3
## cumul_incid_pct
## 1 65.661print("peak infection day when all preval_Inf (children + adults) is at its max: ")## [1] "peak infection day when all preval_Inf (children + adults) is at its max: "output %>%
arrange(-preval_Inf, time) %>%
head(1)## time S1 I1 R1 S2 I2 R2 still_Su1
## 1 39.5 56129.79 59039.8 84830.41 298252.7 219944.9 281802.4 28.065
## preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2
## 1 29.52 42.415 143870.2 71.935 37.282 27.493
## propor_Re2 cumul_incid2 cumul_incid2_pct still_Su preval_Inf propor_Re
## 1 35.225 501747.3 62.718 35.438 27.898 36.663
## cumul_incid cumul_incid_pct
## 1 645617.5 64.562print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S1 I1 R1 S2 I2 R2 still_Su1
## 901 90 9489.424 29.33085 190481.2 63063.21 177.5459 736759.2 4.745
## preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2
## 901 0.015 95.241 190510.6 95.255 7.883 0.022
## propor_Re2 cumul_incid2 cumul_incid2_pct still_Su preval_Inf propor_Re
## 901 92.095 736936.8 92.117 7.255 0.021 92.724
## cumul_incid cumul_incid_pct
## 901 927447.4 92.745print("Total cumulative incidence:")## [1] "Total cumulative incidence:"round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]),0)## [1] 927447round(
((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]))/N,3)*100## [1] 92.7print("Cumulative incidence in children:")## [1] "Cumulative incidence in children:"round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]),0)## [1] 190511round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) / (output$S1[1]+output$I1[1]),3) * 100## [1] 95.3print("Proportion of total infections who are children:")## [1] "Proportion of total infections who are children:"round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) / ((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)])),3) * 100## [1] 20.5print("Cumulative incidence in adults:")## [1] "Cumulative incidence in adults:"round((output$S2[1]+output$I2[1]-output$S2[nrow(output)]),0)## [1] 736937round((output$S2[1]+output$I2[1]-output$S2[nrow(output)]) / (output$S2[1]+output$I2[1]),3) * 100## [1] 92.1print("Proportion of total infections who are adults:")## [1] "Proportion of total infections who are adults:"round((output$S2[1]+output$I2[1]-output$S2[nrow(output)]) / ((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)])),3) * 100## [1] 79.5# part 2 - plots
print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"output %>%
select(time, preval_Inf, preval_Inf1, preval_Inf2) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dashed") +
geom_point(show.legend = FALSE) +
scale_color_manual(values = c("red","blue3","chartreuse"),
labels=c("ALL", "Children","Adults")) +
scale_shape_manual(values = c(4, "o", 3)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(
" gamma ", parameters['gamma'],
" R011 ", round(parameters['R011'],3),
" R012 ", round(parameters['R012'],3),
" R022 ", round(parameters['R022'],3),
" R021 ", round(parameters['R021'],3)
),
color = "Prevalent infections",
subtitle= nicesubtitle) +
theme(legend.position="bottom") 
# part 3 - plots continued
print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"output %>%
select(time, cumul_incid, cumul_incid1, cumul_incid2) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dashed") +
geom_point(show.legend = FALSE) +
scale_color_manual(values = c("red","blue3","chartreuse"),
labels=c("ALL", "Children","Adults")) +
scale_shape_manual(values = c(4, "o", 3)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(
" gamma ", parameters['gamma'],
" R011 ", round(parameters['R011'],3),
" R012 ", round(parameters['R012'],3),
" R022 ", round(parameters['R022'],3),
" R021 ", round(parameters['R021'],3)
),
color = "Cumulative incidence",
subtitle= nicesubtitle) +
theme(legend.position="bottom") 
# part 4 - plots continued
print("Plotting the proprotion of people in each compartment over time")## [1] "Plotting the proprotion of people in each compartment over time"output %>%
select(time, cumul_incid_pct, cumul_incid1_pct, cumul_incid2_pct) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dashed") +
geom_point(show.legend = FALSE) +
scale_color_manual(values = c("red","blue3","chartreuse"),
labels=c("ALL", "Children","Adults")) +
scale_shape_manual(values = c(4, "o", 3)) +
xlab("Time (days)") +
ylab("% of Population") +
labs(title=paste(
" gamma ", parameters['gamma'],
" R011 ", round(parameters['R011'],3),
" R012 ", round(parameters['R012'],3),
" R022 ", round(parameters['R022'],3),
" R021 ", round(parameters['R021'],3)
),
color = "Cumulative incidence %",
subtitle= nicesubtitle) +
theme(legend.position="bottom") 
Additional information: realistic contact patterns
We have used empirical data from the Mossong study to assign values to our contact parameters.
Think more generally about what makes a realistic contact pattern in human populations. It is logical that within a population, the total number of contacts made by children with adults has to be consistent with the total number of contacts made by adults with children. We can check if this is the case in our example by calculating if our contact parameters fulfill this condition:
Average daily number of contacts made by children with adults (\(c_{12}\)) * proportion of children in the population = average daily number of contacts made by adults with children (\(c_{21}\)) * proportion of adults in the population
\[\begin{align} 6 \times 0.2 = 1.2 \\ 1 \times 0.8 = 0.8 \end{align}\]
… which is approximately the same, given that we have rounded the average number of contacts to whole numbers for simplicity.
Extending your model to 3 age groups
Develop the simple age-structured model further by extending the model structure to 3 age groups - to represent children, adults and the elderly in separate sets of SIR compartments. We are modelling the same infection and population as in the previous activity, with the only difference being this additional distinction between adults and the elderly. These 3 age groups experience different forces of infection due to age-specific mixing.
Here is the input information for your model: - The probability of infection per contact, \(b\), is 5% for people of all ages - The total population size is 1 million, and children make up 20% of this population. 15% of the population are elderly. - The population has not been exposed to the infection before.
The contacts between our 3 age groups can be summarised in a contact matrix, as follows: \[\begin{pmatrix} c_{11} & c_{12} & c_{13} \\ c_{21} & c_{22} & c_{23} \\ c_{31} & c_{32} & c_{33} \end{pmatrix}\]where 1 represents children, 2 represents adults and 3 represents the elderly.
The average daily number of contacts between the different age groups in the contact matrix are: \[\begin{pmatrix} 7 & 5 & 1 \\ 2 & 9 & 1 \\ 1 & 3 & 2 \end{pmatrix}\]Question: How many infections occurred in each age group over the course of the epidemic, and what proportion of children, adults and elderly individuals does this represent?
Total incidence after 90 days (% of category):
- ALL: 908,622 (90.9%)
- Children: 190,215 (95.1%)
- Adults: 609,095 (93.7%)
- Elderly: 109,311 (72.9%)
# part 1 initial attempt
nicesubtitle <- "SIR Model age3c5 Age-structured Model c3w2"
print(nicesubtitle)## [1] "SIR Model age3c5 Age-structured Model c3w2"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000 # total population size
N1 <- round(N * .2,0) # children make up 20% of the population
N3 <- round(N * .15,0) # elderly make up 15% of the population
N2 <- N - N1 - N3 # adults make up the rest
duration <- 90 # total number of days - 3 months
tsteps <- 0.1 # chunk into 1/10th day intervals
# children-1 adult-2 elderly-3 contact parameters:
c11 <- 7 # children with other children
c12 <- 5 # children with adults
c13 <- 1 # children with elderly
c21 <- 2 # adults with children
c22 <- 9 # adults with other adults
c23 <- 1 # adults with elderly
c31 <- 1 # elderly with children
c32 <- 3 # elderly with adults
c33 <- 2 # elderly with other elderly
b <- .05 # probability of infection per contact equals 5%, same for children and adults (note: prior exercises used the variable name p for for vaccination rate)
gamma <- 5 ** -1 # recovery rate - people remain infected on average for 5 days
# beta = infection rate day^-1 is a function of b * c)
# with 9 values of c, there will be 9 betas and thus 9 R0s
beta11 <- b * c11 # children with children
beta12 <- b * c12 # children with adults
beta13 <- b * c13 # children with elderly
beta21 <- b * c21 # adults with children
beta22 <- b * c22 # adults with adults
beta23 <- b * c23 # adults with elderly
beta31 <- b * c31 # elderly with children
beta32 <- b * c32 # elderly with adults
beta33 <- b * c33 # elderly with elderly
R011 <- beta11 / gamma
R012 <- beta12 / gamma
R013 <- beta13 / gamma
R021 <- beta21 / gamma
R022 <- beta22 / gamma
R023 <- beta23 / gamma
R031 <- beta31 / gamma
R032 <- beta32 / gamma
R033 <- beta33 / gamma
(parameters <- c(
gamma = gamma, # recovery rate
beta11 = beta11, # infection rate - children w children
beta12 = beta12, # infection rate - children w adults
beta13 = beta13, # infection rate - children w elderly
beta21 = beta21, # infection rate - adults w children
beta22 = beta22, # infection rate - adults w adults
beta23 = beta23, # infection rate - adults w elderly
beta31 = beta31, # infection rate - elderly w children
beta32 = beta32, # infection rate - elderly w adults
beta33 = beta33, # infection rate - elderly w elderly
R011 = R011, # basic repro number - children w children
R012 = R012, # basic repro number - children w adults
R013 = R013, # basic repro number - children w elderly
R021 = R021, # basic repro number - adults w children
R022 = R022, # basic repro number - adults w adults
R023 = R023, # basic repro number - adults w elderly
R031 = R031, # basic repro number - elderly w children
R032 = R032, # basic repro number - elderly w adults
R033 = R033 # basic repro number - elderly w elderly
)) ## gamma beta11 beta12 beta13 beta21 beta22 beta23 beta31 beta32 beta33 R011
## 0.20 0.35 0.25 0.05 0.10 0.45 0.05 0.05 0.15 0.10 1.75
## R012 R013 R021 R022 R023 R031 R032 R033
## 1.25 0.25 0.50 2.25 0.25 0.25 0.75 0.50initial_state_values <- c(S1 = N1-1,
I1 = 1,
R1 = 0,
S2 = N2,
I2 = 0,
R2 = 0,
S3 = N3,
I3 = 0,
R3 = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_age3 <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S1 + I1 + R1 + S2 + I2 + R2 + S3 + I3 + R3
N1 <- S1 + I1 + R1
N2 <- S2 + I2 + R2
N3 <- S3 + I3 + R3
lambda1 <- (beta11 * I1/N1) + (beta12 * I2/N2) + (beta13 * I3/N3)
dS1 <- -(lambda1 * S1)
dI1 <- (lambda1 * S1) -(gamma * I1)
dR1 <- (gamma * I1)
lambda2 <- (beta21 * I1/N1) + (beta22 * I2/N2) + (beta23 * I3/N3)
dS2 <- -(lambda2 * S2)
dI2 <- (lambda2 * S2) -(gamma * I2)
dR2 <- (gamma * I2)
lambda3 <- (beta31 * I1/N1) + (beta32 * I2/N2) + (beta33 * I3/N3)
dS3 <- -(lambda3 * S3)
dI3 <- (lambda3 * S3) -(gamma * I3)
dR3 <- (gamma * I3)
return(list(c(dS1, dI1, dR1, dS2, dI2, dR2, dS3, dI3, dR3)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_age3,
parms = parameters)) %>%
mutate(
still_Su1 = round(S1/(S1+I1+R1),digits=5)*100,
preval_Inf1 = round(I1/(S1+I1+R1),digits=5)*100,
propor_Re1 = round(R1/(S1+I1+R1),digits=5)*100,
cumul_incid1 = round(N1-S1,digits=5)
, cumul_incid1_pct = round(cumul_incid1/N1,digits=5)*100
,
still_Su2 = round(S2/(S2+I2+R2),digits=5)*100,
preval_Inf2 = round(I2/(S2+I2+R2),digits=5)*100,
propor_Re2 = round(R2/(S2+I2+R2),digits=5)*100,
cumul_incid2 = round(N2-S2,digits=5)
, cumul_incid2_pct = round(cumul_incid2/N2,digits=5)*100
,
still_Su3 = round(S3/(S3+I3+R3),digits=5)*100,
preval_Inf3 = round(I3/(S3+I3+R3),digits=5)*100,
propor_Re3 = round(R3/(S3+I3+R3),digits=5)*100,
cumul_incid3 = round(N3-S3,digits=5)
, cumul_incid3_pct = round(cumul_incid3/N3,digits=5)*100
,
still_Su = round((S1+S2+S3)/N,digits=5)*100,
preval_Inf = round((I1+I2+I3)/N,digits=5)*100,
propor_Re = round((R1+R2+R3)/N,digits=5)*100,
cumul_incid = round(N-(S1+S2+S3),digits=5)
, cumul_incid_pct = round(cumul_incid/N,digits=5)*100
)
print("Initial record: ")## [1] "Initial record: "output %>%
arrange(time) %>%
head(1)## time S1 I1 R1 S2 I2 R2 S3 I3 R3 still_Su1 preval_Inf1 propor_Re1
## 1 0 199999 1 0 650000 0 0 150000 0 0 99.999 0 0
## cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 cumul_incid2
## 1 1 0 100 0 0 0
## cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 cumul_incid3
## 1 0 100 0 0 0
## cumul_incid3_pct still_Su preval_Inf propor_Re cumul_incid cumul_incid_pct
## 1 0 100 0 0 1 0print("peak infection day when child I1 is at its max: ")## [1] "peak infection day when child I1 is at its max: "output %>%
arrange(-I1, time) %>%
head(1)## time S1 I1 R1 S2 I2 R2 S3 I3
## 1 36.3 65139.01 60421.35 74439.64 239663.2 188559 221777.7 96911.26 26720.94
## R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct
## 1 26367.8 32.57 30.211 37.22 134861 67.43
## still_Su2 preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3
## 1 36.871 29.009 34.12 410336.8 63.129 64.608
## preval_Inf3 propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 17.814 17.579 53088.74 35.392 40.171 27.57
## propor_Re cumul_incid cumul_incid_pct
## 1 32.259 598286.5 59.829print("peak infection day when adult I2 is at its max: ")## [1] "peak infection day when adult I2 is at its max: "output %>%
arrange(-I2, time) %>%
head(1)## time S1 I1 R1 S2 I2 R2 S3 I3
## 1 36.7 60434.87 60294.62 79270.51 223902.9 189201.8 236895.3 93980.17 27483.03
## R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct
## 1 28536.81 30.217 30.147 39.635 139565.1 69.783
## still_Su2 preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3
## 1 34.447 29.108 36.445 426097.1 65.553 62.653
## preval_Inf3 propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 18.322 19.025 56019.83 37.347 37.832 27.698
## propor_Re cumul_incid cumul_incid_pct
## 1 34.47 621682 62.168print("peak infection day when elderly I3 is at its max: ")## [1] "peak infection day when elderly I3 is at its max: "output %>%
arrange(-I3, time) %>%
head(1)## time S1 I1 R1 S2 I2 R2 S3 I3
## 1 38.9 40463.97 54707.3 104828.7 155189.3 176566.9 318243.8 79448.69 29366.22
## R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct
## 1 41185.09 20.232 27.354 52.414 159536 79.768
## still_Su2 preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3
## 1 23.875 27.164 48.961 494810.7 76.125 52.966
## preval_Inf3 propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 19.577 27.457 70551.31 47.034 27.51 26.064
## propor_Re cumul_incid cumul_incid_pct
## 1 46.426 724898 72.49print("peak infection day when all preval_Inf (children + adults + elderly) is at its max: ")## [1] "peak infection day when all preval_Inf (children + adults + elderly) is at its max: "output %>%
arrange(-preval_Inf, time) %>%
head(1)## time S1 I1 R1 S2 I2 R2 S3 I3
## 1 36.8 59312.82 60211.57 80475.61 220122 189198.6 240679.4 93258.31 27653.5
## R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct
## 1 29088.19 29.656 30.106 40.238 140687.2 70.344
## still_Su2 preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3
## 1 33.865 29.107 37.028 429878 66.135 62.172
## preval_Inf3 propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 18.436 19.392 56741.69 37.828 37.269 27.706
## propor_Re cumul_incid cumul_incid_pct
## 1 35.024 627306.8 62.731print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S1 I1 R1 S2 I2 R2 S3
## 901 90 9784.579 20.14349 190195.3 40904.55 80.02684 609015.4 40688.93
## I3 R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1
## 901 43.53274 109267.5 4.892 0.01 95.098 190215.4
## cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 cumul_incid2
## 901 95.108 6.293 0.012 93.695 609095.5
## cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 cumul_incid3
## 901 93.707 27.126 0.029 72.845 109311.1
## cumul_incid3_pct still_Su preval_Inf propor_Re cumul_incid cumul_incid_pct
## 901 72.874 9.138 0.014 90.848 908621.9 90.862print("Total cumulative incidence:")## [1] "Total cumulative incidence:"round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)]),0)## [1] 908622round(
((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)]))/N,3)*100## [1] 90.9print("Cumulative incidence in children:")## [1] "Cumulative incidence in children:"round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]),0)## [1] 190215round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) / (output$S1[1]+output$I1[1]),3) * 100## [1] 95.1print("Proportion of total infections who are children:")## [1] "Proportion of total infections who are children:"round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) / ((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)])),3) * 100## [1] 20.9print("Cumulative incidence in adults:")## [1] "Cumulative incidence in adults:"round((output$S2[1]+output$I2[1]-output$S2[nrow(output)]),0)## [1] 609095round((output$S2[1]+output$I2[1]-output$S2[nrow(output)]) / (output$S2[1]+output$I2[1]),3) * 100## [1] 93.7print("Proportion of total infections who are adults:")## [1] "Proportion of total infections who are adults:"round((output$S2[1]+output$I2[1]-output$S2[nrow(output)]) / ((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)])),3) * 100## [1] 67print("Cumulative incidence in elderly:")## [1] "Cumulative incidence in elderly:"round((output$S3[1]+output$I3[1]-output$S3[nrow(output)]),0)## [1] 109311round((output$S3[1]+output$I3[1]-output$S3[nrow(output)]) / (output$S3[1]+output$I3[1]),3) * 100## [1] 72.9print("Proportion of total infections who are elderly:")## [1] "Proportion of total infections who are elderly:"round((output$S3[1]+output$I3[1]-output$S3[nrow(output)]) / ((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)])),3) * 100## [1] 12Alternative solution: a scalable approach
The previous approach worked well for an age-structured model because we only have a few age groups and is easy to understand. If we have more than 3 age groups, writing out all compartments and parameters separately quickly becomes repetitive and even unfeasible.
Define S, I and R compartments as vectors within the model function, each of length 3 for the 3 age groups. This allows us to write out the differential equations only once for each epidemiological compartment (S, I, R). (We have to reorder the initial states vector.)
Define the age-specific contact pattern as a matrix. When calculating the force of infection, we multiply the contact matrix with the vector of the proportion of infected people for each age group, giving us 3 age-specific forces of infection stored in the “lambda” vector.
# part 1b subsequent attempt using matrix per provided example
nicesubtitle <- "SIR Model age3c5m Age-structured Model c3w2"
print(nicesubtitle)## [1] "SIR Model age3c5m Age-structured Model c3w2"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000 # total population size
N1 <- round(N * .2,0) # children make up 20% of the population
N3 <- round(N * .15,0) # elderly make up 15% of the population
N2 <- N - N1 - N3 # adults make up the rest
duration <- 90 # total number of days - 3 months
tsteps <- 0.1 # chunk into 1/10th day intervals
# children-1 adult-2 elderly-3 contact parameters:
# Set up an empty contact matrix with rows for each age group and columns for each age group
contact_matrix <- matrix(0,nrow=3,ncol=3)
# Fill in the contract matrix
contact_matrix[1,1] = 7 # children with other children
contact_matrix[1,2] = 5 # children with adults
contact_matrix[1,3] = 1 # children with elderly
contact_matrix[2,1] = 2 # adults with children
contact_matrix[2,2] = 9 # adults with other adults
contact_matrix[2,3] = 1 # adults with elderly
contact_matrix[3,1] = 1 # elderly with children
contact_matrix[3,2] = 3 # elderly with adults
contact_matrix[3,3] = 2 # elderly with other elderly
b <- .05 # probability of infection per contact equals 5%, same for children and adults (note: prior exercises used the variable name p for for vaccination rate)
gamma <- 5 ** -1 # recovery rate - people remain infected on average for 5 days
# beta = infection rate day^-1 is a function of b * c)
# with n values of c, there will be n betas and thus n R0s
# initialize and fill beta matrix
beta_matrix <- matrix(0,nrow=3,ncol=3)
beta_matrix <- b * contact_matrix
# initialize and fill R0 matrix
R0_matrix <- matrix(0,nrow=3,ncol=3)
R0_matrix <- b * contact_matrix / gamma
(parameters <- c(
gamma = gamma, # recovery rate
b = b, # probability of infection per contact
contact_matrix = contact_matrix, # ave daily contacts
beta_matrix = beta_matrix, # infection rate
R0_matrix = R0_matrix # basic repro number
)) ## gamma b contact_matrix1 contact_matrix2 contact_matrix3
## 0.20 0.05 7.00 2.00 1.00
## contact_matrix4 contact_matrix5 contact_matrix6 contact_matrix7 contact_matrix8
## 5.00 9.00 3.00 1.00 1.00
## contact_matrix9 beta_matrix1 beta_matrix2 beta_matrix3 beta_matrix4
## 2.00 0.35 0.10 0.05 0.25
## beta_matrix5 beta_matrix6 beta_matrix7 beta_matrix8 beta_matrix9
## 0.45 0.15 0.05 0.05 0.10
## R0_matrix1 R0_matrix2 R0_matrix3 R0_matrix4 R0_matrix5
## 1.75 0.50 0.25 1.25 2.25
## R0_matrix6 R0_matrix7 R0_matrix8 R0_matrix9
## 0.75 0.25 0.25 0.50# note different order to initial state values to work with matrix
initial_state_values <- c(S1 = N1-1,
S2 = N2,
S3 = N3,
I1 = 1,
I2 = 0,
I3 = 0,
R1 = 0,
R2 = 0,
R3 = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_age3m <- function(time, state, parameters) {
with(as.list(parameters), { # referenced differently
n_agegroups <- 3 # number of age groups
S <- state[1:n_agegroups] # first 3 numbers in the initial_state_values vector
I <- state[(n_agegroups+1):(2*n_agegroups)] # numbers 4 to 6 in the initial_state_values vector
R <- state[(2*n_agegroups+1):(3*n_agegroups)] # numbers 7 to 9 in the initial_state_values vector
N <- S+I+R # people in S, I and R are added separately by age group, so N is also a vector of length 3
# Force of infection acting on susceptible children
lambda <- b * contact_matrix %*% as.matrix(I/N)
# lambda <- beta_matrix %*% as.matrix(I/N)
# %*% is used to multiply matrices in R
# the lambda vector contains the forces of infection for children, adults and the elderly (length 3)
# The differential equations
dS <- -(lambda * S)
dI <- (lambda * S) - (gamma * I)
dR <- (gamma * I)
# Output
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_age3m,
parms = parameters)) %>%
mutate(
still_Su1 = round(S1/(S1+I1+R1),digits=5)*100,
preval_Inf1 = round(I1/(S1+I1+R1),digits=5)*100,
propor_Re1 = round(R1/(S1+I1+R1),digits=5)*100,
cumul_incid1 = round(N1-S1,digits=5)
, cumul_incid1_pct = round(cumul_incid1/N1,digits=5)*100
,
still_Su2 = round(S2/(S2+I2+R2),digits=5)*100,
preval_Inf2 = round(I2/(S2+I2+R2),digits=5)*100,
propor_Re2 = round(R2/(S2+I2+R2),digits=5)*100,
cumul_incid2 = round(N2-S2,digits=5)
, cumul_incid2_pct = round(cumul_incid2/N2,digits=5)*100
,
still_Su3 = round(S3/(S3+I3+R3),digits=5)*100,
preval_Inf3 = round(I3/(S3+I3+R3),digits=5)*100,
propor_Re3 = round(R3/(S3+I3+R3),digits=5)*100,
cumul_incid3 = round(N3-S3,digits=5)
, cumul_incid3_pct = round(cumul_incid3/N3,digits=5)*100
,
still_Su = round((S1+S2+S3)/N,digits=5)*100,
preval_Inf = round((I1+I2+I3)/N,digits=5)*100,
propor_Re = round((R1+R2+R3)/N,digits=5)*100,
cumul_incid = round(N-(S1+S2+S3),digits=5)
, cumul_incid_pct = round(cumul_incid/N,digits=5)*100
)
print("Initial record: ")## [1] "Initial record: "output %>%
arrange(time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2 R3 still_Su1 preval_Inf1 propor_Re1
## 1 0 199999 650000 150000 1 0 0 0 0 0 99.999 0 0
## cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 cumul_incid2
## 1 1 0 100 0 0 0
## cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 cumul_incid3
## 1 0 100 0 0 0
## cumul_incid3_pct still_Su preval_Inf propor_Re cumul_incid cumul_incid_pct
## 1 0 100 0 0 1 0print("peak infection day when child I1 is at its max: ")## [1] "peak infection day when child I1 is at its max: "output %>%
arrange(-I1, time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2
## 1 36.3 65139.01 239663.2 96911.26 60421.35 188559 26720.94 74439.64 221777.7
## R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct
## 1 26367.8 32.57 30.211 37.22 134861 67.43
## still_Su2 preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3
## 1 36.871 29.009 34.12 410336.8 63.129 64.608
## preval_Inf3 propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 17.814 17.579 53088.74 35.392 40.171 27.57
## propor_Re cumul_incid cumul_incid_pct
## 1 32.259 598286.5 59.829print("peak infection day when adult I2 is at its max: ")## [1] "peak infection day when adult I2 is at its max: "output %>%
arrange(-I2, time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2
## 1 36.7 60434.87 223902.9 93980.17 60294.62 189201.8 27483.03 79270.51 236895.3
## R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct
## 1 28536.81 30.217 30.147 39.635 139565.1 69.783
## still_Su2 preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3
## 1 34.447 29.108 36.445 426097.1 65.553 62.653
## preval_Inf3 propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 18.322 19.025 56019.83 37.347 37.832 27.698
## propor_Re cumul_incid cumul_incid_pct
## 1 34.47 621682 62.168print("peak infection day when elderly I3 is at its max: ")## [1] "peak infection day when elderly I3 is at its max: "output %>%
arrange(-I3, time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2
## 1 38.9 40463.97 155189.3 79448.69 54707.3 176566.9 29366.22 104828.7 318243.8
## R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct
## 1 41185.09 20.232 27.354 52.414 159536 79.768
## still_Su2 preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3
## 1 23.875 27.164 48.961 494810.7 76.125 52.966
## preval_Inf3 propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 19.577 27.457 70551.31 47.034 27.51 26.064
## propor_Re cumul_incid cumul_incid_pct
## 1 46.426 724898 72.49print("peak infection day when all preval_Inf (children + adults + elderly) is at its max: ")## [1] "peak infection day when all preval_Inf (children + adults + elderly) is at its max: "output %>%
arrange(-preval_Inf, time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2
## 1 36.8 59312.82 220122 93258.31 60211.57 189198.6 27653.5 80475.61 240679.4
## R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct
## 1 29088.19 29.656 30.106 40.238 140687.2 70.344
## still_Su2 preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3
## 1 33.865 29.107 37.028 429878 66.135 62.172
## preval_Inf3 propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 18.436 19.392 56741.69 37.828 37.269 27.706
## propor_Re cumul_incid cumul_incid_pct
## 1 35.024 627306.8 62.731print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S1 S2 S3 I1 I2 I3 R1
## 901 90 9784.579 40904.55 40688.93 20.14349 80.02684 43.53274 190195.3
## R2 R3 still_Su1 preval_Inf1 propor_Re1 cumul_incid1
## 901 609015.4 109267.5 4.892 0.01 95.098 190215.4
## cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 cumul_incid2
## 901 95.108 6.293 0.012 93.695 609095.5
## cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 cumul_incid3
## 901 93.707 27.126 0.029 72.845 109311.1
## cumul_incid3_pct still_Su preval_Inf propor_Re cumul_incid cumul_incid_pct
## 901 72.874 9.138 0.014 90.848 908621.9 90.862print("Total cumulative incidence:")## [1] "Total cumulative incidence:"round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)]),0)## [1] 908622round(
((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)]))/N,3)*100## [1] 90.9print("Cumulative incidence in children:")## [1] "Cumulative incidence in children:"round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]),0)## [1] 190215round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) / (output$S1[1]+output$I1[1]),3) * 100## [1] 95.1print("Proportion of total infections who are children:")## [1] "Proportion of total infections who are children:"round((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) / ((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)])),3) * 100## [1] 20.9print("Cumulative incidence in adults:")## [1] "Cumulative incidence in adults:"round((output$S2[1]+output$I2[1]-output$S2[nrow(output)]),0)## [1] 609095round((output$S2[1]+output$I2[1]-output$S2[nrow(output)]) / (output$S2[1]+output$I2[1]),3) * 100## [1] 93.7print("Proportion of total infections who are adults:")## [1] "Proportion of total infections who are adults:"round((output$S2[1]+output$I2[1]-output$S2[nrow(output)]) / ((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)])),3) * 100## [1] 67print("Cumulative incidence in elderly:")## [1] "Cumulative incidence in elderly:"round((output$S3[1]+output$I3[1]-output$S3[nrow(output)]),0)## [1] 109311round((output$S3[1]+output$I3[1]-output$S3[nrow(output)]) / (output$S3[1]+output$I3[1]),3) * 100## [1] 72.9print("Proportion of total infections who are elderly:")## [1] "Proportion of total infections who are elderly:"round((output$S3[1]+output$I3[1]-output$S3[nrow(output)]) / ((output$S1[1]+output$I1[1]-output$S1[nrow(output)]) +(output$S2[1]+output$I2[1]-output$S2[nrow(output)]) +(output$S3[1]+output$I3[1]-output$S3[nrow(output)])),3) * 100## [1] 12# part 2 - plots
print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"output %>%
select(time, preval_Inf, preval_Inf1, preval_Inf2, preval_Inf3) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dashed") +
geom_point(show.legend = FALSE) +
scale_color_manual(values = c("red","blue3","chartreuse", "bisque4"),
labels=c("ALL", "Children","Adults","Elderly")) +
scale_shape_manual(values = c(4, "o", 3, "-")) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(
" gamma ", parameters['gamma'],
" R011 ", round(parameters['R0_matrix1'],3),
" R012 ", round(parameters['R0_matrix4'],3),
" R013 ", round(parameters['R0_matrix7'],3),
"\n",
" R021 ", round(parameters['R0_matrix2'],3),
" R022 ", round(parameters['R0_matrix5'],3),
" R023 ", round(parameters['R0_matrix8'],3),
"\n",
" R031 ", round(parameters['R0_matrix3'],3),
" R032 ", round(parameters['R0_matrix6'],3),
" R033 ", round(parameters['R0_matrix9'],3)
),
color = "Prevalent infections",
subtitle= nicesubtitle) +
# labs(title=paste(
# " gamma ", parameters['gamma'],
# " R011 ", round(parameters['R011'],3),
# " R012 ", round(parameters['R012'],3),
# " R013 ", round(parameters['R013'],3),"\n",
# " R021 ", round(parameters['R021'],3),
# " R022 ", round(parameters['R022'],3),
# " R023 ", round(parameters['R023'],3),
# "\n",
# " R031 ", round(parameters['R031'],3),
# " R032 ", round(parameters['R032'],3),
# " R033 ", round(parameters['R033'],3)
# ),
# color = "Cumulative incidence",
# subtitle= nicesubtitle) +
theme(legend.position="bottom") 
# part 3 - plots continued
print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"output %>%
select(time, cumul_incid, cumul_incid1, cumul_incid2, cumul_incid3) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dashed") +
geom_point(show.legend = FALSE) +
scale_color_manual(values = c("red","blue3","chartreuse", "bisque4"),
labels=c("ALL", "Children","Adults","Elderly")) +
scale_shape_manual(values = c(4, "o", 3, "-")) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(
" gamma ", parameters['gamma'],
" R011 ", round(parameters['R0_matrix1'],3),
" R012 ", round(parameters['R0_matrix4'],3),
" R013 ", round(parameters['R0_matrix7'],3),
"\n",
" R021 ", round(parameters['R0_matrix2'],3),
" R022 ", round(parameters['R0_matrix5'],3),
" R023 ", round(parameters['R0_matrix8'],3),
"\n",
" R031 ", round(parameters['R0_matrix3'],3),
" R032 ", round(parameters['R0_matrix6'],3),
" R033 ", round(parameters['R0_matrix9'],3)
),
color = "Cumulative incidence",
subtitle= nicesubtitle) +
theme(legend.position="bottom") 
# part 4 - plots continued
print("Plotting the proprotion of people in each compartment over time")## [1] "Plotting the proprotion of people in each compartment over time"output %>%
select(time, cumul_incid_pct, cumul_incid1_pct, cumul_incid2_pct, cumul_incid3_pct) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dashed") +
geom_point(show.legend = FALSE) +
scale_color_manual(values = c("red","blue3","chartreuse", "bisque4"),
labels=c("ALL", "Children","Adults","Elderly")) +
scale_shape_manual(values = c(4, "o", 3, "-")) +
xlab("Time (days)") +
ylab("% of Population") +
labs(title=paste(
" gamma ", parameters['gamma'],
" R011 ", round(parameters['R0_matrix1'],3),
" R012 ", round(parameters['R0_matrix4'],3),
" R013 ", round(parameters['R0_matrix7'],3),
"\n",
" R021 ", round(parameters['R0_matrix2'],3),
" R022 ", round(parameters['R0_matrix5'],3),
" R023 ", round(parameters['R0_matrix8'],3),
"\n",
" R031 ", round(parameters['R0_matrix3'],3),
" R032 ", round(parameters['R0_matrix6'],3),
" R033 ", round(parameters['R0_matrix9'],3)
),
color = "Cumulative incidence %",
subtitle= nicesubtitle) +
theme(legend.position="bottom") 
Interventions in an age-structured population
Use the age-structured model with 3 age groups to investigate the impact of vaccinating different age groups.
We want to model a vaccine with an all-or-nothing protective effect, which is 100% effective in children and adults, but only 50% effective in the elderly. Like in earlier activities on vaccination in previous weeks, you can model vaccination before the start of the outbreak, by applying the vaccine coverage to the initial conditions.
There are 250,000 vaccine doses available for the entire population. Since the infection causes more severe disease in the elderly, your task is to model which age group to prioritise for vaccination in order to minimise infection in the elderly. Use your model to answer the following questions:
If you can only give the vaccine to one of the 3 age groups, which one would you prioritise to minimise the number of infections in the elderly? Would this also be the best strategy to reduce the overall number of infections in the population?
- If the priority is to minimize the number of infected elderly, then giving the vaccine to the elderly results in 50,023 elderly infections, reducing the age group incidence (-54.2% delta). However, the total cumulative incidence of 842,872 (84.3%) is not the lowest of the 3 options (7% delta only).
- The best option overall would be to give the vaccines to the adult age group, resulting in a total cumulative incidence of 599,453 reducing the total incidence from 90.9% down to 59.9% (-34% delta).
If you distribute the vaccine doses among the 3 age groups in proportion to their population size, which group would benefit the most in terms of the percentage reduction in the cumulative incidence achieved with vaccination? Is the reduction in the total number on infections in the elderly what you would expect given the lower vaccine efficacy in this age group?
- Distributing the vaccines proportionally results in a total cumulative incidence of 620,756 reducing the total incidence from 90.9% down to 62.1% (-31.7% delta).
- Changes in age-group incidence shows close reduction values for both children and adult groups (-31% delta, -32% delta, respectively); adults benefited the most because they consist of a larger age group and also had the highest percentage improvement.
- The elderly group shows 77,408 infections, reducing the age group incidence from 72.9% down to 51.6% (-29% delta). This amount of improvement was initially surprising given the limited (50%) vaccine efficacy for the elderly, but makes sense when we consider that the elderly see twice the amount of children + adults versus other elderly.
Baseline, no vaccination
# part 1
nicesubtitle <- "SIR Model age3c5m Age-structured baseline c3w2"
print(nicesubtitle)## [1] "SIR Model age3c5m Age-structured baseline c3w2"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000 # total population size
N1 <- round(N * .2,0) # children make up 20% of the population
N3 <- round(N * .15,0) # elderly make up 15% of the population
N2 <- N - N1 - N3 # adults make up the rest
Vdose <- 250000 # total number of vaccinations available
Vdose1 <- round(Vdose * .2,0) # children make up 20% of the population
Vdose3 <- round(Vdose * .15,0) # elderly make up 15% of the population
Vdose2 <- Vdose - Vdose1 - Vdose3 # adults make up the rest
duration <- 90 # total number of days - 3 months
tsteps <- 0.1 # chunk into 1/10th day intervals
# children-1 adult-2 elderly-3 contact parameters:
# Set up an empty contact matrix with rows for each age group and columns for each age group
contact_matrix <- matrix(0,nrow=3,ncol=3)
# Fill in the contract matrix
contact_matrix[1,1] = 7 # children with other children
contact_matrix[1,2] = 5 # children with adults
contact_matrix[1,3] = 1 # children with elderly
contact_matrix[2,1] = 2 # adults with children
contact_matrix[2,2] = 9 # adults with other adults
contact_matrix[2,3] = 1 # adults with elderly
contact_matrix[3,1] = 1 # elderly with children
contact_matrix[3,2] = 3 # elderly with adults
contact_matrix[3,3] = 2 # elderly with other elderly
b <- .05 # probability of infection per contact equals 5%, same for children and adults (note: prior exercises used the variable name p for for vaccination rate)
gamma <- 5 ** -1 # recovery rate - people remain infected on average for 5 days
# beta = infection rate day^-1 is a function of b * c)
# with n values of c, there will be n betas and thus n R0s
# initialize and fill beta matrix
beta_matrix <- matrix(0,nrow=3,ncol=3)
beta_matrix <- b * contact_matrix
# initialize and fill R0 matrix
R0_matrix <- matrix(0,nrow=3,ncol=3)
R0_matrix <- b * contact_matrix / gamma
(parameters <- c(
gamma = gamma, # recovery rate
b = b, # probability of infection per contact
contact_matrix = contact_matrix, # ave daily contacts
beta_matrix = beta_matrix, # infection rate
R0_matrix = R0_matrix # basic repro number
)) ## gamma b contact_matrix1 contact_matrix2 contact_matrix3
## 0.20 0.05 7.00 2.00 1.00
## contact_matrix4 contact_matrix5 contact_matrix6 contact_matrix7 contact_matrix8
## 5.00 9.00 3.00 1.00 1.00
## contact_matrix9 beta_matrix1 beta_matrix2 beta_matrix3 beta_matrix4
## 2.00 0.35 0.10 0.05 0.25
## beta_matrix5 beta_matrix6 beta_matrix7 beta_matrix8 beta_matrix9
## 0.45 0.15 0.05 0.05 0.10
## R0_matrix1 R0_matrix2 R0_matrix3 R0_matrix4 R0_matrix5
## 1.75 0.50 0.25 1.25 2.25
## R0_matrix6 R0_matrix7 R0_matrix8 R0_matrix9
## 0.75 0.25 0.25 0.50# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_age3m <- function(time, state, parameters) {
with(as.list(parameters), { # referenced differently
n_agegroups <- 3 # number of age groups
S <- state[1:n_agegroups] # first 3 numbers in the initial_state_values vector
I <- state[(n_agegroups+1):(2*n_agegroups)] # numbers 4 to 6 in the initial_state_values vector
R <- state[(2*n_agegroups+1):(3*n_agegroups)] # numbers 7 to 9 in the initial_state_values vector
N <- S+I+R # people in S, I and R are added separately by age group, so N is also a vector of length 3
# Force of infection acting on susceptible children
lambda <- b * contact_matrix %*% as.matrix(I/N)
# lambda <- beta_matrix %*% as.matrix(I/N)
# %*% is used to multiply matrices in R
# the lambda vector contains the forces of infection for children, adults and the elderly (length 3)
# The differential equations
dS <- -(lambda * S)
dI <- (lambda * S) - (gamma * I)
dR <- (gamma * I)
# Output
return(list(c(dS, dI, dR)))
})
}
# note different order to initial state values to work with matrix
(initial_state_values <- c(S1 = N1-1,
S2 = N2,
S3 = N3,
I1 = 1,
I2 = 0,
I3 = 0,
R1 = 0,
R2 = 0,
R3 = 0))## S1 S2 S3 I1 I2 I3 R1 R2 R3
## 199999 650000 150000 1 0 0 0 0 0# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_age3m,
parms = parameters)) %>%
mutate(
N1 = (S1+I1+R1),
N2 = (S2+I2+R2),
N3 = (S3+I3+R3),
S = (S1+S2+S3),
I = (I1+I2+I3),
R = (R1+R2+R3),
N = (S+I+R),
still_Su1 = round(S1/N1,digits=5)*100,
preval_Inf1 = round(I1/N1,digits=5)*100,
propor_Re1 = round(R1/N1,digits=5)*100,
cumul_incid1 = round(N1-S1,digits=5)
, cumul_incid1_pct = round(cumul_incid1/N1,digits=5)*100
,
still_Su2 = round(S2/N2,digits=5)*100,
preval_Inf2 = round(I2/N2,digits=5)*100,
propor_Re2 = round(R2/N2,digits=5)*100,
cumul_incid2 = round(N2-S2,digits=5)
, cumul_incid2_pct = round(cumul_incid2/N2,digits=5)*100
,
still_Su3 = round(S3/N3,digits=5)*100,
preval_Inf3 = round(I3/N3,digits=5)*100,
propor_Re3 = round(R3/N3,digits=5)*100,
cumul_incid3 = round(N3-S3,digits=5)
, cumul_incid3_pct = round(cumul_incid3/N3,digits=5)*100
,
still_Su = round(S/N,digits=5)*100,
preval_Inf = round(I/N,digits=5)*100,
propor_Re = round(R/N,digits=5)*100,
cumul_incid = round(N-S,digits=5)
, cumul_incid_pct = round(cumul_incid/N,digits=5)*100
)
print("Initial record: ")## [1] "Initial record: "output %>%
arrange(time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2 R3 N1 N2 N3 S I R
## 1 0 199999 650000 150000 1 0 0 0 0 0 2e+05 650000 150000 999999 1 0
## N still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct
## 1 1e+06 99.999 0 0 1 0
## still_Su2 preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3
## 1 100 0 0 0 0 100
## preval_Inf3 propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 0 0 0 0 100 0
## propor_Re cumul_incid cumul_incid_pct
## 1 0 1 0print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S1 S2 S3 I1 I2 I3 R1
## 901 90 9784.579 40904.55 40688.93 20.14349 80.02684 43.53274 190195.3
## R2 R3 N1 N2 N3 S I R N
## 901 609015.4 109267.5 2e+05 650000 150000 91378.05 143.7031 908478.2 1e+06
## still_Su1 preval_Inf1 propor_Re1 cumul_incid1 cumul_incid1_pct still_Su2
## 901 4.892 0.01 95.098 190215.4 95.108 6.293
## preval_Inf2 propor_Re2 cumul_incid2 cumul_incid2_pct still_Su3 preval_Inf3
## 901 0.012 93.695 609095.5 93.707 27.126 0.029
## propor_Re3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf propor_Re
## 901 72.845 109311.1 72.874 9.138 0.014 90.848
## cumul_incid cumul_incid_pct
## 901 908621.9 90.862print(paste0("Total cumulative incidence: ",
round(output$cumul_incid[nrow(output)], 0), " (",
round((output$cumul_incid[nrow(output)]/N), 3)*100, " %)"))## [1] "Total cumulative incidence: 908622 (90.9 %)"print(paste0("Cumulative incidence in children: ",
round(output$cumul_incid1[nrow(output)], 0), " (",
round((output$cumul_incid1[nrow(output)]/N1), 3)*100, " %)"))## [1] "Cumulative incidence in children: 190215 (95.1 %)"print(paste0("Proportion of total infections who are children: ",
round((output$cumul_incid1[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are children: 20.9 %"print(paste0("Cumulative incidence in adults: ",
round(output$cumul_incid2[nrow(output)], 0), " (",
round((output$cumul_incid2[nrow(output)]/N2), 3)*100, " %)"))## [1] "Cumulative incidence in adults: 609095 (93.7 %)"print(paste0("Proportion of total infections who are adults: ",
round((output$cumul_incid2[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are adults: 67 %"print(paste0("Cumulative incidence in elderly: ",
round(output$cumul_incid3[nrow(output)], 0), " (",
round((output$cumul_incid3[nrow(output)]/N3), 3)*100, " %)"))## [1] "Cumulative incidence in elderly: 109311 (72.9 %)"print(paste0("Proportion of total infections who are elderly: ",
round((output$cumul_incid3[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are elderly: 12 %"# Extract cumulative incidence for each run:
(results1a <- data.frame(child_cumul_inc = output$cumul_incid1[nrow(output)],
adult_cumul_inc = output$cumul_incid2[nrow(output)],
elderly_cumul_inc = output$cumul_incid3[nrow(output)],
total_cumul_inc = output$cumul_incid[nrow(output)]))## child_cumul_inc adult_cumul_inc elderly_cumul_inc total_cumul_inc
## 1 190215.4 609095.5 109311.1 908621.9Part 2x3, Vaccination for only for one age group: children
# part 2x
nicesubtitle <- "SIR Model age3c5m Age-structured Vaccine 1 v3 c3w2"
print(nicesubtitle)## [1] "SIR Model age3c5m Age-structured Vaccine 1 v3 c3w2"# Generalize and set individual pvacc_n to the ratio of number of available vaccination doses (250,000) over the total population (1,000,000) based on available 250,000 doses already calculated as:
# Vdose = Vdose1 + Vdose2 + Vdose3
# peff = effective vaccine coverage (in all-or-nothing = vaccination rate p * veff vaccine efficacy)
pvacc1 <- if_else((N1-1) >= Vdose, Vdose/(N1-1), 1) # vaccine coverage in the children using all available
pvacc2 <- 0 # NO vaccine coverage in adults
pvacc3 <- 0 # NO vaccine coverage in elderly
veff1 <- 1.00 # vaccine efficacy in children
veff2 <- 1.00 # vaccine efficacy in adults
veff3 <- 0.50 # vaccine efficacy in the elderly
peff1 <- pvacc1 * veff1 # effective vaccine coverage children
peff2 <- pvacc2 * veff2 # effective vaccine coverage adults
peff3 <- pvacc3 * veff3 # effective vaccine coverage elderly
# using proportional scenario differs use of pvacc and peff
(initial_state_values <- c(
S1 = if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (peff1)),0), 0),
S2 = if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (peff2)),0), 0),
S3 = if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (peff3)),0), 0),
I1 = 1,
I2 = 0,
I3 = 0,
R1 = if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (1-peff1)),0), (N1-1)),
R2 = if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (1-peff2)),0), N2),
R3 = if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (1-peff3)),0), N3)
))## S1 S2 S3 I1 I2 I3 R1 R2 R3
## 0 650000 150000 1 0 0 199999 0 0# vaccination compartment information - double book outside matrix
V1 <- if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (1-peff1)),0), (N1-1))
V2 <- if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (1-peff2)),0), N2)
V3 <- if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (1-peff3)),0), N3)
print(cbind(V1,V2,V3))## V1 V2 V3
## [1,] 199999 0 0# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_age3m,
parms = parameters)) %>%
add_column(V1 = V1, V2 = V2, V3 = V3) %>%
mutate(
N1 = (S1+I1+R1),
N2 = (S2+I2+R2),
N3 = (S3+I3+R3),
S = (S1+S2+S3),
I = (I1+I2+I3),
R = (R1+R2+R3),
V = (V1+V2+V3),
N = (S+I+R),
still_Su1 = round(S1/N1,digits=5)*100,
preval_Inf1 = round(I1/N1,digits=5)*100,
propor_Re1 = round(R1/N1,digits=5)*100,
propor_V1 = round(V1/N1,digits=5)*100,
cumul_incid1 = round(N1-S1-V1,digits=5)
, cumul_incid1_pct = round(cumul_incid1/N1,digits=5)*100
,
still_Su2 = round(S2/N2,digits=5)*100,
preval_Inf2 = round(I2/N2,digits=5)*100,
propor_Re2 = round(R2/N2,digits=5)*100,
propor_V2 = round(V2/N2,digits=5)*100,
cumul_incid2 = round(N2-S2-V2,digits=5)
, cumul_incid2_pct = round(cumul_incid2/N2,digits=5)*100
,
still_Su3 = round(S3/N3,digits=5)*100,
preval_Inf3 = round(I3/N3,digits=5)*100,
propor_Re3 = round(R3/N3,digits=5)*100,
propor_V3 = round(V3/N3,digits=5)*100,
cumul_incid3 = round(N3-S3-V3,digits=5)
, cumul_incid3_pct = round(cumul_incid3/N3,digits=5)*100
,
still_Su = round(S/N,digits=5)*100,
preval_Inf = round(I/N,digits=5)*100,
propor_Re = round(R/N,digits=5)*100,
propor_V = round(V/N,digits=5)*100,
cumul_incid = round(N-S-V,digits=5)
, cumul_incid_pct = round(cumul_incid/N,digits=5)*100
)
print("first record for the run: ")## [1] "first record for the run: "output %>%
arrange(time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2 R3 V1 V2 V3 N1 N2 N3
## 1 0 0 650000 150000 1 0 0 199999 0 0 199999 0 0 2e+05 650000 150000
## S I R V N still_Su1 preval_Inf1 propor_Re1 propor_V1
## 1 8e+05 1 199999 199999 1e+06 0 0 99.999 99.999
## cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 propor_V2
## 1 1 0 100 0 0 0
## cumul_incid2 cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 propor_V3
## 1 0 0 100 0 0 0
## cumul_incid3 cumul_incid3_pct still_Su preval_Inf propor_Re propor_V
## 1 0 0 80 0 20 20
## cumul_incid cumul_incid_pct
## 1 1 0print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S1 S2 S3 I1 I2 I3 R1 R2
## 901 90 0 77202.34 57058.96 1.522998e-08 1745.592 648.2237 2e+05 571052.1
## R3 V1 V2 V3 N1 N2 N3 S I R V
## 901 92292.82 199999 0 0 2e+05 650000 150000 134261.3 2393.816 863344.9 199999
## N still_Su1 preval_Inf1 propor_Re1 propor_V1 cumul_incid1
## 901 1e+06 0 0 100 100 1
## cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 propor_V2 cumul_incid2
## 901 0.001 11.877 0.269 87.854 0 572797.7
## cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 propor_V3 cumul_incid3
## 901 88.123 38.039 0.432 61.529 0 92941.04
## cumul_incid3_pct still_Su preval_Inf propor_Re propor_V cumul_incid
## 901 61.961 13.426 0.239 86.334 20 665739.7
## cumul_incid_pct
## 901 66.574print(paste0("Total cumulative incidence: ",
round(output$cumul_incid[nrow(output)], 0), " (",
round((output$cumul_incid[nrow(output)]/N), 3)*100, " %)"))## [1] "Total cumulative incidence: 665740 (66.6 %)"print(paste0("Cumulative incidence in children: ",
round(output$cumul_incid1[nrow(output)], 0), " (",
round((output$cumul_incid1[nrow(output)]/N1), 3)*100, " %)"))## [1] "Cumulative incidence in children: 1 (0 %)"print(paste0("Proportion of total infections who are children: ",
round((output$cumul_incid1[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are children: 0 %"print(paste0("Cumulative incidence in adults: ",
round(output$cumul_incid2[nrow(output)], 0), " (",
round((output$cumul_incid2[nrow(output)]/N2), 3)*100, " %)"))## [1] "Cumulative incidence in adults: 572798 (88.1 %)"print(paste0("Proportion of total infections who are adults: ",
round((output$cumul_incid2[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are adults: 86 %"print(paste0("Cumulative incidence in elderly: ",
round(output$cumul_incid3[nrow(output)], 0), " (",
round((output$cumul_incid3[nrow(output)]/N3), 3)*100, " %)"))## [1] "Cumulative incidence in elderly: 92941 (62 %)"print(paste0("Proportion of total infections who are elderly: ",
round((output$cumul_incid3[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are elderly: 14 %"# Extract cumulative incidence for each run:
(results2x <- data.frame(child_cumul_inc = output$cumul_incid1[nrow(output)],
adult_cumul_inc = output$cumul_incid2[nrow(output)],
elderly_cumul_inc = output$cumul_incid3[nrow(output)],
total_cumul_inc = output$cumul_incid[nrow(output)]))## child_cumul_inc adult_cumul_inc elderly_cumul_inc total_cumul_inc
## 1 1 572797.7 92941.04 665739.7Part 3x3, Vaccination for only for one age group: adults
The code below does not generate similar results to the class-provided solution. Discrepancy found only with adult vaccine numbers caused by a rounding issue where the solution key used an imprecise rounding of .38 instead of something more realistic such as 0.3846154 which would generate an initial adult vaccination number of 250,000 (all doses) instead of the faulty 247,000.
# part 3x
nicesubtitle <- "SIR Model age3c5m Age-structured Vaccine 2 v3 c3w2"
print(nicesubtitle)## [1] "SIR Model age3c5m Age-structured Vaccine 2 v3 c3w2"# Generalize and set individual pvacc_n to the ratio of number of available vaccination doses (250,000) over the total population (1,000,000) based on available 250,000 doses already calculated as:
# Vdose = Vdose1 + Vdose2 + Vdose3
# peff = effective vaccine coverage (in all-or-nothing = vaccination rate p * veff vaccine efficacy)
pvacc1 <- 0 # NO vaccine coverage in children
pvacc2 <- if_else((N2) >= Vdose, Vdose/(N2), 1) # vaccine coverage in the adults using all available
pvacc3 <- 0 # NO vaccine coverage in elderly
veff1 <- 1.00 # vaccine efficacy in children
veff2 <- 1.00 # vaccine efficacy in adults
veff3 <- 0.50 # vaccine efficacy in the elderly
peff1 <- pvacc1 * veff1 # effective vaccine coverage children
peff2 <- pvacc2 * veff2 # effective vaccine coverage adults
peff3 <- pvacc3 * veff3 # effective vaccine coverage elderly
# using proportional scenario differs use of pvacc and peff
(initial_state_values <- c(
S1 = if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (peff1)),0), 0),
S2 = if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (peff2)),0), 0),
S3 = if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (peff3)),0), 0),
I1 = 1,
I2 = 0,
I3 = 0,
R1 = if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (1-peff1)),0), (N1-1)),
R2 = if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (1-peff2)),0), N2),
R3 = if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (1-peff3)),0), N3)
))## S1 S2 S3 I1 I2 I3 R1 R2 R3
## 199999 400000 150000 1 0 0 0 250000 0# vaccination compartment information - double book outside matrix
V1 <- if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (1-peff1)),0), (N1-1))
V2 <- if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (1-peff2)),0), N2)
V3 <- if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (1-peff3)),0), N3)
print(cbind(V1,V2,V3))## V1 V2 V3
## [1,] 0 250000 0# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_age3m,
parms = parameters)) %>%
add_column(V1 = V1, V2 = V2, V3 = V3) %>%
mutate(
N1 = (S1+I1+R1),
N2 = (S2+I2+R2),
N3 = (S3+I3+R3),
S = (S1+S2+S3),
I = (I1+I2+I3),
R = (R1+R2+R3),
V = (V1+V2+V3),
N = (S+I+R),
still_Su1 = round(S1/N1,digits=5)*100,
preval_Inf1 = round(I1/N1,digits=5)*100,
propor_Re1 = round(R1/N1,digits=5)*100,
propor_V1 = round(V1/N1,digits=5)*100,
cumul_incid1 = round(N1-S1-V1,digits=5)
, cumul_incid1_pct = round(cumul_incid1/N1,digits=5)*100
,
still_Su2 = round(S2/N2,digits=5)*100,
preval_Inf2 = round(I2/N2,digits=5)*100,
propor_Re2 = round(R2/N2,digits=5)*100,
propor_V2 = round(V2/N2,digits=5)*100,
cumul_incid2 = round(N2-S2-V2,digits=5)
, cumul_incid2_pct = round(cumul_incid2/N2,digits=5)*100
,
still_Su3 = round(S3/N3,digits=5)*100,
preval_Inf3 = round(I3/N3,digits=5)*100,
propor_Re3 = round(R3/N3,digits=5)*100,
propor_V3 = round(V3/N3,digits=5)*100,
cumul_incid3 = round(N3-S3-V3,digits=5)
, cumul_incid3_pct = round(cumul_incid3/N3,digits=5)*100
,
still_Su = round(S/N,digits=5)*100,
preval_Inf = round(I/N,digits=5)*100,
propor_Re = round(R/N,digits=5)*100,
propor_V = round(V/N,digits=5)*100,
cumul_incid = round(N-S-V,digits=5)
, cumul_incid_pct = round(cumul_incid/N,digits=5)*100
)
print("first record for the run: ")## [1] "first record for the run: "output %>%
arrange(time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2 R3 V1 V2 V3 N1 N2
## 1 0 199999 4e+05 150000 1 0 0 0 250000 0 0 250000 0 2e+05 650000
## N3 S I R V N still_Su1 preval_Inf1 propor_Re1
## 1 150000 749999 1 250000 250000 1e+06 99.999 0 0
## propor_V1 cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2
## 1 0 1 0 61.538 0 38.462
## propor_V2 cumul_incid2 cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3
## 1 38.462 0 0 100 0 0
## propor_V3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf propor_Re
## 1 0 0 0 75 0 25
## propor_V cumul_incid cumul_incid_pct
## 1 25 1 0print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S1 S2 S3 I1 I2 I3 R1
## 901 90 18900.34 70585.88 61061.16 414.5482 1304.462 644.3864 180685.1
## R2 R3 V1 V2 V3 N1 N2 N3 S I
## 901 578109.7 88294.45 0 250000 0 2e+05 650000 150000 150547.4 2363.397
## R V N still_Su1 preval_Inf1 propor_Re1 propor_V1
## 901 847089.2 250000 1e+06 9.45 0.207 90.343 0
## cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 propor_V2
## 901 181099.7 90.55 10.859 0.201 88.94 38.462
## cumul_incid2 cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 propor_V3
## 901 329414.1 50.679 40.707 0.43 58.863 0
## cumul_incid3 cumul_incid3_pct still_Su preval_Inf propor_Re propor_V
## 901 88938.84 59.293 15.055 0.236 84.709 25
## cumul_incid cumul_incid_pct
## 901 599452.6 59.945print(paste0("Total cumulative incidence: ",
round(output$cumul_incid[nrow(output)], 0), " (",
round((output$cumul_incid[nrow(output)]/N), 3)*100, " %)"))## [1] "Total cumulative incidence: 599453 (59.9 %)"print(paste0("Cumulative incidence in children: ",
round(output$cumul_incid1[nrow(output)], 0), " (",
round((output$cumul_incid1[nrow(output)]/N1), 3)*100, " %)"))## [1] "Cumulative incidence in children: 181100 (90.5 %)"print(paste0("Proportion of total infections who are children: ",
round((output$cumul_incid1[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are children: 30.2 %"print(paste0("Cumulative incidence in adults: ",
round(output$cumul_incid2[nrow(output)], 0), " (",
round((output$cumul_incid2[nrow(output)]/N2), 3)*100, " %)"))## [1] "Cumulative incidence in adults: 329414 (50.7 %)"print(paste0("Proportion of total infections who are adults: ",
round((output$cumul_incid2[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are adults: 55 %"print(paste0("Cumulative incidence in elderly: ",
round(output$cumul_incid3[nrow(output)], 0), " (",
round((output$cumul_incid3[nrow(output)]/N3), 3)*100, " %)"))## [1] "Cumulative incidence in elderly: 88939 (59.3 %)"print(paste0("Proportion of total infections who are elderly: ",
round((output$cumul_incid3[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are elderly: 14.8 %"# Extract cumulative incidence for each run:
(results3x <- data.frame(child_cumul_inc = output$cumul_incid1[nrow(output)],
adult_cumul_inc = output$cumul_incid2[nrow(output)],
elderly_cumul_inc = output$cumul_incid3[nrow(output)],
total_cumul_inc = output$cumul_incid[nrow(output)]))## child_cumul_inc adult_cumul_inc elderly_cumul_inc total_cumul_inc
## 1 181099.7 329414.1 88938.84 599452.6Part 4x3, Vaccination for only for one age group: elderly 50% efficacy
# part 4x
nicesubtitle <- "SIR Model age3c5m Age-structured Vaccine 3 50% eff v3 c3w2"
print(nicesubtitle)## [1] "SIR Model age3c5m Age-structured Vaccine 3 50% eff v3 c3w2"# Generalize and set individual pvacc_n to the ratio of number of available vaccination doses (250,000) over the total population (1,000,000) based on available 250,000 doses already calculated as:
# Vdose = Vdose1 + Vdose2 + Vdose3
# peff = effective vaccine coverage (in all-or-nothing = vaccination rate p * veff vaccine efficacy)
pvacc1 <- 0 # NO vaccine coverage in children
pvacc2 <- 0 # NO vaccine coverage in adults
pvacc3 <- if_else((N3) >= Vdose, Vdose/(N3), 1) # vaccine coverage in the elderly using all available
veff1 <- 1.00 # vaccine efficacy in children
veff2 <- 1.00 # vaccine efficacy in adults
veff3 <- 0.50 # vaccine efficacy in the elderly
peff1 <- pvacc1 * veff1 # effective vaccine coverage children
peff2 <- pvacc2 * veff2 # effective vaccine coverage adults
peff3 <- pvacc3 * veff3 # effective vaccine coverage elderly
# using proportional scenario differs use of pvacc and peff
(initial_state_values <- c(
S1 = if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (peff1)),0), 0),
S2 = if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (peff2)),0), 0),
S3 = if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (peff3)),0), 0),
I1 = 1,
I2 = 0,
I3 = 0,
R1 = if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (1-peff1)),0), (N1-1)),
R2 = if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (1-peff2)),0), N2),
R3 = if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (1-peff3)),0), N3)
))## S1 S2 S3 I1 I2 I3 R1 R2 R3
## 199999 650000 75000 1 0 0 0 0 75000# vaccination compartment information - double book outside matrix
V1 <- if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (1-peff1)),0), (N1-1))
V2 <- if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (1-peff2)),0), N2)
V3 <- if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (1-peff3)),0), N3)
print(cbind(V1,V2,V3))## V1 V2 V3
## [1,] 0 0 75000# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_age3m,
parms = parameters)) %>%
add_column(V1 = V1, V2 = V2, V3 = V3) %>%
mutate(
N1 = (S1+I1+R1),
N2 = (S2+I2+R2),
N3 = (S3+I3+R3),
S = (S1+S2+S3),
I = (I1+I2+I3),
R = (R1+R2+R3),
V = (V1+V2+V3),
N = (S+I+R),
still_Su1 = round(S1/N1,digits=5)*100,
preval_Inf1 = round(I1/N1,digits=5)*100,
propor_Re1 = round(R1/N1,digits=5)*100,
propor_V1 = round(V1/N1,digits=5)*100,
cumul_incid1 = round(N1-S1-V1,digits=5)
, cumul_incid1_pct = round(cumul_incid1/N1,digits=5)*100
,
still_Su2 = round(S2/N2,digits=5)*100,
preval_Inf2 = round(I2/N2,digits=5)*100,
propor_Re2 = round(R2/N2,digits=5)*100,
propor_V2 = round(V2/N2,digits=5)*100,
cumul_incid2 = round(N2-S2-V2,digits=5)
, cumul_incid2_pct = round(cumul_incid2/N2,digits=5)*100
,
still_Su3 = round(S3/N3,digits=5)*100,
preval_Inf3 = round(I3/N3,digits=5)*100,
propor_Re3 = round(R3/N3,digits=5)*100,
propor_V3 = round(V3/N3,digits=5)*100,
cumul_incid3 = round(N3-S3-V3,digits=5)
, cumul_incid3_pct = round(cumul_incid3/N3,digits=5)*100
,
still_Su = round(S/N,digits=5)*100,
preval_Inf = round(I/N,digits=5)*100,
propor_Re = round(R/N,digits=5)*100,
propor_V = round(V/N,digits=5)*100,
cumul_incid = round(N-S-V,digits=5)
, cumul_incid_pct = round(cumul_incid/N,digits=5)*100
)
print("first record for the run: ")## [1] "first record for the run: "output %>%
arrange(time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2 R3 V1 V2 V3 N1 N2 N3
## 1 0 199999 650000 75000 1 0 0 0 0 75000 0 0 75000 2e+05 650000 150000
## S I R V N still_Su1 preval_Inf1 propor_Re1 propor_V1
## 1 924999 1 75000 75000 1e+06 99.999 0 0 0
## cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 propor_V2
## 1 1 0 100 0 0 0
## cumul_incid2 cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 propor_V3
## 1 0 0 50 0 50 50
## cumul_incid3 cumul_incid3_pct still_Su preval_Inf propor_Re propor_V
## 1 0 0 92.5 0 7.5 7.5
## cumul_incid cumul_incid_pct
## 1 1 0print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S1 S2 S3 I1 I2 I3 R1
## 901 90 11029.67 46120.86 24977.38 24.77333 98.34464 23.02433 188945.6
## R2 R3 V1 V2 V3 N1 N2 N3 S I
## 901 603780.8 124999.6 0 0 75000 2e+05 650000 150000 82127.92 146.1423
## R V N still_Su1 preval_Inf1 propor_Re1 propor_V1
## 901 917725.9 75000 1e+06 5.515 0.012 94.473 0
## cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 propor_V2
## 901 188970.3 94.485 7.096 0.015 92.889 0
## cumul_incid2 cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 propor_V3
## 901 603879.1 92.904 16.652 0.015 83.333 50
## cumul_incid3 cumul_incid3_pct still_Su preval_Inf propor_Re propor_V
## 901 50022.62 33.348 8.213 0.015 91.773 7.5
## cumul_incid cumul_incid_pct
## 901 842872.1 84.287print(paste0("Total cumulative incidence: ",
round(output$cumul_incid[nrow(output)], 0), " (",
round((output$cumul_incid[nrow(output)]/N), 3)*100, " %)"))## [1] "Total cumulative incidence: 842872 (84.3 %)"print(paste0("Cumulative incidence in children: ",
round(output$cumul_incid1[nrow(output)], 0), " (",
round((output$cumul_incid1[nrow(output)]/N1), 3)*100, " %)"))## [1] "Cumulative incidence in children: 188970 (94.5 %)"print(paste0("Proportion of total infections who are children: ",
round((output$cumul_incid1[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are children: 22.4 %"print(paste0("Cumulative incidence in adults: ",
round(output$cumul_incid2[nrow(output)], 0), " (",
round((output$cumul_incid2[nrow(output)]/N2), 3)*100, " %)"))## [1] "Cumulative incidence in adults: 603879 (92.9 %)"print(paste0("Proportion of total infections who are adults: ",
round((output$cumul_incid2[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are adults: 71.6 %"print(paste0("Cumulative incidence in elderly: ",
round(output$cumul_incid3[nrow(output)], 0), " (",
round((output$cumul_incid3[nrow(output)]/N3), 3)*100, " %)"))## [1] "Cumulative incidence in elderly: 50023 (33.3 %)"print(paste0("Proportion of total infections who are elderly: ",
round((output$cumul_incid3[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are elderly: 5.9 %"# Extract cumulative incidence for each run:
(results4x <- data.frame(child_cumul_inc = output$cumul_incid1[nrow(output)],
adult_cumul_inc = output$cumul_incid2[nrow(output)],
elderly_cumul_inc = output$cumul_incid3[nrow(output)],
total_cumul_inc = output$cumul_incid[nrow(output)]))## child_cumul_inc adult_cumul_inc elderly_cumul_inc total_cumul_inc
## 1 188970.3 603879.1 50022.62 842872.1Part 5x3, 25% Vaccination in population proportion with elderly 50% efficacy
# part 5x
nicesubtitle <- "SIR Model age3c5m Age-structured Vaccine 4 proportional v3 c3w2"
print(nicesubtitle)## [1] "SIR Model age3c5m Age-structured Vaccine 4 proportional v3 c3w2"# Generalize and set individual pvacc_n to the ratio of number of available vaccination doses (250,000) over the total population (1,000,000) based on available 250,000 doses already calculated as:
# Vdose = Vdose1 + Vdose2 + Vdose3
# peff = effective vaccine coverage (in all-or-nothing = vaccination rate p * veff vaccine efficacy)
pvacc1 <- Vdose1 / N1 # vaccine coverage in children
pvacc2 <- Vdose2 / N2 # vaccine coverage in adults
pvacc3 <- Vdose3 / N3 # vaccine coverage in the elderly
veff1 <- 1.00 # vaccine efficacy in children
veff2 <- 1.00 # vaccine efficacy in adults
veff3 <- 0.50 # vaccine efficacy in the elderly
peff1 <- pvacc1 * veff1 # effective vaccine coverage children
peff2 <- pvacc2 * veff2 # effective vaccine coverage adults
peff3 <- pvacc3 * veff3 # effective vaccine coverage elderly
# using proportional scenario differs use of pvacc and peff
(initial_state_values <- c(
S1 = if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (peff1)),0), 0),
S2 = if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (peff2)),0), 0),
S3 = if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (peff3)),0), 0),
I1 = 1,
I2 = 0,
I3 = 0,
R1 = if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (1-peff1)),0), (N1-1)),
R2 = if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (1-peff2)),0), N2),
R3 = if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (1-peff3)),0), N3)
))## S1 S2 S3 I1 I2 I3 R1 R2 R3
## 149999 487500 131250 1 0 0 50000 162500 18750# vaccination compartment information - double book outside matrix
V1 <- if_else((N1-1) >= pvacc1 * N1, round((N1-1) - ((N1-1) * (1-peff1)),0), (N1-1))
V2 <- if_else((N2) >= pvacc2 * N2, round((N2) - ((N2) * (1-peff2)),0), N2)
V3 <- if_else((N3) >= pvacc3 * N3, round((N3) - ((N3) * (1-peff3)),0), N3)
print(cbind(V1,V2,V3))## V1 V2 V3
## [1,] 50000 162500 18750# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_age3m,
parms = parameters)) %>%
add_column(V1 = V1, V2 = V2, V3 = V3) %>%
mutate(
N1 = (S1+I1+R1),
N2 = (S2+I2+R2),
N3 = (S3+I3+R3),
S = (S1+S2+S3),
I = (I1+I2+I3),
R = (R1+R2+R3),
V = (V1+V2+V3),
N = (S+I+R),
still_Su1 = round(S1/N1,digits=5)*100,
preval_Inf1 = round(I1/N1,digits=5)*100,
propor_Re1 = round(R1/N1,digits=5)*100,
propor_V1 = round(V1/N1,digits=5)*100,
cumul_incid1 = round(N1-S1-V1,digits=5)
, cumul_incid1_pct = round(cumul_incid1/N1,digits=5)*100
,
still_Su2 = round(S2/N2,digits=5)*100,
preval_Inf2 = round(I2/N2,digits=5)*100,
propor_Re2 = round(R2/N2,digits=5)*100,
propor_V2 = round(V2/N2,digits=5)*100,
cumul_incid2 = round(N2-S2-V2,digits=5)
, cumul_incid2_pct = round(cumul_incid2/N2,digits=5)*100
,
still_Su3 = round(S3/N3,digits=5)*100,
preval_Inf3 = round(I3/N3,digits=5)*100,
propor_Re3 = round(R3/N3,digits=5)*100,
propor_V3 = round(V3/N3,digits=5)*100,
cumul_incid3 = round(N3-S3-V3,digits=5)
, cumul_incid3_pct = round(cumul_incid3/N3,digits=5)*100
,
still_Su = round(S/N,digits=5)*100,
preval_Inf = round(I/N,digits=5)*100,
propor_Re = round(R/N,digits=5)*100,
propor_V = round(V/N,digits=5)*100,
cumul_incid = round(N-S-V,digits=5)
, cumul_incid_pct = round(cumul_incid/N,digits=5)*100
)
print("first record for the run: ")## [1] "first record for the run: "output %>%
arrange(time) %>%
head(1)## time S1 S2 S3 I1 I2 I3 R1 R2 R3 V1 V2 V3
## 1 0 149999 487500 131250 1 0 0 50000 162500 18750 50000 162500 18750
## N1 N2 N3 S I R V N still_Su1 preval_Inf1
## 1 2e+05 650000 150000 768749 1 231250 231250 1e+06 74.999 0
## propor_Re1 propor_V1 cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2
## 1 25 25 1 0 75 0
## propor_Re2 propor_V2 cumul_incid2 cumul_incid2_pct still_Su3 preval_Inf3
## 1 25 25 0 0 87.5 0
## propor_Re3 propor_V3 cumul_incid3 cumul_incid3_pct still_Su preval_Inf
## 1 12.5 12.5 0 0 76.875 0
## propor_Re propor_V cumul_incid cumul_incid_pct
## 1 23.125 23.125 1 0print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S1 S2 S3 I1 I2 I3 R1
## 901 90 19214.94 74936.29 53842.36 657.5161 2408.102 869.3974 180127.5
## R2 R3 V1 V2 V3 N1 N2 N3 S I
## 901 572655.6 95288.24 50000 162500 18750 2e+05 650000 150000 147993.6 3935.015
## R V N still_Su1 preval_Inf1 propor_Re1 propor_V1
## 901 848071.4 231250 1e+06 9.607 0.329 90.064 25
## cumul_incid1 cumul_incid1_pct still_Su2 preval_Inf2 propor_Re2 propor_V2
## 901 130785.1 65.393 11.529 0.37 88.101 25
## cumul_incid2 cumul_incid2_pct still_Su3 preval_Inf3 propor_Re3 propor_V3
## 901 412563.7 63.471 35.895 0.58 63.525 12.5
## cumul_incid3 cumul_incid3_pct still_Su preval_Inf propor_Re propor_V
## 901 77407.64 51.605 14.799 0.394 84.807 23.125
## cumul_incid cumul_incid_pct
## 901 620756.4 62.076print(paste0("Total cumulative incidence: ",
round(output$cumul_incid[nrow(output)], 0), " (",
round((output$cumul_incid[nrow(output)]/N), 3)*100, " %)"))## [1] "Total cumulative incidence: 620756 (62.1 %)"print(paste0("Cumulative incidence in children: ",
round(output$cumul_incid1[nrow(output)], 0), " (",
round((output$cumul_incid1[nrow(output)]/N1), 3)*100, " %)"))## [1] "Cumulative incidence in children: 130785 (65.4 %)"print(paste0("Proportion of total infections who are children: ",
round((output$cumul_incid1[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are children: 21.1 %"print(paste0("Cumulative incidence in adults: ",
round(output$cumul_incid2[nrow(output)], 0), " (",
round((output$cumul_incid2[nrow(output)]/N2), 3)*100, " %)"))## [1] "Cumulative incidence in adults: 412564 (63.5 %)"print(paste0("Proportion of total infections who are adults: ",
round((output$cumul_incid2[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are adults: 66.5 %"print(paste0("Cumulative incidence in elderly: ",
round(output$cumul_incid3[nrow(output)], 0), " (",
round((output$cumul_incid3[nrow(output)]/N3), 3)*100, " %)"))## [1] "Cumulative incidence in elderly: 77408 (51.6 %)"print(paste0("Proportion of total infections who are elderly: ",
round((output$cumul_incid3[nrow(output)]/output$cumul_incid[nrow(output)]), 3)*100, " %"))## [1] "Proportion of total infections who are elderly: 12.5 %"# Extract cumulative incidence for each run:
(results5x <- data.frame(child_cumul_inc = output$cumul_incid1[nrow(output)],
adult_cumul_inc = output$cumul_incid2[nrow(output)],
elderly_cumul_inc = output$cumul_incid3[nrow(output)],
total_cumul_inc = output$cumul_incid[nrow(output)]))## child_cumul_inc adult_cumul_inc elderly_cumul_inc total_cumul_inc
## 1 130785.1 412563.7 77407.64 620756.4resultsx <- rbind(results1a, results2x, results3x, results4x, results5x)
resultsx <- resultsx %>%
mutate(
child_cumul_chg = round((child_cumul_inc - resultsx$child_cumul_inc[1])/resultsx$child_cumul_inc[1],3) * 100,
adult_cumul_chg = round((adult_cumul_inc - resultsx$adult_cumul_inc[1])/resultsx$adult_cumul_inc[1],3) * 100,
elderly_cumul_chg = round((elderly_cumul_inc - resultsx$elderly_cumul_inc[1])/resultsx$elderly_cumul_inc[1],3) * 100,
total_cumul_chg = round((total_cumul_inc - resultsx$total_cumul_inc[1])/resultsx$total_cumul_inc[1],3) * 100)
resultsx## child_cumul_inc adult_cumul_inc elderly_cumul_inc total_cumul_inc
## 1 190215.4 609095.5 109311.07 908621.9
## 2 1.0 572797.7 92941.04 665739.7
## 3 181099.7 329414.1 88938.84 599452.6
## 4 188970.3 603879.1 50022.62 842872.1
## 5 130785.1 412563.7 77407.64 620756.4
## child_cumul_chg adult_cumul_chg elderly_cumul_chg total_cumul_chg
## 1 0.0 0.0 0.0 0.0
## 2 -100.0 -6.0 -15.0 -26.7
## 3 -4.8 -45.9 -18.6 -34.0
## 4 -0.7 -0.9 -54.2 -7.2
## 5 -31.2 -32.3 -29.2 -31.7Class-provided solution – interventions in an age-structured population
** Note: Use the revised assignment code above instead of solution provided below which inaccurately caclculates the adult-only scenario.**
# INPUT
# Set up an empty contact matrix with rows for each age group and columns for each age group
contact_matrix <- matrix(0,nrow=3,ncol=3)
# Fill in the contract matrix
contact_matrix[1,1] = 7 # daily number of contacts that children make with each other
contact_matrix[1,2] = 5 # daily number of contacts that children make with adults
contact_matrix[1,3] = 1 # daily number of contacts that children make with the elderly
contact_matrix[2,1] = 2 # daily number of contacts that adults make with children
contact_matrix[2,2] = 9 # daily number of contacts that adults make with each other
contact_matrix[2,3] = 1 # daily number of contacts that adults make with the elderly
contact_matrix[3,1] = 1 # daily number of contacts that elderly people make with children
contact_matrix[3,2] = 3 # daily number of contacts that elderly people make with adults
contact_matrix[3,3] = 2 # daily number of contacts that elderly people make with each other
# The contact_matrix now looks exactly like the one in the activity instructions. We add this matrix as a parameter below.
# Parameters
parameters <- c(b = 0.05, # the probability of infection per contact is 5%
contact_matrix = contact_matrix, # the age-specific average number of daily contacts (defined above)
gamma = 1/5) # the rate of recovery is 1/5 per day
# Run simulation for 3 months
times <- seq(from = 0, to = 90, by = 0.1)
# MODEL FUNCTION
sir_age_model <- function(time, state, parameters) {
with(as.list(parameters), {
n_agegroups <- 3 # number of age groups
S <- state[1:n_agegroups] # assign to S the first 3 numbers in the initial_state_values vector
I <- state[(n_agegroups+1):(2*n_agegroups)] # assign to I numbers 4 to 6 in the initial_state_values vector
R <- state[(2*n_agegroups+1):(3*n_agegroups)] # assign to R numbers 7 to 9 in the initial_state_values vector
N <- S+I+R # people in S, I and R are added separately by age group, so N is also a vector of length 3
# Defining the force of infection
# Force of infection acting on susceptible children
lambda <- b * contact_matrix %*% as.matrix(I/N)
# %*% is used to multiply matrices in R
# the lambda vector contains the forces of infection for children, adults and the elderly (length 3)
# The differential equations
# Rate of change in children:
dS <- -lambda * S
dI <- lambda * S - gamma * I
dR <- gamma * I
# Output
return(list(c(dS, dI, dR)))
})
}Modelling vaccination scenarios
If you can only give the vaccine to one of the 3 age groups, which one would you prioritise to minimise the number of infections in the elderly? Would this also be the best strategy to reduce the overall number of infections in the population?
# 1) Vaccinating only children
# There are more doses of vaccine than children in the population, so the vaccine coverage will be 100% in children
# and 0% in the other groups as per the question
vacc_cov1 <- 1 # vaccine coverage in children
vacc_cov2 <- 0 # vaccine coverage in adults
vacc_cov3 <- 0 # vaccine coverage in the elderly
vacc_eff3 <- 0.5 # vaccine efficacy in the elderly (100% in the other age groups)
# Effective vaccine coverage for each age group:
p1 <- vacc_cov1
p2 <- vacc_cov2
p3 <- vacc_cov3 * vacc_eff3
# Population size in total and for each age group:
N <- 1000000
N1 <- 0.2*N
N2 <- 0.65*N
N3 <- 0.15*N
# Fill in initial state values for a naive population based on effective vaccine coverage:
initial_state_values <- c(S1 = N1-p1*N1,
S2 = N2-p2*N2,
S3 = N3-p3*N3,
I1 = 1, # the outbreak starts with 1 infected person (can be of either age)
I2 = 0,
I3 = 0,
R1 = p1*N1,
R2 = p2*N2,
R3 = p3*N3)
# Run model output
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_age_model,
parms = parameters))
# Calculate cumulative incidence in each age group:
results1 <- data.frame(child_cum_inc = output$S1[1]-output$S1[nrow(output)],
adult_cum_inc = output$S2[1]-output$S2[nrow(output)],
elderly_cum_inc = output$S3[1]-output$S3[nrow(output)],
total_cum_inc = sum(output[1,c("S1", "S2", "S3")])-sum(output[nrow(output),c("S1", "S2", "S3")]))
print(results1)## child_cum_inc adult_cum_inc elderly_cum_inc total_cum_inc
## 1 0 572797.7 92941.04 665738.7Giving all available vaccine doses to children would prevent all infections in children, but still result in 92,941 infections in the elderly, and a total number of infections of 665,739 by the end of the epidemic.
Let’s compare this with the output when giving all the vaccine doses to adults or the elderly instead:
# 2) Vaccinating only adults
# Vaccine coverage in adults = 250,000/650,000
# 0% in the other groups as per the question
vacc_cov1 <- 0 # vaccine coverage in children
vacc_cov2 <- 0.38 # vaccine coverage in adults
vacc_cov3 <- 0 # vaccine coverage in the elderly
vacc_eff3 <- 0.5 # vaccine efficacy in the elderly (100% in the other age groups)
# Effective vaccine coverage for each age group:
p1 <- vacc_cov1
p2 <- vacc_cov2
p3 <- vacc_cov3 * vacc_eff3
# Population size in total and for each age group:
N <- 1000000
N1 <- 0.2*N
N2 <- 0.65*N
N3 <- 0.15*N
# Fill in initial state values for a naive population based on effective vaccine coverage:
initial_state_values <- c(S1 = N1-p1*N1,
S2 = N2-p2*N2,
S3 = N3-p3*N3,
I1 = 1, # the outbreak starts with 1 infected person (can be of either age)
I2 = 0,
I3 = 0,
R1 = p1*N1,
R2 = p2*N2,
R3 = p3*N3)
# Run model output
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_age_model,
parms = parameters))
# Calculate cumulative incidence in each age group:
results2 <- data.frame(child_cum_inc = output$S1[1]-output$S1[nrow(output)],
adult_cum_inc = output$S2[1]-output$S2[nrow(output)],
elderly_cum_inc = output$S3[1]-output$S3[nrow(output)],
total_cum_inc = sum(output[1,c("S1", "S2", "S3")])-sum(output[nrow(output),c("S1", "S2", "S3")]))
# 3) Vaccinating only elderly people
# There are more doses of vaccine than elderly people in the population, so the vaccine coverage will be 100%
# 0% in the other groups as per the question
vacc_cov1 <- 0 # vaccine coverage in children
vacc_cov2 <- 0 # vaccine coverage in adults
vacc_cov3 <- 1 # vaccine coverage in the elderly
vacc_eff3 <- 0.5 # vaccine efficacy in the elderly (100% in the other age groups)
# Effective vaccine coverage for each age group:
p1 <- vacc_cov1
p2 <- vacc_cov2
p3 <- vacc_cov3 * vacc_eff3
# Population size in total and for each age group:
N <- 1000000
N1 <- 0.2*N
N2 <- 0.65*N
N3 <- 0.15*N
# Fill in initial state values for a naive population based on effective vaccine coverage:
initial_state_values <- c(S1 = N1-p1*N1,
S2 = N2-p2*N2,
S3 = N3-p3*N3,
I1 = 1, # the outbreak starts with 1 infected person (can be of either age)
I2 = 0,
I3 = 0,
R1 = p1*N1,
R2 = p2*N2,
R3 = p3*N3)
# Run model output
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_age_model,
parms = parameters))
# Calculate cumulative incidence in each age group:
results3 <- data.frame(child_cum_inc = output$S1[1]-output$S1[nrow(output)],
adult_cum_inc = output$S2[1]-output$S2[nrow(output)],
elderly_cum_inc = output$S3[1]-output$S3[nrow(output)],
total_cum_inc = sum(output[1,c("S1", "S2", "S3")])-sum(output[nrow(output),c("S1", "S2", "S3")]))
print(rbind(results1, results2, results3))## child_cum_inc adult_cum_inc elderly_cum_inc total_cum_inc
## 1 0.0 572797.7 92941.04 665738.7
## 2 181263.6 332793.4 89262.62 603319.6
## 3 188970.3 603879.1 50022.62 842872.0The cumulative incidence in the elderly is lowest if we only vaccinate everyone in the elderly age group (only 50,023 infections), despite the low vaccine efficacy in this age group! However, with this strategy we also get a substantially larger total number of infections than if vaccinating only children or only adults (842,872 infections vs. 665,739 and 603,320 respectively).
The worst strategy for the given question would be to only vaccinate children, since this neither minimises the number of infections in the elderly nor in total.
If you distribute the vaccine doses among the 3 age groups in proportion to their population size, which group would benefit the most in terms of the percentage reduction in the cumulative incidence achieved with vaccination? Is the reduction in the total number on infections in the elderly what you would expect given the lower vaccine efficacy in this age group?
First, we need to calculate the baseline prevalence without vaccination, then the model the vaccine scenario:
# Baseline prevalence (no vaccine)
vacc_cov1 <- 0 # vaccine coverage in children
vacc_cov2 <- 0 # vaccine coverage in adults
vacc_cov3 <- 0 # vaccine coverage in the elderly
vacc_eff3 <- 0.5 # vaccine efficacy in the elderly (100% in the other age groups)
# Effective vaccine coverage for each age group:
p1 <- vacc_cov1
p2 <- vacc_cov2
p3 <- vacc_cov3 * vacc_eff3
# Population size in total and for each age group:
N <- 1000000
N1 <- 0.2*N
N2 <- 0.65*N
N3 <- 0.15*N
# Fill in initial state values for a naive population based on effective vaccine coverage:
initial_state_values <- c(S1 = N1-p1*N1,
S2 = N2-p2*N2,
S3 = N3-p3*N3,
I1 = 1, # the outbreak starts with 1 infected person (can be of either age)
I2 = 0,
I3 = 0,
R1 = p1*N1,
R2 = p2*N2,
R3 = p3*N3)
# Run model output
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_age_model,
parms = parameters))
# Calculate cumulative incidence in each age group:
baseline_results <- data.frame(child_cum_inc = output$S1[1]-output$S1[nrow(output)],
adult_cum_inc = output$S2[1]-output$S2[nrow(output)],
elderly_cum_inc = output$S3[1]-output$S3[nrow(output)],
total_cum_inc = sum(output[1,c("S1", "S2", "S3")])-
sum(output[nrow(output),c("S1", "S2", "S3")]))
# Distributing vaccine doses among age groups proportionally to population size
# 250,000 doses/1 million = 25% coverage in each age group
vacc_cov1 <- 0.25 # vaccine coverage in children
vacc_cov2 <- 0.25 # vaccine coverage in adults
vacc_cov3 <- 0.25 # vaccine coverage in the elderly
vacc_eff3 <- 0.5 # vaccine efficacy in the elderly (100% in the other age groups)
# Effective vaccine coverage for each age group:
p1 <- vacc_cov1
p2 <- vacc_cov2
p3 <- vacc_cov3 * vacc_eff3
# Population size in total and for each age group:
N <- 1000000
N1 <- 0.2*N
N2 <- 0.65*N
N3 <- 0.15*N
# Fill in initial state values for a naive population based on effective vaccine coverage:
initial_state_values <- c(S1 = N1-p1*N1,
S2 = N2-p2*N2,
S3 = N3-p3*N3,
I1 = 1, # the outbreak starts with 1 infected person (can be of either age)
I2 = 0,
I3 = 0,
R1 = p1*N1,
R2 = p2*N2,
R3 = p3*N3)
# Run model output
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_age_model,
parms = parameters))
# Calculate cumulative incidence in each age group:
results4 <- data.frame(child_cum_inc = output$S1[1]-output$S1[nrow(output)],
adult_cum_inc = output$S2[1]-output$S2[nrow(output)],
elderly_cum_inc = output$S3[1]-output$S3[nrow(output)],
total_cum_inc = sum(output[1,c("S1", "S2", "S3")])-sum(output[nrow(output),c("S1", "S2", "S3")]))
print(baseline_results)## child_cum_inc adult_cum_inc elderly_cum_inc total_cum_inc
## 1 190215.4 609095.5 109311.1 908621.9print(results4)## child_cum_inc adult_cum_inc elderly_cum_inc total_cum_inc
## 1 130785 412563.8 77407.67 620756.4# Reduction in prevalence achieved with vaccination:
(baseline_results-results4)/baseline_results## child_cum_inc adult_cum_inc elderly_cum_inc total_cum_inc
## 1 0.3124374 0.3226615 0.2918589 0.3168155Actually, the percentage reduction in prevalence achieved with this vaccination strategy is very similar across the 3 age groups! It is slightly higher in children and adults than in the elderly; the vaccine reduces the cumulative incidence in children and adults by 31-32%, compared to a 29% reduction in the elderly.
At first glance, it might seem counterintuitive that the reduction in incidence in the elderly is nearly as high as for children and adults, despite the vaccine efficacy and therefore the effective vaccine coverage being only half that of the other age groups. However, it makes sense when looking at the contact matrix:
print(contact_matrix)## [,1] [,2] [,3]
## [1,] 7 5 1
## [2,] 2 9 1
## [3,] 1 3 2On average, elderly people in this population make more contacts with children and adults (1+3) than with other elderly people (2) per day, which is why they benefit from a lower proportion of infected children and adults achieved with vaccination as well.
Coding a vector-borne disease (VBD) model
Code a simple model of dengue transmission, using the classic vector-borne disease (VBD) framework. In addition to the human population, VBD models also need to represent the vector to model transmission events. Use this model to explore the effect of varying vector-related parameter values on infection prevalence in humans.
Only include a single serotype within the model. The differential equations for the mosquito population dynamics are:
\[\begin{align} \frac{dS_V}{dt} &= \mu_V N_V - \frac{a b_V}{N_H} S_V I_H - \mu_V S_V \\ \frac{dI_V}{dt} &= \frac{a b_V}{N_H} S_V I_H - \mu_V I_V \end{align}\]
The differential equations for the human host population dynamics are:
\[\begin{align} \frac{dS_H}{dt} &= - \frac{a b_H}{N_H} S_H I_V \\ \frac{dI_H}{dt} &= \frac{a b_H}{N_H} S_H I_V - r I_H \\ \frac{dR_H}{dt} &= r I_H \end{align}\]
Assumptions
- host and vector population: there are no births, background deaths or disease-induced mortality in the host population. Vectors enter the population and die at the same rate, therefore the mosquito population size remains constant over time. Vector survival is independent of infection status, and recruitment of new vectors into the transmission cycle (\(\mu_V N_V\)) depends on the vector population size. There is no heterogeneity in the host or vector population, which means vectors and bites are homogeneously distributed among people.
- transmission dynamics and biting behaviour: transmission occurs only from vector to host and from host to vector (no host-host or vector-vector transmission). The mosquito takes all blood meals from humans (H=1 so this is not represented in the equations).
- natural history of infection: hosts are infectious as soon as they are infected, and can recover. Recovery induces permanent immunity in hosts, i.e. there is no re-infection. Vectors are infectious as soon as they get infected and remain infectious until they die.
Based on information from Nishiura (2006), develop a model of a closed system of 20,000 vectors and 10,000 humans, assuming an infection prevalence of 0.28% in humans and 0.057% in mosquitoes and that everyone else is susceptible.
Parameter values
- Biting rate \(a\) = 1 days\(^{-1}\): Each mosquito bites a human once a day (all bloodmeals are taken from humans).
- Probability of infection from an infected host to a susceptible vector, \(b_V\) = 0.4: When a mosquito bites an infected human, the probability that the mosquito becomes instantaneously infected is 40%.
- Probability of infection from an infected vector to a susceptible host, \(b_H\) = 0.4: When a human gets bitten by an infectious mosquito, the human becomes infected with a probability of 40%.
- Mortality rate of the vector \(\mu_V\) = 0.25 days\(^{-1}\): the Aedes aegypti mosquito has an average life expectancy of 1/0.25 = 4 days.
- Recovery rate of the host, \(r\) = 0.167 days\(^{-1}\): dengue infection in humans lasts 1/0.167 = 6 days on average.
Based on the SIR model, code for a single-serotype dengue model and simulate for a period of 3 months. Plot the number of hosts and vectors in each compartment over this time period.
# part 1
nicesubtitle <- "SIR Model vbd1 vector-borne disease c3w3"
print(nicesubtitle)## [1] "SIR Model vbd1 vector-borne disease c3w3"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
Nh <- 10000 # number of humans
Nv <- 20000 # number of humans
I0h <- round(0.280/100 * Nh,0) # infection prevalence humans
I0v <- round(0.057/100 * Nv,0) # infection prevalence mosquitoes
S0h <- Nh - I0h # remaining humans are all susceptible
S0v <- Nv - I0v # remaining vectors are all susceptible
duration <- 90 # total number of days - 3 months
tsteps <- 1.0 # chunk into 1 day intervals
bite <- 1 # average number of bites per mosquito per day
bh = 0.4 # probability of infection from an infected vector to a susceptible host
bv = 0.4 # probability of infection from an infected host to a susceptible vector
muv = 0.25 # mortality rate of the vector
gammah = 0.167 # recovery rate of the host
# beta = infection rate day^-1 is a function of (bh or bv) * bite)
beta_h <- bh * bite
beta_v <- bv * bite
R0_h <- beta_h / gammah
(parameters <- c(
bite = bite, # mosquito biting rate per day
bh = bh, # prob of infection from infected vector to susceptible host
bv = bv, # prob of infection from infected host to susceptible vector
muv = muv, # mortality rate of the vector
gammah = gammah, # recovery rate of the host
beta_h = beta_h,
beta_v = beta_v,
R0_h = R0_h
)) ## bite bh bv muv gammah beta_h beta_v R0_h
## 1.00000 0.40000 0.40000 0.25000 0.16700 0.40000 0.40000 2.39521# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_vbd1 <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
Nh <- Sh + Ih + Rh
Nv <- Sv + Iv
dSh <- -((bite*bh/Nh) * Sh * Iv)
dIh <- ((bite*bh/Nh) * Sh * Iv) - (gammah * Ih)
dRh <- (gammah * Ih)
dSv <- (muv * Nv) - ((bite*bv/Nh) * Sv * Ih) - (muv * Sv)
dIv <- ((bite*bv/Nh) * Sv * Ih) - (muv * Iv)
return(list(c(dSh, dIh, dRh, dSv, dIv)))
})
}
# initial state values
(initial_state_values <- c(
Sh = Nh-I0h,
Ih = I0h,
Rh = 0,
Sv = Nv-I0v,
Iv = I0v
))## Sh Ih Rh Sv Iv
## 9972 28 0 19989 11# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_vbd1,
parms = parameters)) %>%
mutate(
Nh = (Sh+Ih+Rh),
Nv = (Sv+Iv),
still_Su_h = round(Sh/Nh,digits=5)*100,
preval_Inf_h = round(Ih/Nh,digits=5)*100,
propor_Re_h = round(Rh/Nh,digits=5)*100,
cumul_incid_h = round(Nh-Sh,digits=5)
, cumul_incid_h_pct = round(cumul_incid_h/Nh,digits=5)*100
,
still_Su_v = round(Sv/Nv,digits=5)*100,
preval_Inf_v = round(Iv/Nv,digits=5)*100,
cumul_incid_v = round(Nv-Sv,digits=5)
, cumul_incid_v_pct = round(cumul_incid_v/Nv,digits=5)*100
)
print("first record for the run: ")## [1] "first record for the run: "output %>%
arrange(time) %>%
head(1)## time Sh Ih Rh Sv Iv Nh Nv still_Su_h preval_Inf_h propor_Re_h
## 1 0 9972 28 0 19989 11 10000 20000 99.72 0.28 0
## cumul_incid_h cumul_incid_h_pct still_Su_v preval_Inf_v cumul_incid_v
## 1 28 0.28 99.945 0.055 11
## cumul_incid_v_pct
## 1 0.055print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time Sh Ih Rh Sv Iv Nh Nv still_Su_h
## 91 90 31.06478 0.1230327 9968.812 19998.95 1.047749 10000 20000 0.311
## preval_Inf_h propor_Re_h cumul_incid_h cumul_incid_h_pct still_Su_v
## 91 0.001 99.688 9968.935 99.689 99.995
## preval_Inf_v cumul_incid_v cumul_incid_v_pct
## 91 0.005 1.04775 0.005print(paste0("Cumulative incidence (host): ",
round(output$cumul_incid_h[nrow(output)], 0), " (",
round((output$cumul_incid_h[nrow(output)]/Nh), 3)*100, " %)"))## [1] "Cumulative incidence (host): 9969 (99.7 %)"print(paste0("Cumulative incidence (vector): ",
round(output$cumul_incid_v[nrow(output)], 0), " (",
round((output$cumul_incid_v[nrow(output)]/Nv), 3)*100, " %)"))## [1] "Cumulative incidence (vector): 1 (0 %)"# Extract cumulative incidence for each run:
(results <- data.frame(
h_cumul_inc = output$cumul_incid_h[nrow(output)],
v_cumul_inc = output$cumul_incid_v[nrow(output)])
)## h_cumul_inc v_cumul_inc
## 1 9968.935 1.04775# part 2 plot
print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"(dengueplot <- output %>%
select(time, Sh, Ih, Rh, Sv, Iv) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green", "deepskyblue3", "deeppink3")) +
scale_shape_manual(values = c(0,4,1,"-",".")) +
xlab("Time (days)") +
ylab("Number") +
labs(title=paste(
" bh of", parameters['bh'],
" bv of", parameters['bv'],
" gammah of", parameters['gammah']
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") )
Sensitivity of infection prevalence to the biting rate
Theoretical studies of the behaviour of mathematical models under different parameter assumptions have had an important contribution to our understanding of VBD transmission dynamics and control. Perform a simple univariate sensitivity analysis to explore how varying values of the biting rate affect infection prevalence in humans in our modelled dengue outbreak, keeping all other parameters constant.
The easiest way of visualising this is to plot the infection prevalence over time for given values of the parameter in a single plot. Use a for loop to simulate the model for a set of biting rate values.
# part 3 multiple infection prevalence values comparative plot
print("Plotting multiple infection prevalence values")## [1] "Plotting multiple infection prevalence values"bite_values <- seq(0,1,by=0.1) # a vector of values for the biting rate, ranging from 0 to 1 per day
out_list <- vector(length(bite_values), mode = "list") # store output for each biting rate value
for (i in seq_along(bite_values)) {
out_list[[i]] <- as.data.frame(ode(
y = initial_state_values,
times = times,
func = sir_model_vbd1,
parms = c(parameters[c("bv", "bh", "muv", "gammah")],
bite = bite_values[i])))
}
names(out_list) <- bite_values # rename
# Extract the infected host column from the list and creating a dataframe by time
out_inf <- as_tibble(cbind(time = out_list[[1]]$time,
sapply(out_list,
"[[",
"Ih")))
# plot
out_inf %>%
melt(id = "time") %>%
# Plot the infection prevalence for each biting rate value
ggplot(aes(x = time, y = value, colour = variable, group = variable)) +
geom_line() +
xlab("Time (days)") +
ylab("Human infection prevalence") +
labs(title=paste(
" bh of", parameters['bh'],
" bv of", parameters['bv'],
" gammah of", parameters['gammah']
),
color = "Biting rate \n per day",
subtitle= nicesubtitle)
How do assumptions of the mosquito biting activity affect the human infection prevalence?
For an average number of 0-0.3 bites per mosquito per day, no dengue epidemic occurs. An outbreak only occurs for values of 0.4 and above, with increasing values of \(a\) producing higher peak epidemic sizes and shorter durations.
What might affect the average number of bites a mosquito takes per unit time? Is it realistic to assume that the biting rate stays constant over time in the simulation?
The biting rate, just like other vector-related parameters, can be influenced by a wide range of environmental factors such the temperature and rainfall. Over a short time period, assuming a constant biting rate may be fine, but accounting for seasonality is important over longer timescales (as is turnover in the human host population). The biting rate can also be affected by the host’s behaviour, for example if humans take measures to avoid mosquito bites (e.g. wearing long-sleeve clothing).
How could you extend this model for a more realistic representation of dengue transmission dynamics?
This would depend on the research question, but in general a more realistic model structure could be achieved by representing mosquito population dynamics and the natural history of dengue infection. For example, dengue models usually account for the intrinsic incubation period (latent period in humans) and extrinsic incubation period (latent period in the vector), i.e. the time between exposure to infection and the bite becoming infectious. Only adult mosquitoes are involved in transmission, but different developmental stages from larvae to adult could be represented, particularly if we are interested in studying changes in the mosquito population e.g. in relation to environmental factors. Severe cases of dengue infection occur as a result of secondary infection with a different serotype, so the co-circulation of different serotypes would be important to account for when studying disease burden of dengue in humans.
Coding vector control interventions
Start thinking about the different parameters of the vectorial capacity and how they relate to interventions for vector control.
Vectorial capacity questions
Part 1: Interpreting the vectorial capacity equation
\[\begin{align} V = \frac{m a^2 p^n}{-ln(p)} \end{align}\]
where \(m\) is the ratio of vectors to hosts, \(a\) is the daily biting rate, \(p\) is the probability of the vector surviving through one day and \(n\) is the average duration of the extrinsic incubaton period in days.
Question 1: Can you think of a reason why the biting rate is present more than once in the equation (i.e. it is squared)?
The biting rate is counted twice in the equation because the vector has to take 2 bites to transmit - once to become infected, and a second time to pass it on.
Question 2: How do changes in the vector-to-host ratio change a) the rate at which vectors bite, and b) the rate at which humans are bitten?
Changes in the vector-to-host ratio do not affect a) the rate at which vectors bite, which is defined in a separate parameter (the biting rate \(a\)) and represents the number of bloodmeals taken by each vector per unit time on average, irrespective of how many vectors there are and independent of the given host (humans). However, increasing the vector-to-host ratio will also increase b) the rate at which human hosts are bitten, which is the product \(ma\) (units of bites per day).
Reality is of course more complex, and changes in these parameters are often interacting, so that the biting rate on humans may also depend on changes in host or vector density.
Question 3: What is \(p^n\)?
\(p\) is the probability of the vector surviving through one day and \(n\) is the average duration of the extrinsic incubation period in days, so \(p^n\) is the probability of the vector surviving through the extrinsic incubation period.
Question 4: Why is the duration of the latent period, \(n\), important to account for in models of mosquito-borne diseases?
Because mosquitoes live only for a short time, the latent period of the pathogen is often a considerable proportion of the mosquito’s life expectancy. Therefore, changes in both the life expectancy and incubation period of the pathogen in the vector can have important effects on transmission.
Part 2: Implications for vector control
For Aedes aegypti, the main mosquito transmitting dengue virus, estimated parameter values in the Bello municipality in Colombia are (Catano-Lopez, 2019):
\(m\) = 3, \(a\) = 0.29 days-1, \(p\) = 0.98, \(n\) = 6 days
This gives a vectorial capacity of around 11.
The plots below shows how vectorial capacity changes for varying values of each parameter, while keeping the other parameters constant. The vertical line indicates the Aedes aegypti parameter assumptions above.
# VC calc function
VC_calc <- function(m, bite, p, n) {
VC <- (m * (bite**2) * (p**n))/(-1 * log(p))
return(VC)
}
# VC plot function
VC_plot <- function(nicesubtitle, out_inf, param_values, param_name, param_ref) {
VCplot <- out_inf %>%
melt(id = "param_values") %>%
ggplot(aes(x = param_values, y = value, colour = param_values)) +
geom_line(show.legend = FALSE) +
xlab(" ") +
ylab("Vectorial Capacity") +
labs(title=paste(
param_name,
": ", min(param_values),
" - ", max(param_values)
),
color = param_name,
subtitle= nicesubtitle) +
geom_vline(xintercept=param_ref, linetype="dotted", color="blue")
return(VCplot)
}
# Vectorial Capacity changes
nicesubtitle <- "Vectorial Capacity changes for varying parameter values c3w3"
print(nicesubtitle)## [1] "Vectorial Capacity changes for varying parameter values c3w3"print("initial state values and parameters")## [1] "initial state values and parameters"m <- 3 # ratio of vectors to hosts
bite <- 0.29 # daily biting rate
p <- .98 # probability of the vector surviving through 1 day
n <- 6 # average duration of the extrinsic incubation period, in days
VC <- (m * (bite**2) * (p**n))/(-1 * log(p)) # Vectorial Capacity 11.06278
# plot 1 m Vector-to-host ratio
param_name <- "Vector-to-host ratio"
print(param_name)## [1] "Vector-to-host ratio"param_values <- seq(0, 6 , by= 1)
out_list <- vector(length(param_values), mode = "list")
for (i in seq_along(param_values)) {
out_list[[i]] <- VC_calc(param_values[i], bite, p, n)
}
names(out_list) <- param_values
out_inf <- as_tibble(cbind(param_values, VC=as.numeric(out_list)))
# plot
VC_plot(nicesubtitle, out_inf, param_values, param_name, m)
# plot 2 bite a daily biting rate
param_name <- "Daily biting rate"
print(param_name)## [1] "Daily biting rate"param_values <- seq(0, .6 , by= .1)
out_list <- vector(length(param_values), mode = "list")
for (i in seq_along(param_values)) {
out_list[[i]] <- VC_calc(m, param_values[i], p, n)
}
names(out_list) <- param_values
out_inf <- as_tibble(cbind(param_values, VC=as.numeric(out_list)))
# plot
VC_plot(nicesubtitle, out_inf, param_values, param_name, bite)
# plot 3 p prob of survival through 1 day
param_name <- "Probability of survival through 1 day"
print(param_name)## [1] "Probability of survival through 1 day"param_values <- seq(0, 1 , by= .02)
out_list <- vector(length(param_values), mode = "list")
for (i in seq_along(param_values)) {
out_list[[i]] <- VC_calc(m, bite, param_values[i], n)
}
names(out_list) <- param_values
out_inf <- as_tibble(cbind(param_values, VC=as.numeric(out_list)))
# Replace Inf values with NA for plotting
out_inf$VC <- (ifelse(!is.infinite(out_inf$VC),
out_inf$VC, NA))
# plot
VC_plot(nicesubtitle, out_inf, param_values, param_name, p)## Warning: Removed 1 row(s) containing missing values (geom_path).
# plot 4 n average duration of the extrinsic incubation
param_name <- "Extrinsic incubation period in days"
print(param_name)## [1] "Extrinsic incubation period in days"param_values <- seq(0, 50 , by= 10)
out_list <- vector(length(param_values), mode = "list")
for (i in seq_along(param_values)) {
out_list[[i]] <- VC_calc(m, bite, p, param_values[i])
}
names(out_list) <- param_values
out_inf <- as_tibble(cbind(param_values, VC=as.numeric(out_list)))
# plot
VC_plot(nicesubtitle, out_inf, param_values, param_name, n)
Question 5: Keeping all other parameters constant, by how much (approximately) would you have to change each parameter to reduce the vectorial capacity by half?
Based on the plots (or the equation), to reduce the vectorial capacity to around 5.5 while keeping all other parameter values constant, we would need to: - reduce \(m\) by 50% (linear relationship between \(m\) and \(V\)) - reduce \(a\) from 0.29 to 0.2 (31% reduction) - reduce \(p\) from 0.98 to 0.965 (1.5% reduction) - increase \(n\) from 6 to 40 days (over 600% increase).
Question 6: Based on your answer above, which two parameters would you preferentially target with vector control interventions? In both cases, give an example of an appropriate intervention.
The vectorial capacity is most sensitive to changes in the probability of survival through one day, \(p\), and the biting rate, \(a\): we only need to achieve small reductions in these parameters to achieve a large reduction in the vectorial capacity. Therefore, these might represent the most efficient targets for vector control. Mosquito survival can be reduced using insecticides targeting the adult stage or Wolbachia, whereas the biting rate can be reduced through the use of personal protective interventions (e.g. long-sleeved clothes, insect repellent).
Question 7: When studying the impact of vector control interventions on dengue, what else would we need to account for?
The vectorial capacity only captures the future number of potentially infectious mosquito bites through properties of the mosquito and relies on a number of assumptions (see the Dye 1986 review for details on assumptions). To study the impact of vector control interventions on the transmission of dengue, we would also need to account for the biology and epidemiology of dengue infection in humans, as well as the efficacy and properties of any available interventions, including compliance.
Modelling vector control interventions
The following cover two different ways of modelling implementation of interventions during an outbreak. Previously, we have only represented interventions like vaccination and treatment to be implemented from the start of the simulation (before the beginning of the outbreak). In reality, an infection is often already established or an outbreak has begun before control measures are taken. The following code is for the dengue model but with the addition of a vector control intervention.
Intervention 1 represents an example of an intervention which changes the value of a parameter over time, and Intervention 2 is an example of an intervention which instantaneously affects the population at a specific timepoint. For comparison, here is the output from the dengue model without interventions from the previous activity:
dengueplot +
labs(title="Dengue model")
Intervention 1:
# Intervention 1
nicesubtitle <- "SIR Model vbd2 Dengue Interv1 over time c3w3"
print(nicesubtitle)## [1] "SIR Model vbd2 Dengue Interv1 over time c3w3"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
Nh <- 10000 # number of humans
Nv <- 20000 # number of humans
I0h <- round(0.280/100 * Nh,0) # infection prevalence humans
I0v <- round(0.057/100 * Nv,0) # infection prevalence mosquitoes
S0h <- Nh - I0h # remaining humans are all susceptible
S0v <- Nv - I0v # remaining vectors are all susceptible
duration <- 90 # total number of days - 3 months
tsteps <- 1.0 # chunk into 1 day intervals
bite <- 1 # average number of bites per mosquito per day
bh = 0.4 # probability of infection from an infected vector to a susceptible host
bv = 0.4 # probability of infection from an infected host to a susceptible vector
muv = 0.25 # mortality rate of the vector
gammah = 0.167 # recovery rate of the host
# beta = infection rate day^-1 is a function of (bh or bv) * bite)
beta_h <- bh * bite
beta_v <- bv * bite
R0_h <- beta_h / gammah
(parameters <- c(
bite = bite, # mosquito biting rate per day
bh = bh, # prob of infection from infected vector to susceptible host
bv = bv, # prob of infection from infected host to susceptible vector
muv = muv, # mortality rate of the vector
gammah = gammah, # recovery rate of the host
beta_h = beta_h,
beta_v = beta_v,
R0_h = R0_h
)) ## bite bh bv muv gammah beta_h beta_v R0_h
## 1.00000 0.40000 0.40000 0.25000 0.16700 0.40000 0.40000 2.39521# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_vbd2 <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
Nh <- Sh + Ih + Rh
Nv <- Sv + Iv
if (time <= 15){
bite = 1
} else if (time > 15 & time <= 60){
bite = 0.25
} else if (time > 60) {
bite = 1
}
dSh <- -((bite*bh/Nh) * Sh * Iv)
dIh <- ((bite*bh/Nh) * Sh * Iv) - (gammah * Ih)
dRh <- (gammah * Ih)
dSv <- (muv * Nv) - ((bite*bv/Nh) * Sv * Ih) - (muv * Sv)
dIv <- ((bite*bv/Nh) * Sv * Ih) - (muv * Iv)
return(list(c(dSh, dIh, dRh, dSv, dIv)))
})
}
# initial state values
(initial_state_values <- c(
Sh = Nh-I0h,
Ih = I0h,
Rh = 0,
Sv = Nv-I0v,
Iv = I0v
))## Sh Ih Rh Sv Iv
## 9972 28 0 19989 11# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_vbd2,
parms = parameters)) %>%
mutate(
Nh = (Sh+Ih+Rh),
Nv = (Sv+Iv),
still_Su_h = round(Sh/Nh,digits=5)*100,
preval_Inf_h = round(Ih/Nh,digits=5)*100,
propor_Re_h = round(Rh/Nh,digits=5)*100,
cumul_incid_h = round(Nh-Sh,digits=5)
, cumul_incid_h_pct = round(cumul_incid_h/Nh,digits=5)*100
,
still_Su_v = round(Sv/Nv,digits=5)*100,
preval_Inf_v = round(Iv/Nv,digits=5)*100,
cumul_incid_v = round(Nv-Sv,digits=5)
, cumul_incid_v_pct = round(cumul_incid_v/Nv,digits=5)*100
)
print("first record for the run: ")## [1] "first record for the run: "output %>%
arrange(time) %>%
head(1)## time Sh Ih Rh Sv Iv Nh Nv still_Su_h preval_Inf_h propor_Re_h
## 1 0 9972 28 0 19989 11 10000 20000 99.72 0.28 0
## cumul_incid_h cumul_incid_h_pct still_Su_v preval_Inf_v cumul_incid_v
## 1 28 0.28 99.945 0.055 11
## cumul_incid_v_pct
## 1 0.055print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time Sh Ih Rh Sv Iv Nh Nv still_Su_h
## 91 90 1759.371 987.4055 7253.224 17573.34 2426.656 10000 20000 17.594
## preval_Inf_h propor_Re_h cumul_incid_h cumul_incid_h_pct still_Su_v
## 91 9.874 72.532 8240.629 82.406 87.867
## preval_Inf_v cumul_incid_v cumul_incid_v_pct
## 91 12.133 2426.656 12.133print(paste0("Cumulative incidence (host): ",
round(output$cumul_incid_h[nrow(output)], 0), " (",
round((output$cumul_incid_h[nrow(output)]/Nh), 3)*100, " %)"))## [1] "Cumulative incidence (host): 8241 (82.4 %)"print(paste0("Cumulative incidence (vector): ",
round(output$cumul_incid_v[nrow(output)], 0), " (",
round((output$cumul_incid_v[nrow(output)]/Nv), 3)*100, " %)"))## [1] "Cumulative incidence (vector): 2427 (12.1 %)"# Extract cumulative incidence for each run:
(results <- data.frame(
h_cumul_inc = output$cumul_incid_h[nrow(output)],
v_cumul_inc = output$cumul_incid_v[nrow(output)])
)## h_cumul_inc v_cumul_inc
## 1 8240.629 2426.656# part 2 plot
print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"(dengueplot2 <- output %>%
select(time, Sh, Ih, Rh, Sv, Iv) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green", "deepskyblue3", "deeppink3")) +
scale_shape_manual(values = c(0,4,1,"-",".")) +
xlab("Time (days)") +
ylab("Number") +
labs(title=paste(
" bh of", parameters['bh'],
" bv of", parameters['bv'],
" gammah of", parameters['gammah']
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") )
1a) Based on the code, what effect does the intervention have on the natural history of the infection and how effective is it at this?
The modelled interventions works by reducing the biting rate \(bite\) during a given time period (between days 15 and 60). It has a 75% efficacy at reducing the biting rate (from 1 per day to 0.25 per day).
1b) Give an example of an intervention this might represent and describe in more detail how it is implemented (e.g. when, for how long).
This might represent implementation of a personal protective intervention on the population level, e.g. rolling out an insect repellent scheme to all individuals in the modelled community, which would reduce the biting rate on humans. At the beginning of the simulation and until the timestep corresponding to day 15, \(bite\)=1 days\(^{-1}\). From the plot we can see that a dengue outbreak has started by that time. At day 15, the insect repellent is distributed and assumed to reduce the biting rate to 0.25 per day. Use of insect repellent is maintained for 45 days, but stops at day 60 so the biting rate increases again to its baseline value.
1c) Based on the plot, how does the intervention affect the prevalence of infection in humans?
Compared to the scenario without intervention, the peak of the dengue outbreak occurs slightly earlier and is lower. However, the output at the later timesteps also shows that as soon as the intervention is no longer in effect, the number of infections in humans starts rising again.
Notes:
The dengue vector Aedes aegypti is a day-biting mosquito, which makes personal protective interventions like the one modelled here an important method of prevention. However, as always this simple model makes many simplifying assumptions, such as that uptake and compliance to the use of insect repellent is homogenous across the population.
While we are looking at vector control interventions against dengue infection, the general principle is the same for other types of diseases and control measures. The biting rate in this model, which changes as a result of an intervention, is an example of a time-dependent parameter.
Intervention 2
# Intervention 2
nicesubtitle <- "SIR Model vbd1b Dengue Interv2 event instant effect c3w3"
print(nicesubtitle)## [1] "SIR Model vbd1b Dengue Interv2 event instant effect c3w3"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
Nh <- 10000 # number of humans
Nv <- 20000 # number of humans
I0h <- round(0.280/100 * Nh,0) # infection prevalence humans
I0v <- round(0.057/100 * Nv,0) # infection prevalence mosquitoes
S0h <- Nh - I0h # remaining humans are all susceptible
S0v <- Nv - I0v # remaining vectors are all susceptible
duration <- 90 # total number of days - 3 months
tsteps <- 1.0 # chunk into 1 day intervals
bite <- 1 # average number of bites per mosquito per day
bh = 0.4 # probability of infection from an infected vector to a susceptible host
bv = 0.4 # probability of infection from an infected host to a susceptible vector
muv = 0.25 # mortality rate of the vector
gammah = 0.167 # recovery rate of the host
# beta = infection rate day^-1 is a function of (bh or bv) * bite)
beta_h <- bh * bite
beta_v <- bv * bite
R0_h <- beta_h / gammah
(parameters <- c(
bite = bite, # mosquito biting rate per day
bh = bh, # prob of infection from infected vector to susceptible host
bv = bv, # prob of infection from infected host to susceptible vector
muv = muv, # mortality rate of the vector
gammah = gammah, # recovery rate of the host
beta_h = beta_h,
beta_v = beta_v,
R0_h = R0_h
)) ## bite bh bv muv gammah beta_h beta_v R0_h
## 1.00000 0.40000 0.40000 0.25000 0.16700 0.40000 0.40000 2.39521# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_vbd1 <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
Nh <- Sh + Ih + Rh
Nv <- Sv + Iv
dSh <- -((bite*bh/Nh) * Sh * Iv)
dIh <- ((bite*bh/Nh) * Sh * Iv) - (gammah * Ih)
dRh <- (gammah * Ih)
dSv <- (muv * Nv) - ((bite*bv/Nh) * Sv * Ih) - (muv * Sv)
dIv <- ((bite*bv/Nh) * Sv * Ih) - (muv * Iv)
return(list(c(dSh, dIh, dRh, dSv, dIv)))
})
}
# initial state values
(initial_state_values <- c(
Sh = Nh-I0h,
Ih = I0h,
Rh = 0,
Sv = Nv-I0v,
Iv = I0v
))## Sh Ih Rh Sv Iv
## 9972 28 0 19989 11eventdata <- data.frame(var = c("Sv","Iv"), time = c(15,15),
value = c(0.5,0.5), method = "multiply")
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_vbd1,
parms = parameters,
events = list(data = eventdata))) %>%
mutate(
Nh = (Sh+Ih+Rh),
Nv = (Sv+Iv),
still_Su_h = round(Sh/Nh,digits=5)*100,
preval_Inf_h = round(Ih/Nh,digits=5)*100,
propor_Re_h = round(Rh/Nh,digits=5)*100,
cumul_incid_h = round(Nh-Sh,digits=5)
, cumul_incid_h_pct = round(cumul_incid_h/Nh,digits=5)*100
,
still_Su_v = round(Sv/Nv,digits=5)*100,
preval_Inf_v = round(Iv/Nv,digits=5)*100,
cumul_incid_v = round(Nv-Sv,digits=5)
, cumul_incid_v_pct = round(cumul_incid_v/Nv,digits=5)*100
)
print("first record for the run: ")## [1] "first record for the run: "output %>%
arrange(time) %>%
head(1)## time Sh Ih Rh Sv Iv Nh Nv still_Su_h preval_Inf_h propor_Re_h
## 1 0 9972 28 0 19989 11 10000 20000 99.72 0.28 0
## cumul_incid_h cumul_incid_h_pct still_Su_v preval_Inf_v cumul_incid_v
## 1 28 0.28 99.945 0.055 11
## cumul_incid_v_pct
## 1 0.055print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time Sh Ih Rh Sv Iv Nh Nv still_Su_h
## 91 90 418.6724 1.7317 9579.596 9994.812 5.188421 10000 10000 4.187
## preval_Inf_h propor_Re_h cumul_incid_h cumul_incid_h_pct still_Su_v
## 91 0.017 95.796 9581.328 95.813 99.948
## preval_Inf_v cumul_incid_v cumul_incid_v_pct
## 91 0.052 5.18842 0.052print(paste0("Cumulative incidence (host): ",
round(output$cumul_incid_h[nrow(output)], 0), " (",
round((output$cumul_incid_h[nrow(output)]/Nh), 3)*100, " %)"))## [1] "Cumulative incidence (host): 9581 (95.8 %)"print(paste0("Cumulative incidence (vector): ",
round(output$cumul_incid_v[nrow(output)], 0), " (",
round((output$cumul_incid_v[nrow(output)]/Nv), 3)*100, " %)"))## [1] "Cumulative incidence (vector): 5 (0 %)"# Extract cumulative incidence for each run:
(results <- data.frame(
h_cumul_inc = output$cumul_incid_h[nrow(output)],
v_cumul_inc = output$cumul_incid_v[nrow(output)])
)## h_cumul_inc v_cumul_inc
## 1 9581.328 5.18842# part 2 plot
print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"(dengueplot3 <- output %>%
select(time, Sh, Ih, Rh, Sv, Iv) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green", "deepskyblue3", "deeppink3")) +
scale_shape_manual(values = c(0,4,1,"-",".")) +
xlab("Time (days)") +
ylab("Number") +
labs(title=paste(
" bh of", parameters['bh'],
" bv of", parameters['bv'],
" gammah of", parameters['gammah']
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") )
2a) Based on the code, what effect does the intervention have on the natural history of the infection and how effective is it at this?
Here, an event, defined in the “eventdata” dataframe, is implemented while solving the differential equations (within the ode() command).
“An event occurs when the value of a state variable is suddenly changed, e.g. because a value is added, subtracted, or multiplied. The integration routines cannot deal easily with such state variable changes. Typically these events occur only at specific times. In deSolve, events can be imposed by means of an input data.frame, that specifies at which time and how a certain state variable is altered, or via an event function.”
The “eventdata” dataframe tells us that we are changing the number of susceptible vectors (Sv) and the number of infected vectors (Iv), both at timestep 15. We are multiplying the number in those states at that timepoint by 0.5, which means we are modelling an intervention with 50% efficacy at reducing the vector population in the community.
2b) Give an example of an intervention this might represent and describe in more detail how it is implemented (e.g. when, for how long).
This might represent a fogging intervention, whereby widescale spraying of an insecticide instantly kills a large number of mosquitoes. In the code, this happens 15 days into the simulation, at which point a dengue outbreak has started. In contrast to the continued use of the insect repellent, fogging represents a one-off event in this example.
2c) Based on the plot, how does the intervention affect the prevalence of infection in humans?
The plot shows the abrupt drop in the number of susceptible and infected vectors at day 15, at which point the infection prevalence in humans also starts declining. The peak of the human dengue epidemic is thereby reduced compared to the scenario with no intervention.
Note:
The plot also shows that right after the fogging event, the mosquito population slowly increases again, as new mosquitoes are recruited into the population through maturation of larvae (as defined in the \(\mu_v\) parameter). This is because the mosquito larvae are not killed by fogging.