Facile: a command-line network compiler for systems biology

Input file

Facile's input is simple: a system of chemical reactions. Facile automatically calculates reaction rates using the law of mass action. Therefore, for example, the Michaelis-Menten reaction is specified as

E + S <-> C ; f = 1.6e7; b = 6

C -> P + E ; k = 0.15

The association rate between the enzyme, E, and the substrate, S, is 1.6 × 10⁷ M^-1 s^-1, and the complex, C, dissociates into the product, P, and enzyme with rate 0.15 s^-1. Backward reactions can also be written explicitly, e.g.

E + S <- C ; b = 6

Sources, for example gene expression, and sinks, for example protein degradation, are denoted by null:

null -> A ; s = 10

A -> null ; d = 0.1

implying that protein A is produced at a rate of 10 s^-1 and degrades at a rate of 0.1 s^-1.

Facile also allows explicit specification of reaction velocities, time-varying system parameters, and effective rate expressions, such as Hill functions. A thick arrow, =>, indicates a reaction velocity and overrides the law of mass action. The reaction velocity may be a constant or an expression enclosed in double quotes. For example, the Michaelis-Menten equations could be specified as

E + S < = > C; f = "1.6e7*E*S"; b = "6*C"

C = > P + E; k = "0.15*C"

or as

S = > P; f = "S*V/(K + S)"

in the quasi-steady state approximation. Variables can also be defined, such as

variable V = 0.15*0.2

variable K = (0.15+6)/16

for the quasi-steady state approximation with a total enzyme concentration of 0.2 μ M. Once a variable has been specified, its value is available throughout the rest of the input file.

Time can be included in a reaction rate by using t. For example, the expression of a gene may be periodically modulated by the cell cycle:

null -> A ; s = "10*(1+cos(2*pi*t/2400))"

for a cell cycle of 40 minutes or 2400 seconds.

To run simulations, initial concentrations are also needed. Concentrations are specified beneath the INIT keyword in the input file. Any concentration not specified is set to zero. For the Michaelis-Menten example, we have

INIT

E = 0.2 uM

S = 3

implying that P is initially zero, S is initially 3 M, and E is initially 0.2 μ M. The substrate could also be initialized as S = 3 M.

Output for analytical analysis

Mathematica and Maple are two standard computer algebra packages. They can be used to derive analytical expressions for steady-state concentrations, for the eigenvalues determining the stability of the steady-state, or to evaluate perturbation expansions, such as the method of multiple scales.

The command

facile.pl -M model

converts the textual input for the Michaelis-Menten system shown earlier (and contained in the file model) into Mathematica format. Since Facile automatically applies the law of mass action, it generates the expressions

dEdt = + b C - f E S + k C ;

dSdt = + b C - f E S ;

dCdt = + f E S - b C - k C ;

dPdt = + k C ;

These expressions can be cut-and-pasted or loaded into Mathematica and then algebraically manipulated. For example, they can be solved for the quasi-steady state solution. The -L option generates an equivalent format for Maple.

Output for simulation

Analytical solutions are usually only possible for an approximate form of a model, if at all. Facile provides an ordinary differential equation version of the model for input to Matlab and Octave [18], via the -m option, or XPP, via the -x option. These software packages provide tools for numerically integrating ordinary differential equations and visualizing their solutions. Two files are created for Matlab: one describing the model as a system of ordinary differential equations and another (driver) file for the user. An ode file is created for XPP.

Nevertheless, intracellular molecular numbers are often small, and significant stochastic fluctuations can exist in biochemical networks [19]. Through its SBML output, Facile can interface with many stochastic simulators, but it is also supplied with EasyStoch, a stochastic simulator written in C that supports time-varying kinetic rates. EasyStoch is based on the Gibson-Bruck version [20] of the Gillespie algorithm [21].

Output for bifurcation analysis

A bifurcation analysis is frequently used to show how the qualitative behaviour of a system changes when parameter values are altered. It is particularly important in systems biology where many parameters are unknown and where the dynamical behaviour of the systems is often expected to be robust to parameter changes [22].

