Bifurcation in Cell Cycle Dynamics Regulated by p53.

We study the regulating mechanism of p53 on the properties of cell cycle dynamics in the light of the proposed model of interacting p53 and cell cycle networks via p53. Irradiation (IR) introduce to p53 compel p53 dynamics to suffer different phases, namely oscillating and oscillation death (stabilized) phases. The IR induced p53 dynamics undergo collapse of oscillation with collapse time Δt which depends on IR strength. The stress p53 via IR drive cell cycle molecular species MPF and cyclin dynamics to different states, namely, oscillation death, oscillations of periods, chaotic and sustain oscillation in their bifurcation diagram. We predict that there could be a critical Δtc induced by p53 via IRc, where, if Δt〈Δtc the cell cycle may come back to normal state, otherwise it will go to cell cycle arrest (apoptosis).

Introduction

p53 is well known for its abnormally long stability in response to the stress available against genomic integrity [1]. It conglomerated with its negative inhibitor MDM2 in the nucleus due to their strong interaction [2]. When the cell is in stress condition (due to irradiation, stress inducer molecule etc), p53 concentration level rises which leads to cell cycle arrest until repair or doctoring takes place of the impaired DNA. If the repair is not successful the system goes towards the apoptosis [3-6]. The transcriptional ability of the p53 is kept under controlled level at normal state due to its negative feedback interaction with MDM2 [7]. The hyperbolized concentration of MDM2 helps in degradation of the p53 protein because of its E3-ligase activity, causing adherence of ubiquitin to the lysine rich C-terminal of the p53 molecule [8-10]. Introduction of stress in the system is sensed by the activation of ARF protein, initially situated in nucleolar region in the form of nucleophosmin shifts to the nucleoplasm in its independent and active cast, to mark MDM2 for its degradation, thus assisting the p53 stability [11-13]. Triggering of p53 in response to stress leads to the expression of several downstream genes apart from the MDM2.

p21 protein is one of the most important proteins which is found to be expressed due to p53 accumulation in the cell [14]. p53 acts as a transcription factor for p21. It is also reported that p21 expression is directly proportional to the level of p53 in the system [15]. The role of p21 in controlling G1 phase checkpoint has been widely studied but its role in controlling G2 phase checkpoint is comparatively less studied [16-18]. The G2 phase checkpoint interruption leads to the disruption of cell cycle that leads to halt mitosis [14]. The cyclin-cdk interaction leads to the formation of MPF (Maturation Promoting Factor) [19]. The formation of MPF is very important for transition of G2 phase to mitosis phase [20]. The p21 protein is reported as antagonist for the formation of MPF. Several experimental results suggest that p21 directly interacts with cdk and also with cyclin leading to the inhibition of both cdk as well as cyclin [21]. It is also reported that the interaction of cdk and p21 causes to halt in DNA replication [20, 22].

Cyclin, in cell cycle process, is an important protein which interacts with cyclin dependent kinases and forms MPF. The MPF is responsible for the activation of pRb (Retinoblastoma protein), and helps the liberation of transcription factor E2F from its inhibitory. This E2F maintains the expression profile of genes required to ingress the S-phase of the cell division cycle [23-25]. Further, it is reported by several experimental results that p21 can directly interact with MPF and forms complex and then dissociate [16, 18]. Hence, p53 can able to cross talk with MPF and cyclin through p21.

There have been various experimental and theoretical studies on p53 regulatory network and cell cycle model to understand their regulatory mechanisms and cell fate. p53 – Mdm2 regulatory network has been modeled in order to study the impact of irradiation and change in DNA on cell variability and cell fate [26]. Further, it has also been shown that this DNA damage force the cell to select its fate (DNA repair, cell cycle arrest, apoptosis) via activating p53 [27]. On the other hand, variation in DNA methylation specially in neuronal cells in central nervous system may induce better response to developmental and environmental changes [28]. Moreover, this cell fate in tumor cells can probably be triggered by p53 dependent PUMA accumulation and p53 signal strength [29, 30]. Other method, say, recurrent artificial neural network model has also been implemented to study such network to understand DNA damage responses due to damage signal and parameter modeling to incorporate the changes [31, 32]. Studies in NF-kB model has been done in order to understand how the model system responses to the cellular signal which may trigger to different states like chaos in the dynamics and phase synchronization [33].

