Simulative Model Checking of Steady State and Time-Unbounded Temporal Operators

Abstract

When working with large stochastic models simulation remains the only possible analysis technique. Therefore, simulative model checking is the way to go. While finite time horizon algorithms are well known for probabilistic linear-time temporal logic, we provide an infinite time horizon procedure as well as steady state computation, based on exact stochastic simulation algorithms. We demonstrate the approach on models of the RKIP inhibited ERK pathway and angiogenetic process.

Keywords

• simulative model checking
• stochastic Petri net
• steady state
• unbounded temporal operator
• probabilistic linear-time temporal logic

-

Stochastic modelling of biochemical reaction networks is getting more and more popular. This also increases the demand for efficient analysis of such models. While small and medium-sized models can be analysed numerically, we focus on large or unbounded models. Therefore we use stochastic simulation to overcome the problem of state space explosion.

We use stochastic Petri nets (SPN ) [ 1 ] as modelling paradigm, which gives us a complete formalised and standardised framework, as well as an intuitive way of modelling concurrent behaviour. A number of biochemical species N involved in the biological model are represented as places p1 . . . pN , and the reactions between them refer to the transitions t1 . . . tM . The kinetics of a reaction is defined as possibly state-dependent rate function ht assigned to the transition. Places and transitions are connected via directed arcs. Each arc contains the stoichiometric value of the associated species. The semantics of such an SPN is defined as continuous-time Markov chain (CTMC).

The dynamic behaviour of stochastic models can be analysed in different ways. We showed in [ 2 ] that numerical analysis is currently efficient up to 1× 109 states. Beyond this limit, stochastic simulation remains the only possible technique. Stochastic simulation may be performed with approximate or exact methods. An approximate method is τ -leaping [ 3 ], which generates an approximate realisation of the stochastic process. Its advantage is the ability to jump over several transitions and thus be more efficient in trace generation than exact methods. But for simulative model checking, we need to know the exact occurrences of each transition. Therefore, only exact simulation algorithms are suitable for the purpose of simulative model checking, like Gillespie’s direct method [ 4 ] or the next reaction method [ 5 ] by Gibson & Bruck.

In this paper, we extend the finite time horizon model checking algorithm of probabilistic linear-time temporal logic to an infinite time horizon and provide an algorithm to compute simulatively steady state formulas. 2

Stochastic Simulation

In biochemical reaction networks (with n molecular species and k reactions), the molecular reactions between the species are random processes, because it is impossible to predict the time at which the next reaction will occur. Stochastic modelling has therefore become an important tool to fully understand the system behaviour of such reaction networks.

The stochasticity can be described in a time-dependent manner by the Chemical Master Equation. In probability theory, this identifies the evolution as a continuous-time Markov chain (CTMC), with the integrated master equation obeying a Chapman-Kolmogorov equation. When working with biological systems, it may be infeasible to set up the CTMC as the state space X ⊆ Nn can be very large or even infinite. The largeness of CTMCs makes simulation an important analysis technique: instead of computing the CTMC directly, simulation aims at imitating the CTMC by generating different paths of the CTMC, i.e., a sequence of discrete random variable Xl(t). The discrete random variable Xl(t) describes the number of molecules of species Sl, l ∈ { 1, . . . , n} present at time t. The system state at time t is thus a discrete n-dimensional random vector X(t) = (X1(t), . . . , Xn(t)) ∈ X . Given the system is in state X(t), the probability that a transition/reaction of type j ∈ { 1, . . . , k} will occur in the infinitesimal time interval [ t + τ, t + τ + dτ ) is given by:

P (t + τ, j | X(t))dτ = aj (X(t)) exp (− a0(X(t))τ ) dτ For each reaction j, the rate is given by the propensity function aj , where aj (x)dτ is the conditional probability that a reaction of type j occurs in the infinitesimal time interval [t, t + dτ ), given state X(t) at time t. The sum of the propensities of all possible transitions in the current state X(t) is given by a0(X). Thus, the different (enabled) transitions in the net compete in a race condition and the fastest one determines next state and the time elapsed. In the new state, the race condition is started anew.

