Multi Outbreak ID Control

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 explores how one should implement ID control measures for a scenario with multiple consecutive outbreaks. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab. Before going through this app, you should go through the Reproductive Number apps first. This app closely follows a model and analysis described in [@handel07a], see the "Further Information" tab for this and other related references.

Learning Objectives

• Understand the concept of susceptible threshold level and how it relates to population immunity and the reproductive number.
• Learn what 'overshoot' is and how it happens in infectious disease interventions.
• See how optimal control measures for multiple outbreaks can differ from those for a single outbreak.

The Model {#shinytab2}

Model Overview

For this app, we'll use the basic compartmental SIR model. We include 3 different stages/compartments:

• S - uninfected and susceptible individuals
• I - infected and infectious individuals
• R - recovered/removed individuals

For this app, 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.
• Control measures are applied at some level f (between 0 and 1) which reduces the rate at which new infections occur. Control measures start at time tstart and end at time tend. The parameter f can be non-zero outside of the treatment interval given by tstart and tend. What f does is to reduce the transmission rate by a certain fraction, e.g. f=0.5 means the rate of transmission is reduced to 0.5b. We do not specify the intervention here by which this reduction is achieved, it can be many types of interventions.
• To mimic the process of new infected individuals arriving in a specific location, at intervals tnew, the number of infected increases by 1.

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 = -(1-f) b SI$$ $$\dot I = (1-f)b S I$$ $$\dot R = g I$$

The simulation for this app also has one more (mostly invisible) feature. In previous apps, you learned that an unrealistic feature of models implemented by ordinary differential equations (ODE) is that the number of individuals can drop below 1, i.e., there could be a fraction of infected. For this app, the underlying code is such that if the number of infected drops below 1, it is set to 0. This is in some way a bit of a "hack" to deal with this issue. When you work through the stochastic apps, you will learn better ways of handling this. The advantage of using the "hack" is that we can keep using the ODE model formulation, without stochasticity, which makes things easy.

What to do {#shinytab3}

This app assumes knowledge of the reproductive number concept. If you are not familiar with it, please go through the 'Reproductive Number' apps first.

The tasks below are described in a way that assumes everything is in units of days (rate parameters, therefore, have units of inverse days). If any quantity is not given in those units, you need to convert it first.

#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 = "We'll start with an outbreak in the absence of control. Set the model to 1000 initial susceptible, 1 infected, no recovered. Set the recovery rate, _g_, such that it corresponds to an infectiousness duration of 5 days. Set the infectiousness rate _b_ such that the model has a reproductive number, R~0~, of 4. Keep the intervention level, _f_, at 0 for now. Set _tmax_ and _tnew_ to 100. Values for _tstart_ and _tend_ do not matter since _f_ is 0. Run the simulation. You should get an outbreak with around 20 susceptibles left at the end."
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 2
########################
tid = tid + 1
tasktext = "First without running the simulation and using what you know about R~0~ and how it relates to the susceptible population size, figure out what the number of susceptibles is (assuming every other parameter stays the same) at which you do not get an outbreak anymore. Now test your expectation by trying different values for _S_. You might have to adjust the simulation time to see this better. This value for _S_ provides a threshold below which we don't get an outbreak. This threshold value is often called the _herd immunity_ or _community/population immunity_ level. It corresponds to the threshold value of R~0~=1 which you learned about previously."
nrec = 1 # number of items to record
out_records = c("Number susceptible at which you stop getting an outbreak (population immunity level) based on theory (only use simulation to confirm your expectation)")
out_types = rep("Rounded_Integer",nrec)
out_notes = c("Report as integer")
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 = "In the previous task, you found the threshold/population immunity value for the number of susceptible needed to prevent an outbreak. In the 1st task, you also saw that at the end of the outbreak, there are much fewer susceptible left than the population immunity level. The outbreak _overshoots_ by depleting more susceptibles than the threshold value. Compute the additional number of susceptible that became infected during the outbreak compared to the value that would have prevented an outbreak.  Why is there such a difference? Think about what the value for the reproductive number is during the outbreak at the moment the number of susceptible has dropped to population level immunity. Does the outbreak stop immediately? Why not?"
nrec = 1 # number of items to record
out_records = c("Excess number by which susceptibles drop below herd immunity level")
out_types = rep("Rounded_Integer",nrec)
out_notes = c("Report as integer")
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 = "For a single outbreak, more control is better. We looked at at that in the _Basics of ID Control_ app, where we assumed some type of control lowered the transmission rate. We can repeat that here. Set parameters as in task 1. Set _tmax_ and _tnew_ to 200 (so we can see slow outbreaks if they happen). Set control to start at time 0 (_tstart_) and end at 200 (_tend_). That basically means the control is active for the whole duration of the simulation. Slowly increase control, _f_ (e.g. in steps of 0.05). Run the simulation for each value of _f_. As _f_ increases, you should see smaller outbreaks until you hit a value of _f_ for which you don't see outbreaks anymore. Understand how this value of _f_ relates to the reproductive number. You should be able to figure out the exact value for _f_ at which no outbreak happens anymore based on your knowledge of the reproductive number, and how a reduction in transmission (which is what f does) affects it."
nrec = 2 # number of items to record
out_records = c("Total/cumulative number of infected at end of simulation for f=0.5 (you need to compute this)",
                "Minimum value for f at which no outbreak is possible")
out_types = c("Rounded_Integer",
                "Numeric")
out_notes = c("Report as integer",
              "Report 2 decimals")
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 = "The idea that more control is better also applies if control is only applied during a certain period of the outbreak. Let's explore that by repeating the previous task, but now with control starting 10 days after the beginning of the outbreak and ending 120 days later (tend=130). Slowly increase control, in steps of 0.1, starting at 0. Run the simulation for each value of _f_. Since you don't start control at the beginning, you will always see an outbreak starting, but it will be reduced once control starts, with a higher reduction as _f_ increases."
nrec = 1 # number of items to record
out_records = c("Total/cumulative number of infected at end of simulation for f=0.5 (you need to compute this)")
out_types = rep("Rounded_Integer",nrec)
out_notes = c("Report as integer")
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 = "Doing the previous tasks, you will have noticed that timing of the outbreak matters, the earlier you start the better. Play around a bit with different start and end time for the control and different strengths to explore how they impact the size of the outbreaks you get and how much _overshoot_ you get, i.e. by how much the final number of susceptible is below the herd immunity level."
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 = "As you explored different start- and end-times and strengths for control in the previous task, you might have noticed that occasionally you can get the infected numbers to come back up once control is stopped. Let's explore a bit more how to apply optimal control if there is a chance of multiple outbreaks, either because the pathogen did not get completely wiped out in a specific community, or there are ongoing re-introduction of the pathogen from the outside. Set control to start at time 10 and end at 50. Assume no new infected persons enter the population by setting both total simulation time and _tmax_ and _tend_ to 200 days. All of the other settings as in the previous task. Slowly increase control, _f_ and run the simulation for the different values of _f_. Initially, things look as before. But once you reach a certain level of control, you will see that a second outbreak occurs and the total number of susceptibles drops again. To explore this in more detail, set f=0.6 and f=0.8 and try to understand why in one case you only get a single outbreak and in the other you get 2. To get the values for susceptible at t=50 it is easiest to switch to a plotly plot and read them off. Try to figure out how the number of susceptible when control is stopped relates to _R~0~_ and how that relates to the fact that you do or don't see another outbreak."
nrec = 4 # number of items to record
out_records = c( "Number of susceptible at end of control (t=50) for f=0.6",
                  "Number of susceptible at end of simulation for f=0.6",
                "Number of susceptible at end of control (t=50) for f=0.8",
                "Number of susceptible at end of simulation for f=0.8")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report as 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 8
########################
tid = tid + 1
tasktext = "Run the simulation with the previous settings for _f_ going from 0.5 to 0.7 in steps of 0.025. For each control level, record the number of susceptibles left at the end of the simulation. Based on that, what do you conclude about the impact of different levels of control on the outcome? What seems to be the best level of control, and why?"
nrec = 1 # number of items to record
out_records = c( "Value of _f_ among those you tried for which control is best")
out_types = rep("Numeric",nrec)
out_notes = rep("Report the exact value (three decimals)",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 9
########################
tid = tid + 1
tasktext = "In the previous task, control ends while there are still a few infected around, which can lead to a second outbreak. An alternative scenario is one where control ends after infected are gone, but then a newly infected person enters the population. We can explore this scenario, as well. Using parameter settings as previously, change control to start at time 10 (_tstart_) and end at 110 (_tend_). Set _tmax_ and _tnew_ to 400. Run the simulation for control strength _f=0.8_. You should see a single outbreak with around 475 susceptibles left. You learned above that this value of is not low enough for herd immunity, thus you can get a second outbreak. This is not happening here since the infected dropped to 0, so when control is ended, there is no chance a new outbreak can start. However, that changes if newly infected enter the population. To simulate this, set _tnew=50_, then run the simulation again. A new infected person is no introduced every 50 days. At times 50 and 100, control is still in effect, so they have no impact. But at day 150, there is no control, so this person can spark a new outbreak, which is what you see."
nrec = 2 # number of items to record
out_records = c( "Number of susceptible at end of control (t=110)",
                  "Number of susceptible at end of simulation")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report as 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 10
########################
tid = tid + 1
tasktext = "Open exploration: In a setting like this, where multiple outbreaks are possible and one cannot apply control for a long enough time to drive the disease extinct (or get a vaccine or some other new intervention), the best one can do is implement control to get the susceptible to drop to herd immunity while minimizing the _overshoot_, i.e. the excess drop of susceptible below herd immunity. This can be accomplished in different ways. Play around with start- and end-times for control and control strength to explore ways you can minimize the outbreak. You can also alter transmission or recovery rate to change _R~0~_ and see how that changes results. An option that is not possible with this simulation but could be done in real life is adaptive control by changing the control strength. You could mimic it by running a simulation with control at some level, then use the end values of the simulation as starting values with new control, etc. That's a bit tedious to do through the graphical interface but would not be too hard if you interacted with the simulation function directly through code."
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.
• This app closely follows a model and analysis described in [@handel07a]. Similar ideas applied to drug sensitive and drug-resistant influenza strains are discussed in [@lipsitch07; @alexander2007; @handel09a].

References

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