In this section we briefly describe the implementation of Snoopy’s hybrid simulator by considering the architecture, available algorithms, export and import, and the deployed external libraries.
Architecture
Figure 3 presents the architecture of Snoopy’s hybrid simulator. This architecture consists of three components: the user interface, which comprises the model editor and the simulation dialog; the simulator, which implements the simulation algorithms as well as storing the currently running models and the corresponding result views; and the Snoopy manager, which connects the user interface with the simulation module. Snoopy’s hybrid simulator deploys a simple graphical user interface to permit a rapid configuration of the core simulation procedure. Figure 4 depicts the user interface; in the following, we discuss each of these components.
Model editor
The model editor permits the graphical construction of hybrid models using (coloured) hybrid Petri net notations defined in [11]. Reactions are represented by transitions, while species are denoted by places. More information about hybrid Petri net notations can be found in [11, 19, 29] as well as in Snoopy’s hybrid simulator user manual [21]. In addition to specifying all reactions, the model editor provides other features to configure the model parameters as well as the initial state. The model editor is applied as a pre-step before executing the simulation.
Simulation dialog
The simulation dialog is the user tool to run and manage the simulation. Through it the simulation experiment can be configured and then executed. Moreover, the simulation dialog provides access to the simulation algorithms that are implemented in Snoopy’s hybrid simulator. Once a model has been constructed, a user can access the simulation dialog through Snoopy’s menu bar. The simulation dialog consists of four parts: model configuration, simulator configuration, import/export, and simulation state (compare Fig. 4).
The model configuration section permits to adjust model settings including initial state, reaction rates and other similar parameters. The simulation configuration section deals with specifying the simulation options including the start and end time point of the simulation, the type of the ODE solver, and the hybrid synchronisation method. The import/export section allows users to configure how Snoopy performs any export or import of simulation results. Finally, the simulation state section serves to start and stop the simulation as well as to monitor the simulation state.
The simulation results can be examined using views. Different result views can be defined to explore the simulation output from different perspectives. Each view has its own window to display the results using a dedicated result viewer. Result viewers permit to render the final data using different plotting techniques such as xy-plot. Finally, view curves can be exported into comma separated files (CSV) for further processing.
Simulation model
After a model is constructed, e.g. using the model editor, it can be sent to the simulator for execution. The simulation module takes a copy of the hybrid Petri net model, but ignoring the layout information. Usually, the model is a collection of species, reactions, stoichiometries as well as associated data such as kinetic rates and kinetic rate constants. Nevertheless, this information is mapped in terms of Petri net components. When the simulation model is partitioned into deterministic and stochastic parts, reactions and species are assigned to either of the two regimes. Unlike other implementations (e.g., in [10]), we consider only one version of the species: either discrete or continuous. Later, if a transformation is required (e.g., from number of molecules to concentration, or vice versa), this can be easily done at the position where such a conversion is required (e.g., in the rate equation). Similarly, when considering the simulation output, such a transformation can also be easily applied. However, if a species is manipulated by both a deterministic and a stochastic reaction, it needs to be represented as a continuous species.
Simulator
The simulation algorithms are implemented as a stand-alone, but built-in simulation library. The simulator module reads the model to be simulated and carries out the execution. A number of algorithms, which will be discussed in the next section, are available to execute a hybrid model. Please note, although the core simulator is implemented as a stand-alone library, it can currently not be used as a stand-alone application. However, we are working on this to achieve a stand-alone application (see the “Future improvements” section).
Views
Views are associated with models. Each view is defined over a set of places or transitions of which the dynamic behaviour shall be displayed when the simulation starts. These place/transition sets can be specified by a regular expression. Multiple views can be defined for the same model. A view is also associated with a viewer that displays the selected information. Views can be manipulated or removed after they were initially added to a model.
Snoopy manager
The communication between the user interface and the simulator is done via the Snoopy manager. The Snoopy manager acts as an intermediate agent that sends the GUI command to the simulator and gets the result back to visualise or export them to the chosen file format. As Snoopy is a stand-alone application, the communication between the user interface and the simulation module is done internally and not through a physical communication channel.
Available algorithms
Snoopy’s hybrid simulator encompasses a set of simulation algorithms that together provide a convenient execution of hybrid biological models. The general idea of the hybrid simulation algorithms implemented in Snoopy is as follow. First, the synchronisation module (the hybrid algorithm) prepares the jump equation (see below). Afterwards, the ODE solver numerically integrates the system of ODEs due to the deterministic part until the jump equation is fulfilled. At this point, the synchronisation module switches the control to the stochastic module to select and fire a stochastic reaction. The exact time point of the stochastic event is determined by the jump equation. In what follows, we outline each of these algorithms.
Haseltine and Rawlings algorithm
This is the realisation of the hybrid simulation idea proposed by Haseltine and Rawlings in [9]. According to this method, a system of ODEs is numerically solved until a stochastic event is to occur. The exact occurrence time of the stochastic event is captured through (1).
$$ \int_{t}^{t+\tau}{\sum\limits_{j=0}^{N}{a_{j}^{s}}(\mathbf{x})dt}=-log(r) \,, $$
(1)
where x is the state vector of the model at time t, τ is the firing time of the next slow reaction, r is a random number uniformly distributed in the interval U(0,1), \(a_{j}^{s}(\mathbf {x})\) is the propensity of the j
th slow reaction, and N is the number of stochastic reactions.
In (1) we aim to determine the value of τ. This is achieved by first generating a random number from U(0,1) and then integrating the propensity equations of all slow reactions together with the system of ODEs due to the deterministic part from the current simulation time t until (1) is satisfied. At this point we know that there is a stochastic event which needs to be fired.
Although this method is very accurate, it requires considerable time to switch from stochastic to deterministic simulation [13]. The performance of this method drops rapidly as soon as the number of stochastic events increases. Thus it is suitable only for simple models where the number of potential stochastic events is limited. Moreover, it can produce better results with ODE solvers that do not collect and use history information to advance the numerical integration time.
Accelerated Hybrid Simulation
To overcome the limitation of the Haseltine and Rawlings method, we follow an accelerated approach introduced in [13]. The accelerated algorithm takes advantage of the model structure to boost the overall simulation performance. According to this method, stochastic reactions are classified into two groups: dependent and independent. Dependent reactions affect the system state of the ODE solvers when they occur, while independent reactions have no effects. Therefore, the ODE solver is reinitialised only when a reaction in the dependent group is fired. Thus, the simulation performance becomes better than for the Haseltine and Rawlings method, particularly for bigger models. For instance, in [13] we compared the performance of the Haseltine & Rawlings method and the accelerated approach using three models. We found that there is a notable performance improvement for all three case studies and for certain models; the latter approach is ten times faster than the former one. This save in runtime is mainly due to the reduction of the number of times where the ODE solver is reinitialised. In order to achieve the better performance, the accelerated method approximates the exact capture of the stochastic event occurrence time given by (1) by another equation given in (2).
$$ {\sum\limits_{j=0}^{N}{a_{j}^{s}}(\mathbf{x})}\cdot \Delta \tau =-log(r)\,, $$
(2)
where Δ
τ is the time difference between the occurrence time of the previous event and the current event. Eqs. 1 or (2) has to be satisfied during the integration, when the Haseltine and Rawlings or accelerated method is used, respectively, until the ODE solver stops and returns the control back to the stochastic regime. Please note that although our approach mainly intends to reduce the reinitialisation of ODE solvers employing history information to advance the simulation time (e.g., multi-step ODE solver; see below), it can also be used with single step solvers (e.g., Runge-Kutta) to reduce the frequent recalculation of the step size after each firing of a discrete event.
Improved hybrid rejection-based stochastic simulation
The numerical integration of (1) as well as its approximation in (2) are computationally expensive to be satisfied. Therefore, in [12] a new hybrid simulation method was proposed based on the rejection-based stochastic simulation algorithm (RSSA) introduced in [30] which avoids the calculation of (1) and (2). The RSSA algorithm defines lower and upper bounds of the reaction propensities to minimise the propensity updates. The propensity lower and upper bounds are calculated based on a lower and upper bound of the system state values called fluctuation interval. Propensities are updated only when one or more of the system state entries move completely outside the defined fluctuation interval. The Hybrid Rejection-based Stochastic Simulation Algorithm (HRSSA) exploits this opportunity by switching from the deterministic to the stochastic regime only when the ODE solver reaches the time of a stochastic event or when any of the system state entries is outside the fluctuation interval. In the former case, the discrete regime does not affect the continuous one, while in the latter case the deterministic regime changes the state of the discrete species during the numerical integration. We apply an improved implementation of this method which combines the accelerated and hybrid rejection-based methods. Currently, the improved hybrid rejection stochastic simulation method is tested as the best hybrid algorithm implemented in our tool in terms of performance (compare Table 2). In [12], the performance of the HRSSA algorithm has been compared with state of the art hybrid simulation algorithms using five benchmark models. It turned out that the HRSSA outperforms all competing algorithms.
Dynamic hybrid simulation
The previously discussed simulation approaches are based on static partitioning. Static partitioning adopts a predefined classification of the model reactions into stochastic and deterministic ones. The partitioning itself is usually performed by the user and exploited afterwards by the simulator during the whole simulation process. This approach is constructive for many applications with a clear cut between reactions which have to be simulated stochastically and those which should be simulated deterministically. For instance, in [31] and in many other similar publications that study cell fate, reactions related to the cell nucleus are considered as stochastic, while those happening inside the cytoplasm are considered as deterministic. However, such a clear cut is not always possible to be achieved for all models during the whole simulation period. Reactions can change their state from slow to fast and vice versa during the simulation. For example, in oscillating biological systems, reaction rates also oscillate with respect to time from fast to slow and the other way around. In this case, dynamic partitioning, where reactions are partitioned repeatedly during the simulation, can play a role in speeding up the whole simulation procedure. Our implementation of the dynamic hybrid simulation is based on the improved hybrid rejection method. Using this approach, reactions are repartitioned as soon as any of the state vector entries leaves the fluctuation interval. This will indeed eliminate the need for frequent checks of whether the set of reactions requires repartitioning.
Pure stochastic and pure deterministic simulation
To improve the comfort when simulating biological models with Snoopy’s hybrid simulator, the user has the option to perform a pure stochastic or a pure deterministic simulation of a hybrid model. The direct method [3] is applied to implement the stochastic simulation, while the SUNDIAL CVODE [32] is used to carry out the deterministic simulation. This is a worthwhile feature during the experimentation phase to compare the hybrid results with the pure stochastic and pure deterministic ones. Using this feature, Snoopy’s hybrid simulator ignores any reaction partitioning specified by the user and reads all model reactions as stochastic or deterministic ones, depending on the selected simulation algorithm.
Parallel multi-run simulation
Similar to stochastic simulation, hybrid simulation of biological models might require the execution of multiple runs to calculate average statistics. Snoopy’s hybrid simulator permits the concurrent execution of different runs to take advantage of the existence of multiple cores in the user’s machine.
Simulation of coloured models
A coloured model (such as \({\mathcal {HPN^{C}}}\)) with finite colour sets can be automatically unfolded to an uncoloured model (such as \({\mathcal {HPN}}\)). See [33] for one of the unfolding algorithms deployed in Snoopy. Thus, the simulation of an \({\mathcal {HPN^{C}}}\) model is done on the automatically unfolded \({\mathcal {HPN}}\) model. When the user starts the simulation of an \({\mathcal {HPN^{C}}}\) model, an unfolding dialogue will be triggered, where the user can select an appropriate unfolding method to perform the unfolding. Afterwards, the simulation methods discussed in this section can be used to execute the unfolded model.
To better explain this idea, we consider the \({\mathcal {HPN^{C}}}\) model presented in Fig. 2. First, the discrete subnet, consisting of the two places: closed and open as well as the two transitions: c
h_o
p
e
n and c
h_c
l
o
s
e, is unfolded into three identical subnets (because the colour set chCS consists of three colours). In other modelling scenarios where the unfolded subnets are not identical we can make use of transition guards to imply the required constraints. Similarly, the continuous subnet, consisting of the place Ca and the two transitions C
a_p
u
m
p and diffuse is unfolded into 10,000 identical subsets (because the colour set G
r
i
d2D consists of 100×100 colours). However, because the transition diffuse has a guard expression, only transitions for the colours satisfying this Boolean expression are added. The transition C
a_i
n
f
l
o
w will have only one copy in the unfolded net because the input and output arcs contain a constant expression. Nevertheless, this procedure does not need to be implemented iteratively as we do in this small example. Instead, it can be viewed as a constraint satisfaction problem (CSP) which can be solved by a dedicated CSP solver (e.g. [34]).
Export and import
Snoopy supports the import/export of Petri net models from/to other tools and formats. First, Snoopy imports and exports the (C)ANDL format (cf. [21]) which is a human readable file format used by other software tools (e.g., Marcie [35] and Charlie [36]) which can be employed for a formal analysis of Petri net models (e.g., structure analysis, model checking, etc.). Moreover, Snoopy reads and writes SBML files [37] according to SBML level 2 version 4 by using libSBML [38]. However, we support only a subset of SBML elements that is compatible to our net classes; specifically we do not support any kind of rules or events. Snoopy passes all tests of the SBML Test Suite comprising supported elements. However, the partitioning of hybrid models is lost when exporting to SBML, because SBML has no support for hybrid models yet. The user has to decide whether the exported model has to be treated stochastically or continuously. Furthermore, a coloured model is exported to SBML by first unfolding it into the low level representation and then performing the export.
Implementation language and external libraries
Snoopy’s Hybrid Simulator has been implemented using standard C++. As a component of Snoopy, it adopts wxWidget [39] to implement the graphical user interface under different operating systems. Moreover, the stochastic and deterministic simulation components are implemented in a modular way such that different algorithms can be easily exchanged to execute the stochastic and deterministic regimes. Snoopy’s hybrid simulator adopts internally an external library, SUNDIAL CVODE [32], to solve a system of ODEs. The ODE library provides two main algorithms: one for stiff and one for non-stiff ODEs. Additional ODE solver modules can easily be added in future releases. We also make use of the C++ library Boost [40] to carry out routine tasks such as input parsing and multithreading support.