The experiments on mammalian cells show that p21-cyclin signaling pathway control the decision of cell cycle fate [34]. The other studies in cell cycle dynamics in mammalian cells further show the positive feedback as controlling mechanism of cell cycle regulation [35], role of noise in regulation and exhibition of bifurcation in cell cycle dynamics [36].

Our model incorporates the integration of both p53-Mdm2 regulatory network and cell cycle network in order to study the impact of p53 in deciding the fate of cell cycle dynamics and vice versa. We focus in this work to study and find out the behaviour of different molecular species which are actively involved in the checking of cell cycle at G2 phase regulated by p53. We proposed an integrated model of p53 and cell cycle network to find out the impact of p53 regulator on cell cycle via p21 protein. We organized our work as follows. We hope that the study may open up important behaviors in the dynamics of both p53 and cell cycle oscillators and in the decision making mechanism of cell fate via p53. We explained our proposed model in section II. The result of the large scale simulation of the model is given in section III with discussion. The conclusion based on our results is provided in section IV.

Materials and Methods

Model of cell cycle regulated by p53

We present a model which brings together p53 – MDM2 regulatory network [37] and cell cycle [38-40] via p21 protein (Fig 1) in the light of various theoretical and experimental reports. The model is described briefly as follows. The main component of p53 – MDM2 regulatory network is the feedback loop between p53 and MDM2 [37]. p53 and MDM2 interact to form p53 – MDM2 complex with a rate constant k 17 [37], followed by dissociation of the complex to its respective components with a rate constant k 18 [41, 42]. The transcription rate of MDM2 gene to its mRNA (MDM2 – mRNA) is takes place with rate constant k 20, followed by translation of MDM2 – mRNA to MDM2 with a rate constant k 22 [37, 43] and its (MDM2 – mRNA) self-degradation with a rate constant k 21 [44]. The ubiqiutination of MDM2 protein occurs with rate constant k 23. The p53 synthesis is taken placed with a rate constant k 16, and gets ubiquitinized at the rate constant k 19 [43]. The DNA damage in system is introduced via irradiation with an estimated rate constant of k 24 [37]. Irradiation is reported to be a major cause of DNA damage. The severity of the DNA damage is depended on the dose of exposure of irradiation [34]. The repair of the DNA damage is then occurred at a rate constant k 25 [45, 46]. The activation of ARF due to DNA damage takes place at a rate constant k 26 [37]. Further, this activated form of ARF interacts with MDM2 protein and forms ARF – MDM2 complex with a rate constant k 27 [47]. The degradation of ARF protein is reported to occur at a rate constant k 28 [48]. ARF based degradation of the MDM2 takes place by getting targeted to the complex via proteasome recognition with a much faster rate constant k 29 than individual degradation rates [49]. The p53, being a transcription promoting factor for many of the proteins, also transcribes the gene responsible for the manufacture of p21 protein with a rate constant k 30 as presumed by the approximations made to attain the appropriate oscillations and arrests [14, 50, 51].

Fig 1

The schematic diagram of interaction of p53-Mdm2 reaction network cell cycle oscillator.

The interaction between different molecular species are shown with respect to their rate constant. The blue and black dots indicate creation and decay of the respective molecular species.

The p21 protein is capable of making complex with the cell division promoting factor MPF with a rate constant k 31 [16, 18] with respect to the amount of concerned molecules present in the system [19]. Then the inhibition of MPF, or more appropriately G2 associated Cyclin – Cdk complex, by p21 is approximated with a rate constant k 32 [50, 52]. p21 then gets degraded by the virtue of its half-life in the system with a rate constant k 33 [18, 23]. The cyclin is assumed to translate at the rate constant k 1 [53]. Further, ubiquitin dependent cyclin degradation or protease independent degradation of the cyclin is reported to happen at a rate constant k* [54]. The degradation of the cyclin due to effect of protease activation during cyclin accumulation and interaction between inactive form of MPF with cyclin takes place with a rate constant k 4 [38]. Formation of activated form of MPF (M) occurs due to interaction of cyclin with inactive MPF (M*) with a rate constant k** [38, 55–57]. Further this activated form of MPF (M) converts to inactivated form (M*)with a rate constant k*** [55, 57]. The activated form of MPF(M) interact with inactive protease(X*) to generate activated form of protease (X) with a rate constant k**** [38, 58, 59]. The generation of activated form of cyclin protease (X) occurs due to interaction of cyclin protease with inactive X* with a rate k***** [38, 55]. The activated form of protease (X) can convert into inactive form (X*) with a rate constant k****** [38, 57]. In Fig 1, The blue dots indicates creation and black dots indicates decay of the respective molecular species. The lists of molecular species and biochemical reaction channels involved in this proposed model are listed in Tables 1 and 2 respectively.

