 Software
 Open Access
 Published:
MONALISA for stochastic simulations of Petri net models of biochemical systems
BMC Bioinformaticsvolume 16, Article number: 215 (2015)
The Erratum to this article has been published in BMC Bioinformatics 2015 16:371
Abstract
Background
The concept of Petri nets (PN) is widely used in systems biology and allows modeling of complex biochemical systems like metabolic systems, signal transduction pathways, and gene expression networks. In particular, PN allows the topological analysis based on structural properties, which is important and useful when quantitative (kinetic) data are incomplete or unknown. Knowing the kinetic parameters, the simulation of time evolution of such models can help to study the dynamic behavior of the underlying system. If the number of involved entities (molecules) is low, a stochastic simulation should be preferred against the classical deterministic approach of solving ordinary differential equations. The Stochastic Simulation Algorithm (SSA) is a common method for such simulations. The combination of the qualitative and semiquantitative PN modeling and stochastic analysis techniques provides a valuable approach in the field of systems biology.
Results
Here, we describe the implementation of stochastic analysis in a PN environment. We extended MonaLisa  an opensource software for creation, visualization and analysis of PN  by several stochastic simulation methods. The simulation module offers four simulation modes, among them the stochastic mode with constant firing rates and Gillespie’s algorithm as exact and approximate versions. The simulator is operated by a userfriendly graphical interface and accepts input data such as concentrations and reaction rate constants that are common parameters in the biological context. The key features of the simulation module are visualization of simulation, interactive plotting, export of results into a text file, mathematical expressions for describing simulation parameters, and up to 500 parallel simulations of the same parameter sets. To illustrate the method we discuss a model for insulin receptor recycling as case study.
Conclusions
We present a software that combines the modeling power of Petri nets with stochastic simulation of dynamic processes in a userfriendly environment supported by an intuitive graphical interface. The program offers a valuable alternative to modeling, using ordinary differential equations, especially when simulating singlecell experiments with low molecule counts. The ability to use mathematical expressions provides an additional flexibility in describing the simulation parameters. The opensource distribution allows further extensions by thirdparty developers. The software is crossplatform and is licensed under the Artistic License 2.0.
Background
Chemical and thus biochemical systems, can be modeled at different levels of abstraction. The choice and level of the method depends on the available data. With the development of the omicstechnologies we get many different types of experimental data, ranging from NGS (Next Generation Sequencing) data to image data. It is indispensable to check the quality and completeness of the data to decide which method should be applied for theoretical modeling. In this context, it is essential to provide mathematical formalisms which can handle and combine different levels of abstraction. This can then lead to assurance of several sound analysis techniques and an intuitive graphical interface. Here, we use Petri nets (see next subsection) — a formalism that exhibits all these properties. Stateoftheart tools implementing hybrid Petri nets include commercial tools, e.g., Cell Illustrator [1,2], and colored Petri nets, e.g., CPN [3]. Our aim is to provide an opensource, easily extendable tool for biochemical Petri nets that allow to model the classical P/T networks as well as the Gillespie’s stochastic simulation approach.
If there exists kinetic data in a sufficient amount we can start a quantitative modeling. Mathematical models of (bio)chemical reaction systems are usually formulated in terms of ordinary differential equations (ODE). The simulation of such models is performed numerically, and the result is strictly deterministic. This approach ignores stochastic fluctuations which are important for biological systems operating with low molecule numbers. The concept of Chemical Master Equations (CME) aims to describe stochastic fluctuation in reaction systems. For a more detailed introduction we refer to the textbook by Atkins [4].
The CME is a countable, but infinite, system of firstorder differential equations that determines time evolution of the probabilities of the discrete states of a system. A numerical solution of this equation is infeasible. Alternatively, the Monte Carlo simulation method can be applied to compute valid trajectories through the state space of the system. For this purpose, Gillespie [5] proposed an efficient stochastic simulation algorithm (SSA). For the initial state of a chemical system, the Gillespie’s method simulates the evolution of the number of molecules by estimating when and which reaction would be occurring next. Several implementations of the Gillespie’s method have already been proposed to improve the computational performance. The approximate SSA [6] facilitates a significant speedup factor.
Petri nets
Petri net (PN) is a powerful mathematical concept [7], which is widely applied for modeling systems of chemical reaction, metabolic pathways [8], signaling pathways [9], or gene expression networks [10]. The main idea is the consequent distinction between passive and active parts of the network, enabling for a sound treatment of concurrency. Here, we will give briefly the main definitions which are essential for understanding the idea described in our approach.
Definition 1 (Petri net).
A Petri net (PN) is a directed bipartite graph PN = (P, T, E, f(e), m(P _{ i })) with

