Reproductive Number 2 - Practice

In DSAIDE: Dynamical Systems Approach to Infectious Disease Epidemiology (Ecology/Evolution)

#load various variable definitions that are the same for each app
source('startup_script.R')
currentrmdfile = knitr::current_input() 
appsettings = get_settings(currentrmdfile,appdocdir,packagename)

Overview {#shinytab1}

This app covers more aspects of the reproductive number concept. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab. If you are not familiar with the reproductive number, you should go through the 'Reproductive Number 1' app first.

This app also assumes that you have worked through the apps in the The Basics section.

Learning Objectives

• Know how to estimate the reproductive number from steady state levels of disease.
• Understand how the reproductive number can be used to parameterize models.

The Model {#shinytab2}

Model Overview

For this app, we'll use a basic compartmental SIR model that also includes births, deaths and waning immunity. We allow for 3 different stages/compartments:

• S - uninfected and susceptible individuals
• I - infected and infectious individuals (note that these terms are often used interchangeably, but technically we are talking about someone who is infected and is infectious, i.e. can infect others)
• R - recovered/removed individuals.

In addition to specifying the compartments of a model, we need to specify the dynamics determining the changes for each compartment. Broadly speaking, some processes increase the number of individuals in a given compartment/stage and other processes that lead to a reduction. Those processes are sometimes called inflows and outflows.

For our system, we specify the following processes/flows:

• Susceptible individuals (S) can become infected by infectious individuals (I) at some rate, b. This leads to the susceptible individual leaving the S compartment and entering the I compartment.
• Infected individuals recover and enter the recovered (R) compartment at some rate, g.
• Recovered individuals can lose immunity at rate w and return to the susceptible compartment.
• Natural births (into the S compartment at rate n) and deaths (from all compartments at rate m) are possible.

Model Implementation

The flow diagram and the set of equations which are used to implement this model are as follows:

knitr::include_graphics(here::here('inst/media',appsettings$modelfigname))

$$\dot S =n - b SI - mS + wR$$ $$\dot I = b S I - g I - mI$$ $$\dot R = g I - mR - wR$$

Reproductive number

The app and tasks deal with the reproductive number concept. I recommend reading a bit more about it. I'm following the terminology of my own write-up. You can also check the books listed in the ID introduction app or some of the papers listed in the _Further Resources section of this app._

In the Reproductive Number 1 app, we looked at how one can write down and equation for the reproductive number for the model used there. The model we use here is similar but not quite the same. We can use the basic definition to easily figure out the reproductive number. Recall that it is defined as the average number of new infected (and infectious) individuals caused by one infectious individual.

For the model shown above, one can figure the value out by determining how many new infections are caused by one infected person. A person is infectious for a duration of 1/(g+m), the inverse of all the outflows out of the compartment. During that time they infect others at rate b. Thus the average number of new infections during created in b/(g+m). For the whole population, assuming initially everyone is susceptible, we multiply by the number of initial susceptibles to get $$R_0=\frac{bS_0}{g+m}$$

What to do {#shinytab3}

The tasks below are described in a way that assumes everything is in units of MONTHS (rate parameters, therefore, have units of inverse months). If any quantity is not given in those units, you need to convert it first (e.g. if it says a year, you need to convert it to 12 months).

#this is the running counter for the records which starts at 1 
rc=1

#empty object, will hold all outcomes
alloutcomes = NULL


# Task 1
tid = 1
tasktext = "In _Reproductive Number 1_, the focus was on a single outbreak. But the reproductive number concept is also useful in situations where a disease is endemic. As you learned in the _Patterns of ID_ app (if not, go do it before this one), to achieve an endemic state, we need replenishment of susceptibles through for instance births or waning immunity. We'll start with natural births and deaths. If births (and deaths) are present, you will be able to get multiple outbreaks and endemic states. Set the parameter values such that your hosts have approximately an average lifespan of 41 years (remember to convert to months before taking the inverse). Round to the first significant digit (i.e. 0.00X). Set birth rate such that the population is steady at 1000 in the absence of any infected hosts (if you need a reminder how that's done, revisit the _ID Patterns_ app). 

Now we'll introduce the disease. Set simulation time to around 500 years, 1 infected, _g_=5, _b_=0.015. Run the simulation, make sure you reach a steady state. Compare the steady state values for S (and optional I and R) from the simulation with those predicted from the SIR steady state equations (see _Patterns of ID_ app)."

# Record for task 1
nrec = 2 # number of items to record
out_records = c("The value you used for birth rate",
                    "Number susceptible at steady state")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 

##########################
# Task 2
##########################
tid = tid + 1 #increase task counter
tasktext = "Figure out what the value for the reproductive number is at the endemic steady state. To that end, recall its definition (average number of new infectious persons created by one infectious person) and determine what has to be true if the number of infected neither increases nor decreases. How does that compare to the reproductive number at the peak of an outbreak, which you learned about before?"
nrec = 1 # number of items to record
out_records = c("Value of R at endemic steady state")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


##########################
# Task 3
##########################
tid = tid + 1
tasktext = "At steady state, you have the _R_ value you just figured out in the previous task and a fraction of the population is susceptible, namely  _S~f~_ / _S~0~_ where _S~f~_ is the number of susceptibles at steady state and _S~0~_ is the number of susceptibles at the beginning. You learned that the equation for the reproductive number for this model is _R~0~=bS~0~/(g+m)_. Use it to compute the reproductive number. Compare the value of _R~0~_ at the start (where the fraction of susceptible is 1) with the value of the fraction of susceptible and the R value at steady state. Can you figure out how they relate? It should confirm what you figured out about _R_ at steady state in the previous task."
nrec = 2 # number of items to record
out_records = c("Number susceptible left at steady state",
                "The value for R~0~")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


##########################
# Task 4
##########################
tid = tid + 1
tasktext = "Next double the values for both natural birth and death rates. What do you get for the value of S at steady state? Is that surprising? If it is, take another look at the relation between _R~0~_ and S at steady state, and the equation for _R~0~_."
nrec = 1 # number of items to record
out_records = c("The value of S at steady state (end of simulation)")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 



##########################
# Task 5
##########################
tid = tid + 1
tasktext = "A) Turn off births and deaths and instead assume immunity wanes over 10 months (convert to a rate). Run the simulation with the same transmission rate and recovery rate as task 1. Record the value of S at steady state. How does it compare to the value from the previous task? Is that surprising? Again, look at the equation for _R~0~_. 
B) Now double the initial number of susceptible. What is _R~0~_ now? What do you therefore expect for the number of susceptibles at steady state? Run the simulation to see what you get. If the numbers confuse you at first glance, think carefully about the relation between _R~0~_ and both number and fraction susceptible left at steady state. Also, compare what you find here with the way to get the steady state as discussed in _Patterns of ID_ and revisited in task 1 above." 

nrec = 2 # number of items to record
out_records =  c("Number susceptible at steady state (A)",
                    "Number susceptible at steady state (B)")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


##########################
# Task 6
##########################
tid = tid + 1
tasktext = "One usually obtains R~0~ from the data/literature and uses it to pick values for the transmission parameter, _b_, which is otherwise very hard to estimate. Let's try that. Go online and find (approximate) values for the duration of the infectious period and _R~0~_ for SARS. Start with 1000 susceptible, 1 infected. Assume that for a single outbreak we can ignore births, deaths or waning immunity (all zer0). Use those values and the _R~0~_ equation to compute the transmission parameter, _b_. Then run a simulation with that value for _b_ (for 12 months). Recall from the _Reproductive Number 1_ app how you can use the fraction if suceptibles left at the end of an outbreak to compute _R~0~_. Do that here again to confirm that you get indeed (approximately) the value you started out with."
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


##########################
# Task 7
##########################
tid = tid + 1
tasktext = "Let's repeat the same for a disease that causes an endemic state, tuberculosis (TB). Go online and find estimates for the basic reproductive number of TB. Note that TB is somewhat special: An infected and infectious person (someone with TB disease) might infect many others, but only some of those will go on to develop disease and become infectious. For R~0~, we always need to go from infectious individual to another infectious individual. Also find an estimate for the duration of the infectious period for TB. Then use these quantities to determine parameters _b_ and _g_ and run a simulation with 1000 susceptible and the same birth and death rates as in task 1 (m=0.002 and n=2). Check that the steady state levels of infected/infectious individuals matches with the R~0~ value you chose.

Note: The SIR model is not a good model for TB since for TB, the stage where individuals are infected but not yet infectious is long and important. So to really model TB, one would need to include such details. We are only using TB here for illustrative purposes."

# Record for task 5
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 

#save the fully filled task table to a tsv file
alloutcomes$QuizID = paste0(packagename,"_",appsettings$appid)
alloutcomes$AppTitle = appsettings$apptitle
alloutcomes$AppID = appsettings$appid
#remove a few variables from the data frame
savedoutcomes <- dplyr::select(alloutcomes,QuizID,AppID,AppTitle,TaskID,TaskText,RecordID,Record,Type,Note)     
write.table(savedoutcomes, paste0(appsettings$appid,"_tasktable.tsv"), append = FALSE, sep = "\t", row.names = F, col.names = TRUE)

# Take all the text stored in the table and print the tasks and items to record
write_tasktext(alloutcomes)

Further Information {#shinytab4}

• This app (and all others) are structured such that the Shiny part (the graphical interface you see and the server-side function that goes with it) calls an underlying R script (or several) which runs the simulation for the model of interest and returns the results.
• For this app, the underlying function(s) running the simulation are called r appsettings$simfunction. You can call them directly, without going through the shiny app. Use the help() command for more information on how to use the functions directly. If you go that route, you need to use the results returned from this function and produce useful output (such as a plot) yourself.
• You can also download all simulator functions and modify them for your own purposes. Of course to modify these functions, you'll need to do some coding.
• For examples on using the simulators directly and how to modify them, read the package vignette by typing vignette('DSAIDE') into the R console.
• For the equations mentioned in the tasks, see e.g. [@heffernan2005; @keeling08; @vynnycky10].
• For some information on the reproductive number, see e.g. [@fine2011; @heffernan2005; @roberts2007; @wallinga2007].

References

DSAIDE documentation built on Aug. 24, 2023, 1:07 a.m.