08_StochasticIBM_solutions.qmd

---
title: "08. Stochastic individual-based models (solutions)"
---

## Practical 1. An individual-based SEIR model of SARS-CoV-2 transmission

```{r, eval = FALSE, message = FALSE}
# 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"))
```

## Practical 2. Adding more complex dynamics to the model

```{r, eval = FALSE, message = FALSE}
# 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

```{r, eval = FALSE, message = FALSE}
# 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](08_StochasticIBM_practical.qmd).**