Model Variant 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 exploration of the impact of different model formulations on simulation results. Read about the model in The Model tab. Then, work through the tasks described in the What To Do tab.

Learning Objectives

• Explore the impact of alternative model formulations on outcomes.
• Understand that sometimes, different reasonable model formulations can produce qualitatively different results.

The Model {#shinytab2}

Model Overview

For this model, we track uninfected cells, infected cells and (free) virus. Additionally we model the innate and adaptive immune response - which is of course a very simplified approximation to the real immune response.

We track the following entities by assigning each to a compartment:

• U - uninfected cells
• I - infected cells
• V - (free) virus
• F - innate immune response
• A - adaptive immune response

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. Those processes are sometimes called in-flows and out-flows.

For the purpose of this app, we specify several alternative processes that allow us to explore different model variants by 'turning on and off' specific components of the model.

We specify the following processes/flows:

1. Uninfected cells are produced at rate n, die naturally at rate d~U~ and become infected at rate b.
2. Infected cells die at rate d~I~ and produce virus at rate p.
3. Free virus is removed at rate d~V~ due to any unmodeled processes, or goes on to infect further uninfected cells at rate b.
4. In the absence of virus, the innate response is produced at a rate p~F~ and removed at a rate d~F~. In the presence of virus, the innate response additionally grows according to 3 alternative model formulations: A) Proportional to virus at rate f~1~ and saturating at a maximum level F~max~. B) Proportional to virus at rate f~2~, with a growth rate saturating at high levels of virus, determined by the saturation constant s~V~. C) Proportional to both virus and infected cells at rate f~3~, with a growth rate saturating at high levels of virus and infected cells, determined by the saturation constant s~V~.
5. The innate response can also act on the system in 3 different ways: A) By moving target cells into a "protected" state at rate k~1~ where those cells can not become infected any longer. B) By inducing death of infected cells at rate k~2~. C) By reducing production of virus particles at strength k~3~.
6. The adaptive response growth is also modeled according to 3 alternative model formulations: A) Proportional to the innate response at rate a~1~. B) Proportional to virus at rate a~2~, with a growth rate saturating at high levels of virus, determined by the saturation constant h~V~. C) Proportional to both virus and innate response at rate a~3~, with a growth rate saturating at high levels of virus and innate response, determined by the saturation constant h~V~.
7. The adaptive response can act on the system in 3 ways: A) By killing infected cells at rate k~4~. B) By killing infected cells at rate k~5~, with saturation of the maximum killing rate for high adaptive response levels, determined by the saturation constant s~A~. C) By killing virus at rate k~6~.
8. The adaptive immune response decays at a rate d~A~.

Model Concepts

Both the innate and adaptive response are modeled in a rather abstract manner. We could think of them as some kind of cumulative representation of each arm of the immune response, or alternatively a single dominant innate response component, e.g. interferon for the innate and CD8 T-cells for the adaptive response.

The idea explored in this app and implemented by this model is that results sometimes, but not always, change depending on different (biologically reasonable) ways the immune response is modeled. We can explore those different models by setting certain parameters describing alternative processes to a non-zero value, and all others to zero. We can then study how different model alternatives affect the outcome.

Obviously, the number of alternative models we could make that are biologically reasonable is virtually endless. The better the underlying biology of a given infection is known, the easier it is to pick one model formulation over another. In the end, for most infections, we still don't know enough to pick the "right" model. We often have to choose one or a few reasonable model candidates and hope that they approximate the underlying processes reasonably well.

Model Equations

Implementing the above described processes as a set of ordinary differential equations leads to the following model:

\begin{align} \dot U &= n - d_U U - bVU - k_1FU \ \dot I &= bVU - d_II - k_2FI - k_4 A I - k_5 \frac{A I}{A+s_A} \ \dot V &= \frac{pI}{1+k_3 F} - d_VV - bVU - k_6AV \ \dot F &= p_F - d_F F + f_1 V (F_{max} - F) + f_2 \frac{V}{V+s_V} F + f_3 \frac{VI}{VI+s_V} F \ \dot A &= a_1 F A + a_2\frac{V}{V+h_V}F + a_3 \frac{F V}{ F V + h_V} A - d_A A \end{align}

Model Diagram

The rather messy looking model diagram is here:

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

