Stochastic SIR Model

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 lets you explore a stochastic SIR model. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab.

Learning Objectives

• Get an introduction to stochastic compartmental models
• Understand differences and similarities between stochastic models and their deterministic equivalent
• Understand the concept of reproducibility with random numbers
• Know how the reproductive number affects outbreaks in stochastic 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. This is the same model we used for the Reproductive Number 2 app.

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 that are immune.

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 for the model implemented in this app is:

knitr::include_graphics( system.file(figuredir,appsettings$modelfigname,package=packagename))

In this app, you can run the model both as an deterministic, ordinary differential equation implementation (the kinds of simulations you have mostly seen so far), as well as a stochastic version of that model.

The deterministic model implemented as set of differential equations is given by the following equations:

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

The main focus for this app is a stochastic implementation of the variables and processes described above. This model is not an ordinary differential equation model. It is instead its stochastic equivalent.

We can specify the model by writing down every possible transition/event/reaction that can occur and their propensities (the propensity multiplied with the time step gives the probability that a given event/transition occurs). For our model these are the following:

Event type | Transitions | Propensity | ---------- | ----------- | ---------- | Infection | S => S-1, I => I+1 | bSI | Recovery | I => I-1, R => R+1 | gI | Births | S => S+1 | n | Death of susceptible | S => S-1 | mS | Death of infected | I => I-1 | mI | Death of recovered | R => R-1 | mR | Waning immunity | R => R-1, S => S+1 | w*R |

A note on randomness in computer simulations

This simulation (as well as some of the others) involves using random numbers to introduce stochasticity/uncertainty/noise into the model. This leads to a model that usually more closely reflects the underlying real system. However, in science, we want to be as reproducible as possible. Fortunately, random numbers on a computer are not completely random, but can be reproduced. In practice, this is done by specifying a random number seed, in essence a starting position for the algorithm to produce pseudo-random numbers. As long as the seed is the same, the code should produce the same pseudo-random numbers each time, thus ensuring reproducibility.

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

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 = "Start with 1000 susceptibles, 1 infected, no recovered. Transmission/infection rate of 0.0005, duration of infectious period 5 days, no births and deaths, no waning immunity. Start time 0, time step 0.1, simulation time of 100 days. Set number of simulations and random seed to anything you like. Run the deterministic/ODE model. Hit the _Run_ button multiple times with changed random number seeds. You can also change the number of simulations. Since the ODE model doesn't include any randomness and each time it leads to the same outcome, the random number seed and simulation number inputs are actually ignored (you can check that by looking at the function call for the ODE model for this app, see the _Further Information_ section). 

Now switch to the stochastic model, single simulation, random number seed at 100. Hit the _Run_ button multiple times. What do you expect? What do you see? If you are surprised, re-read _A note on randomness_ in the _Model_ section.

Next, set the number of simulations to 3 and run. You will see the one you already saw, another similar trajectory, and one scenario without an outbreak (note that susceptible stay at their starting value). What happens is that if you choose to run multiple stochastic simulations, the software increases the random number seed by 1 each time (i.e. it runs 3 simulations with random seeds 100, 101, and 102), otherwise you would get exactly the same result multiple times, which is useless. To confirm this, set the number of simulations to 1 and run the model with random seed 100, 101 and 102. Compare to what you got when you run 3 simulations starting at seed 100. 

Try to understand why for some of the stochastic simulations (seed 100 and 102) you get an outbreak, while for others (seed 101) not. Keep playing around with different seeds and different number of simulations until you have a full understanding of what's going on."
nrec = 2 # number of items to record
out_records = c("Final number susceptible, random seed 100, single run of stochastic model",
                "Final number susceptible, random seed 102, single run of stochastic model")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the 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
