SIR 2 model with a dynamic force of infection -restudy

library(deSolve)
library(reshape2)
library(ggplot2)

1 MODEL INPUTS:

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

initial_state_values <- c(S = 1000000-1,  # 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(beta = 0.5,     # the infection rate, which acts on susceptibles
                gamma = 0.25)   # the rate of recovery, which acts on those infected

1.2 TIMESTEPS:

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

1.3 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    
        
    # Calculating the total population size N (the sum of the number of people in each compartment)
      N <- S+I+R
      
    # Defining lambda as a function of beta and I:
      lambda <- beta * I/N
        
    # 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))

output

##     time        S            I            R
## 1      0 999999.0 1.000000e+00 0.000000e+00
## 2      1 999998.4 1.284025e+00 2.840260e-01
## 3      2 999997.7 1.648719e+00 6.487215e-01
## 4      3 999996.8 2.116996e+00 1.117000e+00
## 5      4 999995.6 2.718271e+00 1.718280e+00
## 6      5 999994.0 3.490321e+00 2.490338e+00
## 7      6 999992.0 4.481646e+00 3.481677e+00
## 8      7 999989.5 5.754526e+00 4.754581e+00
## 9      8 999986.2 7.388910e+00 6.389005e+00
## 10     9 999982.0 9.487474e+00 8.487635e+00
## 11    10 999976.6 1.218204e+01 1.118231e+01
## 12    11 999969.7 1.564184e+01 1.464230e+01
## 13    12 999960.8 2.008418e+01 1.908494e+01
## 14    13 999949.4 2.578802e+01 2.478930e+01
## 15    14 999934.8 3.311153e+01 3.211366e+01
## 16    15 999916.0 4.251448e+01 4.151801e+01
## 17    16 999891.8 5.458707e+01 5.359292e+01
## 18    17 999860.8 7.008689e+01 6.909657e+01
## 19    18 999821.0 8.998624e+01 8.900226e+01
## 20    19 999769.9 1.155329e+02 1.145594e+02
## 21    20 999704.3 1.483278e+02 1.473715e+02
## 22    21 999620.1 1.904247e+02 1.894969e+02
## 23    22 999512.0 2.444574e+02 2.435766e+02
## 24    23 999373.2 3.138026e+02 3.129991e+02
## 25    24 999195.1 4.027871e+02 4.021113e+02
## 26    25 998966.6 5.169526e+02 5.164870e+02
## 27    26 998673.3 6.633909e+02 6.632717e+02
## 28    27 998297.2 8.511693e+02 8.516208e+02
## 29    28 997814.9 1.091867e+03 1.093258e+03
## 30    29 997196.6 1.400247e+03 1.403184e+03
## 31    30 996404.3 1.795093e+03 1.800573e+03
## 32    31 995389.8 2.300245e+03 2.309905e+03
## 33    32 994091.8 2.945856e+03 2.962380e+03
## 34    33 992432.4 3.769896e+03 3.797676e+03
## 35    34 990314.0 4.819910e+03 4.866126e+03
## 36    35 987613.7 6.154997e+03 6.231349e+03
## 37    36 984178.7 7.847904e+03 7.973399e+03
## 38    37 979820.5 9.987048e+03 1.019244e+04
## 39    38 974309.1 1.267808e+04 1.301286e+04
## 40    39 967367.9 1.604440e+04 1.658771e+04
## 41    40 958671.4 2.022562e+04 2.110294e+04
## 42    41 947846.2 2.537277e+04 2.678100e+04
## 43    42 934479.0 3.163844e+04 3.388257e+04
## 44    43 918134.9 3.916014e+04 4.270500e+04
## 45    44 898388.8 4.803551e+04 5.357565e+04
## 46    45 874872.7 5.828936e+04 6.683794e+04
## 47    46 847336.1 6.983551e+04 8.282844e+04
## 48    47 815715.3 8.244026e+04 1.018444e+05
## 49    48 780197.5 9.569890e+04 1.241036e+05
## 50    49 741260.7 1.090384e+05 1.497010e+05
## 51    50 699672.3 1.217566e+05 1.785711e+05
## 52    51 656436.1 1.330994e+05 2.104645e+05
## 53    52 612690.8 1.423622e+05 2.449470e+05
## 54    53 569582.4 1.489921e+05 2.814255e+05
## 55    54 528140.7 1.526635e+05 3.191958e+05
## 56    55 489187.4 1.533083e+05 3.575043e+05
## 57    56 453289.5 1.510989e+05 3.956116e+05
## 58    57 420760.2 1.463943e+05 4.328456e+05
## 59    58 391691.2 1.396686e+05 4.686402e+05
## 60    59 366004.1 1.314411e+05 5.025548e+05
## 61    60 343503.9 1.222182e+05 5.342779e+05
## 62    61 323926.4 1.124546e+05 5.636190e+05
## 63    62 306976.1 1.025317e+05 5.904922e+05
## 64    63 292352.6 9.275052e+04 6.148969e+05
## 65    64 279767.3 8.333463e+04 6.368981e+05
## 66    65 268953.6 7.443870e+04 6.566077e+05
## 67    66 259671.3 6.615981e+04 6.741689e+05
## 68    67 251707.5 5.854918e+04 6.897434e+05
## 69    68 244876.2 5.162311e+04 7.035006e+05
## 70    69 239016.3 4.537245e+04 7.156112e+05
## 71    70 233988.6 3.977046e+04 7.262410e+05
## 72    71 229673.7 3.477896e+04 7.355474e+05
## 73    72 225969.3 3.035314e+04 7.436776e+05
## 74    73 222787.9 2.644505e+04 7.507671e+05
## 75    74 220054.6 2.300617e+04 7.569392e+05
## 76    75 217705.6 1.998912e+04 7.623053e+05
## 77    76 215686.1 1.734880e+04 7.669652e+05
## 78    77 213949.3 1.504313e+04 7.710076e+05
## 79    78 212455.3 1.303337e+04 7.745114e+05
## 80    79 211169.7 1.128429e+04 7.775460e+05
## 81    80 210063.3 9.764100e+03 7.801726e+05
## 82    81 209110.9 8.444370e+03 7.824447e+05
## 83    82 208290.9 7.299788e+03 7.844093e+05
## 84    83 207584.7 6.307944e+03 7.861073e+05
## 85    84 206976.5 5.449077e+03 7.875744e+05
## 86    85 206452.7 4.705822e+03 7.888415e+05
## 87    86 206001.4 4.062957e+03 7.899357e+05
## 88    87 205612.6 3.507179e+03 7.908802e+05
## 89    88 205277.6 3.026881e+03 7.916955e+05
## 90    89 204988.9 2.611951e+03 7.923991e+05
## 91    90 204740.2 2.253598e+03 7.930062e+05
## 92    91 204525.8 1.944186e+03 7.935300e+05
## 93    92 204341.1 1.677088e+03 7.939818e+05
## 94    93 204181.8 1.446560e+03 7.943716e+05
## 95    94 204044.6 1.247628e+03 7.947078e+05
## 96    95 203926.3 1.075985e+03 7.949977e+05
## 97    96 203824.4 9.279044e+02 7.952477e+05
## 98    97 203736.5 8.001653e+02 7.954633e+05
## 99    98 203660.8 6.899832e+02 7.956493e+05
## 100   99 203595.5 5.949521e+02 7.958096e+05
## 101  100 203539.2 5.129940e+02 7.959478e+05

1.4 PLOTTING THE OUTPUT

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


# Adding a column for the proportion of the population in each compartment at each timestep
# One way of calculating this is dividing the number in each compartment by the total initial population size
# We can do this in this case because our population is closed, so the population size stays the same
# at every timestep
output_long$proportion <- output_long$value/sum(initial_state_values)

1.4.1 Plot this new column

ggplot(data = output_long,                                               # specify object containing data to plot
       aes(x = time, y = proportion, 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("Proportion of the population") +                                 # add label for y axis
  labs(colour = "Compartment")                                           # add legend title

2 What do you observe when β = 0.5 and γ = 0.25?

An epidemic occurs, reaching a peak 56 days after introduction of the first infectious case, at which point about 15% of the population are infected. By the end of the epidemic, about 80% of the population have been infected and recovered.

2.1 Modelling a scenario where beta drops to 0.1 because an infection control measure is introduced:

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

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

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

# Calculating the proportion in each compartment
output_long$proportion <- output_long$value/sum(initial_state_values)

ggplot(data = output_long,                                               # specify object containing data to plot
       aes(x = time, y = proportion, 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("Proportion of the population") +                                 # add label for y axis
  labs(colour = "Compartment")                                           # add legend title