Biochemical reaction systems are usually modeled by ordinary di erential equations (ODEs). For further analysis, they are often transformed into stochastic Petri nets (SPN), whose state space (or reachability graph) then can be studied to deduce properties. If a biochemical reaction system is in a steady state, from now on called steady situation1, then the rates of the reactions and the concentrations of the species are constant. These concentrations and rates can be established by simulation of the SPN-model. A steady situation1, signi es also that on the model level only a subset of all possible reachable states is pertinent to this situation. It would be of interest to isolate formally and constructively this subset of states. To our knowledge there is no way to achieve this using the SPN-model or the ODE-model. In this article we propose an approach to calculate the part of the state space corresponding to a steady situation1. To do so, we map the SPNmodel onto a Time Petri Net-model (TPN) with the same behaviour as that in the steady situation1 observed in the SPN simulation. Using reduction methods for TPNs we can extract the part of the reachability graph of the SPN-model which is relevant for the steady-situation1. We show that this is exactly the reduced reachability graph of the constructed TPN-model. Finally, the later one can be analyzed qualitatively and quantitatively. In addition, this approach helps for validating the correctness of the calculated (and used) rates in the steady situation1and of the parameters used in the original ODEs, xed by experiments in the wet labs, both being -a priori- subject to a certain degree of uncertainty.
When considering biochemical reaction systems we emphasize the interactions between di erent species during time and we do not take a momentary snapshot of the system. This means that time is an indispensable component in each 1 to avoid confusion, between state in the biochemical system and states or state space in the models, the biological steady state will be called throughout the whole paper steady situation model of such a system. Furthermore, the repetitive occurrence of reactions in the system during a certain time, expressed by their reaction rates, de nes the behaviour of the system. It is obvious, that the rates depend on the concentrations of the species involved in the reactions: the higher the concentration the higher the reaction rate. This is the case until the concentrations of the involved species achieve certain levels. Then the concentrations of the species do no longer change, i.e., the reaction rates stay constant. This situation is the so called steady situation1 in an biochemical reaction system. Finally, the occurrence (or taking place) of biochemical reactions is a stochastical one.
The taking place of a particular reaction can be modeled by an ordinary
di erential equation (ODE), the causal relationship between the interactions
is often modeled by some graph. Both aspects can be represented in a unique
model, by using some variant of Petri Nets, such as e.g. Hybrid Functional Petri
Nets [
Please note that an essential point while constructing an SPN-model is the de nition of its initial marking. It should faithfully map the initial concentrations of all species involved in the reaction system.
Considering the models quoted above, it is not possible by formally analyzing it, to extract the steady situation1 in which the reaction system may stay after some time. The only thing which can be done, and which is done in general, is to simulate the model (over a very high number of runs) until being able to deduce properties concerning the steady situation1, with some remaining uncertainty.
In this paper we are using the uncertain data obtained by simulation of the SPN, and also uncertain parameters estimated by measures and interpolation, and prove analytically if they present in fact those of a steady situation1. For this reason, we rst observe the SPN-model during a high number of runs which yield mean values corresponding to the concentration of species and reaction rates in the steady situation1. In function of these values, we map the SPN onto a Time Petri Net-model (TPN), having the same skeleton, such that the behaviour (that of the observed steady situation1) stays the same.
This works, because TPNs have the same semantics as SPNs with constant
rates, and we dispose of a well established theory [11{13] for studying analytically
their behaviour. In particular, our reduction results concerning state spaces of
TPNs [
This paper is organized as follows: in the next section we recall some basic notions and notations of the used Petri Net classes together with the reduction results of TPN state spaces. In section 3 we introduce the mapping from an SPN in the steady situation1 onto a TPN model. In the subsequent section we illustrate our approach on the core model of the in uence of the Raf-1 Kinase Inhibitor Protein (RKIP) on the Extracellular signal Regulated Kinase (ERK) signalling pathway, chosen as running example, before concluding. 2
In this section we recall the concept of TPNs. After that we introduce some basic notions and fundamental properties, which are important for their quantitative and qualitative evaluation. 2.1
Time Petri Nets (TPN) [
As usual, in this paper, N denotes the set of natural numbers, and Q0+, resp. R0+, the sets of nonnegative rational numbers, resp. real numbers. T denotes the set of all nite words over the alphabet T , l(w) is the length of a given word w.
Some 5-tuple Z = (P; T; v; mo; I) is called a Time Petri net (TPN), if S(Z) := (P; T; v; mo), the skeleton of Z, is a Petri net where P; T are nite sets with P \ T = ;, v : (P T ) [ (T P ) ! N de nes the arcs with their weight, mo : P ! N xes the initial marking, and I : T ! Q0+ (Q0+ [ f1g) is its interval function where 8t 2 T , I(t) = [I1(t); I2(t)] and I1(t) I2(t), spezifying the earliest and latest ring time of t: ef t(t) = I1(t), lf t(t) = I2(t).
As shown in [
A marking m : P ! N can be seen as a vector of size jP j, we refer to it as p-marking. Thus each transition t 2 T induces the p-markings t , t+ and t de ned by t (p) := v(p; t); t+(p) := v(t; p) and t(p) := t+(p) t (p). With these notions the ring rule for TPNs can be de ned. A transition t 2 T is enabled at a marking m i t m (e.g. t (p) m(p) for every place p 2 P ).
The pre-sets and post-sets of a place or transition x are given by x := fy j v(y; x) > 0g and x := fy j v(x; y) > 0g, respectively.
An example for an arbitrary TPN is shown in Fig. 1.
Every possible situation in a given TPN can be described completely by a state z = (m; h), consisting of a p-marking m (the standard marking) and a transition-marking (short: t-marking) h. The t-marking is a transition vector, which describes the current time circumstances in a certain situation. More exactly, each component of the t-marking is either a real number or the sign ]. Thus h(t) can be seen as clock of t. If t is enabled at a marking m, its clock h(t) shows the time elapsed since t became most recently enabled. If t is disabled at m, the clock is switched o (indicated by h(t) = #).
Formally, a pair z = (m; h) with m : P ! N and h : T ! R0+[f#g is called a state of a TPN Z = (P; T; v; mo; I) if 8t 2 T , either (t m and h(t) lf t(t)) or (t 6 m and h(t) = # ): The initial state zo := (mo; ho) of the TPN Z is given by de ning ho as follows 8t 2 T; ho(t) := 0# iiff tt 6 mm00:
Thus, the initial state of Z1, as given in Fig. 1, is z0 = ( (0; 1; 1) ; (0; ]; ]; 0) ). | {z } | {z } p-marking t-marking
The state z = (m; h) is called an integer state, if h(t) is an integer for each enabled transition t in m.
The behaviour of a TPN is de ned by changing from one state into another by ring a transition (without auto-concurrency) or by time elapsing. In a to, denoted by z new clock satis es 8t 2 T; h0(t) := t!o z0 , where the new marking is m0 = m + <8 h#(t) iiff tt 6=6 t0m;0(t + to ) : 0 otherwise.
ring of such
to and the m; t m0
ef t(to) h(to).
Then the state z change into a state z0 = (m0; h0) by the
This de nition implies, that in the case that t0 is still enabled after the ring of t0, it can only re re after at least t new waiting units. To resume, concurrency but no auto-concurrency is possible by the way the evolving of the clocks is de ned.
The state z may also change into a state z0 = (m; h0) by the time elapsing 2 R0+, denoted by z ! z0; where the marking stays the same, but time goes on : 8t 2 T with h(t) 6= # we need h(t) + lf t(t) ) i.e. the time elapsing need to be possible, and the new clock is given 8t 2 T by
h0(t) := h#(t) + iiff tt 6 mm00:
A state z = (m; h) of a TPN Z is called reachable in Z (starting at z0), if there exist states z1; z10; :::; zn; zn0, transitions t1; :::; tn, and times i 2 R0+, for i n, such that z0 !0 z1 t!1 z10 !1 z2 t!2 z20 !2 : : : zn t!n zn0 !n z holds.
The sequence of transitions = t1 : : : tn leading to a reachable state will be called a feasible one (starting at z0) or just a ring sequence of Z. The full sequence ( ) = 0t1 1 : : : tn n is called a (feasible) run of . It shows that in a given TPN the state changes generally consist of alternating series of time elapsing and transition ring. Obviously, for a given run the transition sequence is well de ned, and for a given ring sequence there are in nitely many runs in general.
Eventually, RSZ (z0) is the set of all reachable states in Z starting from an arbitrary state z0. And RSZ := RSZ (z0), that from the initial state is also called the state space of Z.
We may also consider the set of reachable p-markings, also called the p marking space RZ := f m j (m; h) 2 RSZ g in a TPN Z. This is a subset (not necessarily proper) of the reachable markings of the skeleton S(Z). Therefore, a ring sequence in the skeleton S(Z) is not necessarily a ring sequence in Z. The set of p-markings, reachable in Z starting at an arbitrary p-marking m0, is denoted by RZ (m0). A TPN is called bounded, if its set of reachable p-markings is nite, otherwise it is called unbounded.
For di erent reasons the state space of a TPN is in general in nite and dense in terms of the time: the set of reachable p-markings can be in nite or the set of t-markings for a xed p-marking can be in nite or both together. Later on, we consider some approaches for concise state space representations, when RSZ is in nite while RZ is nite.
The de nition of state change by time elapsing can be slightly and
consistently modi ed for the introduction of a reachability graph based on all
reachable essential states for arbitrary TPNs, especially for TPNs including
transitions whose lft s are 1. The set of all reachable essential states for arbitrary
TPNs is de ned as a subset of all reachable integer states of the considered TPN.
We will use the following property (for more details see [
We will formalize the time-dependent notions of measuring the length of runs,
cf. [
Let ( ) be a run of the transition sequence in some TPN Z. The length of the run l ( ) is the sum of all times while executing the run ( ), i.e., n l ( ) := X i=0 i; where n = l( ) and = 0 1 : : : n:
For a given transition sequence in Z, a feasible run ( ) with minimal length will be referred to as minimal run of . Evidently, it satis es : l ( ) := mi0nf l ( 0) j ( 0) is a feasible run of in Z g:
The notion of maximal run can be introduced analogously. It denotes the run with maximal length within all feasible runs of if such an upper bound exists; otherwise it is not de ned.
two states can be found in [
The central problem for the dynamic analysis of a given TPN is the adequate knowledge of its state space. It is important to get a nite description of the in nite state spaces, under the condition that the p-marking space is nite.
It can be shown that - despite the continuous nature of the time intervals - it is su cient to pick up just some \essential" states to determine the entire timed behaviour of the net so that qualitative and quantitative analyses remain possible.
Fig. 2. For some TPN Z. Left hand. Sketched state space of Z: a continuous set.
Right hand. Sketched reduced state space of Z: all reachable integer states.
While the calculation of a single reachable integer state is rather
straightforward, the proof that the knowledge of the integer states is su cient for analyzing
a TPN was quite di cult. Three solutions had been proposed in the past:
considering a global clock [
An edge (z1; z2) labeled (k; t), with k 2 N, in this reduced reachability graph
has the meaning that in state z1 k time units are elapsing before transition t res
leading to state z2. Now, for nding minimal and maximal time paths between
two states/p-markings in a TPN, its reduced reachability graph can be used,
even e ectively. Our algorithms for computing the reduced reachability graph of
a given TPN are implemented in several standard Petri Net tools, like INA [
Stochastic Petri Nets (SPNs) had been introduced at the beginning of the
Eighties, cf. [
Such SPNs are derived from classical PNs by assigning to each transition t a ring rate t. This ring rate speci es a ring delay for the transition. More exactly, the ring delay is a random variable which is distributed exponentially and has t as parameter of the probability density function. In fact, to each transition t a probability density function with parameter t is associated : ft(x; t) =
0; x < 0: ( te tx; x
0; Finally, the ring rate t may be marking-dependent in general. In such a case, we should write t(m), where m is a marking, instead of t. Than the expected value for the ring delay for the transition t in the marking m is t(1m) .
The ring mode is de ned as follows: In a given marking m, each enabled
transition t obtains an instance of the ring delay t(m) from its associated
probability density function. Then a choice is made: the transition with the
minimum ring delay is ring. The ring itself of a transition does not consume
time. The successor marking is than obtained as in the underlying classical PN.
It is well known [
Biochemical reaction networks are mostly described by ordinary di erential
equations (ODEs) or reaction rate equations (RREs), and both can be
converted into each other. Taking account of the rate equations of all reactions in
the systems, ODE like RRE models can be transformed into Continuous Petri
nets or Stochastic Petri nets. More about these transformations can be found,
e.g. in [
A transformation of an ODE model into a Time Petri net model (TPN) using
the reaction rates is shown in [
Independently from the original model, an RRE one or an ODE one, in a rst
step, a timeless Petri net is always derived. This describes the causal relations
between the events in the system. In biochemical systems, these are biochemical
reactions or biochemical signal transductions. Thereafter additional information,
in particular the time parameters, need to be assigned to the Petri net. They
are obtained from the parameters (kinetic rate constants) in the ODEs. Their
values are often determined experimentally. When it is not possible to collect
or identify them in vitro, the parameters are estimated using experimental data
achieved only for some discrete time points. In this case the goal is to estimate the
value of the parameters for each moment so that the values over the time t the
experimental data (cf. [
In this paper we are going to study biochemical reaction systems which possess a steady situation1. This is the case when the system comes in a situation, in which the concentration of all substances stays constant. Usually, in the steady situation1 the concentration of all substances allows that all reactions take place permanently. Now constructing the reachability graph, may be interpreted as considering the path of changes of the single substances. Loosely speaking, we should get a cyclic set of states in the reachability graph corresponding to the steady situation1. The behaviour of the SPN, expressed by Markow chains is isomorphic to the full reachability graph of the underlying PN, i.e., they have the same state space. By convention, we speak in the following of \the reachability graph" of the SPN. But the nodes corresponding to the steady situation1 can not be recognized in this reachability graph, even knowing the rates in the steady situation1. To our knowledge, no method is known until now for separating or extracting the subgraph corresponding to the steady situation1 from the reachability graph of the SPN. Steady state meaning cyclic behaviour, this subgraph (of the reachability graph of the SPN) is supposed to present a cyclic structure (with one or more circles), up to some initiation part. In contrast, knowing the steady situation1 rates we are able to separate the reachable states, we are looking for, using a TPN and its reachability graph. This is due to the reduction results on reachability graphs of TPNs discussed in section 2.3.
Thus, we propose in this paper a methodologie to calculate and verify such set of states which correspond to steady situation1. The starting point will always be an SPN model for a biochemical reaction system, which has a steady situation1. This means that the rate for each enabled transition in each marking is constant. The steady situation1 concentrations and rates can be determined using simulation of the SPN. Examples for which about 10,0000 simulation runs have been done may be considered. These runs has to be merged into one averaged simulation run showing the mean of the concentrations, and thus also of the rates, over the time. We take the expectation values of the steady situation1 rates for our investigations.
Simulation means approximation; thus it is not a priori clear how accurate the determined steady situation1 rates are.
The reciprocal value of the rate is the time which each enabled transition has to wait before it can re. The transition with the minimal waiting time res in an SPN. Consequently, an SPN acts in the steady situation1 exactly like a certain kind of TPN: we propose to construct a TPN, having the same underlying Petri net as the SPN, and where the transitions t will recieve time intervals [at; bt], where at = bt is equal to the above calculated waiting time of t in the SPN in the steady-state.
A qualitative analysis of the TPN can prove whether the subset of all reachable states generate cycle(s) only (up to some initiation part). Furthermore, the time length of these cycles can be computed.
The formal analysis we proposed allows the following interpretations. If the reduced reachability graph of the TPN consists of cycles, then the considered rates achieved by simulation describe a steady situation1, actually. By contraposition we may deduce, that a non cyclic form of the reduced state space indicates severe problems in the set of data used to build the TPN, and by consequence, in the initially established data from experiences in the wet labs.
Additionally the time-length of the cycles can be easily computed and compared with results from the wet labs. Both, the reachability graph of the TPN and the time-length of the cycles are either an indication for the correctness of the models or they invalidate these. Therefore our method o ers a way of validating a complex modeling process of biochemical reaction systems. 4
In this section we will illustrate our approach of analysing a biochemical reaction
system in a steady situation1 along an example introduced in [
In [
Afterwards, simulation studies provides a qualitative validation of the mathematical model compared to experimental results in the wet labs in view of the transient behavior and sensitivity analysis. However, parameter estimation is, as already mentioned, an uncertain factor in such a mathematical model.
Then in [
Additionally, reaction rates are associated to all transitions of this PN [
We are going to investigate the SPN in the simulated steady situation1. First let us have a look on its reachabilty graph depicted on the left hand side of Fig.4. Unfortunately no method exists, to our knowledge, to nd out analytically ERK-PP s9 s8 MEK-PP_ERK
r6 s7 MEK-PP r1
Raf-1Star_RKIP r8 r3
r4 s4 Raf-1Star_RKIP_ERK-PP
s11 r7 s5
ERK r5 r9
r10 s6 RKIP-P s10
RP
r11
RKIP-P_RP
Fig. 3. The Petri net forPUYtRheORYcDoHrOYeM mNYBModCeNSVl oSNCfFthCYOeN RSYCKIFPNT0 pTaNF0thFwNP0aPyNF,0cNoECnSsisting of 11 places
and 11 transitions. The pDlTaPceCPsIsC1T,I..S.C,TsI1S1Bskt-aBn1d-BfoDrCFpDrSotteDiTnrsLoIVr pRErVotein complexes.
Complexes are indicated by aYn uYndeYrscoNre Y bYetwYeenN th0e pNrotYeinYnames, phosphorylated
forms by the su x -P or -PP. The transitions r1, ..., r11 model the reactions. The
preplaces of a transition correspond to the reaction's precursors, and its postplaces to
the reaction's products. The layout follows the suggestions by the graphical notation
used in [
the rate function vi, presented in [
s11
Now, we can calculate the rates v1 v11 of the transitions in the steady situation1 using their rate functions from Table 1 and the data from Table 1 and Table 2. Subsequently, the delay time i for every one of the ten transitions ri is obtained as the reciprocal of the rate in the steady situation1. The resulting values are presented in Table 3.
Finally, the TPN model can be constructed: As skeleton of the TPN model we take the underlying PN of the SPN, i.e. the net given in Fig. 3. To each transition ri transition ri; 1 i 11, a time interval [ i; i] is associated, where i is the calculated delay time, from Table 3. r5
Now the obtained TPN may be analysed. First, we just calculate the reachable p-markings, designed on the right hand side of Fig.4. We observe that due to the time constraints this graph has 9 nodes, i.e., much less p-markings are reachable as in the reachability graph of the SPN, depicted on the left hand side, which is also the reachability graph of the underlying net of the SPN and TPN. We also detect that the right graph is clearly a subgraph of the left one. The complete state space of the TPN is -a priori- in nite, a lot of states may share the same p-marking.
The reduced reachability graph of the considered TPN can now be
constructed, by applying the reduction method described in section 2.3. We did it
with tools INA and Charlie [
Analyzing this reduced reachability graph of the TPN tells us that it consists of the cycle z4; r3; z5; r5; z6; r6; z7; r9; z8; r8; z9; r11; z10; r1 and an initiation path z1; r6; z2; r8; z3; r1:
This path (or panhandle) is caused by the choice of the initial p-marking for the TPN, chosen to be the same as for the SPN. Actually, the initial pmarking for the TPN should be a p-marking which the SPN reaches in the steady situation1. However, the TPN only initiates its behaviour by this path and then comes to the steady situation1, i.e. stays in the cycle. The time-length of the cycle was also calculated, its value is 731 time units.
We also read on this reduced graph that the transitions r2; r4; r7 and r10 will never re. Such transition are called dead. These are the transitions modeling the backward reactions which have rate constants being essentially smaller as the rate constants for the forward reactions. This tell us that in the steady situation1 the backward reactions do never proceed. 5
In this paper we introduce a method for qualitative and quantitative evaluation of an SPN model of biochemical reaction systems in a steady situation1 including validation of all used data. A mathematical model of such a system contains a number of uncertain factors resulting from the estimation of the parameters in the ODEs and the values of the reaction rates in the steady situation1 obtained by simulation of the SPN. The reaction rates and the concentration of the species in the steady situation1 are constant values.
This means that the set of reachable markings in the SPN model in the steady situation1 is nite and they generate a cycle, not necessarily a simple one. But no state reachability analysis of the SPN does allow for isolating those states which correspond to the steady situation1, i.e. does allow to detect the cycle.
Due to the fact of constant values, we are able to propose a mapping from the usual SPN model in the steady situation1 onto a TPN model which has the same behaviour. Contrarily to the SPN model, we can reduce the state space of the TPN to the part we are interested in, consisting of the essential states. The obtained reduced state graph can be further analyzed. Its cyclic or non cyclic form validate or invalidate the used data during the modeling process. The time length of the cycle can be calculated and compared to real time measures, too.
The algorithms for reachability analysis of TPNs, implemented in the tools
[
We will lead some re exions if the initial state for the TPN could be rede ned in a better way by regarding the values for the concentrations of the species in the simulated steady situation1. We are planning to apply the presented method to some other cases of biological or biochemical interaction networks, were more complex steady situations1, with -a priori- non simple cycles.