P and T are disjoint sets, P∪T is the set of all vertices. P is the set of places, T is the set of transitions,

E ⊆((P×T)∪(T×P)) is the set of directed edges,

\(f(e): E \rightarrow \mathbb {N}_{+} \) is the weight function which assigns a positive integer weight to each edge e∈E and

\(m(P_{i}): P \rightarrow \mathbb {N}_{0}\) is the marking which assigns a positive integer number of tokens to each place P _{ i }∈P.
Places (states or species) represent passive and transitions (reactions) active elements. Movable objects are called tokens and are located on the places. A token distribution (marking) represents the number of entities of the corresponding species and defines a system state. The dynamics are modeled by the movements of tokens through firing of transitions, applying firing rules.
Definition 2.
(Pre and postplaces, firing of a transition)

∙T _{ i }:={P _{ j }∈P  (P _{ j }, T _{ i })∈E} is the set of the preplaces of transition T _{ i }, i.e., the set of all places which have an outgoing edge to T _{ i }.

T _{ i }∙:={P _{ j }∈P  (T _{ i }, P _{ j })∈E} is the set of the postplaces of transition T _{ i }, i.e., the set of all places which have an incoming edge from T _{ i }.

Transition T _{ i } is active and can fire iff ∀P _{ j } ∈∙T _{ i }:f(P _{ j },T _{ i })≤m(P _{ j }).
Pre and posttransitions are defined analogously to pre and postplaces.
A transition can only fire, if the number of tokens on all preplaces of the transition is equal to or greater than the weights of the edges between the corresponding places and the transition. In this case, the transition is active or enabled. Firing of a transition consumes tokens from the preplaces and adds tokens to the postplaces, according to the weights of the corresponding edges. We extend the standard definition of PN by introducing constant places. The number of tokens on a constant place is not affected by the firing of transitions but is determined by a userdefined mathematical expression. For instance, pulsatile secretion of insulin in a model of insulin signaling can be described as
where the variable Time stands for the simulated time (in seconds).
In the context of (bio)chemical systems, places represent chemical compounds (e.g. insulin, insulin receptor), complexes (e.g. receptorligand complex) or different states (e.g. inactive or active receptor), and transitions describe chemical reactions. The number of tokens on a place represents the number of molecules of the corresponding biochemical species. In signaling or gene expression PN, the number of tokens can describe the activation state (on/off) or the strength of a response. An overview of the PN formalism and its application to biology has already been earlier described [11,12].
Available software
MONALISA [13] is an opensource software for creation, visualization and analysis of PN. MONALISA implements several analysis techniques such as invariant analysis (implemented in C, [14]), topology properties, knockout analysis, and other, supporting a broad range of file formats like PNML, PNT, SPEED, SBML (read: all levels and versions of SBML Core, write: SBML core Level 2, Version 4), KGML and DAT. MONALISA provides the possibility to easily extend its functionality by new modules. It also supports the annotation with MIRIAM identifiers and SBO terms. An interface to SEDML [15] is planned.
Here, we introduce a module which extends the functionality to applications of dynamic simulations. The applied firing rules determine the strategy to choose the next transition to fire and the time point at which the event takes place. Stochastic firing rules simulate the dynamics of the stochastic kinetic of a mass action reaction system inside a cell. Alternative tools for the stochastic analysis of biochemical PN are Snoopy [16] and VANTED [17]. VANTED is an opensource solution but without any simulation module. Snoopy is not opensource and aims at providing the entire types of PN mainly for the Petri net community. Modeling for biology is one aspect among many others. MonaLisa is focused on biological applications and covers another spectrum than Snoopy, ranging from general network analysis to specific decomposition methods at steadystate as well as nonsteadystate conditions. Moreover, MONALISA includes a simulation mode with multiple parallel simulations. For a more detailed comparison of MONALISA and Snoopy see Table 1. The table compares the different features of both tools, like existing simulation modes, analysis techniques, availability, and supported file formats.
StochPy [18] is another tool for stochastic simulation of biological processes, but is not PNbased. It uses the PySCeS [19] model description language — a textbased model description technique. It provides an integrated statistical output analysis (autocorrelations, propensities, moments, waiting times) and can be easily extended due to its integration with other Python libraries. Implementations of the exact and the approximate SSA are offered, but StochPy provides neither an intuitive graphical user interface nor any PN and topological analysis method.
Implementation
MONALISA is an opensource software (see Additional file 1) distributed under the Artistic License 2.0. It is accessible to modifications and allows adjustment and extension of its functionality by experienced users. A description of all features of the simulation module is given in the documentation file MonaLisa_Documentation.pdf (Additional file 2) provided in the supplementary.
Libraries
The simulation module is implemented in Java, version 1.7, as a plugin of MONALISA. The plotting functionality exploits the JFreeChartlibrary [20]. The integration of the library exp4j [21] develops the option to evaluate userdefined mathematical expressions. The Java class HighQualityRandom [22] implements the standard pseudo random number generator (RNG) Ran defined by Press et al. [23]. This RNG is a good compromise between speed and cryptographic randomness and appropriate for all algorithms described here. For a more detailed discussion of the choice of RNGs for simulations we refer to the textbook of Knuth [24] and the numerous literature published on this topic, see, e.g., [2529].
Modes
PN modeling addresses a broad range of scenarios in systems biology with various goals and requirements. To account for the different needs, we implemented four simulation modes in MONALISA, ranging from simple synchronous and asynchronous modes to the wellestablished Gillespie’s algorithm for stochastic simulation:

