============================================================
             - CHEMICAL METABOLISM SIMULATOR -
============================================================

An open-source chemical reaction modeling platform.

This project is available at:
http://github.com/jpg18/metabolism/


Support
============================================================
This project was partially supported by NSF UBM DUE-0634612.


Dependencies
============================================================

This project uses the following:

   SIMD-oriented Fast Mersenne Twister (SFMT)
      A random number generator. The SFMT libraries have
      been modified for saving and restoring state and are
      included with the source code, but the project page
      can be found at
      http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/.

   Qt
      An API for building graphical user interfaces. Used by
      default but not required (see 'Compilation
      Instructions'). Available at
      http://qt.nokia.com/.

   Qwt
      A Qt library for creating plots. Required for building
      with Qt. Available at
      http://qwt.sourceforge.net/.

   ncurses
      An API for building text user interfaces. Used by
      default but not required (see 'Compilation
      Instructions'). Available at
      http://www.gnu.org/software/ncurses/.

   Boost Iostreams Library
      Provides a framework for defining streams, stream
      buffers and i/o filters. Required for building with
      Qt or ncurses. Available at
      http://www.boost.org/.

   R
      Statistical computing software. Scripts provided for
      analysis. Not required for compiling or running the
      simulation. Available at
      http://www.r-project.org/.

   deSolve
      General solvers for R for initial value problems of
      differential equations. Required by some scripts for
      analysis. The package can be installed in R using:
         install.packages("deSolve")
      The package is also available at
      http://cran.r-project.org/web/packages/deSolve/.

   tikzDevice
      LaTeX graphics for R. Optional parameter 'outputtype'
      or 'uselatex' for each R script allows LaTeX documents
      to be built instead of PDFs or PNGs. The package can
      be installed in R using:
         install.packages("tikzDevice")
      The package is also available at
      http://cran.r-project.org/web/packages/tikzDevice/.


Compilation Instructions
============================================================

Compiles with g++ (Ubuntu 4.4.3-4ubuntu5) 4.4.3,
              gcc (Ubuntu 4.4.3-4ubuntu5) 4.4.3,
              Qt version 4.6.2.

Execute the following for a standard compilation:

   cd metabolism/src
   make

The program created is called "metabolism" (or
"metabolism.app" on Mac).

By default, Qt and ncurses are linked to the program. If
developer's libraries are not available, the program can be
compiled without ncurses, Qt, or both. Alternative make
targets are:

   make metabolism-qt
      for compiling with Qt and without ncurses
   make metabolism-ncurses
      for compiling with ncurses and without Qt
   make metabolism-minimal
      for compiling without either Qt or ncurses


An Example Simulation and Analysis
============================================================

Once compiled, try running the following from the
metabolism/src directory:

   ./metabolism --gui-off --load ../load/michaelismenten.load --iters 2000 

or on Mac:

   metabolism.app/Contents/MacOS/metabolism --gui-off --load ../load/michaelismenten.load --iters 2000

This will run a quick simulation of the Michaelis-Menten
enzyme kinetics system. Data is output to the files
diffusion.out and census.out, and the simulation parameters
are recorded in config.out. To plot the trajectory of the
reaction, run:

   ../scripts/kinetics.R

A PDF file is generated containing plots of the changes in
chemical species concentration against time. View the PDF:

   evince plots.pdf

or on Mac:

   open plots.pdf

Both observed (simulated) and expected trajectories are
plotted for all chemical species. Expected trajectories are
found through numerical integration of the rate laws. Notice
that the simulated reactions occurred at a slightly slower
rate than expected (in particular, note the rate of
"Product" formation). This is because the chemicals need to
be "mixed" or "stirred" to maximize the rate of collisions
between reactants. Now try:

   ./metabolism --gui-off --load ../load/michaelismenten.load --iters 2000 --shuffle
   ../scripts/kinetics.R pdf true config.out census.out plots-shuffle
   evince plots-shuffle.pdf

or on Mac:

   metabolism.app/Contents/MacOS/metabolism --gui-off --load ../load/michaelismenten.load --iters 2000 --shuffle
   ../scripts/kinetics.R pdf true config.out census.out plots-shuffle
   open plots-shuffle.pdf

(The extra parameters passed to the script are explained
in the comments header of the file.) This will run the same
simulation while shuffling the positions of all particles
each iteration, thereby simulating mixing. Now the simulated
reaction rates match the the expected trajectories much
better.


Using the Simulation Platform
============================================================

All command line options are available via --help and can be
used to specify simulation parameters such as the size of
the lattice or the duration of the simulation. By default,
the program produces a configuration file named config.out
which records the simulation parameters. config files can be
loaded to run exact reproductions of simulations using
--load.

To change the types of particles in the lattice or the
reactions that are defined in the reaction table, you will
need to create a custom config file or use one of the
provided files in the directory metabolism/load. Create a
new file and fill it with any parameters you would like to
specify using the format:

   keyword value

with separate lines for each parameter. Parameters that can
be specified and examples of valid values are:

   iters 5000
   seed 1234567890
   x 2000
   y 1000
   shuffle on
   reactions off

In addition to these (which are described in more detail in
--help), particle types, reactions, and conditions for
terminating the simulation can be specified using the
keywords "ele", "rxn", and "extinct". Add a line beginning
with "ele" for each particle type ("element") using the
following format:

   ele name symbol color conc

"name" is a string that names the particle type. "symbol" is
a single character that is used to represent particles of
that type in the ncurses text interface. "color" is a string
that names the color that will be used by the Qt gui to
distinguish particles of that type. The following colors
are used by the default settings and are distinguishable by
users with and without color deficient vision: "teal",
"hotpink", "darkorange", and "yellow". All valid color names
can be found at
http://www.w3.org/TR/SVG/types.html#ColorKeywords.
Finally, "conc" is a real valued number between 0 and 1 that
represents the starting concentration for the particle type.

An example will help clarify how to build a config file. The
Michaelis-Menten enzyme kinetics system consists of three
reactions: E + S -> ES, ES -> E + S, and ES -> E + P. The
particle types for this system could be defined using:

   ele Enzyme E yellow 0.1
   ele Substrate S teal 0.3
   ele ES C darkorange 0.0
   ele Product P hotpink 0.0

(The starting concentrations were chosen arbitarily, and the
concentration for Substrate is 3 times the concentration for
Enzyme; no ES or Product will be initially present in the
lattice.)

When specifying reactions, each unique reaction should
receive its own line beginning with "rxn". Up to two
reactants can be specified (reaction orders higher than
2 are not supported), and up to two products can be
specified. Use the format:

   rxn prob [n] reactant [+] [n] [reactant] -> [n] product [+] [n] [product]

"prob" is the probability between 0 and 1 of the reaction
occurring when the reactants "collide". "reactant" and
"product" are the names of particle types. If you would like
to define a zeroth-order reaction, or if you would like your
reactants to annihilate one another, use the placeholder
symbol "*" for the reactants or products, respectively
(e.g., * -> A, or A -> *). "n" is a stoichiometric
coefficient. Since at most two reactants and two products
can be specified, this number must be either "1" (in which
case it may be omitted) or "2" (in which case there can be
no other reaction participants listed on that side of the
reaction arrow). Additionally, at most two separate
reactions can be specified with the same set of reactants.
All particle types must be listed in a config file before
any reactions that use them.

For our example, the Michaelis-Menten reactions could be
listed using:

   rxn 0.50 Enzyme + Substrate -> ES
   rxn 0.01 ES -> Enzyme + Substrate
   rxn 0.01 ES -> Enzyme + Product

with "prob"'s = 0.5, 0.01, and 0.01.

Finally, the user can optionally define special conditions
under which the simulation should terminate early. A set of
particle types can be specified using:

   extinct n name [...]

"n" is the number of types in the set, and "names"'s are the
names of the particle types. If at any time during the
simulation there exist no particles with types from this
list, the simulation will end. More than one "extinct" line
can be listed in a config file to specify alternative
termination conditions. Types included in these lists must
be defined earlier in the config file.

For our example, the system will no longer undergo any
reactions when all Substrate and ES is depleted, and so we
would not be interested in simulating beyond that point.
Thus, we could use:

   extinct 2 Substrate ES

to force the simulation to end when there are zero Substrate
and zero ES particles in the lattice.

If a config file is not used, or if a config file is used
that lacks any lines that begin with "ele", then the types
"A", "B", "C", and "D" are defined by default. If the config
file lacks any "rxn" lines and the default particle types
are used, the reaction A + B -> C + D (prob = 0.5) is
defined. If the config file does not have any "extinct"
lines and the default particle types are used, then the
simulation will terminate when either all the A particles or
all the B particles are depleted.

============================================================