Noterbook 3 - w1: Solution: SIR model with a constant force of infection

In the last 2 etivities, you have learnt how you can represent differential equations in R through the use of the deSolve package. In this etivity you will use this skill to start coding a simplified SIR model within R. Let’s recap first what an SIR model is:

The differential equations for this system are:

\[\begin{align} \frac{dS}{dt} & = -\lambda S \\ \frac{dI}{dt} & = \lambda S - \gamma I \\ \frac{dR}{dt} & = \gamma I \end{align}\]

We can use the SIR model to describe a disease that can be split into 3 states: susceptible (\(S\)), infected (\(I\)), or recovered (R). All those infected are infectious, and all those recovered are immune, so they cannot get the disease again.

One part of this model should already be very familiar to you! The transition from I to R is what we explored in etivity 1. The new addition is the susceptibles in compartment S. As you learnt in the lecture, depending on how many people in the population are infectious, susceptible people experience a force of infection \(\lambda\) (lambda), which is the transition rate at which they become infected.

In this etivity, we are trying to simulate an outbreak of a new infectious disease that our population of 10\(^{6}\) people has not been exposed to before. This means that we are starting with a single case, everyone else is susceptible to the disease, and no one is yet immune or recovered. This can for example reflect a situation where an infected person introduces a new disease into a geographically isolated population, like on an island, or even when an infections “spill over” from other animals into a human population. In terms of the initial conditions for our model, we can define: S = 10\(^{6}\)-1 = 999999, I = 1 and R = 0.

Now it’s your turn to write and solve this model in R, assuming a constant force of infection of 0.2 days\(^{-1}\). Like in the first coding etivity, the infection has an average duration of 10 days, but this time we want to run the model for 60 days. Try to write the whole code by yourself - but if you have trouble remembering how, just go back to the etivity 1 solution to check. We have also added some hints in the cell below, to help you remember which different parts are required. As you go along, it is good practice to add your own comments to every part, explaining what the code is for. When you’re done writing the code and the comments, run it and print the output to double-check that everything works.

1 Solution: SIR model with a constant force of infection

In this etivity, you had the task of writing an SIR model in R from scratch. Compare your code to the solution given below, but note that yours might look slightly different in terms of the names you have chosen for your variables. That’s fine, as long as the code is easy to follow and produces the correct output!

Remember, the differential equations for the simple SIR model with a constant force of infection look like this:

\[\begin{align} \frac{dS}{dt} & = -\lambda S \\ \frac{dI}{dt} & = \lambda S - \gamma I \\ \frac{dR}{dt} & = \gamma I \end{align}\]

1.1 The input data from the instructions were as follows:

Initial number of people in each compartment:
S = 10\(^6\)-1, I = 1 and R = 0
Parameters:
\(\lambda\) = 0.2 days\(^{-1}\) (this represents a force of infection that’s constant at 0.2)
\(\gamma\) = 0.1 days\(^{-1}\) (corresponding to an average duration of infection of 10 days)
We want to run the model for 60 days.

# LOAD THE PACKAGES:
library(deSolve)
library(reshape2)
library(ggplot2)

# MODEL INPUTS:

# Vector storing the initial number of people in each compartment (at timestep 0)

initial_state_values <- c(S = 999999,  # the whole population we are modelling is susceptible to infection
                          I = 1,       # the epidemic starts with a single infected person
                          R = 0)       # there is no prior immunity in the population

# Vector storing the parameters describing the transition rates in units of days^-1
parameters <- c(lambda = 0.2,  # the force of infection, which acts on susceptibles
                gamma = 0.1)   # the rate of recovery, which acts on those infected

# TIMESTEPS:

# Vector storing the sequence of timesteps to solve the model at
times <- seq(from = 0, to = 60, by = 1)   # from 0 to 60 days in daily intervals

# SIR MODEL FUNCTION: 

# The model function takes as input arguments (in the following order): time, state and parameters
sir_model <- function(time, state, parameters) {  

    with(as.list(c(state, parameters)), {  # tell R to unpack variable names from the state and parameters inputs
        
    # The differential equations
      dS <- -lambda * S               # people move out of (-) the S compartment at a rate lambda (force of infection)
      dI <- lambda * S - gamma * I    # people move into (+) the I compartment from S at a rate lambda, 
                                      # and move out of (-) the I compartment at a rate gamma (recovery)
      dR <- gamma * I                 # people move into (+) the R compartment from I at a rate gamma
      
    # Return the number of people in the S, I and R compartments at each timestep 
    # (in the same order as the input state variables)
    return(list(c(dS, dI, dR))) 
    })
  
}

# MODEL OUTPUT (solving the differential equations):

# Solving the differential equations using the ode integration algorithm
output <- as.data.frame(ode(y = initial_state_values, 
                            times = times, 
                            func = sir_model,
                            parms = parameters))
# Printing the model output returns a dataframe with columns time (containing the times vector), 
# S (containing the number of susceptible people at each timestep),
# I (containing the number of infected people at each timestep) and 
# R (containing the number of recovered people at each timestep).
output