To analyse or understand the behaviour of a biochemical reaction network, many trajectories need to be simulated for a good approximation of the underlying CTMC. Although in principle known a long time before, Gillespie was the first who developed a supporting theory for a stochastic simulation of chemical kinetics [ 4 ]. He presented the Stochastic Simulation Algorithm (SSA; often also called Gillespie’s algorithm), which is a Monte Carlo procedure for numerically generating CTMC. Since Gillespie’s seminal work, several variants and different implementations and optimisations of the SSA have been proposed. Basically, each variant performs the following steps: 1. Initialise time t = t0 and the system’s state X at time t0. 2. Repeat: (a) determine time increment τ ∈ R (b) select next reaction type j depending on the current state X(t) (c) perform state transition imposed by reaction of type j and update state vector X (d) update time t = t + τ .

until simulation time is reached.

The SSA simulates every state transition event, one at a time, and updates the system after each state transition. To determine the time increment τ and to select the next reaction requires to generate random numbers. Different realisations of the CTMC are obtained by different initialisations of the random number generator. Since reliable statements about the system behaviour (variance) can only be made based on many simulations, the usefulness of the simulation approach depends on the simulation time for each individual trajectory. Accelerating simulations is therefore desirable without changing the basic ideas of the algorithm.

Many variants of the SSA aim at reducing the computational cost of selecting the next reaction that will occur. Cao et al. [ 6 ] keep the reactions with larger propensities at the beginning of the list. The position of each reaction in the list is thereby determined after some pre-simulations. McCollum et al. [ 7 ] maintain a loosely sorted order of the reactions as the simulation proceeds. Instead of arranging the reactions in a linear list, Gibson & Bruck [ 5 ] propose to use advanced data structures (trees) to speed up the search for the next reaction that will occur. However, the time to manage the advanced data structures partially compensates the speed-up due to faster search [ 8 ].

Further performance increases of the SSA are obtained if only those propensities are recalculated that actually have changed after a state transition, whereas all others are reused (e.g., [ 8, 5 ]).

An additional speed-up to the SSA is provided by the approximate method τ -leaping [ 3 ], in which time t is advanced by a preselected amount τ and the numbers of firings of the individual transitions during the time interval [ t, t + τ ) are approximated by Poisson random numbers. Thus, instead of (sequentially) tracing every single state transition, several reactions are executed in parallel. With τ -leaping, it is assumed that all propensity functions are approximately constant in [t, t + τ ), which is referred to as the leap condition. To ensure this, it is important to select τ sufficiently small, but also large enough to accelerate simulation. 3

Model Checking

Stochastic simulation produces traces through the state space of the model. For the specification of temporal formulas we define the probabilistic extension of the Linear-time Temporal Logic (LTL) [ 9 ] with numerical constraints [ 10 ], which is called Probabilistic Linear-time Temporal Logic with numerical constraints (PLTLc) [ 11 ]. The grammar of all PLTLc formulas is given in Definition 1. Definition 1. Probabilistic Linear-time Temporal Logic with Constraints ψ := P./ x [φ ] | P =? [φ ]

./ ∈ { <, ≤ , ≥ , >} , x ∈ [ 0, 1 ] φ := X I φ | F I φ | G I φ | φ U I φ | ¬ φ | φ ∧ φ | φ ∨ φ | σ I := [x1, x2] = x ∈ R+| x1 ≤ x ≤ x2 , omit I = [0, ∞ ) σ := ¬ σ | σ ∧ σ | σ ∨ σ | value E value | true | f alse

