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.

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.

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?

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.)

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:

- What does the parameter \(\beta\) represent? Describe it in words.
- Write an R program that approximates the solution to this differential equation numerically (you can adjust the code from above).
- Visualize the solutions for different values of \(\beta\) (or use the Shiny app). What influence does it have?
- What are the equilibrium points of this system? Are they stable or unstable?
- At the start of an epidemic, how do the number of infected people approximately evolve? How does \(\beta\) influence this?
- (Optional:) Analytically derive a solution for this linear differential equation.

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
)
}
```

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\).

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\).

*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?

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?

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?

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} \]

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?

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.

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!

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?

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

- Estimate the parameters of the
**SIR**model (using least squares or any other estimation procedure you can think of) - Visualize the raw data together with the fit. How does it look like?
- Are the parameter estimates you get sensible? Why / why not?