Petri net based modelling and simulation of p16-Cdk4/6-Rb pathway

Abstract

Tumor suppressor gene p16 is of utmost interest in investigation of signal transduction pathways due to its gatekeeper role at the G1/S checkpoint of the cell cycle. Defects in p16 result in uncontrolled cell division which leads to progression of malignancy in an organism. In the present research we focus on p16-Cdk4/6-Rb pathway which is a cornerstone of G1 phase of the cell cycle. We implement Pet net formalism and Cell Illustrator software tool to create model of p16-Cdk4/6-Rb pathway and perform a series of simulations to validate the model.

Keywords

• Replicative senescence
• Cell cycle
• p16-Cdk4/6-Rb pathway
• Hybrid functional Petri net
• Cell Illustrator

-

1.1

Introduction

Biological context

Cell division is a fundamental biological process that is essential to continuity of all living organisms. Cell replication or growth is controlled by a complex network of signals, that control the cell cycle. During the cell cycle cells grow to twice their size, copy their chromosomes, and divide into two new cells. The cell cycle is composed of four distinct phases: G1-phase (gap 1), S-phase (synthesis), G2-phase (gap 2) and M-phase (mythosis) [ 14 ]. Cell cycle checkpoints are used between neighboring phases to monitor and regulate the progress of the cell cycle. A cell cannot proceed to the next phase until otherwise checkpoint requirements have been met.

Tumor suppressor gene p16 plays important role in regulating cell grows and division at checkpoint G1/S [ 34 ]. The p16 gene is major tumor suppressor gene that is responsible for replicative senescence. Cell division is not an in nitely continuous process as cells undergo a nite number of cumulative population doublings [ 17 ]. Most human normal cells permanently stop dividing after a 50-75 cell divisions and enter a state termed cellular or replicative senescence [ 17 ]. Most tumors contain cells that appear to have bypassed this limit and evaded replicative senescence. Immortality, or even an extended replicative lifespan, greatly increases susceptibility to malignant progression because it permits the extensive cell divisions needed to acquire successive mutations. Thus, cellular senescence may act as a barrier to cancer and play an important role in tumor suppression [ 8 ]. Inactivation of tumor suppressor gene p16, which in fact keeps track of replicative senescence, results in uncontrolled cell division, which leads to cancer [ 4 ].

During G1 phase, proteins Cdk4 and Cdk6 form complex with protein CycD, which in turn phosphorylates the Rb protein family. When Rb is phosphorylated by Cdk4/6 it loses its function and releases its target, the E2F family transcription factors, resulting in the initiation of DNA replication [ 31, 32 ]. Otherwise Rb inhibits transcription factor E2F [ 36 ]. E2F is a transcription factor which initiates transcription of genes required for S phase [ 5 ]. In the case of malignant progression action of p16 inhibits binding of Cdk4/6 with CycD which leaves Rb, and other Rb related proteins [ 25, 35 ]. The p16 targets Cdk4 and Cdk6, rather than the CycD, and actually competes with CycD for Cdk binding. Binding of p16 results in changes in conformation of Cdk proteins so that they can no longer bind CycD [ 29 ]. The p16 may also deactivate preassembled Cdk4/6 CycD complex blocking their function [ 29 ].

The proteins and their complexes are involved in natural degradation. In addition, the CycD protein is also tightly regulated by ubiquitin-dependent degradation [ 2, 13, 23 ]. 1.2

Related work

Over the past two decades considerable e orts have been directed towards Petri net based investigation of biological systems. A series of biological phenomena modelled and simulated in terms of Hybrid Functional Petri Net (HFPN) include molecular interactions in the ower developmental network of Arabidopsis thaliana [ 19 ], lac operon gene regulatory mechanism in the glycolytic pathway of Escherichia coli [ 9 ], cell fate speci cation during Caenorhabditis elegans vulval development [ 21 ], antifolate inhibition of folate metabolism [ 3 ], validation of transcriptional activity of the p53 [ 12 ], glycolytic pathway controlled by the lac operon gene [ 11 ], apoptosis signalling pathway [ 26 ], circadian rhythms of Drosophila [ 26 ], switching mechanism of phage [ 26 ].

