MONALISA for stochastic simulations of Petri net models of biochemical systems

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 semi-quantitative 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 open-source 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 user-friendly 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 user-friendly environment supported by an intuitive graphical interface. The program offers a valuable alternative to modeling, using ordinary differential equations, especially when simulating single-cell experiments with low molecule counts. The ability to use mathematical expressions provides an additional flexibility in describing the simulation parameters. The open-source distribution allows further extensions by third-party developers. The software is cross-platform and is licensed under the Artistic License 2.0. Electronic supplementary material The online version of this article (doi:10.1186/s12859-015-0596-y) contains supplementary material, which is available to authorized users.

Now, select a simulation mode from the same-named drop-down menu and click on the "Start simulation"-button. There are four simulation modes implemented, each of which evolves the system in its own way. Simulation modes will be described later on (section 2). If you want to simulate a system of chemical reactions, select "Gillespie SSA" (described in section 2.4).
When the simulation mode becomes activated, it is impossible to edit the net structure. You can adjust font, icon, arrow and edge size with the corresponding field under the "Start simulation"-button. To return to the editing mode, press "End simulation"-button. It will deactivate all simulation controls and stop any running simulation. In the following sections common functionality of simulation mode will be explained using the "Asynchronous"-mode.  One can fire an active transition by clicking on it with the left mouse button. Selecting a place and opening the context menu with the right mouse button allows to input the number of tokens for the place or flagging the place as constant (constant places are described in the next section).

Constant places
One specific feature of the simulation mode in MonaLisa are the constant places. In contrast to "normal" places, the number of tokens on a constant place can not be changed by firing of transitions. According to the standard rules of Petri net, firing of a transition consumes tokens on the pre-places of the transition (the number of consumed tokens is given by the weight of the arc between the pre-place and the transition) and adds tokens to the postplaces, according to the arc weight. The constant places, however, are ignored in the firing step, although they are still considered while determining which transitions are active. This feature can be useful when modeling input conditions which are not affected by the modeled system. For example, in the model of insulin receptor recycling (the model is briefly described in section 4) insulin is modeled as a constant place. It is assumed that insulin concentration is maintained and regulated by processes which are not modeled, and a single cell (which is modeled with the net) does not affect its global concentration.
Furthermore, the number of tokens on a constant place can be described by a complex mathematical expression. For instance, pulsatile insulin secretion can be described as where the variable Time stands for the simulated time (in seconds). For more information how to use mathematical expressions, see section "Using mathematical expressions" 1.10.1.
If a non-constant place carries tokens and is selected to be constant, it becomes a mathematical expression with the value of the token numbers assigned. If a constant place with a mathematical expression is selected to be non-constant, its expression is evaluated and the number of tokens is assigned to the non-constant place.

Performing simulation steps
As mentioned above, you can fire an active transition by clicking on it with the left mouse button. The other way of performing simulations is by specifying the number of simulation steps and letting the algorithm perform them by the selected rules. The controls are located in the area below the "Start-" and the "End simulation"-buttons (figure 3). On top of this area is the name of the selected simulation mode. You can specify the number of simulation steps in the text area "Steps to perform". Alternatively, you can enable the check box "Continuous mode" if the simulation should go on until no active transitions left or is interrupted otherwise.
After the number of steps is selected, click on the "Start simulation sequence"-button. This will start firing sequence of the specified number of steps.
The button turns into the "Stop simulation sequence"-button, clicking on which terminates the ongoing simulation sequence.

History and snapshots
The software tracks the simulation in a history list. The area "History of performed steps" (figure 4) shows which transitions have fired in the last 100 performed steps. The step performed last is highlighted with red background.
If the list is full (shows 100 steps) and a new step is performed, it is written to the end of the list, and the first element is removed. You can navigate through the history by selecting a step with the left mouse button or clicking the buttons " Step back" and " Step forward". The visual output in the NetViewer will be updated to the system's state at the selected step. Firing transitions which are not part of the history (e.g. by selecting a transition from the NetViewer or starting a new simulation sequence) will erase all history steps following the last selected step.
By default, a snapshot of the system is taken after each 100 steps. A snapshot contains the system's state at the corresponding step and the history of the last 100 steps. By selecting a snapshot from the "Snapshots:"-list, its history is loaded into the history area. The consequence is that though the history can only track the last 100 steps, selecting snapshots allows the user to have a complete history of the simulation, giving the possibility to return to every single state occurred during the simulation. Though the shown history can be erased by firing transitions (e.g. if you select the second step from the history and fires a new transition, steps 3-100 will be erased), it does not affect the snapshots.
Snapshots provide a detailed insight into the simulation sequence on the cost of increased memory consumption. For large simulations (several hundreds or thousands of steps) and low memory configurations, it may be advisable to turn snapshots off in the "Preferences"-panel (section 1.9). You can also turn snapshots on and off between the simulation sequences.

