DSAIRM Package Tutorial
2023-08-23
Package description
Dynamical Systems Approaches to Immune Response Modeling
(DSAIRM) is an R package that allows you to explore and learn about
dynamical systems modeling of within-host infection and immune response
dynamics.
The package consists of several simulations/apps that allow for the
simulation and exploration of different topics in within-host infection
and immune response dynamics.
The underlying models are written as compartmental dynamical models,
either deterministic using differential equations (deSolve package) or
stochastic using a Gillespie-type approach (adaptivetau package). A
graphical user interface is wrapped around each simulation/app. The
graphical user interface is written using the functionality of the R
Shiny package. This allows exploration of models and concepts without
the need to write any code. At the same time, the package is structured
in a way that should allow those interested in the actual models and
learning R coding to easily move from one stage to another.
Each app is meant to be fully self-explanatory and contains a
description of the model, a list of tasks the user could try, and
further information on background reading and other resources.
The main way of interaction with the models is through a graphical
user interface. It is not necessary - though possible - to read or write
computer code. The package also provides a path toward moving from the
graphical user interface to direct interaction with pre-written
simulations all the way to modifying the underlying simulations to fit
your needs. These different use cases for the package are described
below.
Intended audience and goal of the package
The main audience are individuals who want to learn how mechanistic
simulation models can be used to study within-host infection dynamics.
While the ideal use of this package is likely as component of a class on
this topic, self-learning might also be possible. Each app contains a
fair amount of documentation and description of the model. By reading
the model description, doing the tasks and probably also reading some of
the papers in the future information section, it might be possible to
teach yourself the beginnings of within-host modeling on your own.
Installation
I assume you have R
installed. I also highly recommend RStudio, though it
is not required. After having opened R/Rstudio, you can install the
package from CRAN by typing the following into the R console (the bottom
left window if you are using RStudio):
install.packages('DSAIRM')
Alternatively, you can install the latest development version (with
usually both more features and more bugs) from github, using the
remotes package. The following commands, typed into the
R console, will get you up and running:
install.packages('remotes') #if not already installed
remotes::install_github('ahgroup/DSAIRM')
The package depends on other packages, which will also be installed
as needed. Package installation is a one-time process, unless you update
or reinstall R itself, or want to install a newer version of DSAIRM.
Using the package through the graphical user interface (Level
1)
The interactive exploration of the models and infectious disease
concepts through the graphical user interface is the main intended use
of this package. The idea is that you start at Level 1, and then
depending on needs and interests, you can decide to move on to the next
levels, described below.
After install, load the package by typing the following into the
R console.
You need to load the package every time you start a new R/Rstudio
session.
If the package loads correctly, you should receive a short greeting.
Now you can open the main menu by typing the following into the
R console:
If successful, you should see a graphical menu open in your browser
that looks something like this:
From the main menu, you can choose the different apps corresponding
to different modeling topics and scenarios. Each app has the inputs you
can set on the left and the outputs (figures and text) on the right.
This is an example of one of the apps:
As seen in the figure, below each app are four tabs with different
sections providing information and instruction. The Overview
and Further Information tabs provide a brief introduction and
pointers to further exploration options, respectively.
The two most important tabs are the Model and What to
do tabs. The Model tab gives a detailed description of the
model underlying the app, shows the model diagram and equations, and
provides any other further information needed to understand what is
going on. The following figure shows an example of this for one of the
apps.
The What to do tab contains a list of suggested tasks you
might want to go through. By doing those tasks, you will learn the topic
a specific app is meant to teach. The following figure shows an example
for part of one of those What to do sections.
Once you finished exploring a specific app, you can return to the
main menu. Eventually, you can exit back to the R console
by clicking the Exit button. You will not need to read or write
any code at this level of exploration and learning. If an when you want
to learn more, you can proceed to level 2.
Directly interacting with the simulation functions (Level 2)
Once you are comfortable interacting with the models and have a good
understanding of the concepts covered by the different apps, it is
possible, without too much effort, to interact with the simulation
models directly. This provides more flexibility at the cost of writing a
limited amount of code.
To facilitate direct interaction and modification of the underlying
simulations, each app is structured in such a way that the underlying
model/simulation is a stand-alone function. For some apps, there are
multiple underlying functions involved. You can call any of these
functions directly, without going through the graphical interface. The
Further Information tab inside each app provides the name of
the corresponding underlying function(s) and how to call their help
files to get information on their use.
Consider as example the first app, called “Basic Bacteria Model”.
This model has 2 underlying simulator functions, one that runs the
discrete-time model called simulate_basicbacteria_discrete
and one that runs the continuous, differential equation model called
simulate_basicbacteria_ode. Assume you want to learn more
about the latter model.
After loading the package (if not already loaded) with
you can learn about the inputs and outputs of the model/function by
looking at its documentation
help('simulate_basicbacteria_ode')
The help file explains that one can run the simulation by specifying
initial number of bacteria and immune response strength, the different
model parameters, as well as some time values. Each parameter has some
default, listed in the help file. Running a simulation model without any
specifications runs it with the defaults:
result <- simulate_basicbacteria_ode()
Calling the simulation function execuctes the underlying dynamical
model (which is described in the Model section of the app).
Most simulation functions produce and return time-series for the
dynamics of each of the variables that are tracked. Results are returned
in a list structure, if a time-series exists, it is generally stored as
ts. The Value section of the help file for each
function explains in detail what the function returns.
You can take the results returned from the simulation function and
further process them. For instance you could plot the time-series for
bacteria as function of time:
plot(result$ts[ , "time"],result$ts[ , "B"],xlab='Time',ylab='Bacteria',type='l')
You can overwrite the default settings by providing your own inputs
for specific parameters, as shown in the next example, which sets the
rate at which bacteria grow, g, and the rate at which the
immune response is induced and grows, r, to values different
than the defaults, while using the defaults for the remainder (this is
equivalent to setting different inputs through the graphical interface
in level 1):
result <- simulate_basicbacteria_ode(g = 0.5, r = 0.002)
plot(result$ts[ , "time"],result$ts[ , "B"],xlab='Time',ylab='Bacteria',type='l')
Any values that you do not specify will be kept at their
defaults.
Note that unless explicitly specified, the models do not have
inherent time units. Instead, those depend on the time units you choose
for the parameters. It is important to ensure that all quantities
(parameters and time settings) have the same time units, e.g. days or
months (or the inverse of those units for the rate parameters).
Not all simulator functions return time-series. For instance the
function underlying the uncertainty and sensitivity app returns specific
outcomes (maximum and final state of variables) for different parameter
values. In every instance, the documentation for the function explains
what is returned, so you know how to take the returned results and
further process them.
The ability to call the simulation functions without using the
graphical interface allows for more flexibilty in exploring the models.
For instance if you wanted to explore the behavior of a model
systematically for different values of a given parameter, this would
need to be done manually through the graphical interface. Calling the
function directly allows you automate this. As an illustration, assume
you wanted to know how the peak bacteria load changes with changes in
the immune response activation rate, r, (while keeping all
other parameters at some specified fixed values). The following lines of
code show how this can be achieved. A loop is run over different
r values, for each r value the simulation is run and
the peak of the bacteria is recorded. At the end, peak of bacteria as
function of immune activation rate is plotted:
#values for immune activation rate, r, for which to run the simulation
rvec = 10^seq(-5,-2,length=20)
#this variable will hold the peak values for each r
Bpeak = rep(0,length(rvec))
for (n in 1:length(rvec))
{
#run the simulation model for different values of r
#for simplicity, all other inputs are kept at their defaults
result <- simulate_basicbacteria_ode(r = rvec[n])
#record max number of bacteria for each value of r
Bpeak[n] <- max(result$ts[,"B"])
}
#plot final result
plot(rvec,Bpeak,type='p',xlab='Immune activation rate',ylab='Peak bacteria',log='xy')
Thus, you can add your own custom code to the existing simulator
functions and with a few lines of extra code analyze and explore many
more questions and scenarios than those accessible through the graphical
interface. This provides a lot more flexibility, but requires writing
some R code to interface with the supplied simulator functions.
Modifying the simulation functions (Level 3)
While level 2 provides you with a fair amount of flexibility and the
ability to investigate the models and ask questions that would be hard
or impossible going through the graphical interface, there is still one
major constraint. Namely this approach only allows use and analysis of
pre-written, existing simulation models. This constraint can be overcome
by directly modifiying and customizing the underlying simulation
functions. As was true for moving from level 1 to 2, this move from
level 2 to what I call level 3 provides further (almost unlimited)
flexibility at the cost of having to write increasingly more
R code.
To make modifying the existing functions easy, copies of all
simulator functions are provided in a folder called
simulatorfunctions which can be downloaded as a zip file from
the main menu. Each function in that folder starts with
simulate_. The Further Information section in each
app tells you which simulator function goes with which app.
The code for each simulator function is (hopefully) well documented.
Some basic to intermediate level of R coding experience is
required to successfully modify the functions. In addition to modifying
the simulator function of interest, you will also have to write code to
interact with your modified function(s) (as described in Level
2).
To provide an example of this approach, assume you have a system and
question for which the basic bacteria ODE model tracking only bacteria
and immune response is reasonable, however the DSAIRM model is not fully
suitable. Instead, for your system, you need a model in which the immune
system growth term includes saturation at some maximum rate when
bacteria numbers are high. This can be accomplished by a change of the
term rBI to rBI/(B+s). (See the documentation for this
app for an explanation of each model term and corresponding equation).
This change leads to a growth at rate rB proportional to the
number of bacteria if B is smaller than some threshold
parameter s, and turns into a growth at fixed maximum rate
r, independent of bacteria numbers, if B is larger
than s.
To implement this, get the code for this model, which is
simulate_basicbacteria_ode.R. (To make things easy, the
name of a .R file containing the code and the name of the
function itself are the same.) After finding the file, make a copy and
rename it (we call it mysimulator.R here).
Then open the file (ideally in RStudio) and make the following
modifications to the code:
old:
simulate_basicbacteria_ode <- function(B = 10, I = 1, g = 1, Bmax = 1e+05, dB = 0.1, k = 1e-06, r = 0.001, dI = 1, tstart = 0, tfinal = 30, dt = 0.05)
new:
mysimulator <- function(B = 10, I = 1, g = 1, Bmax = 1e+05, dB = 0.1, k = 1e-06, r=1e3, dI=1, tstart = 0, tfinal = 30, dt = 0.05, s=1E3)
Note the changed default value for r to ensure the immune
response is sufficiently triggered.
old:
pars = c(g=g,Bmax=Bmax,dB=dB,k=k,r=r,dI=dI)
new:
pars = c(g=g,Bmax=Bmax,dB=dB,k=k,r=r,dI=dI,s=s)
old:
dBdt = g*B*(1-B/Bmax) - dB*B - k*B*I
dIdt = r*B*I - dI*I
new:
dBdt = g*B*(1-B/Bmax) - dB*B - k*B*I
dIdt = r*B*I/(s+B) - dI*I
With these changes made, it is now possible to use the new
mysimulator function to ask specific questions regarding
this model, and thus the underlying biological system of interest for
which this model is a - hopefully useful - approximation. For instance,
you can explore how different values of the saturation parameter
s impact the maximum level of the immune response. This
requires a slight modification of the code shown above in Level
2 as follows:
#initialize the new function
#it needs to be in the same directory as this code
source('mysimulator.R')
#values of saturation parameter to explore
svec = 10^seq(-3,3,length=20)
#this will record the maximum immune response level
Ipeak = rep(0,length(svec))
for (n in 1:length(svec))
{
#run the simulation model for different values of s
#for simplicity, all other inputs are kept at their defaults
result <- mysimulator(s = svec[n])
#record max immune response for each value of s
Ipeak[n] <- max(result$ts[,"I"])
}
A plot showing the results is:
plot(svec,Ipeak,type='p',xlab='Saturation parameter, s',ylab='Peak immune response',log='xy')
Using one of the provided simulator functions as starting point and
modifying it is likely easier than having to write a new model
completely from scratch. Eventually, as you gain more coding experience,
you have (almost) unlimited flexibility and power regarding the models
you can create and questions you can address, of course at the cost of
having to write increasingly more R code. You are
essentially only limited by what can be accomplished using the
R programming language and your ability and interest in
writing code.
Contributing to the package
The package is on GitHub and you can use the usual GitHub process to
file bug reports, send feature requests, contribute updates and bug
fixes, etc. If you have any comments or feedback, I very much welcome
them. Please file a GitHub issue and let
me know what you think.
The package is built in a way that makes it (hopefully) easy for
others to contribute new simulations/apps. To that end, the package
contains this
Markdown file which provides further information on the details of
the package structure. If you plan to develop new apps, or make other
substantial contributions, it’s best to get in touch with me first via
email or GitHub.
Contributors
This R package is developed and maintained by Andreas Handel. A full list of
contributors and a Bibtex entry for the citation can be found
here.
References
