03. Discrete-time deterministic models

In this practical, we will implement a simple SIR model using difference equations.

Practical 1. A simple SIR model

We will begin by implementing the SIR model that was introduced in the slides for this session.

Start by copying the code below into a new R source file in RStudio:

# Times for evaluating the model
time_sir <- seq(0, 20, by = 1)

# Storage for S, I, R at times time_sir
y_sir <- matrix(data = NA,
                nrow = length(time_sir),
                ncol = 3)

# Function to return our update vector with three elements:
update_sir <- function(t, y, parms)
{
    # Get parameter values from parms vector
    beta  <- parms["beta"]
    gamma <- parms["gamma"]
    
    # Get current state variables from y vector
    S <- y[1]
    I <- y[2]
    R <- y[3]

    # FILL IN: out should be a numeric vector with three elements, 
    # with updates to add to S, I, and R respectively
    out <- ...

    return (out)
}

# Set parameter values
parms_sir <- c(beta = 1.3, gamma = 0.23)

# FILL IN: Set initial values for time zero
y_sir[1, ] <- c(..., ..., ...)

# Run each step of the model from the second time step onward
for (i in 2:nrow(y_sir))
{
    # FILL IN: assign new value to y_sir[i, ] using the update_sir function
    y_sir[i, ] <- ...
}

# Plot the number of infectious individuals over time.
# type = "s" makes a "stair-step" line plot suitable for difference equations
plot(x = time_sir, y = y_sir[, 2], type = "s")

Then, fill in the blanks indicated by ... after FILL IN comments, and run the completed code. Use this to answer the questions:

  1. At approximately what time does the peak in infectious population occur and what proportion of the population is infectious?

  2. Approximately what proportion of the population gets infected?

Now, change the mean time spent infectious from 4.35 days to 2 days, keeping the rate of transmission the same.

  1. With these new parameters, at approximately what time does the peak in infectious population occur and what proportion of the population is infectious?

  2. Approximately what proportion of the population gets infected? You may need to change the simulation to run for more than 20 time steps.

Finally, change the mean time spent infectious back to 4.35 days and set the transmission rate to be half what is was before.

  1. At approximately what time does the peak in infectious population occur and what proportion of the population is infectious?

  2. Approximately what proportion of the population gets infected? You may need to change the simulation to run for more than 20 time steps.

Practical 2. Extending the SIR model

Incorporating births and deaths

Working from your code in Practical 1 (use the solution to practical 1 if you are stuck), adapt your SIR model to incorporate birth of new susceptibles, at a rate proportional to the sum of the total population of S + I + R individuals. Balance these new births with deaths occurring in each of the S, I, and R groups, with both the per capita birth rate and the per capita death rate being \(\delta = 0.01\) (\(\delta\) = “delta”).

You will need to add a new parameter, delta = 0.01, and add the new transitions for births and deaths to your update_sir() function.

Visualising for the whole population

Calculate \(N(t) = S(t) + I(t) + R(t)\) the total number of individuals.

  1. Make a plot of \(S(t)\), \(I(t)\), \(R(t)\) and \(N(t)\). Your function \(N(t)\) should be constant at 1 for all values of \(t\). If this is not the case, ensure the model contains births of new S proportional to N, and deaths of each of S, I, and R.

  2. Change the code so the model runs for 100 time steps instead of 20. Compare the number of infectious people over time when delta = 0 (no births and deaths, as before) versus when delta = 0.01. What explains the difference?

  3. Describe the parameters of the model, what they represent, and whether the assumptions they represent are realistic.

Solutions to this practical can be accessed here.