How to Model and Simulate Biological Pathways with Petri Nets A New - - PDF document

SLIDE 1

How to Model and Simulate Biological Pathways with Petri Nets – A New Challenge for Systems Biology – ∗

Satoru Miyano

Human Genome Center, Institute of Medical Science, University of Tokyo miyano@ims.u-tokyo.ac.jp

Hiroshi Matsuno

Faculty of Science, Yamaguchi University matsuno@sci.yamaguchi-u.ac.jp

1 Why Petri nets to model biological pathways ?

1.1 Petri nets – suitable than other mathematical descriptions

In 1999, we surveyed which architecture is suitable when modeling and simulating biopathway for biological and medical scientists. At that time, there were ODE-based attempts for modeling and simulating chemical reactions

• e.g. Gepasi [19], E-Cell[29] - and others - e.g. the Lisp based architecture QSIM [13] and our other work,

the pi-Calculus based architecture, Bio-Calculus [21]. Unfortunately, applications based on these architectures are not acceptable in their fields. This is due to poor GUI interfaces, e.g. lacking biopathway editors, or their architectures themselves. To overcome this situation, we came to the conclusion that an architecture based

• n Petri nets should be suitable because of their intuitive graphical representation and their capabilities for

mathematical analyses.

1.2 Various types of Petri nets have been used for describing biological pathways

Researches on Petri nets have a long history of nearly 40 years from the paper by Dr. Petri [25]. The first attempt to use Petri nets for modeling biological pathways was made by [26], giving a method to represent metabolic

• pathways. Hofest¨

adt [9] expanded this method to model metabolic networks. Subsequently, several enhanced Petri nets have been used to model biological phenomena. Genrich et al. [5] modeled metabolic pathways with a colored Petri net by assigning enzymatic reaction speeds to the transitions, and simulated a chain of these reactions quantitatively. Voss et al. [30] used the colored Petri net in a different way, accomplishing a qualitative analysis of steady state in metabolic pathways. The stochastic Petri net has been applied to model a variety of biological pathways; the ColE1 plasmid replication [7], the response of the σ32 transcription factor to a heat shock [27], and the interaction kinetics of a viral invasion [28]. On the other hand, we have shown that the gene regulatory network of λ phage can be more naturally modeled as a hybrid system of “discrete” and “continuous” dynamics [15] by employing a hybrid Petri net

• architecture. It has also been observed in [6] that biological pathways can be handled as hybrid systems. For

example, protein concentration dynamics, which behave continuously, being coupled with discrete switches. Another example is protein production that is switched on or off depending on the expression of other genes, i.e. the presence or absence of other proteins in sufficient concentrations.

∗This article is distributed at the tutorial in the 25th International Conference on Application and Theory of Petri Nets to be

held in Bologna, Italy, on June 22, 2004.

1

SLIDE 2

normal arc test arc inhibitory arc discrete transition discrete place continuous transition continuous place

Figure 1: Elements of hybrid (functional) Petri net.

2 Hybrid functional Petri net

2.1 Hybrid functional Petri net includes hybrid Petri net, hybrid dynamic net, and functional Petri net

Petri net is a network consisting of place, transition, arc, and token. A place can hold tokens as its content. A transition has arcs coming from places and arcs going out from the transition to some places. A transition with these arcs defines a firing rule in terms of the contents of the places where the arcs are attached. Hybrid Petri net (HPN) [2] has two kinds of places discrete place and continuous place and two kinds of transitions discrete transition and continuous transition. A discrete place and a discrete transition are the same notions as used in the traditional discrete Petri net [26]. A continuous place can hold a nonnegative real number as its content. A continuous transition fires continuously at the speed of a parameter assigned at the continuous transition. The graphical notations of a discrete transition, a discrete place, a continuous transition, and a continuous place are shown in Figure 1, together with three types of arcs. A specific value is assigned to each arc as a weight. When a normal arc is attached to a discrete/continuous transition, w tokens are transferred through the normal arc, in either of normal arcs coming from places or going out to places. An inhibitory arc with weight w enables the transition to fire only if the content of the place at the source of the arc is less than or equal to w. For example, an inhibitory arc can be used to represent repressive activity in gene regulation. A test arc does not consume any content of the place at the source of the arc by firing. For example, a test arc can be used to represent enzyme activity, since the enzyme itself is not consumed. Hybrid dynamic net (HDN) [4] has a similar structure to the HPN, using the same kinds of places and transitions as the HPN. The main difference between HPN and HDN is the firing rule of continuous transition. As we can know from the description about HPN above, for a continuous transition of HPN, the different amounts

• f tokens can be flowed through the two types of arcs, coming from/going out the continuous transition. In

contrast, the definition of HDN does not allow to transfer different amount through these two types of arcs. However, HDN has the following firing feature of continuous transition which HPN does not have; “the speed

• f continuous transition of HDN can be given as a function of values in the places”.

From the above discussion, we can know that each of HPN and HDN has its own feature for the firing mechanism of continuous transition. As a matter of fact, both of these features of HPN and HDN are essentially required for modeling common biological reactions. (See the example of Figure 7 in which four monomers compose one tetramer.) This motivated us to propose hybrid functional Petri net (HFPN) [16] which includes both of these features of HPN and HDN. Moreover, HFPN has the third feature for arcs, that is, a function

• f values of the places can be assigned to any arc. This feature was originated from the functional Petri net

[10] which was introduced in order to realize the calculation of dynamic biological catalytic process on Petri net based biological pathway modeling. Formal definition of the HFPN is given in [22].

2.2 How to use discrete and continuous elements of HFPN

Biological pathways essentially consist of discrete parts such as a genetic switch control and continuous parts such as a metabolic reaction. These discrete and continuous parts can be represented by discrete elements (discrete place and discrete transition) and continuous elements (continuous place and continuous transition) of

• HFPN. For example, a control system turning on or off the gene expression with operator site can be represented

by discrete elements. That is, if the discrete place has a token, the protein necessary for activating the operator site has bound to the operator, that means the gene expression turn on. In addition, using the delay concept

• f the discrete transition, we can easily describe the transcription which happens after a certain time.

2

SLIDE 3

m1 m2 m3 m3 m1 m2

v2 v1

m5 m4 m6

v3

m5 m4 m6

v4

m7

v5

m10 m11 m12

v8

m8 m9

v6 v7

protein A protein A protein A signal cited protein A activation protein A inactivation protein A protein A protein B protein A' protein A protein A' protein A protein B complex AB complex AB (a)Protein complex formation (b)Enzyme reactions (c)Protein degradation (d)Protein State switching (e)Protein transformation active enzyme inhibitor composition protein degradation inactivate activate signal transmit m13 m14

v9

nucleate A extranuclear A (f)Substrate migration move out nuclere decomposition constant concentration competitive inhibition

Figure 2: HFPN components for typical biological reactions. Contents of continuous places represent concen- trations of substances such as proteins. (b) Test arc is used for an enzymic reaction, since enzymes are not consumed by reactions. Two types of enzymic reactions “activating enzyme” and “blocking enzyme” are de-

• scribed. Furthermore, activating enzyme reactions can be classified into two patterns “constant concentration”

and “variable concentration.” Reaction speeds of them are given in Table 1. (e) Discrete elements can be used for state changing reactions instead of constant elements. In this case, when a place gets a number of tokens enough to fire the transition, the reaction proceeds. On the other hand, biological phenomenon such as transcription, translation, and enzymic and metabolic reactions have been treated as events whose conditions change continuously. By modeling the transcription and the translation with continuous elements in the following way, expression levels of mRNA and proteins can be simulated. Continuous places are used for storing concentrations of mRNA and protein. Reaction speeds

• f transcription and translation are assigned at parameters of continuous transitions. In addition, common

formulas for biological reactions such as Michaelis-Menten’s equation can be modeled almost directly with assigning concentrations of substrate and product to continuous places and the formula for the reaction of them to the continuous transition between these two places.

2.3 HFPN components for basic biological reactions

Many biological pathways can be constructed by mainly using the following six basic reactions: (1) protein composition and decomposition, (2) enzymic reaction, (3) protein disintegration, (4) state switching, (5) state changing, and (6) substance migration. HFPN components of these reactions are listed in Figure 2, and reaction speeds of these reactions are summarized in Table 3.1. The variables a, b, c1, c2, d, e, and f are constants. Depending on the speed of biological reactions to be modeled, the appropriate integer in the range of 5 to 20 is assigned to each of six variables. The variables Vmax, Km, and Ki are the maximum speeds, a Michaelis constant, and a blocking constant of a Michaelis-Menten equation [1], respectively. In our model, we let Vmax = 2 and Km = 1, and use a real number in the range of 0.4 to 0.5 for the Ki. 3

SLIDE 4

Table 1: Transition speeds for HPN components of typical biological reactions in Figure 2 Reaction Transition speed composition v1 = m1∗m2

a

protein formation decomposition v2 = m3

a

non-constant concentration v3 = m4∗m5

b

constant concentration v3 = Vmax∗m4

Km+m4

enzymic reactions competitive inhibition v4 =

Vmax∗m5 Km(1+ m6

Ki )+m5

protein degradation v5 = m7

d

protein state switching v6 = m8

