Facile User Manual 0V5x
Contents
- 1 Facile User's Guide
- 1.1 Overview
- 1.2 Installation
- 1.3 Command-Line Parameters
- 1.4 Usage Examples
- 1.5 Input File Description
- 1.6 Simulating a system of ODE(s) over multiple intervals
- 1.7 Simulating Dynamic Rates with EasyStoch
- 1.8 Output File Descriptions
- 1.9 Known Bugs
- 1.10 Future Work
- 2 Programmer's Reference Manual
Facile User's Guide
Overview
Many models in systems biology are expressed as a set of biochemical equations, for example:
Together with the appropriate biochemical parameters and initial conditions, these equations specify the dynamics of the system.
Often, the modeller desires to simulate the dynamics using a deterministic approach, in which case the system must be transformed into a set of ordinary differential equations (ODEs) and then simulated with a tool such as Matlab or XPP. Here, modellers will sometimes use AUTO for bifurcation analysis as well. Another alternative is stochastic simulation using a Gillespie-type algorithm such as EasyStoch, which generates realizations of the system's chemical master equation. Finally, the modeller may prefer an analytical approach using tools such as Mathematica or Maple.
However, to simulate or analyze a model using any of these methods, the model must first be translated into a a machine-readable form. Particularly for larger systems, this mapping process can be tedious, error-prone and time-consuming, as is the creation of scripts customized to the target simulation/analysis platform.
Facile was created to automate this process of translating a system of chemical equations into a machine readable form. Starting from a simple textual description of the biochemical system, Facile can generate simulation and analysis scripts for Matlab, EasyStoch, XPP/AUTO, Mathematica and Maple.
For example, if the user wishes to simulate the system using Matlab, then given the above system and a set of initial conditions, Facile will generate the following ODEs:
The ODEs and initial conditions are encapsulated in a Matlab function and scripts are created to run the simulation, display results in a graphical form, and/or write them to a file.
While the simplest type of simulation would involve specifying a vector of initial conditions and a simulation interval, Facile also has features to help the user describe more complex situations (when this is supported by the target simulation platform). These features include:
- support for time-varying parameters
- support for rate-law or mass-action law equations
- specifying times where parameter values change discontinuously (Matlab only)
- specifying times where dynamic rate parameters must be sampled (EasyStoch only)
Finally, driven by the need to provide a fully reduced system to AUTO for stability analysis, Facile has the capability of identifying conserved moieties and reducing the system to contain only independent variables.
Facile is a command-line and file-oriented tool. As such, it has a multitude of options which confer considerable flexibility to the modeller in terms of integrating Facile into a toolset and simulation methodology. When possible, many configuration parameters specific to the target platform (e.g. specific integration engine to be used in Matlab, initial/final times) can be specified on the command-line, in the input file, or both.
Installation
Facile is straightforward to install in any typical UNIX-like setup. The following installation instructions assume that the standard tools make, tar, gcc are present in the user's path. A standard installation of the Perl scripting language (version 5.8 or higher) is also assumed, and system administrator priviledges may be required to install CPAN Perl modules.
Supported Platforms (UNIX/Windows)
The application was developped in a Linux environment, but should function correctly on any UNIX-like platform (including OS-X) which provides a suitable shell (e.g. bash) along with standard UNIX tools and commands, such as tar and the Perl interpreter. If the accompanying EasyStoch simulator is used, the platform should also include GNU Make, and the GNU C compiler.
On Windows, the easiest solution (and the only one we have tested) is to install CYGWIN. This solution provides a bash shell along with many standard UNIX commands and applications in the Windows environment. When installing, be sure to include the make, gcc, tar and perl packages. CYGWIN Home Page
Installation using tarball
The latest version of Facile is available in a compressed tar format by FTP from www.cnd.mcgill.ca/~swain/downloads/facile
The files can be decompressed with the following commands:
mkdir facile cd facile tar -xvzf facile_RELEASE_0V50.tar.gz
It is NOT necessary or recommended that you copy the scripts into your home directory. Instead, you can add the installation directory to your PATH by adding the following line to your .bashrc file:
export PATH="$PATH:~/toolbox/facile/RELEASE_0V50"
The accompanying EasyStoch stochastic simulator is installed in a similar fashion, with an additional compilation step:
mkdir easystoch cd easystoch tar -xvzf easystoch_RELEASE_0V90.tar.gz cd easystoch_RELEASE_0V90 make
Similarly,
export PATH="$PATH:~/toolbox/easystoch/RELEASE_0V90"
Installation of CPAN modules
If your installation of Perl does not already include them, you may have to install a small number of CPAN modules. For example, the module Class::Std is required, and if it is missing your systems administrator can install it (and any other missing module) as follows:
cpan -i Class::Std
Installing the SBML library (libSBML)
If you wish to use the SBML export feature of Facile, you need to build and install libSBML using the appropriate switches for Perl support, then update Perl's library search path. For example, after downloading libSBML, decompressing and changing to the directory containing the configure script, execute the following commands (we assume administrative privileges):
./configure --with-expat --with-perl make make install
If Perl cannot find the LibSBML module, you will have to update Perl's library search path. For example, add a line similar to one of the following to your shell's initialization script (typically .bashrc), adjusting as appropriate
# e.g. if library was installed in home directory: export PERL5LIB=$PERL5LIB:/home/username/mylibs/lib/perl5/site_perl/5.8.6 # e.g. if library was installed in system directory: export PERL5LIB=$PERL5LIB:/usr/local/lib/perl5/site_perl/5.10/i686-cygwin
You may also have to set the following environment variable for the library to load at runtime:
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib
Note that the SBML export function was tested using v4.3 of libSBML. Earlier versions of libSBML may be incompatible.
More information is available at
http://sbml.org/Software/libSBML/docs/cpp-api/libsbml-installation.html and
http://sbml.org/Software/libSBML/docs/cpp-api/libsbml-accessing.html
N.B. There are no CPAN modules to install for libSBML.
Verifying the Installation
Check that the basic installation is correct by asking for the help page:
facile.pl -h
This will also reveal whether any CPAN modules are missing. If so, install them using the cpan -i command.
Installing downstream tools
Facile merely generates the scripts required for simulation and analysis of your system. It is up to the user to properly install Matlab, XPP, AUTO etc., according to the vendor's instructions. Refer to your systems administrator for help.
Command-Line Parameters
Facile has a sophisticated command-line interface which gives considerable to the user in how the program is used. Some of these options can also be configured in the input file, but command-line arguments override the values given through the CONFIG section. Each option is described below.
#- General options #- #- (-H) --extended_help Extended help screen #- (-h) --help Quick help screen #- #- (-V) --verbose More verbose output #- (-q) --quiet Quiet mode, suppresses some warnings #- (-E) --echo Echo output to STDOUT #- (-o) --prefix (prefix) Specifies prefix for output files #- #- Output type #- #- (-M) --mathematica Mathematica input file (default) #- (-L) --maple Maple input file #- (-m) --matlab Matlab driver script and functions (optionally with -r) #- (-O) --octave Octave driver script and functions (optionally with -r) #- (-y) --python Python/SciPy driver script and functions (optionally with -r) #- (-p) --split Split output into top-lvl script and sub-scripts with #- containing rates and ICs (Matlab, Octave, SciPy only). #- --extern In driver script, don't set external parameters and/or ICs #- to a value in equation file (Matlab/Octave only). This allows #- those parameters and species designated in the CONFIG section #- to be set outside the driver. #- --fast_dydt Use a refactoring algorithm to make dydt function more efficient. #- Caution: compare results w/ and w/o to ensure correctness of #- refactoring result. #- (-x) --xpp XPP input file (optionally with -r) #- (-s) --easystoch EasyStoch simulation input file and Matlab conversion #- script for EasyStoch output. #- (-a) --auto AUTO output, implies (-r) (optionally with -A) #- (-S) --sbml Export to SBML file. #- #- (-r) --reduce Reduce equations using moiety and constraint identification algorithm #- (-R) --matrix Same as -r except that stoichiometry matrix is output to a file. #- #- AUTO related #- #- (-A) --ss_file (file) Reads initial solution from file #- #- Simulation parameters #- #- (-t) --t_final (time) Final integration time 'tf' (Matlab; default NO_TIME_SPECIFIED) #- (-v) --t_sampling (vector) Sampling time list for state vector (Matlab; defaults '[t0 tf]') #- (-k) --t_tick (tick) Progress messages tick interval (Matlab/SciPy) #- (-l) --solver (solver) ODE solver #- For Matlab, overrides matlab_ode_solver configuration variable, #- which defaults to "ode23s". #- For Octave, overrides octave_ode_solver configuration variable, #- which defaults to "lsode". #- (-n) --solver_options (options) ODE solver options #- For Matlab, overrides matlab_ode_options configuration variable. #- For Octave, overrides octave_ode_options configuration variable. #- (-e) --ode_event_times (event_list) Comma-separated list of times where system parameters #- change discontinuously (Matlab) #- (-C) --volume (volume) Specify compartment volume (in L) #- (-P) --plot Plot probes (matlab only). #- #- Miscellaneous #- #- -- run Run the selected target simulator. #-
Usage Examples
facile.pl --run -P -m <input file> # generates Matlab scripts with commented-in plot commands, runs Matlab facile.pl -s <input file> # generates input file for EasyStoch stochastic simulator facile.pl <input file> # defaults to Mathematica output facile.pl -E <input file> # echo Mathematica eqns to screen echo 'A+B->C; f=1' | facile.pl -E - # grab equation from stdin and display
Input File Description
Facile requires an input file specifying the system to be simulated and/or analyzed. The system is specified by giving mass-action and rate law equations, rate constants and initial conditions (ICs). It is divided into two sections, the first giving the rate equations and the second the ICs. Comments are permitted anywhere, either inline or on separate lines. Comments must start with '//' or '#'.
Unless otherwise specified, the first section is an EQN section.
EQN Section
Reactions
Rate parameters are given using units of molarity (mol/L, or M) and seconds, as appropriate. Hence a forward rate constant for a simple binary reaction will have units of M^{-1} s^{-1}, and the rate constant for a unimolecular reaction is given in s^{-1}.
An EQN section is indicated by 'EQN' on a separate line, or is assumed to be the first section. It contains reaction equations whose kinetics follow a mass-action law or an explicit rate law, as indicated by the notation used. All the forms below are permitted, specifying reversible, forward, or backward reactions. Rate constants must be given along with each reaction, as appropriate (i.e. two constants for reversible reactions). For a reversible reaction, the forward rate is specified first. Use a semicolon to separate the equation and multiple rates. The maximum number of terms on either side of each reaction is unlimited, however this has not been thoroughly tested beyond two terms.
# Mass-action law examples A1 + B1 <-> C1 + D1; f1 = 10.0; b1 = 1.5 # reversible A2 + B2 -> C2; f2 = 1.5 # forward G + H <- I b2 = 2.0 # backward M N -> M f3=50 # reactions can also be written without the +'s.
For example, the ODEs for A1 and I will be
and
respectively.
Using the following notation, the user indicates that the reaction does not follow the mass-action law, but rather that the reaction velocity is specified explicitely by the user. To indicate that the rate is not a constant, it may be optionally enclosed in quotes (see Expressions below). Note that careless use may result in non-physical behaviour (e.g. if the reaction velocity is constant, then some node concentrations may become negative). Also note that stochastic simulator (easystoch) does not support this.
# Rate-law examples W + X <=> Y fl="W*X"; bl="2.5*Y" # reversible, quoted expressions X + Y => Z f2=X*Y # forward only, unquoted expression P + Q <= M + N f3="2*M*N" # backward, with 2 terms on each side
The resulting output will, for example, include an ODE for W, , and Z, .
Use the keyword null for a source or a sink. Sources are useful as inputs to the system, and sinks are used to wash away an output product. They can also be combined to clamp a species' concentration to a specific value.
# example of sources and sinks null -> X; src1=10.0 # (source) species X is created at a rate of 10Hz Y -> null; sink1=2.0 # (sink) species Y is destroyed at a rate of 2*[Y] Hz = 2 mol/L/s Z => null; sink2=1.5 # (sink) species Z is destroyed at a rate of 1.5 mol/L/s (n.b. NOT a mass-action law) null -> X; src2=0.5 # (source) duplicate: X is created at rate of src1+src2=10.5 M/s null <-> W; src1; sink1 # (clamp) a source and a sink are combined to clamp the concentration at 5 mol/L
Each numerical rate is assigned to a variable, which can be re-used in several places. A variable can also be assigned to another variable, as shown below.
A + B <-> C; f1=1.3; b1=2.2 # forward and backward rates are specified for bi-directional reaction D + E -> F; f1 # the rate 'f1' is re-used X + Y <- Z; f2=f1 # the rate f2 takes the value of 1.3
Expressions
To support more sophisticated simulations, there is limited support for expressions when specifying rate parameters and/or variables. Rate expressions allow more complex values to be used instead of simple numerical constants. The possible uses of this feature include:
- specifying time-varying sinks and sources, which can be used to construct sophisticated stimulus patterns
- specifying the reaction velocity when it is written in a rate-law form (the velocity could be a function of time or other state variables)
A rate expression may optionally be enclosed in quotes. This was compulsory in releases prior to RELEASE_v23.
The precise effect of using a rate expression depends on the type of simulation performed. The expression is copied verbatim in the associated differential equations for all targets except EasyStoch (see below). Hence any valid expression can be supplied as a rate. The following example is valid for all targets except the stochastic simulator
A + B <-> C; f1=2; b1=1 C -> D; f2=3*f1*A/(1+A)
For a stochastic sim using EasyStoch, the expression is evaluated at t=0 and at each time supplied as a simulation event. Facile evaluates the expression at each event time and includes the resulting rate constant values and their associated event times in the simulation input file of the stochastic simulator. The stochastic simulator will use each evaluated value starting from the specified event time until the another event changes the value of the parameter. Hence only discontinuous, discrete rate value changes are supported. The expression should evaluate to a constant between the supplied event times if results are to be compared to deterministic ODE simulation of the system.
The supplied expression can be a function of time (t) but not of the concentration of any nodes. However, the expression may include other parameters or variables. If it does, their values are substituted in until a numerical value is obtained for the rate. In addition, the function square() can be included in the expression (the definition is the same as that in Matlab), and also the value 'pi'.
This example would work both with Matlab and EasyStoch as a target, since node concentrations do not appear in the rate expression.
A + B <-> C; f1=0.7; b1=0.2 # binary reversible reaction C C -> D; f2=2 # the '+' is optional X -> null; f3=square(2*pi*t)+1 # time-varying sink, T=1 s, A=2Hz Y -> null; f4=2*f3 # expression OK because evaluates to constant
This example works only in Matlab, since node concentrations do appear in the rate expression.
node RNAp # declare RNAp, since does not appear in reactions variable mRNA_synth_rate=100.0 # synthesis rate variable mRNA_degradation=0.01 # degradation rate variable one_day=24*60*60 # one day in seconds variable sr=mRNA_synth_rate*RNAp*(square(2*pi*t/one_day)+1) # mRNA creation during 1/2 the day null -> mRNA; sr mRNA -> null; dr=1.0 # degrades all the time
Variable and Parameter Declarations
A variable or parameter may be declared in the EQN section for use in other expressions. The keywords variable and parameter are synonymous.
For example
parameter c=3 A -> B; k=2*c
Another example is
variable one_day_in_s = 86400 variable freq_CA = 10 # 10Hz variable conc_CA = 5.4e-12 * (1 + sin(2*pi*freq_CA*t/one_day_in_s))
Explicit Node Declarations
Using a node in the context of a reaction specification implicitly declares it. It is also possible to explicitly declare a node without having it appear as a substrate or product. This is useful in situations where the node concentration appears in a rate expression but not in a reaction.
node A, B, C # declare a comma-separated list of nodes
Duplicate reactions
Duplicate reactions are legal and cause extra terms to appear in the rate equations. This is useful when a species has identical and indistinguishable binding sites. For instance, if a protein X has two identical, non-distinguishable binding sites for a ligand L, and the the reaction kinetic parameters with respect to a single binding site are known, one may write:
# File: example_identical_binding_sites.eqn # # example of two identical ligand binding sites # (written as repeated reaction) X + L <-> XL; f_xl = 10.0; b_xl = 1.0 # L binds to first binding site X + L -> XL; f_xl; b_xl # L binds to second binding site, non-reversible! XL + L <-> XLL; f_xl; b_xl # X binds second ligand
Note that only one of the duplicate reactions is written as a reversible reaction, to get the correct backward rate. This yields the following ODEs, where the terms due to the second forward reaction are highlighted:
Another way to write the above system is to realize that the total forward rate for X binding L is double the rate for each individual binding site. This is fine, but is perhaps less convenient when generalizing to multiple binding sites, and it obscures the fact that there are in fact two ligand binding sites. Also, we need to define an extra parameter:
# File: example_identical_binding_sites_alt.eqn # # example of two identical ligand binding sites # (written with adjusted rates) X + L <-> XL; f1_xl = 20.0; b_xl = 1.0 # L binds to first binding site XL + L <-> XLL; f2_xl = 10.0; b_xl # X binds second ligand
PROBE Section
The PROBE section allows the user to specify that certain variables and nodes are to be probed. Any node and variable declared (implicitely or explicitely) in the EQN section can be probed simply by listing them. In addition, new probes can be defined in the same manner as any expression.
An example PROBE section is shown below, showing that nodes E, G and S1 are to be probed, along with variables sourceWf, v1, v3 and v4 (which is defined in the probe section).
PROBE: probe E, v1, v3 probe G probe S1 probe v4="v3 + v1" probe sourceWf
The precise effect of the PROBE section is dependent on the target simulator.
Target | Effect |
Matlab | If the -P switch is used, commands for plotting the probes are included in the driver file. If -P switch is not enabled, the same commands are included but commented out in the driver file. |
Octave | If the -P switch is used, commands for plotting the probes are included in the driver file. If -P switch is not enabled, the same commands are included but commented out in the driver file. |
SciPy | If the -P switch is used, commands for plotting the probes are included in the driver file. If -P switch is not enabled, the same commands are included but commented out in the driver file. |
all others | No effect. |
INIT Section
This section begins with the INIT keyword and contains the initial condition specification. Any nodes not specified are assumed to start at zero.
N.B. Versions of Facile RELEASE_V12 and earlier used the DATA keyword instead of INIT.
INIT A = 1.1uM; B = 2.3uM
Initial conditions are given in micro-molarity (micro-mol/L, suffix uM), molarity (mol/L, suffix M or no suffix), or number of molecules (suffix N). When generating an input file for stochastic simulation, a compartment volume must be specified if ICs are given in M or uM. When generating input files for ODE simulation, a compartment volume must be specified if ICs are given in molarity units.
INIT: A = 1.1M # 1.1 mol/L B = 2.0 # 2.0 mol/L C = 2.3uM # 2.3 micro-mol/L D = 99N # 99 molecules
In the conversion from molarity to # of molecules (or vice-versa), the formula used is
,
where AVOGADRO = 6.022e23.
The units for rate constants specified in the EQN section are assumed to be composed on molarity and seconds (e.g. M^-1s^-1 for binary reactions, s^-1 for unary reactions, etc.). However, this is only relevant if Facile must convert between molarity and number of molecules (e.g. when targeting stochastic simulations, such a conversion is necessary). Otherwise, as long as the user is careful to use the same units in both the EQN and the INIT section of the input file, the results obtained should be correct.
PROMOTER Section
Facile supports the designation and specification of promoter and DNA-protein complexes, which require special treatment during (stochastic) simulations. This feature is intended to be used with a stochastic simulator that supports cell-cycle division. Here is an example input file:
# File: repression.eqn # repression of reporter gene # D: promoter # C: promoter/RNAP complex # T: transcribing RNAP # mR: mRNA # mC1: mRNA/degradosome complex # mC2: mRNA/ribosome complex # mT: translating ribosome # A: protein # R: repressor # Z: repressor/DNA complex # EQN: D <-> C ;f0= 0.42; b0= 0.1 C -> T D ;k0= 0.01 T -> mR ;v0= 0.3 mR -> mC1 ;mf0= 0.114 mC1 -> null ;d0= 0.1 mR <-> mC2 ;mf1= 4.0; mb1= 0.4 mC2 -> mT mR ;k1= 0.3 mT -> A ;v1= 0.048 A -> null ;d1= 2.0e-4 # repression R D <-> Z ;f2= 1.0e7; b2= 0.03 # PROMOTERS: D; promoter_replication_time = 0.4 C; associated_promoter = D Z; associated_promoter = D INIT: D = 1N R = 100N
Here, D is specified to be a promoter with a copy number of 1. The promoter is replicated (copy number changes to 2) when t= 0.4 T, where T is the cell cycle time (T is to be specified when calling the simulation code). C and Z are protein-DNA complexes associated with the promoter D. The initial number of repressor molecules is 100 (specified as # molecules since presumably a stochastic simulation is desired). The rate constant f2 is the binary reaction rate between repressor and DNA (binary rates are always in M^-1 s^-1).
At DNA replication time (here, 0.4 T), the number of free D molecules is doubled and a new free D is created for each associated protein-DNA complex. At cell division (time T), promoters and complexes are distributed appropriately between daughter cells (the details of which depend on the target simulator).
If the stochcelldiv simulator is used (available in the SwainLab toolbox), then at cell division, each protein and mRNA has a probability of 0.5 of remaining in the cell that the simulation is following; the total amount of DNA is halved. In this case, there are 2 DNA molecules just before cell division which exist in either the D or the C or the Z states. Cell volume grows linearly and halves exactly at cell division. (See also: the file ~swainlab/toolbox/GillespieCellCycle/README.pdf for the original documentation of stochcelldiv).
For other target simulators, specification of promoter and protein-DNA complexes has no effect on the output files generated.
Finally, note that there are no consistency checks in Facile to make sure that protein-DNA interactions are appropriately specified in the EQN section of the input file given the designations in the PROMOTER section. In a future version, it may be possible to use a moiety-finding algorithm to verify, e.g., that the DNA is conserved.
MOIETY Section
Facile supports moiety identification and implementation. Moiety identification involves the identification of molecular subgroups that are conserved. For example, the enzymatic reaction E+S<->C->E+P implies two conserved quantities that could be written in terms of the total enzyme E as E_total=E+C and total S, S_total=S+C+P. Once the conserved quantities have been identified, it is possible to implement them. Implementation involves writing one of the state variable in terms of the rest of the variables that appear in the conserved quantity. The ODE of the former state variable will disappear and instead an algebraic expression is defined. Eventually this means that the number of ODEs decreases by the number of identified conserves moieties.
Each conserved quantity will normally involve more than one state variable and it is up to the user to specify which of them will become the dependent variable. The MOIETY section allows the user to specify the state variables that should be made dependent of others via a conserved moiety and/or which variables should remain unconstrained. The format for this section is:
EQN: E + S <-> C; kf=1.0; kb=0.01 C -> E + P kp=1.0 MOIETY: dependent E, S # designate node E and S as the dependent variables
The same result is obtained with
MOIETY: independent C, P # designate node C and P as the independent variables
A combination of dependent and independent nodes is also valid
MOIETY: dependent S independent C
(P.S. can have a comma-separated list on each line -- TODO: test this!!)
Re-ordering of nodes
Note that there is a "side effect" of specifying independent variables. The order of the independent variables as they appear in the MOIETY section will determine the order in which the ODEs will be written. This could be useful for targets that are sensitive to order. For instance, AUTO outputs a file that contains the bifurcation diagram for the first six state variables as they appear in <file.c> (The exact number will depend on several factors, six is just an example). The first six independent variables in the MOIETY section will be the first six state variables in <file.c> and so their solution will be output by AUTO.
BIFURC_PARAM Section
The equations-file <file.c> for AUTO should define the model parameters that may be used as bifurcation parameters. This section allows the user to specify bifurcation parameters for AUTO. The syntax is
BIFURC_PARAM: f1, f2 f3 b2
CONFIG Section
The CONFIG section allows the user to define variables which have special meaning to Facile. Some variables have a corresponding command-line parameter which if specified has precedence over the value given in the CONFIG section.
Variable | Affected Target | Description |
t_final | Matlab, EasyStoch | Total simulation time. The value is assigned to the variable tf, which can be used when specifying the vector of sampling times t_vector. If the variable tf is not included in the t_vector, then t_final has no effect. |
t_vector | Matlab | Vector of sampling times for Matlab and other ODE integrators, and specified in Matlab/Octave syntax. Defaults to "[0 tf]". Another example value is "[0:0.1:tf]", which samples every 0.1 time units. The last element of this vector determines the final simulation time, which is usually the same as tf. |
matlab_ode_solver | Matlab | Selects ode solver used by Matlab (e.g. ode45, ode23s, ode15s, etc.). |
matlab_solver_options{option}=value | Matlab | Matlab solver options. See odeset documentation. |
octave_ode_solver | Octave | Selects ode solver used by Octave (e.g. lsode, etc.). |
octave_solver_options{solver}{option}=value | Octave | Octave solver options. See lsode_options documentation. |
scipy_ode_solver | Octave | Selects ode solver used by SciPy (e.g. odeint or ode.vode). |
scipy_solver_options{solver}{option}=value | SciPy | SciPy solver options. See SciPy documentation. |
cpp_ode_solver | C++ | Selects ode solver used by C++ (only cvode is supported). |
cpp_solver_options{solver}{option}=value | C++ | C++ solver options. See Sundials/CVODE documentation. |
ode_event_times | Matlab | Event times where rate parameters change discontinuously (requiring integration to be restarted). Can also be a tilde, indicating that the system be integrated until a steady-state is reached. See ode_event.m for more details. |
SS_timescale, SS_RelTol | Matlab | Timescale and tolerance for the steady-state check in the integration intervals defined by ode_event_times. Can be a single number applied to all intervals, or up to N+1 numbers corresponding to each integration interval. If the number of values supplied is less than the number of intervals, then the last value supplied is re-used. If no steady-state check is required in the corresponding interval (i.e. no tilde in ode_event_times) then the corresponding value is ignored (i.e. it is a dummy value acting as a placeholder). |
external_parameters | Matlab/Octave | List of parameters that the user will set before running the driver files. Throwing the --extern switch when calling Facile will cause the line(s) which set these parameters in the driver file to be commented out. |
external_initial_values | Matlab/Octave | List of species whose initial condition the user will set before running the driver files. Throwing the --extern switch when calling Facile will cause the line(s) which set these ICs in the driver file to be commented out. |
compartment_volume | EasyStoch | Volume in L. |
easystoch_sample_times | EasyStoch | Times at which rate expressions should be re-sampled and the updated value appended to EasyStoch's simulation input file as a rate change. |
@ xlow | XPP | Undocumented? |
Here is an example CONFIG section:
# File: example_configuration.eqn CONFIG: t_final = 100.0 # simulation time, t_final is assigned to 'tf' variable # Matlab configuration matlab_ode_solver = ode23s # matlab solver for stiff systems matlab_odeset_options{MaxStep} = 0.01 # matlab options for odeset # Octave configuration octave_ode_solver = lsode # matlab solver for stiff systems octave_solver_options{lsode}{maximum step size} = 0.01 # matlab options for odeset # Matlab/Octave/SciPy parameters t_vector = [t0:0.1:tf] # sampling vector ode_event_times = 1.0 2.0 ~ # params change at 1.0 and 2.0, then integrate to steady-state SS_timescale = 0 0 20 0 # timescale of steady-state checks (only 3rd interval matters) SS_RelTol = 1e-3 # tolerance for steady-state check in all integration intervals external_parameters = f1 f2 # set f1 and f2 outside the driver file external_initial_values = X Y # set IC for X and Y species outside the driver file # EasyStoch configuration compartment_volume = 1e-15 # in L easystoch_sample_times = 1.0 2.0 # re-sampled times where new rate value appended to EasyStoch's simulation input file # XPP configuration @ xlow = 0 # ???
Simulating a system of ODE(s) over multiple intervals
A number of situations warrant splitting up the integration of an ODE system into several intervals. First, in the presence of discontinuous changes in system parameters, the integrator may choose a step size which skips over the discontinuous event, resulting in an incorrect solution. Splitting the integration into multiple intervals solves the problem because the solver will start over with an appropriately small initial step size. Second, it is often desired to simulate the system to steady-state and then re-start the simulation under changed conditions. For this reason, Facile provides the means to split the integration into pre-defined intervals that are separated by events.
To do this, the user supplies a vector of N event times via the configuration variable ode_event_times. The integration of the system will then be split into (N + 1) intervals, with integration re-started at each event. An event time may be either
- a positive value indicating an absolute time (Matlab only -- Octave/SciPy support pending)
- a negative value, indicating a time relative to the last event that occurred (Matlab only -- Octave/SciPy support pending)
- a tilde '~', to indicate that the system should be integrated to steady-state (Matlab only -- Octave/SciPy support pending)
These event times should fall between the initial and final times specified by the TV argument. Given a set of event times of the form event_times=[event1 event2 ... eventk ... eventn] and given the configuration variable t_vector = [t0:0.1:tf], the integration will be split into the following intervals:
[t0 event1] [event1 event2] ... [eventn tf]
Integration will always terminate at t=tf, or possibly earlier if the last event is to find the steady-state (i.e. eventn=~) and the stopping condition occurs before tf. In the case where eventn=~, only N intervals are integrated instead of (N+1). When integrating to steady-state in any interval, tf acts as a timeout value for the integration rather than specifying the actual final time.
Positive event times
Specifying a positive event time simply re-starts integration at that absolute time. This is useful when dealing with time-varying parameters in the system of ODEs which change discontinuously at specific times, to prevent the integrator from "stepping over" the discontinuity.
Negative event times
Specifying a negative event time means that the integration should continue until a time of (-eventk) has elapsed since the last event. Thus, eventk=100 indicates that the kth event occurs at absolute time t=100, while eventk=-50 indicates that the kth event occurs 50 time units after the previous event. The previous event can be an absolute time or a tilde (~) that indicates reaching steady-state. Thus it is possible to simulate to steady-state, then wait a certain amount time after steady-state is reached to re-start integration under changed conditions.
Note that the encoding of relative event times with negative numbers requires the user to specify a non-negative initial simulation time and a positive value for tf.
Steady-State events
An event time of ~ indicates that integration of the current interval should continue until a steady-state stopping condition is reached. This stopping condition is derived based on the configuration variables SS_timescales, SS_RelTol and SS_AbsTol. The stopping condition is a component-wise check that evaluates if
- the significant digits of the solution vector y (specified by SS_RelTol) are stable OR
- the magnitude of y has dropped sufficiently below the variable SS_AbsTol to be considered nil
Mathematically, the stopping condition is expressed by the inequality
dy/dt * SS_timescale < max(SS_RelTol * abs(y), SS_AbsTol)
If this is true for for all solution components, the steady- state condition is reached and integration is interrupted.
The simplest situation involving finding the steady-state is to specify a single steady-state event and a timeout. For example, specifying
ode_event_times = ~ t_vector=[0:0.1:1000]
will simulate the system to steady-state or to t=1000, whichever occurs first.
To reliably detect a steady-state condition, the user must set the configuration variables SS_timescales, SS_RelTol and SS_AbsTol appropriately. The SS_timescales variable indicates the expected timescale of the slowest process in the system. This is vector of up to N+1 timescales, allowing you to take into account that changing biochemical parameters may change the settling time of the system. If a partial vector is supplied, the last element is extended as necessary. If empty, all elements default to 100. The SS_RelTol variable, also a vector of up to N+1 elements, specifies the relative tolerance at which a solution component is considered to have reached steady-state over the given timescale. The SS_AbsTol specifies a threshold below which the solution component is considered so small that the relative tolerance should be relaxed.
To illustrate, if SS_RelTol=1e-3 and SS_timescales=10, this means that steady-state is reached when the rate of change of each component of the solution is 0.1% or less over an interval of 10 time units. However, if some components of the solution decay exponentially toward zero then the relative rate of change of these components does not decrease with time and a steady-state will not be detected. That's why it's important to set the SS_AbsTol variable to a threshold value below which the component is considered as being negligibly small and having reached steady-state.
The choice of Matlab integrator (e.g. ode23s vs ode15s) and integration parameters will affect how successful the above stopping condition is at detecting a steady-state of the system. Indeed the integration should be sufficiently precise and without numerical artifacts for reliable success at detecting steady-state. As such, Matlab's RelTol, AbsTol, InitialStep and MaxStep integration parameters (c.f. odeset documentation) should be set appropriately. In particular, the RelTol (determines precision) and AbsTol (a threshold below which solution components can be wildly inaccurate) parameters should be equal or smaller than SS_RelTol and SS_AbsTol respectively.
Two global variables may be accessed from user-specified rate expressions to implement event-dependent or time-varying equations. The variable event_flags is a boolean vector with the same size as event_times and indicates whether or not the corresponding event has occurred. Similarly, the variable event_times indicates at which time the corresponding event has occurred. For steady-state events, this time will be set to 0 until the steady-state is found, at which point it is updated to the correct value. For relative (negative) event times, the time is updated when the interval is integrated.
Example
Here is an example illustrating the use of all 3 types of events, and of the event_flags variable:
EQN: null -> X; k_srcx = 2*(t>=50 && ~event_flags(3)) # 1st stimulus at t=50 and turns off after steady-state null -> Y; k_srcy = 3*(event_flags(4)) # 2nd stimulus 90 after 1st has settled null -> Z; k_srcz = 4*(event_flags(6)) # 3rd stimulus at t=300 X -> null; k_sinkx = 1 Y -> null; k_sinky = 3.3 Z -> null; k_sinkz = 2.2 PROBE: probe X probe Y probe Z INIT: X = 4 Y = 5 Z = 2 CONFIG: ode_event_times = ~ 50.0 ~ -90 ~ 300 ~ SS_timescales = 10 SS_RelTol = 1e-3 SS_AbsTol = 1e-9 # ode15s is faster than ode23s and works on stiff systems matlab_ode_solver = ode15s # set AbsTol to be much lower than SS_AbsTol (this isn't costly), # and also RelTol to yield a couple more sig. digits than SS_RelTol requires matlab_solver_options{InitialStep} = 1e-15 matlab_solver_options{AbsTol} = 1e-48 matlab_solver_options{RelTol} = 1e-5 t_final = 1000.0 t_vector = [t0:0.01:tf] # sampling vector for matlab
Simulating Dynamic Rates with EasyStoch
Facile includes special support for EasyStoch's ability to simulate dynamic rate changes. Indeed, EasyStoch implements the Piecewise-Propensity Next Reaction Method, where rate constants can be varied in a piecewise-constant or piecewise-linear fashion. Through this feature, rate expressions of arbitrary functional form can be approximated. The following example demonstrates this. Here, note the usage of several Facile features:
- the keyword table, which allows the user to specify a list of time/value pairs
- the keyword vector, which allows the user to specify a vector of times at which an expression should be sampled
- the ability to sample a rate expression through the configuration variable easystoch_sample_times
- the ":1" notation associated with rate definitions, which tells Facile which reactions should use a piecewise-linear approximation of the rate expression. The index number of these reactions will be output in the stochastic equation file and can be passed to EasyStoch via the --pwlin switch.
#################################################################### # File: DynRates.eqn # Description: Showcase/test of dynamic rate functionalities # of the EasyStoch stochastic simulator. #################################################################### #---------------- EQN: #---------------- variable f0=1e3 variable fx2 = 2000 variable fsq = square(2*pi*(t-10)/40, 50)+1 # square wave changes at 10,30,50,70... X + Y <-> XY; rateXY = 1e6*0.5*fsq; b0 = 1 # the following tests expression evaluation of easystoch output A1 + B -> A2 + B; f1=1e6; A2 + B -> A3 + B; f2=f0*f0 A3 + B -> A4 + B; f3=3*f2 A4 + B -> A1 + B; f4=1000*fx2*fsq # list of time/value pairs table T1 = (0,1.0, 25,2.2, 55,0.55) D1 <-> E1; f10 = table(T1,t); b10=1 D2 <-> E2; f11 = 1+0.9*sin(2*pi*0.01*t); b11=1 # ':1' notation indicates reaction rates are piecewise-linear D3 <-> E3; f10:1; b10 D4 <-> E4; f11:1; b11=1 #---------------- INIT: #---------------- X = 1e-6; # about 602N in 1e-15 L Y = 2000N; A1 = 1000N; B = 1000N; D1 = 1000N D2 = 2000N D3 = 1000N D4 = 2000N #---------------- CONFIG: #---------------- t_final = 100 compartment_volume = 1e-15 # in L ode_event_times = 10 30 50 70 90 easystoch_sample_times{fsq} = 10 30 50 70 90 easystoch_sample_times{f10} = table(T1,time_list) easystoch_sample_times{f11} = vector(0,15,tf)
To try out this example, first save it in a file named DynRates.eqn, then type
facile.pl -s DynRates.eqn
Next, examine the DynRates.seqn file to determine the index number of the reactions whose rates are marked as piecewise linear. In this file, we see the following:
... # Section IV (dynamic rates/species change): reaction/species index, time of rate change, new value # N.b. the following reactions have piecewise-linear rates according # to the equation file (c.f. EasyStoch's --pwlin switch): # --> 10,12 # Dynamic rate sample times for f10: table(T1,time_list) # Dynamic rate sample times for fsq: 10 30 50 70 90 # Dynamic rate sample times for f11: vector(0,15,tf) 0 1.000000e+01 1.000000e+06 # rateXY = 1e6*0.5*fsq 5 1.000000e+01 4.000000e+06 # f4 = 1000*fx2*fsq 8 1.500000e+01 1.728115e+00 # f11 = 1+0.9*sin(2*pi*0.01*t) ...
We see that the indexes of reactions 10 and 12 should be specified via the --pwlin switch
easystoch --cout --sout --vol0=1e-15 --out_interval=0.1 --pwlin=10,12 --tmax=100 DynRates.seqn > simdata.DynRates.dat
(N.B. The --rate-vary switch was recently renamed to --pwlin)
Admittedly, the above mechanism is a bit clumsy. Future versions of Facile/EasyStoch will remove the need for the user to specify an index number via the --pwlin switch.
Output File Descriptions
Mathematica
These are the files generated by Facile when Mathematica is specified as the target simulator (-M switch).
Note: the Matlab-style exponentiation operator "**" will be substituted with "^" in the output file.
File | Description |
systemname.ma | Mathematica file. |
Matlab / Octave
These are the files generated by Facile when Matlab or Octave is specified as the target simulator (-m or -O switch). The systemname prefix must be specified specified using the -o command-line option.
File | Description |
systemnameDriver.m | Main script to call from within Matlab/Octave. Initializes and runs the ODE simulation. |
systemname_odes.m | Matlab/Octave function that implements dy/dt function of the ODE system. |
systemname_s.m | Matlab/Octave function that converts species string to its array index. |
systemname_r.m | Matlab/Octave function that converts rate string to its array index. |
EasyStoch
These are the files generated by Facile when EasyStoch is specified as the target simulator (-s switch). The systemname prefix is specified specified using the -o command-line option.
File | Description |
systemname.seqn | Input file for the stochastic simulator. |
systemname_convert.m | Matlab script to convert the data file produced by stoch into named vectors. |
XPP
These are the files generated by Facile when XPP is specified as the target simulator (-x switch). The systemname prefix is specified using the -o command-line option.
File | Description |
systemname.ode | XPP file with system definition. |
AUTO
In addition to having the capability of integrating a biological network it may be required to carry out a stability analysis. AUTO is a continuation package suitable for carrying out detailed stability analysis of a set of ODEs.
Facile generates files in AUTO2000 format. Note that Facile deals with the C version of AUTO, that is, AUTO2000. The newer version of AUTO, AUTO-07P, should be able to read C files despite the fact that it compiles in FORTRAN (see AUTO-07P documentation).
In order to run AUTO, it is necessary to supply two files, an equations-file <file.c> and a constants-file <c.file>. (See AUTO documentation). The AUTO support of Facile aims at constructing the equations-file <file.c>. The equations-file contains the definition of the problem---ODEs and parameters --- in a particular format suitable for AUTO. Note that the constants-file <c.file> is a problem-dependent file which should be supplied by the user.
This is the file generated by Facile when AUTO is specified as the target simulator (-a switch). The systemname prefix must be specified using the -o command-line option.
File | Description |
systemname.c | AUTO2000 equations-file with problem definition. |
MOIETY IDENTIFICATION
An important feature of Facile -a is the fact that it implies -r, that is, it reduces the system by identifying and implementing constraints. This is in fact necessary to run AUTO and it motivated writing the moiety identification algorithms implememnted in Facile.
AUTO starting solution, the -A switch
Facile -a may be run with the optional -A switch, which takes a data file <file.dat> as its argument. When AUTO is run for the first time it typically reads a starting solution from the equations-file <file.c>. It is common that the starting solution is a stable node found by integrating the ODEs numerically. When <file.c> is constructed using Facile -a, the starting solution is the initial condition as specified from the input file. However, this initial condition may not be the required starting solution.
The -A switch gives the user the possibility of constructing <file.c> with a starting solution. The starting solution is read from an ascii data file and the values are written in <file.c>. The data file should contain time-series data in columns. Each column is a state variable except for the left-most column which is time. The last row should be a time point when the steady state has been reached.
The -A switch works as follows. Facile reads the last row (steady state data point) of <file.dat>, and it substitutes each column value into the initial solution function of <file.c>. The first column of <file.dat> is time and it is ignored. It is assumed that the order of appearance of the state variables in <file.dat> is exactly the same as the order of the state variables in <file.c>. Making sure this is the case is simple with Facile since the problem definition is the same regardless of target application. As an example, consider using XPP to integrate the ODEs. You may do the following:
- Run Facile with the -rx switches. This will output a reduced system written in XPP, <file.ode>.
facile -rx <file.eqn>
- Integrate <file.ode> with XPP towards a steady state.
- From the XPP data browser, save the data by simply clicking on "Write". (This saves time and all state variable to text file, <file.dat>.)
- Run Facile with the -aA <file.dat> switches.
facile -aA <file.dat> <file.eqn>
Alternatively, if you are using say, matlab, you could follow this procedure:
- Run Facile with the -rm switches. This will output a reduced system ready for matlab.
facile -rm <file.eqn>
- Integrate system with matlab towards a steady state.
- Within the matlab prompt save the time-series data in the aforementioned format.
var=[t y]; %y contains time-series data save '<file.dat>' var -ascii
- Run Facile with the -aA <file.dat> switches.
facile -aA <file.dat> <file.eqn>
Other target integrators should require a similar procedure.
Compliance with C syntax: powers All mathematical expressions in the Facile input file must conform to C syntax when the -a switch is specified. Although this is normally not an issue there is one particular aspect of C which differs from the rest of the target applications such as XPP or Matlab: powers. Powers in C are expressed as pow(a,b) instead of a^b. Facile deals with this by incorporating a script that transforms a^b into pow(a,b). Therefore, it is recommended that expressions involving powers be entered as a^b regardless of target. This way all target applications will read the power expression correctly and the AUTO file will be automatically translated into the C pow(a,b) syntax.
Maple
These are the files generated by Facile when Maple is specified as the target simulator (-L switch).
File | Description |
systemname.maple | Mathematica file. |
SBML
This is the file generated by Facile when SBML output is specified (-S switch).
File | Description |
systemname.sbml | SBML Level 2 (version 4) file. |
Known Bugs
Release | Description |
... | ... |
Future Work
Target Release | Description |
RELEASE_V?? | ... |
... | ... |
Programmer's Reference Manual
This manual is intended for hackers and developers. It provides essential information for understanding how the application is structured, the algorithms used, and the testing methodology.
Design Notes & Implementation Guide
The Facile application is written in the Perl programming language.
[TODO: add basic info about parse, classes, the files which comprise the module, algorithms used, useful references, etc.)
Testcases
The Makefile included with the Facile program has a test target which tests the various features of Facile by generating output for number of test cases. During the tests each output file is compared against its reference version, which is stored in the version control system. This is useful to verify that code changes didn't break anything, and should be run prior to making a new release.
Revision Control (Bazaar)
Bazaar is used for revision control. Source code can be downloaded from http://sourceforge.net/projects/facile
Release Notes
The release notes are in the file RELEASE_HISTORY.TXT