Custom markings and simulation setups
A state of the simulated system is represented by its marking. The program provides a useful option of saving and loading markings (figure 5). You can save the current state at any point by clicking on the "Save marking"-button.
This marking can then be restored by selecting it from the drop-down menu.
The currently selected marking can be removed with the "Delete marking"button. By default, each net is provided with the "empty marking" (no tokens are allocated on the places). Upon saving a marking, mathematical expressions of constant places are evaluated, and the values are stored as numbers. When restoring a marking by selecting it from the menu, only values of non-constant places are adopted, those of constant places are ignored. If you want to load  The option of saving more information about the simulation setup, such as mathematical expressions for constant places or firing rates of transitions, is provided by creating an XML file with setup parameters by clicking on the "Save setup"-button. To load a valid simulation setup, click on the "Load setup"button and choose a file which stores the settings. The structure of the setup files is described in section 3.1. The program tries to load all possible settings from the selected file. The settings of the places and transitions are mapped to the IDs, so it is possible to rename the nodes or create new ones without losing the compatibility with a setup file. The information stored in the setup file depends on the simulation mode and the supported features. Unsupported features are simply ignored upon loading.

Plotting results
The simulation mode of MonaLisa provides plotting functionality with features like zooming or exporting plots to PNG-format files. To plot the results of your simulation, click on the "Show plot"-button (see figure 6). A window called "Results plot" shows up. It plots the number of tokens on the places against the simulated time. It is possible to zoom in or out and to export the plot into a PNG file by right clicking on the image and opening the context menu. You can also zoom in a region by selecting it with the mouse or using the mouse wheel.
Storing results for plotting is memory-consuming, so it might be advisable to disable this feature for long simulations if it is not needed. This can be done in the "Preferences"-window by un-checking the "Built-in plotting"-check box. You can also select the places which results should be plotted with the "Places to plot"-button.

Preferences
To adjust advanced simulation parameters, click on the "Preferences"-button.
This will open a new window with control elements (see figure 8).

Writing simulation results into a text file
To write simulation results into a file, select the "Create log file"-check box. You can now select the destination path of the file by entering it in the text field below or by clicking the "Browse"-button. A text file is created upon first transition firing. For each simulation instance (initiated by starting of a simulation mode), a separate file is created with the name "Sim log -[currentDateAndTime].csv".
The file consists of elements separated by tab stops and line breaks. The first row contains the names of the places. The second row contains the IDs of the nodes of the corresponding places. The third row is a separation placeholder "----". From the fifth row on, the first element is the number of simulated steps, the second element is the simulated time, the third are the names of the transitions fired in this step (separated by a semicolon), and the following elements are the number of tokens on the places. Be aware that elements are separated by tabs only and not by semicolons, since some spreadsheet programs will try to use semicolons as separators, too.

Snapshots and other preferences
You can activate or deactivate snapshots (see section 1.5) by checking or unchecking the "Save snapshots"-check box.
The area below contains further parameters provided by the particular simulation mode. Common parameters are the delay between firings and the update interval. In the "Delay (ms)"-text field you can specify how long the simulator will wait before performing the next step. Zero stands for no delay, i.e. the next firing will occur as soon as the previous was processed. You can use delays for better visual tracking of the simulation progress.   Names of non-constant places are listed in the area to the right and can be inserted into the expression by clicking on them. You can also use the simulated time by selecting the variable "Time". For the simulation modes "synchronous" and "asynchronous", the time equals to the number of simulated steps. For timed PN simulation in the "stochastic" mode, the time is measured in arbitrary units, for the "Gillespie SSA" the units are seconds.
An expression must know which variables it uses to be able to map the names of the places to the number of tokens they carry. Places used as variables are listed in the "Used variables"-list located below the expression-area and are checked. Every time the user selects a place from the right list and inserts it into the expression, a corresponding check box becomes activated. If you enter the name of the place without selecting it from the list, you have to enable it in the "Used variables"-list, otherwise the expression will not be evaluated.
It is, however, strictly advised only to select places which are actually used in the expression and uncheck those which are removed from it. The reason is that the more variables an expression uses, the more expensive it can become to compute a simulation step. Furthermore, it is not allowed to use constant places as variables to prevent an unsolvable reference cycle.

