03_DiscreteDeterministic_solutions.qmd

---
title: "03. Discrete-time deterministic models: Solutions"
---

## Practical 1. A simple SIR model

```{r, eval = T, echo = T}
# 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]

    # Numeric vector with updates to add to S, I, and R respectively
    out <- c(-beta * S * I,
              beta * S * I - gamma * I,
                             gamma * I)

    return (out)
}

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

# Set initial values for time zero
y_sir[1, ] <- c(0.99, 0.01, 0)

# Run each step of the model from the second time step onward
for (i in 2:nrow(y_sir))
{
	# Assign new value to y_sir[i, ] using the update_sir function
    y_sir[i, ] <- y_sir[i - 1, ] +
        update_sir(time_sir[i],
                   y_sir[i - 1, ],
                   parms_sir)
}

# 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")
```

**Questions 1-6** ask:

-  At approximately what time does the peak in infectious population occur and what proportion of the population is infectious? and
-  Approximately what proportion of the population gets infected? 

for three different sets of `beta` (transmission rate) and `gamma` (recovery rate) in the model. The solutions are summarized in the table and figure below:

| Parameters                  | Peak in infectious population | Total proportion infected |
|-----------------------|----------------------------|----------------------|
| `beta = 1.3,  gamma = 0.23` | \(1\) t = 7, I = 0.60         | \(2\) 0.9997              |
| `beta = 1.3,  gamma = 0.5`  | \(3\) t = 8, I = 0.30         | \(4\) 0.943               |
| `beta = 0.65, gamma = 0.23` | \(5\) t = 14, I = 0.31        | \(6\) 0.943               |

```{r, echo = FALSE, eval = TRUE, output = TRUE}

do_sir <- function(beta, gamma)
{
	parms_sir <- c(beta = beta, gamma = gamma)
	y_sir <- matrix(data = NA,
                nrow = length(time_sir),
                ncol = 3)
	y_sir[1, ] <- c(0.99, 0.01, 0)

	# Run each step of the model from the second time step onward
	for (i in 2:nrow(y_sir))
	{
		# Assign new value to y_sir[i, ] using the update_sir function
	    y_sir[i, ] <- y_sir[i - 1, ] +
	        update_sir(time_sir[i],
	                   y_sir[i - 1, ],
	                   parms_sir)
	}
	
	return (y_sir)
}

y_sir_1 <- do_sir(1.3, 0.23)
y_sir_2 <- do_sir(1.3, 0.5)
y_sir_3 <- do_sir(0.65, 0.23)

# Plot all
plot(x = time_sir, y = y_sir_1[, 2], type = "s", col = "red", xlab = "Time (days)", ylab = "Infectious people")
lines(x = time_sir, y = y_sir_2[, 2], type = "s", col = "orange")
lines(x = time_sir, y = y_sir_3[, 2], type = "s", col = "blue")
legend("topright", legend = expression(list(beta==1.3, gamma==0.23), list(beta==1.3, gamma==0.5), list(beta==0.65, gamma==0.23)),
	col = c("red", "orange", "blue"), lty = 1, bty = "n")

```

When we start with $\beta=1.3, \gamma=0.23$ the peak occurs on day 7 with 60% of the population infectious. Almost everyone in the population is infected by the end of the simulation.

When we increase `gamma` from $\gamma=0.23$ to $\gamma=0.5$, this increase in the recovery rate corresponds to a shorter infectious period (2 days instead of 4.35 days). This means people spend less time in the $I$ compartment, which decreases the reproduction number by around $1/2$, and so the peak proportion of the population who are infectious is smaller (by around $1/2$). But since the generation interval is also shorter by around $1/2$, the peak doesn't happen substantially later.

When instead we decrease `beta` from $\beta=1.3$ to $\beta=0.65$, this decrease in the transmission rate reduces the reproduction number, but does not change the generation interval. The peak proportion of the population who are infectious is smaller relative to when $\beta=1.3, \gamma=0.23$. This time, the peak is also shifted later, because the lower reproduction number is not "compensated" in this regard by a shorter generation interval.

For questions (4) and (6), the reproduction number is lower by around the same amount relative to question (2), so the total proportion infected is lower by a similar amount in both (4) and (6).

## Practical 2. Extending the SIR model

### Incorporating births and deaths