tasktext = "Set the random seed to 100, 5 simulations, and choose to run both models. All other parameter values as before. You'll get a plot showing the deterministic model dynamics and 5 runs for the stochastic model (each of the stochastic runs gets consecutive random number seeds, 100, 101,...). You might want to switch to plotly for the figure so you can easier turn on/off specific curves (by clicking on them). Take a look at the final value for the susceptibles for both the deterministic and stochastic models. Does what you see in the figure and the numbers/text underneath the figure make sense? If they don't at first glance, carefully think about the last sentence under the figure and what it means."
nrec = 2 # number of items to record
out_records = c("Final number susceptible, deterministic model",
                "Average final number susceptible, stochastic model")
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 = "We previously discussed that R~0~<1 means you won't get an outbreak, while R~0~>1 means you do get an outbreak. The first statement is true, but the second one is only strictly true in a deterministic model. I've tried to be careful and write that with R~0~>1 you _can_ get an outbreak. But it doesn't always happen, which you just saw in task 1 with the stochastic model. 2 of the 3 runs led to an outbreak, the third did not. Since you were able to get an outbreak, it means we have R~0~>1. Confirm this by computing R~0~ for this model. Still, it can happen that no outbreak occurs. That's because if you start with a single infected, and that person by chance recovers before infecting others, the outbreak is over. Thus, for a stochastic model (which is closer to the real world than a deterministic model), there is a chance you get an outbreak for R~0~>1, but it's not guaranteed. We'll explore that more in the next task."
nrec = 1 # number of items to record
out_records = c("R0 for the model")
out_types = rep("Numeric",nrec)
out_notes = rep("Report to one decimal place",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 = "One can mathematically derive an equation linking the probability that an outbreak occurs, _p_, to the reproductive number (see e.g. [@keeling08]). We won't try to do the math here but instead see if we can use the simulation to figure it out. Set parameters as in task 1. Then run the stochastic model for b = 0.0002, 0.0003, 0.0004, 0.0005, 0.0006, 0.001, 0.002. For each value, compute R~0~ and write it down. Then run 100 stochastic simulations for each _b_ value, with max time 300 and record the number of times you get an outbreak. Plot the relationship (either by hand or using the computer) between _R~0~_ and the probability of getting an outbreak. Note that running 100 stochastic simulations can take a bit of time, so be patient. If your computer is too slow, you might try to only run 50 or less, though with fewer simulations the estimate of the proportion of outbreaks is not as accurate as with more simulations."
nrec = 1 # number of items to record
out_records = c("The relation between R0 and the probability of an outbreak is close to a straight line.")
out_types = rep("Logical",nrec)
out_notes = rep("Report TRUE or FALSE",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 = "There is a similar relation between outbreak probability and the initial number of infected individuals. Instead of doing math, let's again use the model to determine it. Keep _b_ fixed at 0.0003, but now vary the initial number of infected individuals, _I~0~_, from 1 to 10. Everything else as in the previous task. For each value of _I~0~_, run 100 simulations, record the fraction of outbreaks you get for each _I~0~_. Plot the relationship (either by hand or using the computer) between _I~0~_ and the probability of getting an outbreak."
nrec = 1 # number of items to record
out_records = c("The relation between I0 and the probability of an outbreak is close to a straight line.")
out_types = rep("Logical",nrec)
out_notes = rep("Report TRUE or FALSE",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 = "Keep exploring. So far we haven't turned on births and deaths or waning immunity. If you do so, you can get steady states. Note that for the stochastic model, there is technically not a real steady state since the numbers bounce around. Eventually, just by chance, the numbers for the stochastic model will hit 0 and then won't recover from there. However, that could take a very long time (approaching infinity). So in practice the stochastic model often does exhibit something like a steady state, with fluctuations around that state."
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 more details on the equations connecting outbreak probability and R~0~ see e.g. the stochastic dynamics chapters in the IDEMA book or [@keeling08].
• For some further reading concerning stochastic dynamics of ID, see e.g. [@lloyd-smith05; @black1966; @keeling1997].

References

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