Recently hybrid modelling and simulation of biochemical systems have attracted increasing interest. This is motivated by the need of simulating systems which integrate di erent sub-cellular models, and the fact that bio networks themselves are inherently stochastic, however stochastic simulation is time expensive. Compared to other methods of biological modelling, Petri nets are characterized by their intuitive visual representation and executability of biological models. In this paper, we present a hybrid Petri net class that incorporates both continuous and stochastic capabilities. The presented class is intended to model and simulate hybrid biological systems such that they contain some parts which are simulated deterministically while other parts are simulated stochastically.
Computer simulation is an essential tool for studying biochemical systems. The deterministic approach (continuous simulation) is the traditional way of simulating biochemical pathways. In this approach, reactions and their in uence on the concentrations of the involved species are represented by a set of ordinary di erential equations (ODEs). The changes in reactants and products are obtained through solving the resulting ODEs using numerical integration algorithms. While this approach has the advantage of a well established mathematical basis and strong documentation, it lacks to capture the phenomena which occur due to the underlying discreteness and random uctuation in molecular numbers [Pah09],[LCP+08], especially in situations where the number of molecules is few.
Stochastic simulation [Gil76] provides a very natural way of simulating biochemical pathways, since it can successfully capture the uctuations of the underlaying model. Furthermore it deals correctly with the problem of extremely low number of molecules [ACT05]. In stochastic simulation, species are no longer represented as continuous concentrations which change continuously with time, instead they are represented as discrete entities such that their dynamics can be simulated using the machinery of Markov process theory. In [SYS+02] an example is given comparing deterministic versus stochastic modeling using a simple model of the intracellular kinetics of a generic virus.
A major drawback of the stochastic simulation is that it is computationally expensive, when it comes to simulate larger biological models [Pah09],[LCP+08], [ACT05], especially when there are a large number of molecules of some chemical species. The reason behind this problem comes from the fact that we have to simulate every reaction event when we use stochastic simulation to simulate biological systems [LCP+08]. This drawback motivates scientists to search for other methods to enhance the capability of the stochastic approach. Hybrid simulation is one of these methods.
Hybrid simulation [ACT05],[Kie+04],[Rue+07] of biochemical system using both deterministic and stochastic approaches has been recently introduced to take the advantage of capturing the randomness and uctuation of the discrete stochastic model and allows at the same time a reasonable computation time. This goal is achieved by simulating fast reactions deterministically, while simulating slow reaction stochastically. While this method provides a promising approach for simulating biochemical models, there are some open questions which need to be solved [Pah09].
Petri nets provide a very useful way of modelling biochemical pathways [RML93],[BGH+08],[HGD08],[Mat+03] since they provide an intuitive approach of transforming the biological model into a graphical representation which coincides with the qualitative description of this model. Furthermore, they can be easily transformed later for quantitative simulation.
Continuous Petri nets are used in biological modelling to introduce an easy way of modelling complex biological pathways and simultaneously hide the mathematical complexities of the underlying ODE. Contrary, in stochastic Petri nets and their simulation, transitions re with exponentially distributed random waiting time.
Hybrid Petri nets [AD98] incorporate both continuous and discrete capabilities and can be used to model systems which contain both discrete and continuous elements. Many various of hybrid Petri nets have been introduced during the last two decades, with di erent modeling goals. Some examples can be found in [Mat+03],[TK93] and [PB09]. An overview of continuous, discrete and hybrid Petri nets can be found in [DA10] .
In this paper, we introduce the de nition of a hybrid continuous-stochastic Petri net, HPN, and integrate it into Snoopy [HRR+08],[RMH10], a tool to design and animate or simulate hierarchical graphs, among them the qualitative, stochastic and continuous Petri nets, which incorporate the modeling capabilities of the previously introduced stochastic and continuous classes [GH06], [GHL07],[HLG+09], The new net class HPN, is intended to model biological pathways that require hybrid simulation, such that the resulting Petri net can be simulated deterministically and stochastically based on the model speci cation.
This paper is organized as follow: Firstly we brie y review the motivations of using continuous and stochastic Petri nets to model biochemical reactions. Then we introduce our hybrid stochastic-continuous Petri net class, by rstly presenting a formal de nition as well as the connectivity rules between its elements. The illustration of the modeling capabilities of HPN to model biological systems is then demonstrated using two examples. At the end we conclude by a summary and autlook of future work. 2
The tight analogy between Petri net and biochemical reactions makes it a natural choice to model these reactions [RML93],[HGD08]. Being bipartite, concurrency, and stochasticity are common properties shared by Petri nets and biochemical interactions. Qualitative Petri net [HGD08] can be used to analyze the biochemical systems qualitatively, while stochastic and continuous Petri nets are used to simulate them quantitatively. Before we discuss the various aspects of the hybrid stochastic-continuous Petri net, we provide a short overview of continuous and stochastic Petri nets as well as how they can be used to model biological systems. Detailed discussion can be found in [BGH+08] and [HGD08], and for a general introduction to Petri net see [DA10] and [Mur89].
Continuous Petri nets provide a way for modeling systems in which states change continuously with time. In this class of Petri nets, places contain nonnegative real values and transitions re continuously with time. In systems biology, continuous Petri nets provide a very useful way of representing ODEs. Preplaces of the transitions represent reactants species and the marking of these places represents species' concentrations. Each transition is associated with a rate function which de nes the kinetic rate. The corresponding ODE which represents the reaction which is modeled by this transition can be generated using (1) [GH06]. dp dt = X f (t; p)v(t)
X f (p; t)v(t) t2 p t2p (1) where v(t): is the rate function and f (t; p): is the weight connecting transition t with place p and p; p are the pre- and post-transitions of place p, respectively. Note that place names are read as real variables.
The resulting system of ordinary di erential equations of all places describes the changes with respects to time in all biochemical species. Our HPN supports the same functionality as the aforementioned continuous Petri net.
In contrast to continuous Petri nets, stochastic Petri nets preserve the discrete state description. The biochemical models are simulated stochastically by associating a probability-distributed ring rate (waiting time) with each transition. This means that there is a time which has to elapse before an enabled transition t 2 T res [HLG+09], where T is the set of all stochastic transitions. The probability density function of the exponentially distributed random variable, xt, which represents the waiting time, is given by (2) fxt ( ) = t(m):e t(m) ; t 0 (2) where t(m) is a marking dependent kinetic rate which is associated with each stochastic transition. t(m) is equivalent to the propensity of the reaction t, a(xi), of the stochastic simulation algorithms which are presented in [Gil76].
Because of the deterministic nature of continuous Petri nets, the concentration of particular species will have the same values at each time point for repeated experiments, which is the main di erence between simulation of stochastic and continuous biological models, and hence for Petri nets as well. In a typical execution of stochastic Petri nets, each transition gets its own local timer. When a particular transition becomes enabled, the local timer is set to an initial value which is computed by means of the corresponding probability distribution. The local timer is then decremented at a constant speed and the transition will re when the time reaches zero. A race will take place in the case of con ict between more than one enabled transition.
To extend the modeling capabilities of stochastic Petri nets (SPN) in biological system, two extensions, general stochastic petri nets (GSP Nbio) and deterministic stochastic petri nets (DSP Nbio), of SPN are introduced in [HLG+09]. These extensions add inhibitor and read arcs and deterministically time-delayed transitions to stochastic Petri nets.
In the following section, we present the merging of stochastic Petri nets (using the extended version) and continuous one, to produce a hybrid continuousstochastic Petri nets which are capable of modeling and simulating hybrid biochemical reactions. 3
In this section we describe the hybrid continuous stochastic Petri nets capable of modeling systems which consist of discrete and continuous parts. The discrete parts may be considered as a set of reactions which involves species with low number of molecules such that it is adequate to simulate them in a discrete way. On the other hand, continuous elements of this class can represent a set of reactions which involves species with large number of molecules, which are computationally too expensive to be simulated stochastically. Continuous and stochastic Petri nets complement each other. We get modelling power of uctuation and discreteness, when using the stochastic simulation and at the same time we can simulate the computationally expensive parts deterministically using ODEs solvers.
Generally speaking, biochemical systems can involve reactions from more than one type of biological networks, for example regulatory, metabolic or transduction pathways. Incorporation of reactions which belong to distinct (biological) networks, tend to result in sti systems. This follows from the fact that regulatory network's species may contain a few number of molecules, while metabolic networks' species may contain a large number of molecules [Kie+04]. In our hybrid Petri nets, reactions which involves species with a small number of molecules are represented by discrete entities, so that they can be simulated stochastically, while reactions which include a large number of molecules are represented by continuous entities, so that they can be simulated deterministically. The connection between the discrete and continuous parts takes place using either special arcs (read, inhibitor, or equal arcs) or in some cases using the standard arcs based on the de ned connection rules.
In the rest of this section, we will discuss in more detail the newly introduced hybrid continuous-stochastic Petri nets in terms of the graphical representation of its elements as well as the ring rules and connectivity between the continuous and stochastic parts. 3.1
As expected, HCSPN contains two types of places: discrete and continuous. Discrete places (single line circle) contain integer numbers which represent for example the number of molecules in a given species. On the other hand, continuous places - which are represented by shaded line circle - contain real numbers which represent the concentration of a given species. This means that we can combine the power of the previously discussed continuous and stochastic Petri nets together in one class. HCSPN contains a variety of transition types: continuous, stochastic, deterministic, immediate, and scheduled transitions [HLG+09]. Continuous transitions - shaded line square - re continuously in the same way like in continuous Petri nets. Their semantics are governed by ordinary di erential equations. Their ODEs de ne the changes in the transitions' pre- and post-places.
Stochastic transitions which are drawn in Snoopy as a square, re randomly with an exponential random distribution delay. The user can specify a set of ring rate functions, which determine the random ring delay. Deterministic (time delay) transitions - black square - re after a speci ed time delay, immediate transitions - black bar - re with zero delay, and they have higher priority in the case of a con icts with other transitions. They may carry weights which specify the relative ring frequency in the case of con icts between more than one immediate transition. Scheduled transitions - grey square - re at a userspeci ed time point or time interval.
The connection between those two types of nodes (places and transitions), takes place using a set of di erent arcs. HCSPN contains ve types of edges: standard, inhibitor, read, equal and reset arcs. Standard edges connect transitions with places or vice versa . They can be continuous, i.e carry real value weights (or in the biochemical context stoichiometry), or discrete i.e carry nonnegative integer value weights. Special arcs like inhibitor, read, equal and reset arcs provide only connection from places to transitions, but not vice versa. The connection rules and their underlying semantics are given below. Fig. 1 provides a graphical illustration of those elements. While this graphical notation is the default one, they can be easily customized using our Petri nets editing tool, Snoopy.
Discrete
Continuous
Stochastic
Continuous
Immediate Deterministic <?> _SimStart,?,_SimEnd
Scheduled Discrete or Continuous Inhibitory Read
Equal
Reset
{ P = fPcont [ Pdiscg whereby Pcont is the set of continuous places to which nonnegative real values can be assigned and Pdisc is the set of discrete places to which nonnegative integer values can be assigned. { T = Tcont [ Tstoch [ Tim [ Ttimed [ Tscheduled with: 1. Tcont, the set of continuous transitions, which re continuously over time. 2. Tstoch, the set of stochastic transitions, which re stochastically with exponentially distributed waiting time. 3. Ttimed, the set of deterministic transitions, which re with a deterministic time delay. 4. Tscheduled, the set of scheduled transitions, which re at prede ned ring time points. 5. Tim, the set of immediate transitions, which re with waiting time zero and it has higher priority compared to other transitions. { A = fAcont [ Adisc [ Ainhibit [ Aread [ Aequal; [Aresetg, is the set of directed edges, whereby: 1. Acont : ((Pcont T ) [ (T Pcont)) ! IR0 : de nes the set of continuous, directed arcs, weighted by nonnegative real values. 2. Adisc : ((P T ) [ (T P )) ! IN0 : de nes the set of discrete, directed arcs, weighted by nonnegative integer values. 3. Aread : (P T ) ! IR+if P 2 Pcont or Aread : (P T ) ! IN+if P 2 Pdisc, de nes the set of read arcs. 4. Aequal : (P T ) ! IR0+if P 2 Pcont or Aequal : (P T ) ! IN0+if P 2
Pdisc, de nes the set of equal arcs. 5. Ainhibit : (P T ) ! IR+ [ f0+gif P 2 Pcont or Ainhibit : (P T ) ! IN+if P 2 Pdisc, de nes the set of inhibits arcs, where 0+ means very small positive real number but not zero. 6. Areset : (P Tdiscrete) de nes the set of reset arcs, where Tdiscrete =
Tstoch [ Tim [ Ttimed [ Tscheduled is the set of discrete transitions. { V is a set of functions ff,g,d,wg where : 1. f : Tcont ! Hc is a function which assigns a rate function hc to each continuous transition t 2 Tcont, such that : fhct jhct : IRj0 tj ! IR+; t 2 Tcontg is the set of all rates functions and f (t) = hct ; 8t 2 Tcont. 2. g : Tstoch ! Hs is a function which assigns a stochastic hazard function hst to each transition t 2 Tstoch, whereby fhst jhst : INj0 tj ! IR+; t 2 Tstochg is the set of all stochastic hazard functions and g(t) = hst 8t 2 Tstoch . 3. d : Ttimed ! IR+, is a function which assigns a constant time to each deterministic transitions representing the waiting time. 4. w : Tim ! Hw is a function which assigns a weight function hw to each immediate transition t 2 Tim, such that : fhwt jhwt : INj0 tj ! IR+; t 2 Timg is the set of all weight functions and w(t) = hwt ; 8t 2 Tim { m0 = fmcont [ mdiscg : is the set of initial marking for both the continuous (Pcont) and discrete places (Pdisc), whereby mcont 2 IR0+jPcontj, mdisc 2 IN0+jPdiscj.
A critical question arises when considering the mixing between discrete and continuous elements: how are these two di erent parts connected with each other? Fig. 2, provides a graphical illustration of how the connection between di erent elements of the introduced HCSPN takes place. Note that other discrete transitions (immediate, deterministic and scheduled transitions) follow the same connection rules as stochastic transitions.
Firstly, we will consider the connection between continuous transitions and the other elements of the HCSPN. Continuous transitions can be connected with continuous places in both directions using continuous arcs (i.e arc with real value weight). This means that continuous places can be pre- and post-places of continuous transitions. These connections represent deterministic, biological interaction. According to the previous formal de nition, each continuous transition takes a rate function. This rate function represents the kinetics of the deterministic reaction. Like in continuous Petri net, the ring of this transition can be represented as an ODE. The continuous transition can be connected also with a discrete or continuous places, but only by one of the special arcs (inhibitor, read, equal). Read arcs allow to specify positive side conditions, while inhibitor arcs allow to specify negative side conditions. It is worth being mention, that the markings of the transition preplaces connected by these special arcs do not change when the transition res. This type of connection allows a connection between the discrete and continuous parts of the biochemical model.
Discrete places are not allowed to be connected with continuous transitions using standard arcs, because the rings of continuous transitions are governed by an ODE which requires real values in the pre- and post-places. Discrete transitions (stochastic, deterministic, immediate and scheduled) can be connected with discrete or continuous places in both directions using standard arcs. However, the arc's weight should be considered, i.e the connection between discrete transitions and discrete places takes place using arcs with nonnegative integer numbers, while the connection between continuous place and discrete transitions is weighted by nonnegative real numbers. The general rule to determine the weight type of the arcs is the type of the transition's pre/post places.
The connection between continuous places and discrete transitions will result in a model like discussed in [TK93], in which the changes in the continuous places are governed by ring of stochastic transitions. Discrete transitions can also have discrete or continuous places as the transition pre-places using the special arcs. 3.3
Due to the use of both stochastic and continuous parts in HCSPN, we have now two di erent clocks: one for the continuous parts and the other for the stochastic ones. The ODEs solver which represents the semantics of the continuous Petri net evolves deterministically with approximate time steps, while the stochastic transitions re stochastically with exact time steps. Because we intend to use HCSPN to simulate biochemical reactions, we provide a synchronization mechanism between the stochastic and continuous Petri nets, since some species (places) may belong simultaneously to both continuous and stochastic Petri nets due to the partition of the reactions. In this part of the paper we propose a Petri net interpreted synchronization algorithm based on the algorithm presented in [ACT+04].
Many synchronization algorithms are used in the literature to synchronize between the deterministic regime and the stochastic one in hybrid simulation of biochemical reactions; some of them can be found in [Pah09,ACT05,Kie+04,Rue+07]. We opted to use the algorithm in [ACT+04], since it has a rigid mathematical basis for the synchronization of the two di erent clocks.
The algorithm which is presented here is based on the direct method [Gil76], see [ACT05,ACT+04] for other variations based on the rst and second reaction method. The algorithm is based on the function f ( jt). f ( jt) will decide when we can switch from the continuous world to the stochastic one. We rstly draw an exponentially distributed random variable and initialize f ( jt) = 0, then we start to simulate the continuous transitions using the ODE solvers. During the 0.1 or 0.1 or 0.1 or 0.1 or 0.1 or 0.1 or 0.1 or 0.1 or 0.1
Continuous transition continuous simulation, f ( jt) will be increased according to the time evolution of the ODE presented in (3) d
f ( jt) = dt
X j2Tstoch gj (m( ); ) (3) where gj (m( ); ) is the rate function, which is associated with each stochastic transition and was de ned in the aforementioned formal de nition of the HCSPN, and m( ) is the current marking of the transition's pre-places. We repeat the continuous simulation until time = s such that f ( jt) = . The mathematical derivation which is presented in [ACT+04] proves that a stochastic event will occur at time = s, which means that we can execute the stochastic simulation at that time. Then we update the current marking according to the red transitions using the arcs' weights which connect this red transition with their pre-places and then we advance the simulation time. The previous steps are then repeated until we reach the end of simulation time. In the following we present the algorithm in a more formal way. 1. Start by the initial marking m0 and the initial time t = t0; 2. Generate an exponentially distributed random variable . 3. Set g( jt) = 0 and simulate the continuous transitions using the ODE solver starting at time = t and progress g( jt) according to equation (3) Until time = s such that g( jt) = . 4. Perform the stochastic simulation using the discrete transitions. 5. Update the current marking m(t) according to the red transitions. 6. Repeat steps 2-5 until we reach the end of simulation time . 4
In this section, we demonstrate by examples how the HCSPN is used to model biological systems. The two examples which are presented here are: the genes operons model and the modeling of the role of LL-6R in regulation of early haematopoiesis. 4.1
In this example, we model two genes operons using the HCSPN class. The original model can be found in [MDN+00]. The HCSPN in Fig. 3 describes the transcription of an operon containing two genes. The two genes are represented by two discrete places, Gene1, Gene2, respectively. The transcription of Gene one is represented by the transition transcriptionG1, which is a stochastic transition. This transition is associated with a ring rate function, which determines when this transition res. After the transcription took place, an amount of concentration which represents the mRNA of Gene one is added to the continuous place mRN A1. This concentration value is equal to the rate function of the continuous transition, transcriptionG1, multiplied by the weight of the arc connecting transition transcriptionG1 with place mRN A1. The concentration of the mRNA of Gene1 can be degradated continuously, when transition deg1 res, if the value of the place mRN A1 is greater than zero. A process called translation can take place depending on the concentration of mRNA. However this process does not change the concentration's value of the mRN A1 value. So we choose to connect them using a read arc.
After the translation process took place, the protein of Gene one which is represented by the continuous place P rotein1 can be degraded, when the transition labeled Degprotein1 res. A similar story can happens to Gene two after the polymerase of the RNA of Gene one into Gene two. The ring rate functions of the stochastic transitions and the rates of the continuous transitions can be speci ed by the user by selecting between a set of kinetic rate functions among them is the mass action kinetics. This example demonstrate by a simple way the modeling power of the HCSPN in system biology. RNA PolyGen1
Gene1 Gene2
TranscriptionG1 mRNA1 0
TranslationG1 After we presented a simple example to illustrate the di erent elements of the HCSPN class, in this section we present a more realistic biological example, modeling the role of a speci c cytokine, interleukin-6, in the regulation of early hoematopoiesis [TTC+06]. Fig.4 shows the modeling of this pathway using the HCSPN Petri net. Haematopoiesis is a complex phenomena beadings to the continuous production of all types of mature blood cells. The use of hybrid Petri nets to model the regulation of early haematopoiesis is motivated by the need of discrete elements for modeling the cellular evaluation, as well as continuous elements to model molecular interactions [TTC+06].
Consequently, the model of the IL-6R regulation of the early haematopoiesis consists of two submodels: the cellular submodel and the molecular one. In the former the three di erent cells types, equiescent, permissive, and committed cells are modeled by three discrete places, Pq, Pp, and C, respectively. Deterministic transitions are used to model the biological processes which take place between these cells types. In the later submodel, continuous places model the molecules involved in the regulation of the haematopoiesis by IL-6, while biological processes are modeled using continuous transition. The bright gray arcs represent the positive feedback loop involving the sLL-6R. Note that in the cellular submodel, arcs weight equal to one are not displayed.
The resulting hybrid Petri net model can be simulated ( continuously and stochastically). Because there are no stochastic transitions in this model, the stochastic simulation is simpli ed to simulate the ring of the discrete Petri net submodel. 5
In this paper we have presented our research in progress of de ning and implementing a hybrid continuous stochastic Petri net class which includes both T14
0
T7 sIL_6R
T11
T9
IL_6 0
T21 discrete and continuous modeling capabilities of biochemical interactions. The presented class is intended to model systems which are sti , i.e contain some species with high number of molecules as well as species with low number.
Snoopy supports the export of drawn models to many other tools. For the hybrid class it can be exported to Modelica's hybrid Petri net library [PB09] for further simulation.
Our hybrid model is based on xed partitioning of the biochemical system, i.e. the reactions are initially divided into discrete and continuous parts. Further extension of this work aims to permit the dynamic partitioning of the reactions during the simulation based on some criterias like the number of molecules in each species or the reaction propensity.
Mostafa Herajy is supported by the GERLS (German Egyptian Research Long Term Scholarships) program, which is administered by the DAAD in close cooperation with the MHESR and German universities.