Supported syntax
The evaluation of a mathematical expression is based on the free Java library exp4j [2]. It supports the numerical input in standard and scientific notations (1*10^2 = 1E2) and various operators and functions, which are listed in table 1. Additionally, the variable "pi" can be used for the number π = 3.14159 . . . . where [condition] is composed of the three parts: The value of the first conditional expression whose conditions hold is returned.
If a case has no conditions, it is automatically valid (e.g. "5+3" is a simple case without conditions). If none of the described cases is valid, 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: if Time <= (5*60) then 1000; if Time <= (10*60) then 100.
As long as the seconds' count stays below 300, the first case holds and 1000 is returned. If the number of simulated seconds is greater than 5*60, but less than 10*60, the second case is evaluated. After that, no case holds and 0 is returned.

Different simulation modes
This section describes the four simulation modes which are implemented in MonaLisa in detail.

Asynchronous
Asynchronous is the common mode of a PN simulation. In this mode, one transition is fired per simulation step, and the firing does not consume any time. All active transitions have the same probability to be chosen to fire.
The setup file stores the information about the number of tokens on nonconstant places and the mathematical expressions of the constant onces. The "Time"-column in the log file represents the number of simulated steps.

Synchronous
In the synchronous mode, multiple active transitions can fire in one step simultaneously. By default, the simulator tries to fire all active transitions at once.

Stochastic
Stochastic PN are a specific form of transition-timed PN. In a transition-timed net, a transition must wait a defined time before it can fire after it was activated.
Generally, the units of time are not specified but can be modeled as seconds or other common time units. This technique allows to model systems with time-consuming processes, therefore timed PN are widely used for modeling biochemical reactions.
In stochastic PN, waiting times dt are chosen stochastically from the probability density function with the firing rate r(T i ) of transition T i . The waiting times are exponentially distributed. A firing rate expresses how often (in average) a transition will fire in a single time unit. A transition with the firing rate r(T i ) = 1 will fire once every time unit, a transition with the firing rate r(T i ) = 0.5 once every second time unit and a transition with r(T i ) = 2 will fire twice in one time unit.
To specify the firing rates, click on the "Firing rates"-button in the "Sto-  The firing rates of the transitions can be dependent on the current marking.
This can be activated in the "Preferences"-window by selecting the "Marking dependent firing rates"-check box. In this case, the firing rate of transition T i is multiplied by the number of possible firings q(T i ). q(T i ) describes how many firings the transition can perform at current marking before the tokens on the pre-places are depleted. The probability density function of the waiting times evolves to To mimic the stochastic behavior, waiting times are drawn using a random number generator (RNG). Each simulation should be unique and the waiting times, and by that the sequence of the fired transitions, slightly different. This randomness is given by the "pseudo random" number sequences generated by the RNG. However, the sequence produced by the RNG is not really random and depends on the "seed", a 48-bit long integer, of the generator. Two simulations of the same parameters started from the same seed will always produce identical results. The seed can be changed by clicking the "Set RNG seed"-button in the "Preferences"-window (figure 11). The current seed is shown in the appearing window. It is important to understand that simulation results can only be identical if their seeds were equal before any simulation step was performed.
The  is computed automatically upon clicking on the "Save"-button.
Finally, rate constants must be defined for all reactions in the second table.
The units are mole per liter and seconds; the concrete unit depends on the order of the reaction. For zero order reactions, the unit is M × s −1 , for the first order s −1 and for the second order M −1 × s −1 . New settings are applied by clicking on the "Save"-button. Closing the window without saving will discard all changes.
In the "Preferences"-window you can set the seed for the random number generator. See 2.3 for more details on the random number generator and seeds.
Furthermore, the user can limit the number of parallel simulation runs. This The exact algorithm should be used for simulating systems with low molecule count, where stochastic effects have more significance.