c1 v7 = m9 c2

protein transformation v8 = m11

e

substrate migration v9 = m13

f

4

SLIDE 5

3 lac Operon Gene Regulatory Mechanism and Glycolytic Pathway

This section demonstrates how we can model biopathways by using HFPN. As an example, we adopt lac operon gene regulatory mechanism and glycolytic pathways in E. coli and describe HFPNs with reference to biological facts in the literature [1, 14, 31]. The HFPN model of “transcription control switch” (Figure 4) is firstly described, then this model is going to be extended with adding “positive regulation” (Figure 5), “negative regulation” (Figure 6), and “hydrolyzing lactose to glucose and galactose” (Figure 8). For the transitions, the following parameters are summarized in Tables 2 and 3;

• type (discrete or continuous),
• delay time (discrete) or firing speed (continuous),
• variable of the place inputted to the transition,
• weight and type of the input arc to the transition,
• variable of the place outputted from the transition, and
• weight of the output arc from the transition (for discrete transitions only).

(Note that, since the transition T 59 is the continuous transition defined in Definition 1, the above parameters except variables of places are not described in the corresponding columns of Table 2.) In the following, we will present the strategy for determining the parameters in Tables 2 and 3. Delays and Firing Speeds of Transitions

• Default delay time of discrete transition is “1”.

If necessary, appropriate delay time is assigned to a discrete transition depending on the underlying biological facts.

• Basically, the firing speed of a continuous transition is given by the following formula;

mX a , where mX is the variable of a place incoming to the transition through an arc and the value a is a positive

• integer. The value a represents the reaction speed of reaction and 10 is chosen as the default value of a

in our modeling. In the case of high (low) degradation speed, a smaller (larger) number should be used for the the value a.

• At all continuous places representing concentrations of substances, continuous transitions are attached for

modeling the degradation of substance together with its degradation rate. Basically, the following two types of degradation rates are used; – high speed degradation rate: mX/10, and – low speed degradation rate: mX/10000 Weights of Arcs

• Default weights of the arcs from/to a discrete transition is “1”. If necessary, appropriate weight(s) are

assigned to the arc(s) depending on the underlying biological facts.

• Default weight of the arc to a continuous transition is “0”. If necessary, appropriate weight is assigned to

the arc depending on the underlying biological facts. Note that, from the definition of HFPN [16], it does not need to assign a weight to the arc outgoing from a continuous transition. The delay of discrete transition and/or the firing speed of continuos transition of HFPN and the weights of these transitions are going to be tuned up at the parameter tuning stage in the modeling. The places which have initial values greater than zero are listed in Table 4. 5

SLIDE 6

Figure 3: lac operon: The enzyme β-galactosidase, the product of the lacZ gene, hydrolyzes lactose to glucose and galactose. The lacZ gene is transcribed in a single mRNA along with two other genes, lacY and lacA. lacY encodes the permease that brings lactose into the cell, and lacA encodes an acetylase that is believed to detoxify thioglactosides, which, along with lactose, are transported into cells by lacY. The “operator” lies within the “promoter”, and the “CAP site” lies just upstream of the promoter.

3.1 Transcription Control Switch

Figure 3 shows the structure of lac operon. In the absence of lactose, repressor is bound to the operator. The repressor prevents RNA polymerase from starting RNA synthesis at the promoter. On the other hand, in the presence of lactose and the absence of glucose, the catabolite activator protein (CAP) is bound to the CAP site. Since the CAP helps RNA polymerase to bind the promoter, the transcription of the lac operon begins in this situation. This regulation mechanism can be expressed by the Petri net of Figure 4 consisting of only discrete elements.

m2 CAP site m3

• prator

m1 Promoter T63 T64 T65 T42 T43

Figure 4: HFPN representation of the control mechanism of lac operon transcription switch. Discrete elements are used for representing switching mechanism. The place “promoter” (m1) represents the status of the transcription of lac operon. That is, if this place contains token(s), the lac operon is transcribing, but if this place has no token, transcription of the operon does not begin. The rates of releasing CAP and repressor from the DNA, are assigned to the transitions T 42 and T 43 as delay times, respectively. The production rates of CAP and repressor are also assigned to the transitions T 63 and T 64, respectively, as the delay times. Each time the transition T 65 fires, the place “promoter” gets one token. This transition fires when the both

• f following conditions hold;
• the place “CAP site” (m2) contains token(s) and,

– This represents that the protein CAP is bound to the CAP site.

• the place “operator” (m3) has no token.

– This represents that the repressor is not bound to the operator.

3.2 Positive regulation

The DNA binding of CAP is regulated by glucose to ensure that the transcription of lac operon begins only in the absence of glucose. In fact, glucose starvation induces an increase in the intracellular levels of cyclic AMP (cAMP). Transcription of lac operon is activated with the help of CAP to which cAMP binds. When glucose is plentiful, cAMP levels drop; cAMP therefore dissociates from CAP, which reverts to an inactive form that 6

SLIDE 7

m4 CAP m5 cAMP m11 AMP m12 ADP m6 glucose m2 CAP site m3

• prator

m1 Promoter T45 T63 T64 T65 T42 T43 T19 T20 T21 T1 T17 T80 T79 T81 T82

100 5

glycolytic Pathway (Figure 8)

Figure 5: Positive regulation mechanism: This figure properly contains the HFPN model of Figure 4. can no longer bind DNA. This regulatory mechanism by CAP and cAMP is called positive regulation of the lac

• peron. Figure 5 shows an HFPN model which represents positive regulation of the lac operon.

Continuous places are used for representing the concentrations of the substances CAP, cAMP, AMP, ADP, and glucose. Tokens in the places “CAP” (m4) and “cAMP” (m5) should not be consumed by the firing of the transition T 63, since both of CAP and cAMP are not lost by forming a complex of these two substances. Then, two test arcs are used from the places “CAP” and “cAMP” to the transition T 63. The weight of the arc from the place cAMP to the transition T 63 is set to 100, while the weight of the arc from the place “CAP” to that transition is 1, which was determined by hand-tuning with referring the simulation results. After both of the concentrations of CAP and cAMP exceed the thresholds which are given to these two arcs as weights, the transition T 63 can fire, transferring a token from the transition T 63 to the place “CAP site.” In general, reactions among cAMP, AMP, and ADP are reversible. Then, between the places “cAMP” and “AMP” (the places “AMP” and “ADP”), the transition T 80 (the transition T 82) is described for representing the reversible reaction together with the transition T 79 (the transition T 81). To the places “cAMP”, “AMP”, and “ADP”, 1, 200, and 200 are assigned as initial values, respectively. Recall that when glucose is plentiful, cAMP levels drop. This phenomenon is represented by the inhibitory arc from the place “glucose” to the transition T 80. When the concentration in the place “glucose” exceeds the threshold given at this inhibitory arc, the transition T 80 stops its firing. In this model, since we supposed that the CAP is produced continuously, this production mechanism is modeled with the place “CAP” and the transitions T 45 and T 1. 5 is set to the place “CAP” as an initial value, since CAP is produced by the production mechanism which is independent of the mechanism being described here. Finally, the transitions T 1, T 17, T 19, T 20, and T 21 are used for representing the natural degradation of the corresponding substances.

3.3 Negative regulation

In the presence of lactose, a small sugar molecule called allolactose is formed in the cell. Allolactose binds to the repressor protein, and when it reaches a high enough concentration, transcription is turned on by decreasing the affinity of the repressor protein for the operator site. The repressor protein is the product of the lacI gene which is lied upstream of the lac operon. Actually, after forming a tetramer, the repressor protein can bind the

• perator site.

By adding this negative regulation mechanism to Figure 5, Figure 6 is obtained. Since, in the literature, it is described that repressor should be produced enough prior to other substances productions, parameters relating to this negative regulation are set to faster value than other parameters. Discrete place is used for representing the promoter site of the gene lacI. When this discrete place “repressor 7

SLIDE 8

m4 CAP m13 repressor promoter m14 mRNA repressor m15 repressor m16 4 repressor m17 DNA bind m29 lactose

• utside of a cell

m9 lactose m8 arolactose m10 complex m5 cAMP m11 AMP m12 ADP m6 glucose m18 not bind m7

• prator bind

m2 CAP site m3

• prator

m1 Promoter T45 T57 T58 T59 T61 T62 T63 T64 T65 T42 T43 T4 T5 T6 T19 T20 T21 T1 T2 T3 T60 T13 T14 T17 T15 T18 T80 T79 T81 T82 T78 T46 T77 T76 T94

1 4 1 1 4 100 5 glycolytic Pathway (Figure 8)

Figure 6: Negative regulation mechanism: This figure properly contains the HFPN model of Figure 5. promoter” gets token(s), the transcription of the lacI gene begins. We can determine the transcription frequency by the delay rate of the transition T 46. Transcription and translation mechanisms are modeled by the places “mRNA repressor” (m14) and “repressor” (m15) and the transitions T 57, T 2, T 58, and T 3. The reaction composing monomers (Figure 7 (a)) to a tetramer can be represented by HDN as is shown in Figure 7 (b). (HPN can not model this type of reaction, since, any reaction speed should be realized by assigning a function of values of the places to the continuous transition.) By comparing this representation and the representation of Figure 7 (c), we can recognize that HFPN allows us to represent such reaction naturally and intuitively. This function of HFPN is used in the modeling of tetramer formation constituted by two places “repressor” and “4 repressor” (m16) and three transitions T 3, T 59, and T 60. Actually, the following functions are assigned for the input/output arcs to/from the transition T 59 as the speed of flows;

