FERN - a Java framework for stochastic simulation and evaluation of reaction networks

Abstract

Background: Stochastic simulation can be used to illustrate the development of biological systems over time and the stochastic nature of these processes. Currently available programs for stochastic simulation, however, are limited in that they either a) do not provide the most efficient simulation algorithms and are difficult to extend, b) cannot be easily integrated into other applications or c) do not allow to monitor and intervene during the simulation process in an easy and intuitive way. Thus, in order to use stochastic simulation in innovative high-level modeling and analysis approaches more flexible tools are necessary.

Results: In this article, we present FERN (Framework for Evaluation of Reaction Networks), a Java framework for the efficient simulation of chemical reaction networks. FERN is subdivided into three layers for network representation, simulation and visualization of the simulation results each of which can be easily extended. It provides efficient and accurate state-of-the-art stochastic simulation algorithms for well-mixed chemical systems and a powerful observer system, which makes it possible to track and control the simulation progress on every level. To illustrate how FERN can be easily integrated into other systems biology applications, plugins to Cytoscape and CellDesigner are included. These plugins make it possible to run simulations and to observe the simulation progress in a reaction network in real-time from within the Cytoscape or CellDesigner environment.

Conclusion: FERN addresses shortcomings of currently available stochastic simulation programs in several ways. First, it provides a broad range of efficient and accurate algorithms both for exact and approximate stochastic simulation and a simple interface for extending to new algorithms. FERN's implementations are considerably faster than the C implementations of gillespie2 or the Java implementations of ISBJava. Second, it can be used in a straightforward way both as a stand-alone program and within new systems biology applications. Finally, complex scenarios requiring intervention during the simulation progress can be modelled easily with FERN.

Figure 1

Stochastic simulation. This figure shows the flow of one simulation step. On the left-hand side the flow for the original Gillespie algorithm can be seen. On the right-hand side, we illustrate how the different steps are modified by the Gibson-Bruck, enhanced Gillespie and tau-leaping algorithms. Here, U(0, 1) denotes the uniform distribution on the range of 0 to 1 and a_μ the reaction propensity for reaction μ.

Figure 2

FERN design. The figure illustrates the overall design of FERN into three layers. Each layer is represented by one interface or abstract class: Interface Network and abstract classes Simulator and Observer.

Figure 3

Example program. This figure shows a small example on how to use FERN for running a simulation on a network. First, a network is loaded from an SBML file and then a simulator is created. In the next step, an observer is created and registered with the simulator. In this example, the observer records the current amount of molecule X every second of simulated time. The SBML events are registered with the simulator and the simulation is started to run for 50 seconds. Finally, the recorded results for X are printed.

Figure 4

UML diagram of Network related classes. This figure shows an UML diagram for the Network interface and related interfaces and classes.

Figure 5

Runtime Comparisons. The EGF signaling pathway described by Lee et al. [38] was simulated with the three exact methods provided by FERN (original and enhanced Gillespie algorithm and the next reaction method by Gibson and Bruck) for a simulated time of 800 seconds both with an SBML network using expression trees to represent MathML expressions and a FernML network. For each combination of network type and stochastic simulation algorithm, 1,000 simulations were performed and the average runtime in milliseconds was calculated. The same simulations were performed with the Gillespie and Gibson-Bruck algorithms of ISBJava. All results were obtained on one processor of an Intel Core2Duo with 2.4 GHz. Standard errors in all cases were < 1.5 milliseconds.

Figure 6

Results for the LacZ model. Average results of 1,000 simulations are shown for the LacZ protein over ten bacterial generations (red). After each generation (2100 s) the number molecules for each species was divided by 2 to simulate cell division. The blue line shows a linear fit to the increasing LacZ concentration during the first generation. This yields a rate of protein synthesis of 21s^-1.