E ∈ { <, ≤ , ≥ , >, =, 6=} value := value ∼ value | P lace | $V ariable | Int | Real | f unction ∼ ∈ {

+, − , ∗ , /} The probability operator P has two different modes. If it is used with the question mark as P=? [φ ] then it will return the probability P r(φ ) that φ is true. In the second case, P./ x [φ ] returns true, if P r(φ ) ./ x is fulfilled, false otherwise. In simulative model checking we compute a confidence interval ( c.i.); in consequence of that, we have to introduce an additional return value in the second case. For simplicity, we assume the c.i. to have a lower and an upper bound Bl, Bu ∈ R≥ 0, such that the probability P r(φ ), which is not known in our case, is Bl ≤ P r(φ ) ≤ Bu.

 true  f alse

if x ./ [Bl, Bu] ∧ x 6∈ [Bl, Bu] P./ x [φ ] =

if x 6./ [Bl, Bu] ∧ x 6∈ [Bl, Bu]  unknown if x ∈ [Bl, Bu] The operators ¬ , ∧ , ∨ are the standard boolean operators not, and, or. Whereas X , F , G , U denote the temporal operators NEXT, FINALLY, GLOBALLY and UNTIL. The NEXT operator (X I φ ) refers to true in the next state and within the time interval I. The UNTIL operator (φ 1 U I φ 2) indicates that a state where φ 2 holds is eventually reached within the time interval I, while φ 1 continuously holds. The FINALLY operator (F I φ ) means that at some point within the time interval I a state where φ holds is eventually reached. Whereas the GLOBALLY operator (G I φ ) refers to the condition φ continuously holding true within the time interval I. The latter two are syntactic sugar, as they rely on the equivalences F φ ≡ true U φ and G φ ≡ ¬ F ¬ φ .

A trace T fulfils a linear-time temporal logic formula φ if the following | = relations hold:

T | = X φ ⇒ ⇐ T | = F φ ⇒ ⇐ ∃

T | = G φ ⇒ ⇐ ∀ T | = φ 1 U φ 2 ⇒ ⇐ ∃

T | = ¬ φ ⇒ ⇐ T | = φ 1 ∧ φ 2 ⇒ ⇐ T | = φ 1 ∨ φ 2 ⇒ ⇐ T | = v1 E v2 ⇒ ⇐

T (1) | = φ i ∈ N : T (i) | = φ i ∈ N : T (i) | = φ T 6| = φ T | = φ 1 ∧ T | = φ 2 T | = φ 1 ∨ T | = φ 2 evalState(v1, T (i)) E evalState(v2, T (i))

i ∈ N : T (i) | = φ 2 and ∀ j ∈ N ∧ j < i : T (j) | = φ 1 The function evalState(v, T (i)) assigns a numerical value to the expression v by looking up the tokens that each place x ∈ P (v) has in state T (i) of trace T .

In the next sections we present an algorithm to compute steady state formulas, and afterwards an algorithm to compute time-unbounded temporal operators in a simulative manner. Time-bounded algorithms for simulative model checking are well known, i.e., [ 10 ]. 3.1

Steady State Computation

In steady state simulation, the measures of interest are defined as limits, as the length of the simulation goes to infinity. There is no natural event to terminate the simulation, so the length of the simulation is made large enough to get “good” estimates of the quantities of interest. Steady-state simulation generally poses two problems: 1. The existence of a transient phase may cause the estimate to be biased. 2. The simulation runs are long, and usually one cannot afford to carry out many independent simulations.

These are several methods that allow to cope with these problems to some extent. Among them are: the batch means method, the method of independent replicas, and the regeneration method. Each of these methods has its advantages and disadvantages. In our implementation we use a sample batch means algorithm to compute the steady state.

We choose Skart [ 12 ], which is an automated sequential procedure for on-thefly steady state simulation output analysis, because it is specifically designed to handle observation-based statistics and usually requires a smaller initial sample size compared with other well-known simulation analysis procedures [ 12 ]. This algorithm partitions a long simulation run into batches, computes an average statistics for each batch and constructs an interval estimate using the batch means. Based on this interval estimate Skart decides whether a steady state is reached or more samples were needed. A detailed description of the algorithm is given in [ 12 ].

