Stochastic Basic Virus Model

In DSAIRM: Dynamical Systems Approach to Immune Response Modeling

#load various variable definitions that are the same for each app
source('startup_script.R')
sapply(files_to_source, source) #source some helper files defined in the files_to_source variable
currentrmdfile = knitr::current_input()  #get current file name
appsettings = get_settings(currentrmdfile,appdocdir,packagename) #get settings for current app

Overview {#shinytab1}

This app allows exploration of a stochastic model that is almost identical to the deterministic model introduced in the Basic Virus Model app. In fact, both models are being run so comparison is possible. Read about the model in The Model tab. Then, work through the tasks described in the What To Do tab.

Learning Objectives

• Be familiar with stochastic compartmental models.
• Understand the advantages and disadvantages of stochastic models.

The Model {#shinytab2}

Model Overview

This model consists of 3 compartments and is almost identical to the deterministic basic virus infection model. We track the following entities, by assigning each to a compartment:

• U - uninfected cells
• I - infected cells
• V - (free) virus

In addition to specifying the compartments of a model, we need to specify the dynamics determining the changes for each compartment. Broadly speaking, there are processes that increase the numbers in a given compartment/stage, and processes that lead to a reduction.

For this model, we consider the following processes:

1. Uninfected cells are produced at some rate n and naturally die at some rate d~U~.
2. Virus infects cells at rate b.
3. Infected cells produce new virus at rate p and die at rate d~I~.
4. Free virus is removed at rate d~V~ or goes on to infect further uninfected cells.

Model Diagram

The diagram illustrating this compartmental model is shown in the figure (with $g=1$ for this model).

knitr::include_graphics(path = paste0("../media/",appsettings$modelfigname))

Model Equations

If we were to implement this model as a continuous-time, deterministic model, it would have the following set of ordinary differential equations:

$$ \begin{aligned} \dot U & = n - d_U U - bUV \ \dot I & = bUV - d_I I \ \dot V & = pI - d_V V - b UV \end{aligned} $$

Model Concepts

Instead of implementing a continuous-time, deterministic model via ordinary differential equations, we use a stochastic model here. For such a model, the differential equation formulation is not valid. One can write down an equivalent formulation as a stochastic model by specifying every possible process (also called transition/event/reaction) that can occur and their propensities (the propensity multiplied with the time step gives the probability that a given process/event/transition occurs). For our model these are the following:

Event type | Transitions | Propensity | ---------- | ----------- | ---------- | Production of U | U => U+1 | nU | death/removal of U | U => U-1 | d~U~U | infection | U => U-1, V => V-1, I => I+1 | b*UV | death if I | I => I-1 | d~I~I | production of V | V => V+1 | pI | removal of V | V => V-1 | d~V~V |

Notes

• The one difference to the deterministic virus model is that we do not include an additional conversion factor to account for experimental units. The reason is that this stochastic model makes the assumption that all entities/variables are in discrete units, i.e., the number of cells and number of (infectious) virions, and can only change by integer numbers. Because of that, our model requires that virus is in units of actual virion numbers and we don't need to convert. That also means that the parameter g of the deterministic model - which we run for comparison - is set to 1.

• The adaptivetau package/function we use to run the stochastic model does not allow specification of starting time and time step. It always starts at t=0 and each time step is chosen by the simulation routine based on the adaptivetau algorithm. Therefore, to allow comparison with the ODE, the start time for the ODE is fixed to 0. Also, the time step at which data is returned is fixed at tfinal/1000.

What To Do {#shinytab3}

The tasks below are described in a way that assumes everything is in units of 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 = "Run the simulation with the pre-set parameter values. It runs the deterministic model and you should see an infection happening. Switch to the stochastic model (under _Models to run_) with the same parameter settings and run the stochastic model. You should see no infection. Change the random number seed to 122 and re-run. You should get an infection. 

\nNow, set _Models to run_ to both. You should get infections from each model; however, they look somewhat different. The infection produced by the stochastic model takes off faster. It also leads to slightly more infected cells. Figure out how you can determine the total number of infected cells caused by the infection. Hint: The maximum number of infected cells is not the same as the total number. Another hint: Contemplate what happens to the uninfected cells, and therefore what the final value of the uninfected cells might tell you about the infected cells."
nrec = 2 # number of items to record
out_records = c("Total number of infected cells for ODE model",
                "Total number of infected cells for stochastic model (_rngseed_ = 122)")
out_types = c("Rounded_Integer", "Integer")
out_notes = c("Report the rounded integer", "Report the 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 = "With the settings as they just were, increase the number of simulations (for the stochastic model) to 2, re-run the simulation. Notice that you see two distinct outcomes for the two stochastic simulations. Also, look at the message below the plot (what the app calls an outbreak is what I call an infection here, basically a within-host outbreak). Based on that, make sense of the fact that the average total number of infected cells for the deterministic simulation is much higher than for the stochastic (which contrasts to the finding from the previous task)."
nrec = 1 # number of items to record
out_records = c("Average total number of infected cells for stochastic model")
out_types = rep("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 3
#########################
tid = tid + 1
tasktext = "Set number of simulations to 4, keep everything else the same. You'll get 2 outbreaks. The maximum value for the virus for the stochastic simulations is around half that of the deterministic model. Based on the text below the plot, you know that the stochastic results are the averages over all runs. Thus, you are adding two values that are similar to the ODE model to 2 values that are 0. Reporting the average (mean) over situations that both do and don't lead to infection might not be too insightful. It might be better to report the number/fraction of times you get an infection, and then just the average for some quantities for those simulations where an infection occurs. For these settings, is it true that the average of the virus peak for the stochastic simulations is above that of the ODE? You might want to switch to plotly and zoom in on the peaks to determine this." 
nrec = 1 # number of items to record
out_records = c("TRUE or FALSE: The average virus peak for the stochastic simulations that lead to an infection is above the ODE value.")
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 4
#########################
tid = tid + 1
tasktext = "Set _Model to run_ to stochastic, number of simulations 20. Run the simulation. The message under the plot should inform you that for  runs, an infection/outbreak occured. The reason we sometimes do not get an infection is that we start with a single virion, which sometimes, just by chance, goes extinct before the infection can take off. This does never happen for the deterministic model. If the infection starts with more virions, the chance of all of them going extinct is much less. Test that by starting with 5 initial virions, re-run the simulation. You should find 13 outbreaks/infections happening. Repeat for 10 virions. You can explore different random seeds and more simulations to convince yourself that an increase in initial virus leading to more infections is a robust result (and doesn't depend on specific choices of the random seed or number of samples)."
nrec = 1 # number of items to record
out_records = c("Number of outbreaks/infections when starting with 10 virions")
out_types = rep("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 5
#########################
tid = tid + 1
tasktext = "The starting number of virions has a strong impact on the chance an infection occurs. Some model parameters do as well. If the virus is fitter (e.g. higher virus production rate, faster rate of infection), the chances that an infection occurs also increases. Let's try that. 

\nReset all inputs. Turn on stochastic simulations. Set to 20 simulations. Run the model. You should now get 4 outbreaks/infections (the random number seed is different from above). Now set p=10. Re-run. The number of infections should have increased. 

\nReturn to p=5 and set b=0.0002. You will again get more outbreaks/infections. Again, play with the random number seed and the number of simulations to convince yourself that this is a general finding."

nrec = 1 # number of items to record
out_records = c("Number of outbreaks/infections for b=0.0002.")
out_types = rep("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 6
#########################
tid = tid + 1
tasktext = "With the model settings we explored so far, stochasticity/randomness matters in deciding if an infection occurs or not. Once it has started, the resulting time-series look fairly smooth and not too noisy. The reason for that is that the numbers of virus grows quickly, and there are a lot of cells; thus, in general, the law of large numbers applies and the overall behavior of the system is fairly well approximated by a deterministic model. We can change this by reducing the number of cells. Reset all inputs. Then set U=1000. To make up for this change, increase b to b=0.002. Run 20 simulations of the stochastic model. You will see that the time-series look much noisier, even at the peak. You should find that 4 simulations produce an infection. Record the average number of peak/maximum infected cells.

The important take-home message from this is that for small numbers (e.g., at the start or end of an infection, or when overall numbers are large), stochasticity is often important. Once numbers are fairly large (100s or more), then deterministic models are usually doing a good job describing the system dynamics."
nrec = 1 # number of items to record
out_records = c("Average number peak/max infected cells (round to integer).")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Round to 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 = "Keep exploring. You can for instance turn on births and deaths of uninfected cells to explore chronic/steady state conditions and see how stochasticity impacts those."
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 functions running the simulation are called r appsettings$simfunction[1] and r appsettings$simfunction[2]. 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('DSAIRM') into the R console.

References

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