# Individual-based SARS-CoV-2 transmission model, practical 1
library(ggplot2)
## Model parameters
beta <- 0.5 # Transmission parameter
delta <- 1 / 2.5 # Rate of transitioning out of latent state
gamma <- 1 / 5 # Rate of transitioning out of infectious state
omega <- 1 / 180 # Rate of waning immunity
dt <- 1 # Time step of simulation (1 day)
days <- 365 # Duration of simulation (365 days)
steps <- days / dt # Total number of time steps
n <- 1000 # Population size
## Data frame to store simulation results
results <- data.frame(ts = 1:steps, S = 0, E = 0, I = 0, R = 0)
## Initialize simulation
# Set the seed for the pseudorandom number generator, for reproducibility
set.seed(12345)
# Since this is an individual-based model, we track the properties of all n
# individuals in the simulation. One kind of property we can track is a state,
# such as S (susceptible), E (exposed), I (infectious), or R (recovered). We
# will store each individual's state as a string, either "S", "E", "I", or "R".
state <- rep("S", n) # Each individual's state: start with all susceptible
state[1:10] <- "E" # Start 10 individuals in the "exposed" state
## Run simulation
# We'll use the built-in function txtProgressBar to track the simulation's
# progress. Really helps for planning coffee breaks! It needs to know the
# minimum and maximum values to expect, and style = 3 tells it to report the
# percentage complete.
bar <- txtProgressBar(min = 1, max = steps, style = 3)
# Loop over each time step . . .
for (ts in 1:steps) {
# Calculate the force of infection
lambda <- beta * sum(state == "I") / n
# Loop through each individual . . .
for (i in 1:n) {
if (state[i] == "S") {
# Transition S -> E (infection) at rate lambda
if (runif(1) < 1 - exp(-lambda * dt)) {
state[i] <- "E"
}
} else if (state[i] == "E") {
# Transition E -> I (latent to infectious) at rate delta
if (runif(1) < 1 - exp(-delta * dt)) {
state[i] <- "I"
}
} else if (state[i] == "I") {
# Transition I -> R (infectious to recovered) at rate gamma
if (runif(1) < 1 - exp(-gamma * dt)) {
state[i] <- "R"
}
} else if (state[i] == "R") {
# Transition R -> S (waning of immunity) at rate omega
if (runif(1) < 1 - exp(-omega * dt)) {
state[i] <- "S"
}
}
}
# Save population state for this time step
results[ts, "S"] <- sum(state == "S")
results[ts, "E"] <- sum(state == "E")
results[ts, "I"] <- sum(state == "I")
results[ts, "R"] <- sum(state == "R")
# Update progress bar; close progress bar if we are finished
setTxtProgressBar(bar, ts)
if (ts == steps) {
close(bar)
}
}
## Plot simulation results
ggplot(results) +
geom_line(aes(x = ts, y = S, colour = "S")) +
geom_line(aes(x = ts, y = E, colour = "E")) +
geom_line(aes(x = ts, y = I, colour = "I")) +
geom_line(aes(x = ts, y = R, colour = "R"))08. Stochastic individual-based models (solutions)
Practical 1. An individual-based SEIR model of SARS-CoV-2 transmission
Practical 2. Adding more complex dynamics to the model
# Individual-based SARS-CoV-2 transmission model, practical 2
library(ggplot2)
## Model parameters
beta <- 0.5 # Transmission parameter
iota <- 1e-5 # Importation rate
wane <- 0.05 # Rate of antibody waning
dt <- 1 # Time step of simulation (1 day)
days <- 365 * 2 # Duration of simulation (2 years)
steps <- days / dt # Total number of time steps
n <- 1000 # Population size
## Some helper functions
# Calculates infectiousness as a function of state and age: zero if state is
# not "I"; nonzero if state is "I", and slightly decreasing with age
infectiousness <- function(state, age) {
ifelse(state == "I", 1.25 - age / 160, 0)
}
# Calculates susceptibility of individuals with antibody level(s) ab
susceptibility <- function(ab) {
pnorm(ab, 5, 1, lower.tail = FALSE)
}
# Generates n random delays from the latent-period distribution
# (approximately 2 days, on average)
latent_delay <- function(n) {
rlnorm(n, meanlog = 0.5, sdlog = 0.6)
}
# Generates n random delays from the infectious-period distribution
# (approximately 5 days, on average)
infectious_delay <- function(n) {
rlnorm(n, meanlog = 1.5, sdlog = 0.5)
}
# Generates n random increments to antibody levels following recovery
ab_increment <- function(n) {
rnorm(n, mean = 12, sd = 2)
}
## Data frame to store simulation results
results <- data.frame(ts = 1:steps, S = 0, E = 0, I = 0, AMeanU = 0, AMeanV = 0)
## Initialize simulation
# Set the seed for the pseudorandom number generator, for reproducibility
set.seed(12345)
# Initialize state variables
state <- rep("S", n) # Each individual's state: start with all susceptible
age <- runif(n, 0, 80) # Each individual's age: random distribution from 0 to 80
delay <- rep(0, n) # Delay for latent and infectious periods
antib <- rep(0, n) # Antibody concentration for each individual
vacc <- rep(FALSE, n) # Vaccinated status
state[1:10] <- "E" # Start 10 individuals in the "exposed" state
## Run simulation
# Initialize progress bar
bar <- txtProgressBar(min = 1, max = steps, style = 3)
# Loop over each time step . . .
for (ts in 1:steps) {
# Calculate the force of infection
lambda <- beta * sum(infectiousness(state, age)) / n + iota
# Loop through each individual . . .
for (i in 1:n) {
# Update individual i's non-state variables
# Time remaining in latent/infectious periods
delay[i] <- delay[i] - dt
# Antibody waning
antib[i] <- antib[i] - wane * dt
# Vaccination at time step 300 for over-40s
if ((ts == 300) && (age[i] >= 40)) {
vacc[i] <- TRUE
antib[i] <- antib[i] + 2 * ab_increment(1)
}
# Update individual i's state
if (state[i] == "S") {
# Transition S -> E (infection) at rate lambda
if (runif(1) < 1 - exp(-lambda * dt)) {
if (runif(1) < susceptibility(antib[i])) {
state[i] <- "E"
delay[i] <- latent_delay(1)
}
}
} else if (state[i] == "E") {
# Transition E -> I (latent to infectious)
if (delay[i] < 0) {
state[i] <- "I"
delay[i] <- infectious_delay(1)
}
} else if (state[i] == "I") {
# Transition I -> S (infectious to susceptible)
if (delay[i] < 0) {
state[i] <- "S"
antib[i] <- antib[i] + ab_increment(1)
}
}
}
# Save population state for this time step
results[ts, "S"] <- sum(state == "S")
results[ts, "E"] <- sum(state == "E")
results[ts, "I"] <- sum(state == "I")
results[ts, "AMeanU"] <- mean(antib[!vacc])
results[ts, "AMeanV"] <- mean(antib[vacc])
# Update progress bar; close progress bar if we are finished
setTxtProgressBar(bar, ts)
if (ts == steps) {
close(bar)
}
}
## Plot simulation results
ggplot(results) +
geom_line(aes(x = ts, y = S, colour = "S")) +
geom_line(aes(x = ts, y = E, colour = "E")) +
geom_line(aes(x = ts, y = I, colour = "I"))
ggplot(results) +
geom_line(aes(x = ts, y = AMeanU, colour = "Unvaccinated")) +
geom_line(aes(x = ts, y = AMeanV, colour = "Vaccinated")) +
labs(x = "Time step", y = "Mean antibody level")Practical 3. Optimizing the model to run faster
# Individual-based SARS-CoV-2 transmission model, practical 3
library(ggplot2)
## Model parameters
beta <- 0.5 # Transmission parameter
iota <- 1e-5 # Importation rate
wane <- 0.05 # Rate of antibody waning
dt <- 1 # Time step of simulation (1 day)
days <- 365 * 4 # Duration of simulation (4 years)
steps <- days / dt # Total number of time steps
n <- 5000 # Population size
## Some helper functions
# Calculates infectiousness as a function of state and age: zero if state is
# not "I"; nonzero if state is "I", and slightly decreasing with age
infectiousness <- function(state, age) {
ifelse(state == "I", 1.25 - age / 160, 0)
}
# Calculates susceptibility of individuals with antibody level(s) ab
susceptibility <- function(ab) {
pnorm(ab, 5, 1, lower.tail = FALSE)
}
# Generates n random delays from the latent-period distribution
# (approximately 2 days, on average)
latent_delay <- function(n) {
rlnorm(n, meanlog = 0.5, sdlog = 0.6)
}
# Generates n random delays from the infectious-period distribution
# (approximately 5 days, on average)
infectious_delay <- function(n) {
rlnorm(n, meanlog = 1.5, sdlog = 0.5)
}
# Generates n random increments to antibody levels following recovery
ab_increment <- function(n) {
rnorm(n, mean = 12, sd = 2)
}
## Data frame to store simulation results
results <- data.frame(ts = 1:steps, S = 0, E = 0, I = 0, AMeanU = 0, AMeanV = 0)
## Initialize simulation
# Set the seed for the pseudorandom number generator, for reproducibility
set.seed(12345)
# Initialize state variables
state <- rep("S", n) # Each individual's state: start with all susceptible
age <- runif(n, 0, 80) # Each individual's age: random distribution from 0 to 80
delay <- rep(0, n) # Delay for latent and infectious periods
antib <- rep(0, n) # Antibody concentration for each individual
vacc <- rep(FALSE, n) # Vaccinated status
state[1:10] <- "E" # Start 10 individuals in the "exposed" state
## Run simulation
# Initialize progress bar
bar <- txtProgressBar(min = 1, max = steps, style = 3)
# Loop over each time step . . .
for (ts in 1:steps) {
# Calculate the force of infection
lambda <- beta * sum(infectiousness(state, age)) / n + iota
##### NOTE - There is no inner loop over individuals anymore!
# Update non-state variables (for all individuals simultaneously)
# Time remaining in latent/infectious periods
delay <- delay - dt
# Antibody waning
antib <- antib - wane * dt
# Vaccination at time step 300 for over-40s
if (ts == 300) {
vacc[age >= 40] <- TRUE
antib[vacc] <- antib[vacc] + 2 * ab_increment(sum(vacc))
}
# Update state variables (for all individuals simultaneously)
##### trE selects all individuals who will transition states from S to E.
trE <- (state == "S") & (runif(n) < 1 - exp(-lambda * dt)) &
(runif(n) < susceptibility(antib))
trI <- (state == "E") & (delay < 0)
trS <- (state == "I") & (delay < 0)
# Do state transitions
# transition S -> E
state[trE] <- "E"
delay[trE] <- latent_delay(sum(trE))
# transition E -> I
state[trI] <- "I"
delay[trI] <- infectious_delay(sum(trI))
# transition I -> S
state[trS] <- "S"
antib[trS] <- antib[trS] + ab_increment(sum(trS))
# Save population state for this time step
results[ts, "S"] <- sum(state == "S")
results[ts, "E"] <- sum(state == "E")
results[ts, "I"] <- sum(state == "I")
results[ts, "AMeanU"] <- mean(antib[!vacc])
results[ts, "AMeanV"] <- mean(antib[vacc])
# Update progress bar; close progress bar if we are finished
setTxtProgressBar(bar, ts)
if (ts == steps) {
close(bar)
}
}
## Plot simulation results
ggplot(results) +
geom_line(aes(x = ts, y = S, colour = "S")) +
geom_line(aes(x = ts, y = E, colour = "E")) +
geom_line(aes(x = ts, y = I, colour = "I"))
ggplot(results) +
geom_line(aes(x = ts, y = AMeanU, colour = "Unvaccinated")) +
geom_line(aes(x = ts, y = AMeanV, colour = "Vaccinated")) +
labs(x = "Time step", y = "Mean antibody level")Return to the practical here.