We extend PLTLc with the steady state operator S. Definition 2 states the syntax of it. The return values are quite the same as for the probability operator P. Inside of S are only state formulas allowed, i.e., no temporal operators. Definition 2. Extension of PLTLc with steady state operator S. ψ := P./ x [φ ] | P =? [φ ] | S ./ x [σ ] | S =? [σ ]

./ ∈ { <, ≤ , ≥ , >} , x ∈ [ 0, 1 ] Steady-state formulas are computed with Algorithm 1. At first the simulation trace is created until the steady state is reached (line 4–8). To get an unbiased result, we cut off the first n states, which bias the steady state (line 9). The steady state probability is now the ratio To/Ts of the occupation time To of the fulfilling states and the simulation time Ts (line 11–19). But this gives correct results only for those Petri nets, where the reachability graph consists of only one strongly connected component. The complexity of this decision is the same as for constructing the reachability graph. To solve this problem, one has to make several steady state computations and average the values. In that way the individual steady state estimates are weighted according to the strongly connected components.

Algorithm 1 Steady state computation for one simulation run . state si, sojourn time ti in si 3.2

Verification of Time-Unbounded Until

In Algorithm 2 we present the algorithm for checking time-unbounded until formulas. It is nearly the same as for time-bounded formulas except that one has to stop the simulation trace at some time point. The decision of doing that is the crucial part of the algorithm. We assume that reaching the steady state is a reasonable stopping criteria (line 38). A system in a steady state has numerous properties that are unchanging in time. This implies that for any property p of the system, the partial derivative with respect to time is zero: ∂p

= 0 ∂t If a system is in steady state, then the recently observed behaviour of the system will continue into the future. In stochastic systems, the probabilities that various states will be repeated will remain constant. 4

MARCIE: An Implementation

MARCIE [ 13 ] is a tool for analysing generalised stochastic Petri nets (GSPN ), supporting qualitative and quantitative analyses including model checking capabilities. Particular features are symbolic state space analysis including efficient saturation-based state space generation, evaluation of standard Petri net properties as well as Computational Tree Logic model checking. Further it offers symbolic Continuous Stochastic Logic model checking and permits to compute expectations for rewards which can be added to the core GSPN . Most of MARCIE’s features are realised on top of an Interval Decision Diagram (IDD) implementation [ 14 ]. IDDs are used to efficiently encode interval logic functions representing marking sets of bounded Petri nets. Thus, MARCIE falls into the category of symbolic analysis tools. However, it additionally comprises approximative and simulative engines, which work explicitly, to support also stochastic analysis of very large and unbounded nets. It includes two exact simulation algorithms, firstly Gillespie’s direct method [ 4 ], and secondly the next reaction method by Gibson & Bruck [ 5 ].

Parallelising the simulation is a good way to speed-up the computation. We use the Message Passing Interface (MPI) to develop a portable and scalable simulation tool for large-scale models. The desired number of simulations runs is equally distributed to the worker processes. Therefore the run-time decreases with the number of workers.

MARCIE is completely written in C++, and makes use of the libraries GMP, pthreads, flex/bison and boost. It comprises about 45,000 lines of source code. MARCIE is available for non-commercial use; we provide statically linked binaries for Linux and Mac OS X. The tool, the manual and a benchmark suite can be found on our website http://www-dssz.informatik.tu-cottbus.de/marcie.html. MARCIE itself comes with a textual user interface. Options and input files can also be specified by a generic Graphical User Interface (GUI), written in Java, which can be easily configured by means of a XML description. The GUI is part of our Petri net analyser Charlie [ 15 ]. Algorithm 2 Unbounded Until for one simulation run 5