Using fast simulation mode
To use the fast simulation mode, start the "Gillespie SSA" from the "Simulation"-tab of the ToolBar first. Enter the parameters of the system, like the volume, the initial concentrations of the compounds and the reaction rate constants, by clicking on the "Enter simulation data"-button or by loading an existing setup from an XML file with the "Load setup"-button. You can also set the seed of the Random Number Generator with the "Set RNG seed"-button from the "Preferences"-window. Any other setting, including the number of steps, snapshots or places to plot, will have no effect on the fast simulation mode. Clicking on the "Fast simulation mode"-button will bring the window of the fast simulation mode to the foreground (figure 13).
One important point that needs to be understood is that each new created fast simulation mode starts with the same seed of the RNG. That means that if you start two fast simulation modes from the same Gillespie SSA with the same settings, the outcome of the simulations will be identical. If you want a diversity in the results, you have to use different seeds. However, you must not change the seed each time by going to the preferences. The button "New random" in the fast mode window creates a new RNG with a new random seed.
In the upper text field you can specify the name and the location of the file for the results. The file consists of the elements which are separated by tab stops and line breaks. The first element of the first row is called " Step", the second In the drop-box "Select algorithm:" you can select the Exact SSA or the Approximate SSA. The exact algorithm simulates an occurrence of one reaction per step (described in [3]). This is the precise version of the SSA and provides a very close insight on simulation course, but is very time intensive for the larger systems or the high molecule counts. The approximate algorithm tries to speed up the simulation by guessing the number of reaction occurrences in a defined time step and executing several reactions per one simulation step (described in [4]). In most cases, the approximate algorithm provides a measurable speed increase compared to the exact one without significant precision loss on higher molecule counts. As a rule of thumb, one could say that for "low" molecule counts (up to 1000 for the most species) the exact method should be used. For the higher numbers of reactants the approximate method should be applied.
Some comparisons of the methods can be found in the section 4.3.
The next drop-box allows to set up the Number of parallel simulations in the range from 1 to 500. Parallel simulation runs start with the same input data but different seeds and so allow the examination of the stochastic effects.
The results are written in separate files, for each run an index is appended to the name specified in the "Output"-field. Note: The number of runs which are effectively executed parallel is limited to the number of cores available to the system (usually it would correspond to the number of cores of the machine).
Remaining runs are queued and executed as soon as some of other runs are finished. Runs are considered as "finished" if the time given in "Time span to simulate" is simulated. You can also limit the number of parallel simulation runs (disregarding from which "fast simulation mode"-instance they are executed) in the "Preferences"-window.
In the text field "Time span to simulate" you can specify which time should be simulated. The value consists of four parts separated by a ":", which are (from left to right) days, hours, minutes and seconds. The simulation will run until the specified time is reached and will stop after that. If the field is left with the value "000:00:00:00", the simulation will run as long as there are reactions which can occur, or the user stops the simulation manually (identical to the "Continuous mode").
By default, each simulation step is written into the output file. For long term simulations the number of steps can be very large, and the size of the output file becomes unmanageable. In the text area "Write results into file interval:" you can specify the interval of simulation time between writing the results. For example, if you select "00:00:00:01", the results will be written upon every simulated second. You can even change this value during the simulation (you need to stop the simulation, change the interval value and start the simulation again). If you wish, you can first simulate the system for 5 minutes (type "000:00:05:00" in the "Time span to simulate"-text field) with every step protocoled (leave the "Write results into file interval:"-text field with "00:00:00:00"). Let us assume that critical and rapid reactions happen during this time, and after that, the system is developing towards its steady state very slowly. You can now select the output interval to one second or even one minute and let the simulation run for further hours or days. In this case, you have both, a high resolution of the critical processes in the first 5 minutes and a manageable output file size because of the coarse resolution of the slower steady state development afterwards.
After having set up the simulation time, the update interval and the number of parallel runs, click on the "Start simulation sequence"-button. You can now see up to 500 tabs (depending on the selected number of the simulation runs) with the information about the simulation progress (figure 13), and the button's capture changes to "Stop simulation sequence". Clicking on it will terminate all simulation runs. The information in the big text field is updated every ten seconds and shows the number of the simulated steps and the simulated time. After the defined time has been simulated, the simulation stops and the button turns back to "Start simulation sequence".
Each tab has three control buttons. The button "Apply state to PN" writes the numbers of the compound molecules to the marking of the PN. This is the convenient way to apply simulation results to the PN. The button "Places to plot" opens a window where you can select which places' values should be displayed in the plot. The plot is created with the "Show plot"-button and uses the values from the output file. That means that only results which were written into the file will be displayed.
It is possible to change the parameters of the system while a fast mode simulation is running and to start a new instance of the fast mode. You can even run further simulations of Gillespie SSA. It is, however, not allowed to end the Gillespie SSA and start another mode.