1.
Asynchronous: One randomly chosen transition fires per simulation step, without any time consumption. All active transitions have the same probability to fire.

2.
Synchronous: Multiple active transitions fire simultaneously per step. By default, the simulator tries to fire all active transitions at once. Transitions which share preplaces are called concurrent and compete for the tokens of the shared preplaces. If the number of tokens on such a place is not sufficient for all enabled posttransitions, transitions to fire are randomly chosen with equal probability until all tokens are consumed. The simulator also allows to select the fraction of active transitions to fire. If the fraction is less than 100%, transitions are randomly selected with equal probability in each step.

3.
Stochastic: An enabled transition T _{ i } has to wait a defined time dt before it can fire. The waiting time dt is simulated for each transition as an exponential distributed random variable E x p(λ), a detailed description is available in the textbook of D.J Wilkinson [30]. The parameter λ is the firing rate specified by the modeler for the particular transition. Among all enabled transitions only the transition with the lowest waiting time fires. If several transitions have the lowest waiting time, one of them is selected randomly. To speed up the simulation, waiting times are recomputed only for transitions with changed number of tokens on the preplaces. Waiting times of the posttransitions of constant places are recomputed in every step.

4.
Stochastic simulation algorithm: This algorithm implements the exact and approximate SSA. The exact SSA implements the direct method of Gillespie [5,30], applying an internal data structure similar to the “dependency graph” introduced by Gibson and Bruck [31] for better performance.
The exact SSA consists of four steps:

i
Computing the rates of all reactions based on the stochastic rate constants c(T _{ i }) and the number of reactant molecules. We refer to the textbook of Wilkinson [30] for a detailed description.