Case Studies

In this section we demonstrate our approach on the models of the RKIP inhibited ERK pathway and angiogenetic process. All Petri nets were modelled with Snoopy [ 16 ] and analysed with MARCIE [ 13 ]. The experiments were carried out on a machine with 4x AMD Opteron™ 6276 with 2.3 GHz and 256GB RAM running CentOS 6. 5.1

RKIP inhibited ERK pathway

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 [ 17 ]. Later on, it was discussed as qualitative and continuous Petri nets in [ 18 ], and as three related Petri net models in [ 19 ]. The stochastic Petri net SPN ERK comprises 11 places and 11 transitions connected by 34 arcs and is shown in Fig. 1. All Raf1Star

RKIP ERKpp r1

r2

Raf1Star_RKIP r8 r3 r4 r11 r10 RP MEKpp_ERK

Raf1Star_RKIP_ERKpp

RKIPp_RP r6 r7 r5 r9 MEKpp

ERK

RKIPp Raf1Star + RKIP r1↔,r2 Raf1Star RKIP

r3,r4 Raf1Star RKIP ERKpp Raf1Star RKIP + ERKpp ↔

Raf1Star RKIP ERKpp →r5 Raf1Star + ERK + RKIPp

ERKR+KIMPEpK+ppRPr6↔r,r9↔7,rM10ERKKpIpPEpRRKP →rr→811ERRKKIpPp++RMPEKpp transition rate functions use mass action kinetics with the original parameter values from [ 17 ]. The model is scalable by the initial amount of tokens in the places RKIP, MEKpp, ERK and RP. The more initial tokens on each of these places, the bigger the state space of the Petri net. Table 1 shows the number of reachable states for different initial markings.

In any case such a state was reached, therefore the probability of the formula is 1, see Table 2.

This is the same property as in [ 2 ], so we can verify the correctness of the results. The results in Table 3 show first that the resulting confidence interval covers the probability computed by the Jacobi method in [ 2 ]. Second the algorithm scales nearly linear with the number of worker processes. A very interesting behaviour regards the relationship between the state space size and the total run-time of the computation. One could expect an increase of the run-time, but it stays the same. This is a result of the level semantics described in [ 20 ]. Angiogenesis is a complex phenomenon that goes from a molecular level to macroscopic events. This Petri net models a part of the signal transduction pathway involved in the angiogenetic process and was originally published in [ 21 ]. The stochastic Petri net SPN ANG comprises 39 places and 64 transitions connected by 185 arcs.

The model is scalable by the initial amount of tokens in the places Akt, DAG, Gab1, KdStar, Pip2, P3k, Pg and Pten. The more initial tokens on each of these places, the bigger the state space of the Petri net. The number of reachable states for different initial markings are shown in Table 4. As in the previous case study we check for reachability first. Now we want to know the probability of eventually reaching a state where no tokens reside on place Akt:

P=? [F [Akt = 0]] .

In contrast to SPN ERK , Table 5 shows that the probability ranges from about 0.44 (N = 1) to 0.9 (N = 6). That means a state where no tokens lay on place Akt is not always reached, because the CTMC consists of several strongly connected components and in some of them such a state does not exist. Secondly we compute the steady state probability of being in a state that has no tokens on place Akt:

S=? [Akt = 0] .

The results in Table 6 show that the steady state probability is nearly the same as in the reachability case as the overall steady state probability consists of two parts, first the probability of reaching a strongly connected component and second the steady state probability inside these component. The result means the steady state probability inside a strongly connected component, where a state exists with Akt = 0, is almost 1. That’s why the overall steady state probability almost coincides with the reachability probability. To compute the transient probability of the formula P=?[φ 1 U φ 2] in state s means to compute the probability distribution starting in s and making states absorbing, which satisfy ¬ φ 1 ∨ φ 2. The resulting linear system of equations can be solved numerically by iterative methods like Gauss-Seidel or Jacobi. There are several tools available that support such solvers, among them MARCIE [ 13 ]. The drawback of numerical solvers is their restriction to bounded CTMCs and the complexity is typically O (n) and in worst case O n2 . On the other hand they compute an “exact” result. The same methods were used to compute the steady state distribution of bounded CTMCs for computing the steady state probability of formulas like S=?[φ ].

