DRAFT Notes: Developing the SIR Model
Infectious Disease Modelling Specialization, Imperial College London
analysisinsightdata
Thursday, September 09, 2021
Reference: Greek Letters et al.
S = number of susceptible individuals in the population
I = number of infected individuals in the population
R = number of recovered/removed individuals in the population
N = total number of individuals in the population (S + I + R)
c = average number of contacts a susceptible makes per time
p = probability of a susceptible becoming infected on contact with infected person
beta β = infection rate (function of c * p)
gamma γ = recovery rate
lambda λ = rate of movement from exposed to infectious (beta * I/N)
R0 R-naught = basic reproduction number of infectious agents (beta / gamma)
Reff Rt = Effective Reproduction Number (R0 * S/N)
pvacc p = vaccination rate
pvacc_threshold pc = critical vaccination threshold (1 - (1/R0))
sigma σ = waning immunity rate
mu μ = background mortality rate
birth b = immigration or birth rate (same as mu if stable population)
M = number of deceased (if used)
—
Modelling an infected cohort
Find out how long it takes for a cohort of infected people to recover. You need to keep track of 2 populations: those that are infected (compartment \(I\)), and those that have recovered (compartment \(R\)). Infected people recover at a rate \(\gamma\) (gamma). The differential equations describing this are:
\[\begin{align} \frac{dI}{dt} & = -\gamma I \\ \frac{dR}{dt} & = \gamma I \end{align}\]
We are looking at a cohort of 10\(^6\) currently infected people, and no one has recovered so far. The average duration of infection is 10 days. The question we want to answer is how many people will recover from the infection over a 4-week period.
nicesubtitle <- "Modelling an infected cohort"
print(nicesubtitle)## [1] "Modelling an infected cohort"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_number_infected <- 1000000 # the initial infected population
# size
initial_number_recovered <- 0 # the initial number of people in
# the recovered state
recovery_rate <- 10 ** (-1) # the rate of recovery gamma,
# in units of days^-1
follow_up_duration <- 7 * 4 # the duration to run the model for,
# in units of days
# The initial state values are stored as a vector
# and each value is assigned a name.
initial_state_values <- c(I = initial_number_infected,
R = initial_number_recovered)
# Parameters are also stored as a vector with assigned names and values.
parameters <- c(gamma = recovery_rate)
# In this case we only have one parameter, gamma.
# The times vector creates a sequence of timepoints at which we want to
# calculate the number of people in the I and R compartment.
times <- seq(from = 0, to = follow_up_duration, by = 1)
initial_state_values## I R
## 1e+06 0e+00parameters## gamma
## 0.1times## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## [26] 25 26 27 28cohort_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
dI <- -(parameters['gamma']) * state['I']
dR <- parameters['gamma'] * state['I']
return(list(c(dI, dR)))
})
}
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = cohort_model,
parms = parameters))Question: Based on the output, how many people have recovered after 4 weeks? What proportion of the total population does this correspond to?
print("recovered")## [1] "recovered"output['R'][29,]## [1] 939189.9print("proportion of total population")## [1] "proportion of total population"output['R'][29,] / initial_number_infected## [1] 0.9391899Plot your model with time on the x axis and the number of infected and recovered people on the y axis.
# First turn the output dataset into a long format,
# so that the number in each compartment at each timestep
# are all in the same column
output_long <- melt(as.data.frame(output), id = "time")
# Plot the number of people in each compartment over time
ggplot(data = output_long, # specify object containing data to plot
aes(x=time, y=value,color=variable)) + # assign columns to axes and groups
geom_line() + # represent data as lines
xlab("Time (days)")+ # add label for x axis
ylab("Number of people")+ # add label for y axis
labs(title=paste("gamma of",recovery_rate))
Question: Based on the plot, at what timepoint were infected and recovered individuals equal in number?
print("roughly right before 7 days")## [1] "roughly right before 7 days"output %>%
mutate(differs = abs(I - R),
proportion = R/(I + R)) %>%
arrange(differs, time) %>%
head(1)## time I R differs proportion
## 1 7 496585.3 503414.7 6829.4 0.5034147Vary \(\gamma\) to see how it affects the output. Change gamma (\(\gamma\)) to correspond to an average infectious period of:
a) 2 days
b) 20 days.
What is the recovery rate in these 2 cases?
print("average infectious period of 2 days")## [1] "average infectious period of 2 days"(recovery_rate <- 2 ** (-1))## [1] 0.5parameters <- c(gamma = recovery_rate)
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = cohort_model,
parms = parameters))
# Since the initial number in each compartment and
# the timesteps of interest haven't changed,
# these are the only parts of the code we need to rerun.
# Now, copy-paste your plot code from above here to visualise the output.
melt(as.data.frame(output), id = "time") %>%
ggplot(aes(x=time, y=value,color=variable)) +
geom_line() +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste("gamma of",recovery_rate))
print("closest to all recovered")## [1] "closest to all recovered"output %>%
mutate(differs = abs(I - R),
proportion = R/(I + R)) %>%
arrange(-R, time) %>%
head(1)## time I R differs proportion
## 1 28 0.8315196 999999.2 999998.3 0.9999992print("average infectious period of 20 days")## [1] "average infectious period of 20 days"(recovery_rate <- 20 ** (-1))## [1] 0.05parameters <- c(gamma = recovery_rate)
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = cohort_model,
parms = parameters))
# Since the initial number in each compartment and
# the timesteps of interest haven't changed,
# these are the only parts of the code we need to rerun.
# Now, copy-paste your plot code from above here to visualise the output.
melt(as.data.frame(output), id = "time") %>%
ggplot(aes(x=time, y=value,color=variable)) +
geom_line() +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste("gamma of",recovery_rate))
print("closest to all recovered")## [1] "closest to all recovered"output %>%
mutate(differs = abs(I - R),
proportion = R/(I + R)) %>%
arrange(-R, time) %>%
head(1)## time I R differs proportion
## 1 28 246597 753403 506806.1 0.753403Question: What changes do you observe in the transition to the recovered compartment if \(\gamma\) is higher or lower? For example, how long does it take for everyone to recover in both cases?
- With a gamma of 0.5 (recovery of 2 days), virtually everyone (99.9%) has recovered by the 28th day.
- With a gamma of 0.05 (recovery of 20 days), only about 75% have recovered by the 28th day.
—
Simulating competing hazards
We will focus on adding disease-induced mortality to this model in order to explore the concept of competing hazards. You will also calculate the case fatality ratio, and compare it against the result.
The model we want to specify has 3 compartments: \(I\) (infected), \(R\) (recovered) and \(M\) (dead).
The differential equations for the model look like this: \[\begin{align} \frac{dI}{dt} & = -\gamma I -\mu I \\ \frac{dR}{dt} & = \gamma I \\ \frac{dM}{dt} & = \mu I \end{align}\]
Question: what do \(\gamma\) (gamma) and \(\mu\) (mu) represent?
- \(\gamma\) (gamma) represents the recovery rate (people who recover).
- \(\mu\) (mu) is often used to denote mortality rate (people who have died).
- Case fatality rate is mu divided by mu + gamma.
- The proportion who people who recover and survive is given by gamma divided by mu + gamma.
This corresponds to the following model diagram:
The equation describing the rate of change in the recovered (\(R\)) compartment (second line) is not affected by this addition. However, we need a new equation describing the rate of change in the deceased (\(M\)) compartment (third line). People move into this compartment from the infected compartment (\(I\)) at a rate \(\mu\) - this transition also needs to be added in the rate of change in the infected compartment \(I\) (first line).
Incorporate the new compartment and transition. At the start, there are 10\(^6\) infected people. No one has recovered or died yet. The recovery rate \(\gamma\) is 0.1 days\(^{-1}\) and the mortality rate \(\mu\) is 0.2 days\(^{-1}\). We want to model the course of the infection over 4 weeks.
Question: after 4 weeks, do you expect more people to have recovered or more people to have died, and why?
- Because the rates of death exceeds the rate of recovery, I expect more people to die than recover.
nicesubtitle <- "Simulating competing hazards v1"
print(nicesubtitle)## [1] "Simulating competing hazards v1"print("initial state values and parameters")## [1] "initial state values and parameters"initial_state_values <- c(I = 1000000,
R = 0,
M = 0)
parameters <- c(gamma = 0.1,
mu = 0.2)
times <- seq(from = 0, to = (7 * 4), by = 1)
# Look back at the code in the previous activity if you cannot
# remember how to define these vectors.
cohort_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
dI <- (-gamma * I) - (mu * I)
dR <- gamma * I
dM <- mu * I
return(list(c(dI, dR, dM)))
})
}
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = cohort_model,
parms = parameters))
melt(as.data.frame(output), 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","green","black")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste("gamma of", parameters['gamma'],
" mu of", parameters['mu']),color = "Compartment",
subtitle= nicesubtitle)
Question: based on the model output, what proportion of the initially infected cohort died before recovering over the 4 week period?
- About 67%.
# Answer: deaths/initial number ~approx .67
output[output$time == 28,"M"] / output[output$time == 0,"I"]## [1] 0.6665168Question: now use the competing hazards formula to calculate the case fatality rate. Does this agree with your answer to the previous question?
- Yes, it is about 67%.
# Answer: Number is close
# Case fatality rate is mu divided by mu + gamma ~approx .67
parameters['mu']/(parameters['mu']+parameters['gamma'])## mu
## 0.6666667Question: Which value of \(\mu\) do you need to get a case fatality rate of 50% assuming \(\gamma\) stays fixed?
- 0.1 which is derived from:
.5 = mu/(mu + gamma);
.5mu + .5gamma = mu;
.5gamma = mu -.5mu;
.5gamma = .5mu;
gamma = mu# Which value of mu do you need to get a case fatality rate of 50% assuming gamma stays fixed?
# Answer: .5 = mu/(mu + gamma);
# .5mu + .5gamma = mu;
# .5gamma = mu -.5mu;
# .5gamma = .5mu;
# gamma = mu
parameters['mu'] <- parameters['gamma']
parameters## gamma mu
## 0.1 0.1Simulate the model using using the \(\mu\) value that you have just calculated, and verify that the code gives you a CFR of 50%.
nicesubtitle <- "Simulating competing hazards v2"
print(nicesubtitle)## [1] "Simulating competing hazards v2"output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = cohort_model,
parms = parameters))
melt(as.data.frame(output), 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","green","black")) +
scale_shape_manual(values = c(0, 4, 1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste("gamma of", parameters['gamma'],
" mu of", parameters['mu']),color = "Compartment",
subtitle= nicesubtitle)
print("deaths/initial number ~approx .498")## [1] "deaths/initial number ~approx .498"output[output$time == 28,"M"] / output[output$time == 0,"I"]## [1] 0.4981511print("Case fatality rate is mu divided by mu + gamma is .5")## [1] "Case fatality rate is mu divided by mu + gamma is .5"parameters['mu']/(parameters['mu']+parameters['gamma'])## mu
## 0.5—
SIR Model with a Constant Force of Infection
Start coding a simplified SIR model within R.
SIR model: 
The differential equations for this system are:
\[\begin{align} \frac{dS}{dt} & = -\lambda S \\ \frac{dI}{dt} & = \lambda S - \gamma I \\ \frac{dR}{dt} & = \gamma I \end{align}\]
We can use the SIR model to describe a disease that can be split into 3 states: susceptible (\(S\)), infected (\(I\)), or recovered (R). All those infected are infectious, and all those recovered are immune, so they cannot get the disease again.
One part of this model was already discussed: the transition from I to R. The new addition is the susceptibles in compartment S. Depending on how many people in the population are infectious, susceptible people experience a force of infection \(\lambda\) (lambda), which is the transition rate at which they become infected.
We are trying to simulate an outbreak of a new infectious disease that our population of 10\(^{6}\) people has not been exposed to before. This means that we are starting with a single case, everyone else is susceptible to the disease, and no one is yet immune or recovered. This can for example reflect a situation where an infected person introduces a new disease into a geographically isolated population, like on an island, or even when an infections “spill over” from other animals into a human population. In terms of the initial conditions for our model, we can define: S = 10\(^{6}\)-1 = 999999, I = 1 and R = 0.
Solve this model assuming a constant force of infection of 0.2 days\(^{-1}\). The infection has an average duration of 10 days, but we want to run the model for 60 days.
nicesubtitle <- "SIR Model with a Constant Force of Infection"
print(nicesubtitle)## [1] "SIR Model with a Constant Force of Infection"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_state_values <- c(S = (1000000 - 1),
I = 1,
R = 0)
parameters <- c(gamma = 10 **-1,
lambda = 0.2)
# TIMESTEPS:
times <- seq(from = 0, to = 60, by = 1)
# SIR MODEL FUNCTION
# We are renaming this to sir_model.
# Remember to include the input arguments,
# differential equations and output objects here.
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
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))Plot the number of people in each compartment over time. You should see that at the peak of the epidemic, around 500000 people are infected.
print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"melt(as.data.frame(output), 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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste("gamma of", parameters['gamma'],
" lambda of", parameters['lambda']),color = "Compartment",
subtitle= nicesubtitle)
Question: based on the plot, describe the pattern of the epidemic over the 2 month period. How does the number of people in the susceptible, infected and recovered compartments change over time? After how many days does the epidemic reach its peak? When does it end?
- Peak infection right about 7 days.
- End of infection visually ends at 60 (max) period with approx. 995,000 recovered.
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
## 1 7 246596.7 499976.7 253426.6print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R
## 61 60 6.144158 4945.213 995048.6print("at day 60, proportion = recovered/initial number ~approx 99.5%") ## [1] "at day 60, proportion = recovered/initial number ~approx 99.5%"output[output$time == 60,"R"] / output[output$time == 0,"S"]## [1] 0.9950496print("attempt to figure out end of I by extending time, and day 145 looks like the closest to absolute recovery")## [1] "attempt to figure out end of I by extending time, and day 145 looks like the closest to absolute recovery"(as.data.frame(ode(y = initial_state_values,
times = seq(from = 0, to = 150, by = 1),
func = sir_model,
parms = parameters))) %>%
mutate(differs = abs(initial_state_values['S'] - R)) %>%
arrange(differs, time) %>%
head(1)## time S I R differs
## 1 145 2.543637e-07 1.008694 999999 0.008694556- Susceptible (S) decreases over time.
- Infected (I) increases until it peaks close to 500,000 by day 7, then decreases.
- Recovered (R) increases over time and is 99.5% recovered at day 60, and if we extend the timeline it goes closest to 100% recovered at day 145.
—
SIR Model with a Dynamic Force of Infection
The previous activity introduced the concept of the force of infection. However, it made the unrealistic assumption that the force of infection can be expressed as a constant rate. In reality \(\lambda\) changes over time, depending on the number of infected people in the population.
The diagram for a simple SIR model looks as follows, with \(\lambda\) defined as a function of the infection rate, \(\beta\) (beta), and the proportion of the population that is infectious, \(I/N\).
Question: Write out the differential equations for an SIR model with a dynamic force of infection. Which assumptions are you making in this model structure?
The differential equations for an SIR model with a dynamic force of infection are:
\[\begin{align} \frac{dS}{dt} & = -\beta \frac{I}{N} S \\ \\ \frac{dI}{dt} & = \beta \frac{I}{N} S - \gamma I \\ \\ \frac{dR}{dt} & = \gamma I \end{align}\]
Some assumptions inherent in this model structure are:
- a homogeneous population - everyone in the same compartment is subject to the same hazards
- a well-mixed population - all susceptible people have the same risk as getting infected, dependent on the number of infected people
- a closed population - there are no births or deaths, so the population size stays constant
Simulate this SIR model and extend it to capture a dynamic force of infection. Our assumptions for the initial conditions (\(S = 999999\), \(I = 1\) and \(R = 0\)), simulation period (60 days) and \(\gamma\) (0.1 days\(^{-1}\)) remain the same as in the previous activity, and the daily infection rate \(\beta\) equals 1.
Plot the curves of infection, susceptibility and recovery over time.
Question: After how many days does the epidemic peak? What is the peak prevalence?
- Peak infection right about day 18.
- Peak prevalence of approximately 67%.
nicesubtitle <- "SIR Model with a Dynamic Force of Infection"
print(nicesubtitle)## [1] "SIR Model with a Dynamic Force of Infection"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_state_values <- c(S = 1000000 - 1,
I = 1,
R = 0)
parameters <- c(gamma = .1,
beta = 1)
# TIMESTEPS:
times <- seq(from = 0, to = 60, by = 1)
# SIR MODEL FUNCTION
sir_model2 <- 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_model2,
parms = parameters))
# PLOT
melt(as.data.frame(output), 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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste("gamma of", parameters['gamma'],
" beta of", parameters['beta']),color = "Compartment",
subtitle= nicesubtitle)
# Q: After how many days does the epidemic peak? What is the peak prevalence?
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
## 1 18 100188.4 669741.4 230070.2print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R
## 61 60 51.14455 11864.48 988084.4Explore how your epidemic curve changes under different assumptions for \(\beta\) and \(\gamma\).
Question: How does the pattern of the epidemic change with different values for \(\beta\) or \(\gamma\)? For example, how do these different parameter combinations affect when the peak of the epidemic happens, how many people are infected at the peak, and when the epidemic ends?
- Holding beta constant, increasing gamma results in a flatter infection curve with a later peak.
- Holding gamma constant, increasing beta results in a steeper infection curve with an earlier peak.
nicesubtitle <- "SIR Model with a Dynamic Force of Infection v2"
print(nicesubtitle)## [1] "SIR Model with a Dynamic Force of Infection v2"print("initial state values and parameters")## [1] "initial state values and parameters"sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
# tell R to look for variable names within the
# state and parameters objects
N <- S+I+R
# New: defining lambda as a function of beta and I:
lambda <- beta * I/N
# The differential equations
dS <- -lambda * S # people move out of (-) the
# S compartment at a rate lambda
# (force of infection)
dI <- lambda * S - gamma * I # people move into (+) the I compartment
# from S at a rate lambda,
# and move out of (-) the I compartment
# at a rate gamma (recovery)
dR <- gamma * I # people move into (+) the R compartment
# from I at a rate gamma
# Return the number of people in the S, I and R compartments at each
# timestep
# (in the same order as the input state variables)
return(list(c(dS, dI, dR)))
})
}
run_sir_model <-
function(beta, gamma, S0 = 999999, I0 = 1, R0 = 1, duration = 60) {
initial_state_values <- c(S = S0, # nearly the whole population we are
# modelling is susceptible to infection
I = I0, # the epidemic starts with a single
# infected person
R = R0) # there is no prior immunity in the
# population
parameters <- c(beta = beta, # the rate of infection,
gamma = gamma) # the rate of recovery,
# which acts on those infected
times <- seq(from = 0, to = duration, by = 1)
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters))
print("Plotting the proportion of people in each compartment over time")
plot(x = output$time, # time on the x axis
y = output$S, # the number of susceptible people at
col = "blue", # each timestep on the y axis
type = "l", # type = "l" tells R we want lines
# rather than points
ylim = c(0,(S0+I0+R0)), # the limits of the y axis
# (from 0 to the total number of
# people)
xlab = "Time (days)",
ylab = "Number of people") # add axis labels
lines(x = output$time, # add the number of
y = output$I, # infected people at each
col = "red") # timestep on the y axis
lines(x = output$time, # number of recovered
y = output$R, # people at each
col = "green") # timestep on the y axis
legend(x = "right", # add a legend on the right-hand
# side of the plot
legend = c("S", "I", "R"), # labels S, I and R for blue,
col = c("blue", "red", "green"), # red, green lines respectively
lty = c(1,1)) # both lines are
# solid linetype (lty = 1)
title(main = paste("beta =", beta, # main title
", gamma =", gamma))
}def.par <- par(no.readonly = TRUE)
# Prepare the format for the plot display
# par(mfrow = c(2,4)) # Show 2*4 separate plots in one page
# (2 rows, 4 per row)
par(mfrow = c(2,2))
# Run the model and plot the output with different beta and gamma values.
# You might also have to change the time that you run the model for,
# by changing the value of the duration argument,
# in units of days (for example, to reproduce your model above,
# this would be duration = 60).
run_sir_model(beta = 1,
gamma = .1,
duration = 60)## [1] "Plotting the proportion of people in each compartment over time"run_sir_model(beta = 1,
gamma = .5,
duration = 60)## [1] "Plotting the proportion of people in each compartment over time"run_sir_model(beta = 3,
gamma = .1,
duration = 60) ## [1] "Plotting the proportion of people in each compartment over time"run_sir_model(beta = 3,
gamma = .5,
duration = 60)## [1] "Plotting the proportion of people in each compartment over time"
# Repeat the run_sir_model() command here for every parameter combination
# you want to try, and adapt the par() command above depending
# on how many plots you want to look at
par(def.par) #- reset to defaultDifferent recovery rates affect the epidemic just as much as different forces of infection. With \(\beta\) held constant at 1, an increasing value for \(\gamma\) tends to lead to a later and lower peak of infected people, and an earlier rise in the recovered curve. If people can stay infected for a long time before recovering (\(\gamma\) = 0.025, corresponding to an average duration of infection of 40 days), the number of infected people stays high over a longer period and declines slowly - the epidemic flattens out. In contrast, if recovery happens very quickly after infection (\(\gamma\) = 0.5), there is a only small peak in the prevalence of infection and the epidemic dies out quickly. If \(\gamma\) is as large as 1, no epidemic takes place after the introduction of 1 infected case.
—
SIR dynamics with varying parameters
The roles of \(\beta\) and \(\gamma\).
Question: Imagine a disease where every person infects 1 person on average, every 2 days, and is infectious for 4 days. What are the values for \(\beta\) and \(\gamma\)?
- beta of 0.5 and gamma of 0.25
So far, when looking at the dynamics of an epidemic, we have always plotted the number of people in each compartment. Sometimes however it is more informative to look at the proportion of the population in each compartment rather than the actual number, particularly if the population size changes over time. Although we are modelling an epidemic in a closed population with constant population size in this activity, it is a good opportunity to practice plotting the prevalence as a proportion of the total population who are in the \(S\), \(I\) and \(R\) compartments.
Model the epidemic with the \(\beta\) and \(\gamma\) values you calculated above, assuming introduction of a single infected person in a totally susceptible population of 1 million. Then plot the proportion of susceptible, infected and recovered people over time.
Question: What do you observe?
- Susceptible decreases then hovers around 20% after 120 days
- Infection peaks at 15.3% on day 55 then goes down
- Recovered increases then hovers around 80% after 120 days
nicesubtitle <- "SIR dynamics with varying parameters v1"
print(nicesubtitle)## [1] "SIR dynamics with varying parameters v1"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_state_values <- c(S = 1000000 - 1,
I = 1,
R = 0)
(parameters <- c(beta = 1/2,
gamma = 1/4))## beta gamma
## 0.50 0.25# TIMESTEPS:
times <- seq(from = 0, to = 120, by = 1)
# 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)
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
## 1 55 489187.4 153308.3 357504.3 48.919 15.331 35.75print("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
## 121 120 203205.6 26.39732 796768 20.321 0.003 79.677print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste("gamma of", parameters['gamma'],
" beta of", parameters['beta']),color = "Compartment",
subtitle= nicesubtitle)
Imagine an intervention is introduced to control infection, for example infected people are isolated so that they cannot spread infection. As a result, \(\beta\) drops to 0.1. Model the epidemic under this new scenario.
Question: What do you observe?
- beta of 0.1 and gamma of 0.25
- No epidemic
- Susceptible stays around 100% after 120 days
- Infection stays around 0% after 120 days
- Recovered stays around 0% after 120 days
nicesubtitle <- "SIR dynamics with varying parameters v2"
print(nicesubtitle)## [1] "SIR dynamics with varying parameters v2"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_state_values <- c(S = 1000000 - 1,
I = 1,
R = 0)
(parameters <- c(beta = .1,
gamma = 1/4))## beta gamma
## 0.10 0.25# TIMESTEPS:
times <- seq(from = 0, to = 700, by = 1)
# 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)
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
## 1 0 999999 1 0 100 0 0print("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
## 701 700 999998.3 9.971222e-44 1.666665 100 0 0print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste("gamma of", parameters['gamma'],
" beta of", parameters['beta']),color = "Compartment",
subtitle= nicesubtitle)
Now, vary \(\gamma\) with \(\beta\) held constant at 0.1 days\(^{-1}\). Notice that there are some values where the initial infected case fails to cause an epidemic.
Question: Assuming \(\beta\) equals 0.1 days\(^{-1}\), what value of \(\gamma\) do you need in order to get an epidemic? In real life, what could give rise to such a change in \(\gamma\)?
As the infection rate becomes lower, i.e. the infection spreads more slowly, you have to follow the population over a longer period to observe a visible change in the proportion infected. You should see that different combinations of infection and recovery rates can give rise to very different infection dynamics.
- beta of 0.1 and gamma of 0.01
- Recovery rate gamma might drop if the pathogen mutates, making it more difficult to recover
- Susceptible decreases then hovers around <1% after 218 days
- Infection peaks at 67% on day 180 then goes down
- Recovered increases then hovers around 99% after 618 days
nicesubtitle <- "SIR dynamics with varying parameters v3"
print(nicesubtitle)## [1] "SIR dynamics with varying parameters v3"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_state_values <- c(S = 1000000 - 1,
I = 1,
R = 0)
(parameters <- c(beta = .1,
gamma = .01))## beta gamma
## 0.10 0.01# TIMESTEPS:
times <- seq(from = 0, to = 700, by = 1)
# 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)
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
## 1 180 100188.8 669741.4 230069.8 10.019 66.974 23.007print("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
## 701 700 47.44886 4366.841 995585.7 0.005 0.437 99.559print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste("gamma of", parameters['gamma'],
" beta of", parameters['beta']),color = "Compartment",
subtitle= nicesubtitle)
Question: Can you think of a condition involving \(\beta\) and \(\gamma\) that is necessary for an epidemic?
- For an epidemic to occur, recovery rate gamma needs to be lower than infection rate beta; beta/gamma > 1.
—
Simulating the effective reproduction number Reff
In the simple SIR model, we get an epidemic only if \(\beta\)/\(\gamma\) > 1, i.e. if the average number of secondary cases caused by a single infected case in a totally susceptible population is greater than 1. As susceptibility declines over the course of the epidemic, the effective reproduction number Reff determines the shape of the epidemic curve as it reflects the amount of immunity in the population at any given time.
In a simple homogenous SIR model, Reff is directly related to the proportion of the population that is susceptible: \[\begin{align} R_{eff} = R_{0} \frac{S}{N} \end{align}\]
Model an epidemic and study the connection between the behaviour of Reff and the epidemic curve. We are looking at a closed fully susceptible population, into which a single infected person is introduced. Chose your combination of parameters so that R0 > 1. Solve the differential equations, and use the output to calculate Reff at each timestep. Generate 2 plots for comparison: one of the proportion of susceptible, infected and recovered individuals over time, and one of the effective reproduction number over time.
nicesubtitle <- "Simulating the effective reproduction number R eff v1"
print(nicesubtitle)## [1] "Simulating the effective reproduction number R eff v1"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_state_values <- c(S = 1000000-1,
I = 1,
R = 0)
duration <- 150
beta <- 1/2
gamma <- 1/4
R0 <- beta / gamma
(parameters <- c(beta = beta,
gamma = gamma,
R0 = R0))## beta gamma R0
## 0.50 0.25 2.00# TIMESTEPS:
times <- seq(from = 0, to = duration, by = 1)
# 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):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 55 489187.4 153308.3 357504.3 48.919 15.331 35.75 0.9783747print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(-Reff, time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 55 489187.4 153308.3 357504.3 48.919 15.331 35.75 0.9783747print("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
## 151 150 203187.7 0.3076275 796812 20.319 0 79.681 0.4063754print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0']),
color = "Compartment",
subtitle= nicesubtitle)
print("2 y-axis plot, where S and R eff overlap")## [1] "2 y-axis plot, where S and R eff overlap"output2 <- output %>%
mutate(Reff_ = Reff / parameters['R0']) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = ifelse(variable == "Reff_", value, value / sum(initial_state_values)))
ggplot(output2, aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line(linetype="dotdash", aes(y=proportion)) +
scale_y_continuous(sec.axis= sec_axis(~.* parameters['R0'], name="R effective")) +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","orange")) +
scale_shape_manual(values = c(0,4,1,6)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Question: How does Reff vary over the course of the epidemic? What do you notice about the connection between the change in Reff and the epidemic curve over time? In relation to Reff, when does the epidemic peak and start to decline?
- R eff varies proportional to Susceptible (S) over time.
- Peak infection is on Day 55 where (I) is 15%.
- For this run, the point when R eff goes below 1 is also on day 55, where R eff = 0.978. (Note that inflection points can occur throughout the smallest time period, and sometimes the data row with peak infection might not be the exact same row showing R eff crossing below 1.)
Use a daily infection rate of 0.4 and a daily recovery rate of 0.1 to get an R0 of 4. Model this epidemic over the course of about 3 months (100 days).
- With R0 increasing to 4 (from 2), peak infection and point when R eff goes below 1 occurs slightly earlier (days 50-51), and reaches a higher peak of about 40% (versus about 15%).
nicesubtitle <- "Simulating the effective reproduction number R eff v2"
print(nicesubtitle)## [1] "Simulating the effective reproduction number R eff v2"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_state_values <- c(S = 1000000-1,
I = 1,
R = 0)
duration <- 100
beta <- 0.4
gamma <- .1
R0 <- beta / gamma
(parameters <- c(beta = beta,
gamma = gamma,
R0 = R0))## beta gamma R0
## 0.4 0.1 4.0# TIMESTEPS:
times <- seq(from = 0, to = duration, by = 1)
# 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):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 50 258083.8 403299 338617.2 25.808 40.33 33.862 1.032335print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(-Reff, time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 51 219671.7 401423.5 378904.8 21.967 40.142 37.89 0.8786866print("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 20473.64 7372.583 972153.8 2.047 0.737 97.215 0.08189457print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0']),
color = "Compartment",
subtitle= nicesubtitle)
print("2 y-axis plot, where S and R eff overlap")## [1] "2 y-axis plot, where S and R eff overlap"output2 <- output %>%
mutate(Reff_ = Reff / parameters['R0']) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = ifelse(variable == "Reff_", value, value / sum(initial_state_values)))
ggplot(output2, aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line(linetype="dotdash", aes(y=proportion)) +
scale_y_continuous(sec.axis= sec_axis(~.* parameters['R0'], name="R effective")) +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","orange")) +
scale_shape_manual(values = c(0,4,1,6)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
—
Modelling population turnover
We want to incorporate births and deaths into the basic SIR model. The structure for this model looks like this:
The differential equations for an SIR model with population turnover look like this:
\[\begin{align} \frac{dS}{dt} & = -\beta \frac{I}{N} S - \mu S + bN \\ \frac{dI}{dt} & = \beta \frac{I}{N} S - \gamma I - \mu I \\ \frac{dR}{dt} & = \gamma I - \mu R \end{align}\]
This structure assumes that every individual in each compartment experiences the same background mortality \(\mu\) (there is no additional mortality from the infection for example, and we make no distinction by age). Those who have died no longer contribute to infection (a sensible assumption for many diseases). Babies are all born at a rate b into the susceptible compartment.
Note that, as always, we calculate the people dying in each of the compartments by multiplying the number in that compartment by the rate \(\mu\). Even though babies are all born into the same compartment, the birth rate still depends on the population in all of the compartments, hence why we need to multiply b by the total population size N. In the following activity, we choose a value of b to allow a constant population size, so that all deaths are balanced by births.
Modelling an acute disease epidemic in a fully susceptible human population
Parameters:
\(\beta\) = 0.4 days\(^{-1}\) = 0.4 * 365 years\(^{-1}\)
\(\gamma\) = 0.2 days\(^{-1}\) = 0.2 * 365 years\(^{-1}\)
\(\mu\) = 1/70 years\(^{-1}\)
b = 1/70 years\(^{-1}\)
Extend the SIR model to capture this population turnover process. In the first instance, we are looking at an acute disease epidemic introduced into a fully susceptible human population. The infection and recovery rates are 0.4 and 0.2 days\(^{-1}\) respectively. We can calculate the background mortality rate based on the average human lifespan - let’s assume that this is 70 years, so mu = 1/70 years\(^{-1}\), or 1/(70x365) days\(^{-1}\). In this example, we are also assuming that the population size stays constant over time at 10\(^{6}\) - to achieve that, we set the birth rate b = mu.
Code this model and run it first for 1 year, and plot only the number of infected people over time at this point. Confirm that you observe your regular SIR epidemic curve before proceeding.
nicesubtitle <- "Modelling population turnover v1, 365 days"
print(nicesubtitle)## [1] "Modelling population turnover v1, 365 days"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_state_values <- c(S = 1000000-1,
I = 1,
R = 0)
duration <- 365
tsteps <- 1
beta <- 0.4
gamma <- 0.2
mu <- 1/(70 * 365)
birth <- mu
R0 <- beta / gamma
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
mu = mu,
birth = birth
))## beta gamma R0 mu birth
## 4.000000e-01 2.000000e-01 2.000000e+00 3.913894e-05 3.913894e-05# 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) -(S * mu) +(N * birth)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 68 512728.7 153245.8 334025.5 51.273 15.325 33.403 1.025457print("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
## 1 69 482248.6 153076.3 364675.1 48.225 15.308 36.468 0.9644972print("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
## 366 365 212026.6 3.463514e-10 787973.4 21.203 0 78.797 0.4240533print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle)
Now, change the duration of the model run to several generations, for example 400 years.
A note on timesteps:
So far, we have only modelled short-term disease dynamics, where it was intuitive to express the timesteps in units of days. In this example, the simulation timescale is much longer, so it might be easier to change the units to years (although the result is also correct if we model 400 years as 146000 days). Now, we need to start thinking about the interval of the timesteps, in the code below.
nicesubtitle <- "Modelling population turnover v2, \n 400 years, daily intervals step2"
print(nicesubtitle)## [1] "Modelling population turnover v2, \n 400 years, daily intervals step2"print("initial state values and parameters")## [1] "initial state values and parameters"# times <- seq(from = 0, to = 400, by = 1/365) # from 0 to 400 YEARS in DAILY intervals
# MODEL INPUTS:
initial_state_values <- c(S = 1000000-1,
I = 1,
R = 0)
duration <- 365 * 400
tsteps <- 2
beta <- 0.4
gamma <- 0.2
mu <- 1/(70 * 365)
birth <- mu
R0 <- beta / gamma
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
mu = mu,
birth = birth
))## beta gamma R0 mu birth
## 4.000000e-01 2.000000e-01 2.000000e+00 3.913894e-05 3.913894e-05# 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) -(S * mu) +(N * birth)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 68 512728.3 153246 334025.7 51.273 15.325 33.403 1.025457print("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
## 1 70 453778.1 151113.9 395108 45.378 15.111 39.511 0.9075561print("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
## 73001 146000 503432.5 69.27998 496498.2 50.343 0.007 49.65 1.006865print("Plot Infection only - 400 days: ")## [1] "Plot Infection only - 400 days: "output %>%
filter(time <= 400) %>%
select(time, I) %>%
ggplot(aes(x=time, y=I)) +
geom_line(linetype="dotdash", color="red") +
geom_jitter(color="red", shape=4, show.legend = FALSE) +
xlab("Time (days)") +
ylab("I Infected") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle)
print("Plot Infection only - all days: ")## [1] "Plot Infection only - all days: "output %>%
select(time, I) %>%
ggplot(aes(x=time, y=I)) +
geom_line(linetype="dotdash", color="red") +
geom_jitter(color="red", shape=4, show.legend = FALSE) +
xlab("Time (days)") +
ylab("I Infected") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle)
print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle)
Before investigating the disease dynamics, play around with the timesteps at which you are solving the model in your code above and plot the output. In the cell below, it is daily - also try solving it only at every 2, 3, 4 and 5 days. Keep in mind that as we are running the model for longer now, it might take your computer longer to produce the output.
Question: How does changing the interval of the timesteps to solve the model at influence the output? Does the plot look correct in each case? If not, at what resolution of timesteps do you get erroneous results, and why?
- Extending duration Time to much longer periods reveals stairstep patterns in visualizations, owing to the compressed space to plot but which may follow natural up and down patterns under limited space.
- Larger intervals of timesteps might cause the routine not to resolve properly. At the margins of large timesteps, the visualizations reveal stairstep patterns instead of smooth curves.
nicesubtitle <- "Modelling population turnover v3, \n 400 years, yearly intervals"
print(nicesubtitle)## [1] "Modelling population turnover v3, \n 400 years, yearly intervals"print("initial state values and parameters")## [1] "initial state values and parameters"# times <- seq(from = 0, to = 400, by = 1/365) # from 0 to 400 YEARS in DAILY intervals
# MODEL INPUTS:
initial_state_values <- c(S = 1000000-1,
I = 1,
R = 0)
duration <- 400
tsteps <- 1/365
beta <- 0.4 * 365
gamma <- 0.2 * 365
mu <- 1/(70 * 1)
birth <- mu
R0 <- beta / gamma
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
mu = mu,
birth = birth
))## beta gamma R0 mu birth
## 146.00000000 73.00000000 2.00000000 0.01428571 0.01428571# 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) -(S * mu) +(N * birth)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 0.1863014 512728.7 153245.8 334025.5 51.273 15.325 33.403 1.025457print("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
## 1 0.1890411 482248.6 153076.3 364675.1 48.225 15.308 36.468 0.9644972print("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
## 146001 400 503898.4 80.8277 496020.8 50.39 0.008 49.602 1.007797print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (years)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (years)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle)
Finally, run and plot the number of infected people again with the timestep that seems most sensible based on your findings, to study the long-term disease dynamics.
Question: What do you observe about the long-term disease dynamics under these assumptions? Can you explain why this pattern occurs based on what you have learnt in the last weeks?
- The infection rate cycles over time, with the peaks getting progressively smaller.
- Infections require Susceptibles, which can get depleted and need to be replenished through births .
Now make a plot that also shows the number of susceptible and recovered individuals, as a proportion of the total population.
Question: Does this confirm your hypothesis? Also think about how it relates to the effective reproduction number.
- Yes, the Infections go up as Susceptibles increase, though the peaks for both S and I seem to also decrease propotionally over time.
- Reff tracks with the Susceptibles compartment.
# Redo using daily interval and plot n=400 then all
nicesubtitle <- "Modelling population turnover v4, redo"
print(nicesubtitle)## [1] "Modelling population turnover v4, redo"print("initial state values and parameters")## [1] "initial state values and parameters"nicesubtitle <- "Modelling population turnover"
# MODEL INPUTS:
initial_state_values <- c(S = 1000000-1,
I = 1,
R = 0)
duration <- 365 * 400
tsteps <- 2
beta <- 0.4
gamma <- 0.2
mu <- 1/(70 * 365)
birth <- mu
R0 <- beta / gamma
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
mu = mu,
birth = birth
))## beta gamma R0 mu birth
## 4.000000e-01 2.000000e-01 2.000000e+00 3.913894e-05 3.913894e-05# 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) -(S * mu) +(N * birth)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 68 512728.3 153246 334025.7 51.273 15.325 33.403 1.025457# point when R eff goes below 1 - day / where R eff = ?
output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 70 453778.1 151113.9 395108 45.378 15.111 39.511 0.9075561print("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
## 73001 146000 503432.5 69.27998 496498.2 50.343 0.007 49.65 1.006865#
# Limit to 400 days to see short term details
#
nicesubtitle <- "Modelling population turnover v4, redo n=400"
print(nicesubtitle)## [1] "Modelling population turnover v4, redo n=400"print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"output %>%
filter(time <= 400) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
# Plotting the proportion of people in each compartment over time
output %>%
filter(time <= 400) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
filter(time <= 400) %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle)
# 2 y-axis plot, where S and R eff overlap
output2 <- output %>%
filter(time <= 400) %>%
mutate(Reff_ = Reff / parameters['R0']) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = ifelse(variable == "Reff_", value, value / sum(initial_state_values)))
ggplot(output2, aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line(linetype="dotdash", aes(y=proportion)) +
scale_y_continuous(sec.axis= sec_axis(~.* parameters['R0'], name="R effective")) +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","orange")) +
scale_shape_manual(values = c(0,4,1,6)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
#
# Open up to long time period
#
nicesubtitle <- "Modelling population turnover v4, redo n=ALL"
print(nicesubtitle)## [1] "Modelling population turnover v4, redo n=ALL"print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
# Plotting the proportion of people in each compartment over time
output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle)
# 2 y-axis plot, where S and R eff overlap
output2 <- output %>%
mutate(Reff_ = Reff / parameters['R0']) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = ifelse(variable == "Reff_", value, value / sum(initial_state_values)))
ggplot(output2, aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line(linetype="dotdash", aes(y=proportion)) +
scale_y_continuous(sec.axis= sec_axis(~.* parameters['R0'], name="R effective")) +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","orange")) +
scale_shape_manual(values = c(0,4,1,6)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
In the second example, we are modelling a similar acute disease, but this time we are changing the mortality and birth rates to represent a population with a much more rapid turnover. The infection parameters are the same as before, but we are assuming the average lifespan is 4 weeks this time.
Code this scenario below and run the model for 1 year, and plot the prevalence of susceptibility, infection and recovery over time.
Question: How do the disease dynamics compare to the previous example? Why does this occur (what is different compared to the disease in the human population)?
- With the short life span scenario, we do not observe the same epidemic cycle pattern.
- Also, Infections do not go to zero over time.
# short life people scenario
nicesubtitle <- "Modelling population turnover v5, short life span"
print(nicesubtitle)## [1] "Modelling population turnover v5, short life span"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
initial_state_values <- c(S = 1000000-1,
I = 1,
R = 0)
duration <- 365
tsteps <- 1
beta <- 0.4
gamma <- 0.2
mu <- 1/(7 * 4)
birth <- mu
R0 <- beta / gamma
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
mu = mu,
birth = birth
))## beta gamma R0 mu birth
## 0.40000000 0.20000000 2.00000000 0.03571429 0.03571429# 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) -(S * mu) +(N * birth)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 83 591760.7 131429.3 276810 59.176 13.143 27.681 1.183521print("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
## 1 93 498195.6 105147.1 396657.3 49.82 10.515 39.666 0.9963913print("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
## 366 365 589257.4 62237.4 348505.2 58.926 6.224 34.851 1.178515print("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) %>%
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")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
# Plotting the proportion of people in each compartment over time
output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle)
# 2 y-axis plot, where S and R eff overlap
output2 <- output %>%
mutate(Reff_ = Reff / parameters['R0']) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = ifelse(variable == "Reff_", value, value / sum(initial_state_values)))
ggplot(output2, aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line(linetype="dotdash", aes(y=proportion)) +
scale_y_continuous(sec.axis= sec_axis(~.* parameters['R0'], name="R effective")) +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","orange")) +
scale_shape_manual(values = c(0,4,1,6)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Additional questions
The examples we have modelled here can be compared to real-world diseases. An acute disease that causes a pattern of epidemic cycles in human populations like in the first example is measles. The second example on the other hand is more comparable to endemic swine flu on pig farms, where pigs enter and leave the farm in a matter of weeks (hence the ‘average lifespan’ that you modelled, of 4 weeks).
However, the basic SIR model is very simple and does not account for many of the factors that can affect disease dynamics in reality. For example, measles notifications data from the pre-vaccination era suggest that epidemic oscillations did not flatten out over time like in our model. This tells us that there are other processes that sustain these cycles, that we have not included in the model.
Question: Can you think of other drivers for epidemic cycles?
- mutation
- lowered immunity level
- lack of public health standards or consistency in treatment
- cycles of nature coinciding with dormancy/activity of disease
This activity was also your first look at a simple demographic process with births and deaths. As you might have noticed, the population size in this example stayed constant over time. However, if we were to model the demography of a country over time, an increasing population size is usually more realistic.
Question: How would you change the birth and death rate to model a growing population? What other factors could drive population growth?
- a growing population can be caused by a increasing birth rate; can be caused by improvements in public health measures, technology/medicine, incomes/social services, etc.
- a growing population can be caused by a decreasing death rate; can be caused by improvements in public health measures, technology/medicine, incomes/social services, etc.
- increased immigration and decreased outflows can also result in a growing population
In this activity, we also made the assumption that all babies are born susceptible. Of course, this is not always the case in reality: maternal antibodies can confer immunity to newborns, neonatal vaccination can make newborns immune to infection, and some infections can be transmitted from mothers to their children, like HIV.
Question: How would you adapt the differential equations to model a disease where a proportion p of babies born to infected mothers are infected at birth? Write them out on paper.
- the current function calculates the current susceptible (S) to include the (b) birth term multiplied by S; we do not see a separate term coming out of the infected (I) so that may be an option.
- There is also the thought that b and mu might already factor in these scenarios.
—
A simple model for vaccination
We are going back to the simple SIR model, without births or deaths, to look at the effect of vaccination in a very simple way – we are assuming it already happened before we run our model! By changing the initial conditions, we can prepare the population so that it has received a certain coverage of vaccination. Even though this is very simplistic, it will allow us to study some important effects of vaccination on the infection dynamics.
We are starting with the transmission and recovery parameters beta = 0.4 days\(^{-1}\) and gamma = 0.1 days \(^{-1}\). To incorporate immunity from vaccination in the model, we assume that a proportion p of the total population starts in the recovered compartment, representing the vaccine coverage and assuming the vaccine is perfectly effective. Again, we assume the epidemic starts with a single infected case introduced into the population.
Model this scenario for a duration of 2 years, assuming that the vaccine coverage is 50%, and plot the prevalence in each compartment over time. Confirm that you observe an epidemic peaking at around 125 days after introduction of the infectious person in the population. Also have a look at the proportion susceptible and recovered, to double-check this looks like what you would expect given your initial conditions.
Modelling a disease where \(\beta\) = 0.4 days\(^{-1}\), \(\gamma\) = 0.1 days \(^{-1}\) and the vaccine coverage p = 0.5
Question: Does everyone in the population need to be vaccinated in order to prevent an epidemic? What do you observe if you model the infection dynamics with different values for p? Can you explain why?
- No, it is not the case every time that everyone needs to be vaccinated to prevent an outbreak. As seen in the second run for this scenario, using the critical vaccination threshold of 75% (instead of the first run value of 50%) is enough to prevent an outbreak.
nicesubtitle <- "A simple model for vaccination v1, R0 of 4 and 50% vacc"
print(nicesubtitle)## [1] "A simple model for vaccination v1, R0 of 4 and 50% vacc"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000
duration <- 365 * 2
tsteps <- 1
beta <- 0.4
gamma <- 0.1
pvacc <- .50
R0 <- beta / gamma
# Reff threshold calc
pvacc_thresh <- 1 - (1/R0)
initial_state_values <- c(S = (N-1) * (1 - pvacc),
I = 1,
R = (N-1) * pvacc )
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## beta gamma R0 pvacc pvacc_thresh
## 0.40 0.10 4.00 0.50 0.75# 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):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 130 248893.4 76711.5 674395.1 24.889 7.671 67.44 0.9955737print("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
## 1 130 248893.4 76711.5 674395.1 24.889 7.671 67.44 0.9955737print("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
## 731 730 101593.4 8.665453e-11 898406.6 10.159 0 89.841 0.4063737print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
Hopefully, you have seen that once a certain proportion of the population is immune, no epidemic occurs. This is called the critical vaccination threshold or herd immunity threshold.
nicesubtitle <- "A simple model for vaccination v2, R0 of 4 and calc vacc threshold"
print(nicesubtitle)## [1] "A simple model for vaccination v2, R0 of 4 and calc vacc threshold"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000
duration <- 365 * 2
tsteps <- 1
beta <- 0.4
gamma <- 0.1
pvacc <- .50
R0 <- beta / gamma
# Reff threshold calc
pvacc_thresh <- 1 - (1/R0)
# override pvacc with pvacc_thresh
pvacc <- pvacc_thresh
initial_state_values <- c(S = (N-1) * (1 - pvacc),
I = 1,
R = (N-1) * pvacc )
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## beta gamma R0 pvacc pvacc_thresh
## 0.40 0.10 4.00 0.75 0.75# 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):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 0 249999.8 1 749999.2 25 0 75 0.999999print("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
## 1 0 249999.8 1 749999.2 25 0 75 0.999999print("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
## 731 730 249927 0.9893463 750072 24.993 0 75.007 0.9997081print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
- When the p vacc rate was 50%, Infection peaked at approx. 7.7% on Day 130.
- When the p vacc rate was increased to the calculated threshold of 1-(1/R0) = 75%, Infection never had a chance to increase.
- Some trial and error runs show that lowering the p vacc rate lead to a quicker ramp up in Infection, as well as a higher peak.
Investigate how the herd immunity threshold changes if we are modelling a disease with a different infection and recovery rate.
Question: What proportion of the population needs to be vaccinated in order to prevent an epidemic if beta = 0.4 and gamma = 0.2 days\(^{-1}\)? What if beta = 0.6 and gamma = 0.1 days\(^{-1}\)?
if beta = 0.4 and gamma = 0.2 days\(^{-1}\), trying 10% vaccination rate (v3), then critical vaccination threshold of 83% afterwards.
No, it is not the case every time that everyone needs to be vaccinated to prevent an outbreak. In the prior scenario, critical vaccination threshold of .75 was enough. For the following example, the cvt is about .83.
nicesubtitle <- "A simple model for vaccination v3, lower R0 and 10% vacc"
print(nicesubtitle)## [1] "A simple model for vaccination v3, lower R0 and 10% vacc"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000
duration <- 365 * 2
tsteps <- 1
beta <- 0.4
gamma <- 0.2
pvacc <- .10 # try 10 percent
R0 <- beta / gamma
# Reff threshold calc
pvacc_thresh <- 1 - (1/R0)
initial_state_values <- c(S = (N-1) * (1 - pvacc),
I = 1,
R = (N-1) * pvacc )
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## beta gamma R0 pvacc pvacc_thresh
## 0.4 0.2 2.0 0.1 0.5# 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):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 83 496363.7 106094 397542.3 49.636 10.609 39.754 0.9927275print("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
## 1 83 496363.7 106094 397542.3 49.636 10.609 39.754 0.9927275print("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
## 731 730 240812.5 1.350074e-23 759187.5 24.081 0 75.919 0.481625print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
Use critical vaccination threshold:
nicesubtitle <- "A simple model for vaccination v4, lower R0 and calc vacc threshold"
print(nicesubtitle)## [1] "A simple model for vaccination v4, lower R0 and calc vacc threshold"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000
duration <- 365 * 2
tsteps <- 1
beta <- 0.4
gamma <- 0.2
pvacc <- .0 # to be overwritten below with calculated threshold
R0 <- beta / gamma
# Reff threshold calc
pvacc_thresh <- 1 - (1/R0)
# override pvacc with pvacc_thresh
pvacc <- pvacc_thresh
initial_state_values <- c(S = (N-1) * (1 - pvacc),
I = 1,
R = (N-1) * pvacc )
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## beta gamma R0 pvacc pvacc_thresh
## 0.4 0.2 2.0 0.5 0.5# 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):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 0 499999.5 1 499999.5 50 0 50 0.999999print("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
## 1 0 499999.5 1 499999.5 50 0 50 0.999999print("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
## 731 730 499854.6 0.9788434 500144.5 49.985 0 50.014 0.9997091print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
Changed parameters using 10% vaccination rate –>
- if beta = 0.6 and gamma = 0.1 days\(^{-1}\):
nicesubtitle <- "A simple model for vaccination v5, higher R0 and 10% vacc"
print(nicesubtitle)## [1] "A simple model for vaccination v5, higher R0 and 10% vacc"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000
duration <- 365 * 2
tsteps <- 1
beta <- 0.6
gamma <- 0.1
pvacc <- .10 # try 10 percent
R0 <- beta / gamma
# Reff threshold calc
pvacc_thresh <- 1 - (1/R0)
initial_state_values <- c(S = (N-1) * (1 - pvacc),
I = 1,
R = (N-1) * pvacc )
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## beta gamma R0 pvacc pvacc_thresh
## 0.6000000 0.1000000 6.0000000 0.1000000 0.8333333# 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):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 35 163728.4 452241 384030.6 16.373 45.224 38.403 0.9823704print("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
## 1 35 163728.4 452241 384030.6 16.373 45.224 38.403 0.9823704print("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
## 731 730 4167.877 1.352032e-23 995832.1 0.417 0 99.583
## Reff
## 731 0.02500726print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
Changed parameters using critical vaccination threshold –>
- if beta = 0.6 and gamma = 0.1 days\(^{-1}\):
nicesubtitle <- "A simple model for vaccination v6, higher R0 and calc vacc threshold"
print(nicesubtitle)## [1] "A simple model for vaccination v6, higher R0 and calc vacc threshold"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000
duration <- 365 * 2
tsteps <- 1
beta <- 0.6
gamma <- 0.1
pvacc <- .0 # to be overwritten below with calculated threshold
R0 <- beta / gamma
# Reff threshold calc
pvacc_thresh <- 1 - (1/R0)
# override pvacc with pvacc_thresh
pvacc <- pvacc_thresh
initial_state_values <- c(S = (N-1) * (1 - pvacc),
I = 1,
R = (N-1) * pvacc )
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## beta gamma R0 pvacc pvacc_thresh
## 0.6000000 0.1000000 6.0000000 0.8333333 0.8333333# 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):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 0 166666.5 1 833332.5 16.667 0 83.333 0.999999print("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
## 1 0 166666.5 1 833332.5 16.667 0 83.333 0.999999print("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
## 731 730 166593.9 0.9841127 833405.1 16.659 0 83.341 0.9995634print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n p vacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
As you can see, the proportion of the population that needs to be vaccinated varies with different infection-related parameters. Think about why that is so in the context of herd immunity, and what is different between the scenarios you have just modelled.
Question: Remember that vaccination changes the effective reproduction number, by reducing the number of people who are susceptible. Based on your answers to the previous questions, can you use the formula for the effective reproduction number Reff to derive a formula for calculating the critical vaccination threshold?
“In mathematical modelling terms, herd immunity is just the same as saying that Reff < 1. We can derive the herd immunity threshold by solving the formula for Reff for p when Reff = 1:”
\[\begin{align} R_{eff} & = R_{0} \frac{S}{N} \\ R_{eff} & = R_{0} (1-p) \\ p = 1-\frac{1}{R_{0}} \end{align}\]
“Remember, we can calculate R0 by dividing \(\beta\) by \(\gamma\).”
# this is where we get the following formula
pvacc_threshold <- 1 - (1/R0)—
Neonatal vaccination to reduce prevalence of an endemic disease in livestock
A neonatal vaccine has recently become available against an endemic viral infection in livestock. The Ministry of Agriculture has asked you to model the impact of different vaccination scenarios to inform the design of a farm animal vaccination programme. They provide you with the following information:
- The infection rate of the disease in this population is 1 day\(^{-1}\) and the average duration of the infectious period is 20 days.
- The average lifespan of the population is 3 years.
- The total population size in the country is 300000 and has not varied over time. The disease is endemic but thought to be relatively rare in this population.
- The vaccine, which has perfect efficacy, is given to a proportion p_vacc of newborn animals (that is, it is a neonatal vaccine). It is currently not known whether immunity to the disease and to the vaccine is life-long or wanes over time.
Setting up the model:
The SIR structure needs to be extended to incorporate vaccinated births (going into the R compartment), unvaccinated births (going into the S compartment), deaths, and waning immunity. As we are modelling an endemic infection, the initial conditions for the population don’t matter as long as they add up to 300000 according to the instructions (this is because all initial conditions will end up at the same endemic equilibrium, given enough time). For the baseline scenario, we are assuming no vaccine coverage (p_vacc = 0) and no waning of immunity (\(\sigma\) = 0).
nicesubtitle <- "Neonatal livestock vaccination v2 run for baseline "
print(nicesubtitle)## [1] "Neonatal livestock vaccination v2 run for baseline "print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
# modelling endemic infection - initial conditions for population won't matter
endemic_I <- round(N * .1) # assume 1% "small" endemic infected
endemic_R <- round(N * .6) # assume 39% endemic recovered
duration <- 365 * 6 # 3 year life span so track for 6 years
tsteps <- 2 # chunk into 2 day intervals
beta <- 1 # rate given as 1 day^-1
gamma <- 1/20 # 20 days
pvacc <- 0 # assume no vaccine coverage
sigma <- 0 # no waning immunity
mu <- 1/(3 * 365) # 3 year life span
birth <- mu # stable population so set to same as mu
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
initial_state_values <- c(S = (N-(endemic_I + endemic_R)),
I = endemic_I,
R = endemic_R)
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0, # basic reproduction number
pvacc_thresh = pvacc_thresh, # vaccination threshold
mu = mu, # background mortality rate
birth = birth, # birth rate
pvacc = pvacc, # vaccine rate
sigma = sigma # immunity waning rate
))## beta gamma R0 pvacc_thresh mu birth
## 1.00000e+00 5.00000e-02 2.00000e+01 9.50000e-01 9.13242e-04 9.13242e-04
## pvacc sigma
## 0.00000e+00 0.00000e+00# 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) -(S * mu) +(N * birth * (1-pvacc)) +(sigma * R)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu) +(N * birth * (pvacc)) -(sigma * R)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 10 11876.03 78307.91 209816.1 3.959 26.103 69.939 0.7917352print("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
## 1 10 11876.03 78307.91 209816.1 3.959 26.103 69.939 0.7917352print("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
## 1096 2190 15273.97 5107.193 279618.8 5.091 1.702 93.206 1.018265print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
What is the endemic prevalence of the disease currently (the baseline prevalence), assuming permanent immunity?
- The baseline Infected rate appears to be 1.7%.
# uses values from output v2 run above
print("See equilibrium values assumed at 2 years: ")## [1] "See equilibrium values assumed at 2 years: "output[output$time == 365*2,]## time S I R still_Su preval_Inf propor_Re Reff
## 366 730 15250.79 5122.388 279626.8 5.084 1.707 93.209 1.016719print("See equilibrium values assumed at 3 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 3 years (+1 because step is 2): "output[output$time == 365*3+1,]## time S I R still_Su preval_Inf propor_Re Reff
## 549 1096 15275.12 5107.069 279617.8 5.092 1.702 93.206 1.018341# save SIR values at 3 years for use as initial values in next run
endem_S <- round(output[output$time == 365*3+1,2])
endem_I <- round(output[output$time == 365*3+1,3])
endem_R <- round(output[output$time == 365*3+1,4])
endem_N <- endem_S + endem_I + endem_RWhat proportion of newborn animals would you need to vaccinate to reduce the prevalence by half, assuming life-long immunity?
- Setting the vaccination rate at 0.50 reduces 3 or 5 year prevalence to 0.805% which is close to the target of .85%.
- Setting the vaccination rate at 0.48 reduces 5 year prevalence to 0.841% which is even closer to the target of .85%.
nicesubtitle <- "Neonatal livestock vaccination v3 reduce I by 50%"
print(nicesubtitle)## [1] "Neonatal livestock vaccination v3 reduce I by 50%"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
# modelling endemic infection - initial conditions for population won't matter
endemic_I <- endem_I # take from last run baseline at year 3
endemic_R <- endem_R # take from last run baseline at year 3
duration <- 365 * 6 # 3 year life span so track for 6 years
tsteps <- 2 # chunk into 2 day intervals
beta <- 1 # rate given as 1 day^-1
gamma <- 1/20 # 20 days
# trying to reduce prevalence from 1.7% down to 0.85% (half) via vaccination
pvacc <- 0.48 # 50% gets us 0.805%, so let us try 48% vaccination rate
sigma <- 0 # no waning immunity
mu <- 1/(3 * 365) # 3 year life span
birth <- mu # stable population so set to same as mu
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
initial_state_values <- c(S = (N-(endemic_I + endemic_R)),
I = endemic_I,
R = endemic_R)
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0,
pvacc_thresh = pvacc_thresh, # vaccination threshold
mu = mu, # background mortality rate
birth = birth, # birth rate
pvacc = pvacc, # vaccine rate
sigma = sigma # immunity waning rate
))## beta gamma R0 pvacc_thresh mu birth
## 1.00000e+00 5.00000e-02 2.00000e+01 9.50000e-01 9.13242e-04 9.13242e-04
## pvacc sigma
## 4.80000e-01 0.00000e+00# 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) -(S * mu) +(N * birth * (1-pvacc)) +(sigma * R)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu) +(N * birth * (pvacc)) -(sigma * R)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 0 15275 5107 279618 5.092 1.702 93.206 1.018333print("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
## 1 4 14768.5 5089.632 280141.9 4.923 1.697 93.381 0.9845667print("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
## 1096 2190 15273.94 2524.305 282201.8 5.091 0.841 94.067 1.018263print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
print("See equilibrium values assumed at 3 years (+1 because step is 2): ") ## [1] "See equilibrium values assumed at 3 years (+1 because step is 2): "output[output$time == 365*3+1,]## time S I R still_Su preval_Inf propor_Re Reff
## 549 1096 15276.8 2512.547 282210.7 5.092 0.838 94.07 1.018453print("See equilibrium values assumed at 5 years (+1 because step is 2): ") ## [1] "See equilibrium values assumed at 5 years (+1 because step is 2): "output[output$time == 365*5+1,]## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15274.72 2524.406 282200.9 5.092 0.841 94.067 1.018315Would it be possible to eliminate the disease from the population using neonatal vaccination under the assumption of lifelong immunity?
- Apparently, yes. We can use the “critical vaccination threshold” calculated based on R0 which is 1-(1/R0) == .95.
nicesubtitle <- "Neonatal livestock vaccination v4 48% vacc, no wane"
print(nicesubtitle)## [1] "Neonatal livestock vaccination v4 48% vacc, no wane"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
# modelling endemic infection - initial conditions for population won't matter
endemic_I <- endem_I # take from last run baseline at year 3
endemic_R <- endem_R # take from last run baseline at year 3
duration <- 365 * 6 # 3 year life span so track for 6 years
tsteps <- 2 # chunk into 2 day intervals
beta <- 1 # rate given as 1 day^-1
gamma <- 1/20 # 20 days
# trying to reduce prevalence from 1.7% down to 0.85% (half) via vaccination
pvacc <- 0.48 # placeholder, see below for threshold overwrite
sigma <- 0 # no waning immunity
mu <- 1/(3 * 365) # 3 year life span
birth <- mu # stable population so set to same as mu
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
pvacc <- pvacc_thresh # overwrite in attempt to eliminate
initial_state_values <- c(S = (N-(endemic_I + endemic_R)),
I = endemic_I,
R = endemic_R)
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0,
pvacc_thresh = pvacc_thresh, # vaccination threshold
mu = mu, # background mortality rate
birth = birth, # birth rate
pvacc = pvacc, # vaccine rate
sigma = sigma # immunity waning rate
))## beta gamma R0 pvacc_thresh mu birth
## 1.00000e+00 5.00000e-02 2.00000e+01 9.50000e-01 9.13242e-04 9.13242e-04
## pvacc sigma
## 9.50000e-01 0.00000e+00# 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) -(S * mu) +(N * birth * (1-pvacc)) +(sigma * R)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu) +(N * birth * (pvacc)) -(sigma * R)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 0 15275 5107 279618 5.092 1.702 93.206 1.018333print("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
## 1 2 14763.96 5098.295 280137.7 4.921 1.699 93.379 0.9842637print("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
## 1096 2190 13699.03 1.413955e-10 286301 4.566 0 95.434 0.9132689print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
print("See equilibrium values assumed at 3 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 3 years (+1 because step is 2): "(out_v4_y3 <- output[output$time == 365*3+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 549 1096 11466.83 1.326657e-06 288533.2 3.822 0 96.178 0.7644556print("See equilibrium values assumed at 5 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 5 years (+1 because step is 2): "(out_v4_y5 <- output[output$time == 365*5+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 13186.01 1.282425e-09 286814 4.395 0 95.605 0.8790674If the average duration of immunity is only 1 year, how would this impact the proportional reduction in the prevalence with the vaccine coverage you obtained above compared to the baseline?
- Without waning immunity, a vaccination rate of .95 is enough to handle the scenario with R0 of 20 and eliminate infections.
- With 1-year waning immunity, no vaccination results in a baseline prevalence of 6.5%.
- With 1-year waning immunity, a vaccination rate of .95 results in a baseline prevalence of 4.8%.
- The reduction in prevalence between no vaccination and a vaccination rate of .95 With 1-year waning immunity is 25%.
nicesubtitle <- "Neonatal livestock vaccination v5 no vacc, 1 yr wane"
print(nicesubtitle)## [1] "Neonatal livestock vaccination v5 no vacc, 1 yr wane"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
# modelling endemic infection - initial conditions for population won't matter
endemic_I <- endem_I # take from last run baseline at year 3
endemic_R <- endem_R # take from last run baseline at year 3
duration <- 365 * 6 # 3 year life span so track for 6 years
tsteps <- 2 # chunk into 2 day intervals
beta <- 1 # rate given as 1 day^-1
gamma <- 1/20 # 20 days
# trying to reduce prevalence from 1.7% down to 0.85% (half) via vaccination
pvacc <- 0.48 # placeholder, see below for zero out overwrite
sigma <- 1/(1 * 365) # average duration of immunity 1 year
mu <- 1/(3 * 365) # 3 year life span
birth <- mu # stable population so set to same as mu
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
pvacc <- 0 # no vacc scenario
initial_state_values <- c(S = (N-(endemic_I + endemic_R)),
I = endemic_I,
R = endemic_R)
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0,
pvacc_thresh = pvacc_thresh, # vaccination threshold
mu = mu, # background mortality rate
birth = birth, # birth rate
pvacc = pvacc, # vaccine rate
sigma = sigma # immunity waning rate
))## beta gamma R0 pvacc_thresh mu birth
## 1.000000000 0.050000000 20.000000000 0.950000000 0.000913242 0.000913242
## pvacc sigma
## 0.000000000 0.002739726# 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) -(S * mu) +(N * birth * (1-pvacc)) +(sigma * R)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu) +(N * birth * (pvacc)) -(sigma * R)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 64 15208.95 23433.61 261357.4 5.07 7.811 87.119 1.01393print("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
## 1 66 14819.61 23392.2 261788.2 4.94 7.797 87.263 0.9879743print("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
## 1096 2190 15273.97 19385.6 265340.4 5.091 6.462 88.447 1.018265print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
print("See equilibrium values assumed at 3 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 3 years (+1 because step is 2): "(out_v5_y3 <- output[output$time == 365*3+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 549 1096 15273.97 19385.6 265340.4 5.091 6.462 88.447 1.018265print("See equilibrium values assumed at 5 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 5 years (+1 because step is 2): "(out_v5_y5 <- output[output$time == 365*5+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15273.97 19385.6 265340.4 5.091 6.462 88.447 1.018265nicesubtitle <- "Neonatal livestock vaccination v6 48% vacc, 1 yr wane"
print(nicesubtitle)## [1] "Neonatal livestock vaccination v6 48% vacc, 1 yr wane"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
# modelling endemic infection - initial conditions for population won't matter
endemic_I <- endem_I # take from last run baseline at year 3
endemic_R <- endem_R # take from last run baseline at year 3
duration <- 365 * 6 # 3 year life span so track for 6 years
tsteps <- 2 # chunk into 2 day intervals
beta <- 1 # rate given as 1 day^-1
gamma <- 1/20 # 20 days
# trying to reduce prevalence from 1.7% down to 0.85% (half) via vaccination
pvacc <- 0.48 # placeholder, see below for threshold overwrite
sigma <- 1/(1 * 365) # average duration of immunity 1 year
mu <- 1/(3 * 365) # 3 year life span
birth <- mu # stable population so set to same as mu
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
pvacc <- pvacc_thresh # overwrite in attempt to eliminate
initial_state_values <- c(S = (N-(endemic_I + endemic_R)),
I = endemic_I,
R = endemic_R)
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0,
pvacc_thresh = pvacc_thresh, # vaccination threshold
mu = mu, # background mortality rate
birth = birth, # birth rate
pvacc = pvacc, # vaccine rate
sigma = sigma # immunity waning rate
))## beta gamma R0 pvacc_thresh mu birth
## 1.000000000 0.050000000 20.000000000 0.950000000 0.000913242 0.000913242
## pvacc sigma
## 0.950000000 0.002739726# 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) -(S * mu) +(N * birth * (1-pvacc)) +(sigma * R)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu) +(N * birth * (pvacc)) -(sigma * R)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 72 15190.56 17455.19 267354.2 5.064 5.818 89.118 1.012704print("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
## 1 74 14906.54 17428.61 267664.8 4.969 5.81 89.222 0.9937694print("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
## 1096 2190 15273.97 14534.54 270191.5 5.091 4.845 90.064 1.018265print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
print("See equilibrium values assumed at 3 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 3 years (+1 because step is 2): "(out_v6_y3 <- output[output$time == 365*3+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 549 1096 15273.97 14534.54 270191.5 5.091 4.845 90.064 1.018265print("See equilibrium values assumed at 5 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 5 years (+1 because step is 2): "(out_v6_y5 <- output[output$time == 365*5+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15273.97 14534.54 270191.5 5.091 4.845 90.064 1.018265print("R0 of 20, 1-year waning immunity, no vacc rate, year 5")## [1] "R0 of 20, 1-year waning immunity, no vacc rate, year 5"out_v5_y5## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15273.97 19385.6 265340.4 5.091 6.462 88.447 1.018265print("R0 of 20, 1-year waning immunity, a vacc rate of .95, year 5")## [1] "R0 of 20, 1-year waning immunity, a vacc rate of .95, year 5"out_v6_y5## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15273.97 14534.54 270191.5 5.091 4.845 90.064 1.018265print("reduction in prevalence with 95% neonatal vaccination rate, year 5")## [1] "reduction in prevalence with 95% neonatal vaccination rate, year 5"1 - out_v6_y5[3] / out_v5_y5[3]## I
## 914 0.2502406Would it be possible to eliminate the disease from the population using neonatal vaccination under these assumptions?
- Even with 100% vaccination, the baseline prevalence ends up at about 4.8%.
nicesubtitle <- "Neonatal livestock vaccination v7 100% vacc, 1 yr wane"
print(nicesubtitle)## [1] "Neonatal livestock vaccination v7 100% vacc, 1 yr wane"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
# modelling endemic infection - initial conditions for population won't matter
endemic_I <- endem_I # take from last run baseline at year 3
endemic_R <- endem_R # take from last run baseline at year 3
duration <- 365 * 6 # 3 year life span so track for 6 years
tsteps <- 2 # chunk into 2 day intervals
beta <- 1 # rate given as 1 day^-1
gamma <- 1/20 # 20 days
# trying to reduce prevalence from 1.7% down to 0.85% (half) via vaccination
pvacc <- 1.0 # try 100% vacc rate and comment out threshold below
sigma <- 1/(1 * 365) # average duration of immunity 1 year
mu <- 1/(3 * 365) # 3 year life span
birth <- mu # stable population so set to same as mu
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
# pvacc <- pvacc_thresh # no overwrite
initial_state_values <- c(S = (N-(endemic_I + endemic_R)),
I = endemic_I,
R = endemic_R)
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0,
pvacc_thresh = pvacc_thresh, # vaccination threshold
mu = mu, # background mortality rate
birth = birth, # birth rate
pvacc = pvacc, # vaccine rate
sigma = sigma # immunity waning rate
))## beta gamma R0 pvacc_thresh mu birth
## 1.000000000 0.050000000 20.000000000 0.950000000 0.000913242 0.000913242
## pvacc sigma
## 1.000000000 0.002739726# 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) -(S * mu) +(N * birth * (1-pvacc)) +(sigma * R)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu) +(N * birth * (pvacc)) -(sigma * R)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 72 15269.74 17133.96 267596.3 5.09 5.711 89.199 1.017982print("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
## 1 74 14982.38 17116.72 267900.9 4.994 5.706 89.3 0.9988251print("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
## 1096 2190 15273.97 14279.22 270446.8 5.091 4.76 90.149 1.018265print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
print("See equilibrium values assumed at 3 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 3 years (+1 because step is 2): "(out_v7_y3 <- output[output$time == 365*3+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 549 1096 15273.97 14279.22 270446.8 5.091 4.76 90.149 1.018265print("See equilibrium values assumed at 5 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 5 years (+1 because step is 2): "(out_v7_y5 <- output[output$time == 365*5+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15273.97 14279.22 270446.8 5.091 4.76 90.149 1.018265If an adjuvant (a vaccine promoter) was given along with the vaccine, that would extend the duration of immunity to 2.5 years on average, what vaccine coverage would be needed to reduce the baseline prevalence by half? Would it be possible to eliminate the disease from the population under these assumptions using neonatal vaccination?
- With 2.5-year waning immunity, no vaccination results in a baseline prevalence of 3.7%.
- With 2.5-year waning immunity, a 100% vaccination rate results in a baseline prevalence of 1.9%.
- The reduction in prevalence between no vaccination and a 100% vaccination rate with 2.5-year waning immunity is 47%.
nicesubtitle <- "Neonatal livestock vaccination v8 no vacc, 2.5 yr wane"
print(nicesubtitle)## [1] "Neonatal livestock vaccination v8 no vacc, 2.5 yr wane"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
# modelling endemic infection - initial conditions for population won't matter
endemic_I <- endem_I # take from last run baseline at year 3
endemic_R <- endem_R # take from last run baseline at year 3
duration <- 365 * 6 # 3 year life span so track for 6 years
tsteps <- 2 # chunk into 2 day intervals
beta <- 1 # rate given as 1 day^-1
gamma <- 1/20 # 20 days
# trying to reduce prevalence from 1.7% down to 0.85% (half) via vaccination
pvacc <- 0.0 # no vacc, no threshold overwrite below
sigma <- 1/(2.5 * 365) # average duration of immunity 2.5 years
mu <- 1/(3 * 365) # 3 year life span
birth <- mu # stable population so set to same as mu
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
# pvacc <- pvacc_thresh # no overwrite
initial_state_values <- c(S = (N-(endemic_I + endemic_R)),
I = endemic_I,
R = endemic_R)
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0,
pvacc_thresh = pvacc_thresh, # vaccination threshold
mu = mu, # background mortality rate
birth = birth, # birth rate
pvacc = pvacc, # vaccine rate
sigma = sigma # immunity waning rate
))## beta gamma R0 pvacc_thresh mu birth
## 1.00000e+00 5.00000e-02 2.00000e+01 9.50000e-01 9.13242e-04 9.13242e-04
## pvacc sigma
## 0.00000e+00 1.09589e-03# 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) -(S * mu) +(N * birth * (1-pvacc)) +(sigma * R)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu) +(N * birth * (pvacc)) -(sigma * R)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 80 15331.97 13028.96 271639.1 5.111 4.343 90.546 1.022131print("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
## 1 84 14937.56 13003.69 272058.7 4.979 4.335 90.686 0.9958376print("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
## 1096 2190 15273.97 10999.07 273727 5.091 3.666 91.242 1.018265print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
print("See equilibrium values assumed at 3 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 3 years (+1 because step is 2): "(out_v8_y3 <- output[output$time == 365*3+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 549 1096 15273.97 10999.07 273727 5.091 3.666 91.242 1.018265print("See equilibrium values assumed at 5 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 5 years (+1 because step is 2): "(out_v8_y5 <- output[output$time == 365*5+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15273.97 10999.07 273727 5.091 3.666 91.242 1.018265nicesubtitle <- "Neonatal livestock vaccination v9 100% vacc, 2.5 yr wane"
print(nicesubtitle)## [1] "Neonatal livestock vaccination v9 100% vacc, 2.5 yr wane"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
# modelling endemic infection - initial conditions for population won't matter
endemic_I <- endem_I # take from last run baseline at year 3
endemic_R <- endem_R # take from last run baseline at year 3
duration <- 365 * 6 # 3 year life span so track for 6 years
tsteps <- 2 # chunk into 2 day intervals
beta <- 1 # rate given as 1 day^-1
gamma <- 1/20 # 20 days
# trying to reduce prevalence from 1.7% down to 0.85% (half) via vaccination
pvacc <- 1.0 # 100% vacc, no threshold overwrite below
sigma <- 1/(2.5 * 365) # average duration of immunity 2.5 years
mu <- 1/(3 * 365) # 3 year life span
birth <- mu # stable population so set to same as mu
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
# pvacc <- pvacc_thresh # no overwrite
initial_state_values <- c(S = (N-(endemic_I + endemic_R)),
I = endemic_I,
R = endemic_R)
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0,
pvacc_thresh = pvacc_thresh, # vaccination threshold
mu = mu, # background mortality rate
birth = birth, # birth rate
pvacc = pvacc, # vaccine rate
sigma = sigma # immunity waning rate
))## beta gamma R0 pvacc_thresh mu birth
## 1.00000e+00 5.00000e-02 2.00000e+01 9.50000e-01 9.13242e-04 9.13242e-04
## pvacc sigma
## 1.00000e+00 1.09589e-03# 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) -(S * mu) +(N * birth * (1-pvacc)) +(sigma * R)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu) +(N * birth * (pvacc)) -(sigma * R)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 106 15273.5 5950.337 278776.2 5.091 1.983 92.925 1.018233print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## [1] time S I R still_Su preval_Inf propor_Re
## [8] Reff
## <0 rows> (or 0-length row.names)print("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
## 1096 2190 15273.97 5731.295 278994.7 5.091 1.91 92.998 1.018265print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
print("See equilibrium values assumed at 3 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 3 years (+1 because step is 2): "(out_v9_y3 <- output[output$time == 365*3+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 549 1096 15273.98 5731.291 278994.7 5.091 1.91 92.998 1.018265print("See equilibrium values assumed at 5 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 5 years (+1 because step is 2): "(out_v9_y5 <- output[output$time == 365*5+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15273.97 5731.295 278994.7 5.091 1.91 92.998 1.018265print("R0 of 20, 2.5-year waning immunity, no vacc rate, year 5")## [1] "R0 of 20, 2.5-year waning immunity, no vacc rate, year 5"out_v8_y5## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15273.97 10999.07 273727 5.091 3.666 91.242 1.018265print("R0 of 20, 2.5-year waning immunity, a 100% vacc rate, year 5")## [1] "R0 of 20, 2.5-year waning immunity, a 100% vacc rate, year 5"out_v9_y5## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15273.97 5731.295 278994.7 5.091 1.91 92.998 1.018265print("reduction in prevalence with 95% neonatal vaccination rate, year 5")## [1] "reduction in prevalence with 95% neonatal vaccination rate, year 5"1 - out_v9_y5[3] / out_v8_y5[3]## I
## 914 0.4789293Reducing baseline prevalence by half (actually only 43%)
We saw that with a 1-year waning immunity, no vaccination results in a baseline prevalence of 6.5%. Using the vaccine adjuvant, we now have a 2.5-year waning immunity; with no vaccination the baseline prevalence drops to 3.7% (43% drop). Adding 100% vaccination results in a baseline prevalence of 1.9%, which appears to be the best scenario we’ve tested so far.
Based on your results, what overall recommendation would you give to the Minister?
Assuming safety studies are mature and acceptable, I recommend neonatal vaccination ideally at 100% along with the vaccine adjuvant while continuing to set up parallel studies that can track other interventions that could help drive the prevalence further down. Also, long term studies could provide more insight.
We still need more research and studies to determine other interventions. Even a hypothetical magic adjuvant that would extend immunity to 12 years will still result in about 0.3% prevalence at 5 years, which is not zero.
nicesubtitle <- "Neonatal livestock vaccination v10 100% vacc, 12 yr wane"
print(nicesubtitle)## [1] "Neonatal livestock vaccination v10 100% vacc, 12 yr wane"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
# modelling endemic infection - initial conditions for population won't matter
endemic_I <- endem_I # take from last run baseline at year 3
endemic_R <- endem_R # take from last run baseline at year 3
duration <- 365 * 6 # 3 year life span so track for 6 years
tsteps <- 2 # chunk into 2 day intervals
beta <- 1 # rate given as 1 day^-1
gamma <- 1/20 # 20 days
# trying to reduce prevalence from 1.7% down to 0.85% (half) via vaccination
pvacc <- 1.0 # 100% vacc, no threshold overwrite below
sigma <- 1/(12 * 365) # assume magic adjuvant for average duration of immunity 12 years
mu <- 1/(3 * 365) # 3 year life span
birth <- mu # stable population so set to same as mu
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
# pvacc <- pvacc_thresh # no overwrite
initial_state_values <- c(S = (N-(endemic_I + endemic_R)),
I = endemic_I,
R = endemic_R)
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0,
pvacc_thresh = pvacc_thresh, # vaccination threshold
mu = mu, # background mortality rate
birth = birth, # birth rate
pvacc = pvacc, # vaccine rate
sigma = sigma # immunity waning rate
))## beta gamma R0 pvacc_thresh mu birth
## 1.000000e+00 5.000000e-02 2.000000e+01 9.500000e-01 9.132420e-04 9.132420e-04
## pvacc sigma
## 1.000000e+00 2.283105e-04# 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) -(S * mu) +(N * birth * (1-pvacc)) +(sigma * R)
dI <- (lambda * S) -(gamma * I) -(I * mu)
dR <- (gamma * I) -(R * mu) +(N * birth * (pvacc)) -(sigma * R)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
# tidyverse pipe to add fields including Reff
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)) )
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
## 1 0 15275 5107 279618 5.092 1.702 93.206 1.018333print("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
## 1 2 14862.5 5099.979 280037.5 4.954 1.7 93.346 0.9908334print("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
## 1096 2190 15246.59 985.4795 283767.9 5.082 0.328 94.589 1.016439print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
print("Plotting R eff")## [1] "Plotting R eff"output %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" = R0 of", parameters['R0'],
" \n mu of", round(parameters['mu'],5),"/ b of", round(parameters['birth'],5)," / step ", tsteps,
" \n pvacc of", parameters['pvacc']),
color = "Compartment",
subtitle= nicesubtitle)
print("See equilibrium values assumed at 3 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 3 years (+1 because step is 2): "(out_v10_y3 <- output[output$time == 365*3+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 549 1096 14600.78 1088.564 284310.7 4.867 0.363 94.77 0.9733851print("See equilibrium values assumed at 5 years (+1 because step is 2): ")## [1] "See equilibrium values assumed at 5 years (+1 because step is 2): "(out_v10_y5 <- output[output$time == 365*5+1,])## time S I R still_Su preval_Inf propor_Re Reff
## 914 1826 15402.47 978.0193 283619.5 5.134 0.326 94.54 1.026831