• for the input arc to T 59 :

4×m15 5

, and

• for the output arc from T 59: m15

5 .

Note that the speed of the input arc is four times as fast as the speed of the output arc. For the repressor forming tetramer, we considered that about 96% of them bind to operator site, about 0.0333% of them bind to other DNA sites, and about 0.0001% of them do not bind to DNA. These percentages are determined from the description in [14]. The places “operator bind” (m7), “DNA bind” (m17), and “not bind” (m18) represent the amount of these repressors. According to the binding rate described above, the firing speeds of the transitions T 60, T 61, and T 62 were given as follows; T 60 : 96 × m16 100 T 61 : 33 × m16 10000 T 62 : m16 10000. We separate the concentration of lactose to two places “lactose outside of a cell” (m29) and “lactose” (m9) for the convenience of describing the function of the lacY gene in the next subsection. Concentration of the 8

SLIDE 9

I1 I2 I2 I1

T01 T01 T01 T01

HDN

Reaction decomposing monomer to tetramer

HFPN v1 v2

v01=v2=

[I1] 5

v1=4*v2

monomer tetramer

(a) (b) (c) Figure 7: HDN and HFPN representations of the reaction composing monomers to a tetramer. The speed of v1 is four times as fast as the speed of v2. v01 is the firing speed of the transition T01. [I1] ([I2]) represents the content of the place I1 (I2). HPN can not model this type of reaction, since, any reaction speed should be realized by assigning a function of values of the places to the continuous transition. arolactose is represented by the place “arolactose” (m8) whose accumulation rates are given at the transitions T 94 and T 76. It is known [14] that arolactose is produced from the lactose existing outside of the cell as well as the lactose inside the cell. Since it is natural to consider that the production speed of arolactose is faster than the passing rate of arolactose through the cell membrane, the speed of transition T 76 (m9/5) is set to be faster than the speed of transition T 94 (m29/10). The negative regulation in Figure 6 is realized in the following way;

• The place “operator” gets token if

– the concentration of the place “operator bind” exceeds the threshold 1 given at the test arc to the transition T 64 as a weight, and – the concentration of the place “arolactose” does not exceed the threshold 4 given at the inhibitory arc to the transition T 64. Four molecules of allolactose need to bind with one molecule of tetramer repressor. Accordingly the function

• f input arc from the place ”allolactose” to the transition T 78 should be four times faster than that of input

arc from the place ”4 repressor”. In summary, we gave formulas,

• 4 × m8 × m16 to the arc from the place ”allolactose” to T 78,
• m8 × m16 to the arc from the place ”4 repressor” to T 78, and
• m8 × m16 to given at the arc from the T 78 to the place ”complex”.

The transition T 78 gives the complex forming late of the arolactose and the tetramer of repressor proteins. The place “complex” (m10) represents the concentration of the complex. Arolactose can also release the tetramer

• f repressor protein from the operator site by forming complex with it. The transition T 77 and the arcs from/to

the transition model this mechanism. Discrete transitions are used for the transitions T 77 and T 78, since only discrete transition is available for the arc from the discrete place operator (m3) according to Definition 1. Small delay time (0.1) is assigned to the transition T 77, since arolactose should bind to repressor immediately after it is produced. The transitions T 4, T 5, T 6, T 13, T 14, T 15 and T 18 are used for representing the natural degradation of the corresponding substances. 9

SLIDE 10

3.4 Hydrolyzing Lactose to Glucose and Galactose

The lac operon consists of three genes, lacZ, lacY, and lacA. The lacZ gene is transcribed in a single mRNA along with two other genes, lacY and lacA. lacY encodes the permease that brings lactose into the cell, and lacA encodes an acetylase that is believed to detoxify thiogalactosides, which, along with lactose, are transported into cells by lacY. The lac operon transcription and translation mechanism is described in Figure 8 together with the effects

• f two products of the genes lacZ and lacY on hydrolyzing lactose to glucose and galactose. The effect of the

gene lacA is not included in this figure, since the lacA gene does not work in both of lac operon gene regulatory mechanism and glycolytic pathway.

1

m4 CAP m13 repressor promoter m14 mRNA repressor m15 repressor m16 4 repressor m17 DNA bind m29 lactose

• utside of a cell

m9 lactose m8 arolactose m10 complex m5 cAMP m11 AMP m12 ADP m30 galactose m6 glucose m18 not bind m7

• prator bind

m2 CAP site m3

• prator

m21 X1 m19 mRNA lacZ m23 mRNA lacY m27 mRNA lacA m20 LacZ m24 LacY LacA m22 X2 m25 X3 m26 X4 m1 Promoter T45 T57 T58 T59 T61 T62 T63 T64 T65 T42 T43 T66 T72 T75 T4 T5 T6 T19 T16 T20 T21 T69 T67 T70 T73 T1 T2 T3 T60 T13 T14 T17 T15 T18 T68 T71 T7 T8 T10 T9 T11 T80 T79 T81 T82 T78 T46 T77 T76 T94 T74 m28 T12

glycolytic Pathway (Figure 8) 1 1 4 100 4 5

Figure 8: HFPN model of lac operon gene regulatory mechanism. This figure properly contains the HFPN model in Figure 6. The places “mRNA lacZ”, “mRNA lacY”, and “mRNA lacA” represent the concentrations of mRNAs transcribed from the genes lacZ, lacY, and lacA, respectively, and the transcription rates are given at the discrete transitions T 66, T 69, and T 71 as delay times. The places “LacZ”, “LacY”, and “LacA” represent the concentrations of proteins translated from the lacZ, 10

SLIDE 11

lacY, and lacA mRNAs, and the translation rates are given at the continuous transitions T 67, T 70, and T 93. As is shown in Table 4, we give the initial values 5, 2.5, 1 of proteins LacZ, LacY, and LacA, respectively, to the corresponding places. These values are chosen according to the production ratios of LacZ, LacY, and LacA proteins [14]. Actually, the formulas m19, m23

2 , and m27 5

are assigned to the transitions T 67, T 70, and T 73 as the speeds, according to the fact that the proteins of lacZ, lacY, and lacA are produced in the ratios 1 : 1

2 : 1 5.

Degradation rates of mRNAs (proteins) are assigned to the transitions T 7, T 9, and T 11 (T8, T 10, and T 12). In this model, for representing the lac operon DNA, the following discrete elements,

• discrete transitions T 66, T 68, T 69, T 71, and T 72, and
• discrete places X1 (m21), X2 (m22), X3 (m25), and X4 (m26),

are used. The discrete places X1 (X3) represents the Boolean status of the transcription of lacZ gene (lacY gene). That is, each time transcription of lacZ (lacY) is finished, the place X1 (X3) gets token. At the discrete transition T 68 (T71), the delay time 0.051 (0.065) required for RNA polymerase moving from the end of lacZ gene to the beginning of lacY gene (the end of lacY gene to the beginning of lacA gene)) is assigned. The delay times 3.075 and 1.254 are assigned to the transitions T 66 and T 69, respectively, according to the fact that the length of lacZ gene (lacY gene) is 3075bp (1254bp). The delay time 0.682 at the transition T 72 represents the length of lacA gene (the length of lacA gene is 682bp). The lengths of the genes are obtained from the website [32]. We can know the transcription status of the gene lacY (lacA) by observing whether the discrete places X2 (X4) contains token(s) or not. Recall that product of the gene lacZ is an enzyme which hydrolyzes lactose to glucose and galactose. This reaction is represented by using the places “lactose”, “galactose” (m30), and glucose (m6), and the transitions T 75, T 16, and T 17. In our modeling, 20 is assigned to each of the places “lactose” and “lactose outside of a cell” as an initial value. Test arc is used from the place lacZ to the transition T 75, since by the reaction of enzyme, no substance is consumed. We considered that the production rates of glucose and lactose depend on both of the concentration of lactose and the concentration of product of lacZ gene. The following formula representing speed of the transition T 74 reflects this idea; m24 × m29 m29 + m24 × 10. We mentioned that the gene lacY encodes the permease that brings lactose into the cell. In Figure 8, this function is realized with the places “LacY” (m29), “lactose outside of a cell” (m29), and “lactose” (m9), and the transitions T 74, T 13, and T 14. Since the product of lacY gene is an enzyme, test arc is used from the place “LacY” to the transition T 74, and the speed of this transition is given by the following formula according to the same idea above; m20 × m9 m9 + m20 × 10. The weight 2.5 (5) of the arc from the place LacY (m24) (LacZ (m20)) to the transition T 77 (T75) corre- sponds to the basal concentration of LacY (LacZ) presented in Table 4.

3.5 Glycolytic Pathway

