In ahgroup/DSAIRM: Dynamical Systems Approach to Immune Response Modeling
#load various variable definitions that are the same for each app
source('startup_script.R')
#source helper functions (defined in the startup script)
#needs to happen inside each Rmd file since knitr starts a new environment
sapply(files_to_source, source)
#load the settings file for the current app
#so we can automatically include figure, list the functions in the further information section
#and use other information specific to the current app for the task table generation
currentrmdfile = knitr::current_input()
#currentrmdfile = "basicbacteria_documentation.Rmd" #for debugging
appsettings = get_settings(currentrmdfile,appdocdir,packagename)
Overview {#shinytab1}
This app allows exploration of a very basic bacteria (or other extracellular pathogens) infection model, with one compartment for bacteria and one compartment for the immune response. The main goal for this app is to get you familiar with simulation models, as well as the overall setup and ideas behind using these simulations, and how to run them. Read about the model in The Model tab. Then, work through the tasks described in the What To Do tab.
Learning Objectives
- Become familiar with simulation models.
- Learn about steady states and how to compute them.
- Understand how model parameters impact outcomes.
The Model {#shinytab2}
Model Overview
All the models we consider are compartmental models. Compartmental means that we place the agents/players/entities/individuals of interest into distinct compartments/categories. We then only track the total number of individuals in each of these compartments. In this simple model, we track bacteria and some (fairly abstract) notion of immune response strength, using the following notation:
- B - bacteria
- I - immune response
In addition to specifying the compartments or variables of a model, we need to specify the dynamics determining the changes for each compartment. Broadly speaking, there are processes or flows 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 our system, we specify the following processes/flows:
- Bacteria grow/divide at some maximum rate, g, and saturate as they approach some maximum carrying capacity, B~max~.
- Bacteria die at a natural death rate, d~B~.
- Bacteria are killed by the immune response at some rate, k.
- The immune response grows proportional to the number of bacteria and itself at some rate, r.
- The immune response decays at some rate, d~I~.
Model Diagram
A very good way to describe compartmental models and move from a verbal description toward a mathematical/computational formulation is by using diagrams. Being able to go back and forth between verbal description, diagram and mathematical/computational model is a crucial skill when building models.
The diagram for a compartmental model is often called flow diagram. It consists of a box for each compartment/variable (here B and I), and arrows pointing in and out of boxes to describe flows and interactions. For the model described above, the flow diagram looks as follows:
knitr::include_graphics(path = paste0("../media/",appsettings$modelfigname))
For the diagrams in this R package, solid arrows indicate physical flows, i.e. movement from a compartment to another (e.g. bacteria moving out of the compartment because of death, or bacteria increasing in the compartment due to growth), while dashed arrows indicate interactions without physical flow (e.g. infected cells killing bacteria). This notation is not universal and it is common in the literature to see no distinction made between these 2 types of flows and only solid arrows being used.
Model Equations
Next, we need to implement this verbal model and flow diagram in a way that it can be run on a computer. While we already specified some of the model features, e.g. that we only track compartments/total numbers and not individual bacteria, both the verbal model and diagram do not yet fully specify the underlying mathematical model. In fact, there are several ways the diagram could be translated into a mathematical/computer model. For this app, we consider two approaches, namely discrete-time deterministic and continuous-time deterministic models. In some of the other apps, you will see another variant, namely continuous-time stochastic model implementations.
Implementation I
The most common way to implement compartmental models is as a continuous-time, deterministic process, formulated as a set of ordinary differential equations (ODEs). Each compartment/variable gets an equation. The right side of each equations specifies the processes going on in the system and how they change the numbers in each compartment via inflows and outflows. For the model described above, the equations look like this:
\begin{align}
\dot B &= g B (1-\frac{B}{B_{max}}) - d_B B - k BI \
\dot I &= r B I - d_I I
\end{align}
We are using here the short-hand notation where a dot over the variable indicates differentiation with respect to time, and we also do not explicitly indicate that all variables depend on time. This is the most common notation. A more detailed, completely equivalent notation is:
\begin{align}
\frac{d}{dt}B(t) &= gB(t)(1-\frac{B(t)}{B_{max}})-d_B B(t) - k B(t)I(t) \
\frac{d}{dt} I(t) &= r B(t) I(t) - d_I I(t)
\end{align}
This notation specifies explicitly that B and I are functions of time, and that the differentiation is with respect to time. Since this is always the case in our models, it is not necessary to be that explicit and the more compact notation shown above is the one you will see in this package and it's also the most common throughout the literature.
Implementation II
Continuous time models implemented as ordinary differential equations are the most common types of models. However, other implementations of the above model are possible. One alternative formulation is a discrete-time deterministic equivalent to the ODE model. For such an implementation, the equations are:
\begin{align}
B_{t+dt} &= B_t + dt \left( gB_t (1-\frac{B}{B_{max}}) - d_B B_t - k I_t B_t \right) \
I_{t+dt} &= I_t + dt \left( r I_t B_t - d_I I_t \right)
\end{align}
In words, the number of bacteria and immune response at a time step dt in the future is given by the number at the current time, t, plus/minus the various growth and death/removal processes. The latter need to be multiplied by the time step, since less of these events can happen if the time step is smaller. As the time-step gets small, this discrete-time model approximates the continuous-time model above. In fact, when we implement a continuous-time model on a computer, the underlying simulator runs a smart version of a discrete-time model and makes sure the steps taken are so small that the numerical simulation is a good approximation of the continuous-time model. If you want to learn more about that, you can check out the deSolve R package documentation, which we use to run our simulations.
Model Concepts
In general, the entities that change in our model (i.e. here B and I) are called variables; they are variable and change during the simulation. To run a simulation, we need to specify the starting values for each variable. Those are often called initial conditions.
In contrast, the quantities that are usually fixed for a given scenario are called parameters. For this model, those are g, d~B~, k, r and d~I~. Values for parameters are chosen to match the known biology of a specific disease/scenario we want to model. All parameters need to be in the same time units, e.g. per day or per hour.
You might find for some parameter settings that numbers for bacteria/immune response/virus/cells/etc. drop below 1 (and possibly rebound later). You'll see that in one of the tasks below. If the units of our model are assumed to be in number of pathogen/cells, this of course is biologically not reasonable. Numbers less than 1 are an artifact of the underlying differential equation model. All ODE models (and the discrete time model we consider here) have that problem. A hack-y solution is to monitor the simulation and if a quantity drops below 1, set it to 0 or stop the simulation. A cleaner solution is to treat all entities (virus, cell) as discrete units and allow them to only change in integer steps (in a probabilistic manner). This approach will be discussed in the apps that use stochastic models. For most apps, neither approach is taken, so you'll see numbers less than 1. That's okay for the purpose we use the apps here. Just be careful if you use these kinds of models in your research you might need to pay attention to this issue.
Notes
-
This model is of course very simple. However, as you will see, even with just 2 entities and their interactions, we can get interesting dynamics. This model could be extended by introducing further compartments. For instance we could include separate compartments for the innate and adaptive immune response. You will encounter more detailed models in some of the other apps.
-
This and most other simulations/apps do not have natural time units (unless specifically stated). You could, therefore, assume that your model runs in units of days or minutes/hours/weeks/months..., based on what's most suitable for the system you want to study. You have to make sure that all your parameters are in the right time units. Always make sure to check if a given simulation can handle different time units or assumes specific ones. In general, I specify for each What To Do section what kind of units I have in mind for the different tasks.
-
Some of the tasks in the What To Do section (and in future apps) are fairly open-ended. The idea is that these tasks give you something to get started, but you should feel free to explore the simulations any way you want. Play with them, query them, go through iterations of thinking about what you expect, observing it, and, if discrepancies occur, figuring out why. Essentially, the best way to use these apps is to do your own science/research with them.
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 cell lives for a week, you need to convert it to 7 days, then take the inverse if you want a death rate, which is what usually goes into the model).
#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 100 initial bacteria and an initial level of 1 for the immune response. Assume bacteria grow at a rate of 1 per day (i.e., g = 1), the carrying capacity is 10^5^ and they live for about 2 days (remember that the inverse of the lifespan is the rate of death). Assume the immune response is activated and grows at a rate of 10^-4^, kills bacteria at a rate of 10^-4^ and decays at a rate of 2 per day. Set simulation start time to 0 and final time of 100 (which we assume to be days). Set the time step to 0.01. As you'll see, if we run an ODE model, the time step is only relevant for plotting, not for the underlying model simulation (while this is not the case for the discrete time model). Only run the continuous time ODE model. Plot both x- and y-axes on a linear scale (i.e, no log scales for now). You can stick with ggplot as the plot engine. Run the simulation, see what you get. You should see some oscillations and then the system settles down, with bacteria and immune response at the end of the simulation at around 19996 and 3003, respectively. Change the time step to 0.05, re-run the simulation. Record the number of bacteria at the end of the simulation (Hint: not much changes)."
nrec = 1 # number of items to record
out_records = c("Number of bacteria at end of simulation, _dt_=0.05")
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 = "Play around with the plot settings. Switch the plotting to have x-axis, y-axis or both plotted on a log scale. Also try plotly as the plot engine. Leave all other settings as before. You should see that while the look of the plot changes, the underlying numbers do not. This is something to be aware of when you see plots in papers or produce your own. The best plot to use is the one that shows results of interest in the clearest form. Usually, the x-axis is linear and the y-axis is either linear or logarithmic. One nice feature about plotly is that the plot is interactive and you can read off numbers. Use this feature to determine the day at which the bacteria and immune response have their second peak. (*Hint: day 32 for bacteria, a bit later for the immune response.*)"
nrec = 1 # number of items to record
out_records = c("Day of second peak for the immune response")
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 = "Go back to both linear scales for plotting. Keep all variable and parameter settings as before. Set _Models to run_ to _both_. This runs and shows both the continuous-time and discrete-time models. Start with a time step of 0.01. Run the simulation. You should see the results from the 2 models essentially on top of each other and barely distinguishable. You should see that the maximum number of bacteria is 49,513 for the ODE model, and very close for the discrete-time model."
nrec = 1 # number of items to record
out_records = c("Maximum number of bacteria for the discrete-time 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 4
#########################
tid = tid + 1
tasktext = "Now try different values for the time step, _dt_. Leave all other settings as before. You should notice that as _dt_ gets larger, the continuous-time model results remain the same, but the discrete-time models change and start moving away from those of the continous-time model. At a time step above 0.1, the results start to look very different. Somewhere above a time step of 0.5, it becomes so large that for these parameter settings, the simulation 'crashes' and you get an error message."
nrec = 1 # number of items to record
out_records = c("Final value for **B**, continous model, _dt_=0.1",
"Final value for **B**, discrete-time model, _dt_=0.1")
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 = "Now we'll explore how model parameters impact outcomes. Set _dt_ to 0.01, change the simulation time, *tfinal*, to 200 days. Also change the bacteria growth rate from 1 to 2 per day. Leave all other settings as before (you can keep running the discrete model, or switch back to ODE only). Run the simulation. Switch back and forth between growth rates of 1 and 2 and examine how bacteria and immune response dynamics change. Then, also try a growth rate of 3."
nrec = 2 # number of items to record
out_records = c("Final value for **B**, continous model, *g* = 2",
"Final value for **B**, continous model, *g* = 3")
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 = tid + 1
tasktext = "In the previous task, you might have been surprised to find that the bacteria growth rate has no impact on the number of bacteria at the final, steady state. This is an indication that even a simple 2-variable model can lead to interesting, and maybe non-intuitive results.
\nFor a model as simple as the one we have here, one can mathematically compute the steady state, i.e., the state at which bacteria and immune response don't change further. To do so, realize that _no change_ means both left hand sides of the ODE model are zero. Equivalently, for the discrete time model, it means $B_{t+dt} = B_t$ and the same for **I**. We'll focus on the ODE model here. With the right side of the equations being zero, the model turns into just 2 algebraic equations, namely $0= g B (1-B/B_{max}) - d_B B - k BI$ and $0 = r B I - d_I I$. You can now solve this such that you end up with two equations of the form $B = XX$ and $I = YY$ where *XX* and *YY* are some combinations of model parameters. Let's do that for the second equation.
\nWe rewrite the equation as $r B I = d_I I$, then divide by **I** and *r* to arrive at $B = d_I/r$. This shows that indeed, the number of bacteria at the steady state does not depend on the growth rate *g*. I'll let you do the same for the first equation to get **I** at steady state (note that at some poin in the solving process, you'll have to insert the steady state value of **B** we just found into the equation). You should end up with $I = (r B_{max} (g-d_B) - d_I g)/(k r B_{max})$. Based on this, we expect that doubling *d~I~* will lead to an increase (doubling) in $B$ at steady state, and a decrease in **I**. Let's test that.
\nSet everything back as in task 1 (you can use the _Reset Inputs_ button), but increase the simulation time to 200 days (to ensure we reach steady state). Run the model. Record **B** and **I** at the end. Then set $d_I=3$, run again, and again record values for **B** and **I**. Do it again for $d_I=4$."
nrec = 6 # number of items to record
out_records = c("Final value for **B**, continous model, _d~I~_ = 2",
"Final value for **B**, continous model, _d~I~_ = 3",
"Final value for **B**, continous model, _d~I~_ = 4",
"Final value for **I**, continous model, _d~I~_ = 2",
"Final value for **I**, continous model, _d~I~_ = 3",
"Final value for **I**, continous model, _d~I~_ = 4")
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 7
#########################
tid = tid + 1
tasktext = "Again, set everything back as in task 1 (you can use the _Reset Inputs_ button). Set the y-axis to log for the plot. Change immune response decay rate, _d~I~_, to 0.2. You should get a plot where both immune response and bacteria numbers drop below 1 and take on fractional values. Contemplate what that means. If not clear, re-read the bullet points at the start of this page. This is one of the draw-backs of ODE based models, and we'll revisit this topic in the stochastic apps."
nrec = 2 # number of items to record
out_records = c("Minimum value of **B** (report as X.YZ)",
"Minimum value of **I** (report as X.YZ)")
out_types = rep("Numeric",nrec)
out_notes = rep("Round to two significant digits, report as non-scientic notation (as X.YZ)",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 = "Go wild! Change any inputs you want to change, see how it affects the results. To start building intuition, it is best to change one quantity at a time. Before you run a simulation with new settings, contemplate what you think might happen. Then see if it does. Always go through an iterative cylce: Think about expectations -> run simulation -> update understanding (__Do Science/Research__)."
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.
The model we are exploring here belongs to a class of well-studied models in ecology known as predator-prey models. If you want to learn more about these kinds of models, see e.g. [@otto11]. The models in those references are described in the context of ecology, but results transfer to within-host situation. For somewhat similar within-host models (applied to TB and antibiotic resistance), see [@antia96; @handel09b].
References
ahgroup/DSAIRM documentation built on Oct. 4, 2023, 6:50 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#load various variable definitions that are the same for each app source('startup_script.R') #source helper functions (defined in the startup script) #needs to happen inside each Rmd file since knitr starts a new environment sapply(files_to_source, source) #load the settings file for the current app #so we can automatically include figure, list the functions in the further information section #and use other information specific to the current app for the task table generation currentrmdfile = knitr::current_input() #currentrmdfile = "basicbacteria_documentation.Rmd" #for debugging appsettings = get_settings(currentrmdfile,appdocdir,packagename)
Overview {#shinytab1}
This app allows exploration of a very basic bacteria (or other extracellular pathogens) infection model, with one compartment for bacteria and one compartment for the immune response. The main goal for this app is to get you familiar with simulation models, as well as the overall setup and ideas behind using these simulations, and how to run them. Read about the model in The Model tab. Then, work through the tasks described in the What To Do tab.
Learning Objectives
- Become familiar with simulation models.
- Learn about steady states and how to compute them.
- Understand how model parameters impact outcomes.
The Model {#shinytab2}
Model Overview
All the models we consider are compartmental models. Compartmental means that we place the agents/players/entities/individuals of interest into distinct compartments/categories. We then only track the total number of individuals in each of these compartments. In this simple model, we track bacteria and some (fairly abstract) notion of immune response strength, using the following notation:
- B - bacteria
- I - immune response
In addition to specifying the compartments or variables of a model, we need to specify the dynamics determining the changes for each compartment. Broadly speaking, there are processes or flows 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 our system, we specify the following processes/flows:
- Bacteria grow/divide at some maximum rate, g, and saturate as they approach some maximum carrying capacity, B~max~.
- Bacteria die at a natural death rate, d~B~.
- Bacteria are killed by the immune response at some rate, k.
- The immune response grows proportional to the number of bacteria and itself at some rate, r.
- The immune response decays at some rate, d~I~.
Model Diagram
A very good way to describe compartmental models and move from a verbal description toward a mathematical/computational formulation is by using diagrams. Being able to go back and forth between verbal description, diagram and mathematical/computational model is a crucial skill when building models. The diagram for a compartmental model is often called flow diagram. It consists of a box for each compartment/variable (here B and I), and arrows pointing in and out of boxes to describe flows and interactions. For the model described above, the flow diagram looks as follows:
knitr::include_graphics(path = paste0("../media/",appsettings$modelfigname))
For the diagrams in this R package, solid arrows indicate physical flows, i.e. movement from a compartment to another (e.g. bacteria moving out of the compartment because of death, or bacteria increasing in the compartment due to growth), while dashed arrows indicate interactions without physical flow (e.g. infected cells killing bacteria). This notation is not universal and it is common in the literature to see no distinction made between these 2 types of flows and only solid arrows being used.
Model Equations
Next, we need to implement this verbal model and flow diagram in a way that it can be run on a computer. While we already specified some of the model features, e.g. that we only track compartments/total numbers and not individual bacteria, both the verbal model and diagram do not yet fully specify the underlying mathematical model. In fact, there are several ways the diagram could be translated into a mathematical/computer model. For this app, we consider two approaches, namely discrete-time deterministic and continuous-time deterministic models. In some of the other apps, you will see another variant, namely continuous-time stochastic model implementations.
Implementation I
The most common way to implement compartmental models is as a continuous-time, deterministic process, formulated as a set of ordinary differential equations (ODEs). Each compartment/variable gets an equation. The right side of each equations specifies the processes going on in the system and how they change the numbers in each compartment via inflows and outflows. For the model described above, the equations look like this:
\begin{align} \dot B &= g B (1-\frac{B}{B_{max}}) - d_B B - k BI \ \dot I &= r B I - d_I I \end{align}
We are using here the short-hand notation where a dot over the variable indicates differentiation with respect to time, and we also do not explicitly indicate that all variables depend on time. This is the most common notation. A more detailed, completely equivalent notation is:
\begin{align} \frac{d}{dt}B(t) &= gB(t)(1-\frac{B(t)}{B_{max}})-d_B B(t) - k B(t)I(t) \ \frac{d}{dt} I(t) &= r B(t) I(t) - d_I I(t) \end{align}
This notation specifies explicitly that B and I are functions of time, and that the differentiation is with respect to time. Since this is always the case in our models, it is not necessary to be that explicit and the more compact notation shown above is the one you will see in this package and it's also the most common throughout the literature.
Implementation II
Continuous time models implemented as ordinary differential equations are the most common types of models. However, other implementations of the above model are possible. One alternative formulation is a discrete-time deterministic equivalent to the ODE model. For such an implementation, the equations are:
\begin{align} B_{t+dt} &= B_t + dt \left( gB_t (1-\frac{B}{B_{max}}) - d_B B_t - k I_t B_t \right) \ I_{t+dt} &= I_t + dt \left( r I_t B_t - d_I I_t \right) \end{align}
In words, the number of bacteria and immune response at a time step dt in the future is given by the number at the current time, t, plus/minus the various growth and death/removal processes. The latter need to be multiplied by the time step, since less of these events can happen if the time step is smaller. As the time-step gets small, this discrete-time model approximates the continuous-time model above. In fact, when we implement a continuous-time model on a computer, the underlying simulator runs a smart version of a discrete-time model and makes sure the steps taken are so small that the numerical simulation is a good approximation of the continuous-time model. If you want to learn more about that, you can check out the deSolve R package documentation, which we use to run our simulations.
Model Concepts
In general, the entities that change in our model (i.e. here B and I) are called variables; they are variable and change during the simulation. To run a simulation, we need to specify the starting values for each variable. Those are often called initial conditions.
In contrast, the quantities that are usually fixed for a given scenario are called parameters. For this model, those are g, d~B~, k, r and d~I~. Values for parameters are chosen to match the known biology of a specific disease/scenario we want to model. All parameters need to be in the same time units, e.g. per day or per hour.
You might find for some parameter settings that numbers for bacteria/immune response/virus/cells/etc. drop below 1 (and possibly rebound later). You'll see that in one of the tasks below. If the units of our model are assumed to be in number of pathogen/cells, this of course is biologically not reasonable. Numbers less than 1 are an artifact of the underlying differential equation model. All ODE models (and the discrete time model we consider here) have that problem. A hack-y solution is to monitor the simulation and if a quantity drops below 1, set it to 0 or stop the simulation. A cleaner solution is to treat all entities (virus, cell) as discrete units and allow them to only change in integer steps (in a probabilistic manner). This approach will be discussed in the apps that use stochastic models. For most apps, neither approach is taken, so you'll see numbers less than 1. That's okay for the purpose we use the apps here. Just be careful if you use these kinds of models in your research you might need to pay attention to this issue.
Notes
-
This model is of course very simple. However, as you will see, even with just 2 entities and their interactions, we can get interesting dynamics. This model could be extended by introducing further compartments. For instance we could include separate compartments for the innate and adaptive immune response. You will encounter more detailed models in some of the other apps.
-
This and most other simulations/apps do not have natural time units (unless specifically stated). You could, therefore, assume that your model runs in units of days or minutes/hours/weeks/months..., based on what's most suitable for the system you want to study. You have to make sure that all your parameters are in the right time units. Always make sure to check if a given simulation can handle different time units or assumes specific ones. In general, I specify for each What To Do section what kind of units I have in mind for the different tasks.
-
Some of the tasks in the What To Do section (and in future apps) are fairly open-ended. The idea is that these tasks give you something to get started, but you should feel free to explore the simulations any way you want. Play with them, query them, go through iterations of thinking about what you expect, observing it, and, if discrepancies occur, figuring out why. Essentially, the best way to use these apps is to do your own science/research with them.
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 cell lives for a week, you need to convert it to 7 days, then take the inverse if you want a death rate, which is what usually goes into the model).
#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 100 initial bacteria and an initial level of 1 for the immune response. Assume bacteria grow at a rate of 1 per day (i.e., g = 1), the carrying capacity is 10^5^ and they live for about 2 days (remember that the inverse of the lifespan is the rate of death). Assume the immune response is activated and grows at a rate of 10^-4^, kills bacteria at a rate of 10^-4^ and decays at a rate of 2 per day. Set simulation start time to 0 and final time of 100 (which we assume to be days). Set the time step to 0.01. As you'll see, if we run an ODE model, the time step is only relevant for plotting, not for the underlying model simulation (while this is not the case for the discrete time model). Only run the continuous time ODE model. Plot both x- and y-axes on a linear scale (i.e, no log scales for now). You can stick with ggplot as the plot engine. Run the simulation, see what you get. You should see some oscillations and then the system settles down, with bacteria and immune response at the end of the simulation at around 19996 and 3003, respectively. Change the time step to 0.05, re-run the simulation. Record the number of bacteria at the end of the simulation (Hint: not much changes)." nrec = 1 # number of items to record out_records = c("Number of bacteria at end of simulation, _dt_=0.05") 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 = "Play around with the plot settings. Switch the plotting to have x-axis, y-axis or both plotted on a log scale. Also try plotly as the plot engine. Leave all other settings as before. You should see that while the look of the plot changes, the underlying numbers do not. This is something to be aware of when you see plots in papers or produce your own. The best plot to use is the one that shows results of interest in the clearest form. Usually, the x-axis is linear and the y-axis is either linear or logarithmic. One nice feature about plotly is that the plot is interactive and you can read off numbers. Use this feature to determine the day at which the bacteria and immune response have their second peak. (*Hint: day 32 for bacteria, a bit later for the immune response.*)" nrec = 1 # number of items to record out_records = c("Day of second peak for the immune response") 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 = "Go back to both linear scales for plotting. Keep all variable and parameter settings as before. Set _Models to run_ to _both_. This runs and shows both the continuous-time and discrete-time models. Start with a time step of 0.01. Run the simulation. You should see the results from the 2 models essentially on top of each other and barely distinguishable. You should see that the maximum number of bacteria is 49,513 for the ODE model, and very close for the discrete-time model." nrec = 1 # number of items to record out_records = c("Maximum number of bacteria for the discrete-time 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 4 ######################### tid = tid + 1 tasktext = "Now try different values for the time step, _dt_. Leave all other settings as before. You should notice that as _dt_ gets larger, the continuous-time model results remain the same, but the discrete-time models change and start moving away from those of the continous-time model. At a time step above 0.1, the results start to look very different. Somewhere above a time step of 0.5, it becomes so large that for these parameter settings, the simulation 'crashes' and you get an error message." nrec = 1 # number of items to record out_records = c("Final value for **B**, continous model, _dt_=0.1", "Final value for **B**, discrete-time model, _dt_=0.1") 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 = "Now we'll explore how model parameters impact outcomes. Set _dt_ to 0.01, change the simulation time, *tfinal*, to 200 days. Also change the bacteria growth rate from 1 to 2 per day. Leave all other settings as before (you can keep running the discrete model, or switch back to ODE only). Run the simulation. Switch back and forth between growth rates of 1 and 2 and examine how bacteria and immune response dynamics change. Then, also try a growth rate of 3." nrec = 2 # number of items to record out_records = c("Final value for **B**, continous model, *g* = 2", "Final value for **B**, continous model, *g* = 3") 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 = tid + 1 tasktext = "In the previous task, you might have been surprised to find that the bacteria growth rate has no impact on the number of bacteria at the final, steady state. This is an indication that even a simple 2-variable model can lead to interesting, and maybe non-intuitive results. \nFor a model as simple as the one we have here, one can mathematically compute the steady state, i.e., the state at which bacteria and immune response don't change further. To do so, realize that _no change_ means both left hand sides of the ODE model are zero. Equivalently, for the discrete time model, it means $B_{t+dt} = B_t$ and the same for **I**. We'll focus on the ODE model here. With the right side of the equations being zero, the model turns into just 2 algebraic equations, namely $0= g B (1-B/B_{max}) - d_B B - k BI$ and $0 = r B I - d_I I$. You can now solve this such that you end up with two equations of the form $B = XX$ and $I = YY$ where *XX* and *YY* are some combinations of model parameters. Let's do that for the second equation. \nWe rewrite the equation as $r B I = d_I I$, then divide by **I** and *r* to arrive at $B = d_I/r$. This shows that indeed, the number of bacteria at the steady state does not depend on the growth rate *g*. I'll let you do the same for the first equation to get **I** at steady state (note that at some poin in the solving process, you'll have to insert the steady state value of **B** we just found into the equation). You should end up with $I = (r B_{max} (g-d_B) - d_I g)/(k r B_{max})$. Based on this, we expect that doubling *d~I~* will lead to an increase (doubling) in $B$ at steady state, and a decrease in **I**. Let's test that. \nSet everything back as in task 1 (you can use the _Reset Inputs_ button), but increase the simulation time to 200 days (to ensure we reach steady state). Run the model. Record **B** and **I** at the end. Then set $d_I=3$, run again, and again record values for **B** and **I**. Do it again for $d_I=4$." nrec = 6 # number of items to record out_records = c("Final value for **B**, continous model, _d~I~_ = 2", "Final value for **B**, continous model, _d~I~_ = 3", "Final value for **B**, continous model, _d~I~_ = 4", "Final value for **I**, continous model, _d~I~_ = 2", "Final value for **I**, continous model, _d~I~_ = 3", "Final value for **I**, continous model, _d~I~_ = 4") 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 7 ######################### tid = tid + 1 tasktext = "Again, set everything back as in task 1 (you can use the _Reset Inputs_ button). Set the y-axis to log for the plot. Change immune response decay rate, _d~I~_, to 0.2. You should get a plot where both immune response and bacteria numbers drop below 1 and take on fractional values. Contemplate what that means. If not clear, re-read the bullet points at the start of this page. This is one of the draw-backs of ODE based models, and we'll revisit this topic in the stochastic apps." nrec = 2 # number of items to record out_records = c("Minimum value of **B** (report as X.YZ)", "Minimum value of **I** (report as X.YZ)") out_types = rep("Numeric",nrec) out_notes = rep("Round to two significant digits, report as non-scientic notation (as X.YZ)",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 = "Go wild! Change any inputs you want to change, see how it affects the results. To start building intuition, it is best to change one quantity at a time. Before you run a simulation with new settings, contemplate what you think might happen. Then see if it does. Always go through an iterative cylce: Think about expectations -> run simulation -> update understanding (__Do Science/Research__)." 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.
The model we are exploring here belongs to a class of well-studied models in ecology known as predator-prey models. If you want to learn more about these kinds of models, see e.g. [@otto11]. The models in those references are described in the context of ecology, but results transfer to within-host situation. For somewhat similar within-host models (applied to TB and antibiotic resistance), see [@antia96; @handel09b].
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.