AUTO is a standard tool for numerical bifurcation analysis [8]. Like any bifurcation analysis tool, AUTO requires all mass conservations in the system to be explicitly written and the model to be presented in a reduced form.

Facile automatically finds any mass conservations and incorporates them into a reduced system model which is only a function of the independent variables (those chemical species that are not related by mass conservations). For example, consider again the Michaelis-Menten system. It has two mass conservations: the total amount of enzyme is conserved, in either the E or C form, and the total amount of substrate is conserved, in either the S, C, or P form. Running Facile with the -rM option, produces a Mathematica version of the model in the reduced form:

E = - C + E_tot ;

S = - P - C + S_tot ;

dCdt = + f E S - b C - k C ;

dPdt = + k C ;

where E_tot is the total amount of enzyme and S_tot is the total amount of substrate. Facile determines the numerical value of these two constants from the initial conditions specified in the input file.

Once the mass conservations have been identified, there is not a unique reduced form for the system. For example, C and P could be related to E_tot and S_tot in the Michaelis-Menten example, and E and S would then be the independent variables. In the input file to Facile, the user can specify the chemical species that should be independent and those that should be dependent. Any chemical species not specified will have its dependency assigned automatically. The MOIETY section in the input file is used for these specifications. For example,

MOIETY

dependent E, S

ensures that the Michaelis-Menten equations are reduced to the form above. So, too, does

MOIETY

independent C, P

A combination of independent and dependent variables is also valid.

Although XPPAUT does allow bifurcation analysis (and is supported by the -x option in Facile), the complex biochemical networks studied in systems biology usually require the stand-alone version of AUTO because it allows the user to tune more of AUTO's parameters. Via its -a option, Facile produces a C version of the reduced model in the format expected by AUTO. AUTO also needs to know the parameters to be varied for the bifurcation analysis. These are specified in the BIFURC_PARAM section of the Facile input file and are incorporated into the C version of the reduced model. For example, if the user wished to investigate whether the system bifurcates as the association rate, f, and dissociation rate, b, for the Michaelis-Menten reaction are varied, Facile requires

BIFURC_PARAM

f, b

in its input file.

As a starting point for generating the bifurcation diagram, bifurcation analysis software uses the concentrations of the chemical species in a system at an attractor of the system. Usually a starting point is found by integrating the system to steady-state. Facile with its -A option will automatically load the values of the steady-state concentrations from a text file, generated by Matlab or XPP, and incorporate them into the C version of the model required by AUTO.

Facile greatly speeds up bifurcation analysis of biochemical models. Here we illustrate the steps that a user could follow to set up a bifurcation analysis. Facile can be used to generate both the simulation file and the C AUTO file. The simulation file is useful for finding the initial point for the bifurcation analysis and for exploring the dynamics of the model in different regions of the bifurcation diagram. As an example, we will use XPP as the integrator. We specify the model as a system of chemical reactions in the text file model. For XPP output, we run facile.pl -rx model at the command-line to generate model.ode. This XPP file is a reduced form of the model in the form of ordinary differential equations: rates are automatically calculated by the law of mass action and any mass conservations are found and used to reduce the model to the minimum number of independent variables. We then integrate the model to a steady-state and use the XPP data browser to save the corresponding time course data to a space-delimited text file, output.dat. Alternatively, we could run XPP without an interface with the command xppaut model.ode -silent. This command runs XPP and saves the time series generated in output.dat. By running facile.pl -aA output.dat model, we create a C version of the model, which includes the steady-state concentrations taken from the last time point of output.dat. We then directly load into AUTO the model.c file created. Matlab could also be used to integrate the model equations (with the -rm option in Facile). In summary, the single-line command facile.pl -rx model; xppaut model.ode -silent; facile.pl -aA output.dat model builds both the C model expected by AUTO and an XPP simulation file.

SBML export

SBML export is enabled by the -S option and uses the libSBML package, which must also be installed.

Miscellaneous

As mentioned, variables or parameters can be defined in the Facile input file using the variable command. The input file also has a CONFIG section where if desired the user can specify the time interval to run simulations, the particular differential equation solver to be used by Matlab, and XPP configuration commands. See the online Facile manual.