We model the glycolytic pathway with HFPN that was first modeled with the discrete Petri net by Reddy [26]. In the glycolytic pathway, a glucose is converted into two molecules of pyruvate, where the enzymatic reactions shown in Figure 9 and Table 5 are involved. (a) The process converting a glucose to two glyceraldehyde 3-phospahte ((1)-(5)), where two ATP molecules are invested to produce glucose 6-phosphate and fructose 6-phosphate. (b) The process in which the aldehyde group of each glyceralydehyde 3-phosphate molecule is oxidized to a carboxylic acid ((6)-(7)), where the energy from this reaction is utilized for the synthesis of ATP from ADP and inorganic phosphate. (c) The process converting 3-phosphoglycerate to pyruvate ((8)-(10)), where two ATP molecules that were invested in the process (a) are recovered. 11

SLIDE 12

Figure 9: A part of glycolytic pathway participating glucose (Gluc), glucose 6-phosphate (G6P), fructose 6- phosphate (F6P), fructose 1,6- diphosphate (FBP), dihydroxyacetone phosphate (DHAP), glyceraldehyde 3- phosphate (GAP), 1,3-diphosphoglycerate (1,3-BPG), 3-phosphoglycerate (3PG), 2-phosphoglycerate (2PG), phosphoenolpyruvate (PEP), and pyruvate (Pry). More details about this pathway shall be explained in the following modeling process. Figure 10 shows the HFPN modeling of the glycolytic pathway. Main pathway from glucose to pyruvate acid HFPN modeling starts with creating continuous places corresponding to glucose (m6), intermediates (m31– m39), and pyruvate (m40). Default continuous places are introduced at this step though these places are meant to represent the concentrations of the corresponding substrates. Then, by following the pathway in Figure 9, we put continuous transitions (T83–T 93) together with normal arcs between two consecutive places in the pathway. These transitions and arcs shall represent the reactions, but default transitions and arcs are initially introduced without any parameter tuning. By considering the natural degradation of substrates, we put continuous transitions T 17 for glucose and T 23–T 31 for intermediates with normal arcs. Since these transitions are used for representing degradation and do not produce any, no arcs are going out from these transitions. By taking into account the fact that natural degradation is very slow in glycolysis, each of the firing speed of these transitions is given by the formula mX/10000, where X = 17, 23, 24, . . ., 31. Productions of ATP from ADP and NADH from phosphoric acid Then we consider ADP, ATP, Pi (phosphoric acid) and NADH. In the pathway of Figure 9, two ADP molecules and two Pi’s are invested for producing four ATP molecules and two NADH molecules. Continuous places “ADP” (m12), “ATP” (m51), “Pi” (m52, initial value=200), and “NADH” (m53) are created to represent ADP, ATP, Pi and NADH. We attach continuous transitions (T21 and T 22) representing the natural degradation

• f ADP and ATP. The firing speeds of them are set to be very slow (T21 : m21/10000, T 22 : m22/10000) as

is set for intermediates in the above. For the process of ATP→ADP in the reaction (1), the normal arc from the place “ATP” to the transition T 83 and the normal arc from the transition T 83 to the place “ADP” are

• introduced. In the same way, normal arcs connected to T 85 are introduced for the process of ATP→ADP in

the reaction (3). Similarly, for representing the reactions (7) and (10), the transitions T 90 and T 93 are used,

• respectively. The reaction (6) Pi→NADH is realized with the transition T 89. Since it can be considered that

ATP and ADP are degraded slowly, the formulas m22/10000 and m21/10000 are assigned to the transitions T 22 and T 21 as degradation speed, respectively. Installing enzyme reactions Places having valuables m41, m42,...,m50 represent concentrations of enzymes. Initial values of these valu- ables are set to 5. The production rates of these enzymes are assigned to the transitions T 47,T 48,...,T 56 whose 12

SLIDE 13

m6

glucose

m12

ADP

m51

ATP

m31

G6P

m32

F6P

m33

FBP

m35

GAP

m36

1,3-BPG

m37

3PG

m38

2PG

m39

PEP

m40

pyruvic acid

m52

Pi

m53 NADH m49

enolase

m50

pyruvate kinase

m48

phosphoglycerate mutase

m47

phosphoglycerate kinase

m46 glyceraldehyde-3-phosphate dehydrogenase m45 triosephosphate isomerase m44

aldolase

m43

phosphofructokinase

m42

phosphoglucose isomerase

m41

hexokinase

m34 DHAP T17 T21 T22 T83 T23 T84 T24 T85 T25 T86 T27 T88 T28 T89 T29 T90 T30 T91 T31 T92 T47 T32 T48 T33 T49 T34 T50 T35 T26 T87 T51 T36 T52 T37 T53 T38 T54 T39 T55 T40 T56 T41

lac operon gene regulatory network (Figure 6)

Vmax

Pv

Km

Pk

Figure 10: HFPN model of glycolytic pathway. Michaelis-Menten’s equation is applied at the reactions in the main pathway from Glucose to pyruvic acid. Test arcs are used for representing enzyme reactions. 13

SLIDE 14

speeds are set to 1. The transitions T 32, T 33,...,T 41 have the firing speeds of mX/10s (X = 41, 42, ..., 50), which represent the speeds of natural degradation of these enzymes. Test arcs are used for enzyme reactions, since an enzyme itself is not consumed in the reaction. To each of these test arcs, the value 3 is chosen as a weight of the arc. Reaction speeds in the main pathway For representing the reaction speeds in the main pathway, we adopt Michaelis-Menten equation such as; Vmax[S] Km + [S], where [S] is the concentration of substrate, Vmax is the maximum speed of the reaction, and Km is a Michaelis

• constant. In our model, we let Vmax = 1 and Km = 1
• 2. The two independent places Pv and Pk in Figure 10

represent these valuables Vmax and Km. The values of Vmax and Km can be easily manipulated by changing the contents of the places Pv and Pk, respectively. Michaelis-Menten equation is used for representing each of the firing speeds of the transitions T 83, T 84,..., and T 93. For example, for the transition T 84, the formula Vmaxm31 Km + m31 is used. ADP molecules and two Pi’s are invested for producing four ATP molecules and two NADH molecules. It is known that

• 1. one fructose 1,6-diphosphate (FBP) molecule is invested for producing two glyceraldehyde 3-phosphate

(GAP) molecules, and

• 2. one dihydroxyacetone phosphate (DHAP) is invested for producing two GAP molecules.

These reactions are realized by assigning the following formulas to the arcs inputted to/outputted from the transitions T 86 and 88 as firing speeds;

• 1. for the reaction from FBP to GAP;

the arc from FBP to T 86 = Vmaxm33 Km + m33, the arc from T 86 to GAP = 2 × Vmaxm33 Km + m33 ,

• 2. for the reaction from DHAP to GAP;

the arc from DHAP to T 88 = Vmaxm34 Km + m34, the arc from T 88 to GAP = 2 × Vmaxm34 Km + m34 . This demonstrates the effectiveness of representing biopathways by HFPN.

3.6 Simulation Results of Mutants by Genomic Object Net

Simulations for five mutants were carried out by Genomic Object Net [34] which are listed below together with the methods for realizing these mutants in the HFPN model of Figure 9. lacZ− the mutant which can not produce β-galactosidase,

• delete the transition T 69,

lacY− the mutant which can not produce β-galactoside permease,

• delete the transition T 70,

lacI− the mutant in which 4 repressor monomers can not constitute one active repressor tetramer, 14

SLIDE 15

• delete the transition T 59,

lacIs the mutant to which arolactose can not bind,

• delete the transitions T 77 and T 78 together with the inhibitory arc from the place “arolactose” to

the transition T 64, lacI−d the mutant which can not bind to the DNA,

• delete the transitions T 60 and T 61.

Behavior of the wild type In both of Figures 11 and 12, the concentration behaviors of lactose (outside of a cell), lactose, glucose, LacZ (β-galactosidase), and LacY (β-galactoside permease) are described. From the beginning, glucose is degraded, since it is consumed by glycolytic pathway. At the time 55 the glucose is disappeared, the transcription of lac

• peron begins, producing LacZ protein (time=60), and successively, LacY protein is begun to produce. With

comparing the concentration behavior of lactose and lactose (outside of a cell), we can recognize that the LacY protein works well. At the time 65 the concentration of LacZ exceeds 10, decomposition of lactose to glucose and galactose starts, increasing the concentration of glucose again. Just after the glucose is completely consumed again, the transcription of lac operon is stopped, keeping the concentration of LacZ and LacY proteins at some levels (from the assumption that the degradation speed of these proteins are very show) [14].

100 200 300 t 10 20 30 40 50 100 200 300 t 0,05 0,1 0,15 100 200 300 t 10 20 30 40 50 100 200 300 t 5 10 15 100 200 300 t 10 20 30 40 50 100 200 300 t 10 20 30 100 200 300 t 20 40 60

glucose lactose (inside

• f cell)

Lac Z Lac Y lactose (outside

• f cell)

Wild type lac I - lac I s lac I -d

100 200 300 t 2,433 2,467 2,5 100 200 300 t 4,85 4,9 4,95 5 100 200 300 t 49,7 49,8 49,9 50 100 200 300 t 10 20 30 40 50 100 200 300 t 5 10 15 100 200 300 t 100 200 100 200 300 t 100 200 300 400 100 200 300 t 10 20 30 40 50 100 200 300 t 10 20 30 40 50 100 200 300 t 5 10 15 100 200 300 t 10 20 30 40 50 100 200 300 t 100 200 100 200 300 t 100 200 300 400

