Figures
Abstract
Experimental observations performed in the p53-Mdm2 network, one of the key protein modules involved in the control of proliferation of abnormal cells in mammals, revealed the existence of two frequencies of oscillations of p53 and Mdm2 in irradiated cells depending on the irradiation dose. These observations raised the question of the existence of birhythmicity, i.e. the coexistence of two oscillatory regimes for the same external conditions, in the p53-Mdm2 network which would be at the origin of these two distinct frequencies. A theoretical answer has been recently suggested by Ouattara, Abou-Jaoudé and Kaufman who proposed a 3-dimensional differential model showing birhythmicity to reproduce the two frequencies experimentally observed. The aim of this work is to analyze the mechanisms at the origin of the birhythmic behavior through a theoretical analysis of this differential model. To do so, we reduced this model, in a first step, into a 3-dimensional piecewise linear differential model where the Hill functions have been approximated by step functions, and, in a second step, into a 2-dimensional piecewise linear differential model by setting one autonomous variable as a constant in each domain of the phase space. We find that two features related to the phase space structure of the system are at the origin of the birhythmic behavior: the existence of two embedded cycles in the transition graph of the reduced models; the presence of a bypass in the orbit of the large amplitude oscillatory regime of low frequency. Based on this analysis, an experimental strategy is proposed to test the existence of birhythmicity in the p53-Mdm2 network. From a methodological point of view, this approach greatly facilitates the computational analysis of complex oscillatory behavior and could represent a valuable tool to explore mathematical models of biological rhythms showing sufficiently steep nonlinearities.
Citation: Abou-Jaoudé W, Chaves M, Gouzé J-L (2011) A Theoretical Exploration of Birhythmicity in the p53-Mdm2 Network. PLoS ONE 6(2): e17075. https://doi.org/10.1371/journal.pone.0017075
Editor: Magnus Rattray, University of Sheffield, United Kingdom
Received: October 5, 2010; Accepted: January 20, 2011; Published: February 14, 2011
Copyright: © 2011 Abou-Jaoudé et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was supported in part by the INRIA-INSERM initiative action ColAge. The authors thank INRIA for funding the post-doctoral stay of WAJ. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. No additional external funding received for this study.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Periodic phenomena are encountered at all levels of biological organization, with periods ranging from fractions of a second to years [1]. At the intracellular level, periodic phenomena have been reported in various biochemical systems such as calcium signalling, circadian rythms, cell cycle, glycolysis, cAMP signaling in Dictyostelium [1], [2] or the p53-Mdm2 network [3]. Most of the time, biochemical oscillations display a simple pattern with a single oscillatory regime of stable period and amplitude. However, rhythmic processes can sometimes present a more complex behavior. One mode of complex oscillatory behavior is the coexistence of two simultaneously stable oscillatory regimes for the same external conditions. This phenomenon, called birhythmicity [1], [4], is the counterpart of bistability for oscillatory dynamics. Such a behavior has been observed in a number of chemical oscillatory systems [5], [6] but, although some studies suggest its occurrence in the heart and the neuronal system [7], birhythmicity has not yet been firmly observed experimentally in biological systems.
Recent experiments performed in the p53-Mdm2 network, one of the key protein module involved in the control of proliferation of abnormal cells in mammals [8], [9], [10], reported two oscillatory regimes of p53 and Mdm2 in irradiated cells [11]: a low-frequency oscillatory regime at low irradiation dose with a period of about 10 h and high-frequency oscillations at high irradiation dose with a period of about 6 h (Figure 1). This observation raised the question of the existence of birhythmicity in the p53-Mdm2 network which would be at the origin of the two oscillatory regimes experimentally observed as a function of the irradiation dose.
Power spectrum of nuclear Mdm2-YFP fluorescence dynamics in individual cells. Top: an example of a cell showing fluctuations with a characteristic frequency of ∼10 hours (exposed to 0.3Gy of gamma irradiation), and the power spectrum of the signal (by Fourier transform). Bottom: an example of a cell with multiple oscillations with a period of ∼6 hours (exposed to 5Gy), and the power spectrum of the signal (right). Reprinted by permission from Macmillan Publishers Ltd: Molecular Systems Biology, advance online publication, 2006 (doi: 10.1038/msb4100068).
A theoretical answer to this question has been recently proposed by Ouattara, Abou-Jaoudé and Kaufman [12], [13]. In the framework of a simple 3-dimensional differential model of the p53-network, they showed that this system could display birhythmicity for a certain range of irradiation dose [12]. The simultaneous presence of two distinct periodic orbits allowed the model to reproduce in particular (i) the two frequencies experimentally observed, (ii) the increase of the fraction of cells oscillating with a high frequency when the irradiation dose increases and (iii) the changes in the oscillation frequency that have been observed for some cells after irradiation [12], [13].
Following this work, the aim of this paper is to investigate the mechanisms at the origin of birhythmicity in Ouattara, Abou-Jaoudé and Kaufman's model (OAK model). As this 3-dimensional continuous non-linear differential model is difficult to analyze, in a first step, we approximated it by a 3-dimensional piecewise linear differential model where the Hill functions have been replaced by step functions and, in a second step, by a 2-dimensional piecewise linear differential model by setting one autonomous variable as a constant in each domain of the phase space delimited by the thresholds of the step functions. Analyzing the dynamics of the system in the framework of a piecewise linear differential description allowed us to reveal the phase space structure of the two oscillatory regimes composing the birhythmic behavior.
Following this analysis, we find that two features related to the phase space structure of the system are at the origin of birhythmicity observed in the OAK model: (1) the existence of two embedded cycles in the transition graph of the piecewise differential models; (2) the presence of a bypass containing two folds in the orbit of the large amplitude oscillatory regime of low frequency. Calculation of a first return map further permitted to give an interpretation of the role of these two features in the emergence of the birhythmic phenomenon. Finally, from this analysis, an experimental strategy is proposed to test the existence of birhythmicity in the p53-Mdm2 network.
Results
Birhythmicity in the OAK Model
The model of the p53-Mdm2 core network proposed by Ouattara, Abou-Jaoudé and Kaufman is a three-variable model, inspired by the work of Ciliberto et al. [14], which describes the interaction between p53, cytoplasmic Mdm2 and nuclear Mdm2. We briefly recall the biological data from which this model has been elaborated (see [12], for a more detailed description). Nuclear Mdm2 down regulates p53 by accelerating its degradation through ubiquitination [15]–[17] and by blocking its functional activity [18]–[20]. p53 up regulates cytoplasmic Mdm2 level by enhancing the transcription of gene MDM2 [21], [22]. p53 inhibits the translocation of Mdm2 from the cytoplasm to the nucleus [23]. These interactions are summarized by the influence diagram shown in Figure 2. Importantly, the degradation rate of nuclear Mdm2 (dMn) depends on DNA damage caused upon cell irradiation [24], [25] and increases when DNA damage increases [13].
Normal arrows correspond to positive interactions, blunt arrows to negative interactions. P, Mc and Mn represent p53, cytoplasmic Mdm2 and nuclear Mdm2 respectively.
To describe the dynamics of this network, Ouattara, Abou-Jaoudé and Kaufman built a three-dimensional differential model (which will be called the OAK Model from now on) to represent the temporal evolution of the level of p53, cytoplasmic and nuclear Mdm2. We recall the equations of the OAK Model in Text S1. For appropriate parameter values, a bifurcation analysis of this model as a function of dMn showed the existence of a domain of birhythmicity separating two oscillatory regimes of different amplitude, period and mean values (see Figure 2b in [13], and Figure S1) which allowed to reproduce the two characteristic periods of p53 oscillations observed experimentally by Geva-Zatorsky et al. as a function of cell irradiation dose (Figure 1).
A three-dimensional reduced model
Starting from the OAK model, we first performed an analysis of the different terms involved in the equations of this model to characterize their contribution to the oscillatory dynamics of the system, in particular to the birhythmic behavior. Our analysis suggests that the terms modeling the translocation process of Mdm2 from the nucleus to the cytoplasm and the bilinear term modeling Mdm2-mediated acceleration of p53 degradation would have little influence on the oscillatory dynamics (see Text S2). Suppressing these terms, the OAK Model becomes:(Model 1)where P, Mc and Mn represent the concentration of p53, cytoplasmic Mdm2 and nuclear Mdm2 respectively.
Numerical simulations of Model 1 show that, for appropriate parameter settings, the bifurcation diagram as a function of dMn is very similar to the bifurcation picture of the OAK Model if the non-linear terms are slightly steeper (i.e. higher Hill coefficients n) (Figure 3A, and Figure S1). As dMn increases, Model 1 displays successively a stable steady state with a low level of p53, the coexistence of this stable steady state with two unstable steady states, a stable oscillatory state, and finally a stable steady state with high levels of p53. In the oscillatory domain, Model 1 shows two limit cycles, a large amplitude oscillatory regime of low frequency for low dMn values and a small amplitude oscillatory regime of high frequency for high dMn values, separated by a birhythmic region where these two regimes coexist (Figure 3A). In the domain of birhythmicity, numerical simulations moreover show that the phase portrait of the two limit cycles for the OAK Model is very similar to the one for Model 1 (Figure 3B, and Figure S1). In the remainder of this work, Model 1 will thus be used as the starting model to investigate the mechanisms leading to the emergence of the birhythmic behavior in the OAK Model.
(A) Bifurcation diagram of p53 level as a function of dMn for Model 1. Solid lines (resp. dashed lines) represent the stable (resp. unstable) equilibrium points. Bold (resp. white) dots are the maxima and minima of the stable (resp. unstable) limit cycles. The system shows a birhythmic domain for 1.75 h−1<dMn<2.11 h−1. (B) Projection of the two oscillatory regimes on the plane (Mn,P) for Model 1 for dMn = 1.9 h−1. The thresholds KMc and KMn related to P are indicated in solid lines. The parameter values of Model 1 are the same as for the OAK Model (KMn = 0.1 nM) except: n = 6, dP = 2.5 h−1 and KMc = 0.4 nM (see legend of Figure S1 for the parameter values of the OAK Model).
We next qualitatively analyzed the main differences between the two oscillatory regimes for Model 1. Our analysis shows significant differences between the shape of the orbits of the two oscillatory regimes in the phase space (Figure 3B). Indeed:
- the small amplitude oscillatory regime of short period has a regular shape. Moreover, p53 level stands well above KMn threshold all along the trajectory of the cycle. The Hill function
modeling the inhibition of Mdm2 translocation into the nucleus by p53 has thus no significant impact along this cycle
.
- on the contrary, the large amplitude oscillatory regime of long period contains a long “tail” in the phase space domain where p53 is below KMn threshold and its orbit crosses the two thresholds, KMn and KMc, related to P.
Finally, we performed an analytical study of the dynamics of the system by investigating general properties concerning the existence of limit cycles. Contrary to 2-dimensional differential systems where the Poincaré-Bendixson theorem gives general conditions for the existence of periodic orbits, there is no such equivalent theorem concerning the existence of periodic orbits for 3-dimensional differential system except for particular systems such as competitive or cooperative systems [26]. In the case of Model 1, an analysis of the sign of the elements of the Jacobian Matrix interestingly shows that this model is indeed equivalent to a competitive system (Text S3). According to Hirsch [27], a Poincaré-Bendixson theorem in dimension 3 holds. Moreover, for the parameter values for which the system displays birhythmicity, a numerical study of the equilibrium points suggests that there is only one equilibrium point which is unstable and for which the stable manifold is 1-dimensional (Text S3 and Figure S2). Applying the Poincaré-Bendixson theorem for 3-dimensional competitive systems, it can be proved that the system admits a periodic orbit for the parameter values for which the system displays birhythmicity [26]. Under some supplementary assumptions satisfied by our system, we can moreover show that there exists a stable periodic orbit [28]. However, these theorems do not tell how many periodic orbits the system has and thus do not provide a mathematical proof of the existence of birhythmicity in our model.
From continuous to discrete model: approximation of the Hill functions by step functions
Although we can derive some general properties concerning the dynamics of the 3-dimensional differential Model 1, it is still difficult to get insight into the dynamics of the system, in particular into the birhythmic behavior.
However, the equations of Model 1 involve several highly non-linear terms of the Hill form:
for increasing Hill functions;
for decreasing Hill functions.
We thus approximated Model 1 by replacing these terms by step functions with two levels (Figure 4):for increasing Hill functions and:
for decreasing Hill functions, where ν represents the threshold of the step functions.
Hill function which appears in the equation of Mc in Model 1 (in black) and its approximation into the step function
defined as:
if
and
if
(in red) for n = 6 and KMc = 0.4 nM. n is the Hill coefficient and characterizes the steepness of the Hill function.
Such an approximation has been originally proposed by Glass and Kaufman [29] for biochemical networks and has been since widely used to model genetic systems [30]–[33]. This approximation leads to a piecewise linear differential model whose equations are:(Model 2)
The space of variables can thus be decomposed into 6 domains (D11, D21, D12, D22, D13, D23) delimited by the threshold values of the step functions: KP, KMn and KMc (Figure 5A). The equations of evolution in each domain of the phase space are detailed in Table S1. In each domain, the equations are affine and stable and one can calculate a so-called target equilibrium point of the domain towards which the system will tend (Table S2). The target equilibrium point is an analogous of the focal point defined in a class of piecewise linear diagonal models [34] (see next section). If a target equilibrium point of a domain Dij does not belong to its domain, then the system starting from Dij will leave Dij at some time as it will reach sooner or later a boundary of the domain. If a target equilibrium point of a domain Dij belongs to Dij, it corresponds to a stable equilibrium point of the system. Importantly, the equation of Mn depends on Mc. It follows that the sign of the derivative of Mn can change in each domain of the phase space according to Mc. The transition graph, which is the graph defining the possible transitions between the different domains [35], can thus not be directly derived from the position of the target equilibrium points (see next section).
(A) Subdivision of the phase space for Model 2 in 6 domains delimited by the thresholds KP, KMn and KMc. (B) Graph of the transitions followed by the two oscillatory regimes composing birhythmicity shown in Figure 6. The small amplitude oscillatory regime of short period crosses domains D22, D12, D13 and D23 successively (in red). The large amplitude oscillatory regime of long period crosses domains: D22, D21, D11, D12, D13 and D23 successively (in green).
For appropriate parameter settings, numerical simulations show that the bifurcation diagram as a function of dMn is similar to the bifurcation picture of Model 1: a stable steady state of low p53 values for low dMn values corresponding to the target equilibrium point of domain D21 (P = 0, ,
), a stable steady state of high p53 levels for high dMn values corresponding to the target equilibrium point of domain D13 (
,
,
) and an oscillatory regime for intermediate dMn values (see Text S4 and Figure 6). In the oscillatory region, there exists a range of dMn values for which the system displays birhythmicity with the coexistence of two oscillatory regimes of different amplitude and period (Figure 6):
- a small amplitude oscillatory regime of short period, crossing domains D12, D22, D13 and D23, in which the KMn threshold is not functional (i.e. P>KMn along the orbit). This periodic orbit corresponds to the small amplitude oscillatory regime appearing in Model 1 (Figure 3B);
- a large amplitude oscillatory regime of long period passing through all the domains of the phase space. This periodic orbit corresponds to the large amplitude oscillatory regime of Model 1 (Figure 3B) and contains a broad outward excursion in domain D21, in the form of a “tail”. As also observed for Model 1, this “tail” contains two folds where the sign of dMn/dt changes: a fold in the transition between D22 and D21 and a fold inside domain D21.
Projection of the two oscillatory regimes composing the birhythmic behavior of Model 2 in the plane (Mn,P). The dashed lines represent the thresholds of the step functions: KMn and KMc for P, KP for Mn. The parameter values are the same as for Figure 3 (KMn = 0.1 nM, KP = 0.2 nM) except KMc = 0.6 nM. The period of the large amplitude oscillatory regime is significantly longer than the period of the small amplitude oscillatory regime (see Figure 8A and 8C).
The sequence of transitions followed by the small and the large amplitude limit cycles are shown in Figure 5B. As emphasized in the previous section, the orbits of the two limit cycles are significantly different. Using Model 2, we can clearly separate the trajectories of the two limit cycles by decomposing them into a sequence of transitions in the state space. Each of the two orbits evolves in a dynamical regime with distinct qualitative properties and distinct biological consequences. Starting from the threshold delimiting domain D12 and D13 (i.e. p53 levels in its highest threshold KMc and nuclear Mdm2 levels is below its threshold KP), the two oscillatory regimes can be described as follows. As p53 level goes above the KMc threshold, p53 starts to activate the expression of gene MDM2 and thus the production of cytoplasmic Mdm2. Mdm2 is then translocated into the nucleus where it accumulates. When nuclear Mdm2 level gets beyond its threshold KP, it starts to inhibit p53 activity (domain D23). p53 level then decreases. As its level goes below its threshold KMc, p53 no longer activates the expression of MDM2. Cytoplasmic Mdm2 thus starts to decrease leading to a decrease of nuclear Mdm2 level (domain D22). As p53 and nuclear Mdm2 levels both decrease, two alternatives can be chosen by the system leading to either the small or the large amplitude oscillatory regime, which corresponds to the branching point in the graph shown in Figure 5B:
- either nuclear Mdm2 level first goes below its threshold KP. Nuclear Mdm2 then no longer inhibits p53 and p53 level starts to increase again (domain D12) until it reaches its threshold KMc which closes the trajectory of the small amplitude limit cycle.
- or p53 level first goes below its threshold KMn. p53 then no longer inhibits Mdm2 entry into the nucleus (domain D21) and nuclear Mdm2 starts to rise again, leading to the first fold observed in the transition between D22 and D21 in the large amplitude limit cycle. However, at the same time, cytoplasmic Mdm2 decreases because p53 level is too low to activate MDM2 transcription. When cytoplasmic Mdm2 level gets too low to feed nuclear Mdm2, nuclear Mdm2 level then starts to decrease leading to the second fold observed in the large amplitude limit cycle in domain D21. As nuclear Mdm2 level goes below its threshold KP, nuclear Mdm2 no longer inhibits p53 activity whose level starts to increase again (domain D11 and D12) until it reaches its threshold KMc which closes the trajectory of the large amplitude limit cycle.
The analysis of the orbits of the two oscillatory regimes thus reveals two key qualitative features which characterize the birhythmic behavior:
- the choice between two transitions when the system reaches domain D22 which induces a separation in the phase space between the two oscillatory regimes. This branching point gives rise to two embedded cycles in the graph of transitions representing the two oscillatory regimes in the phase space (Figure 5B). This result is in accordance with the logical analysis of the p53-Mdm2 network which brought out an oscillatory regime composed of several small and large amplitude embedded cycles in the transition graph of the system (see Figure 3 in [12]);
- the presence of a “tail” containing two folds in the orbit of the high amplitude oscillatory regime (Figure 6, domain D21) whose emergence can be interpreted in terms of the competition between two opposite processes: Mdm2 translocation from the cytoplasm to the nucleus on one hand and p53-mediated inhibition of Mdm2 nuclear entry on the other hand.
Analysis of the basic mechanisms for the emergence of birhythmicity
A two-dimensional reduced model.
In order to get more insight into the basic mechanisms leading to birhythmicity, we looked for a further approximation of Model 2 which would preserve the birhythmic behavior as well as the characteristics of the orbits of the two oscillatory regimes underlined by the previous analysis. In Model 2, in each domain, the evolution of Mn depends on Mc whereas the evolution of Mc does not depend on the other variables of the system (Table S1). We thus considered Mc as a forcing external parameter applied on the evolution of Mn. In this respect, we set Mc as a constant Mcij in each domain Dij of the phase space under some constraints on the Mcij values which will be detailed further.
This approximation leads to a 2-dimensional piecewise linear differential model with the phase space delimited by the thresholds KMc, KMn and KP. The evolution of each variable in each domain is of the form:where x represents the level of P or Mn, kij is a constant which depends on the domain Dij of the phase space and d is the degradation rate of P or Mn. This type of piecewise diagonal linear systems belongs to a class of dynamical systems proposed originally by Glass and Kaufman [29]. Interestingly, these systems have mathematical properties which favor the qualitative analysis of their dynamics [35]–[38]. In particular, one can derive the so-called transition graph, which describes all the possible transitions between the different domains, from the position of the target equilibrium points of each domain (called the focal points) [31], [35], [36], [38].
In order to reproduce the two folds forming the outward excursion which characterizes the large amplitude oscillatory regime in Model 1 and Model 2, we added another threshold K (K<KMn) for p53 level. As the evolution of Mn in each domain is now monotone, the introduction of this new threshold allows changing the sign of the derivative of Mn as the system crosses threshold K and thus recovering in particular the fold observed in domain D21 in Model 2 (Figure 6). The space of variables can thus be decomposed into 8 domains (D11, D21, D12, D22, D13, D23, D14, D24), defined by the thresholds KMc, KMn plus the additional threshold K for p53, and KP for nuclear Mdm2 (Figure 7). The equations of evolution in each domain of the phase space are detailed in Table S3 (Model 3).
(A) Subdivision of the phase space for Model 3 in 8 domains delimited by the thresholds KP, KMn, KMc and the additional threshold K (in red). (B) Transition graph of Model 3. The graph contains a branching point in domain D23 and two embedded cycles. From this domain, the system can either go to domain D13 or domain D22 depending on the initial conditions.
Moreover, in order to keep the basic characteristics of the birhythmic behavior shown in Model 2, we imposed the following constraints on the parameters values of Model 3 (see Table S4 and Text S5 for more details):
- constraint (1): the transition graph of Model 3 contains the sequence of transitions corresponding to the two oscillatory regimes in Model 2 (Figures 5 and 7).
- constraint (2): the setting of the parameters Mcij in each domain Dij of the phase space of Model 3 has to be in accordance with the evolution of Mc in the two oscillatory regimes from one domain to another in Model 2 (Figure 8).
Temporal simulation (h) of the concentration (in nM) of p53 (P), cytoplasmic Mdm2 (Mc) and nuclear Mdm2 (Mn) for Model 2 (A and C) and Model 3 (B and D) in the small amplitude short period (A and B) and the large amplitude long period (C and D) oscillatory regime. For Model 3, Mc has been set as a constant Mcij in each domain Dij of the phase space following constraint (2) (see text). The parameter values are indicated in Figure 6 for Model 2 and in Figure 9 for Model 3.
Finally, in order to simplify the analytical calculation of the return map (see section on the analysis of the first return map), the values of the degradation rate, dMn and dP, have been chosen to be the same. In this case, the trajectories in each domain of the phase space are straight lines, which greatly reduce the explicit computation of the return map. However, considering different degradation rate values does not change the main results of our analysis. Indeed, we can still numerically recover birhythmicity and the basic characteritics of this behavior for nonequal degradation rates (not shown).
Numerical simulations of the 2-dimensional model.
For appropriate parameter settings respecting the constraints stated above, numerical simulations show that the system exhibits a birhythmic behavior whose phase portrait is similar to the projection in the plane (Mn,P) of the phase portrait of the birhythmic behavior observed in Model 2 (Figure 9). In accordance with constraint (1), the small amplitude limit cycle crosses domains D14, D24, D23 and D13 whereas the large amplitude limit cycle passes through all the domains of the phase space. As shown in the transition graph (Figure 7B), when the system reaches domain D23, it can either go to D13 or D22. This branching point thus induces a separatrix in D23, passing through the threshold intersection (Mn = KP, P = KMn) and the focal point of D23, which delimits the basins of attraction of the two oscillatory regimes of the birhythmic behavior (Figure 10A and 10C, blue curve). Finally, the large amplitude cycle contains a “tail” with two folds which appear at two domain transitions: D23 to D22 when the level of P crosses KMn and D22 to D21 when the level of P level crosses the additional threshold K. These folds reproduce the two folds observed in the large amplitude limit cycle in Model 2 (see Figure 6).
Simulation of the two oscillatory regimes in the phase space for Model 3. The dashed lines represent the thresholds KMc, KMn and K for P, KP for Mn. The parameter values are the same as for Figure 6 (KMn = 0.1 nM) except KMc = 0.4 nM, K = 0.05 nM, KP = 2 nM, dP = 3 h−1, Mc11 = Mc21 = Mc12 = 0, Mc22 = 5 nM, Mc13 = 9 nM, Mc23 = 11.3 nM, Mc14 = Mc24 = 25 nM. Note that, since the degradation constants dMn and dP have the same values, the trajectories in each domain are straight lines.
(A and C) In red, simulation in the phase plane for the initial conditions Mn = 2 nM and P = 0.515 nM (x>xD, panel A) or P = 0.401 nM (x<xD, panel C) for Model 3. The dashed lines represent the thresholds KMc = 0.4 nM, KMn = 0.1 nM and K = 0.05 nM for P, KP = 2 nM for Mn. (B and D) In red, simulation of the orbit of the first return map from and to the [0, x) axis for the initials conditions of panel A (panel B) and panel C (panel D). The first return map F is shown in green and has been derived analytically (see Text S6). The parameter values are indicated in Figure 9. The discontinuity in F at x = xD∼0.077 arises when the trajectory hits the threshold intersection point (Mn = KP = 2, P = KMn = 0.1) (blue curve). The fixed points of F correspond to the two limit cycles of the system shown in Figure 9. For x>xD (resp. x<xD), the trajectory tends to the large (resp. small) amplitude oscillatory regime corresponding to the fixed point x = x2∼0.093 (resp. x = x1∼0.048).
Analysis of a first return map of the 2-dimensional model.
We next looked for a proof of the existence of birhythmicity by analyzing the first return map from and to the boundary between the domains D14 and D24 ([0, x) axis in Figure 10) which is crossed by both oscillatory regimes. The return map (or Poincaré map) of a surface S (here the [0, x) axis) is a mapping of S obtained by following the trajectory between two consecutive intersections with S [39]. This mapping is described by the equation:where xn (resp. xn+1) represents the nth (resp n+1th) intersection of the trajectory of the system with S, and F represents the first return map. In this description, the periodic orbits crossing S correspond to fixed points of the function F (i.e. points x for which
).
Here, an analytical expression of the first return map of the [0, x) axis can be calculated (see Text S6). The results of our analysis show that the return map (Figure 10B and 10D, green curve) presents two strictly positive fixed points which correspond to the two oscillatory regimes shown by the numerical simulations in Figure 9, proving the existence of birhythmicity in Model 3. Moreover and surprisingly, a discontinuity point appears which separates the [0, x) axis in two intervals: the interval [0, xD] for which the system tends to the small amplitude oscillatory regime (x = x1) (Figure 10C and 10D); the interval x>xD for which the system tends to the large amplitude oscillatory regime (x = x2) (Figure 10A and 10B). For the critical point x = xD, the trajectory hits the threshold intersection point (Mn = KP, P = KMn) (Figure 10A and 10C, blue curve), drawing the separatrix of the system in domain D23. Therefore, the branching point of the transition graph, which induces this separatrix, is at the origin of the discontinuity point observed in the return map.
We then investigated the role of the outward excursion which characterizes the large amplitude oscillatory regime in the emergence of the birhythmicity by using the return map description. To do so, we suppressed this “tail”, while keeping the two embedded cycles in the transition graph, by adding a transition from domain D22 to domain D12 in the transition graph (Figure 11A). We then calculated the first return map in the [0, x) axis for this modified Model 3 (see Text S6 for details of the calculation). The results show that the large amplitude limit cycle disappeared (Figure 11B). Indeed, the first return map ceases to be discontinuous and instead presents a non smooth point for the same xD value than before (Text S6). The loss of the discontinuity point accompanies the loss of the fixed point, x = x2, corresponding to the large amplitude oscillatory regime. Therefore, the “tail” is here necessary for the emergence of the birhythmic behavior. More precisely, the suppression of the “tail” leads to the disappearance of the discontinuity of the first return map which has for consequence the disappearance of the large amplitude limit cycle.
(A) Transition graph of the modified model 3. (B) Graph of the corresponding first return map (in green) from and to the [0, x) axis (see Figure 10). The parameter values are indicated in Figure 9, except Mc22 = 1. The focal point of D22 is now in D11 inducing the additional transition from D22 to D12 (arrow in red). The fixed points of the first return function correspond to limit cycles of the system. The system presents a unique limit cycle (x = x1∼0.048) corresponding to the small amplitude oscillatory regime (see Figure 10C and 10D). The first return map shows a non smooth point at xD∼0.077 (see Text S6).
An experimental strategy to reveal birhythmicity in the p53-Mdm2 network
From the previous analysis, an experimental strategy can be proposed to test the existence of the birhythmic behavior predicted by the OAK Model. The experiments performed on irradiated cells by Geva-Zatorsky et al. showed the existence of two frequencies of p53 oscillations as a function of the irradiation dose but the co-existence between these two types of oscillations for the same irradiation dose has not been firmly observed experimentally so far [11]. The 2-dimensional piecewise linear reduced system (Model 3) provides a model of the p53-Mdm2 network from which experimental strategies can be implemented to distinguish between existence of a single oscillatory regime (mono-rhythmicity) whose frequency increases with the irradiation dose, or coexistence of two oscillatory regimes for the same irradiation dose (birhythmicity).
Looking at the phase portrait of the oscillatory regimes for Model 3 (Figure 9), it's straightforward to predict that, starting from the small amplitude limit cycle of short period, an increase in p53 concentration will lead to a shift of the system to the other oscillatory regime. Indeed, numerical simulations of this model show that applying a transitory pulse of p53 when the system oscillates in the short period oscillatory regime leads to a shift to long period oscillations (Figure 12). It has to be noted that these type of strategies in which a transient pulse of an inducer is applied to provoke a permanent switch from one regime to another have already been implemented to reveal bistable behavior in genetic networks [40]–[42]. We thus propose the following experimental test: introduce an inducible P53 gene into the cells. Then apply to the cells a set of increasing irradiation doses for which most of the cells are oscillating with a high frequency. For each dose, apply a transitory pulse of the inducer (which will induce a transitory expression of the inducible P53 gene) and wait until the distribution of the periods among the cell population reaches a stationary regime. Has the distribution shifted to long period p53 and Mdm2 oscillations? If yes, this experimental observation is consistent with the existence of a birhythmic behavior in the system. However, in order to show the coexistence of two distinct period distributions for the same dMn values, the DNA repair process triggered by cell irradiation has to be slow (i.e. dMn decreases slowly as DNA damage is repaired) compared to the time required to reach the stationary regime after a p53 pulse induction. If this assumption is not fulfilled, a more elaborated experiment, in which the repair process would be slowed down, has to be designed.
Temporal simulation of p53 level for Model 3. A pulse of p53 is applied at t = 1.2 h (dashed line). Before applying the pulse, the system is oscillating in the small amplitude oscillatory regime. The p53 pulse induces a shift from the small amplitude short period to the large amplitude long period limit cycle. The parameter values are indicated in Figure 9.
Discussion
The aim of this work has been to analyze the mechanisms at the origin of the birhythmic behavior observed in Ouattara, Abou-Jaoudé and Kaufman's differential model (OAK Model) of the p53-Mdm2 network. To do so, the Hill functions have been firstly approximated into step functions leading to a 3-dimensional piecewise linear differential model in which the birhythmic behavior is conserved. An autonomous variable of the system has then been set as a constant in each domain defined by the thresholds of the step functions to obtain a 2-dimensional piecewise linear differential model also displaying birhythmicity.
Converting the sharp non-linearities of the model into discontinuous step functions revealed the phase space structure of the birhythmic behavior. This approximation coupled with a reduction of the dimensionality of the system, from a 3-dimensional to a 2-dimensional system, and then from 2-dimensional to a 1-dimensional first return map description, yielded a reduced model from which the basic qualitative features observed in the phase space leading to birhythmicity have been extracted: (1) the presence of two embedded cycles in the transition graph of the piecewise linear models, (2) the presence of a “tail” in the orbit of the large amplitude oscillatory regime of long period.
By calculating the first return map, we have theoretically shown the existence of birhythmicity in the framework of the 2-dimensional piecewise linear model. This map also elucidates the role of these two features in the emergence of the birhythmic phenomenon: (1) is associated with a non smooth point of the return map whereas (2) leads to a discontinuity in the return map, which is at the origin of the emergence of the birhythmic behavior. In general, on a methodological level, our analysis shows that the mechanisms leading to complex dynamical behavior such as birhythmicity can be studied by developing qualitative models of the system. One major advantage of this approach is to greatly facilitate the computational analysis of the model. In this regard, this method could be implemented to analyze the dynamics of other mathematical models of biological rhythms showing sufficiently steep nonlinearities.
Finally, in the light of our analysis, we proposed an experimental strategy to test the existence of birhythmicity in the p53-Mdm2 network. This strategy relies on a key characteristic of multistable systems which is to convert a transient stimulus into a permanent response of the cell. However, one may wonder what the physiological interest of birhythmicity in the p53 response to DNA damage is. One possible advantage could be to enhance the response of the system in contrast to mono-rhythmic systems which may require large perturbations and long time to induce a significant frequency change. Birhythmicity would thus allow a better flexibility of the p53 response by permitting the cell to quickly switch from one oscillatory regime to another after DNA damage.
Materials and Methods
The simulations are performed with XPPAUT (http://www.math.pitt.edu/~bard/xpp/xpp.html) for Model 1 and with Matlab for Model 2 and 3.
Supporting Information
Table S1.
Equations of evolution for Model 2. The domains Dij of the phase space are delimited by the threshold values of the step functions: KP, KMn and KMc (Figure 5A).
https://doi.org/10.1371/journal.pone.0017075.s001
(DOC)
Table S2.
Target equilibrium points for Model 2. In each domain Dij, the equations of evolution are linear and the Jacobian matrix is triangular with negative elements in the diagonal (see Table S1). In each domain Dij of the phase space, the system will thus tend towards the target equilibrium of Dij.
https://doi.org/10.1371/journal.pone.0017075.s002
(DOC)
Table S3.
Equations of evolution for Model 3. The domains Dij of the phase space are delimited by the threshold values of the step functions KP, KMn, KMc and the additional threshold K. Mc has been set as a constant Mcij in each domain Dij.
https://doi.org/10.1371/journal.pone.0017075.s003
(DOC)
Table S4.
Conditions on the parameter values of Model 3 to respect constraint (1). The focal points have been chosen such that the transition graph contains the two embedded cycles composing the graph of transitions of Model 2 (Figure 6): one cycle crossing respectively domains D14, D24, D23 and D13; one cycle crossing respectively domains D14, D24, D23, D22, D21, D11, D12 and D13 (Figure 7). The conditions on the parameter values of Model 3 to respect constraint (1) can then be directly derived.
https://doi.org/10.1371/journal.pone.0017075.s004
(DOC)
Text S2.
Estimation of the terms modeling the translocation of nuclear Mdm2 to the cytoplasm and the Mdm2-mediated acceleration of p53 degradation in the OAK Model.
https://doi.org/10.1371/journal.pone.0017075.s006
(DOC)
Text S4.
Comparison between the bifurcation diagrams of Model 1 and Model 2.
https://doi.org/10.1371/journal.pone.0017075.s008
(DOC)
Text S5.
Conditions on the parameter values of Model 3 to respect constraint (2).
https://doi.org/10.1371/journal.pone.0017075.s009
(DOC)
Text S6.
Analytical expression of the first return map for Model 3.
https://doi.org/10.1371/journal.pone.0017075.s010
(DOC)
Figure S1.
Bifurcation diagram and projection of the phase portrait of birhythmicity in the plane (Mn,P) for the OAK model. (A) Bifurcation diagrams of p53 level as a function of dMn for the OAK Model (Ouattara et al., 2010). Solid lines (resp. dashed lines) represent the stable (resp. unstable) equilibrium points. Bold (resp. white) dots are the maxima and minima of the stable (resp. unstable) limit cycles. The system shows a birhythmic domain for 1.09 h-1<dMn<1.3h−1. (B) Projection of the two oscillatory regimes on the plane (Mn,P) for dMn = 1.2 h−1. The thresholds, KMc and KMn, related to P are indicated in dashed lines. The parameter values for the bifurcation diagram are KMn = 0.1 nM, KMc = 0.6 nM, kP = 5 nM.h−1, kMc = 0.1 nM.h−1, k'Mc = 1.2 nM.h−1, dP = 0.1 h−1, d'P = 2.3 nM−1.h−1, dMc = 0.6 h−1, kin = 0.45 h−1, k'in = 0.4 h−1, kout = 0.045 h−1, KP = 0.2 nM and Vr = 10.
https://doi.org/10.1371/journal.pone.0017075.s011
(DOC)
Figure S2.
Study of the number of equilibrium points for Model 1. Graph of the rational function R (left) and a zoom of this graph (right) for the parameter values indicated in Figure 3. The roots of R give the equilibrium points for Model 1. The equilibrium points are included in the interval [0; ].
https://doi.org/10.1371/journal.pone.0017075.s012
(DOC)
Acknowledgments
The authors would like to thank Professor Marcelle Kaufman for critical reading of the manuscript.
Author Contributions
Conceived and designed the experiments: WA-J MC J-LG. Performed the experiments: WA-J. Analyzed the data: WA-J. Contributed reagents/materials/analysis tools: WA-J MC J-LG. Wrote the paper: WA-J MC J-LG.
References
- 1.
Goldbeter A (1996) Biochemical Oscillations and Cellular Rythms. The molecular basis of periodic and chaotic behavior. Cambridge University Press. 605 p.
- 2. Goldbeter A (2002) Computational approaches to cellular rhythms. Nature 420: 238–245.
- 3. Bar-Or RL, Maya R, Segel LA, Alon U, Levine AJ, et al. (2000) Generation of oscillations by the p53–mdm2 feedback loop: a theoretical and experimental study. Proc Natl Acad Sci USA 97: 11250–11255.
- 4. Decroly O, Goldbeter A (1982) Birhythmicity, chaos, and other patterns of temporal self-organization in a multiply regulated biochemical system. Proc Natl Acad Sci U S A 79: 6917–6921.
- 5. Alamgir M, Epstein I (1983) Birythmicity and compound oscillation in coupled chemical oscillators: chlorite-bromate-iodide system. J Am Chem Soc 105: 2500–2501.
- 6. Citri O, Epstein I (1988) Mechanistic study of a coupled chemical oscillator: the bromate-chlorite-iodide reaction. Journal of Physical Chemistry 92: 1865–1871.
- 7. Hounsgaard J, Hultborn H, Jespersen B, Kiehn O (1988) Bistability of alpha-motoneurones in the decerebrate cat and in the acute spinal cat after intravenous 5-hydroxytryptophan. J Physiol 405: 345–367.
- 8. Ventura A, Kirsch DG, McLaughlin ME, Tuveson DA, Grimm J, et al. (2007) Restoration of p53 function leads to tumour regression in vivo. Nature 445: 661–665.
- 9. Vogelstein B, Lane D, Levine AJ (2000) Surfing the p53 network. Nature 408: 307–310.
- 10. Vousden KH, Lane DP (2007) p53 in health and disease. Nat Rev Mol Cell Biol 8: 275–283.
- 11.
Geva-Zatorsky N, Rosenfeld N, Itzkovitz S, Milo R, Sigal A, et al. (2006) Oscillations and variability in the p53 system. Mol Syst Biol. 2. Available: http://www.nature.com/msb/journal/v2/n1/full/msb4100068.html. Accessed 13 June 2006.
- 12. Abou-Jaoudé W, Ouattara DA, Kaufman M (2009) From structure to dynamics: frequency tuning in the p53-Mdm2 network I. Logical approach. J Theor Biol 258: 561–577.
- 13. Ouattara DA, Abou-Jaoudé W, Kaufman M (2010) From structure to dynamics: frequency tuning in the p53-Mdm2 network. II Differential and stochastic approaches. J Theor Biol 264: 1177–1189.
- 14. Ciliberto A, Novak B, Tyson JJ (2005) Steady states and oscillations in the p53/mdm2 network. Cell Cycle 4: 488–493.
- 15. Brooks CL, Gu W (2006) p53 ubiquitination: Mdm2 and beyond. Mol Cell 21: 307–315.
- 16. Fang S, Jensen JP, Ludwig RL, Vousden KH, Weissman AM (2000) Mdm2 is a RING finger-dependent ubiquitin protein ligase for itself and p53. J Biol Chem 275: 8945–8951.
- 17. Inoue T, Geyer RK, Howard D, Yu ZK, Maki CG (2001) MDM2 can promote the ubiquitination, nuclear export, and degradation of p53 in the absence of direct binding. J Biol Chem 276: 45255–45260.
- 18. Chen J, Lin J, Levine AJ (1995) Regulation of transcription functions of the p53 tumor suppressor by the mdm-2 oncogene. Mol Med 1: 142–152.
- 19. Kruse J-P, Gu W (2009) Modes of p53 regulation. Cell 137: 609–622.
- 20. Oliner JD, Pietenpol JA, Thiagalingam S, Gyuris J, Kinzler KW, et al. (1993) Oncoprotein MDM2 conceals the activation domain of tumour suppressor p53. Nature 362: 857–860.
- 21. Barak Y, Juven T, Haffner R, Oren M (1993) Mdm2 expression is induced by wild type p53 activity. EMBO J 12: 461–468.
- 22. Freedman DA, Wu L, Levine AJ (1999) Functions of the MDM2 oncoprotein. Cell Mol Life Sci 55: 96–107.
- 23. Gottlieb TM, Leal JFM, Seger R, Taya Y, Oren M (2002) Cross-talk between Akt, p53 and Mdm2: possible implications for the regulation of apoptosis. Oncogene 21: 1299–1303.
- 24. Bakkenist CJ, Kastan MB (2003) DNA damage activates ATM through intermolecular autophosphorylation and dimer dissociation. Nature 421: 499–506.
- 25. Stommel JM, Wahl GM (2005) A new twist in the feedback loop: stress-activated MDM2 destabilization is required for p53 activation. Cell Cycle 4: 411–417.
- 26.
Hirsch MW, Smith HL (2003) Monotone systems, a mini-review, in Positive Systems. Proceedings of the First Multidisciplinary Symposium on Positive Systems (POSTA 2003). In: Benvenuti L, De Santis A, Farina L, editors. Lecture Notes on Control and Information Sciences vol 294. Heidelberg: Springer Verlag.
- 27. Hirsch MW (1988) Systems of differential equations which are competitive or cooperative: III. Competing species. Nonlinearity 1: 51–71.
- 28. Zhu H-R, Smith H (1994) Stable periodic orbits for a class of three dimensional competitive systems. J Diff Eqn 110: 143–156.
- 29. Glass L, Kauffman SA (1973) The logical analysis of continuous, non-linear biochemical control networks. J Theor Biol 39: 103–129.
- 30. Omholt SW, Kefang X, Andersen O, Plahte E (1998) Description and analysis of switchlike regulatory networks exemplified by a model of cellular iron homeostasis. J Theor Biol 195: 339–350.
- 31. de Jong H, Gouzé J-L, Hernandez C, Page M, Sari T, et al. (2004) Qualitative simulation of genetic regulatory networks using piecewise-linear models. Bull Math Biol 66: 301–340.
- 32. Ropers D, de Jong H, Page M, Schneider D, Geiselmann J (2006) Qualitative simulation of the carbon starvation response in Escherichia coli. Biosystems 84: 124–152.
- 33.
Dayarian A, Chaves M, Sontag ED, Sengupta AM (2009) Shape, size, and robustness: feasible regions in the parameter space of biochemical networks. PLoS Comput Biol. 5. Available:http://www.ploscompbiol.org/article/info%3Adoi%2F10.1371%2Fjournal.pcbi.1000256. Accessed 2 January 2009.
- 34. Glass L, Pasternack JS (1978) Stable oscillations in mathematical models of biological control systems. J Math Biol 6: 207–223.
- 35. Snoussi E (1989) Qualitative dynamics of piecewise-linear differential equations: a discrete mapping approach. Dyn Stabil Syst 4: 189–207.
- 36. Glass L (1975) Classification of biological networks by their qualitative dynamics. J Theor Biol 54: 85–107.
- 37. Edwards R (2000) Analysis of continuous-time switching networks. Physica D 146: 165–199.
- 38. Casey R, de Jong H, Gouzé , J-L (2006) Piecewise-linear models of genetic regulatory networks: equilibria and their stability. J Math Biol 52: 27–56.
- 39.
Strogatz H (1994) Nonlinear Dynamics and Chaos. With applications to Physics, Biology, Chemistry, and Engineering. Westview. 498 p.
- 40. Novick A, Weiner M (1957) Enzyme induction as an all-or-none phenomenon. Proc Natl Acad Sci U S A 43: 553–566.
- 41. Gardner TS, Cantor CR, Collins JJ (2000) Construction of a genetic toggle switch in Escherichia coli. Nature 403: 339–342.
- 42. Xiong W, Ferrell JE (2003) A positive-feedback-based bistable ‘memory module’ that governs a cell fate decision. Nature 426: 460–465.