Table 1

List of molecular species.

S.No.	Species Name	Description	Notation
1.	Cyclin	Unbounded Cyclin protein	x 1
2.	MPF	Maturation promotion factor	x 2
3.	Cyclin – Protease	Unbounded Cyclin Protease	x 3
4.	p53	Unbounded p53 protein	x 4
5.	Mdm2	Unbounded Mdm2 protein	x 4
6.	Mdm2_p53	Mdm2 with p53 complex	x 6
7.	Mdm2_mRNA	Mdm2 messenger mRNA	x 7
8.	ARF	Unbounded ARF protein	x 8
9.	ARF_Mdm2	ARF_Mdm2 complex	x 9
10.	IR	Irradiation	x 10
11.	DamDNA	Damaged DNA	x 11
12.	p21	p21 protein	x 12
13.	p21_M	p21 and M complex	x 13

Table 2

List of Chemical Reactions, Rate constants and their values.

S.No.	Biochemical reaction	Description	Rate Constant	Values of Rate Constant	Ref.
1	ϕ→k_1x_1	Synthesis of Cyclin	k 1	0.000416667 × 10−2 sec −1	[38–40]
2	x_1→k^*ϕ	Decay of Cyclin	k^*⟨x_1⟩,where,k^*=k_2x_1x_3k_3+x_1	k 2 = 0.004166667sec −1, k 3 = 0.02sec −1	[38, 54]
3	x_1→k_4ϕ	Cyclin decay	k 4⟨x 1⟩	0.0000167sec −1	[38]
4	ϕ→k^**x_2	Creation of MPF	k^**,where,k^**=k_5(1−x_2)k_6+(1−x_2),k_5=k_14x_1k_13+x_1	k 6 = 0.01, k 13 = 0.5, k 14 = 0.00sec −1	[16, 38, 57]
5	x_2→k^***ϕ	Decay of MPF	k^***⟨x_2⟩,where,k^***=k_7x_2k_8x_2	k 7 = 0.0025sec −1, k 8 = 0.01sec −1	[16, 55, 57]
6	x_2+x_12→k_31x_13	Formation of MPF_p21 complex	k 31⟨x 2⟩⟨x 12⟩	0.0001mol −1 sec −1	[16, 18, 19]
7	ϕ→k_7x_3	Activation of protease molecule	k^****,where,k^****=k_9(1−x_3)k_10+(1−x_3),k_9=x_2k_15	k 10 = 0.01, k 15 = 0.001667	[38, 58, 59]
8	x_3→k^*****ϕ	Inactivation of protease molecule	k^*****⟨x_3⟩,where,k^*****=k_11x_3k_12+x_3	k 11 = 0.0008333, k 12 = 0.01	[38, 57]
9	ϕ→k_16x_4	creation of p53	k 16	0.078	[37, 42, 43]
10	x_4+x_5→k_17x_6	synthesis of p53_MDM2 complex	k 17⟨x 4⟩⟨x 5⟩	1.155 × 10−3 mol −1 sec −1	[37, 42]
11	x_6→k_18x_4+x_5	Dissociation of p53_MDM2 complex	k 18⟨x 6⟩	1.155 × 10−5 sec −1	[37, 41, 42]
12	x_6→k_19x_5	ubiquitination of p53	k 19⟨x 6⟩	8.25 × 10−4 sec −1	[37, 42, 43]
13	x_4→k_20x_4+x_7	creation of MDM2_mRNA	k 20⟨x 4⟩	1.0 × 10−4 sec −1	[37, 42, 44]
14	x_7→k_21ϕ	decay of MDM2_mRNA	k 21⟨x 7⟩	1.0 × 10−4 sec −1	[37, 42, 44]
15	x_7→k_22x_5+x_7	synthesis of MDM2	k 22⟨x 7⟩	4.95 × 10−4 sec −1	[37, 42, 43]
15	x_5→k_23ϕ	decay of MDM2	k 23⟨x 5⟩	4.33 × 10−4 sec −1	[37, 42, 43]
16	x_10→k_24x_11	creation of DNA damage	k 24⟨x 10⟩	1.0sec −1	[37, 45]
17	x_11→k_25ϕ	recovery of damaged DNA	k 25⟨x 11⟩	2.0 × 10−5 sec −1	[37, 46]
18	x_11→k_26x_8	Activation of ARF	k 26⟨x 11⟩	3.3 × 10−5 sec −1	[37]
19	x_5+x_8→k_27x_9	synthesis of MDM2_ARF complex	k 27⟨x 5⟩⟨x 8⟩	0.01mol −1 sec −1	[37, 47]
20	x_8→k_28ϕ	decay of ARF	k 28⟨x 9⟩⟨x 8⟩	0.001sec −1	[37, 48]
21	x_9→k_29x_8	degradation of MDM2	k 29⟨x 9⟩	0.001sec −1	[37, 49]
22	x_4→k_30x_4+x_12	synthesis of p21	k 30⟨x 4⟩	0.001sec −1	[14, 50, 51]
23	x_2+x_12→k_30x_13	synthesis of p21_MPF complex	k 31⟨x 4⟩	0.0001sec −1	[14, 50, 51]
24	x_13→k_32x_12	dissociation of p21_MPF complex	k 32⟨x 13⟩	0.002sec −1	[16, 18, 50]
25	x_12→k_33ϕ	decay of p21 complex	k 33⟨x 12⟩	0.005sec −1	[16, 50]