Since this is already a fairly complex model, it is worth spending some time mapping each component in the equations to the arrows and boxes in the diagram and the written description. By going back and forth between verbal formulation and mathematical notation, you will start to build some expertise in going from one to the other.

Notes

• It is recommended that before you work your way through this app and tasks, you first explore and do the tasks of the apps in the Basic Virus Models section, especially the Basic Virus Model app. Otherwise, this model and app will likely be too overwhelming.

What To Do {#shinytab3}

For the tasks below, it is assumed that the model is run 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 = "We'll begin with the basic virus model by turning off any immune response related component. Set all parameter values related to the immune response such that there is no __F__ and __A__ present at any time during the simulation. Set initial conditions to 10^5^ uninfected cells, no infected cells, and 10 virions. Set the infection rate to 10^-5^, no production or death of uninfected cells, lifespan of infected cells and virus of 1 day and 6 hours, respectively. The rate of virus production should be 100 per day. You should get a single acute viral infection with no immune response present and a maximum of 70217 infected cells."
nrec = 1 # number of items to record
out_records = c("Day at which the viral load peaks")
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 = tid + 1
tasktext = "Let's explore different mechanisms for the innate response. Keep initial level of innate response at 0, but set innate production and removal rate to 1. Set all innate growth parameters (the _f~i~_) and all innate actions (_k~1~_, _k~2~_, _k~3~_) to 0. Run the simulation, confirm that you get an innate response that settles at a steady state (balance between production and removal), but that there is no further increase and that the rest of the dynamics doesn't change (i.e. you should get the same maximum number of infected cells as above)."
nrec = 2 # number of items to record
out_records = c("Maximum number of infected cells", "Value for innate response at steady state")
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 = "Now, let's explore the different types of innate response induction by playing with the _f~i~_ parameters and the saturation parameter, _s~V~_. First, set _f~1~_ to 10^-5^ and _F~max~_ to 1000. Relate what you see in the plots to the equations so you get an idea of how different terms in the equations behave. Continue to keep the _k~i~_ at 0. This means the rest of the variables should not change. Make ten-fold changes to _f~1~_ and _F~max~_. Try to predict what will happen and then test your predictions by running the simulation."
nrec = 6 # number of items to record
out_records = c("Peak of **F** with _f~1~_ = 10^-5^ and _F~max~_ = 1000", "Day at which **F** peaks with _f~1~_ = 10^-5^ and _F~MAX~_ = 1000", "Peak of **F** with _f~1~_ = 10^-6^ and _F~max~_ = 1000", "Day at which **F** peaks with _f~1~_ = 10^-6^ and _F~max~_ = 1000", "Peak of **F** with _f~1~_ = 10^-5^ and _F~max~_ = 10000", "Day at which **F** peaks with _f~1~_ = 10^-5^ and _F~max~_ = 10000")
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 = tid + 1
tasktext = "Let's switch to the alternative innate response induction proceses. Set _f~1~_ back to zero. Then, set _f~2~_ to 3 and _s~V~_ to 10^5^. Play around with different values for _f~2~_ and _s~V~_, each time trying to predict how changes in the parameters will affect the model.  

\nNext, set _f~2~_ to zero and set _f~3~_ to 2. For any given value of _s~V~_, how does switching between innate response induction via _f~2~_ differ from _f~3~_? (*Hint: take another look at the model equations.*) Play around with different values for _f~3~_ and _s~V~_, again trying to predict the output and then running the simulation to confirm. Remember, each time you switch processes, 'turn off' the alternative ones."
nrec = 4 # number of items to record
out_records = c("Peak of **F** with _f~2~_ = 3 and _s~V~_ = 10^5^", "Peak of **F** with _f~2~_ = 3 and _s~V~_ = 10^7^", "Peak of **F** with _f~3~_ = 2 and _s~V~_ = 10^5^", "Peak of **F** with _f~3~_ = 2 and _s~V~_ = 10^7^")
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 = tid + 1
tasktext = "Explore the alternative representations of innate response induction by switching modeled process. Do you think that any of the modeled processes is more ***biologically reasonable*** than the others? Why or why not? Are you able to recreate similar dynamics regardless of which process is modeled? What happens if you turn more than one of the alternative processes on at once; is *this* ***biologically reasonable***? Can you think of improvements to the parameterization?"
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = rep("",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 = "Now, let's explore what happens when you have non-zero _k~i~_. First, set the innate induction parameters as follows: _f~3~_ = 5, _f~1~_ = _f~2~_ = 0, and _s~V~_ = 10^5^ (this way we can explore the _k~i~_ effects using the same innate activation processes). Try different values for the _k~i~_ parameters. At first, turn on one process at a time; after, mix and match. 

\nYou'll find that the action of the innate response now impacts the other variables, which in turn can impact further innate activation. Some of the resulting dynamics can get complex. Pay attention to how different processes of innate activation and innate action do or don't produce different overall dynamics."
nrec = 4 # number of items to record
out_records = c("Peak of **I** with _k~1~_ = 0.001, other _k~i~_ = 0", "Peak of **I** with _k~2~_ = 0.001, other _k~i~_ = 0", "Peak of **I** with _k~3~_ = 0.001, other _k~i~_ = 0", "TRUE or FALSE: you can track the extent of the infection (i.e., how many total cells were infected) by calculating the difference between the number of ***uninfected cells*** at the beginning and the end of the simulation for all non-zero values of _k~1~_, _k~2~_, and _k~3~_.")
out_types = c(rep("Rounded_Integer",3), "Logical")
out_notes = c(rep("Report the rounded integer",3), "Report TRUE or FALSE")
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 = "Finally, turn on the adaptive response. Note that in this model, there is no adaptive response without innate response (check it by setting the innate response to 0 while having non-zero adaptive growth rates, _a~i~_). Also, note that the alternative adaptive growth processes for _a~1~_ and _a~3~_ further need an initial non-zero value for the adaptive immune response to be able to grow (again, check it for yourself either by reviewing the model equations or running the simulation). 

\nTurn the innate response back on as in Task 6 with _k~2~_ = 0.001 and _k~1~_ = _k~3~_ = 0. Explore how different non-zero adaptive growth rates, _a~i~_, affect the adaptive response dynamics. Start by leaving the adaptive action parameters at 0. Then, fiddle with those parameters as well and set them to non-zero values. Also, pay attention to the impact of the saturation constants."
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = rep("",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 = "Now that we have played with each of the immune response processes, let's see if we can put together a bit more complicated model. 

\nSet _d~F~_ = 10, the innate response growth parameters to _f~3~_ = 15, _f~1~_ = _f~2~_ = 0, and _s~V~_ = 10^5^, and the innate response action parameters to _k~2~_ = 0.001 and _k~1~_ = _k~3~_ = 0. Set the adaptive response initial value to 1, adaptive response growth parameters to _a~3~_ = 5, _a~1~_ = _a~2~_ = 0, and _h~V~_ = 10^7^, and the adaptive response action parameters to _k~4~_ = 0, _k~5~_= 0.01, _k~6~_ = 0.1, and _s~A~_ = 1000. Can you intuit what the model output will look like? Run the simulation for 100 days. You should have about 3 uninfected cells remaining at the end of the simulation.

\nNow, let's add some processes for the birth and death of uninfected cells. Set the rate of uninfected cell production, _n_, to 10000 and the rate of natural death of uninfected cells to a value that would, in the absence of infection, equate to a constant 10^5^ number of uninfected cells throughout the simulation. If you have difficulty determining the corresponding value for _d~U~_, you can review steady states in the ***Basic Virus*** app. Run the simulation. You should get a pretty interesting pattern. How would you interpret this? Does this model, parameter or dynamics, seem ***biologically reasonable***?"
nrec = 1 # number of items to record
out_records = c("Peak of **A**")
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 = "As you'll notice, some of the specific model choices (i.e. specific parameters/terms in the model being turned on or off) lead to similar results, while other times results are quite different. For a specific system under study, some ways to formulate the immune response and the dynamics one gets might be better than others. Think about a system you are familiar with and consider which (if any) of the possible mechanisms implemented in this model might best describe that system. Keep playing with the model and see if you can recreate the system!"
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = rep("",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 

#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 running the simulation is 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('DSAIRM') into the R console.

Explorations of different models and their impacts on outcomes for acute virus infections in general can be found in e.g. [@li14] and a more detailed discussion of previously published models and comparison to data for influenza can be found in [@dobrovolny13].

References

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