part 2 of session 4, ODEs

03_DiscreteDeterministic_solutions.qmd

@@ -54,7 +54,12 @@ for (i in 2:nrow(y_sir))
 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:
+**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 |
 |-----------------------|----------------------------|----------------------|
@@ -90,10 +95,9 @@ 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")
+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")
-title(xlab = "Time (days)", ylab = "Infectious people")
 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")
 
@@ -103,7 +107,7 @@ When we start with $\beta=1.3, \gamma=0.23$ the peak occurs on day 7 with 60% of
 
 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 $\gamma=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.
+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).
 

04_ODEs_practical.qmd

@@ -112,4 +112,262 @@ In this model, once infected, susceptible individuals move to the exposed class.
 
 We will assume that susceptible individuals are vaccinated at a rate $v = 0.01$. The vaccine is 100% effective, so once vaccinated, individuals cannot become infected. HINT: you will need to create a new class $V$, and you can assume that the initial number of vaccinated individuals is 0.
 
+## Practical 2: Advanced use of `deSolve`
+
+In this practical, we'll be exploring some more advanced uses of the `deSolve` package.
+
+Copy and paste the following R Script into RStudio and follow along with the instructions in the comments:
+
+```{r, eval = FALSE}
+######################################################
+#              ODEs in R: Practical 2                #
+######################################################
+
+# Load in the deSolve package
+library(deSolve)
+# If the package is not installed, install using the install.packages() function
+
+########## PART I. TIME DEPENDENT TRANSMISSION RATE
+
+# The code below is for an SI model with a time dependent transmission rate.
+# The transmission rate is a function of the maximum value of beta (beta_max)
+# and the period of the cycle in days.
+
+# Define model function 
+SI_seasonal_model <- function(t, state, parms)
+{
+    # Get variables
+    S <- state["S"]
+    I <- state["I"]
+    N <- S + I
+    # Get parameters
+    beta_max <- parms["beta_max"]
+    period <- parms["period"]
+    # Calculate time-dependent transmission rate
+    beta <- beta_max / 2 * (1 + sin(2*pi*t / period))
+    # Define differential equations
+    dS <- -(beta * S * I) / N
+    dI <- (beta * S * I) / N
+    res <- list(c(dS, dI))
+    return(res)
+}
+
+# Define parameter values
+parameters <- c(beta_max = 0.4, period = 10)
+
+# Define time to solve equations
+times <- seq(from = 0, to = 50, by = 1)
+
+# 1. Plot the equation of beta against a vector of times to understand the time
+#    dependent pattern. 
+# How many days does it take to complete a full cycle? 
+# How does this relate to the period?
+
+##### YOUR CODE GOES HERE #####
+
+# 2. Now solve the model using the function ode and plot the output 
+# HINT: you need to write the state vector and use the ode function
+# using func = SI_seasonal_model
+
+##### YOUR CODE GOES HERE #####
+
+# 3. Change the period to 50 days and plot the output. What does the model output look like?
+# Can you explain why?
+
+
+
+########## PART II. USING EVENTS IN DESOLVE
+
+## deSolve can also be used to include 'events'. 'events' are triggered by 
+# some specified change in the system. 
+
+# For example, assume an SI model with births represents infection in a livestock popuation.
+# If more than half of the target herd size becomes infected, the infected animals are culled
+# at a daily rate of 0.5. 
+
+# As we are using an open population model (i.e. a model with births and/or deaths), 
+# we have two additional parameters governing births: the birth rate b, and the 
+# target herd size (or carrying capacity), K
+
+# Define model function 
+SI_open_model <- function(times, state, parms)
+{
+    ## Define variables
+    S <- state["S"]
+    I <- state["I"]
+    N <- S + I
+    # Extract parameters
+    beta <- parms["beta"]
+    K <- parms["K"]
+    b <- parms["b"]
+    # Define differential equations
+    dS <- b * N * (K - N) / K - ( beta * S * I) / N
+    dI <- (beta * S * I) / N 
+    res <- list(c(dS, dI))
+    return(res)
+}
+
+# Define time to solve equations
+times <- seq(from = 0, to = 100, by = 1)
+
+# Define initial conditions
+N <- 100
+I_0 <- 1
+S_0 <- N - I_0
+state <- c( S = S_0, I = I_0)
+
+# 4. Using beta = 0 (no infection risk), K = 100, b = 0.1, and an entirely susceptible 
+# population (I_0 = 0), investigate how the population grows with S_0 = 1, 50 and 100.
+
+# What size does the population grow to? Why is this?
+# Answer: 
+
+# How do you increase this threshold?
+# Answer:
+
+# Define parameter values
+# K is our target herd size, b is the birth rate)
+parameters <- c(beta = 0, K = 100, b = 0.1)
+
+# To include an event we need to specify two functions: 
+# the root function, and the event function.
+
+# The root function is used to trigger the event 
+root <- function(times, state, parms)
+{
+    # Get variables
+    S <- state["S"]
+    I <- state["I"]
+    N <- S + I
+    # Get parameters
+    K <- parms["K"]
+    # Our condition : more than half of the target herd size becomes infected
+    condition <- I <  K * 0.5 # This is a logical condition (TRUE/FALSE)
+    return (as.numeric(condition)) # Make this numeric, event occurs if root == 0
+}
+
+# The event function describes what happens if the event is triggered
+event_I_cull <- function(times, state, parms) 
+{
+    # Get variables
+    I <- state["I"]
+
+    # Extract parameters
+    tau <- parms["tau"]
+  
+    I <- I * (1 - tau) # Cull the infected population
+  
+    state["I"] <- I # Record new value of I
+  
+    return (state)
+}
+
+# We add the culling rate tau to our parameter vector 
+parameters <- c(beta = 0.1, K = 100, tau = 0.5, b = 0.1)
+
+# Solve equations
+output_raw <- ode(y = state, times = times, func = SI_open_model, parms = parameters,
+                  events = list(func = event_I_cull, root = TRUE), rootfun = root)
+# Convert to data frame for easy extraction of columns
+output <- as.data.frame(output_raw)
+
+# Plot output
+par(mfrow = c(1, 1))
+plot(output$time, output$S, type = "l", col = "blue", lwd = 2, ylim = c(0, N),
+     xlab = "Time", ylab = "Number")
+lines(output$time, output$I, lwd = 2, col = "red")
+legend("topright", legend = c("Susceptible", "Infected"),
+       lty = 1, col = c("blue", "red"), lwd = 2, bty = "n")
+
+# 5. What happens to the infection dynamics when the infected animals are culled?
+
+# Answer:
+
+
+# 6. Assume now that when an infected herd is culled, the same proportion of 
+# animals is ADDED to the susceptible population. 
+
+# HINT: You will need to change the event function to include additions to the 
+# S state.
+
+
+##### YOUR CODE GOES HERE #####
+
+
+# 7. What happens to the infection dynamics when the infected animals are culled?
+# How is this different to when only infected animals are culled?
+
+# Answer: 
+
+
+
+
+########## PART III. USING RCPP
+
+## Here we will code our differential equations using Rcpp to compare the speed
+# of solving the model.
+
+library(Rcpp)
+
+## The SIR Rcpp version can be compiled as follows:
+cppFunction('
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// This is a the SIR model coded in C++
+// Note that in Rcpp, you must declare all variables 
+
+// [[Rcpp::export]]
+List SIR_cpp_model(NumericVector t, NumericVector state, NumericVector parms)
+{
+    // Get variables
+    double S = state["S"];
+    double I = state["I"];
+    double R = state["R"];
+    double N = S + I + R;
+    
+    // Get parameters
+    double beta = parms["beta"];
+    double gamma = parms["gamma"];
+    
+    // Define differential equations
+    double dS = - (beta * S * I) / N;
+    double dI = (beta * S * I) / N - gamma * I;
+    double dR = gamma * I;
+    
+    NumericVector res_vec = NumericVector::create(dS, dI, dR);
+    
+    List res = List::create(res_vec);
+    
+    return res;
+}
+')
+
+# This can be solved in R using deSolve as follows
+
+# Define parameter values
+parameters <- c(beta = 0.4, gamma = 0.1)
+
+# Define initial conditions
+N <- 100
+I_0 <- 1
+S_0 <- N - I_0
+R_0 <- 0
+state <- c(S = S_0, I = I_0, R = R_0)
+
+# Solve equations
+output_raw <- ode(y = state, times = times, func = SIR_cpp_model, parms = parameters)
+# Convert to data frame for easy extraction of columns
+output <- as.data.frame(output_raw)
+
+# plot results
+par(mfrow = c(1, 1))
+plot(output$time, output$S, type = "l", col = "blue", lwd = 2, ylim = c(0, N),
+     xlab = "Time", ylab = "Number")
+lines(output$time, output$I, lwd = 2, col = "red")
+lines(output$time, output$R, lwd = 2, col = "green")
+legend( "topright", legend = c("Susceptible", "Infected", "Recovered"),
+        bg = rgb(1, 1, 1), lty = rep(1, 2), col = c("blue", "red", "green"), lwd = 2, bty = "n")
+```
+
 **Solutions to this practical can be accessed [here](04_ODEs_solutions.qmd).**