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.

| 0: System Requirements
|__ 0.0: Operating System
|__ 0.1: Compilers
| 1: Tutorial with Simple Model
| 2: Compilation and Model Building
|__ 2.0: Species and Reactions
|____ 2.0.0: Declaring Species
|____ 2.0.1: Declaring Reactions
|____ 2.0.2: Defining Reaction Rate Formulas
|____ 2.0.3: Defining Reaction Inputs and Outputs
|__ 2.1: Compiling and Generating Parameter Templates
|__ 2.2: Parameters
|____ 2.2.0: CSV Parser Specifications
|____ 2.2.1: Parameter Sets
|____ 2.2.2: Perturbations
|____ 2.2.3: Gradients
|____ 2.2.4: Analyses XML
| 3: Running the Simulation
|__ 3.0: Description of the Simulation
|____ 3.0.0: Preamble
|____ 3.0.1: Deterministic
|____ 3.0.2: Stochastic
|__ 3.1: Input
|____ 3.1.0: Required Files
|____ 3.1.1: Optional Files
|____ 3.1.2: Command Line Arguments
|__ 3.2: Output
|____ 3.2.0: Simulation Log
|____ 3.2.1: Analysis
|______ 3.2.1.0: Output Destination
|______ 3.2.1.1: Basic Analysis
|______ 3.2.1.2: Oscillation Analysis
|______ 3.2.1.3: Concentration Check
| 4: 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.

Back to Top

step-by-step instructions on whole process using a simple model

Back to Top

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.

1. Empty cells, blank rows, and whitespace

   To illustrate, the following two examples are equivalent.

   3.14, , 2001, -2.18,
   
   41,       2.22e-22
   

   3.14,2001,-2.18,41,2.22e-22
   

2. Comments, i.e. rows that begin with a #

   # I am a comment! Below is the data.
   9182, 667
   

3. Any cell that contains a character which is not a number, ., +, -, or e.

   Only the following scientific notation is supported:

   -0.314e+1, 3.00e+8, 6.63e-34
   

   These would be considered invalid:

   3.33*10^4, 1.3E-12, 4.4x10^2
   

   Often 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

add instructions here on analyses.xml

Back to Top

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.

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.

Back to Top

Copyright (C) 2016-2017 John Stratton ([email protected]), Ahmet Ay ([email protected])

This software was developed with significant contributions from several students over the years: Yecheng Yang, Nikhil Lonberg, and Kirk Lange.

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.

Back to Top