Statistical model checking [ 22–24 ] is a quite similar approach to simulative model checking, but differs in some details. Hypothesis testing, i.e., sequential probability ratio test (SPRT), has good performance compared to the computation of point estimates, but it can only check formulas like P./ x . In the end, the user gets a result of true or false and has no idea of the scale of the estimated probability.

The Monte Carlo Model Checker MC2 [ 11 ] computes a point estimate of a Probabilistic LTL logic (with numerical constraints) formula to hold for a model. MC2 does not include any simulation engine but works offline by taking a set of sampled trajectories generated by any simulation or ODE solver software.

Last not least, a combination of simulation and reachability analysis were used to compute time-unbounded formulas in [ 22, 25 ]. But this approach suffers from the same restrictions of bounded state spaces as the numerical methods. 7

Conclusions

In this paper we presented an infinite time horizon model checking algorithm plus steady state operator for probabilistic linear-time temporal logic. We verified the results of the simulative approach against the numerical solutions of the Jacobi and Gauss-Seidel methods. We proved the efficiency of our algorithm and the scalability by using several worker processes through MPI.

As our algorithm is based on stochastic simulation, its run-time does not directly depend on the size of the state space, as for the numerical methods, but on the rate functions of the transitions and the structural size of the Petri net. That is the greater the sum of the transitions rates, the smaller the time steps are, and the more simulation steps need to be done to reach a certain time point.

The main drawback of simulation-based methods remains. The achieved accuracy depends on the number of simulation runs and grows exponentially with the expected accuracy. Therefore methods of choice for bounded and mediumsized models are numerical, otherwise simulation plays out its strength.

References

1. Marsan , M.A. , Balbo , G. , Conte , G. , Donatelli , S. , Franceschinis , G.: Modelling with Generalized Stochastic Petri Nets . Wiley Series in Parallel Computing, John Wiley and Sons ( 1995 ) 2nd Edition .

2. Heiner , M. , Rohr , C. , Schwarick , M. , Streif , S.: A comparative study of stochastic analysis techniques . In: Proc. CMSB 2010 , ACM ( 2010 ) 96 - 106

3. Gillespie , D.T.: Approximate accelerated stochastic simulation of chemically reacting systems . The Journal of Chemical Physics 115 ( 4 ) ( 2001 ) 1716 - 1733

4. Gillespie , D. : Exact stochastic simulation of coupled chemical reactions . The Journal of Physical Chemistry 81 ( 25 ) ( 1977 ) 2340 - 2361

5. Gibson , M.A. , Bruck , J.: Efficient exact stochastic simulation of chemical systems with many species and many channels . The Journal of Physical Chemistry A 104 ( 2000 ) 1876 - 1889

6. Cao , Y. , Gillespie , D.T. , Petzold , L.R. : Adaptive explicit-implicit tau-leaping method with automatic tau selection . J Chem Phys 126 ( 22 ) ( Jun 2007 ) 224101

7. McCollum , J.M. , Peterson , G.D. , Cox , C.D. , Simpson , M.L. , Samatova , N.F. : The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior . Comput. Biol. Chem . 30 ( 1 ) ( 2006 ) 39 - 49

8. Cao , Y. , Li , H. , Petzold , L. : Efficient formulation of the stochastic simulation algorithm for chemically reacting systems . J Chem Phys 121 ( 9 ) ( Sep 2004 ) 4059 - 4067

