Facile User Manual 0V3x

From Facile
Jump to: navigation, search

Facile User's Guide

Overview

Modelling biological systems frequently requires the translation of a set of biochemical equations into a machine-readable and simulatable form. Typically, one writes down a set of biochemical equations for the system being modelled, such as

With the appropriate rate constants and initial conditions, these equations specify the dynamics of the system. Depending on the system being modelled, one may chose to simulate the chemical master equation of the system using a Gillespie algorithm, or map the system into a set of ordinary differential equations (ODEs) and simulate with a tool such as Matlab or XPP, or analyze the system with tools such as Mathematica, Maple or AUTO. In any case, and 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 Mathematica, Matlab, XPP, EasyStoch, AUTO 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 and display results in a graphical form (?), 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

Finally, driven by the need to provide a fully reduced system to AUTO for stability analysis, Facile has the capability of identifying conserved moities 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 user in terms of integrating Facile into his/her 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_V30.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_V30"

The accompanying EasyStoch stochastic simulator is installed in a similar fashion, with an additional compilation step:

mkdir easystoch
cd easystoch
tar -xvzf easystoch_RELEASE_V6.tar.gz
cd easystoch_RELEASE_V6
make easystoch

Similarly,

export PATH="$PATH:~/toolbox/easystoch/RELEASE_V6"

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.


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 rates 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 Declarations

A variable may be declared in the EQN section for use in other expressions.

For example

variable 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. In this release only Matlab is supported.


Target Effect
Matlab The probes are plotted, if the -P switch is set.
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 of them are used to specify parameters in the input file instead of on the command-line. The variables are given in this example:

# File: example_configuration.eqn

CONFIG:
compartment_volume = 1e-15                       # in L
event_times = 1.0 2.0 ~                          # N event times indicating where to restart integration
SS_timescale  = 10 10 20 20                      # timescale of N+1 integration intervals for steady-state check
SS_RelTol = 1e-3                                 # tolerance for steady-state check in each integration interval
t_final = 100.0                                  # simulation time
t_vector = [t0:0.1:tf]                           # sampling vector for matlab
matlab_ode_solver = ode23s                       # matlab solver
matlab_odeset_options = odeset('MaxStep',0.01)   # matlab options for odeset

The SS_timescale and SS_RelTol variables can be a single value or a list of values corresponding to each integration interval. In the latter case, dummy values must be specified for intervals where no steady-state check is performed. If the number of values supplied is less than the number of intervals, then the last value supplied is re-used.

Also, XPP configuration commands are allowed:

CONFIG:
@ xlow = 0

Parameters specified on the command-line have precedence if they are also specified in the CONFIG section.

Simulating 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 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 or EasyStoch)
  • a negative value, indicating a time relative to the last event that occurred (Matlab only)
  • a tilde '~', to indicate that the system should be integrated to steady-state (Matlab only)

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 and SS_RelTol, and also on the Matlab's AbsTol integration parameter (c.f. odeset documentation). 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 AbsTol to be considered nil

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

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 and SS_RelTol 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. For example, 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.

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:

event_times = ~ 50.0 ~ -90 ~ 300 ~
SS_timescales = 10
SS_RelTol = 1e-3
t_final = 1000.0
t_vector = [t0:0.01:tf]        # sampling vector for matlab

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. Each option is described below.

#- General options
#-
#- (-H) --extended_help     Extended help screen
#- (-h) --quick_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)
#- (-p) --split             Split output into top-lvl script and sub-scripts with
#-                          containing rates and ICs (Matlab only).
#- (-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
#-
#- (-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)
#- (-l) --solver     (solver)    ODE solver (Matlab; defaults to ode23s)
#- (-n) --solver_options (option_string) Option string for odeset (Matlab)
#- (-e) --events   (event_list)  Comma-separated list of times where system parameters
#-                               change discontinuously (Matlab, EasyStoch)
#- (-C) --volume   (volume)      Specify compartment volume (in L)
#- (-P) --plot                   Plot probes (matlab only).

Usage Examples

facile.pl <input file>               # defaults to Mathematica output
facile.pl -E <input file>            # echo Mathematica eqns to screen
facile.pl -m <input file>            # generates Matlab scripts
facile.pl -s <input file>            # output for stochastic simulator
echo 'A+B->C; f=1' | facile.pl -E -  # grab equation from stdin and display

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

These are the files generated by Facile when Matlab is specified as the target simulator (-m switch). The systemname prefix must be specified specified using the -o command-line option.

File Description
run_systemname A bash shell script that starts matlab and runs the driver script (not implemented).
systemnameDriver.m Main script to call from within matlab. Initializes and runs the ODE simulation.
systemname.m Matlab function that implements the ODE system.
systemname_s.m Matlab function that converts species string to its array index.
systemname_r.m Matlab 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_simulation_input 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