Complex ID control Scenarios - 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}

For this module, we will explore a fairly complex model that allows the study of different types of interventions. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab.

Learning Objectives

• Assess how control applied to different modes of transmission impacts outcomes

The Model {#shinytab2}

Model Overview

This model is fairly big and has many parts that can be turned on or off depending on parameter settings. The model allows for 3 types of transmission: direct, through an environmental stage, and through a vector stage. The (human) host is modeled in some detail, the environment and vectors are modeled with 1 and 2 compartments. The following compartments are included:

• S - uninfected and susceptible individuals
• P - individuals who are infected and do not yet show symptoms. Those individuals can potentially be infectious
• A - individuals who are infected and do not show symptoms. Those individuals can potentially be infectious
• I - individuals who are infected and show symptoms. Those individuals are likely infectious, but the model allows to adjust this, including no infectiousness.
• R - recovered/removed individuals. Those individuals have recovered and are immune. They can lose their immunity in this model.
• D - individuals who have died due to the disease.
• E - pathogen in the environment
• S~V~ - susceptible vectors
• I~V~ - infected/infectious vectors

The included processes/mechanisms are the following:

• Susceptible individuals (S) can become infected by pre-symptomatic (P), asymptomatic (A) or symptomatic (I) hosts. The rates at which infections from the different types of infected individuals (P, A and I) occur are governed by 3 parameters, b~P~, b~A~, and b~I~.
• Susceptible individuals (S) can also become infected by contact with the environment or infected vectors, at rates b~E~ and b~v~.
• Susceptible vectors (S~v~) can become infected by contact with symptomatic hosts at rate b~h~.
• All infected hosts first enter the presymptomatic stage. They remain there for some time (determined by rate g~P~, the inverse of which is the average time spent in the presymptomatic stage). A fraction f of presymptomatic hosts move into the asymptomatic category, and the rest become symptomatic infected hosts.
• Asymptomatic infected hosts recover after some time (specified by the rate g~A~). Similarly, the rate g~I~ determines the duration the symptomatic hosts stay in the symptomatic state. For symptomatic hosts, two outcomes are possible. Either recovery or death. The parameter d determines the fraction of hosts that die due to disease.
• Recovered individuals are initially immune to reinfection. They can lose their immunity at rate w and return to the susceptible compartment.
• Symptomatic and asymptomatic hosts shed pathogen into the environment at rates p~A~ and p~I~. The pathogen in the environment decays at rate c.
• New susceptible hosts and vectors enter the system (are born) at rates n~h~ and n~v~. Mortality (death unrelated to disease) for hosts and vectors occurs at rates m~h~ and m~v~.

Note that we only track people that die due to the disease in our D compartment. All hosts dying due to other causes just "exit the system" and we don't further keep track of them (though we could add another compartment to "collect" and track all individuals who died from non-disease-related causes.)

Also, note that we made several simplifications to keep the model from getting too complex. For instance, presymptomatic individuals do not shed into the environment, and only symptomatic hosts are assumed to be able to infect vectors. Further details relaxing these assumptions could, of course, be included, at the expense of a larger and more complex model.

Model Comments

In the Environmental Transmission app, I used P for pathogen in the environment, but here I'm using E since P is used for pre-symptomatic. One could as well switch things around and use E for exposed (instead of my P) and then P for pathogen in environment. I decided to leave it as is. It's good to get used to different notation.

I also decided for simplicity to not give the human compartments (or the environmental pathogen) any subscripts, so no subscript means human, if a compartment has a subscript it refers to vectors (i.e. S~v~ and I~v~).

Another point worth mentioning are the transmission terms. I generally use a single subscript to describe transmission from a group, e.g. b~A~ for transmission/infectiousness of asymptomatic. If there are multiple groups that can be susceptible and infectious, a common notation is to start with the receiving group first, then the sending/transmitting group, e.g. if susceptible adults (A) are infected by children (C), most authors write b~AC~. I follow this convention. Note however that it is equally ok to use b~AC~ to mean that infected adults transmit to susceptible children. I actually like this sender first perspective better, and used it originally, but switched to stick with the convention used in the main introductory textbooks on this topic.

In general you need to read papers/model descriptions carefully, and hopefully the authors do a good job explaining exactly what is meant. Such that there is no confusion. Just read carefully every time and don't jump to conclusions based on what you have seen before or what you think it means.

Model Implementation

The flow diagram and equations describe the model implemented in this app. Note that births and natural deaths are not drawn to keep the diagram from getting too cluttered.

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

$$\dot S = n_h - S (b_P P + b_A A + b_I I + b_E E + b_v I_v) + wR - m_h S $$ $$\dot P = S (b_P P + b_A A + b_I I + b_E E + b_v I_v) - g_P P - m_h P$$ $$\dot A = f g_P P - g_A A - m_h A$$ $$\dot I = (1-f) g_P P - g_I I - m_h I $$ $$\dot R = g_A A + (1-d) g_I I - wR - m_h R$$ $$\dot D = d g_I I $$ $$\dot E = p_I I + p_A A - c E $$ $$\dot S_v = n_v - b_h I S_v - m_v S_v $$ $$\dot I_v = b_h I S_v - m_v I_v $$

What to do {#shinytab3}

Notes

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

Some of the simulations might take a few seconds to run. Be patient.

#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 = "Set the model parameters such that it corresponds to the following setting: 1000 susceptible hosts and vectors, 1 initially symptomatic host, all other compartments start at 0. Assume that there are no births or natural deaths occurring for either hosts or vectors. Assume that only symptomatic individuals transmit, at a rate of 0.002. All other transmission rates should be 0. Assume that the duration of the symptomatic period is 1 month long, the duration of the presymptomatic period is half a month long. Assume that there are no asymptomatic infections. You can, therefore, set the rate of recovery of asymptomatics to anything, it doesn't matter because nobody will be asymptomatic. Assume that no environmental shedding and decay occurs. Assume nobody dies due to disease, and immunity does not wane.  Simulation duration approximately 5 years. With parameters set to correspond to the scenario just described, run the simulation and ensure you get a single outbreak with ~20% susceptibles left at the end."
nrec = 1 # number of items to record
out_records = c("Peak/Maximum number of infected symptomatic (I).")
out_types = c("Rounded_Integer")
out_notes = c("Report the rounded 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 2
########################
tid = tid + 1
tasktext = "Let's now assume that 50% of infected hosts will be asymptomatic (_f_) and that the duration of the asymptomatic stage (_gA_) is the same as the symptomatic stage (_gI_). Also, assume that asymptomatic infected are half as infectious as symptomatic infected, and that pre-symptomatic are as infectious as symptomatic. Run the simulation, you should get an outbreak with around 107 susceptibles left. Now, we will envision a scenario where we can only detect and isolate individuals that show symptoms. Assume that isolation reduces transmission rates by half. Implement such a scenario, run the simulation. 

If you are up for a small, optional, challenge, compute R~0~ for this scenario. It will be the sum of transmission coming from all 3 infected compartments (P, A, I). You will also need to make sure you multiply _A_ and _I_ by their respective fractions (_f_ and _1-f_) otherwise you would be double-counting. Other than that it follows the standard procedure you learned about and you can compare what you get from doing the math with the final size equation. You can then use your R~0~ equation for the next following tasks."
nrec = 1 # number of items to record
out_records = c("Number susceptibles remaining at end of isolation simulation.")
out_types = c("Rounded_Integer")
out_notes = c("Report the rounded 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 = "Set parameter values back to what they were at the start of the previous task. Now, we envision a scenario where we can isolate everyone who has become infected, independent of symptom status. We assume this reduces transmission by half for anyone isolated. Implement such a scenario, again run the simulation and record the number of susceptibles left at the end."
nrec = 1 # number of items to record
out_records = c("Number susceptibles remaining at end of simulation.")
out_types = c("Rounded_Integer")
out_notes = c("Report the rounded 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 = "Set parameter values back to what they were at the start of the previous task before you started changing them. Now assume that we can administer a drug. It will likely only be given to symptomatic individuals. First assume that the drug reduces the infectiousness of symptomatic individuals by half. Run the simulation. Convince yourself that this has the same effect as assuming isolation of of symptomatic individuals above. Now assume that the drug also reduces the duration of the symptomatic period from a month to half a month. What do you expect this shortening of the symptomatic period will do? Run the simulation to check."
nrec = 1 # number of items to record
out_records = c("Number susceptibles remaining at end of simulation.")
out_types = c("Rounded_Integer")
out_notes = c("Report the rounded 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 5
########################
tid = tid + 1
tasktext = "Set parameter values back to what they were at the start of the previous task. Then change transmission terms such that bP=0.001, bA=0.002, bI=0.001. Can you think of a pathogen where asymptomatic individuals might transmit more than symptomatic individuals? Some STI might fall into this category. Pathogens that require close contact could also be of that type if symptomatic are at home sick and can't spread the disease. Run the simulation with those settings, you should see an outbreak with around 203 susceptible left. Now assume that we can administer a drug which reduces symptoms in people, but doesn't change their transmission rates. We mimic this by assuming the fraction that are asymptomatic increases to 75%. What do you expect this will do to the overall outbreak size? Run the simulation to check."
nrec = 1 # number of items to record
out_records = c("Number susceptibles remaining at end of simulation.")
out_types = c("Rounded_Integer")
out_notes = c("Report the rounded 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 = "Now let's consider a pathogen that has some environmental transmission component. Set everything as in task 1. Then turn on environmental shedding by symptomatic individuals by setting pI=2 and set decay of pathogen in environment to c=1. Run the simulation. You should see E come up and go back down similar to the I curve. Make a note of the maximum value for E, it should be around 194. Since we haven't turned on transmission from the environment yet, the outbreak size should be the same as in task 1. Now set bE=0.0005. You will see a larger outbreak (less susceptible hosts remaining) and also higher levels of E. Why do you see a larger outbreak and more E now? Think about how the pathogen moves through the system.

Again, if you are up for an optional challenge and work on your R~0~ skills, determine R~0~ for this scenario. You can use the discussion in the _Environmental Transmission_ app to guide you. Since here we have transmission both through environment and direct transmission, R~0~ will be the sum of transmission coming from all 3 infected compartments (P, A, I) and the environement E. Though for this task, we have no transmission from P and A, so the sum is only the transmission terms from I and E."
nrec = 2 # number of items to record
out_records = c("Number susceptibles remaining at end of simulation.",
                "Max/peak of E during 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 7
########################
tid = tid + 1
tasktext = "Let's repeat the previous task, with only environmental transmission possible. Start by setting transmission bI=0 and bE=0, leave everything else unchanged. Run the simulation. Nothing should happen. The single infected person sheds a bit into the environment, but since we haven't turned on transmission from the environment yet, nothing else happens. If you zoom into the plot (with plotly) to the beginning, you see I drop to 0 and E come up a little bit and then drop as well. Now set bE=0.002. You will see an outbreak, now driven entirely through environmental transmission (compared to a mix of direct and environmental transmission in our previous task)."
nrec = 2 # number of items to record
out_records = c("Number susceptibles remaining at end of simulation.",
                "Max/peak of E during 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 8
########################
tid = tid + 1
tasktext = "Time to apply control. Keep settings as you just had them, and assume an intervention is able to remove pathogen faster from the environment (e.g. through better waste removal/treatment). Double the rate of pathogen decay in the environment."
nrec = 2 # number of items to record
out_records = c("Number susceptibles remaining at end of simulation.",
                "Max/peak of E during 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 9
########################
tid = tid + 1
tasktext = "Set pathogen decay rate back to what it was at the beginning of the previous task. Now we assume that we have an intervention that reduces shedding of pathogen by infected (e.g. through better access to sanitation). To mimic this, halve the rate of environmental shedding by symptomatic individuals. Run the simulation."
nrec = 2 # number of items to record
out_records = c("Number susceptibles remaining at end of simulation.",
                "Max/peak of E during 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 10
########################
tid = tid + 1
tasktext = "As you compare the previous 2 interventions (increase in _c_ or reduction in _pI_), which one was better? Explore a bit more how different choices of those two parameter affect the outbreak. If you computed R~0~ for this scenario, you can directly read off from the equation how those two parameters affect R~0~ and thus outbreak size."
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 11
########################
tid = tid + 1
tasktext = "Now we'll switch to vector-borne transmission. Set everything as in task 1. Also set initial number of susceptible and infected vectors to 1000 and 1. Assume that direct transmission between hosts does not occur, so all those rates should be 0. Same for environmental transmission. Set transmission from host to vector and vector to host to 0.002. Assume that vectors (say mosquitoes) live for half a month. Set vector birth rate such that vector population balances at 1000 in the absence of any infections. Run the simulation for 10 years. You should end up with around 228 susceptibles at the end of the simulation. Now consider some vector control measures. Assume we sprayed against mosquitoes and it reduced the population size by 90%. Set the initial susceptible vector population to that reduced value. Run the simulation. You should find that this does not have much of a long-term impact because the vector popoulation quickly rebounds."
nrec = 2 # number of items to record
out_records = c("Number susceptible hosts remaining at end of simulation with control.",
                "Number susceptible vectors remaining at end of simulation with control.")
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 12
########################
tid = tid + 1
tasktext = "Instead of killing vector populations, we now consider a reduction in transmission, e.g. due to the use of bed nets. Set everything as at the beginning of the previous task. Assume that an intervention reduces transmission _to_ vectors to 0.0015, leave the transmission from vectors unchanged. Run the simulation, observe. Now assume that an intervention only reduces transmission _from_ vectors to 0.0015, the rate _to_ vectors now remains at 0.002. Run the simulation, observe. Finally, assume that the intervention reduces transmission both to and from vectors to 0.0015. Run the simulation, observe."
nrec = 2 # number of items to record
out_records = c("Number susceptible hosts remaining at end of simulation with both transmission rates reduced.",
                "Number susceptible vectors remaining at end of simulation with both transmission rates reduced.")
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 13
########################
tid = tid + 1
tasktext = "While you just worked through a lot of tasks, there are still model parameters and thus scenarios that we haven't touched on. We didn't explore impact of disease-induced death, waning immunity, or presence of all three modes of transmission (though I'm not sure what ID that would apply to). All of those settings, and the impact of different interventions, could be explored. If you have a favorite pathogen, you could try to set parameters in this model to mimic it and explore different interventions and their impact." 
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 some additional reading, see e.g. [@klepac2013, @klepac2015, @kirsch2016, @tognotti2013].

References

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