```{r, echo = T, eval = T, output = T}
# 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"]
    delta <- parms["delta"] ## NEW: delta parameter
    
    # Get current state variables from y vector
    S <- y[1]
    I <- y[2]
    R <- y[3]
    N <- S + I + R ## NEW: calculation of N

    # Numeric vector with updates to add to S, I, and R respectively
    ## NEW: transitions for birth and death with parameter delta
    out <- c(-S * delta - beta * S * I + N * delta,
             -I * delta + beta * S * I - gamma * I,
             -R * delta                + gamma * I)

    return (out)
}

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

# Set initial values for time zero
y_sir[1, ] <- c(0.99, 0.01, 0)

# Run each step of the model from the second time step onward
for (i in 2:nrow(y_sir))
{
	# Assign new value to y_sir[i, ] using the update_sir function
    y_sir[i, ] <- y_sir[i - 1, ] +
        update_sir(time_sir[i],
                   y_sir[i - 1, ],
                   parms_sir)
}

# 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")
```

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

**Answer.** Here is one possible method:

```{r, eval = T, output = T}
# Add total population N(t) to y_sir matrix
y_sir <- cbind(y_sir, rowSums(y_sir))

# Plot each model component
par(mfrow = c(2, 2))
components = c("S", "I", "R", "N")

for (i in 1:ncol(y_sir)) {
    plot(y_sir[, i] ~ time_sir, type = "s", xlab = "Time (years)", ylab = components[i])
}
par(mfrow = c(1, 1))
```

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?

```{r, eval = T, output = F, echo = F}
# Setup
time_sir <- seq(0, 100, by = 1)
y_sir <- matrix(data = NA, nrow = length(time_sir), ncol = 3)

# First run
parms_sir <- c(beta = 1.3, gamma = 0.23, delta = 0.01)
y_sir[1, ] <- c(0.99, 0.01, 0)
for (i in 2:nrow(y_sir)) {
    y_sir[i, ] <- y_sir[i - 1, ] + update_sir(time_sir[i], y_sir[i - 1, ], parms_sir)
}
y_sir_01 <- y_sir

# Second run
parms_sir <- c(beta = 1.3, gamma = 0.23, delta = 0)
y_sir[1, ] <- c(0.99, 0.01, 0)
for (i in 2:nrow(y_sir)) {
    y_sir[i, ] <- y_sir[i - 1, ] + update_sir(time_sir[i], y_sir[i - 1, ], parms_sir)
}
y_sir_00 <- y_sir
```

**Answer.** When the birth rate is 0.01, we eventually see a second outbreak occur around year 50. This is because new susceptibles enter the population via birth, which eventually furnishes enough susceptibles to make the effective reproduction number $R_t>1$ and allow a new outbreak to occur.

Here is one way of plotting infectious individuals over time for `delta = 0` versus `delta = 0.01`, assuming that the solution `y_sir` for `delta = 0` is stored in the new variable `y_sir_00`, and `y_sir` for `delta = 0.01` is stored in the new variable `y_sir_01`:

```{r, eval = T, output = T, echo = T}
plot(time_sir, y_sir_00[, 2], col = "black", type = "s", xlab = "Time (years)", ylab = "I(t)")
lines(time_sir, y_sir_01[, 2], col = "red", type = "s", xlab = "Time (years)")
legend("topright", c("delta = 0", "delta = 0.01"), col = c("black", "red"), lty = 1, title = "Birth rate")
```

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

**Answer:**

The transmission rate `beta` represents the rate at which an individual makes "effective contact" with another individual, where effective contact results in a new infection if the contacted individual is infectious and the contacting individual is susceptible.

The recovery rate `gamma` represents the recovery rate, such that `1/gamma` is the infectious period.

The birth rate `delta` captures the rate of new births and also the rate of deaths.

There's an implicit assumption in the model that transmission is not passed to newborns, i.e., only susceptibles are born. This is likely a reasonable assumption to make for many diseases.

As we are dealing with proportions of the total population, it's reasonable to keep N(t) constant, but the birth and death rates may not be reasonable.

Additionally, we assume that the birth rate is proportional to the entire population, rather than restricting our attention to individuals of reproductive age. We also assume that the disease itself does not cause a loss of life expectancy.

**Click [here](03_DiscreteDeterministic_practical.qmd) to return to the practical.**