Introduction

These are exercises for the workshop Epidemiological Modelling: An Introduction, given at the QHELP (Higher Education Learning Platform for Quantitative Thinking) seminar in Leuven, March 2022. They can be done in any order, depending on your experience and interest, but if you have no strong preference I suggest you tackle them in sequence. The slides, including the Shiny apps, are available from https://r.qhelp.eu/qhelp/Presentations/Epidemiology-Workshop.

Differential equations

We have the following linear differential equation:

\[ \tag{1} \frac{dx}{dt} = f(x) = r \cdot x \enspace . \]

Recall Euler’s approximation:

\[ \tag{2} x(t + \Delta_t) = x(t) + \Delta_t \cdot f(x(t)) \enspace . \]

Recall also that the solution to this differential equation is given by:

\[ \tag{3} x(t) = x_0 e^{r t} \enspace . \]

As a shortcut, we write \(\dot{x} \equiv \frac{dx}{dt}\). Define an equilibrium point \(x^{\star}\) as a point where \(\dot{x} = 0\). An equilibrium is stable if the system returns to it after perturbations, and unstable if it does not.

Numerical approximation

Using pen and paper

Make sure you understand the idea behind the numerical approximation. It helps to do some calculations by hand. Approximate the solution to this differential equation with \(x(0) = 1\) for \(\Delta_t = 0.50\) and \(\Delta_t = 0.25\) until you reach \(x(3)\), assuming \(r = 1\). What differences in value do you observe for \(x(1), x(2)\), and \(x(3)\)? Why is that?

Using R

Now write an R program that approximates the solution to Equation (1). As a starting point, you might be helped by this code:

solve_de <- function(f, x0, delta_t = 0.01, nr_iter = 1000) {
  
  x <- rep(x0, nr_iter + 1)
  
  for (i in seq(1, nr_iter)) {
    x[i + 1] <- x[i] + delta_t * f(x[i])
  }
  
  time <- seq(0, nr_iter * delta_t, delta_t)
  res <- cbind(time, x)
  res
}

First, make sure that your hand calculations from the previous section are correct.

Second, simulate from the system using \(\Delta_t \in \{0.01, 0.50, 1\}\) until the same 5 units of time pass with \(x_0 = 1\). (Note: the argument ‘nr_iter’ in the function changes depending on \(\Delta_t\)!). Visualize the solutions in the same plot together with the analytical solution. How does the quality of the approximation to the true function differ? What explains this?

Further exercises:

  • What influence does the parameter \(r\) have? Study cases where \(r < 0\) and \(r > 0\).
  • What are the equilibrium points of this system? Are they stable or unstable?
  • (Optional:) Rederive the solution given in Equation (3).

(If you are interested in learning more about linear differential equations, you can check out this blog post.)

SI model

Recall the basic SI model:

\[ \begin{aligned} \frac{dS}{dt} &= -\beta I \cdot S \\[0.50em] \frac{dI}{dt} &= \beta I \cdot S \enspace , \end{aligned} \]

where we assume that \(S + I = 1\) and which therefore reduces to:

\[ \tag{4} \frac{dI}{dt} = \beta I \cdot (1 - I) \]

Here are a few exercises:

SIR model

Recall the basic SIR model:

\[ \begin{aligned} \frac{dS}{dt} &= -\beta I S \\[0.50em] \frac{dI}{dt} &= \beta I S \mathbin{\color{black}{-}} \color{black}{\gamma I} \\[0.50em] \color{black}{\frac{dR}{dt}} &\mathbin{\color{black}{=}} \color{black}{\gamma I} \end{aligned} \]

Define \(R_0 \equiv \frac{\beta}{\gamma}\). This following code might help you for the exercises below:

library('RColorBrewer')

solve_SIR <- function(S0, I0, beta = 1, gamma = 1, delta_t = 0.01, nr_iter = 1000) {
  
  res <- matrix(NA, nrow = nr_iter + 1, ncol = 4)
  colnames(res) <- c('S', 'I', 'R', 'Time')
  
  res[1, ] <- c(S0, I0, 1 - S0 - I0, 0)
  
  dS <- function(S, I) -beta * I * S
  dI <- function(S, I)  beta * I * S - gamma * I
  
  for (i in seq(1, nr_iter)) {
    ???
  }
  
  res[, 3] <- 1 - res[, 1] - res[, 2]
  res
}