ii
The next firing time dt is computed by
$$dt = \frac{ln(1U_{1})}{\text{sum of reaction rates}} $$from a uniformly distributed (in the interval [0,1]) random number U _{1}.

iii
A reaction i is chosen such that the following equation is satisfied
$$\sum_{j=1}^{i1}r(T_{j}) < U_{2} \ \text{sum of all rates} \leq \sum_{j=1}^{i}r(T_{j}) $$for a second uniformly distributed random number U _{2}.

iv
The number of molecules is updated according to the chosen reaction i and the time is increased by dt.
The approximate SSA integrates the τleaping algorithm of Gillespie and Petzold [32,33]. In brief, the algorithm chooses a time interval τ and decides how many times each reaction will occur in this period. The time interval should be small enough so that the reaction rates do not change significantly during this interval. To avoid negative populations, reactions that do not have sufficient reactants to fire at least 20 times are considered as critical and are simulated in an exact way only.
At the beginning of each simulation step critical reactions are determined. The algorithm computes two firing times – one for the noncritical reactions τ _{1} according to [33] and a second time τ _{2} for the critical reactions.
If τ _{1}<τ _{2}, no critical reaction occurs. For each noncritical reaction, the number of its occurrences is generated as a random Poisson variable with the mean value r(T _{ i }) τ _{1}, and the reaction is executed according to the chosen number of times. If τ _{2}<τ _{1}, one critical reaction is chosen like in the exact SSA and the number of occurrences of the noncritical reactions are chosen as for the first case.

i
Converting input data
In contrast to ODEbased approaches, the stochastic methods operate with numbers of molecules instead of concentrations. Therefore, the continuous biological parameters like concentrations of chemical compounds and the reaction rate constants must be converted to moleculebased data. So, the user has to provide the volume of the simulated environment as an additional parameter, and a preprocessor converts concentrations and mass action reaction rate constants to the appropriate numbers of molecules and stochastic rate constants, respectively:
where n _{ X } is the number of molecules of compound X, [X] the concentration of compound X, N _{ A }=6·10^{23} the Avogadro constant and V the reaction volume.
The conversion of the mass action rate constant k to the stochastic rate constant c depends on the order of the reaction. The mass action reaction rate constant is given in units of M s ^{−1} for a zero order reaction, in units of s ^{−1} for a first order reaction, in units of M ^{−1} s ^{−1} for a second order reaction, and in general, in units of M ^{−(x−1)} s ^{−1} for a reaction of the order of x. The mass action reaction constant k gives the reaction rate for standardized concentrations of 1 M (i.e., N _{ A } molecules in 1 l reaction volume). To get the reaction rate constant for the number of molecules we have to divide k by (V N _{ A })^{(x−1)} and get
A substance (preplace) P _{ j } may contribute with high stoichiometric order (i.e., weight f(P _{ j },T _{ i })) to the reaction T _{ i }. The combinatorial factor
accounts for the stoichiometric factors of the species P∈∙T _{ i }. Finally, we get the general formula
This equation can be simplified to:

c=k V N _{ A } for zeroorder reaction,

c=k for firstorder reaction,

\(c = \frac {k}{V \ N_{A}}\) for secondorder reaction of the form A+B→C, and

