FERN is an object oriented library implemented in Java (see Additional file 1). Although it consists of more than 100 classes and interfaces, most classes are just implementations of one of three major interfaces and abstract classes (see Figure 2):
-
1.
The interface Network provides the network structure of the model.
-
2.
The abstract class Simulator performs simulations on a Network. It additionally calls the registered observers during the simulation run.
-
3.
The abstract class Observer traces the simulation progress and creates the simulation output.
A simple simulation can be performed in only five lines of code, one line for each of: loading a network file, creating a simulator, creating and registering an observer, running the simulations and printing the results (see Figure 3). More complex examples for using FERN can be found in the FERN distribution. In the following the three layers of FERN are described in more detail.
Networks
The interface Network describes the network's structure, i.e. the reactions and species in the networks. For this purpose, each reaction and each species is described by an integer value. Furthermore, the network stores basic information like species names and their initial molecule numbers. For the simulation more information is necessary which is stored in three additional classes (see Figure 4):
-
The AmountManager controls the amount of each molecular species during the course of a simulation.
-
The AnnotationManager can store additional annotations for the network, its species and reactions.
-
The PropensityCalculator calculates the propensities for the reactions by the specified kinetic laws.
There are three types of implementations of the Network interface:
-
Readers which can read network data from files (e.g. FernMLNetwork, SBMLNetwork)
-
Wrappers which redirect method calls to existing network classes (e.g. CytoscapeNetworkWrapper)
-
Evolution algorithms which create networks from scratch by certain rules (e.g. AutocatalyticNetwork)
For each network, stochastic simulations can be performed with all implemented simulation algorithms.
Import and export of networks
FERN supports two formats for loading and exporting networks: the SBML format [33] as well as the simpler but also XML based FernML format. For reading and writing the SBML format, FERN uses the Java bindings of the C library (libSBML) available at http://www.sbml.org. Thus, it can be easily adapted to new developments of the SBML format. Currently, SBML version 2 levels 1–3 are supported.
From the model loaded by libSBML from the SBML file, a FERN SBMLNetwork is created using the list of compartments, species, reactions, parameters and events in the model. Events have to be registered with a simulator by the SBMLNetwork if they are to be triggered during the simulation (see Figure 3 for an example). Triggering of events is handled by specific observers.
Currently, the SBMLNetwork class uses only the features of SBML necessary for the simulation of the network. It supports MathML to define complex reaction mechanisms but not rules, constraints or function definitions. If these features are required they can be incorporated easily by extending the SBMLNetwork class and loading these features from the SBML model created by libSBML. Since many systems biology applications support SBML (e.g. CellDesigner [34]), the SBML format can be used as an interchange format between FERN and these other applications.
SBML is a powerful format which can provide lots of information about a model. In contrast, FernML stores only the topology of the reaction network, optional annotations and the simulation parameters (see Additional file 2 and Additional file 3). This results in a much more simplified input format. More complex aspects, such as volume change due to cell growth and division, can then be modeled in Java using the FERN library in a straightforward way (see Results for an example). As a consequence, arbitrarily complex models can be designed.
Since FernML supports only the reaction rate equations used by Gillespie [5], the propensities can be recalculated at each step efficiently by a few arithmetic operations. SBML uses MathML to store the kinetics of a reaction. This allows for more complex reaction mechanisms and is particularly useful if the model cannot be formulated exclusively with first or higher order rate equations. To evaluate MathML expressions, FERN creates expression trees from them which have to be evaluated every time a propensity is calculated. Since this is one of the essential steps of SSAs, the simulation of an SBML network in FERN can be significantly slower than the simulation of the same network as a FernML network (see Figure 5). Thus, if only simple reaction rate equations are used, an SBML network should be converted to a FernML network using the provided conversion methods before performing the simulation.
FERN is not restricted to the input formats currently available. Any new input format can be easily included by implementing the Network interface or extending the AbstractNetworkImpl class.
Simulation algorithms
FERN provides implementations for three exact stochastic simulation algorithms, three state-of-the-art tau leaping procedures (see [8, 10]) and a hybrid method combining SSA and tau-leaping [19]. The exact SSAs implemented include the original direct method of Gillespie [5], the next reaction method of Gibson and Bruck [6] and an enhanced version of the direct method. This enhanced method uses the dependency graph technique of the next reaction method to only update the propensity functions which are affected by the firing of a reaction. Apart from this improvement, it is identical to the direct method.
The tau-leaping algorithms are all based on the modified tau-leaping procedure proposed by Cao et al. [9] which avoids the problem of negative populations observed for the original tau-leaping procedure. This method switches to an exact SSA (in our implementations the enhanced Gillespie) for a few steps if the selected τ becomes too small. The three implementations differ only in the way the error is bounded (see [10] for details). The error is bounded either by the sum of all propensity functions (TauLeapingAbsoluteBoundSimulator), the relative changes in the individual propensity functions (TauLeapingRelativeBoundSimulator) or the relative changes in the molecular populations (TauLeapingSpeciesPopulationBoundSimulator).
Furthermore, FERN implements the hybrid method by Puchalka and Kierzek [19] which partitions the system during the simulation into slow reactions which involve only small molecule numbers and fast reactions which involve large molecule numbers. The slow reactions are then simulated using an exact SSA while the fast reactions are simulated with tau-leaping. This algorithm was chosen over other hybrid methods for two reasons. First it uses only stochastic simulation algorithms, i.e. exact SSA and tau-leaping, and no further assumptions such as quasi-steady state. Second, the partitioning of the system is performed dynamically according to the state of the system and updated after each reaction step. Our implementation of the hybrid method uses our more efficient enhanced Gillespie algorithm (see Figure 5) instead of the Gibson and Bruck algorithm used by Puchalka and Kierzek. On the model of LacZ and LacY gene expression by Kierzek [30], the hybrid method speeds up the runtime by a factor of 98 compared to the enhanced Gillespie algorithm.
Each simulation algorithm is implemented by extending the abstract class Simulator or one of its subclasses. A simulation consists of the following steps (see also Additional file 3). First, the data structures are initialized and the simulation is started by passing a simulation controller implementing the SimulationController interface. The simulation controller decides after each step if the simulation can continue. The most basic one is the DefaultController which lets the simulation run until a given simulated time is reached.
In each step, the behavior of the simulator depends on the simulation algorithm implemented (see Figure 1). The direct and enhanced Gillespie algorithms draw the time interval τ till the next reaction from an exponential distribution. The reaction to be fired is then drawn with a probability proportional to its reaction propensity. For this purpose, a random variable r2 is first drawn from the uniform distribution between 0 and 1. The corresponding reaction μ is then identified via a linear search such that . The Gibson-Bruck algorithm generates τ
k
for each reaction R
k
and at each step fires the reaction with the minimum τ
k
obtained from a priority queue. Tau-leaping methods also choose a time interval τ but in this case several reactions can be fired during this interval (for more details see [9]).
After each step, the molecule numbers for reactants and products and the propensity functions for the reactions are recalculated. Here, reactants and products of a reaction can be identified efficiently from the adjacency list for this reaction stored in the network structure. Propensities are updated efficiently for all simulation algorithms apart from the original Gillespie algorithm using a dependency graph which stores for each reaction all reactions whose propensity is changed by the firing of this reaction.
Future developments of the algorithms can easily be included into FERN by extending one of the SSA implementations or the original Simulator class. In the same way ODE solvers or simulators for spatial models which are not provided by FERN can be integrated.
Observer system
FERN uses observers to trace the simulation progress and react to events. For this purpose, each observer has to implement functions which describe its response at specific time points of the simulations. Such responses may occur either at the beginning or the end of a simulation, before each step, after a reaction is fired or when a certain time is reached. In order to be notified of these events, observers have to be registered with the simulator.
Observer implementations are provided for tracing the molecule numbers for some species in arbitrary intervals, recording the firings of reactions, computing distributions of molecule numbers at a certain time over many simulation runs as well as many others. Several observers can be registered for a simulation at the same time and most of them can also handle repeated simulation runs, e.g. to create average curves or curves containing all trajectories for the individual simulation runs.
Visualization
In general, the observers use gnuplot to present their results. Once gnuplot is installed on a system and accessible e.g. via the path variable, the Gnuplot class makes it possible to easily create plots and retrieve them as Image objects, save them as files or present them online in a window. Plots can be customized using appropriate gnuplot commands.
Furthermore, FERN was used to implement Cytoscape [35] and CellDesigner [34] plugins for running and visualizing the simulations from within the Cytoscape or CellDesigner environments (for more details see Results).
Stochastics
An important feature of FERN is that random number generation is handled by the singleton class Stochastics. Accordingly, only one instance of this class is instantiated during a FERN run and all calls for random numbers are referred to this instance. This has several advantages. First, the underlying random number generator can be easily replaced if faster and better random number generators are developed. Currently, the Mersenne Twister implementation of the Colt Project is used http://dsd.lbl.gov/~hoschek/colt/. Second, by setting the seed value for the random number generator explicitly, the simulation can be made deterministic and e.g. interesting trajectories can be reproduced. Third, it is possible to count the number of random number generations necessary for different implementations of SSAs.