In [ 16 ] the authors proposed a hybrid Petri net model of cell cycle. The model comprises both stochastic and deterministic approaches. In this model, stochasticity is used to capture change of the cell size and e ect of noises. This model is centered upon interactions between complexes CycB-Cdk1, Cdh1-APC, and monomers Cdc14 and Cdc20 [ 33 ]. The study expands macro-level understanding of cell cycle control. However, this study does not provide any insights into understanding quantitative behavior of biological components involved in the cell cycle regulation. Indeed, cell cycle regulation is a complex biological mechanism that consists of hundreds of biological components, processes and pathways. It is hard if not impossible to perform quantitative analysis of cell cycle regulation based on modest size model.

Contributions

The present research exploits HFPN to create a model of p16-Cdk4/6-Rb pathway, which is a cornerstone of cell cycle regulation at G1/S checkpoint. We combine biological facts described in Subsection 1.1 and quantitative knowledge on reaction rates provided in [ 11, 12 ] in a HFPN model. Then we use Cell Illustrator software to perform simulation-based model checking to validate the HFPN model. Simulation-based model checking in general provides interesting biological insights which could be used for future wet-lab experiments [ 22 ]. Once the model validated it can be used for obtaining broader understanding of cell cycle regulation.

The manuscript is organized as follows. Section 2 provides a succinct background on HFPN. In Section 3 we develop a HFPN model of p16-Cdk4/6-Rb pathway, and explain relationship between HFPN objects and their biological counterparts. Section 4 presents and analyzes the simulation results. Finally, conclusions are outlined in Section 5. 2

Hybrid Functional Petri Net Biological systems are characterized by interaction of di erent structured processes. A continuous process is used to represent a biological reaction, at which a real number called the reaction speed or reaction rate is assigned as a parameter. Concentration change of the biological components or substrates after the biological reaction is completed is also represented as a real number. Promotion/inhibition mechanisms and checking for presence of this or that biological component or phenomenon are typical discrete processes. Change of quantity in a discrete process is usually expressed by integers or Boolean values.

When modelling biological pathways it is desirable to use a modelling framework that combines both continuous and discrete processes. Related software tools are consequently expected to comprise di erent structured data types including real numbers, integers, Boolean, etc. HFPN [ 26, 21 ] was originally proposed for modelling and simulating biological systems employing hybrid structure and dedicated software Cell Illustrator [ 11, 27 ] provides suitable platform for visualization and simulation of HFPN models.

While modelling with HFPN, the researchers prefer to use terminology that is slightly di erent than the traditional one [ 28 ]. In order to ensure compliance with the biological content Petri net objects such as place, transition, arc and token are respectively renamed as entity, process, connector and quantity. To increase the readability of the paper below we provide a brief description of HFPN model elements. For more detailed information on this issue the readers are referred to [ 10 ].

In context of HFPN an entity is an abstract object that represents biological component or substrate such as DNA, mRNA, protein, enzyme, complex of proteins, etc. Each entity is assigned a numeric value called quantity, which stands for concentration of related substrate. Variables are used to carry concentration values. A process is another abstract object that is used to model biological reaction or phenomenon like transcription, translation, binding, nuclear export/import, ubiquitination and natural degradation. A process de nes the change rate of entity value and establishes interactions among entities. Rate of change is expressed as a formula. The entities and processes are classi ed as being discrete, continuous and generic. A discrete entity is quanti ed by integers. A discrete process causes integervalued change of a quantity. A continuous entity is quanti ed by real numbers, and consequently continuous process causes change of a quantity according to reaction rate formula, which is also represented by real numbers. A generic entity contains structured data type composed of di erent structured data such as Boolean, double and integer. A generic process handles structured data assigned to associated entities. In HFPN we distinguish between process connector, inhibitory connector and association connector. A process connector is adjacent from input entity to a process or from process to its output entity. Weight parameter is used to specify an activation threshold. Process connectors ensure ow of tokens in the model. An inhibitory connector is used to inhibit a process. Inhibitory connectors are integral elements of biological models with competing processes. An association connector establishes adjacency relation between speci ed entity and process under circumstance that occurrence of related process does not cause concentration change. An association connector is often used in modelling of enzymatic and catalytic reactions. 9 5 c 4 1 m 7 c2 ld o h s e r h t 9 1 m 2 3 0 c m P _ B R m * 4 1 m 5 1