Implementation details
This section provides information about the implementation of some features and may be useful to advanced users.

Structure of the XML file for simulation setups
Simulation setups can be stored in an XML file, which has the following structure: • <SimulationSetup>: root element, required.
• <volume>: Element which stores the volume of the simulated system for the Gillespie SSA mode in a child text node.
• <places>: Sub-elements store information about the places of the PN.   • <detReactionRateConstant>: Rate constant of the reaction which is represented by this place. Used by the Gillespie SSA.
• <expressionText>: Text of the mathematical expression which describes the reaction rate constant, stored in a child text node.
• <variable>: Element of a variable of the mathematical expression. Each variable of a mathematical expression must be described by a <variable>-element! Has attributes name (a string with the name of the variable) and placeId (the ID of the place which represents the variable).
Variables pi and Time need not to be specified.

Random number generator
The Pseudo-Random Number Generator (PRNG) used by simulation modes is described in [7]. The Java-implementation is taken from http://www.javamex.
com/tutorials/random_numbers/numerical_recipes.shtml. By default the System.nanoTime()-value is used as the seed. Setting the new seed is implemented as a creation of a new Random-object with the given seed.

Drawing waiting times in stochastic mode
At each simulation step of the stochastic mode a transition with the lowest waiting time is fired. Waiting times are computed according to the firing rates of the transitions. First, a uniformly distributed random number rand in the range between 0 and 1 is drawn using the RNG. If a transition T i has a reaction rate r(T i ), the waiting time dt(T i ) of that transition is If the firing rates should be marking dependent, the static firing rate is

Gillespie SSA implementation
The Gillespie SSA mode implements the Stochastic Simulation Algorithm (SSA) of coupled chemical systems which was first described by Gillespie [3]. It describes the development of a system in a physically plausible way by simulating the occurrences of the chemical reactions utilizing the knowledge of the mass action kinetic. It is an alternative to the widely used deterministic modeling and simulations, which model chemical reactions by a set of Ordinary Differential Equations (ODE).

Converting input data
In contrast to conventional ODE-based approaches, the SSA operates with the (1) Because only reactions of order zero, one and two can be simulated by SSA, equation 1 can be simplified to:

Performing a simulation step
and the reaction rate of the reaction T i is Secondly, after the rates of all reactions are estimated, their sum is computed.
If this sum is zero (which means that no reaction can occur), the simulation is stopped. Otherwise, two uniformly distributed random variables in the range between 0 and 1, U 1 and U 2 are generated. U 1 is used for calculating the next firing time the equation: time of next firing = − ln(1 − U 1 ) sum of reaction rates .
Thirdly, the next occurring reaction is chosen using the random variable U 2 .
This is done in a while-loop, which iterates over the reactions and sums up their rates. When the sum of the reactions' rates is equal to or greater than the sum of all reactions' rates multiplied by U 2 , the algorithm stops the execution of the loop and choses the last processed reaction. In other words, the reaction T i is chosen to satisfy the following equation: Think of a wheel of fortune where reactions have areas of the size according to their reaction rates.

Exact SSA of the fast mode
The fast mode was designed to provide high performance simulation of millions of steps in a reasonable time. Amongst other methods, this is achieved by a complete abstraction from the PN structure using some optimizations which will be described here. The algorithm itself is the same one as used in Gillespie For each reaction, a list of its educts and a list of its products are stored, too.
Finally, reactions which have at least one educt represented by a constant place are stored in a separate list. This sort of a "dependency graph" allows to limit the number of reactions whose rates must be recomputed at each step.
In the very first step of the simulation, rates of all reactions are computed. In the next steps, only the rates of the reactions whose educt's number was changed in the previous step are recomputed. Also the rates of all reactions which have educts represented by constant places are recomputed, as the number of tokens on a constant place could be changed without being affected by a reaction directly (e.g. time-dependent places).