The biochemical reaction network shown in Fig 1 are represented by the twenty five reaction channels listed in Table 2, which are participated by thirteen molecular species (Table 1) defined by a vector at any instant of time t, x(t) = {x 1(t), x 2(t), …, x (t)}, where, T is the transpose of the vector and N = 13. The variables are the concentrations of the molecular species. The time evolution of these variables can be translated from the twenty five reaction channels into the following set of nonlinear ordinary differential equations (ODE) based on Mass action law of chemical kinetics, where, the expressions for M* and X* in the Fig 1 are given by, M* = 1-x 10 and X* = 1-x 11. The set of coupled ODEs can be solved using Runge Kutta method of standard numerical integration algorithm [60].

Results and Discussion

We numerically simulate the proposed model and the results demonstrate new phenomena in bifurcation diagram which may be significant to correlate with various experimental situations. The interaction of p53 regulatory network and cell cycle network highlights different form of signal processing between non-identical networks which could be the way of regulating one another. We study the complicated way of this interaction in order to understand some of the basic mechanisms of network interaction.

Dynamics of p53 driven by irradiation

We first present the spatio-temporal behaviour of p53 upon exposure of irradiation in Fig 2. The p53 dynamics maintains minimum concentration level at IR = 0 (normal condition). As IR dose increases p53 start showing damped oscillatory behaviour (Fig 2 second and third panels) indicating stressed behaviour of p53. The increase in IR dose induces increase in time to attain stability of p53 dynamics (amplitude death) indicating increase in unstability of p53 dynamics (Fig 2 third panel). This could be due to the fact that the increase in IR dose may cause high DNA damage leading to more stress in p53.

Fig 2

Plot shows the temporal variation in the concentration and oscillatory pattern of p53 protein due to the effect of various exposure of IR (Gy) i.e (0,0.1,1,5,10) in left panels.

Similarly, temporal variation in the concentration and oscillatory pattern of MDM2 protein due to the effect of various exposure of IR (Gy) i.e (0,0.1,1,5,10) are shown in right panels.

However, if the IR dose is comparatively strong (IR = 5), the damage within the DNA is also high which may cause the collapse of the p53 oscillatory behaviour (Fig 2 fourth panel) and then repaired back the DNA damage to come back to p53 oscillatory condition. We also found that the time of collapse (Δt) increases as IR dose increases (Fig 2 fifth panel) and it becomes difficult to repair back the DNA damage. In general p53 will collapse forever and will not be recovered back if Δt → ∞ (probable case of apoptosis). However, in real situation, one probably can define a critical Δt such that, if Δt⟨Δt , p53 could come back after DNA repair, and otherwise it will go to apoptosis. Nevertheless, it is very difficult to find out this Δt .

