In this paper, we present an approach to improve the efficiency of stochastic simulation for large and dense networks. We are using stochastic Petri nets as modelling framework. The underlying continuoustime Markov chain (CTMC) is converted to an equivalent discrete-time Markov chain (DTMC), this itself gains no efficiency. We improve the efficiency via discrete-time leaps, even though this results in an approximate method. The discrete-time leaps are done by applying the maximum firing rule, this reduces drastically the number of steps. The presented algorithm is implemented in our modelling and simulation tool Snoopy, as well as in our advanced analysis and model checking tool MARCIE. We demonstrate the approach on models of different sizes and complexities.
Throughout the past decades typical biological models increased in their size and
complexity, because of advances in systems and molecular biology, in particular
through the high-throughput omic technologies. This progress gained momentum
for the arise of multi-scale modelling [
We present an approach to improve the efficiency of stochastic simulation for large and dense networks. We are using stochastic Petri nets as modelling framework, because of the well defined mathematical background and the nice graphical representation. The underlying continuous-time Markov chain (CTMC) is converted to an equivalent discrete-time Markov chain (DTMC), this itself gains no efficiency. We improve the efficiency via discrete-time leaps. The discrete-time leaps are done by applying the maximum firing rule, this reduces drastically the needed number of steps. The algorithm is called discrete-time leap method for the simulation of stochastic Petri nets.
We start by giving the mandatory definitions and continue with the presentation of the discrete-time leap method. Afterwards we demonstrate our approach on some case-studies of different sizes and complexities. 2
A Petri net is a 4-tuple PN = (P, T, f, m0) with a finite set of places P = { p1, p2, . . . , pm} , defining the state of the net and a finite set of transitions T = { t1, t2, . . . , tn} , defining the state changes in the net. The sets of places and transitions are disjoint and not empty. Places and transitions are connected via directed arcs and defined by the arc weight (multiplicity) function f : ((P × T ) ∪ (T × P )) → N. We define any marking m as a mapping m : P → N0. It maps the set of places onto the set of natural numbers, where m(p) defines the number of tokens in place p ∈ P . If appropriate and the context precludes confusion, we treat a place like a variable and write simply p instead of m(p). An initial marking is denoted by m0. The set of pre-places (input places) of transition t is defined as • t = { p ∈ P | f (p, t) > 0} . The set of post-places (output places) of transition t is defined as t• = { p ∈ P | f (t, p) > 0} .
A transition t ∈ T is called enabled in marking m, if ∀ p ∈ • t : m(p) ≥ f (p, t). If a transition t ∈ T is enabled in m, it may fire. The firing of a transition t occurs in two parts. First, the transition removes the amount of tokens f (p, t) from each of its pre-places • t. Second, it adds the amount of tokens f (t, p) to each of its post-places t• . The change in the marking induced by firing transition t is denoted by Δt = { p ∈ • t ∪ t• : Δt (p) = f (t, p) − f (p, t)} . (1) When t in m fires, a new marking m0 = m + Δt is reached. This is denoted by m→− − t m0. The firing itself does not consume any time and takes place atomically. This defines the standard firing rule .
A transition t may still be enabled after firing once, i.e., it is possible to fire transition t in a sequence like this m→− − t m→ 0− − t m→00− − t . .→ .− − t mn until it is not enabled any more. The number of occurrences of transition t in such a sequence starting in marking m is named enabledness (concurrency) degree and is defined as edt = minp∈ • t m(p) f (p, t) .
The standard firing rule given above can be extended by the enabledness degree in the following way m0 = m + Δt · edt(m). (2) Usually referred to as maximum (auto-concurrency) firing rule.
A stochastic Petri net (SPN ) builds on PN , but transitions have an exponentially distributed firing delay, characterised by the firing rate λ . A transition may lose its enabledness while waiting for the delay to expire. The firing itself does not consume time and follows the standard Petri net firing rule. The firing rates are typically transition-specific and marking-dependent. They are defined by stochastic firing rate functions, also known as propensity or hazard functions. The mapping V : T → H, where H is the set of hazard functions, associates to each transition a function ht from H. The sum of all transition rates in a marking m is denoted by E(m) = Pt∈ T ht(m).
The semantics of a SPN is a continuous time Markov chain (CTMC). It is a stochastic process with an exponential probability distribution that has the Markov property. That means, its future behaviour depends only on its current state and not on former behaviour. The denfiition of a CTMC is comparable with the reachability graph, except that the arcs are labelled with the stochastic firing rate of the transition. A CTMC of a SPN is a tuple CTMCSPN = (RSP N (m0), Q, m0) with RSP N (m0) denoting the state space of the underlying net, the transition rate matrix Q = RSP N (m0) × R SP N (m0) → R0+ and m0 the initial state. Now, the probability of a transition t enabled in state m to fire (which results in state m0) within τ time units is 1 − e− Q(m,m0)· τ . (4)
When working with biological systems modelled as SPN , stochastic
simulation is a suitable analysis technique. Although in principle known a long time
before, Gillespie was the first who developed a supporting theory for stochastic
simulation of chemical kinetics [
Starting from the initial marking m0, one has to repeatedly fire transitions. In order to fire a transition, one must answer two questions: 1. When will the next transition fire? 2. Which transition will fire next? The enabled transitions in the net compete in a race condition. The fastest one determines the next marking and the simulation time elapses. In the new marking, the race condition starts anew.
The SSA creates a sequence of discrete random variables X(τ ). The discrete random variable Xp(τ ) describes the number of tokens on place p ∈ P at time τ . The system state (marking) at time τ is thus a discrete n-dimensional random vector X(τ ) = (Xp1 (τ ), . . . , Xpn (τ )) ∈ X . The time Δτ to the next transition is an exponentially distributed random variable with mean 1/E(m); the probability density function (pdf) is
The next transition to fire is a discrete random variable with probability mass function (pmf):
P (Δτ
| m) = E(m) · e− E(m)· Δτ .
P (t | m) = ht(m)/E(m). (5)
Given the system is in state X(τ ), the probability that a transition t ∈ T will occur in the time interval [τ, τ + Δτ ) is given by:
P (Δτ, t | m) = E(m) · e− E(m)· Δτ
· ht(m)/E(m) = ht(m) · e− E(m)· Δτ . 3
The discrete-time leap method (δ -leaping for short) aims at the simulation of stochastic Petri nets used for modelling biochemical reaction networks. In an attempt to decrease the time complexity of the simulation algorithm, we exploit the uniformization of the underlying CTMC in combination with the maximum firing rule.
The idea of converting a CTMC into a DTMC goes back to [
Generating paths through the DTMC is as expensive as for the CTMC and we would not gain any efficiency by doing it in an exact way. That’s why, we are leaping over several states. That means, all enabled transitions that are not mutually exclusive, are forced to fire within one leap. When the net is filled up with tokens, every transition will fire within every leap. Furthermore, we embody the maximum firing rule and let each transition fire concurrently to itself. (7) (8) (9) 3.1
How often a transition is allowed to fire concurrently depends on its enabledness degree and is determined randomly at each step.
firing rate u random[0, enablness degree]
The construction of the DTMC induces that the times between transitions are all exponentially distributed. Hence, these times are randomized by a Poisson process. Since for the uniformized DTMC the number of transitions in any time interval of length δ has a Poisson distribution with rate λ . The Poisson distribution with rate λ is an approximation of the bounded discrete binomial distribution with two parameters k and pr according to the Poisson limit theorem: λ = k · pr.
The binomial distribution is used to model the number of successes in a sequence of k independent yes/no experiments with a probability pr to succeed. In our case, the enabledness degree corresponds to the sequence’s length k = edt. The success probability pr is deduced from equation (4), because of the exponentially distributed times between transitions. Given out of edt maximum firings and a firing rate ht, we compute the probability pr for transition t in marking m for δ units of time as follows pr = (1 − e− ehdtt((mm)) · δ 0 edt(m) > 0 otherwise.
(10)
This leads us to a good approximation of the stochastic simulation results.
We compare the simulation results of δ -leaping, with δ = 0.1, against the direct
method [
P1 P2
P1
P2 P1 P3 P2
P1 P3 P2 2
The second order reaction P 1 + P 2 → P 3 uses mass-action kinetics with rate constant 0.1. The results shown in Fig. 2 are quite as good as in Fig. 1.
We discussed the firing of single and independent transitions, but a typical model consists of much more than that. There are two types of subnets that need some attention: conflicts and sequences.
A conflict exists inside a Petri net, if two or more transitions share a pre-place, see Fig. 3. If the amount of tokens on this place is just as much as the arc weights, one has to determine, which transition will fire. This is not an issue in the stochastic simulation algorithm. In an exact stochastic simulation only one transition is selected at a time and there is no need for further conflict resolution. As in the standard coniflct resolution for Petri nets, we have to choose a transition non-deterministically. The standard conflict resolution proposes a non-deterministic selection according to a uniform distribution on the number of affected transitions. But this is not the best solution for us, because this does not pay attention to the firing rates of the transitions. We are aiming at a weighted non-deterministic selection, which accounts for the firing rates of transitions as well.
P1 P2 P3
P1 P2
P3 100 80 60
The handling of transition sequences, see Fig. 4, is closely related to the conflict resolution. We embody the maximum firing rule and force every transition to fire (if enabled) in one time step. We shuffle the transitions and let them fire (if enabled) sequentially to approximate the stochastic behaviour.
Luckily, we can treat both issues, transitions in conflict or in sequence, in
one solution. We generate a weighted random sequence of all transitions t ∈ T in
each step and let them fire (if enabled) sequentially. The serial firing precludes
additional attempts to solve conflicts. Our algorithm is based on the modern
version of the Fisher–Yates shuffle [
P2 heuristic method, but it provided the best results in our experiments. wt =
X f (p, t)f(p,t) p∈ • t
We get quite meaningful results using our weighted random shuffle algorithm for the simulation of conflicts (Fig. 3) and sequences (Fig. 4). The results are very close to stochastic simulation. The heuristic method, used to compute the transition weights, needs further attention in future work for non mass-action kinetics. 3.3
So far, we discussed the essential steps of the discrete-time leap method. Here, we integrate the various steps in one algorithm, which is given in detail in Algorithm 2. Each simulation run starts with the initialization phase (Algorithm 2, line 2–4), where the simulation time, the marking, and the transition sequence are set. The
Require: transition sequence T , transition weights W , | T | = | W | = n Ensure: shuffled transition sequence 1: procedure weightedRandomShuffle(transitions T , weights W ) 2: for i = n; i > 1; i ← i − 1 do 3: r ← uniformRandomInteger([1, i − 1]) 4: t1 ← Tn− i, t2 ← Tn− i+r 5: w ← Wt2 /(Wt1 + Wt2 ) 6: if bernoulliSampling(w) = 1 then 7: swap(Tn− i, Tn− i+r) 8: end if 9: end for 10: return T 11: end procedure
next step is the generation of the weighted random sequence of transitions using the weightedRandomShuffle algorithm (Algorithm 2, line 6). After that for each transition the following steps are done (Algorithm 2, line 7–14), compute the enabledness degree k, the firing rate h, and pick a random number f according to the binomial distribution defined by k and equation (10), finally the transition fires f -times. The simulation time is elapsed by δ time unit. All these steps are performed until the end of the simulation time interval τ max is reached. Just like for stochastic simulation, one has to do several simulation runs, in order to get reasonable results. The overall result is the mean of all runs. The presented results of the discrete-time leap method approximate the stochastic simulation algorithm quite well. This might not be true for all kinds of biological mechanisms in a stochastic Petri net model, as it is the case for, e.g., a simplified birth-death process, shown in Fig. 5. The result of the stochastic simulation using mass-action kinetics with a rate constant of 1 shows a steady state of P1 at the initial marking. Transitions T1 and T2 are both likely to fire and thus, fire alternately in average. The token consumed by T1 is produced by T2 and vice versa.
Due to the maximum rfiing rule, the result of the discrete-time leap method differs. Let us assume an initial marking of 100 tokens in P1 and each transition consumes 50% of the tokens in each firing for the purpose of demonstration. There exist two possible firing sequences S1 = [T1 , T2 ] and S2 = [T2 , T1 ]. The T1 P1
T2 2 (a) PN 100 80 60 execution of S1 results in the following path m0(10 0→)−− − T1 m1(5 0→)−− − T2 m2(75) and S2 generates m0(10 0→)−− − T2 m1(15 0→)−− − T1 m2(75). Both sequences decrease the amount of tokens on P1. So P1 approaches zero tokens in the long term in any case. This contradicts the exact stochastic simulation results. Equal results can be obtained by adapting the rate constant of the stochastic transition, i.e., setting the rate constant of T2 to 0.5. 4
In this section, we demonstrate our approach on models of the RKIP inhibited
ERK pathway and the mitogen-activated protein kinase, and on two genome scale
metabolic models of E.coli K-12. All Petri nets were modelled with Snoopy [
We compared the discrete-time leap method using δ = 1 with stochastic
simulation using Gillespie’s direct method [
This model shows the influence of the raf kinase inhibitor protein (RKIP) on
the extracellular signal regulated kinase (ERK) signalling pathway. A model of
non-linear ordinary differential equations was originally published in [
The model is scalable by the initial amount of tokens N in the places RKIP,
MEKpp, ERK and RP. All transition rate functions of SPN ERK use mass action
kinetics with rate constants taken from [
We performed experiments with different values of N = { 10, 100, 1000} . The simulation results are given in Fig. 6 and the run-times are provided in Table 1.
Raf1Star
RKIP
Raf1Star_RKIP
ERKPP
MEKPP_ERK
Raf1Star_RKIP_ERKPP
RKIPP_RP
MEKPP
ERK
RKIPP
RP
Raf1Star
RKIP
Raf1Star_RKIP
ERKPP
MEKPP_ERK
Raf1Star_RKIP_ERKPP
RKIPP_RP
MEKPP
ERK
RKIPP
RP
20 40 Time 60 80 100
For all instances of N the results of direct method and δ -leaping match quite
well. For N = 10 the simulation run-time for 100 000 simulation runs for direct
method is lower than for δ -leaping. The simulation run-time increases for both
algorithms with N = 100, but the discrete-time leap method is a little faster now.
For N = 1000 the results follow the ones before, the run-time for direct method
increases by a factor of 10, whereas the run-time for δ -leaping increases by only
30%.
The mitogen-activated protein kinase (MAPK) is the core of the ERK/MAPK
pathway that can, for example, carry cell division and differentiation signals from
the cell membrane to the nucleus. The model was published in [
The Petri net SPN MAP K comprises 22 places and 30 transitions connected
by 90 arcs. The model is scalable by the initial amount of tokens in six places.
All transition rate functions of SPN MAP K use mass action kinetics with rate
constants taken from [
We performed experiments with the following 3 different values of N =
{ 1, 10, 100} . The simulation results for all three instances of N are given in Fig. 7.
They match quite well for the given places, as well as for the others. The run-time
behaviour for this model develops in a comparable way to SP NERK , see Table 1.
The run-time of the direct method in the rfist instance N = 1 is lower than for
δ -leaping, but it increases much faster in the other instances of N .
20 40 Time 60 80 100
This reduced model has been developed to illustrate the basic structure of the
whole genome metabolism of E.coli K-12. The reduction was originally done by
hand [
The Petri net comprises 93 places and 172 transitions connected by 589 arcs. The model is scalable by the initial amount of tokens in 12 places belonging to a place invariant. It has some places with a high connectivity (Fig. 9(a)), e.g., M_h_c(52), M_h2o_c(27) and M_h_e(26). These places are involved in many reactions and changing their values leads to many updates of transition rates. This has a strong influence on the overall speed of the stochastic simulation.
We performed experiments with the following 3 different values of N = { 50, 100, 500} . The simulation results in Fig. 8 and the run-times in Table 1 are quite interesting. The trajectories of M_adp_c, M_amp_c and M_atp_c correlate quite well in the direct method and δ -leaping. The traces of M_accoa_c, M_coa_c and M_succoa_c show some similarities, but M_accoa_c reaches a non zero steady state for N = 500 in direct method, whereas it goes to zero in δ -leaping. The reason for this, might be the different firing rules. The plots of M_q8_c and M_q8h2_c differ the most.
The model is still under investigation for its structural correctness and there
are no kinetic rate constants available; so divergent simulation results have
been expected, as well as the run-time variations differing in several orders of
magnitude, which speak for themselves. Here, the direct method is clearly out of
race.
2000
4000
6000
substrain MG1655, and is one of the 55 GEM models published by Monk et al. [
160 140 120 100 eokn 80 T 60 40 20 00 160 140 120 100 eokn 80 T 60 40 20 00 102 ity v itc e n n co eam101 s h it w s ce a l p f o # 100 genome scale metabolic model (b).
The Petri net comprises 2046 places and 3703 transitions connected by 13001 arcs. The model is scalable by the initial amount of tokens in 101 of 2046 places. The model is supposed to be rather dense. This is confirmed by looking at the place connectivity in Fig. 9(b). The top six places at the connectivity ranking are M_h_c with 1198 arcs, M_h2o_c with 617 arcs, M_atp_c with 402 arcs, M_h_p with 369 arcs, M_pi_c with 339 arcs and M_adp_c with 314 arcs. Places with such high connectivity have a large impact on the performance of stochastic simulation, because the greater connectivity of a place, the more transitions change the number of tokens on this place and the more transition rates have to be evaluated each time this happens.
We performed experiments with the following 3 different values of N = { 50, 100, 500} . The simulation results given in Fig. 10, as well as the run-times nT11116024600000 MMMMM_____ttttsttdhuddemfccbAmeacCapdpPp___c__cccc eok80 40 20 00 2000 4000 Time 6000 8000 10000 in Table 1 are quite interesting. The plots of the randomly chosen places are comparable to some extent, but others differ a lot. This is not surprising, because the model is under investigation and there are no kinetic parameters available. We were not able to finish one stochastic simulation run within 40 days for N = 500. The simulation run-times for δ -leaping are moderate for such a big and dense model. The scaling parameter N has little influence on the simulation run-time of δ -leaping. 5
We presented the discrete-time leap method for the simulation of stochastic Petri nets. It converts the underlying CTMC into a stochastically equivalent DTMC. Generating paths through the DTMC is as expensive as for the CTMC and we would not gain any efficiency by doing it in an exact way. That’s why, we are leaping over several states. The discrete time model and the maximum firing rule in combination with binomial sampling and weighted random shuffling of the transitions make an efficient simulation algorithm, that leads to comparable results. The weighted random shuffle is in need for further improvements to obtain even better approximations of the stochastic simulation.
We can conclude that the discrete-time leap method computes reasonable results and it has a very good run-time performance especially for larger and dense networks. It is less sensitive to higher number of tokens and thus higher transition rates than stochastic simulation algorithms in terms of run-time. This recommends δ -leaping for models, which stochastic simulation is not capable of simulating in reasonable time, e.g., genome scale metabolic models. Furthermore, it is suitable for in silico experiments, where the differences between modified models are of interest, such as knock-out scenarios of certain species or reactions. For particular models it might be necessary to adapt the kinetic rate constants in order to obtain similar results, see Section 3.4. Here, further investigations are needed if we can compute the required adaptations from the net structure.
The discrete-time leap method is implemented in our tools Snoopy [