##    time            S          I          R
## 1     0 9.999990e+05      1.000      0.000
## 2     1 8.187299e+05 172214.061   9056.004
## 3     2 6.703194e+05 296821.937  32858.688
## 4     3 5.488111e+05 384013.528  67175.386
## 5     4 4.493285e+05 441982.410 108689.084
## 6     5 3.678791e+05 477302.608 154818.340
## 7     6 3.011939e+05 495234.957 203571.160
## 8     7 2.465967e+05 499976.743 253426.573
## 9     8 2.018963e+05 494864.916 303238.803
## 10    9 1.652987e+05 482541.541 352159.774
## 11   10 1.353351e+05 465088.345 399576.571
## 12   11 1.108030e+05 444135.873 445061.148
## 13   12 9.071779e+04 420952.533 488329.674
## 14   13 7.427343e+04 396516.440 529210.125
## 15   14 6.080994e+04 371573.802 567616.261
## 16   15 4.978696e+04 346686.177 603526.866
## 17   16 4.076211e+04 322268.614 636969.278
## 18   17 3.337319e+04 298620.488 668006.326
## 19   18 2.732365e+04 275950.307 696726.044
## 20   19 2.237071e+04 254395.666 723233.625
## 21   20 1.831558e+04 234039.257 747645.159
## 22   21 1.499553e+04 214921.669 770082.801
## 23   22 1.227730e+04 197051.601 790671.099
## 24   23 1.005180e+04 180413.979 809534.219
## 25   24 8.229718e+03 164976.376 826793.906
## 26   25 6.737922e+03 150694.067 842568.011
## 27   26 5.516543e+03 137513.992 856969.465
## 28   27 4.516563e+03 125377.829 870105.609
## 29   28 3.697848e+03 114224.364 882077.788
## 30   29 3.027542e+03 103991.298 892981.160
## 31   30 2.478741e+03  94616.601 902904.658
## 32   31 2.029421e+03  86039.514 911931.065
## 33   32 1.661549e+03  78201.265 920137.185
## 34   33 1.360361e+03  71045.572 927594.066
## 35   34 1.113770e+03  64518.965 934367.266
## 36   35 9.118772e+02  58570.980 940517.143
## 37   36 7.465818e+02  53154.252 946099.167
## 38   37 6.112494e+02  48224.527 951164.223
## 39   38 5.004486e+02  43740.622 955758.929
## 40   39 4.097326e+02  39664.335 959925.932
## 41   40 3.354606e+02  35960.336 963704.203
## 42   41 2.746519e+02  32596.029 967129.319
## 43   42 2.248659e+02  29541.405 970233.729
## 44   43 1.841046e+02  26768.894 973047.002
## 45   44 1.507321e+02  24253.202 975596.066
## 46   45 1.234090e+02  21971.163 977905.428
## 47   46 1.010387e+02  19901.583 979997.378
## 48   47 8.272349e+01  18025.097 981892.179
## 49   48 6.772825e+01  16324.028 983608.243
## 50   49 5.545120e+01  14782.255 985162.293
## 51   50 4.539959e+01  13385.087 986569.513
## 52   51 3.717004e+01  12119.146 987843.684
## 53   52 3.043225e+01  10972.258 988997.310
## 54   53 2.491581e+01   9933.350 990041.734
## 55   54 2.039934e+01   8992.358 990987.243
## 56   55 1.670156e+01   8140.135 991843.163
## 57   56 1.367408e+01   7368.375 992617.951
## 58   57 1.119539e+01   6669.536 993319.268
## 59   58 9.166009e+00   6036.774 993954.060
## 60   59 7.504492e+00   5463.877 994528.618
## 61   60 6.144158e+00   4945.213 995048.643

#Plotting the output:

output_long <- melt(as.data.frame(output), id = "time")                  # turn output dataset into long format

ggplot(data = output_long,                                               # specify object containing data to plot
       aes(x = time, y = value, colour = variable, group = variable)) +  # assign columns to axes and groups
  geom_line() +                                                          # represent data as lines
  xlab("Time (days)")+                                                   # add label for x axis
  ylab("Number of people") +                                             # add label for y axis
  labs(colour = "Compartment")                                           # add legend title

In the following cell, plot the number of people in each compartment over time. You should see that at the peak of the epidemic, around 500000 people are infected

1.1.1 Based on the plot, describe the pattern of the epidemic over the 2 month period. How does the number of people in the susceptible, infected and recovered compartment change over time? After how many days does the epidemic reach its peak? After how many days does it end?

The number of infected people quickly increases, reaching a peak of 500000 infected people after around 7 days, before steadily decreasing again. The number of recovered people starts to rise shortly after the first people become infected. It increases steadily (but less sharply than the curve of infected people) until the whole population has become immune - by day 53, 99% are in the R compartment, and nearly no susceptible people remain after 60 days.