Similarly, we also present the plots of temporal variation of the concentration of MDM2 due to exposure of irradiation in right panels of Fig 2. We observed similar kind of behaviour as obtained in case of p53 protein dynamics. This is probably due to intercorrelation between p53 and MDM2 in the system via feedback mechanism. It is also noted that corresponding variations in the behaviours of both p53 and MDM2 (as observed by comparing panels in Fig 2) are due to their positive as well as negative feedback regulations prescribed to them.

Phase diagram of p53 compelled by IR

We simulate the maxima of p53 amplitudes after removing the transients as a function of IR (Fig 3) to capture the different phases namely oscillation and oscillation death regimes. The behaviour of Δt as a function of IR follows the functional form with the values of A = 6778 and B = 0.00887 (fitting values of the function to the data) (Fig 3 inset). The separation between two phases oscillation death and oscillating regimes are clearly visible after the IR ∼ 3.45 and Δt increases as IR increases.

Fig 3

Plot for showing the impact of IR on p53 maxima.

Different p53 maxima observed at different values of IR (Gy) with respect to time. The p53 maxima verses IR dose is shown at left hand side inset and also IR dose verses time is shown in right hand inset.

Generally as Δt → ∞ when IR → ∞, but numerically we approximately found that after IR = R ∼ 11 Δt become Δt ∼ 79 hours and becomes constant (Fig 3 inset). This means that for any Δt⟨Δt , the p53 can able to recover back to normal stable state by repairing DNA damage, otherwise, the system can’t able to come back to normal state, but will go to apoptosis.

Bifurcation in Cyclin regulated by p53

Since cell cycle and p53 regulatory networks are interacted through p21 (Fig 1), the temporal behaviour of cyclin can be regulated by p53 via IR and p21. When IR = 0, the two networks work in normal condition, leaving p53 dynamics at low level (stabilized state) (Fig 2 upper panel) and sustain oscillation in cyclin dynamics (Fig 4 upper left panel). As IR increases, p53 will get activated through DNA damage giving oscillatory behaviour affecting the dynamics of cyclin. When IR = 0.1, the cyclin dynamics shows chaotic behaviour upto t = 145 hours, and then the dynamics becomes sustain oscillation (Fig 4 second left panel and upper right panel). The chaotic behaviour in cyclin dynamics could due to the sudden activation in p53 dynamics due to IR irradiation.

Fig 4

Plot shows the temporal variation in the oscillatory pattern of cyclin due to the effect of various exposure of IR (Gy) i.e (0,0.1,1,5,10) at left side panels and their corresponding bifurcation diagram are shown at right panels.

Now as IR increases (IR = 0.5), we get various situations in the cyclin dynamics, namely, the emergence of period two (for t ∼ [10-40] hours), period 3 (for t ∼ [40-85] hours), chaotic regime (for t ∼ [85-175] hours) and sustain oscillation regime (for t⟩175 hours) (Fig 4 second right upper panel). Further, as IR increases the emergence of oscillation death regime started to exist in the cyclin dynamics (Fig 4 fourth right panel onwards) and the oscillation death regime become larger. Further increase in IR compels the period 2 and 3 regimes to vanish after some value of IR (IR⟩9) and the chaotic regime becomes larger.

The perturbation induced by p53 through IR to the cyclin via p21 clearly induces cyclin dynamics to various states shown by the bifurcation diagram (Fig 4 right panels). We also notice that as one decrease or increase to cross over to sustain oscillation, the state just before it is chaotic regime. The emergence of oscillation death regime starts from IR⟩3 and then switches to sustain oscillation after sometime. This oscillation death regime corresponds to the collapse time due to strong sudden DNA damage. Once the DNA damage is recovered it comes back to sustain oscillation. If the IR is very large then oscillation death regime is large enough that DNA damage can not be repaired back halting the cell cycle permanantly and goes to apoptosis.

Dynamics of MPF regulated by p53