9. Pnueli , A. : The temporal logic of programs . In: Proceedings of the 18th IEEE Symposium on the Foundations of Computer Science , IEEE Computer Society Press ( 1977 ) 46 - 57

10. Fages , F. , Rizk , A. : On the analysis of numerical data time series in temporal logic . In: Proc. CMSB 2007 , LNCS/LNBI 4695, Springer ( 2007 ) 48 - 63

11. Donaldson , R. , Gilbert , D. : A Monte Carlo model checker for probabilistic LTL with numerical constraints . Technical report , University of Glasgow, Dep. of CS ( 2008 )

12. Tafazzoli , A. , Wilson, J.R. , Lada , E.K. , Steiger , N.M.: Skart: A skewness- and autoregression-adjusted batch-means procedure for simulation analysis . In: Winter Simulation Conference . ( 2008 ) 387 - 395

13. Schwarick , M. , Rohr , C. , Heiner , M. : Marcie - model checking and reachability analysis done efficiently . In: Proc. 8th International Conference on Quantitative Evaluation of SysTems (QEST 2011 ), IEEE CS Press ( September 2011 ) 91 - 100

14. Tovchigrechko , A. : Model Checking Using Interval Decision Diagrams . PhD thesis , BTU Cottbus, Dep. of CS ( 2008 )

15. Franzke , A. : A concept for redesigning Charlie . Technical report, BTU Cottbus, Dep. of CS ( 2008 )

16. Rohr , C. , Marwan , W. , Heiner , M. : Snoopy-a unifying Petri net framework to investigate biomolecular networks . Bioinformatics 26 ( 7 ) ( 2010 ) 974 - 975

17. Cho , K.H. , Shin , S.Y. , Kim , H.W. , Wolkenhauer , O. , McFerran , B. , Kolch , W. : Mathematical modeling of the influence of RKIP on the ERK signaling pathway . In: Proc. CMSB 2003 , LNCS 2602, Springer ( 2003 ) 127 - 141

18. Gilbert , D. , Heiner , M. : From Petri nets to differential equations - an integrative approach for biochemical network analysis . In: Proc. ICATPN 2006 , LNCS 4024, Springer ( 2006 ) 181 - 200

19. Heiner , M. , Donaldson , R. , Gilbert , D. In: Petri Nets for Systems Biology , in Iyengar, M.S. (ed.), Symbolic Systems Biology: Theory and Methods . Jones and Bartlett Publishers, Inc. ( 2010 )

20. Calder , M. , Duguid , A. , Gilmore , S. , Hillston , J.: Stronger computational modelling of signalling pathways using both continuous and discrete-state methods . In: Proc. CMSB 2006 , LNBI 4210, Springer ( 2006 ) 63 - 78

21. Napione , L. , Manini , D. , Cordero , F. , Horvath , A. , Picco , A. , Pierro , M.D. , Pavan , S. , Sereno , M. , Veglio , A. , Bussolino , F. , Balbo , G.: On the Use of Stochastic Petri Nets in the Analysis of Signal Transduction Pathways for Angiogenesis Process . In: Proc. CMSB 2009 , LNCS/LNBI 5688, Springer ( 2009 ) 281 - 295

22. Younes , H. , Clarke , E. , Zuliani , P. : Statistical verification of probabilistic properties with unbounded until . In: Formal Methods: Foundations and Applications . Volume 6527 of Lecture Notes in Computer Science. Springer ( 2011 ) 144 - 160

23. Basu , S. , Ghosh , A.P. , He , R. : Approximate model checking of PCTL involving unbounded path properties . In: ICFEM. ( 2009 ) 326 - 346

24. Ballarini , P. , Forlin , M. , Mazza , T. , Prandi , D. : Efficient parallel statistical model checking of biochemical networks . In: Proc. PDMC . ( 2009 ) 47 - 61

25. Zapreev , I.S.: Model checking Markov chains : techniques and tools . PhD thesis , University of Twente, Enschede (March 2008 )