Approximate SSA of fast mode
The approximate SSA algorithm is based on the so-called τ -leaping method which was originally proposed by Gillespie [4], for the implementation see [5].
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 that the reaction rates do not change significantly during the interval. To avoid negative populations, reactions which can fire maximal 20 times are considered as critical and simulated in an exact way only.
At the beginning of each simulation step critical reactions are determined.
The algorithm then computes two firing times -the one of the non-critical reactions τ 1 and that of the critical reactions τ 2 . τ 1 is computed according to equation 24 of part 3.2 of the section "Stochastic Simulation for Biochemical Systems" in [5], using the set of non-critical reactions. τ 2 is computed like the waiting time in the exact SSA using the set of critical reactions only.
If τ 1 < τ 2 , no critical reaction will occur in the next step.

Insulin receptor model
Insulin is an important hormone which regulates various processes, amongst others the glucose intake by the cells. The answer is mediated by the IR which is located in the cell surface membranes. Impairments in the key components of this system may cause such diseases as the metabolic syndrome or Type 2 Diabetes mellitus.
The PN structure which models the system is depicted in figure 14. The PN consists of 7 places and 11 transitions. The unbound inactive insulin receptor is modeled by the place IR, the free insulin by the place Insulin. We consider that only one cell is simulated and insulin concentration is regulated by external processes, therefore this place is constant. Insulin can bind to the free receptor (transition IR Insulin Binding). The receptor-ligand complex (place IR I) can dissociate (transition IR I Insulin Release), and the insulin receptor molecule becomes available for further bindings. Alternatively, the unbound receptor can be internalized to the cytoplasm (transition IR Inter) and is added to the internal receptor pool (place IR Int).
As soon as insulin is bound to the receptor, the receptor-ligand complex can be phosphorylated and thus activated (transition IR I Phos). The bound phos-

Simulating the model
We choose the "Gillespie SSA"-mode for the simulation of this model. The parameters are adapted from [8]. The volume of the system is set to 1E-9 l.
Initial concentration of the membrane located insulin receptor (place IR) is 9E-of insulin should be applied for 24 hours. The mathematical expression for the place Insulin reads if Time <= 86400 then 1E-9 and for the place IR Total IR + IR_I + IR_Int + IR_I_P_Int + IR_I_P.
The (deterministic) reaction rate constants are shown in figure 15. The rate of IR synthesis is dependent on the concentration of the receptor in the cytosol.
If the concentration decreases below the steady state value of 1E-13 M, the accelerated synthesis rate is used. This is described by the expression if IR_Int+IR_I_P_Int >= 1E-13 then 2.78E-19; 1.67E-18.

Comparing Exact and Approximate SSA
To analyze the performance's differences between the exact and the approximate algorithms, two different simulations were performed with each of the algorithms.
The first simulation was performed on the above described model and parameters. Insulin concentration was set to 1E-7 M, and the system was simulated for

Elements of the GUI
This section lists and explains all elements of the user interface that are provided by the simulation module.

"Simulation"-tab of the ToolBar
• "Simulation mode:"-drop-down menu: Selects the mode by which the system will be simulated. Once a simulation has been started, the mode can not be changed. To change the simulation mode, the running simulation must be stopped. No parallel simulations of the different modes are allowed. Following modes can be selected (for more details, see " Stepby-step guide" 1): • Gillespie SSA (default value): Simulation of the natural behavior of (bio-)chemical systems. Chemical compounds or states are represented by places and reactions by transitions. Each chemical reaction has a reaction rate constant, chemical compounds have initial concentrations. Simulation mimics the development of the system by "performing" reactions, which take place depending on their rate constants and compound concentrations.
• Asynchronous: Simulation mode in which one transition is fired at each step. The transition is randomly selected from the set of active transitions.
• Synchronous: Simulation mode in which multiple active transitions are fired simultaneously at each step. Concurrent transitions are randomly selected.
• Stochastic: Simulation mode in which one transition is fired at each step. The probability of a transition to be selected is given by its firing rate. Firing rates can be dependent on the activation state of the transitions (transitions, which are activated repeatedly, have higher firing rates).
• "Start simulation"-button: Starts the selected simulation mode. Disables the "Control"-tab of the ToolBar. Disables net editing. Enables control elements of the simulation mode. Changes the behavior of the mouse in the NetViewer to the simulation mode (selecting an active transition fires it; selecting a place and right-clicking on it enables the token number input).
• "End simulation"-button: Stops the current simulation and switches the controls back to the editing mode. Disables the simulation controls.
• "Steps to perform:"-text area (default value 1): Specifies the number of steps the simulator will try to perform. It will stop earlier if no active transition exists.
• "Continuous mode"-check box: If selected, simulation steps will be performed until no active transitions left or interrupted otherwise. A firing sequence can be interrupted by pressing the "Stop simulation sequence"or the "End simulation"-button or closing the NetViewer.
• "Start simulation sequence"-button: Begins the execution of a specified number of simulation steps.
• "Statistics"-button: Opens a window with information about the current simulation state, including the number of performed steps, the total number of fired transitions and the numbers of firings per transition.
• "Show plot"-button: Opens a window with an XY-plot of simulation results.
• "Preferences"-button: Opens a window for adjusting various simulation parameters.
• "History of performed steps"-list: List of the 100 last performed simulation steps. Selecting a step will set the system state (marking) corresponding to the step and update the visual output. The step whose state is currently shown is highlighted by red background. Initiating a new firing sequence will erase all following history entries.
• " Step back"-button: If possible (not possible if the currently selected step is the first one), performs a simulation step from the history previous to the selected one.
• " Step forward"-button: If possible (not possible if the currently selected step is the last one), performs a simulation step from the history following the selected one.
• "Snapshots"-list: List of the simulation snapshots. A snapshot stores the last 100 steps of the simulation and is created after every 100 steps by default. It can be deactivated in the "Preferences"-panel. •

