EPI 569 – Homework 2: Modeling a directly-transmitted pathogen

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

Background 1

• A new pathogen has emerged in the Republic of Zubrowka. Here’s what we know so far:
  • All individuals are susceptible` and become infected by contact with an infected individual
  • People are infectious for a while and then become immune
  • Immunity eventually decays.
• We can partition our population of size N into susceptible (S), infected (I), and immune (R) people, the above-mentioned state variables, so N = S + I + R.

Background 2

• Then, consistent with the lecture notation, we write these as difference equations:
• In these equations,
  • \(\lambda\) is the force (or rate) of infection per susceptible person. Which is:
  • contact rate, \(\alpha\), multipled by
  • probability of infection on contact with an infectious person \(\beta\), multipled by
  • the probability that a randomly encountered person is infectious, \(\frac{I_{t-1}}{N}\).

Background 3

• Other symbols:
  • \(\mu\) is the birth and death rate, maintaining constant population size;
  • \(\omega\) is the rate at which immunity is lost;
  • \(\sigma\) is the recovery rate and
  • \(\delta\) t is the interval between successive calculations or time-step.

Question 1 (1pt)

• Draw a flow diagram of this transmission process

• ANSWER

• QUESTION 1 & 2 ANSWERS ARE IN THIS PIC: ## Question 2 (0.5pt)

• Rewrite the difference equations above as a system of differential equations. You may use the Equation editor in Word, markdown in R or another program of your choosing. You can also hand write and insert a picture.

• Go through these equations to ensure that you understand what each term implies about the epidemiology of infection. Review the formulae in the spreadsheet. Make sure that they match the difference equations above.

• ANSWER

Background 4

• There are two verions of the model. One is an Excel spreadsheet that uses difference equations. In that version, you can easily see what’s going on the model. Second, there is an R version of the model (below) that uses differential equations. Add the parameter values to the spreadsheet or R script. Use values that seem sensible to you. Remember to keep all in the same units (days). For now, set \(\omega\) to 0 (i.e. lifelong immunity) and \(\mu\) to 0 (i.e. no births/ deaths/ migrations).

• Hint: In models with constant per capita rates, the time spent in a compartment is exponentially distributed, so the mean ‘residence times’ are reciprocals of the exit rates or vice versa. For example, \(\sigma\) is the reciprocal of the duration of clinical illness.

Background 5 – Model parameters

• You’ll use the parameter values below for questions 4 and 6

Parameter	Symbol	Value
Daily contact rate	\(\alpha\)	5.0
Transmision probability	\(\beta\)	0.1
Recovery rate	\(\sigma\)	0.2
Loss of immunty rate	\(\omega\)	0
Birth and death rate	\(\mu\)	0

Question 3 (0.5pt)

• In completely susceptible populations, another way to express the average number of infections per infectious person or basic reproduction number is \(R_{0}=\frac{\alpha \beta}{\mu +\sigma}\)
• Can you explain this in words?
• Hint: What are \(\alpha \beta\) and \(\frac{1}{\mu +\sigma}\), and their product?
• ANSWER R0 {R naught} is a measure of how contagious a disease is with nothing to stop it. The equation lays out how fast it spreads (the contact rate multiplied by the probability of infection) divided by, or over, the changes in the susceptible population which is caluclated by adding the balance of the birth and death rate to the recovery rate.

Question 4 (0.5pt)

• Calculate the basic reproductive number for your model.
• ANSWER R0 = (5*0.1)/(0+0.2) = 2.5

Question 5 (0.5pt)

• If the initial conditions are \(S_{0}\) = 100, \(I_{0}\) = 0, and \(R_{0}\) = 0, what does this say about the composition of the population at t = 0? * ANSWER This means the entire population is entirely susceptible and there are no infected or recovered people. This is a population that has never ancountered a novel pathogen. Additionally, the population has never been infected byand has not recovered from a disease so immunity or length thereof is unknown. R0 is zero because we do not know the parameters by which we would calculate this measure.

Question 6 (1pt)

• Now, introduce an infectious person by changing \(I_{0}\) = 1 and the values from the parameter table abive into the model spreadsheet. Calculate the number of individuals in S, I, and R compartments from t=1 to t=100.

• When does the epidemic peak? When does it end?

• Does everyone get eventually infected? If not, why

  • Hint: Consider the effective reproductive number \(R_{e}\)?

Question 6 (1pt)

• ANSWER
• Based on the numbers in the excel spread sheet, it appears that the epidemic peaks right around 18 days. It looks like it ends at 78 days.
• In a real population, not everyone gets infected. There will be some with natural immunity and some who had a previous infection that conferred immunity to this novel pathogen. The Re (R effective) is a measure of how this natural immunity impacts the spread of disease. Also, herd immunity will take over and a remaining percent of the population will be “protected” from infection.

# Define parameters -- THESE ARE THE MODEL PARAMETERS
parms <- c(alpha =5,        # alpha = daily contacts
           beta=0.1,        # beta = probability of infection on contact
           sigma=0.2,       # sigma = rate of recovery per day
           mu = 0.0,        # mu =  per capita birth and death rate
           omega = 0.0)     # omega = rate of immune loss per day

# Initial conditions --  THESE ARE THE CONDITIONS AT THE START OF THE SIMULATION
init <- c(S=99,           # number initilly susceptible
          I=1,            # number initally infectious
          R=0)            # initially immune or "recovered"

# Define model equations -- do not change -- or change with care!
# These are the model equations.  They are written as a function called sir_ode.  
# "parms" and "init" input the parameters and inital conditions into the equations 
sir_ode <- function(times,init,parms){
  with(as.list(c(parms,init)), {
    # ODEs
    dS <- mu*(S+I+R) + omega*R -alpha*beta*I*S/(S+I+R) - mu*S 
    dI <- alpha*beta*I*S/(S+I+R)-sigma*I - mu*I
    dR <- sigma*I  - omega*R - mu*R
    list(c(dS,dI,dR))
  })
}

# This creates the output from model equations.  
#If you want to run the model for longer, change the second term eg: seq(0,200,...)
times <- seq(0,100,length.out=100)
sir_out <- lsoda(init,times,sir_ode,parms)
sir_out_long <- melt(as.data.frame(sir_out),"time")

#Plotting the model output
ggplot(sir_out_long,aes(x=time,y=value,colour=variable,group=variable))+
  geom_line(lwd=2)+             #Add line
  xlab("Time")+ylab("Number")   #Add labels

Question 7 (0.5pt)

• Now, let’s say that some resident ofZubrowka had returned from a neighboring country and had immunity from exposure there. Lower \(S_{0}\) to reflect this, but keep the total population size, N, the same.

• Is there a threshold at which an epidemic does not occur? If so, what is it and why?
  *Hint: Again, consider the effective reproductive number RE

• ANSWER

• we would need Re = 1 for an epidemic not to occur. Re = R0*n (where n is the susceptible portion of the population), therefore the percent of susceptible population needed to not have an epidemic n=R/Re = 1/2.5 = 40%. Thus S=40 would be needed, and we would need 59 people to have been exposed while visiting the neigboring country.

Question 8 (0.5pt)

• Set the parameters back to the original values and now add some simple demography. Assume that people are entering and leaving the population reflected in a balanced death and birth rate.
  

• Set the rate to 0.1 per day.

• Also try other values. What happens now? How does changing the birth/death rate affect \(R_{0}\)? How does it affect the equilibrium and why?

• ANSWER

• If birth occurs more than death, then the susceptible population will increase and an epidemic threshold may not be met. If death occurs more than birth (say from this infection), then eventually the population will decline.

• IF Beta/ mu > 1 is an outbreak, and Beta/mu<1 is steady state or not an outbreak. Assuming beta is steady, then the larger mu is, the smaller R0 is. The smaller mu is, the larger R0 is. All else being equal, an increase in population will reduce R0, changing the creating a different trajectory and equilibrium.

Question 9 (0.5pt)

• Now, we will also allow immunity to wane by changing \(\omega\) to a value greater than 0.

• Choose parameter values that keep \(R_{0}\) bigger than 1 (e.g. (\(\omega\) = 0.05/day). How does adding this extra process affect the dynamics?

• Hint: You may wish to increase the increment between successive calculations (i.e. to 5 days), \(\delta\)t, to increase the span of time charted in the spreadsheet version.

• ANSWER

• Changing omega to 0.5 changes the loss of immunity from no loss of immunity to some loss of immunity. This keeps the immunity level high because people are cycling between the immune population and the susceptible population and almost never entering the recovered population and the disease reaches equilibrium at 25 days and becomes endemic.

Question 10 (0.5pt)

• \(R_{e}\) is defined as the average number of secondary infections produced by an infectious person at any given time E. What is the value of \(R_{e}\) when the infection has attained the endemic equilibrium? How does the duration of immunity affect \(R_{0}\)? What does this mean? (0.5pt)

• ANSWER

• when an epidemic reaches equilibrium

• Re is equal to 1 when a pathogen reaches endemic equilibrium. That means each infected person infects one additional person. A reduced duration of immunity can possibly increase the Re because there are more susceptible people to be infected and no recovered people taken out of the equations. This means that people enter into the Susceptible bucket sooner and increase transmission pressure, increasing Re, possibly increasing the likelihood of an outbreak.

Extensions

• The model(s) you have used so far in this exercise, even though simple, capture quite well many observed patterns of infection in human host populations. Nevertheless, they might be extended to include a variety of complexities.
• Consider how you would model interventions in this model, how would you model:
  • social distancing?
  • vaccination?
    
  • demographic dynamics (i.e. when birth != death rate)?
  • an infection with a carrier state?

Extra credit (0.5pt)

• Which of these complexities can be incorporated into the model without structural changes to the model (i.e. changing parameters only and not the equations)
• ANSWER
• Social distancing - for social distancing you would change Beta, the contact rate. You would decrease the contact rate as social distancing is implemented. This will increase Lambda, the force of infection. No new variables.
• vaccination - these people would be removed from the number in the “S” bucket and moved over to the “R” bucket so no new variables would be needed.
• demographic dynamics (i.e. when birth != death rate)? - when the birth is not equal to the death rate, the value for mu changes, no new variable is added to the equation.
• an infection with a carrier state - Because the carrier state occurs post-infection, this would be a SIRC Model. additional variables include:
• rho = the proportion of acute infections that become carriers.
• Tau = length of time in carrier state.
• C = number of people in carrier state