A generic library for simulating networks of ordinary and delay differential equations for systems modeling.
Contributors to this file should be aware of Adam Pritchard's Markdown Cheatsheet.
- System Requirements 0. Operating System
- Tutorial with Simple Model
- Compilation and Model Building 0. Species and Reactions 0. Declaring Species 1. Declaring Reactions 2. Defining Reaction Rate Formulas 3. Defining Reaction Inputs and Outputs
- Running the Simulation 0. Description of the Simulation 0. Preamble 1. Deterministic 2. Stochastic
- Authorship and License
linux, mac. windows if-y
CMake version 2.8+ required, along with a supported build manager (make is the only tested system currently.) A C++ compiler with support for at least the C++11 standard is required. In order to compile the CUDA accelerated code, both a CUDA 6.0+ compiler and NVIDIA GPU hardware with "Compute Capability 3.0+" are needed.
Create a sub-directory of your main directory called lemonade-build and navigate to it:
mkdir lemonade-build
cd lemonade-build
Next, make five new empty files within the lemonade-build directory: $ touch specie_list.hpp reactions_list.hpp reaction_deltas.hpp model_impl.hpp param_set.csv
Open specie_list.hpp with your text editor of choice and add the following code:
#ifndef CRITICAL_SPECIE
#define CRITICAL_SPECIE SPECIE
#define UNDO_CRITICAL_SPECIE_DEF
#endif
CRITICAL_SPECIE(lemon)
CRITICAL_SPECIE(sugar)
SPECIE(lemonade)
SPECIE(money)
SPECIE(good_lemonade)
SPECIE(customer_satisfaction)
#ifdef UNDO_CRITICAL_SPECIE_DEF
#undef CRITICAL_SPECIE
#undef UNDO_CRITICAL_SPECIE_DEF
#endif
Save and exit from this file.
Open reactions_list.hpp and add the following code:
#ifndef DELAY_REACTION
#define DELAY_REACTION REACTION
#define UNDO_DELAY_REACTION_DEF
#endif
REACTION(make_lemon)
REACTION(make_sugar)
REACTION(make_lemonade)
REACTION(customer_purchase)
REACTION(buy_lemons)
REACTION(buy_sugar)
REACTION(make_quality)
REACTION(quality_to_satisfaction)
REACTION(satisfaction_to_tips)
#ifdef UNDO_DELAY_REACTION_DEF
#undef DELAY_REACTION
#undef UNDO_DELAY_REACTION_DEF
#endif
Save and exit from this file.
Open reaction_deltas.hpp and add the following code:
#include "utility/common_utils.hpp"
#include "core/reaction.hpp"
#include <cstddef>
STATIC_VAR int num_deltas_make_lemonade = 3;
STATIC_VAR int deltas_make_lemonade[] = {-2,-3,1};
STATIC_VAR specie_id delta_ids_make_lemonade[] = {lemon,sugar,lemonade};
STATIC_VAR int num_deltas_make_lemon = 1;
STATIC_VAR int deltas_make_lemon[] = {1};
STATIC_VAR specie_id delta_ids_make_lemon[] = {lemon};
STATIC_VAR int num_deltas_make_sugar = 1;
STATIC_VAR int deltas_make_sugar[] = {1};
STATIC_VAR specie_id delta_ids_make_sugar[] = {sugar};
STATIC_VAR int num_deltas_customer_purchase = 2;
STATIC_VAR int deltas_customer_purchase[] = {-1,2};
STATIC_VAR specie_id delta_ids_customer_purchase[] = {lemonade,money};
STATIC_VAR int num_deltas_buy_lemons = 2;
STATIC_VAR int deltas_buy_lemons[] = {-3,2};
STATIC_VAR specie_id delta_ids_buy_lemons[] = {money,lemon};
STATIC_VAR int num_deltas_buy_sugar = 2;
STATIC_VAR int deltas_buy_sugar[] = {-2,25};
STATIC_VAR specie_id delta_ids_buy_sugar[] = {money,sugar};
STATIC_VAR int num_deltas_make_quality = 2;
STATIC_VAR int deltas_make_quality[] = {-4,1};
STATIC_VAR specie_id delta_ids_make_quality[] = {lemonade, good_lemonade};
STATIC_VAR int num_deltas_quality_to_satisfaction = 3;
STATIC_VAR int deltas_quality_to_satisfaction[] = {-2,-1,1};
STATIC_VAR specie_id delta_ids_quality_to_satisfaction[] = {good_lemonade, lemonade, customer_satisfaction};
STATIC_VAR int num_deltas_satisfaction_to_tips = 2;
STATIC_VAR int deltas_satisfaction_to_tips[] = {-1,4};
STATIC_VAR specie_id delta_ids_satisfaction_to_tips[] = {customer_satisfaction, money};
Save and exit from this file.
Open model_impl.hpp and add the following code:
#ifndef MODEL_IMPL_H
#define MODEL_IMPL_H
#include "core/reaction.hpp"
#include "core/specie.hpp"
#include "core/model.hpp"
#include "sim/base.hpp"
#include <cstddef>
template<>
template<class Ctxt>
RATETYPE reaction<make_lemonade>::active_rate(const Ctxt& c) {
return c.getRate(make_lemonade) * c.getCon(lemon)*c.getCon(sugar);
}
template<>
template<class Ctxt>
RATETYPE reaction<make_lemon>::active_rate(const Ctxt& c) {
return (c.getRate(make_lemon) * (c.getCon(lemon) - (c.getCon(lemon) *c.getCritVal(rcrit_lemon))));
}
template<>
template<class Ctxt>
RATETYPE reaction<make_sugar>::active_rate(const Ctxt& c) {
return (c.getRate(make_sugar) * (c.getCon(sugar) - c.getCritVal(rcrit_sugar)));
}
template<>
template<class Ctxt>
RATETYPE reaction<customer_purchase>::active_rate(const Ctxt& c) {
return c.getRate(customer_purchase);
}
template<>
template<class Ctxt>
RATETYPE reaction<buy_lemons>::active_rate(const Ctxt& c) {
return c.getRate(buy_lemons);
}
template<>
template<class Ctxt>
Real reaction<buy_sugar>::active_rate (const Ctxt& c) {
return c.getRate(buy_sugar);
}
template<>
template<class Ctxt>
Real reaction<make_quality>::active_rate (const Ctxt& c) {
return c.getRate(make_quality)*c.getCon(lemon)*c.getCon(sugar);
}
template<>
template<class Ctxt>
Real reaction<quality_to_satisfaction>::active_rate (const Ctxt& c) {
return c.getRate(quality_to_satisfaction);
}
template<>
template<class Ctxt>
Real reaction<satisfaction_to_tips>::active_rate (const Ctxt& c) {
return c.getRate(satisfaction_to_tips);
}
#endif // MODEL_IMPL_H
Save and exit from this file.
Open param_set.csv and add the following code:
make_lemonade,customer_purchase,buy_lemons,buy_sugar,make_lemon,make_sugar,make_quality,quality_to_satisfaction,satisfaction_to_tips,rcrit_lemon,rcrit_sugar,
2.25,0.75,1.655,1.237,4.12,12.12,0.657,4.3,2.7,12,19,
Save and exit from this file.
Reopen terminal, and navigate to the lemonade_build directory. Then, run the following commands:
cmake <path_to_DENSE>/DENSE && make;
./simulation -p ./param_set.csv -c 100 -w 100 -t 10 -u 0.5
This will run a stochastic simulation of a lemonade stand for 10 minutes, fetching data every 30 seconds. To run the same simulation with deterministic algorithms, add the -s 0.5 flag to the above. To print the output to a seperate file, add -e "output.csv".
Declare species in specie_list.hpp. List the specie names between the two sets of C++ macros (the lines that begin with #) in the same format as below. The following example lists two normal species, alpha and charlie, and one critical specie, bravo.
A critical specie is a specie whose effect on a particular reaction is bounded by a threshhold. Imagine a specie that doesn't have an effect on a reaction rate until it reaches a certain concentration or might stop having a linear relationship with the rate after a certain conc level. Whether this critical value is an upper or lower bound is decided in the rate equations in model_impl.hpp but the value itself is inputted alongside delays and initial rates (see 3.1) in the parameter set. To establish the pairing of this specie and its critical value, it is declared as a critical specie.
SPECIE(alpha)
CRITICAL_SPECIE(bravo)
SPECIE(charlie)
Declare reactions in reactions_list.hpp. List the reaction names between the two sets of C++ macros (the lines that begin with #) in the same format as below. The following example lists one delay reaction, alpha_degredation, and three normal reactions, bravo_synthesis, alpha_synthesis, and bravo_degredation. While this particular reaction naming scheme is not required, it can be helpful.
REACTION(alpha_synthesis)
REACTION(bravo_synthesis)
DELAY_REACTION(alpha_degredation)
REACTION(bravo_degredation)
Define all of reaction rate functions in model_impl.hpp.
For example, if a reaction is enumerated alpha_synthesis, it should be declared as a
function like this:
RATETYPE reaction<alpha_synthesis>::active_rate(const Ctxt& c) const { return 6.0; }
Or, for a more interesting reaction rate, you might do something like:
RATETYPE reaction<bravo_degredation>::active_rate(const Ctxt& c) const {
return c.getRate(bravo_degredation) * c.getCon(bravo) * c.neighbors.calculateNeighborAvg(charlie);
}
Refer to the Context API in the following section for instructions on how to get delays and critical values for more complex reaction rate functions.
Contexts are iterators over the concentration levels of all species in all cells. Use them to get the conc values of specific species that factor into reaction rate equations.
To get the concentration of a specie where c is the context object and bravo is the specie's enumeration:
c.getCon(bravo)
To get the delay time of a particular delay reaction that is enumerated as alpha_degredation and is properly identified as a delay reaction in reactions_list.hpp (see 2.0.1: Declaring Reactions):
RATETYPE delay_time = c.getDelay(dreact_alpha_degredation);
To get the past concentration of alpha where delay_time, as specified in the previous example, is the delay time for R_ONE:
c.getCon(alpha, delay_time);
To get average concentration of charlie in that cell and its surrounding cells:
c.calculateNeighborAvg(charlie)
To get the past average concentration of SPECIE in that cell and its surround cells:
c.calculateNeighborAvg(charlie, delay_time)
To get the critical value of a critical specie that is enumerated as bravo and is properly identified as a critical specie in specie_list.hpp (see 2.0.0: Declaring Species):
c.getCritVal(rcrit_bravo)
Define each reaction's reactants and products in reaction_deltas.hpp.
Say a reaction enumerated as R_ONE has the following chemical formula:
2A + B --> C
The proper way to define that reaction's state change vector is as follows:
STATIC_VAR int num_deltas_R_ONE = 3;
STATIC_VAR int deltas_R_ONE[] = {-2, -1, 1};
STATIC_VAR specie_id delta_ids_R_ONE[] = {A, B, C};
Running make after having initialized CMake in the desired directory will automatically run csv_gen as the simulation is being compiled. csv_gen will generate *_template.csv files formatted for the directory's particular model. The easiest way to fill these out is with an Excel-like program such as LibreOffice Calc. Remember to always save changes using the original *.csv file extension. Changes should also be saved in a file name different from the one automatically generated so that there is no chance csv_gen will overwrite your settings.
At its core, CSV files contain numerical values seperated by commas. Listed below are three categories of characters/strings that the simulation's CSV parser DOES NOT parse.
-
Empty cells, blank rows, and whitespace
To illustrate, the following two examples are equivalent.
3.14, , 2001, -2.18, 41, 2.22e-223.14,2001,-2.18,41,2.22e-22 -
Comments, i.e. rows that begin with a
## I am a comment! Below is the data. 9182, 667 -
Any cell that contains a character which is not a number,
.,+,-, ore.Only the following scientific notation is supported:
-0.314e+1, 3.00e+8, 6.63e-34These would be considered invalid:
3.33*10^4, 1.3E-12, 4.4x10^2Often times cells which do not contain numbers are intended to be column headers. These are not parsed by the simulation, and can technically be modified by users as they wish.
flying_pig_synthesis, plutonium_synthesis, cultural_degredation, 21.12021, 33, 101.123,It is futile, however, to add/remove/modify the column headers with the expectation of changing the program's behavior. Data must be entered in the default order if it is to be parsed properly.
The parameter set template is named param_sets_template.csv by default. Parameter set files can contain more than one set per file (each being on their own line). When a file is loaded into the simulation, all sets are initialized and executed in parallel.
Below is an example of a parameter sets file that contains three sets:
alpha_synthesis, alpha_degredation,
332, 101,
301, 120,
9.99e+99, 1.0e-99,
The perturbations template is named param_pert_template.csv by default. Perturbation files should only contain one set of perturbations. Only this one perturbations set is applied to all parameter sets when a simulation set is being run.
Use 0 to indicate that a reaction should not have perturbations. In the example below, alpha_synthesis has a perturbation while alpha_degredation does not.
alpha_synthesis, alpha_degredation,
0.05, 0,
The gradients template is named param_grad_template.csv by default. Gradient files should only contain one set of gradients. Only this one gradients set is applied to all parameter sets when a simulation set is being run.
Use 0 under all four columns of a reaction to indicate that it should not have a gradient. In the first example, alpha_synthesis does not have a gradient, while in the second example, alpha_degredation does.
alpha_synthesis_x1, alpha_synthesis_y1, alpha_synthesis_x2, alpha_synthesis_y2,
0, 0, 0, 0,
alpha_synthesis_x1, alpha_synthesis_y1, alpha_synthesis_x2, alpha_synthesis_y2,
2, 0.8, 5, 2.52,
alpha_degredation's gradient is between cell columns 2 and 5, with a multiplier of 0.8 starting at column 2, linearly increasing to 2.52 by column 5.
Gradient Suffixes Chart
| Suffix | Meaning |
|---|---|
| x1 | start column |
| y1 | start multiplier |
| x2 | end column |
| y2 | end multiplier |
add instructions here on analyses.xml
The application supports simulation of well-stirred chemical systems in a single enclosed environment and in multicellular networks. The simulation uses delay-differential equations to advance a given model and estimate concentrations for given species updated by given reactions. The size of dt varies depending on which simulation algorithm is run.
The Deterministic Simulation Algorithm uses rate reaction equations to approximate the concentration levels of every specie over a constant, user-specified, time step. Files concerning this simulation type are found in source/sim/determ. The user can turn on deterministic simulation by specifying a time step in the command line (see 3.1.2).
The Stochastic Simulation Algorithm loosely follows Dan Gillespie's tau-leaping process. The dt is calculated from a random variable as the time until the next reaction event occurs. Molecular populations are treated as whole numbers and results are non-deterministic unless a random seed is provided by the user in the command line (see 3.1.2). The algorithm is much more performance intensive than the deterministic algorithm and is most ideal for smaller tissue sizes and shorter simulation durations.
After the simulation has been compiled, the only file required to run a deterministic or stochastic simulation is a filled out param_sets_template.csv. It is suggested that this file be renamed to param_sets.csv upon completion.
In order to take advantage of perturbations and gradients, param_pert_template.csv and param_grad_template.csv need to be filled out. Rename these to param_pert.csv and param_grad.csv or something similar upon completion.
Similarly, analyses must also be configured in a seperate file. A blank analyses.xml should be included in the /template_model directory.
The table below can also be accessed by running simulation either without any command line arguments or with any of the following flags: -h, --help, --usage.
Short and long flags are equivalent; either can be used to get the same program behavior. Short flags must be preceeded by - while long flags must be preceeded by --. As examples: -h and --help. Arguments that have a parameter require additional text to proceed the flag itself like so: -p ../param_sets.csv, -b 0.05, and -o "alpha, bravo".
RATETYPE is set to double (double-precision floating-point) by default and can be changed in source/util/common_utils.hpp.
| Short | Long | Parameter | Description |
|---|---|---|---|
h |
help or usage |
none | Print information about all command line arguments. |
n |
no-color |
none | Disable color in the terminal. |
p |
param-sets |
string |
Relative file location and name of the parameter sets *.csv. ../param_sets.csv, for example. |
g |
gradients |
string |
Enables gradients and specifies the relative file location and name of the gradients *.csv. ../param_grad.csv, for example. |
b |
perturbations |
string |
Enables perturbations and specifies the relative file location and name of the perturbations *.csv. ../param_pert.csv, for example. |
b |
perturbations |
RATETYPE |
Enables perturbations and specifices a global perturbation factor to be applied to ALL reactions. The `-b |
e |
data-export |
string |
Relative file location and name of the output of the logged data *.csv. ../data_out.csv, for example. |
o |
specie-option |
string |
Specify which species the logged data csv should output. This argument is only useful if `-e |
i |
data-import |
string |
Relative file location and name of *.csv data to import into the analyses. ../data_in.csv, for example. Using this flag runs only analysis. |
a |
analysis |
string |
Relative file location and name of the analysis config *.xml file. ../analyses.xml, for example. USING THIS ARGUMENT IMPLICITLY TOGGLES ANALYSIS. |
c |
cell-total |
int |
Total number of cells to simulate. |
w |
total-width |
int |
Width of tissue to simulate. Height is inferred by dividing cell-total by total-width. |
r |
rand-seed |
int |
Set the stochastic simulation random number generator seed. |
s |
step-size |
RATETYPE |
Time increment by which the deterministic simulation progresses. USING THIS ARGUMENT IMPLICITLY SWITCHES THE SIMULATION FROM STOCHASTIC TO DETERMINISTIC. |
t |
time-total |
RATETYPE |
Amount of simulated minutes the simulation should execute. |
u |
anlys-intvl |
RATETYPE |
Analysis AND file writing interval. How frequently (in units of simulated minutes) data is fetched from simulation for analysis and/or file writing. |
v |
time-col |
none | Toggles whether file output includes a time column. Convenient for making graphs in Excel-like programs but slows down file writing. Time could be inferred without this column through the row number and the analysis interval. |
how is data_out.csv (could be diff name, -e) formatted?
TODO LATER... will change next week based on what we do with analysis log
Basic Analysis calculates the average concentration level of each specie over a given time interval for each cell of a given set and across all of the selected cells. The object also calculates minimum and maximum concentration levels for each specie across the set of cells and for each cell.
MAKE EDITS DEPENDING ON IMPLEMENTATION OF ANALYSIS LOG
Oscillation Analysis identifies the local extrema of a given local range over a given time interval for a given set of cells. The object also calculates the average period and amplitude of these oscillations for each cell in the given set. MAKE EDITS DEPENDING ON IMPLEMENTATION OF ANALYSIS LOG
Concentration Check allows the user to abort simulation prematurely if a concentration level of a given specie (or for all species) escapes the bounds of a given lower and upper value for any given set of cells and time interval.
Copyright (C) 2016-2019 John Stratton ([email protected]), Ahmet Ay ([email protected])
This software was developed with significant contributions from many students over the years, including: Yecheng Yang, Nikhil Lonberg, Kirk Lange, Maxwell Brown, Yiwen Xiang, Jeremy Davis, Jack Taylor, Nick McClellan.
Copying and distribution of this file, with or without modification, are permitted in any medium without royalty provided the copyright notice and this notice are preserved. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.