"Simulation preferences"-window
• "Create log file"-check box (default disabled): If selected, simulation results are written into a tab stop separated textfile which is stored in the location defined in the text field below. See section 1.9.1 for details on the file structure.
• "Save snapshots"-check box (default enabled): If selected, a snapshot of the system state, containing the last 100 steps, is created every 100th step.
• "Built-in plotting"-check box (default enabled): If selected, a plot of the simulation results can be created with the "Show plot"-button. If disabled, the "Show plot"-button is inactive.
• "Places to plot"-button: Enables selecting places whose results will appear in the plot.
• "Delay (ms)"-text field (default 0): Time in milliseconds between the firing steps. If 0, firing steps are executed with no delay, otherwise the simulator will hold for the defined time before performing the next step.
• "Update interval (firings)"-text field (default 1): The interval between updating the visual output in the NetViewer (i.e. number of tokens on places and transition states), given in simulated steps. If 1, visual output is updated after each simulation step. If 0, the NetViewer is only updated when the simulation sequence stops. Disabling visual update or increasing the interval may improve simulation speed.

NetViewer in simulation mode
• Token numbers are drawn inside places (as red dots for token numbers below 5 or as a digits for token numbers above or equal to 5).
• Black transition: Transition is inactive and cannot fire.
• Black transition with green border: Transition is active and can fire.
• Black transition with red border: Inactive transition which fired in the last step.
• Green transition with red border: Active transition which fired in the last step.
• Left clicking on an active transition fires it.
• Left clicking on a place followed by right clicking: opens a pop-up menu: • "Set the number of tokens for this species": Opens a window where the number of tokens for the place can be entered. Selecting Cancel or closing the window discards the changes.
• "Constant token numbers"-check box: If selected, the place is constant, i.e. the number of tokens on this place will not be modified by the simulation.

Synchronous mode -specific controls
• "Fire at once (% of active transitions)" in the preferences menu The value is rounded up.

Stochastic mode -specific controls
• "Simulated time"-text label: shows the simulated time (in arbitrary units). Updated every time when a firing sequence is completed.
• "Firing rates" -button: Opens a window where the firing rates of the transitions can be specified. A firing rate can be altered by double-clicking on it and confirming a new value using the "Enter"-key.
• "Marking dependent firing rates"-check box in the preferences menu with the activation degree q(T i ).
• "Set RNG seed"-button: Opens a window for setting the seed of the random number generator (RNG) used for drawing the waiting times. Two simulations started from identical marking, with identical firing rates and identical RNG seed will always produce identical results. Precisely, a new RNG is generated each time when the "OK"-button of the seed input window is pressed. The number which is shown in the seed input window by its creation is the currently used seed.

