Basic Bacteria Exploration

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 you to explore the effect of specific model parameters on some outcomes of interest for the Basic Bacteria Model. I assume you worked your way through that app. If not, do that first. Then, read more about the scenario in this app in The Model tab. Then, work through the tasks described in the What To Do tab.

Learning Objectives

• Learn the concept of model exploration for a range of parameters.

The Model {#shinytab2}

Model Overview

The model used here is the one introduced in the Basic Bacteria Model app. If you haven't done so, check out and explore that app first. Here, we are only considering the continuous time, ordinary differential equation implementation of the model.

Model Diagram

As mentioned, the underlying model is the continuous-time, ODE version of the model from the Basic Bacteria Model app. Details can be found there. For ease of reference, the flow diagram is shown again here:

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

Model Equations

Set of ODE equations:

\begin{align} \dot B & = gB(1-\frac{B}{B_{max}})-d_B B - k BI \ \dot I & = r B I - d_I I \end{align}

Model Concepts

What's different here compared to the Basic Bacteria Model app is that instead of running the simulation once for a given choice of parameter values and looking at the resulting time-series, this app lets you explore the impact of each model parameter on some outcomes.

For the Outcomes of Interest, we consider the maximum number of bacteria and immune response and their values at the end of the simulation (when the system may have settled down to a steady state). You learned about steady states in the Basic Bacteria Model app. To ensure a steady state, you need to run the simulation for a long enough time. A steady state corresponds biologically to a chronic infection condition. In this app, those 4 steady-state outcomes of interest are labeled Bpeak, Ipeak, Bsteady and Isteady.

The Inputs of Interest are the different model parameters. You can choose one of the model parameters to be varied between some minimum and maximum value. The other parameter values remain fixed. For each parameter value, the model is run and the 4 outcomes of interest computed. The resulting plot shows how those 4 outcomes of interest vary with the parameter you investigated.

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

#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 initial bacteria to 100 and immune response to 10. Set bacteria growth rate to 2. The value does actually not matter since we'll be exploring the simulation for different values of _g_. Assume that bacteria live for a day, set the death rate accordingly (remember to convert duration to rate). Set carrying capacity to 10^5^. Set both the rate at which bacteria are killed by the immune response (IR) and the rate at which IR is activated and grows to 10^-4^. Set immune response decay rate to 2. Choose _g_ as the parameter to vary, go from 2 to 10, do 10 different parameter values/samples, linear spacing, and no log scales for plotting. Set maximum simulation time to 300. Run simulation, check to ensure all simulations reach steady state (see output message below plot). Then, reduce simulation time to 20 days; note some simulations did not reach steady state and how that affects results."
nrec = 1 # number of items to record
out_records = c("Number of simulations that did not reach steady state for _tfinal_=20.")
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 2
#########################
tid = tid + 1
tasktext = "Set simulation time back to 300, leave everything else the same. Play around with the number of parameter samples. A higher number gives you a more detailed plot, but it also means the simulation runs more often and it takes longer. Next, explore how changing the minimum and maximum of the parameter range you scan over changes the results. Then, explore how spacing of parameter values (linear versus log) changes results. Compare and contrast that with plotting results on a linear versus log scale. 

\nThink a bit about the results you see in the plot for the different outcomes. Would you have expected that as the bacteria growth rate increases, the peak bacteria load first increases, but then levels off? What about the fact that the steady state value for the bacteria does not change as the growth rate changes? And that the immune response peak and steady state both increase, even though the parameter you vary is the bacteria growth rate. Is that surprising? At this point, you might want to go back to the steady state tasks of the Basic Bacteria app and compare what you found there with what you see here. You should expect the mathematical expressions and the numerical results to agree. You can use the mathematical expression to check the result you are asked to report and read off from the plot."
nrec = 1 # number of items to record
out_records = c("__Bsteady__ for pretty much all values of _g_")
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 = "As you learned in the ***Basic Bacteria Model*** app, this model is simple enough that we can compute the steady states mathematically, without needing to run simulations. However, even for a simple model like this, one can't write down similarly simple equations for the maximum of __B__ and __I__. You can take a look at the model equations and diagram and see if by thinking through them, you can understand why you got the results for the peak values you find when running the simulations. With enough practice, it is often possible to intuit specific results based on the underlying equations, but that requires a good bit of modeling practice. As models get more complicated, even experienced modelers can often not intuit what model behavior one should expect as specific parameters vary. Often the only way to find out is by actually running the simulations. Let's explore this a bit for a different parameter. 

\nReset all inputs to their defaults. Then, set the parameter to vary to be *d~I~*. Choose 1 to 10 as the range for the parameter, 20 samples. For the __Bsteady__ and __Isteady__ values, you can compare the results with the mathematical expressions. For the peak values, you can try to intuit what you might expect _B~peak~_ and _I~peak~_ do as you increase _d~I~_ based on looking at the model equations. You can then check by running the simulations. You'll see that _I~peak~_ stays at 10 (the initial value) above _d~I~ = 5_. Try to understand what is going on here. Also, as _d~I~_ increases,  _B~steady~_ approaches _B~peak~_ and both reach and remain at the same value. Try to understand why that is. (Hint: At high _d~I~_, the immune response is so small it has no more effect in the dynamics of the bacteria.)."
nrec = 1 # number of items to record
out_records = c("__Bsteady__ for values of _d~I~_ > 5 (get it from plotly plot)")
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 4
#########################
tid = tid + 1
tasktext = "Do the same explorations you did above for any other parameter you want to investigate. Note that here we focus on a single parameter at a time. With some coding effort, we could change the underlying simulation to loop over say 2 parameters and produce outcomes for sets of parameter values, e.g. _B~peak~_ as a function of _g_ and _r_. The results could be plotted as a 2-dimensional heatmap. While this could be extended to more than 2 parameters, it will become hard to visualize and long to run. If there are many parameters that could change, a different approach is useful; you'll learn more about this in the ***Uncertainty and Sensitivity Analysis*** app."
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}

For this app, the underlying function running the simulation is called r appsettings$simfunction. That function repeatedly calls r appsettings$underlying_function.

This app (and all others) are structured such that the Shiny part (the graphical interface you are using) calls one or several underlying R functions which run the simulation for the model of interest and return the results. 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.

Some 'real world' examples where models have been used to explore outcomes of interest as function of parameters for simple models applied to TB, see e.g. [@antia96] (specifically Figure 2B) or Malaria, see e.g. [@kochin10].

References

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