T 1 . 0 * 6 1 1 3 c 0 8 1 m N g e d n o i t a l y r o h p s o h P n io r s e i t iv n d u l o l c e C d e t a l y r ho ien p t s o o r h p p P _ B R F2 lxe E p _ P m D co D n il c y c x _6 lep K D m _C co 4 K D C F 2 x E e _ l P p D m _ o B c R m B R Connector c4 c25 In this section we provide step-by-step explanation on how HFPN model of p16-Cdk4/6-Rb pathway is created according to the biological facts provided in Subsection 1.1, and describe relationship between HFPN objects and their biological counterparts.

The entities used in the model are detailed in Table 1. The entities represent mRNAs, nuclear and cytoplasmic proteins, protein complexes, phosphate, ubiquitin, mutation and G1-dysfunction. A variable associated with a continuous entity quanti es concentration of speci ed substrate. To ensure continual phosphorylation of Rb we assume that there exist su cient amount of phosphate and Rb DP E2F concentration. This is why variables m15 and m16 are initially set to 100. Likewise, m25 is set to 100 to guarantee continual ubiquitination of CycD. The initial concentration of mRNAs and consequently protein monomers and their complexes are set to 0 since simulation starts with transcription of related mRNAs. The entities G1-dysfunction and mutation are used to indicate boolean status or presence/absence of corresponding events. The entity nr div counts the number of cell divisions.

The processes used in the present research include transcription, translation, nuclear import/export, binding, ubiquitination, phosphorylation and natural degradation. Relationship between processes and biological phenomena are illustrated in Table 2 and Table 3. It was reported that mutations in the p16 binding site result in diminished capability of p16 binding to Cdk4/6. This particularly leads to loss of function of p16 as an inhibitor of Cdk4/6-CycD complex. In this model, boolean status of mutation is controlled by T16 and m20. Assignment m20==1 constitutes presence of mutation, consequently leading to occurrence of T16 which in deed arrests p16 in cytoplasm. Otherwise T3 occurs generating nuclear import of p16. Likewise, the presence/absence of dysfunction in the G1 phase is controlled by entity G1-dysfunction and variable m22. Assignment m22==1 indicates the presence of dysfunction in the G1 phase. Next p16 acts as inhibitor of Cdk4/6-CycD complex. We use two Boolean variables with total of four distinct combinations. The rules set for associated connectors and processes depend on four distinct combinations of two Boolean variables m20, which represents the mutation in p16, and m22, which stands for the dysfunction in G1 phase. Occurrence of transitions T3, T14, T16, and T17 respectively depend on the rules on connectors c4, c25, c34, and c36. For instance, T3, nuclear import of p16, occurs if there is no mutation in p16. That is, T3 can re only if m20==0. All connectors together with their ring styles, ring scripts, and connector types are described in Table 4. A snapshot of HFPN model is illustrated in Fig. 1.