We present the temporal behaviour of MPF regulated by p53 as a function of IR (Fig 5) which induces at different states in MPF shown by bifurcation diagrams. The impact upon the MPF due to p53 via IR is not a direct phenomenon but through p21 molecule in the network. Various studies reported that p21 directly interact with cyclin dependent kinases, which has very important role in the formation of maturation promoting factor (MPF). The interaction of p21 with cdk leads to less availability of cdk due to the formation of MPF. Moreover, various experimental results also reported that p21 directly interacts with MPF [16, 18]. It is observed that an IR = 0, the MPF dynamics shows sustain oscillatory behaviour indicating no impact of p53. Further, as IR dose increases the oscillatory behaviour of MPF is abruptly changed inducing different states of MPF as we obtained in the case of cyclin. The increases in IR dose induce different states oscillation death, period 1, 2, 3, chaotic and sustain oscillation regimes indicated by the bifurcation diagram for various IR values. Moreover, as IR increases the width of oscillation death [16] regime also increases and if IR is not strong enough the DNA can able to repair back otherwise the system will go to apoptosis.

Fig 5

Plot shows the temporal variation in the oscillatory pattern of MPF(Maturation Promoting Factor) due to the effect of various exposure of IR (Gy) i.e (0,0.1,1,5,10) at left side panels and their corresponding bifurcation diagram are shown at right panels.

Bifurcation in MPF and Cyclin

We study the regulation of cell cycle dynamics by p53 via IR. The maxima values of MPF (MPF ) and cyclin (Cyc ) as a function of IR are calculated for a range of time in the range [0, 50] hours (Fig 6). It is observed that for low IR dose, MPF exhibits chaotic behaviour. Howeover, if IR dose is comparatively high, MPF becomes almost constant. If the value of IR is moderate, period 1, 2, 3 etc are exhibited in the bifurcation diagram. This indicates that MPF is p53 dependent via IR and p53 controls the MPF behaviour in the system.

Fig 6

Plot shows the impact of various IR dose (in Gy) on MPF maxima (at upper panel) as well as Cyclin maxima (at lower panel).

Similarly, one can also observe the IR dependent maxima of cyclin Cyc in the bifurcation diagram (Fig 6 lower panel). The moderate values of IR induce different periods in Cyc . Excess values of IR show different behaviour in Cyc .

Conclusion

We study the way how p53, one of the largest hubs in cellular network, regulates and controls cell cycle dynamics. We studied the behaviour of different molecules which are highly involved in the checking of cell cycle at G2 phase driven by p53 via IR. The simulation results of the model provided us to understand the biological phenomenon and mechanism of cell cycle arrest due to DNA damage faced by the cell due to the irradiation. The results we got are closely in agreement with the previous experimental reports [16, 17]. Our study suggests that the temporal dynamics of molecular species involved in cell cycle, considered in the model, are controlled by p53. The role of p21 protein in the delay of G2 phase was considered as a cross-talk between p53 regulatory network and cell cycle. The sudden irradiation to the system with high dose induces collapse of the system due to DNA damage, leading to cell cycle arrest. The cell cycle is resumed again to normal situation by repairing back the DNA damage. Moreover, the time of recovery from cell cycle arrest and then resumption of oscillation depends on the amount of dose of IR exposed to the system.

During the process of regulation of cell cyle by p53 via IR we observed the emergence of different periods (1, 2, 3 etc) in the bifurcation diagram of oscillatory dynamics of cell cycle variables (MPF and Cyc ) which may have various information of certain biological significance. Further, the dynamics of these variables switched to various states, namely, chaotic, oscillation death (stabilized state), bifurcating to various periods of oscillation and sustain oscillation states during the process of time evolution. These states could be the different phases of the variables to self-recover back to its normal condition from the sudden stress given to the system. However, how these complicated states are used by the system dynamics when the system is perturbed need to be investigated further.

The study also demonstrates the mechanism of cell cycle arrest induced by perturbed p53 via IR indicated by collapse of the oscillation (oscillation death) for certain interval of time (Δt). This collapse time is a function of strength of the perturbation imparted to the system. Our study shows that there is a minimum value of IR = R , below which the system comes back to its normal state, otherwise the system will go to apoptosis. Our findings will probably be useful for the further study on the impact of p53 on cell cycle checking at G2 phase and related dynamics.