plot_epidemic <- function(res, main = '', legend_x = 50) {
  cols <- brewer.pal(3, 'Set1')
  time <- res[, 4]
  
  matplot(
    time, res[, 1:3], type = 'l', col = cols, axes = FALSE, lty = 1, lwd = 2,
    ylab = 'Subpopulations', xlab = 'Time t', xlim = c(0, tail(time, 1)),
    ylim = c(0, 1), main = main, cex.main = 1.75, cex.lab = 1.5,
    xaxs = 'i', yaxs = 'i'
  )
  
  axis(1, cex.axis = 1.25)
  axis(2, las = 2, cex.axis = 1.25)
  legend(
    legend_x, 0.65, col = cols, legend = c('S', 'I', 'R'),
    lty = 1, lwd = 2, bty = 'n', cex = 1.5
  )
}

SIR dynamics

Here are a few exercises:

  • What does the parameter \(\gamma\) represent? Describe it in words.
  • Finalize the R program that approximates the solution to this differential equation numerically.
  • Suppose the units of choice is days (i.e., \(\Delta_t = 1\) is one day). Simulate the system with \(\beta = 0.50\), \(\gamma = 0.20\), and \(I_0 = 0.01\) for 100 days and visualize the results. What do you observe?
  • Simulate from the system once using \(\beta = 2\) and \(\gamma = 1\) for 2000 iterations and once using \(\beta = 1\) and \(\gamma = 0.50\) using 1000 iterations, both with \(\Delta_t = 0.01\). Visualize the results. What do you observe and why?
  • At the start of the epidemic were \(S = 1\), how do the infectious fraction approximately evolve with time? How do \(\gamma\) and \(R_0\) influence this?

Here are some exercises with the Shiny app for the SIR dynamics. Use the Shiny app and explain what you observe and why:

  • Fix \(I_0 = 1\) and \(\beta = 1\) and vary \(\gamma\) from \(0\) to \(2\). (Optional: Study the differential equation for \(I\) when \(S\) is very close to zero and derive its approximate solution. What do you observe?)
  • Fix \(I_0 = 0.20\) and \(\beta = 1\) and vary \(\gamma\) from \(0\) to \(1.50\).

SIR vector field

Here are another set of exercises using the Shiny app for the SIR vector field. Use the Shiny app and explain what you observe and why. (Always enforce that \(S + I + R = 1\) when changing sliders):

  • Setting \(\beta = 1\) and \(\gamma = 0\), play around with values for \(I_0\) and \(S_0\). Where is the trajectory being pulled towards? Why? In terms of values for \(S\), \(I\), and \(R\), what are the equilibrium points for this system?
  • Now set \(\gamma = 0.50\), \(I_0 = 0.20\), and \(S_0 = 0.80\). Where is the trajectory being pulled towards? Do the number of infected people increase or decrease? How does this change as you adjust \(S_0\) towards \(S_0 = 0.40\)? At what values for \(S_0\) does the number of infected do not increase anymore?
  • Make sense of the observations above by noting that outbreaks (i.e., increases in \(I\)) are possible whenever \(R_t > 1\) and using the relationships between \(R_t\) and \(R_0\).

SIR nullclines

Nullclines allow us to organize the vector field in a systematic way; these are the lines where one part of the system does not change. For example, the \(I\)-nullclines are given by the lines for which \(\frac{dI}{dt} = 0\). Correspondingly, the \(S\)-nullclines are given by the lines for which \(\frac{dS}{dt} = 0\).

Find the \(S\)-nullclines, using:

\[ \frac{dS}{dt} = -\beta I S = 0 \]

And then find the \(I\)-nullclines, using:

\[ \frac{dI}{dt} = \beta I S - \gamma I = 0 \]

In both cases, your solution should be two nullclines. There is one nullcline that is especially interesting. Which one is it, and why? What does it say?

Threshold phenomena

Here are some more exercises:

  • What is the basic reproductive number and what is the effective reproductive number?
  • Let \(p\) denote the proportion of people being immune, and recall that \(R_t = R_0 \cdot S\). At what threshold of immunity \(p_c\) do outbreaks become impossible?
  • Compute the herd immunity threshold for pathogens with \(R_0 \in \{0.50, 1, 2, 5, 10, 20\}\). What do you observe and why?
  • As a function of \(R_0\), what fraction of the populaton needs to be susceptible in order for outbreaks to be possible?

Size of the outbreak