Figure 11: Simulation results of lac repressor mutant. Concentrations of proteins LacZ and LacY keep grow- ing, since the mutants lacI− and lacI−d lose abilities to bind at the operator site. In the mutant lacIs, the transcription of lac operon does not begin, since repressor can not remove from the operator site. Behavior of the lac repressor mutant Figures 11 shows the simulation results of the mutants lacI−, lacIs, and lacI−d obtained from Genomic Object

• Net. In the lacI− and lacI−d mutants, LacZ protein and LacY protein are produced, while these proteins are not

produced in the lacIs mutant. Furthermore, in the lacI− and lacI−d mutants, the concentrations of LacZ and LacY proteins keep growing (except the period of glucose re-production), even after stopping the decomposition

• f lactose. Note that these simulation results support the biological experimental results [14].

15

SLIDE 16

lac Z - lac Y -

100 200 300 t 10 20 30 40 50 100 200 300 t 5 10 15 100 200 300 t 10 20 30 40 50 100 200 300 t 10 20 30 100 200 300 t 20 40 60

glucose lactose (inside

• f cell)

Lac Z Lac Y lactose (outside

• f cell)

Wild type

100 200 300 t 10 20 30 40 50 100 200 300 t 49,85 49,9 49,95 50 100 200 300 t 10 20 30 40 50 100 200 300 t

• 1
• 0,5

0,5 1 100 200 300 t 10 20 30 40 50 100 200 300 t

• 1
• 0,5

0,5 1 100 200 300 t 100 200 300 100 200 300 t 200 400 600 100 200 300 t

• 1
• 0,5

0,5 1 100 200 300 t 10 20 30 40 50

Figure 12: Simulation results of lacZ− and lacY− mutants. We can see that the glucose is not produced in the lacZ− mutant and the lactose can not enter inside of a cell in the lacY− mutant. 16

SLIDE 17

Behavior of the lac operon mutant Figures 12 shows the simulation results of the mutants lacZ− and lacY− obtained from Genomic Object

• Net. From this figure, we can observe that, in the lacZ− mutant, once glucose is completely consumed, it is

never produced again. On the other hand, in the lacY− mutant, the concentration of lactose (inside of a cell) never grows, since a lactose can not go through a cell membrane in this mutant. Note that these observations from the simulation results support the experimental results [14]. 17

SLIDE 18

