Question 1

  1. Differential equations for Model 3

“\(\frac{dS}{dt} = \Lambda - \beta*\frac{I}{N}*S -\mu*S\)” “\(\frac{dE}{dt} = \beta*\frac{I}{N}*S - \epsilon*E - \mu*E\)” “\(\frac{dI}{dt} = \epsilon*E - \gamma*I -\mu*I -\alpha*I\)” “\(\frac{dR}{dt} = \gamma*I -\mu*R\)”

Question 2

• Births occur at a constant rate and all individuals are susceptible • Natural deaths occur at a constant rate and natural deaths affect all compartments • Homogenous population: individuals in the population are believed to have the same characteristics (Mishra et al., 2011) • The population mixes homogeneously, all individuals have an equal probability of contact • I’m adding this from the Density-dependent vs. Frequency-dependent Disease Transmission blog post but not sure if it’s needed - we assume that transmission in frequency-dependent and therefore proportional to the proportion of infectious individuals in the population • Recovered individuals become completely and permanently immune, therefore there are no reinfections • There is no vaccine or other intervention methods in this model

Question 3

Beta = 0.34 Recovery rate = 1/7 = 0.1429 mu = 1 / (365 * 82) = 0.00003 per day

R0 = Beta/Recovery Rate + Per Capita Birth Rate R0 = (0.34/ 0.143 + 0.00003) R0 = 2.38

Therefore, the basic reproductive number for ‘wild-type’ SARS-CoV-2 in a Canadian population is 2.38.

Question 4

Using the example code provided in Part 2, code this simple SARS-CoV-2 model in R. There are a couple of things that are slightly different than in the example so you will need to pay close attention to the model time, and initial conditions since Table 1 tells you that the model timestep here is days. Run your model for 3 months (90 days).

require(deSolve)

 seirmod=function(t, y, parms){
   #Pulling state variables from y vector
   S=y[1]
   E=y[2]
   I=y[3]
   R=y[4]
   #Pulling parameter values from parms vector
   beta <- parms["beta"]
   eps <- parms["eps"]
   mu <- parms["mu"]
   gamma <- parms["gamma"]
   alpha <- parms["alpha"]
   Lambda <- parms["Lambda"]
   N <- parms["N"]
   #Defining the equations
   dS = Lambda - beta * S * I / N - mu * S
   dE = beta * S * I / N - eps * E - mu * E
   dI = eps * E - gamma * I - mu * I - alpha * I
   dR = gamma * I - mu * R
   res=c(dS, dE, dI, dR)
   #Return list of gradients
   list(res)
 }

Setting the new simulation times below:

times  = seq(0, 90, by=1/10) #I think this is 10 simulation run-throughs per day?
parms  = c(beta = 0.34, eps = 1/2.3, gamma = 1/7, mu = 0.0003, alpha = 0.006, Lambda = 0.0003, N = 1)

INitial conditions:

start = c(S=0.999, E = 0, I=0.001, R = 0)

Running the model:

out = ode(y = start, times = times, func = seirmod, 
     parms = parms)
out=as.data.frame(out)
head(round(out, 3))
##   time     S E     I R
## 1  0.0 0.999 0 0.001 0
## 2  0.1 0.999 0 0.001 0
## 3  0.2 0.999 0 0.001 0
## 4  0.3 0.999 0 0.001 0
## 5  0.4 0.999 0 0.001 0
## 6  0.5 0.999 0 0.001 0

Plotting the initial results:

plot(x = out$time, y = out$S, ylab = "Fraction", 
      xlab = "Time", type = "l")
 lines(x = out$time, y = out$I, col = "red")
 lines(x = out$time, y = out$R, col = "green")

#note that this figure is generated using simple plotting and not ggplot.

Question 5

Plotting an updated figure:

Question 6

Question 7