\(c = \frac {2 k}{V \ N_{A}}\) for secondorder reaction of the form 2A→B.
Mathematical expressions
Reaction rate constants and the number of tokens on constant places can be defined by mathematical expressions which may contain names of nonconstant places and the simulated time as variables.
Supported syntax
The evaluation of a mathematical expression is based on the free Java library exp4j [21]. It supports numerical input in standard and scientific notations. Supported operators and functions are listed in Table 2.
Conditional expression
It is possible to use multiple conditional expressions to define the reaction rate constant or the number of tokens on a constant place. The syntax has the form
Each condition is composed of the three parts:
with mathematical expressions [expression 1] and [expression 2]. The operator ([operator]) is one of “ =”, “ <”, “ >”, “ <=”, or “ >=”, describing comparisons. Multiple conditional expressions can be separated by a semicolon. The value of the first conditional expression whose conditions are satisfied is returned. If none of the described cases can be satisfied, 0 is returned. For example, the number of insulin molecules should be 1000 for the first 5 minutes, 100 for the next 5 minutes and 0 afterwards:
As long as the reaction time is below 300, the first conditional expression is valid, and the value 1000 is returned. If the reaction time exceeds 5·60, the second conditional expression becomes valid, and the number of insulin molecules is set to 100. After the reaction time of 10·60 seconds the number of insulin molecules drops down to zero.
Mathematical expressions are evaluated after each simulation step. In combination with timedependent expressions the Gillespie algorithm may lead to discretization artifacts in the case of the combination of a very large waiting time with a fast variation of the conditional expressions on a similar time scale. The discretization artifacts are negligible for standard applications. In further versions of MONALISA, we will implement the rejection method [34] to provide a more robust algorithm for this non–homogeneous Poisson process.
Results and discussion
MONALISA provides a PN editor, including visualization and a special graphical user interface for the simulation mode (see Figure 1). The simulation process is visualized by showing the current number of tokens on the places, highlighting with different colors active transitions and transitions which fired last. The user may change interactively the number of tokens on places and select manually a transition to fire.
Four different simulation modes can be chosen  asynchronous, synchronous, stochastic and the SSA. Various useful features are available:

History of fired steps,

Statistic information (number of fired steps, number of firings each transition has performed),

Saving and loading of marking,

Saving and loading of the simulation set–ups, including the marking, firing rates, reaction rate constants, mathematical expressions, and other information in XML format,

Plot of simulation results,

Export of simulation results into a tabstop separated file,