In this exercise, you should run the model for a sufficiently long time (\(10,000\) iterations with \(\Delta_t = 0.01\) should be enough) so that the values for \(S\), \(I\), and \(R\) do not change anymore. The model has then reached equilibrium. Fix \(\gamma = 0.50\), and run the model using values \(R_0 \in \{0.01, 0.05, \ldots, 4.95, 5\}\) and \((S_0, I_0) = (0.999, 0.001)\). The size of the outbreak is given by the number of people who have become infected and ultimately end up in the \(R\) compartment.

Plot the fraction of people recovered across the values for \(R_0\). What do you observe?

SIRS model

Recall the basic SIRS model:

\[ \begin{aligned} \frac{dS}{dt} &= -\beta I S \mathbin{\color{black}{+}} \color{black}{\mu R} \\[0.50em] \frac{dI}{dt} &= \beta I S - \gamma I \\[0.50em] \frac{dR}{dt} &= \gamma I \mathbin{\color{black}{-}} \color{black}{\mu R} \end{aligned} \]

SIRS dynamics

Here are a few exercises:

  • What does the parameter \(\mu\) represent? Describe it in words.
  • Extend the SIR R program to approximates the solution to this differential equation numerically.
  • Suppose the units of choice is days (i.e., \(\Delta_t = 1\) is one day). Simulate the system with \(\beta = 0.50\), \(\gamma = 0.20\), \(\mu = 0.02\), and \(I_0 = 0.01\) for 100 days and visualize the result (you can use the ‘plot_epidemic’ from above). What do you observe?
  • Simulate it for 500 days. What do you observe now?

Here are some exercises with the Shiny app for the SIRS dynamics. Look at the system for 400 days. Use the Shiny app and explain what you observe and why:

  • Fix \(I_0 = 0.01\), \(\beta = 0.50\), \(\gamma = 0.20\), and vary \(\mu\) from \(0\) to \(0.05\). What effect does this have? Why?
  • Fix \(I_0 = 0.01\), \(\beta = 0.50\), \(\mu = 0.02\), and vary \(\gamma\) from \(0\) to \(0.50\). What effect does this have? Why? What are especially important effects when changing \(\gamma\)?
  • Fix \(I_0 = 0.01\), \(\gamma = 0.20\), \(\mu = 0.02\), and vary \(\beta\) from \(0.50\) to \(1\). What effect does this have? Why?
  • Fix \(I_0 = 0.01\), \(\beta = 0.50\), \(\gamma = 0.20\), and \(\mu = 0.40\). What dynamics do unfold? Why?

SIRS vector field

Open ended: Play around with the Shiny app for the SIRS vector field. Fixing all other parameters, vary one parameter at a tiime. Make sure that you understand what is happening, and why.

SIRS nullclines

The \(S\)-nullcline is the line where:

\[ \frac{dS}{dt} = -\beta I S + \mu R = 0 \]

The \(I\)-nullcline is the line where:

\[ \frac{dI}{dt} = \beta I S - \gamma I = 0 \]

Find the nullclines! Recall that at equilibrium, both lines intersect. Answer the following questions:

  • How will the equilibrium shift as \(\beta\) increases?
  • How will the equilibrium shift as \(\mu\) increases?

Check your answers using the SIRS vector field Shiny app!

Equilibria of the SIRS model

An equilibrium (or fixed point) is a point \((S^{\star}, I^{\star})\) where

\[ \begin{aligned} \frac{dS}{dt} &= -\beta I S + \mu (1 - S - I) &= 0 \\ \frac{dI}{dt} &= \beta I S - \gamma I &= 0 \end{aligned} \]

For which points \((S^{\star}, I^{\star})\) does this occur? What happens to the fixed points as \(R_0\) changes?

Parameter estimation

Use the covidregionaldata package to get COVID-19 case reports for a country of your choice. Select the period that coincides with the first wave, starting when the first case was reported. For Austria, this can be done like this:

library('covidregionaldata')

aut <- get_national_data(countries = 'Austria')
start <- which(aut$cases_new > 0)[1]
end <- which(aut$date == '2020-05-15') # By visual inspection

austria <- aut[seq(start, end), ]
plot(
  austria$date, austria$cases_new,
  xlab = 'Date', ylab = 'Reported Cases', pch = 20, axes = FALSE,
  ylim = c(0, 1200), main = '1st COVID-19 Wave in Austria', font.main = 1
)
axis.Date(
  1,
  at = c('2020-02-26','2020-03-15', '2020-04-15', '2020-05-15'),
  format = '%m/%d/%y'
)
axis(2)

Here are the exercises