Table 2: Transitions in Figure 8. All transitions in this figure are listed in “Name” column. The symbol D or C in “Type” column represents the type of transition, discrete transition or continuous transition, respectively. In “Delay/Speed” column, the firing speed of continuous transition or the delay time of discrete transition is described depending on the type of transition. The column “From”, which represents incoming arc(s) to a transition, is divided into three sub-columns, “variable” (variable name(s) of the place(s) attached to the incoming arc(s)), “weight” (weight of the incoming arc(s)), and “type” (N, T, I, and F represent normal, test, inhibitory arcs, and functional arcs respectively). The column “To”, which represents outgoing arc(s) from a transition, is divided into two sub-columns, “variable” (variable name(s) of the place(s) attached to the outgoing arc(s) and “weight” (weight of the outgoing arc(s)). Name Type Delay / Speed From To Comment variable weight type variable weight T 1 C m4/10 m4 N — —

degradation rate of CAP

T 2 C m14/10 m14 N — —

degradation rate

• f mRNA repressor

T 3 C m15/10 m15 N — —

degradation rate

• f repressor

T 4 C m17/10 m17 N — —

degradation rate of repressor binding to DNA

T 5 C m18/10 m18 N — —

degradation rate of repressor not binding to DNA

T 6 C m7/10 m7 N — —

degradation rate of repressor binding to operator region

T 7 C m19/10 m19 N — —

degradation rate of lacZ mRNA

T 8 C m20/10000 m20 N — —

degradation rate of LacZ

T 9 C m23/10 m23 N — —

degradation rate of lacY mRNA

T 10 C m24/10 m24 N — —

degradation rate of LacY

T 11 C m27/10 m27 N — —

degradation rate of lacA mRNA

T 12 C m28/10000 m28 N — —

degradation rate of LacA

T 13 C m29/10000 m29 N — —

degradation rate

• f lactose outside of cell

T 14 C m9/10000 m9 N — —

degradation rate of lactose

T 15 C m8/2 m8 N — —

degradation rate of arolactose

T 16 C m30/10000 m30 N — —

degradation rate of galactose

T 17 C m17/10000 m17 N — —

degradation rate of glucose

T 18 C m10/10000 m10 N — —

degradation rate of complex

T 19 C m5/10000 m5 N — —

degradation rate of cAMP

T 20 C m11/10000 m11 N — —

degradation rate of AMP

T 21 C m12/10 m12 N — —

degradation rate of ADP

T 42 D 1 m2 1 N — —

CAP releasing rate

T 43 D 1 m3 1 N — —

repressor releasing rate

T 45 C 1 — — N m4 —

CAP production rate

T 46 D 1 — — — m13 1

activation of repressor gene

T 57 D 1.082 m13 1 N m14 1

transcription rate of repressor

T 58 C m14 m14 N m15 —

translation rate of repressor mRNA

T 59 C — m15 — F† m16 —

conformation rate of repressor

T 60 C 96 ∗ m16/100 m16 N m7 —

repressor binding rate to

• perator

T 61 C 399 ∗ m16/10000 m16 N m17 —

repressor binding rate to the DNA other than repressor site

T 62 C m16/1000 m16 N m18 —

rate of repressor which do not bind any DNA

18

SLIDE 19

Table 3: Transitions in Figure 8 (continued) Name Type Delay / Speed From To Comment variable weight type variable weight T 63 D 1 m4 1 T m2 1 binding rate of CAP m5 100 T to the CAP site T 64 D 1 m7 1 T m3 1 binding rate of repressor m8 4 I the operon T 65 D 1 m2 1 T m1 1 logical operation of the places m3 1 I “CAP site” and “operator” T 66 D 3.075 m1 1 T m19 1 transription rate of lacZ m21 1 T 67 C m19 m19 1 T m20 — translation rate of lacZ T 68 D 0.051 m21 1 N m22 1 moving rate of RNA polymerase T 69 D 1.254 m22 1 N m23 1 transription rate of lacY m25 1 T 70 C m23/2 m23 1 T m24 — translation rate of lacY T 71 D 0.065 m25 1 N m26 1 moving rate of RNA polymerase T 72 D 0.682 m26 1 N m27 1 transription rate of lacA T 73 C m27/5 m27 1 T m28 — translation rate of lacA T 74 C

m24×m29 m29+m24×10

m29 N m9 — transforming rate of lactose into m24 2.5 T a cell T 75 C

m20×m9 m9+m20×10

m9 N m30 — decomposing rate of lactose m20 5 T m6 — to galactose and glucose T 76 C m9/5 m9 1 T m8 — producing rate of arolactose from lactose inside of a cell T 77 D 0.5 m3 1 N m10 1 conformation rate of m8 1 N repressor and arolactose T 78 C m8 × m16 m8 4 N m10 1 conformation rate of m16 1 N repressor and arolactose T 79 C m5/10 m5 N m11 — reaction rate: cAMP to AMP T 80 C m11/10 m11 N m5 — reaction rate: AMP to cAMP m6 5 I T 81 C m11/10 m11 N m12 — reaction rate: AMP to ADP T 82 C m12/10 m12 N m11 — reaction rate: ADP to AMP T 94 C m29/10 m29 T m8 — producing rate of arolactose from lactose outside of a cell † During m15 > 0, the amount flow into the transition T 59 in the speed (4×m15)/5 from the place “repressor” and the amount flow from the transition T 59 to the place “4 repressor” in the speed m15/5. 19

SLIDE 20

Table 4: The places having non zero initial values in Figure 8. Name variable initial value Comment CAP m4 5 concentration of CAP cAMP m5 100 concentration of cAMP glucose m6 50 concentration of glucose AMP m11 200 concentration of AMP ADP m12 200 concentration of ADP LacZ m20 5 concentration of LacX LacY m24 2.5 concentration of LacY LacA m28 1 concentration of LacA lactose outside of a cell m29 50 lactose outside of a cell Table 5: Mapping between numbers in Figure 9 and metabolites in reactions Index Enzyme / Reaction Index Enzyme / Reaction 1 hexokinase 2 phosphoglucose isomerase 3 phosphofructokinase 4 aldolase 5 triosephosphate isomerase 6 glyceraldehyde-3-phosphate dehydrogenase 7 phosphoglycerate kinase 8 phosphoglycerate mutase 9 enolase 10 pyruvate kinase 11 lactate dehydrogenase 20

SLIDE 21

4 Modeling Cell Cycle of Fission Yeast with the HFPN

Figure 13 shows a whole HFPN model of the cell cycle regulation pathway of fission yeast. Four blocks of the figure correspond to Figure 14, Figure 15, Figure 16, and Figure 17. Small part at the central-left side of this figure, which is not involved in any of these four blocks, represents the phases of cell cycle. That is, when the place G1 has token(s), the cell cycle is in G1 phase. The places S, G1, and G2 have similar meaning to the place

• G1. If the places off M and off S have token(s), the cell cycle is not in the phases M and S, respectively. When

the place DNA replication (DNA replication end) gets token(s), DNA replication begins (ends). For expressing the reaction that the protein Cig1 is produced only in the G1/S phase, the test arc is used from the place G1 to the transition representing the production of the protein Cig1 (refer to Figure 17). It is known that the S phage does not begin unless the protein Cdc18 is broken down [23, 20]. The test arc from the place representing the breakdown of Cdc18 to the transition DT15 is used to express this. Firing the transition DT15 means that the cell cycle changes from the G1 phase to the S phase.

T T T Cdc13 Cdc2

T1 v1 v9 v11 v10

v2 v3

v7 v8 v6 d1 d2 T3 Cdc13/Cdc2 Cdc13/Cdc2+2P

Wee1 active_Wee1 d d

d Cdc13/Cdc2+3P d active_MPF lamin lamina d d3 v16 v14 v15 v13 v12 Cdc13_are_broken _by_proteasome v5 v4

active_Myt1 stimulate not_stimulate Myt1 Nim1 d d T T T T T T T d

v17

d d d d d d d d d dT d d d PP2A Pyp3 active_Cdc25 Cdc25 Crk1 Mcs2 Crk1/Mon1 Mon1

v18

Crk1/Mon1/Mcs2 active_CAK CAK MAT1

v19

CAKK

v20 v21 v22 v60 d

DT1 DT2

v23 v24

mitosis_start_prophase

DT4

prometaphase

DT5

metaphase DT6 anaphase

DT7

telophase

DT8

break_mitosis

DT9 Cell division count DT3 DT10

DT11 DT12 DT13 DT14 ketugou bunnkai

DT15 DT16 d d cohesin +ubiquitin cohesin_are_broken _by_proteasome active_Cut1 Cut2_are_broken _by_proteasome cohesin d d d v25 v26 v27 v28 v29 v30 v31 v32 v33 v34 Cut1/Cut2+ubiquitin Cut1/Cut2 signal_Slp1/APC/C T Td d d d Cut1 Cut2 T d Slp1 T d APC/C T d Srw1 d Srw1/APC/C Slp1/APC/C d signal_Srw1/APC/C

2 2 2 1 2 1 1 1

T T v35 v36 d d Cig1 Cig2 d d Cig1/Cdc2 Cig2/Cdc2 v23 v39 v40 v38 v41 d

Rum1_are_broken _by_proteasome

v42 v43 v44 v44 d

Cdc18_are_broken _by_proteasome

T

T4

v37 Rum1 Cdc18 T d d Cdc18+P

SCF_ Pop1/Pop2

d Rum1+P

SCF/Cdc18 +ubiquitin SCF/Cdc18 +ubiquitin

d d d d

G1 DNA_ replication_start

• ff_S

rum1_translation_signal S DNA_ replication_end G2 M

• ff_M

DT16 DT15 DT19 DT17 DT20 DT21 DT18

repli

s

initiator_signal Rad17_Sp/Rfc2,3,4,5 Rfc1,2,3,4,5 pr St DNA damage? DNA_replication_delay DT22 DT23 DT26 DT24 T d Rfc2 T d Rad1 T d Rad24_Sp T T d d d T d Hus1 T d Rad9_Sp T d Rfc1 T d Rfc5 T d Chk1 T d Rad17_sp d d v46 v45 v51 v56 v55 Rad1/Hus1 /Rad9_Sp d v47 v52 v53 v58 v48 T d Rfc3 d active_Chk1+P d d signal_Rad17_Sp /Rfc2,3,4,5 signal_Rfc1,2,3,4,5 d T d Rfc4 d Cdc25+P d Cdc25/Rad24_Sp signal_Rad1/ Hus1/Rad9_Sp T pol_alpha d v57 v59 d DNA_replication? DT27 d signal _pol_alpha G2 STOP? DT25 active _Rad3+P Rad3 v50 v49 v54 d extranuclear Cdc25 /Rad24_Sp Cds1 active_Cds1 DT28 S+1 1 pr2 St2 DT30 1 repli+1

Figure 13: A whole HFPN model of the budding yeast cell cycle regulation

4.1 MPF regulation

Figure 14 shows the HFPN model of MPF regulation mechanism of budding yeast. A biological picture describ- ing this MPF regulation is found in the URL [33]. Protein complex formation (Figure 2 (a)) It is known that the proteins Cdc2 and Cdc13 form a complex [24]. Complex formations can be modeled by constructing an HFPN model described in Figure 2 (a). Based on this, the complex formation is constructed with the places Cdc2, Cyc13, and Cdc2/Cdc13 at the left side in Figure 14. The reaction speed of this complex formation is given by [Cdc2]∗[Cdc13]

20

with referring Table 1 1. The transitions d1, d2, and d3 are used for representing degradation of the proteins and the complex, which are attached to the places Cdc13, Cdc2, and Cdc13/Cdc2,

• respectively. According to Table 1, the degradation speeds are given by the formula

[P] 100, where [P] denotes the

content of a place P at which the transition for degradation is attached. The transitions T1 and T3 are used for productions of the proteins Cdc13 and Cdc2, respectively. The speeds of these transitions are set to 1.

1[P] represents the content of a place P.

21

SLIDE 22

The remaining part of Figure 14 were modeled by same manner as described above. That is, a continuous place of HFPN is placed for each substance, and places are connected by transitions and arcs whose parameters (weights of arcs and speeds of transitions) are determined based on the corresponding biological reactions. At transitions labeled with the symbol T in Figure 14, production rates of the corresponding proteins are assigned. In the following, transitions labeled with the symbols T or d have the same meaning as above. The transitions v1, v17, v18, and v19 of Figure 14 are used for complex formations as shown in Figure 2 (a). Enzymic Reactions (Figure 2 (b)) Transitions v2, v3, v5, v7, v8, v9, v10, v11, v20, v21, v22, v23, v24, and v60 are classified into the following three categories according to the enzymic reaction types as shown in Table 1. (Non-constant concentration) Although proteins MPF, CAK, and Cdc25 works as enzyme, activation level of these proteins are not always constant. Based on this fact, the formulas for non-constant concentration in Table 1 are used for the transitions v8, v9, v10, and v60. (Constant concentration) The formula for constant concentration in Table 1 was used for the transitions v2, v3, v5, v7, v11, v20, v22, and v23, since concentrations of kinases corresponding to these transitions keeps almost constant level. The place Pyp3 represents the protein Pyp3 which back up the protein Cdc25 [18]. (Competitive inhibition) For the transitions v21 and v24, the formulas for competitive inhibition in Table 1 were used from the following reasons. The transition v21 represents the reaction that inhibits the protein Cdc25 from becoming non-active by the protein PP2A which keeps almost constant concentration level throughout the cell cycle [12]. The transition v24 represents the reaction that convert the lamin into the lamina as the active MPF is broken down [1]. Protein degradation (Figure 2 (c)) At transitions labeled with the symbol d in Figure 14, degradation rates of the corresponding proteins are

• assigned. For the variable d, we should choose an integer not greater than any other variables, since the speed
• f protein disintegration should be slower than any other reactions. Accordingly, 100 is used for the variable d

in our model. Protein state switching (Figure 2 (d)) State of some proteins will be switched between active and non-active state depending on the concentration

• f other protein. Each of the following combinations of transitions (v4, v5), (v21, v22), (v23, v24), and (v6, v7,

v8) constitutes the state switching shown in Figure 2 (d). Basically, transition speeds are assigned according to the reaction speeds given in Table 1. However, if some enzymic reaction relates to a state switching of protein, a formula for the enzymic reaction will be used. For example, for the protein Wee1, it can be considered that both activation by the MPF and non-activation by the Nim1 are enzymic reactions. Thus, formulas for enzymic reaction are used for the transitions v7 and v8. Protein transformation (Figure 2 (e)) The state of a protein will be transformed by reactions such as phosphorylation and ubiquitination [1]. For example, a part of HFPN in Figure 14 consisting of the places Cdc13/Cdc2+2P, Cdc13/Cdc2+3P, and active CAK and the transition v9 represents the reaction which transforms the protein complex Cdc13/Cdc2+2P to Cdc13/Cdc2+3P by the phosphorylation with the active CAK. Substance migration (Figure 2 (f)) This reaction is not included in Figure 14, but used in the model of checkpoint mechanism of Figure 15. Refer to the next subsection. Usage of discrete elements The transitions v5 and v7 are enzymic reactions with the protein Num1 and the protein complex active MPF. Although it is known that each of these reactions does not occur without a stimulation, the detail mechanism

• f the stimulation is still unclear. In such case, discrete elements are used. At the top-left part of Figure 14,

two discrete places and two discrete transitions constitute a part which can represents two states stimulate and not stimulate. At each of the transitions DT1 and DT2, the time for state changing is assigned. 22

SLIDE 23

T T T Cdc13 Cdc2

T1 v1 v9 v11 v10

v2 v3

v7 v8 v6 d1 d2 T3

Cdc13/Cdc2 Cdc13/Cdc2+2P

Wee1 active_Wee1 d d

d

Cdc13/Cdc2+3P

d

active_MPF lamin lamina

d d3 v16 v14 v15 v13 v12

Cdc13_are_broken _by_proteasome v5 v4

active_Myt1 stimulate not_stimulate Myt1 Nim1 d d T T T T T T T d

v17

d d d d d d d d d dT d d d PP2A Pyp3 active_Cdc25 Cdc25 Crk1 Mcs2 Crk1/Mon1 Mon1

v18

Crk1/Mon1/Mcs2 active_CAK CAK MAT1

v19

CAKK

v20 v21 v22 v60 d

DT1 DT2

v23 v24

Figure 14: HFPN model for MPF regulation

4.2 Checkpoint mechanism

The concept of checkpoint was defined by Hartwell et al. in 1989 [8] as follows. “The events of the cell cycle

• f most organisms are ordered into dependent pathways in which the initiation of late events is dependent on

the completion of early events. Control mechanisms enforcing dependency in the cell cycle are here called checkpoints.” Several checkpoint systems are known to work in the cell cycle. Among them, we will focus on the system

• f “S/M checkpoint” which coordinates the DNA synthesis and the beginning of chromosome separation. The

steps of S/M checkpoint system are the following;

• 1. the checkpoint mechanism is started by an initiator signal such as ultraviolet rays radiation,
• 2. the sensor senses an abnormal status in a cell triggered by the initiator signal, and
• 3. the transducer mediates the abnormal status to the effector or the target.

In Figure 15, two checkpoint mechanisms of fission yeast “DNA damaging checkpoint” and “DNA replication checkpoint” are described together with the part of MPF activation mechanism of Figure 14. In the following,

• nly the mechanism for the DNA damaging checkpoint is explained, since these two mechanisms have similar

structures. The initiator signal is generated at the beginning of S-phase only when the DNA is damaged. This logic of initiator signal is realized with discrete elements at the top-left corner of Figure 15. The content of the discrete place pr represents the DNA damaging status (1:damaged, 0:not damaged). Note that the formula S + 1 is assigned at the inhibitory arc attached to the place pr. Since the content of the discrete place is always set at 1, the place initiator signal can get token only when the discrete place S in Figure 13 has at least one token, where the place S indicates whether the cell cycle is in S-phase ([S]=1) or not ([S]=0). When the concentration of the protein complex Rad17Sp, Rfc2, -3, -4, -5 (place Rad17 Sp/Rfc2,3,4,5) exceeds some fixed level after detecting the initiator signal (transition v46 and the place signal Rad17 Sp/Rfc2,3,4,5), the signal of DNA damaging is passed to the protein complex Rad1/Hus1/Rad9Sp (the places Rad1/Hus1/Rad9 Sp 23

SLIDE 24 T T d d d PP2A

active _Cdc25

Fig.3

Cdc25 v21 v22

initiator_signal

Rad17_Sp/Rfc2,3,4,5

Rfc1,2,3,4,5 pr St DNA damage? DNA_replication_delay DT22 DT23 DT26 DT24 T d Rfc2 T d Rad1 T d Rad24_Sp T T d d d T d Hus1 T d Rad9_Sp T d Rfc1 T d Rfc5 T d Chk1 T d Rad17_sp d d

v46 v45 v51 v56 v55

Rad1/Hus1 /Rad9_Sp d

v47 v52 v53 v58 v48

T d Rfc3 d active_Chk1+P d d signal_Rad17_Sp /Rfc2,3,4,5 signal_Rfc1,2,3,4,5 d T d Rfc4 d Cdc25+P d Cdc25/Rad24_Sp signal_Rad1/ Hus1/Rad9_Sp T pol_alpha d

v57 v59

d DNA_replication? DT27 d signal _pol_alpha G2 STOP? DT25 active _Rad3+P Rad3

v50 v49 v54

d extranuclear Cdc25 /Rad24_Sp Cds1 active_Cds1 DT28

S+1

1

pr2 St2 DT30

1 repli+1

Figure 15: HFPN model of checkpoint mechanism and signal Rad1/Hus1/ Rad9 Sp, and the transition v47). By transmitting this signal to the protein kinase Rad3 (the place Rad3), this kinase is activated through the phosphorylation (the place active Rad3+P). The activated Rad3 phosphorylates the protein Chk1 (the place Chk1), activating this protein (the place active Chk1+P). Note that the protein complex Rad1/Hus1/Rad9Sp and the proteins Rad3 and Chk1 are transducers. The effector in this pathway is the protein Cdc25, which is inactivated by the activated Chk1 (the place active Chk1+P) through the phosphorylation (the transition v52). The inactivated Cdc25 (the place Cdc25+P) is captured by the mobilizing factor Rad24Sp (the place Rad24 Sp. Since the captured Cdc25 (the place Cdc25/Rad24 Sp) is transported to the outside of nucleus, separating from the MPF2, it lose the ability to dephosphorize the Tyr15 site. Eventually, the cell cycle is stopped at G2-phase. By watching the discrete place DNA damage? (G2 STOP?), we can know the status of the signal propagation of DNA damaging (the status of cell cycle termination at G2-phase).

4.3 HFPN models for the M-phase progress and the CKI effect on MPF

Chromosome segregation and cell division occur in M-phase in the order of dramatic events, prometaphase, metaphase, and anaphase. On the other hand, CDK Inhibitor (CKI) is a protein which works as a brake for the MPF, preventing the MPF being out of control from the cell cycle regulation. Rum1 is one of CKIs in yeast cell cycle, which is started to produce at the anaphase of cell cycle, being completely degraded in the S-phase. That is, accumulation of the Rum1 suppresses the MPF activation during the G1-phase. Although both are the important processes in the cell cycle, we only present the HFPN models as shown in Figure 16 and Figure 17 due to the limitation of the space for this paper. Refer to the webpage [33] for the details of them.

5 Simulations by Genomic Object Net

After describing an HFPN of the biological pathway to be modeled, parameters of transition speed/delay and initial values of places have to be determined based on the biological knowledge and/or the facts described in biological literature. In general, many trial and error processes are required until appropriate parameters for

2This transportation corresponds to the reaction of “substance migration” of Figure 2 (f)

24

SLIDE 25 Cdc2 v11 v10 d2 T3 Cdc13/Cdc2+3P d

active_MPF

Fig.3

lamin

lamina d v16 v14 v15 v13 v12 Cdc13_are_broken _by_proteasome T T dT d d d PP2A Pyp3 active_Cdc25 Cdc25 v21 v22 v60 d v23 v24

mitosis_start_prophase

DT4

prometaphase

DT5

metaphase DT6 anaphase

DT7

telophase

DT8

break_mitosis

DT9 Cell division count DT3 DT10

DT11 DT12 DT13 DT14 ketugou bunnkai DT15 DT16 d d cohesin +ubiquitin cohesin_are_broken _by_proteasome active_Cut1 Cut2_are_broken _by_proteasome cohesin d d d

v25 v26 v27 v28 v29 v30 v31 v32 v33 v34

Cut1/Cut2+ubiquitin Cut1/Cut2 signal_Slp1/APC/C T Td d d d Cut1 Cut2 T d Slp1 T d APC/C T d Srw1 d Srw1/APC/C Slp1/APC/C d signal_Srw1/APC/C

2 2 2 1 2 1 1 1

Figure 16: HFPN model for M-phase progress simulation are determined. Since GON provides the GUI specially designed for biological pathway modeling, we can perform these processes very easily and smoothly. The constructed HFPN model of fission yeast cell cycle was simulated by using GON. With showing the simulation result of a known behavior of DNA-damaged fission yeast, the appropriateness of the constructed HFPN model is verified. In addition, we will present a new biological hypothesis obtained from the simulation with changing expression level of the protein Rum1.

5.1 Verification of the constructed model: observing behaviors of DNA-damaged fission yeast

Figure 18 shows concentration behaviors of the proteins and protein complexes related to MPF activation, where DNA-damage occurs during a cell cycle. That is, as observed in Figure 18 (a), the amount of active Cdc25 is reduced rapidly due to the DNA-damage which was realized by putting a token in the place pr in Figure 15 at the time around 550. The signal triggered by this action propagates to the transition v52 which reduces the amount of the place active Cdc25. The mechanism of this signal propagation was explained in the subsection 4.2 with Figure 15. Figure 18 (b) shows the behaviors of the places Cdc13, active MPF and Cdc13/Cdc2+3P in Figure 14. With comparing Figure 18 (a), we can observe that both the protein Cdc13 and the MPF stop oscillating at the time of the DNA-damage. This observation reflects the biological fact that activation of the protein complex Cdc13/Cdc2 is controlled by the Cdc25 activation [1]. This simulation result supports the appropriateness of the constructed HFPN model.

5.2 Simulations of mutants: producing hypotheses for further biological experi- ments

Figure 19 shows the simulation results of Rum1 mutants. The simulation result of (b) is obtained by removing the arc from the transition T 4 to the place Rum1 in Figure 17. Furthermore, by adding the transition with the arc going into the place Rum1, the simulation results of (c) and (d) are obtained. The values 2.0 and 1.0 are assigned to the transition for the case of weak over expression (c) and strong over expression (d), respectively. Behaviors of the protein Rum1 for the figures (a)-(d) are presented in the website [33]. 25

SLIDE 26

Cdc13

Cdc2

T1 v1 d1 d2 T3

Cdc13/Cdc2

Fig.3

d3

T T v35 v36 d d

Cig1 Cig2

d d

Cig1/Cdc2 Cig2/Cdc2

v23 v39 v40 v38 v41 d

Rum1_are_broken _by_proteasome

v42 v43 v44 v44 d

Cdc18_are_broken _by_proteasome

T T4 v37

Rum1 Cdc18

T d d

Cdc18+P

SCF_ Pop1/Pop2

d

Rum1+P

SCF/Cdc18 +ubiquitin SCF/Cdc18 +ubiquitin

d d d d

Figure 17: HFPN model for MPF activating process and CKI effect It is known from the biological experiments that cell cycle does not stop even when the gene rum1 is deactivated [3]. The simulation result of (b) supports this fact. In addition, it can be observed that the period

• f cell cycle of (b) is shorter than that of (a) of wild type. This observation supports the fact that the protein

Rum1 controls the period of cell cycle by breaking down the complex Cdc2/Cdc13 [17]. Biological experiments in [11] elucidate the phenomenon that over-expressing of the protein Rum1 causes the arrest of cell cycle. The simulation result of (d) supports this biological phenomenon. On the other had, the period of simulation result of (c) is longer than that of (d), suggesting that the period of cell cycle becomes longer if the over-expressing rate is decreased to some level. Note that this suggestion has not been confirmed by any biological experiment yet. The above argument highlights that, with the help of computer simulations, biologists can effectively produce hypotheses which will guide them in performing the further biological experiments. 26

SLIDE 27

• (a)
• (b)

Figure 18: Abnormal behaviors of DNA-damaged fission yeast 27

SLIDE 28

• (a) wild type
• (b) deactivated
• (c) over expressed weakly
• (d) over expressed strongly

Figure 19: Simulation results of Rum1 mutants 28

SLIDE 29

References

[1] Alberts, B., Johnson, A., Lewis, J., Raff, M., Roberts, K., and Walter, P. (2002). Molecular Biology of the Cell, Fourth Edition,” Garland Science, 2002. [2] Alla, H, David, R. (1998). Continuous and hybrid Petri nets, J. Circuits Syst Comput, 8, 159–88. [3] Benito, J., Martin-Castellanos, C., and Moreno, S. (1998). Regulation of the G1 phase of the cell cycle by periodic stabilization and degradation of the p25rum1 CDK inhibitor, EMBO J., vol.17, no.2, 482–497. [4] Drath. R. (1998). Hybrid object nets: an object oriented concept for modelling complex hybrid systems, Proceedings of the Hybrid Dynamical Systems, 3rd International Conference on Automation of Mixed Pro- cesses, ADPM’98, 437–42. [5] Genrich, H., Kuffner, R. and Voss, K. (2001). Executable Petri net models for the analysis of metabolic pathways, International Journal on Software Tools for Technology Transfer, 3(4), 394–404. [6] Ghosh, R. and Tomlin, C. (2001). Lateral inhibition through delta-notch signaling: A piecewise affine hybrid model, Proc. 4th International Workshop on Hybrid Systems:Computation and Control, Lecture Notes in Computer Science 2034, 232–246. [7] Goss, P.J.E. and Peccoud, J. (1998). Quantitative modeling of stochastic systems in molecular biology by using Stochastic Petri nets, Proc. Natl. Acad. Sci. USA 95, 6750–6755. [8] Hartwell, L.H. and Weinert, T.A. (1989). Checkpoints: controls that ensure the order of cell cycle events, Science, vol.246, 629–634, 1989. [9] Hofest¨ adt, R. (1994). A Petri net application to model metabolic processes, SAMS 16, 113–122. [10] Hofest¨ adt, R and Thelen, S. (1998). Quantitative modeling of biochemical networks, In Silico Biology 1, 39–53. [11] Jallepalli, P.V. and Kelly, T.J. (1996). Rum1 and Cdc18 link inhibition of cyclin-dependent kinase to the initiation of DNA replication in Schizosaccharomyces pombe, Genes Dev., vol.10, no.5, 541–552. [12] Kinoshita, N., Ohkura, H., Yanagida, M. (1990). Distinct essential roles of type1 and 2A protein phos- phatases in the control of the fission yeast cell division cycle, Cell, vol.63, no.2, 405–415. [13] Kuipers, B.J., Shults, B. (1994). Reasoning in logic about continuous systems, Principles of Knowledge Rep- resentation and Reasoning: Proceedings of the Fourth International Conference KR94, Morgan Kaufmann 391–402. [14] Lewin, B. (1997) Genes VI, Oxford University Press and Cell Press. [15] Matsuno, H., Doi, A., Nagasaki, M. and Miyano, S. (2000). Hybrid Petri net representation of gene regu- latory network. Pacific Symposium on Biocomputing 2000, 338–349. [16] Matsuno, H., Tanaka, Y., Aoshima, H., Doi, A., Matsui, M., Miyano, S. (2003). Biopathways Representation and Simulation on Hybrid Functional Petri Net, In Silico Biology, 3(3), 389–404. [17] Matsuoka,K., Kiyokawa, N., Taguchi, T., Matsui, J., Suzuki, T., Mimori, K., Nakajima, H., Takenouchi, H., Weiran, T., Katagiri, Y.U., and Fujimoto, J. (2002). Rum1, an inhibitor of cyclin-dependent kinase in fission yeast, is negatively regulated by mitogen-activated protein kinase-mediated phosphorylation at Ser and Thr residues, Eur. J. Biochem,, vol.269, 3511–3521. [18] Millar, J.B., Lenaers, G., and Russell, P. (1992). Pyp3 PTPase acts as a mitotic inducer in fission yeast, EMBO J., vol.11, no.13, 4933–4941. [19] Mendes, P. (1993). GEPASI: a software package for modelling the dynamics, steady states and control of biochemical and other systems, Comput, Appl. Biosci. 9, 563–571. [20] Muzi-Falconi, M., Brown, G.W., and Kelly, T.J. (1996). cdc18+ regulates initiation of DNA replication in Schizosaccharomyces pombe, Proc. Natl. Acad. Sci. U S.A., vol.93, issue.4, 1566–1570. 29

SLIDE 30

[21] Nagasaki, M., Onami, S., Miyano, S., Kitano, H. (1999). Bio-calculus: Its Concept and Molecular Interac- tion, Genome Informatics 1999, 10, 133–143. [22] Nagasaki, M., Doi, A., Matsuno, H., Miyano, S. (2004). A Versatile Petri Net Based Architecture for Modeling and Simulation of Complex Biological Processes, Genome Informatics, 15(1), in press. [23] Nishitani, H. and Nurse, P. (1995). p65cdc18 plays a major role controlling the initiation of DNA replication in fission yeast, Cell, vol.83, no.3, 397–405. [24] Nurse, P. (1994) “Ordering S phase and M phase in the cell cycle,” Cell, vol.79, no.4, 547–550. [25] Petri, C.A. (1962) “Kommunikation mit Automaten,” Bonn: Institut fur Instrumentelle Mathematik, Schriften des IIM Nr.3, also, English translation, “Communication with Automata” New York : Griff- iss Air Force Base, Tech. Rep. RADC-TR-65-377, vol.1, Suppl. 1, 1966. [26] Reddy, V.N., Mavrovouniotis, M.L. and Liebman, M.N. (1993). Petri net representations in metabolic pathways, Proc. First ISMB, 328–336. [27] Srivastava, R., Peterson, M.S., and Bentley, W.E. (2001). Stochastic kinetic analysis of the Escherichia coli stress circuit using σ32-targeted antisense, Biotechnol. Bioeng. 75(1), 120–129. [28] Srivastava, R., You, L., Summers, J., and Yin, J. (2002). Stochastic versus deterministic modeling of intracellular viral kinetics, J. Theor. Biol. 218, 309–321. [29] Tomita, M., Hashimoto, K., Takahashi, K., Shimizu, T.S., Matsuzaki, Y., Miyoshi, F., Saito, K., Tanida, S., Yugi, K., Venter, J.C., Hutchison, C.A.3rd. (1999). E-CELL: software environment for whole-cell simulation, Bioinformatics, 15, 72–84. [30] Voss, K., Heiner, M., and Koch, I. (2003). Steady state analysis of metabolic pathways using Petri nets, In Silico Biology, 3(3), 367–387. [31] Watson, J.D., Hopkins, N.H., Roberts, J.W., Steitz, J.A. and Weiner A.M. (1987) Molecular biology of the gene, fourth edition. The Benjamin/Cummings Publishing Company Inc. [32] Escherichia coli GenoList browser, http : //genolist.pasteur.fr/Colibri/genome.cgi [33] http : //genome.ib.sci.yamaguchi-u.ac.jp/˜gon/yeast/ [34] http : //www.GenomicObject.Net/ 30

SLIDE 31

Acknowledgments

The authors would like to thank Dr. Masao Nagasaki, Mr. Atsushi Doi, Ms. Mika Matsui,

• Ms. Sachie Fujita, and Ms. Rie Yamane who have been working together with us in Genomic

Object Net project.

31