Characterizing Infection States - Documentation

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 allows you to explore a model that has several more compartments than the basic SIR model discussed in the Introduction to ID app. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab.

The Model {#shinytab2}

Model Overview

In the Introduction to ID app, you explored a simple 3-compartment model, the basic SIR model. The model for this app has a few additional compartments, which allows us to include more details/realism into our model. We again focus on tracking individuals with regard to their infection/disease status. For this model, we track the following compartments/stages:

• S - susceptible, uninfected individuals
• P - presymptomatic individuals who are infected and do not yet show symptoms. Those individuals can potentially be infectious.
• A - asymptomatic, infected individuals. 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.
• D - individuals who have died due to the disease.

Of course, as with the basic SIR model, we could include further details by extending the number of compartments. In general, for each additional feature you want to track, the existing number of compartments needs to be replicated by the discrete categories you have. For gender, one would need to have 2x the compartments. Similarly if one wanted to stratify according to young/medium/old age, 3x the compartments are required, etc.

In addition to specifying the compartments of a model, we need to specify the dynamics determining the changes for each compartment. In general, more compartments leads to more processes and more parameters governing these processes.

For this model, we include the following processes:

• Susceptible individuals (S) can become infected by pre-symptomatic (P), asymptomatic (A) or symptomatic (I) hosts. Different types of infected individuals (P, A and I) transmit to other individuals and infect them at rates governed by the parameters, b~P~, b~A~, and b~I~. In other words, those b~i~ parameters determine how infectious an individual in stages P, A and I is.
• All infected individuals 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, infected symptomatic hosts leave that stage at rate g~I~. For symptomatic hosts, two outcomes are possible, either recovery or death. The parameter d determines the fraction of hosts that die.
• Recovered individuals are immune to reinfection.

Model Implementation

The flow diagram and the set of ordinary differential equations (ODE) which are used to implement this model are as follows:

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

$$\dot S = -S (b_P P + b_A A + b_I I)$$ $$\dot P = S (b_P P + b_A A + b_I I) - g_P P$$ $$\dot A = f g_P P - g_A A$$ $$\dot I = (1-f) g_P P - g_I I$$ $$\dot R = g_A A + (1-d) g_I I$$ $$\dot D = d g_I I$$

What to do {#shinytab3}

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 (e.g. if it says a week, you need to convert it to 7 days).

#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 = "Start with 1000 susceptibles, 1 initially infected presymptomatic host, _P~0~_, simulation duration of 200 days. Assume that only symptomatic individuals transmit, at rate _b~I~_ = 0.001. Assume that the duration of the presymptomatic, asymptomatic and symptomatic periods are all 5 days long. (Hint: The parameters _g~P~_, _g~A~_, and _g~I~_ are the inverses of these periods.) Assume that there are no asymptomatic infections, _f_=0, and nobody dies due to disease _d_=0. Make sure the results agree with model settings. For instance you shouldn't get any deaths, and no asymptomatic infected. If you did it correctly, your outbreak should end with around 7 susceptibles remaining. From the graph, contemplate how you would estimate the day at which the outbreak peaks. What's the problem? How would you solve it? Rerun the simulation, with the same input settings (i.e. don't change anything, just hit the Run Simulation button again). Do you see any changes in the results? Why (not)? (This will be different once we move to the stochastic models)."

# Record for task 1
nrec = 2 # number of items to record
out_records = c("Total number of recovered (at end of simulation)",
                    "Number susceptible left at 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 2
tid = 2
tasktext = "Assume now that after the pre-symptomatic stage, half of the infected move on to become symptomatic while the other half is asymptomatic. Set the model parameter controlling the fraction who become asymptomatic accordingly. Don't change any other assumption. Do you expect to get a change in the outbreak size? Why? Compute the total number of infected as you did for the _Basic SIR_ app. To do so, look at the diagram for the model and make sure you understand how individuals can flow through the system. That information will help you compute the total (cumulative) number infected (both symptomatic and asymptomatic)."

# Record for task 2
nrec = 2 # number of items to record
out_records = c("Number susceptible left at end of simulation",
                    "Total/cumulative number who got infected")
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 = 3
tasktext = "Now assume that asymptomatic individuals transmit at the same rate as symptomatic individuals. Change _b~A~_ accordingly. Leave everything else as before. How do you expect the results to change? (Try to make as precise/quantitative a prediction as you can). Compare your results to the first task and make sure you understand why you get the results you do."

# Record for task 3
nrec = 3 # number of items to record
out_records = c("Number recovered at end of simulation",
                    "Number susceptible left at end of simulation",
                    "Total deaths")
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 = 4
tasktext = "Next, let's assume that half the symptomatic infected die. Change the corresponding parameter in the model to reflect this new assumption. Leave everything else as in task 3. How do you expect the results to change?"

# Record for task 4
nrec = 3 # number of items to record
out_records = c("Number of recovered at end of simulation",
                    "Number susceptible left at end of simulation",
                    "Total deaths")
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 = 5
tasktext = "Set all 3 transmission parameters (_b~P~_, _b~A~_, and _b~I~_) to 0.0005, leave everything else as before. Note that this lowers the transmission rate of the asymptomatic and symptomatic groups by a factor of 2, but now pre-symptomatic individuals are also infectious. What do you expect will happen to the peak of _P_ and _I_ as well as the total number of individuals who got infected (i.e., _R_ + _D_ at the end of simulation), compared to the previous task?" 

# Record for task 5
nrec = 4 # number of items to record
out_records = c("Max number (peak) of pre-symptomatic (P)",
                    "Max number (peak) of symptomatic infected (I)",
                    "Number susceptible left at end of simulation",
                    "Total deaths")
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 = 6
tasktext = "Further explore how changes in the infectiousness of the different groups (_b~P~_, _b~A~_, _b~I~_) and the average time a person spends in each of those states (_g~P~_, _g~A~_, _g~I~_) affects the infection dynamics. Every time, think about what you expect to get, then run the simulation, compare your expectations with the results. Then make sense of it."

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 = 7
tasktext = "Using the settings from task 5, but now assume that 90% of infections are asymptomatic. How do things change?"
nrec = 4 # number of items to record
out_records = c("Max number of pre-symptomatic (P)",
                    "Max number of symptomatic infected (I)",
                    "Number susceptible left at end of simulation",
                    "Total deaths")
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 = 8
tasktext = "Further explore how changes in the fraction becoming asymptomatic and fraction dying does (or does not) affect the infection dynamics. Every time, think about what you expect to get, then run the simulation, compare your expectations with the results. Then make sense of it."

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.
• Some further information covering the concepts this app addresses can be found in [@keeling07, @vynnycky10, @kajita07, @milwid2016].

References

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