Application of mathematical expressions to describe rate constants and the number of tokens.
The concept of constant places is an assistant feature of MONALISA. Constant places are ignored in the firing step, though they contribute to the transitions’ activation states. Constant places are useful for modeling boundary conditions, external factors or inhibitory effects.
For biochemical systems, the stochastic and the Gillespie’s algorithm (SSA) are the most interesting simulation modes. Each stochastic simulation of identical initial parameters will have a different outcome. To reproduce a simulation, the user may want to set the seed of the Pseudo Random Number Generator (PRNG). Repetition of simulations with a fixed seed will always produce identical results.
Apart from the exact SSA, the simulator implements the approximate SSA with improved performance for systems with fast reactions or large molecule counts [6]. Up to 500 parallel simulations of identical initial parameters can be performed and take full advantage of multicore processors.
Case study of insulin receptor (IR) recycling
We want to explain modeling and simulating of the complex biochemical system of insulin receptor activation and recycling. This case study intends to demonstrate a possible workflow but not a thorough analysis of the IR system.
Insulin is an important hormone that regulates various cellular processes, among others the intake of glucose by such tissues as fat, muscle and brain or glucose production by the liver. The response of a cell to an elevated insulin level is mediated by the insulin receptor, which is mainly located in the cell surface membranes. Impairments in the key components of the IR regulatory system can cause diseases such as the metabolic syndrome or type 2 Diabetes mellitus.
Figure 2 gives an overview of the processes in the model. Figure 3 depicts the corresponding PN. It consists of 6 places and 11 transitions. The unbound inactive insulin receptor is modeled by the place IR, the free insulin by the place Insulin. We consider only a single cell. The insulin concentration is regulated by external processes. Consequently, we choose a constant place to model the timedependent external insulin level. Insulin can bind (transition Bind_Insulin) to the free receptor. The receptorligand complex (place IR_I) can dissociate again (transition Diss_Insulin). Alternatively to the binding of insulin, the free receptor can be internalized to the cytoplasm (transition Inter_IR) and is added to the internal receptor pool (place IR_Int).
As soon as insulin is bound to the receptor, the receptorligand complex can be phosphorylated and becomes activated (transition Phos_IR_I). The bound phosphorylated receptor is represented by the place IR_I_P. The phosphorylation is reversible, the reverse action (transition Dephos_IR_I_P) correlates with the dissociation of insulin. Alternatively, the activated receptorligand complex can be internalized (transition Inter_IR_I_P). The phosphorylated complex in cytosol (place IR_I_P_Int) can be either transported back to the membrane (transition Deinter_IR_I_P) or deactivated in the cell (transition Dephos_IR_I_P). The inactive unbound receptor in the cytosol can be transported back to the membrane (transition Deinter_IR) or be destroyed by degradation (transition Degradation). A process of synthesis sustains the level of IR.
Simulation of the model
The model is simulated using the Gillespie’s algorithm. Model parameters (initial concentrations of IR and reaction rate constants) are adapted from [35]. The volume of the system is set to 1E9 l. Initial concentrations of the places are listed in Table 3 and the mass action reaction rate constants in Table 4. The rate of IR synthesis depends on the concentration of the receptor in the cytosol. If the concentration decreases below the steadystate value of 1E13 M, an increased synthesis rate is used. This is described by the expression
We simulated two days with a typical 24 h insulin profile of a healthy person, whereat the mean basal insulin concentration is about 6·10^{−11} M [36]. Meal intake, simulated at 09:00 h, 13:00 h and 18:00 h, triggers the rapid increase of the insulin concentration to approximately 3.6·10^{−10} M [36]. The hormone concentration returns to the basal level within the next 3 hours. This profile is described by a mathematical expression with four cases (one conditionfree case for the basal concentration and three timedependent cases for the meal intakes) and a negative exponential function which describes the decrease of insulin concentration after meal intake. The expression can be found in the supplementary text file “example_MathExp” (Additional file 3).
Results of a simulation of two days are plotted in Figure 4. The figure shows on top the number of insulin molecules plotted against time. Whereas the number of molecules of different IR states is depicted at the bottom. The PN (Additional file 4) is provided as the MONALISAproject file IR_Model.mlproject. The simulation setup file “IR_Model_params.xml”, (Additional file 5) can be found in the supplement.
Conclusions
MONALISA combines the powerful modeling concept of Petri net formalism with stochastic simulation abilities implementing a simulation module. The intuitive graphical user interface allows to focus on modeling and simulation. Different simulation modes make the software suitable for working with PN in various areas, especially for biological tasks like analysis and simulation of metabolic systems, gene regulatory networks or signal transduction pathways. Constant places and the usage of mathematical expressions for describing simulation parameters offer an improved flexibility compared to other tools, allowing for modeling of nonstandard kinetics and complex relationships to the external environment. Useful features, like builtin plotting, navigation through the simulation history, export of the simulation setups to XMLfiles and setting the seed of the random number generator, can help to easily adjust parameters to follow and assess simulation results. We demonstrate a small case study that models the prominent motif of the recycling of the insulin receptor. The case study presents the abilities of MONALISA to model and simulate complex biochemical systems and demonstrates the applicability of the software for studying biological processes.
Availability and requirements
Project name: Simulation mode for MONALISA  a tool for development, visualization and analysis of Petri netsProject home page: http://www.bioinformatik.unifrankfurt.de/tools/monalisa/ Operating system(s): Linux, Windows, MacOSProgramming language: Java, C (for invariantscomputation)Other requirements: Java 1.7 or higherLicense: Artistic License 2.0Any restrictions to use by nonacademics: None
Abbreviations
 ODE:

Ordinary differential equation
 PN:

Petri net
 CME:

Chemical master equation
 SSA:

Stochastic simulation algorithm
 PRNG:

Pseudo random number generator
 IR:

Insulin receptor
 r _{ i } :

Rate of reaction i
 k _{ i } :

Deterministic rate constant of reaction i
 c _{ i } :

Stochastic rate constant of reaction i
 M:

Molar concentration = moles per liter
 V :

Volume
 l:

Liter
 n _{ X } :

The number of molecules of compound X
 [X]:

Concentration of the compound X (or of the compound represented by the place X)
 N _{ A } :

Avogadro constant = 6·10^{23}
References
 1