A net fragment bound to T15 is shown in Fig. 2. This fragment reveals the structural basis for phosphorylation of Rb. Other than connector rules, the phosphorylation of Rb (T15) has its activity rule as: (m20==0 && m22==0 jj (m20==1 && m22==0) jj (m20==1 && m22 ==1). Here, the rst statement part is for the case when p16 is not mutated, and there is no dysfunction in the G1 phase. It is known that replicative senescence should occur if there is no mutation and dysfunction in a cell, which means that the cell stops dividing after 50 divisions [ 17 ]. In our model, the m17 is de ned as a counter which keeps track the number of divisions, and in the case of no mutation and no dysfunction, it is checked whether the counter is less than 50 or not. If it is not, the cell should stop dividing, which means that RB should not be phosphorylated after 50 divisions. The other two statement parts in the activity rule of T15 are the cases when p16 is mutated. If p16 is mutated, then the replicative senescence will not occur and the cell will divide continually leading to progression of malignancy.

Process rates are chosen in accordance with the reaction speeds for speci c reaction types adopted in [ 11, 12 ]. Process rate for transcription is set to 1 to ensure continual mRNA production. The process rates are listed in Table 2. 4

Simulations and Results In this research, simulations were carried out using Cell Illustrator 5.0 (professional version) that is licensed to Eastern Mediterranean University. Simulation results for concentration behaviour of nuclear and cytoplasmic proteins and their complexes are illustrated in Fig. 3-5. We performed simulations for the following four cases: 1. The p16 is not mutated and there is no dysfunction in the G1 phase (m20==0 && m22==0). 2. The p16 is not mutated and there is dysfunction in the G1 phase (m20==0 && m22==1). 3. The p16 is mutated and there is no dysfunction in the G1 phase (m20==1 && m22==0). 4. The p16 is mutated and there is dysfunction in the G1 phase (m20==1 && m22==1).

It is generally assumed that p16 is transported to the nucleus and acts as a CKI to regulate the G1/S cell cycle checkpoint. This phenomenon has been reported in normal cells where the protein was mainly found in the nucleus but not in the cytoplasm [ 5 ]. This fact is supported by the simulation results that are illustrated in Fig.3. For all four cases the concentration of p16(C) is at level 17.5 after almost 50 pt (Petri net time), at which the steady state starts. On the other hand, if there is no dysfunction in the G1 phase (m22==0) and if p16 is not mutated (m20==0) p16(N) is at level 175, that is, almost 10 times more than that of in cytoplasm. It should be noticed that small oscillations in the p16(C) graphs are result of natural degradation which is 10 times slower than the translation process. Mutation in p16 arrests it in cytoplasm. This is why when p16 is mutated its concentration in nucleus is constantly 0.

Healthy and functioning p16 protein forms a complex with Cdk4/6 if it detects a dysfunction. Simulation results, that are illustrated in Fig. 3 and Fig. 4, have shown that p16 Cdk4/6 concentration in cytoplasm and nucleus are respectively at level 125 and 15, i.e. p16 Cdk4/6 concentration in cytoplasm is almost 8 times more than that in nucleus, indicating that p16 Cdk4/6 is accumulated in cytoplasm rather than in nucleus. We were not able to nd an experimental result to compare this nding with. The reasonable explanation for this fact however could be the di erence between reaction rates of nuclear export and binding, i.e., the former is 10 times faster than the latter.

It was reported in [ 24 ] that levels of Cdk proteins in cells vary little throughout the cell cycle. Simulation results for change of Cdk4 and Cdk6 concentrations in nucleus and cytoplasm are shown in Fig. 4-5. These results fully agree with this fact, in sense that concentration of Cdk4 and Cdk6 in nucleus and cytoplasm are respectively at the level 12 and 17 throughout the simulations. This fact remains true even for Cdk4/6 (Fig. 5). 5

Conclusions

The present research explores interaction between HFPN and biological processes, to the bene t of both elds. On the one hand we adopt HFPN for modelling and simulation of speci c biological pathways, and consequently expand the list of HFPN applications. On the other hand, through modelling and simulating with HFPN we obtain broader understanding of cell cycle regulation.

The fact that in normal cells p16 protein is mainly accumulated in the nucleus but not in the cytoplasm [ 5 ] is con rmed by simulation results. The simulation results have shown that the p16 CDK4/6 protein complex is accumulated in cytoplasm rather than in nucleus. We were not able to nd an experimental result to compare this nding with. The simulation results are in agreement with the fact that levels of Cdk proteins in cells vary little throughout the cell cycle [ 24 ]. . 6 / 4 K D C 6 1 p d n a ) N ( 6 1 p , ) C ( 6 1 p r o f s t l u s e r n o i t a l u m i S g i 1 = = 2 2 1 = = 0 2 0 = = 2 2 1 = = 0 2 1 = = 2 2 0 = = 0 2 0 = = 2 2 0 = = 0 2 . ) N ( 4 K D C d n a ) C ( 4 K D C r o f s t l u s e r n o i t a l u m i S g i 1 = = 2 2 1 = = 0 2 0 = = 2 2 1 = = 0 2 1 = = 2 2 0 = = 0 2 0 = = 2 2 0 = = 0 2 , ) C ( 6 K D C r o f s t l u s e r n o i t a l u m i

References

1. Agherbi , H. , Gaussmann-Wenger , A. , Verthuy , C. , Chasson , L. , Serrano , M. , Djabali , M.: Polycomb mediated epigenetic silencing and replication timing at the INK4a/ARF locus during senescence . PLoS One 4 ( 5 ), e5622 ( 2009 ).

2. Alao , J.P.: The regulation of cyclin D1 degradation: roles in cancer development and the potential for therapeutic invention . Molecular Cancer 6 ( 24 ) ( 2007 ).

3. Assaraf , Y. G. , Ifergana , I. , Kadryb , W. N. , PinteRb R . Y.: Computer modelling of antifolate inhibition of folate metabolism using hybrid functional petri nets . Journal of Theoretical Biology 240 , 637 - 647 ( 2006 ).

4. Baker , D.J. , Perez-Terzic , C. , Jin , F. , Pitel , K. , Niederlander , N.J. , Jeganathan , K. , Yamada , S. , Reyes , S. , Rowe , L. , Hiddinga , H.J. , Eberhardt , N.L. , Terzic , A., van Deursen, J.M.: Opposing roles for p16Ink4a and p19Arf in senescence and ageing caused by BubR1 insu ciency . Nat Cell Biol 10 ( 7 ), 825 - 836 ( 2008 ).

5. Bartkova , J. , Lukas , J. , Guldberg , P. , Alsner , J. , Kirkin , A.F. , Zeuthen , J. , Bartek , J.: The p16-cyclin D/Cdk4-pRb pathway as a functional unit frequently altered in melanoma pathogenesis . Cancer Res , 56 ( 23 ), 5475 - 5483 ( 1996 ).

6. Batt , G. , Ropers , D., de Jong, H., Geiselmann , J. , Mateescu , R. , Page , M. , Schneider , D. : Validation of qualitative models of genetic regulatory networks by model checking: analysis of the nutritional stress response in Escherichia coli . Bioinformatics 21 ( Suppl 1 ) 19 - 28 ( 2005 ).

7. Calder , M. , Vyshemirsky , V. , Gilbert , D. , Orton , R. : Analysis of signalling pathways using the prism model checker . In: Proceeding of Computational Methods in Systems Biology , 3 -5 April, Edinburgh, pp. 179 - 190 ( 2005 ).

8. Campisi , J.: Cellular senescence as a tumorsuppressor mechanism . Trends Cell Biol 10 , S27 - S31 ( 2001 ).

9. Castellini , A. , Franco , G. , Manca , V. : Hybrid functional Petri nets as MP systems . Nat Comput 9 ( 1 ), 61 - 81 ( 2010 ).

10. Cell IllustratorTM: User Guide . The University of Tokyo, 2002 -2010 Human Genome Center , Institute of Medical Science.

11. Doi , A. , Fujita , S. , Matsuno , H. , Nagasaki , M. , Miyano , S. : Constructing biological pathway models with Hybrid Functional Petri Nets . In Silico Biol 4 ( 3 ), 271 - 291 ( 2004 ).

12. Doi , A. , Nagasaki , M. , Matsuno , H.,, Miyano , S. : Simulation-based validation of the p53 transcriptional activity with hybrid functional Petri net . In Silico Biology 6 ( 1-2 ), 1 - 13 ( 2006 ).

13. Germain , D. , Russell , A. , Thompson , A. , Hendley , J.: Ubiquitination of free cyclin D1 is independent of phosphorylation on threonine 286 . J Biol Chem 275 ( 16 ), 12074 - 12079 ( 2000 ).

14. Enoch , T. , Nurse , P. : Coupling M phase and S phase: controls maintaining the dependence of mitosis on chromosome replication . Cell 65 ( 6 ), 921 - 923 ( 1991 ).

15. Fisher , J., Piterman , N. , Hajnal , A. , Henzinger , T.A. : Predictive modeling of signaling crosstalk during C. elegans vulval development . PLoS Comput Biol 3 ( 5 ), e92 ( 2007 ).

16. Herajy , M. , Schwarick , M.: A hybrid Petri net model of the eukaryotic cell cycle . In: Proc. 3th International Workshop on Biological Processes and Petri Nets , 29 - 43 , 2012 .

17. Hay ick , L.: The Limited in Vitro Lifetime of Human Diploid Cell Strains . Exp Cell Res 37 , 614 - 636 ( 1965 ).

18. Heath , J. , Kwiatkowska , M. , Norman , G. , Parker , D. , Tymchyshyn , O. : Probabilistic model checking of complex biological pathways . Theor Comput Sci 391 ( 3 ), 239 - 257 ( 2008 ).

19. Kaufmann , K. , Nagasaki , M. , Jauregui , R.: Modelling the molecular interactions in the ower developmental network of Arabidopsis thaliana . Annals of Botany 107 ( 9 ), 1545 - 1556 ( 2011 ).

20. Kwiatkowska , M.Z. , Norman , G. , Parker , D. : Using probabilistic model checking in systems biology . SIGMETRICS Performance Evaluation Review , 35 ( 4 ), 14 - 21 ( 2008 ).

21. Li , C. , Nagasaki , M. , Ueno , K. , Miyano , S. : Simulation-based model checking approach to cell fate speci cation during Caenorhabditis elegans vulval development by hybrid functional Petri net with extension . BMC Systems Biology 3 ( 42 ) ( 2009 ).

22. Li , C. , Nagasaki M ~., Koh , C.H. , Miyano , S. : Online model checking approach based parameter estimation to a neuronal fate decision simulation model in Caenorhabditis elegans with hybrid functional Petri net with extension . Mol Biosyst . 7 ( 5 ), 1576 - 1592 , ( 2011 ).

23. Lin , D.I. , Barbash , O. , Kumar , K.G. , Weber , J.D. , Harper , J.W. , KleinSzanto , A.J. , Rustgi , A. , Fuchs , S.Y. , Diehl , J.A. : Phosphorylation-dependent ubiquitination of cyclin D1 by the SCF(FBX4-alphaB crystallin) complex . Mol Cell 24 ( 3 ), 355 - 366 ( 2006 ).

24. Malumbres , M. , BaRbacid , M.: Cell cycle, CDKs and cancer: a changing paradigm . Nat Rev Cancer , 9 ( 3 ), 153 - 166 ( 2009 ).

25. Matheu , A. , Maraver , A. , Collado , M. , Garcia-Cao , I. , Canamero , M. , Borras , C. , Flores , J.M. , Klatt , P. , Vina , J. , Serrano , M. , Anti-aging activity of the Ink4/Arf locus . Aging Cell 8 ( 2 ), 152 - 161 ( 2009 ).

26. Matsuno , H,, Tanaka , Y. , Aoshima , H. , Doi , A. , Matsui , M. , Miyano , S. : Biopathways representation and simulation on Hybrid Functional Petri Nets . In Silico Biol 3 ( 3 ), 389 - 404 ( 2003 ).

27. Nagasaki , M. , Doi , A. , Matsuno , H. , Miyano , S. : Genomic object net: I. A platform for modelling and simulating biopathways . Appl Bioinformatics 2 ( 3 ), 181 - 184 ( 2004 ).

28. Reisig , W.: Petri Nets: an introduction . EATCS, Monographs on Theoretical Computer Science . Springer, New York, 1985 .

29. Russo A.A. , Tong , L. , Lee , J.O. , Je rey , P.D., Pavletich , N.P. : Structural basis for inhibition of the cyclin-dependent kinase Cdk6 by the tumour suppressor p16INK4a . Nature 395 ( 6699 ), 237 - 243 ( 1998 ).

30. Sherr C.J.: Cancer cell cycle . Science 274 , 1672 - 1677 ( 1996 ).

31. Stevaux O. , Dyson N.J.: A revised picture of the E2F transcriptional network and RB function . Curr Opin Cell Biol 14 , 684 - 691 ( 2002 ).

32. Trimarchi J.M. , Lees J.A. : Sibling rivalry in the E2F family . Nat Rev Mol Cell Biol 3 , 11 - 20 ( 2002 ).

33. Tyson , J. , Novak , B. : A Systems Biology View of the Cell Cycle Control Mechanisms . Elsevier, San Diego, CA, ( 2011 ).

34. Vidal A. , Ko A. : Cell-cycle inhibitors: three families united by a common cause . Gene 247 ( 1-2 ), 1 - 15 ( 2000 ).

35. Walkley , C.R. , Orkin , S.H. : RB is dispensable for self-renewal and multilineage di erentiation of adult hematopoietic stem cells . Proc Natl Acad Sci USA , 103 ( 24 ), 9057 - 9062 ( 2006 ).

36. Weinberg , R.A. , The retinoblastoma protein and cell cycle control . Cell 81 ( 3 ), 323 - 330 ( 1995 ).