Gillespie SSA -specific controls
• "Simulated time"-text label: shows the simulated time (in seconds).
Updates every time when a firing sequence is completed.
• "Enter simulation data"-button: Opens a window where the user can enter the parameters of the simulated system, including the volume of the system (in liter), the initial concentrations (in M) of the compounds and the rate constants of the reactions (in M × s −1 , s −1 or M −1 × s −1 ). Reaction rate constants and concentrations of the constant places are defined as mathematical expressions. Scientific notation is allowed.
• "Fast simulation mode"-button: Starts the faster implementation of the Stochastic Simulation Algorithm. This mode uses the PN based description of the system but does not update its state during the simulation.
It uses the input data such as volume, initial concentrations of compounds, reaction rate constants and the RNG seed of the Gillespie SSA mode. Up to 500 parallel simulations with different RNG seeds are allowed, and the results are written to separate text files. The fast mode can use the exact stochastic simulation algorithm or the approximate one for simulation speed up without significant accuracy loss.
• "Set RNG seed"-button: Opens a window for setting the seed of the random number generator (RNG) used for drawing the waiting times. Two simulations started from identical marking, with identical firing rates and identical RNG seed will always produce identical results. Precisely, a new RNG is generated each time when the "OK"-button of the seed input window is pressed. The number which is shown in the seed input window by its creation is the currently used seed.
• "Limit number of threads to:"-check box: If selected, only the specified number of simulation runs of the fast simulation mode can be executed in parallel. If more runs are started, they will be queued and executed as resources become available.

Fast simulation mode
• "Output"-text field: Location and name of the CSV file in which simulation results will be written. If several simulation runs are started, the index of the run is appended to the defined name.
• "Start simulation sequence"-button: Starts the defined number of simulation runs and executes them until the specified time is simulated or the simulation is aborted by the user.
• "Stop simulation sequence"-button: Sends a termination request to all simulation runs. When all runs have successfully terminated their work, the button turns back to "Start simulation sequence".
• "Select algorithm:"-drop-box: Allows to select which algorithm should be used for the simulation: • Exact SSA (default): In each simulation step, the time of the next occurrence and the reaction which will occur are computed. Only one reaction per step is simulated.
• Approximate SSA: Estimates how often each reaction will occur in a specific time interval (time intervals selected dynamically). Speeds up the simulation for the large molecule amounts, but can be less precise on low molecule numbers.
• "Time span to simulate:"-text field: Defines the time which should be simulated (days:hours:minutes:seconds). The default value is "000:00:00:00", which means that the simulation will continue until it is interrupted by the user or no reaction can occur.
• "Write results into file interval:"-text field: Interval of simulated time at which the results will be written into the file. The default value is "00:00:00:00", which means that every step is protocoled.
• "Number of parallel simulations"-field: Allows to select up to 500 simulation runs, which will be started with the equal data but different seeds of the RNG.
• "New random"-button: Every fast mode simulation which is started with the "Fast simulation mode"-button has the same seed of the RNG.
This produces identical results if the initial parameters are identical and the seed was not altered. Pressing the "Next random"-button sets a new, random seed for the RNG of the fast simulation instance.
• "Apply state to PN"-button in the run tab: Writes the actual number of compound molecules into the marking of the PN.
• "Places to plot"-button: Opens a window where the user can select which places' values should be plotted.
• "Show plot"-button: Opens a window with an XY-plot of the simulation results of the current run. The values are taken from the output file, so the "Write results into file interval:"-text field has the influence on which points will appear in the plot. A: Generally, scientific simulations consume much resources. However, you can speed up the process by disabling visual components and following some rules: • Try to avoid mathematical expressions where these are not necessary.
• Try to use as few variables in the mathematical expressions as possible. Check that no unused variables are selected in the "Used variables"-list (see 1.10.1).
• Disable updating the visual output (number of tokens on the places in NetViewer) or increase the update interval in the "Preferences"window (see 1.9) (speeds up simulation significantly).
• Disable built-in plotting or reduce the number of plotted places in the "Preferences"-window (can consume a lot of memory for thousands of simulation steps).
Of course, disabling features is annoying and you would like to have great performance with full feature set. Unfortunately, these features consume computational power at various amounts, so you have to decide for yourself which functions to use in order to find the best balance between the functionality and the performance.
8 Changelog • 04.03.2014 • Fixed a bug which appeared when loading a saved marking • Several bugfixes. • Mathematical expressions have a new function for integer division div(x,y).
• Number of parallel simulation runs increased to 500. The actual number of parallel threads is limited to 50, remaining runs are queued and executed upon finishing previous runs.
• Performance of parallel simulation runs improved.
• Confirmation-dialog when trying to close Fast simulation mode with running simulations added.

Acknowledgments
Ina Koch, Jörg Ackermann, Klaus Lindauer for supervising. Jens Einloft for support of MonaLisa and cooperation on implementation of simulator. Jennifer