Nagasaki M, Doi A, Matsuno H, Miyano S. Genomic Object Net: I. A platform for modelling and simulating biopathways. Appl Bioinform. 2003; 2(3):181–184.
 2
Nagasaki M, Saito A, Jeong E, Li C, Kojima K, Ikeda E, et al.Cell Illustrator 4.0: A computational platform for systems biology. In silico Biol. 2010; 10(1):5–26.
 3
Ratzer A, Wells L, Lassen H, Laursen M, Frank J, Stissing M, et al.CPN Tools for editing, simulating, and analysing coloured Petri nets. In: Applications and Theory of Petri Nets 2003: 24th International Conference, ICATPN 2003. Heidelberg, New York: Springer: 2003. p. 450–62.
 4
Atkins P, De Paula J, Vol. Bd. 1. Atkins’ Physical Chemistry. Cambridge, UK: Oxford University Press; 2002.
 5
Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977; 81(25):2340–361.
 6
Gillespie DT. Approximate accelerated stochastic simulation of chemically reacting systems. J Chem Phys. 2001; 115:1716.
 7
Reisig W. Petrinetze: Eine Einführung, 2nd edn. Berlin Heidelberg: In German, Springer Verlag; 1986.
 8
Junker B, Klukas C, Schreiber F. VANTED: A system for advanced data analysis and visualization in the context of biological networks. BMC Bioinf. 2006; 7(1):109.
 9
Sackmann A, Heiner M, Koch I. Application of Petri net based analysis techniques to signal transduction pathways. BMC Bioinf. 2006; 7(1):482.
 10
Grunwald S, Speer A, Ackermann J, Koch I. Petri net modelling of gene regulation of the Duchenne muscular dystrophy. Biosystems. 2008; 92(2):189–205.
 11
Koch I, Heiner M. Petri nets. In: Analysis of Biological Networks. Weinheim: John Wiley & Sons, Inc.: 2007. p. 139–179.
 12
Koch I, Reisig W, Schreiber F. Modeling in Systems Biology: The Petri Net Approach. Heidelberg, New York: Springer; 2010.
 13
Einloft J, Ackermann J, Nöthen J, Koch I. MonaLisavisualization and analysis of functional modules in biochemical networks. Bioinformatics. 2013; 29(11):1469–1470.
 14
Koch I, Ackermann J. On functional module detection in metabolic networks. Metabolites. 2013; 3(3):673–700.
 15
Waltemath D, Adams R, Bergmann F, Hucka M, Kolpakov F, Miller A, et al.Reproducible computational biology experiments with SEDML  The Simulation Experiment Description Markup Language. BMC Syst Biol. 2011; 5(1):198.
 16
Heiner M, Herajy M, Liu F, Rohr C, Schwarick M. Snoopy – A Unifying Petri Net Tool. In: Application and Theory of Petri Nets. Heidelberg: Springer: 2012. p. 398–407.
 17
Junker BH, Klukas C, Schreiber F. VANTED: A system for advanced data analysis and visualization in the context of biological networks. BMC Bioinf. 2006; 7(1):109.
 18
Maarleveld TR, Olivier BG, Bruggeman FJ. StochPy: A Comprehensive, UserFriendly Tool for Simulating Stochastic Biological Processes. PLoS ONE. 2013; 8(11):79345.
 19
Olivier BG, Rohwer JM, Hofmeyr JHS. Modelling cellular systems with PySCeS. Bioinformatics. 2005; 21(4):560–1.
 20
JFreeChart. http://www.jfree.org/jfreechart/. Accessed 21.05.2015.
 21
Asseg F. Exp4j. http://www.objecthunter.net/exp4j/index.html. Accessed 21.05.2015.
 22
Javamex – A Java implementation of the Numerical Recipies Random Number generator. http://www.javamex.com/tutorials/random_numbers/. Accessed 21.05.2015.
 23
Press WH. Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge: Cambridge University Press; 2007.
 24
Knuth DE. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. Menlo Park California: AddisonWesley Publishing Company; 1997.
 25
Marsaglia G. Random numbers fall mainly in the planes. Proc Natl Acad Sci U S A. 1968; 61(1):25.
 26
Marsaglia G. Diehard: A Battery of Tests of Randomness. http://www.stat.fsu.edu/pub/diehard/. Accessed 21.05.2015.
 27
Ackermann J, Tangen U, Bödekker B, Breyer J, Stoll E, McCaskill J. Parallel random number generator for inexpensive configurable hardware cells. Comput Phys Commun. 2001; 140(3):293–302.
 28
Brown RG, Eddelbuettel D, Bauer D. Dieharder: A Random Number Test Suite. http://www.phy.duke.edu/~rgb/General/dieharder.php/. Accessed 21.05.2015.
 29
Pawlikowski K, Jeong HD, Lee JS. On credibility of simulation studies of telecommunication networks. Commun Mag IEEE. 2002; 40(1):132–9.
 30
Wilkinson DJ. Stochastic Modelling for Systems Biology. Boca Raton: CRC Press; 2011.
 31
Gibson MA, Bruck J. Efficient exact stochastic simulation of chemical systems with many species and many channels. J Phys Chem A. 2000; 104(9):1876–89.
 32
Cao Y, Gillespie DT, Petzold LR. Avoiding negative populations in explicit Poisson tauleaping. J Chem Phys. 2005; 123(5):054104.
 33
Gillespie DT, Petzold LR. Improved leapsize selection for accelerated stochastic simulation. J Chem Phys. 2003; 119(16):8229–234.
 34
Lewis PA, Shedler GS. Simulation of nonhomogeneous Poisson processes by thinning. Naval Res Logist Q. 1979; 26(3):403–413.
 35
Sedaghat AR, Sherman A, Quon MJ. A mathematical model of metabolic insulin signaling pathways. Am J Physiol Endocrinol Metab. 2002; 283(5):1084–1101.
 36
Polonsky K, Given B, Van Cauter E. Twentyfourhour profiles and pulsatile patterns of insulin secretion in normal and obese subjects. J Clin Investig. 1988; 81(2):442.
Acknowledgments
PB would like to thank Sanofi Aventis Deutschland GmbH for financial support as well as Jennifer Scheidel and Leonie Amstein for testing the simulation module and providing valuable feedback. This work was partly supported by the LOEWE program of the State of Hesse (Germany) and by the Cluster of Excellence Frankfurt Macromolecular Complexes.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
PB developed the IR model and the simulation module and drafted the manuscript. JE participated in the implementation of the simulator. KL and JA contributed to the development of the IR model. IK has supervised the project. All authors read and approved the final manuscript.
An erratum to this article is available at http://dx.doi.org/10.1186/s1285901507257.
Additional files
Additional file 1
Current version of MONALISA. The current version of MonaLisa, including the implemented simulation module, is provided as executable.jarfile as well as source files. To execute the compiled file, open a console (command line under Windows or terminal under UNIXlike systems), navigate to the directory where the MonaLisa.jar is located and type java jar MonaLisa.jar.
Additional file 2
Documentation of the Simulation Mode of MONALISA. Documentation of the simulation mode.
Additional file 3
Mathematical expression of an insulin profile. The text file “example_MathExp” contains the mathematical expression which is used for modeling a 24h profile of insulin concentration in healthy patients with three meal intakes.
Additional file 4
MONALISA project of the insulin receptor recycling model. The Petri net of the IR recycling model is provided as the MonaLisaproject file “IR_Model.mlproject”. To open it, start MonaLisa and select “Open project”.
Additional file 5
Simulation setup. Simulation setup, which includes initial concentrations of compounds and reaction rate constants, is stored in the XMLfile “IR_Model_params.xml”. To load it, open the IR_Model  project file, start the “Gillespie SSA”mode from the “Simulation”tab and click the “Load setup”button.
Rights and permissions
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Stochastic simulation algorithm
 Gillespie’s algorithm
 Petri net
 MonaLisa
 Insulin receptor recycling