Oscillations of the p53-Akt network: implications on cell survival and death.

Intracellular protein levels of p53 and MDM2 have been shown to oscillate in response to ionizing radiation (IR), but the physiological significance of these oscillations remains unclear. The p53-MDM2 negative feedback loop -- the putative cause of the oscillations -- is embedded in a network involving a mutual antagonism (or positive feedback loop) between p53 and AKT. We have shown earlier that this p53-AKT network predicts an all-or-none switching behavior between a pro-survival cellular state (low p53 and high AKT levels) and a pro-apoptotic state (high p53 and low AKT levels). Here, we show that upon exposure to IR, the p53-AKT network can also reproduce the experimentally observed p53 and MDM2 oscillations. The present work is based on the hypothesis that the physiological significance of the experimentally observed oscillations could be found in their role in regulating the switching behavior of the p53-AKT network between pro-survival and pro-apoptotic states. It is shown here that these oscillations are associated with a significant decrease in the threshold level of IR at which switching from a pro-survival to a pro-apoptotic state occurs. Moreover, oscillations in p53 protein levels induce higher levels of expression of p53-target genes compared to non-oscillatory p53, and thus influence cell-fate decisions between cell cycle arrest/DNA damage repair versus apoptosis.

Introduction

In response to DNA damage from exposure to ultraviolet (UV) or ionizing radiation (IR), damped oscillations in the levels of p53 and MDM2 proteins have been observed in populations of cells such as fibroblasts (human and mouse) and MCF7 breast cancer cells [1], [2]. Isolated single cells exhibit oscillations that are more sustained, sometimes lasting more than three days, depending on the persistence of DNA damage [3]. The frequencies of these oscillations generally increase with increasing IR intensity, and oscillation periods are typically in the range of 4 to 7 hrs; unlike this narrow range in the observed periods, the amplitudes of p53 oscillations are quite variable. Because it is the number of oscillations rather than their amplitudes that is observed to be dependent on the dose of IR, p53 oscillations have been described as ‘digital’. MDM2 oscillations are coupled with the p53 oscillations, with corresponding peaks separated by 1.5 to 2.5 hrs. Similar observations have also been observed in vivo in mouse intestine and spleen [4] where damped p53 pulses with periods of 4.5 to 6 hrs are reported.

The putative mechanistic origin of the oscillations is the negative feedback loop between p53 and MDM2 [1], [3], [5]–[9]. Transcriptional activity of p53 induces MDM2 expression [10], and MDM2 ubiquitinates p53 thereby marking the latter for degradation [11]–[17]. DNA damage, such as double strand breaks (DSBs), triggers a cascade of intracellular signaling pathways that lead to post-translational modifications (e.g., phosphorylation) of both p53 and MDM2, rendering p53 more stable with enhanced transcriptional activity [18]–[20] and MDM2 less stable with shortened half life [21]. The net result is that DNA damage-induced protein modifications prolong the time for MDM2 to cause p53 degradation; this time delay is necessary to generate p53-MDM2 oscillations [1], [3], [6], [9].

The p53-MDM2 loop is, of course, not isolated from other molecular interactions in a cell. For example, as depicted in Figure 1A, the p53-MDM2 loop is embedded in a larger network that includes AKT, a gene associated with cell survival signaling pathways. (There are several isoforms of AKT, but among these, we will specifically refer to AKT1 as AKT in this paper). The significance of this p53-AKT network (Figure 1A) is highlighted by the following facts: it involves two known tumor suppressor genes (p53 and PTEN) [22]–[26], and two oncogenes (MDM2 and AKT) [27]–[29]. Note that the p53-AKT network is composed of two feedback loops: the mutual antagonism between p53 and AKT which is a positive feedback loop (edges 1 to 5), and the p53-MDM2 negative loop (edges 5 and 6). The mutual antagonism between p53 and AKT is significant because p53 and AKT have opposite effects on apoptosis. p53 induces the expression of pro-apoptotic genes such as PUMA, BAX, NOXA, FAS and BAD [30]–[34], while AKT inhibits (by phosphorylating) pro-apoptotic proteins such as BAD [35], [36], CASPASE-9 [37] and FKHRL1 [38]. Note that BAD is a common target of p53 and AKT [39], which may further explain BAD's crucial role in apoptosis [40]. These considerations led us to hypothesize that the p53-AKT network is an important module that controls pathways to cell survival or death [41], [42].

Figure 1

(A) Qualitative network showing mutual antagonism between AKT and p53. An arrow means a pathway that leads to activation or upregulation; a hammerhead represents inhibition or downregulation. AKT is antagonized by p53 via edges 1–3, and p53 is antagonized by AKT via edges 4–5. A p53-MDM2 negative feedback loop is shown by edges 5 and 6. (B) The oscillatory p53-AKT model (the Model). The rate expressions vi's are given in Table 1. Note: mdm2 and pten denote mRNAs of MDM2 and PTEN gene, respectively. Mechanisms of the p53-AKT mutual antagonism: the network is activated by two cellular events namely growth factor stimulation and cellular stress. The former leads to phosphorylation of PIP2 (phosphatidyl inositol-4,5-bisphosphate) to PIP3 (Phosphatidyl inositol 3,4,5. triphosphate) and thereby PIP3-mediated activation of AKT [43]. Active AKT (AKTa) phosphorylates MDM2 leading to MDM2a translocation to the nucleus [44] where it inhibits p53, and thereby promotes cell survival [45]–[50]. Under cellular stress however, transcriptionally active p53 induces PTEN [51]–[54] among other target genes [55]. PTEN antagonizes AKT-mediated cell survival by deactivating PIP3 through efficient dephosphorylation [56]. p53-dependent transcription of PTEN upon irradiation-induced DNA damage has been reported in MCF-7 (human breast cancer epithelial), A172 (human glioblastoma), MEF (mouse embryonic fibroblasts) and tissues of mouse small intestine, colon, kidney and liver [57]–[59] where p53 oscillations are also observed [1], [3]–[5], [60].

Table 1

Kinetic rate equations and parameters of the Model.

Differential Equations			Kinetic Parameters
			Chosen values		Range
d[p53]/dt	=	v0 - v1 - v2	k0,basal = 0.02		0.005–0.2
d[AKTa]/dt	=	v8 - vm8	kd = 0.02		0.02–0.2
d[MDM2a]/dt	=	v6 - vQ6 - v9	k2 = 0.0625		0.0184–0.092
d[MDM2]/dt	=	v10 - v3 - v6+vQ6	j2 = 0.01		0.03–0.3
d[mdm2]/dt	=	v4+v5 - v7	k5 = 0.0375		0.024
d[PTEN]/dt	=	v14 - v15	j5 = 0.5		∼1
d[pten]/dt	=	v11+v12 - v13	k6 = 22		22 (Ref. [81])
d[PIP3]/dt	=	v16 - vm16	j6 = 0.6		0.5±0.1 (Ref. [81])
			kQ6 = 0.5		0.0000297–2.92
[AKT_Tot]	=	[AKT]+[AKTa]	jQ6 = 0.1		0.00238–2.23
[PIP_Tot]	=	[PIP2]+[PIP3]	k8 = 20		∼20
			j8 = 0.1		∼0.1
v0	=	k0	km8 = 0.2		0.0000297–2.92
v1	=	kd * [p53]	jm8 = 0.1		∼0.1
v2	=	k2 * [MDM2a] * [p53]/(j2+[p53])	k10 = 0.02	(New)	(Ref. [6])
v3	=	dMDM2 * [MDM2]	pmdm2 = 0.0009	(New)	(Ref. [6])
v4	=	pmdm2	dmdm2 = 0.01		(Ref. [6])
v5	=	k5 * [p53]n2/(j5n2+[p53]n2)	dMDM2 = 0.005		0.0028–0.0347
v6	=	k6 * [AKTa] * [MDM2]/(j6+[MDM2])	dMDM2a,basal = 0.005		0.0028–0.0347
vQ6	=	kQ6 * [MDM2a]/(jQ6+[MDM2a])	k11 = 0.006		∼0.006
v7	=	dmdm2 * [mdm2]	j11 = 2		>1
v8	=	k8 * [PIP3] * [AKT]/(j8+[AKT])	dPTEN = 0.0054		0.0025–0.0083
vm8	=	km8 * [AKTa]/(jm8+[AKTa])	ppten = 0.0009	(New)	–
v9	=	dMDM2a * [MDM2a]	dpten = 0.01	(New)	–
v10	=	k10 * [mdm2]	k14 = 0.02		–
v11	=	k11 * [p53]n1/(j11n1+[p53]n1)	k16 = 0.15		∼0.15
v12	=	ppten	j16 = 0.1		∼0.1
v13	=	dpten * [pten]	km16 = 73		73±4.4
v14	=	k14 * [pten]	jm16 = 0.5		0.1–1
v15	=	dPTEN * [PTEN]	n1, n2 = 4
v16	=	k16 * [PIP2]/(j16+[PIP2])	k0, IR = 0.0025	(New)
vm16	=	km16 * [PTEN] * [PIP3]/(jm16+[PIP3])	dMDM2a,IR = 0.0025	(New)
			[AKT_Tot] = 1
Note:			[PIP_Tot] = 1
k0	=	k0,basal+k0,IR *ρ * exp(-λ*t)
dMDM2a	=	dMDM2a,basal+dMDM2a,IR *ρ * exp(-λ*t)

New kinetic rate equations and parameters (relative to those used in Ref. [41]) are indicated. Time is in minutes; concentration is in µM.

The major issue addressed in this paper is the impact of the oscillatory dynamics of the p53-MDM2 loop on the antagonism between pro-apoptotic pathways (promoted by high p53 levels) and pro-survival pathways (promoted by high AKT actvity). We have shown earlier [41] that the p53-AKT network predicts an all-or-none (bistable) switching behavior between pro-survival and pro-apoptotic states, and will now show how this switch is affected by p53-MDM2 oscillations.

Results

Model Formulation

The model mechanism considered here is shown in Figure 1B; it will be referred to from hereon as the Model. It is an extension of a model we previously studied (model Q3 in Ref. [41]) which exhibits a bistable switch between pro-survival and pro-apoptotic cellular states. The said model extension involves two additional variables, namely, the mRNA transcripts of the MDM2 and PTEN genes; these transcripts are symbolized by mdm2 and pten, respectively, in Figure 1B. The inclusion of mdm2 introduces a time delay in the negative feedback loop between p53 and MDM2 proteins, and increases the potential for oscillations. The dynamics of the Model is described by the set of ordinary differential equations (ODEs) and parameter values that are listed in Table 1.

Post-translational modifications of p53 and MDM2 due to DNA damage are assumed to be reflected in an increase of the rate parameter k (for p53 synthesis and activation) and of the rate parameter d (for degradation of active MDM2). Another assumption made is that k and d are directly proportional to IR intensity (ρ), that is, k and d. The justification of this assumption is based on the many known pathways associated with IR-induced DNA damage that upregulate p53 and down-regulate MDM2 (see Figure S1 of Supplemental Data for details).

Steady States and Oscillations

The steady states of the Model are determined by equating all the right-hand sides of the ODEs in Table 1 to zero, and solving for the roots of the resulting system of algebraic equations (see ). The steady states of p53 and active MDM2 (MDM2a) as functions of ρ are shown in Figure 2. The Model exhibits multiplicity of steady states within a range of ρ (6 to 20.8 Gy – in Figure 2). Within this range of ρ, the high-p53 (correspondingly, low-MDM2a) steady states are all locally stable nodes, and the middle branch of steady states are all unstable saddle points. Decreasing from the right-knee of the steady state curves to ρ ∼16 Gy, the dynamics about the low-p53 (high-Mdm2a) branch of steady states is damped oscillatory (stable spiral). The steady states enveloped by the gray dotted curves in Figure 2 exhibit unstable spirals that lead to sustained oscillations (limit cycles). The peaks and troughs of these periodic oscillations are indicated by the dots above and below the steady states, respectively. It is interesting to note that the steady state at ρ = 0 Gy is a stable spiral (damped oscillations); this is reminiscent of reported experimental observations of damped oscillations in some cell populations that are not exposed to DNA damage-causing radiation [3].

Figure 2

Steady states of (A) p53 and (B) active MDM2 (MDM2a) as a function of the intensity of ionizing radiation, ρ.

Local stability of these states is indicated as either stable (black curve: stable nodes; dark gray curve: stable spirals) or unstable (broken black curve: saddle points; light gray curve: limit cycles). Limit cycle oscillations whose amplitudes are shown by the dots arise from the unstable steady states (light gray curve) in the lower (A) or upper (B) branch of the steady state curves.

Generally, we found that the frequencies of the oscillations in p53, mdm2, MDM2 and MDM2a increase with increasing ρ (Figure S2 of Supplemental Data), in accord with experimental observations [3]. Furthermore, the oscillation frequencies of p53, mdm2, MDM2 and MDM2a are identical. The Model predicts oscillation periods of 3.5 to 5.2 hours, well within the reported range of experimental values [3].

The Model also reproduces experimental measurements of the time separation or delay (1.2–1.9 hrs) between the peaks of MDM2a and p53 oscillations. Generally, these time delays are not sensitive to ρ (Figure S2 of the Supplemental Data). The Model also predicts that the amplitudes of the oscillations of active AKT, PIP3, pten mRNA and PTEN will be difficult to measure experimentally because the concentration differences between crests and troughs are less than 4% of the mean concentration (data not shown) – this could explain why oscillations in these molecules have not been reported so far.

Cells exposed to increasing IR intensities

Computer experiments using the Model were performed to simulate the behavior of a cell that is exposed to a pulse of IR with fixed intensity ρ. For each computer simulation, the initial cellular levels of the proteins and transcripts are those of the steady states of a cell that is not exposed to IR (that is, the steady state of the Model with ρ = 0 Gy). As we increase ρ in the range where limit cycle oscillations are exhibited, the system initially displays a high-amplitude oscillation that settles down to the unique limit cycle with smaller amplitude at each ρ. Initially large amplitude p53 oscillations are observed in experiments with single cells [3]. The peaks of these initial oscillations are shown as black squares in Figure 3, and superimposed with the steady-state curve of mdm2. This figure shows that the initial peaks increase with ρ and that there is a particular ρ (symbolized by ρ*) where the system crosses a boundary surface associated with the unstable steady states (the dotted middle branch of steady states) and gets attracted to the upper branch of stable steady states. Since ρ* is less than the value of ρ corresponding to the right-knee of the steady state curve, we refer to ρ* as an early-switching point.

Figure 3

The steady-state curve of mdm2 mRNA is plotted as a function of the intensity of ionizing radiation, ρ.

The black squares indicate the first peaks of the mdm2 oscillations when the Model is simulated at different values of ρ, with each simulation having identical initial conditions (equal to the steady states of the system at ρ  =  0 Gy). The vertical line at ρ = ρ* = 13.8 Gy marks the value where the peak of mdm2 oscillation is attracted to the upper branch of steady states.

Effect of cell-cell variation

In our computer simulations, cell-to-cell variability is reflected in the distribution of initial concentrations of all the molecular species in the Model. For every molecular species, the initial concentration is varied around its steady state concentration at ρ = 0 Gy, with deviations ranging from 0 up to 200% of the steady state (this range is divided equally into 100 intervals and then randomly permutated to generate 101 Latin sets of initial conditions in each Latin hypercube sample). A total of 10 independent Latin hypercube samplings are carried out. For each Latin set of initial concentrations, trajectories are computed for each ρ within the range where there are multiple steady states (6 to 20.8 Gy). The percentages of sets of initial concentrations that lead to high-p53 steady states are given in Figure 4. Note that 10–20% of initial conditions can lead to an early-switching point at ρ equal to 8 Gy. At 18 Gy about 90% of the cells are predicted by the Model to make this transition. Thus, switching from the oscillatory state to the high-p53 state can be induced by cell-cell variations in their initial concentrations at ρ that are significantly lower than that corresponding to the right-knee of the steady-state curve (see Figure 2). Nevertheless, the propensity for early switching to high-p53 states becomes more deterministic (less dependent on cell-cell variations) as ρ increases.

Figure 4

Percentage of Latin sets of initial concentrations leading to high-p53 states.

The initial concentrations of all molecular species are varied simultaneously, with deviations ranging from 0 up to 200% of the steady state concentration at ρ  =  0 Gy (range is divided equally into 100 intervals). Hypercube sampling, which involves random permutations of these intervals, is used to select 101 unique Latin sets of initial conditions. The percentage of these sets of initial conditions leading to high-p53 state at various selected ρ is computed. A total of 10 independent hypercube samplings are performed and the result of each hypercube sample is shown in the figure as a bar.

Effect of cell-type variation

Measured basal concentrations of p53 in seven different cell lines range from about 1×104 to 22×104 molecules per cell [6], suggesting a cell type-dependent rate of p53 synthesis, and therefore a significant variation in the value of the parameter k of the Model. In addition, several reports [61]–[63] have documented cell-type specific post-translational modifications of various DNA-binding domains of p53 – suggesting that the model parameter j (associated with p53's affinity to the promoters of target genes) is also subject to wide variations depending on cell types. We investigated how these cell-type specific differences affect the behavior of the system by varying the parameters k and j. For each combination of k and j, a steady-state curve is computed. We sampled 4896 combinations of k and j , and determined their corresponding steady-state curves. All steady-state curves with multiple steady states exhibit either limit cycles or damped oscillations, or both, at the low-p53 steady state branch. Altogether, three general types of bifurcation curves or switching behavior are obtained among the 4896 parameter combinations (Figure S3 of the Supplemental Data):

Monostable (observed in 9% of the steady state curves). The only state of the system is a high-p53 stable node and no oscillation is manifested. This occurs at relatively high p53 production rate and fast dissociation rate of p53 from promoter site.

Early switch (observed in 87% of the steady state curves). This is the predominant type observed.

Saddle-node switch (observed in 4% of the steady state curves). For this type, the system switches to high-p53 state exactly at the right-knee of the curve.

Thus, the early switching phenomenon is likely to arise in cell types that exhibit p53 oscillations.

Consequences of p53 oscillations on expression of target genes

Since p53 is a transcription factor with many target genes [64], it is of interest to determine what the effects of p53 oscillations are on the expression levels of these genes. As an example, we considered a simplistic model for the expression of a representative p53-target gene X. This model includes the rate equations v11, v13, v14 and v15 shown in Figure 1B. The time-course of an oscillatory p53 is estimated by , where M is the mean p53 level between the crest and trough, A is the oscillation amplitude, and P is the oscillation period. The values of the parameters M, A, and P depend on ρ, and these values are estimated from the oscillations generated in Figure 2. The other rate expressions and all the kinetic parameters of this simplistic model are identical to those used in the Model (Figure 1B).

Figure 5A shows examples of temporal profiles of Xp (protein of p53-target gene X) expression for both non-oscillatory and oscillatory p53. As clearly shown, the rate of Xp synthesis is significantly increased when p53 oscillates. Figure 5B shows that the steady state of Xp is increased due to p53 oscillations for the entire range of ρ where limit cycles are manifested, despite the fact that the mean level of p53 oscillations is equal to that of non-oscillatory p53. We observe that the steady state level of Xp depends, not on the period P of p53 oscillations (data not shown), but on the mean (M) and amplitude (A) of the oscillations. The steady state levels of Xp increase with A, as depicted in Figure 5C.

Figure 5

Consequences of p53 oscillations on gene expression.

(A) Time courses of Xp expression induced by both non-oscillatory (gray) and oscillatory p53 (black) shown in the upper panel (this oscillation profile of p53 is similar to the limit cycle of the Model at ρ  = 12Gy). The level of the non-oscillatory p53 is equal to the mean of the oscillating values of p53. (B) Steady state level of Xp induced by both non-oscillating (gray curve) and oscillating (black curve) p53 for the entire range of ρ where limit cycle exists in the Model. The black curve is obtained by taking the average of the crest and trough of Xp at steady state. Similarly, the oscillation profiles of p53 at each ρ are estimated from the limit cycle of the Model. The level of non-oscillatory p53 is equal to the mean of the oscillatory p53. (C) Steady state level of Xp induced as a function of p53 oscillation amplitude. Parameter values used in the simulations (see text): the mean (M) and period (P) of oscillatory p53 are fixed at 0.5 µM and 4 hr, respectively. The trend shown in (C) is similar to those obtained using other values of p53 oscillation amplitudes (A = 0 to 0.5 µM), means (M = 0.1 to 2 µM) and periods (P = 4 to 7 hrs).

Many examples of the target protein Xp can be cited from the literature. For example, the cyclin-dependent kinase inhibitor p21 is a transcriptional target of p53 and has been observed to oscillate in response to p53-MDM2 oscillations upon irradiation in a wild-type human cell population [2]. In addition, cDNA microarray data analysis of p53 target genes revealed 51 genes whose mRNA levels oscillate [65].

To investigate whether p53 oscillations could also induce higher levels of expression of target genes with slow expression kinetics, simulations similar to Figure 5B are performed for each of the following four cases: (i) p53 binding affinity to gene promoter is decreased 2.5 fold, (ii) p53 induction of genes occurs 8 hrs post-irradiation, (iii) p53 transcriptional rate is decreased 6 fold, and (iv) target protein degradation rate is decreased 4 fold. Case (i) assumes that p53 binding to a promoter affects the timing of gene transcription. Lower p53 binding affinity to promoters has been invoked by others [66], [67] to explain the rapid induction of cell cycle genes (within 2 hrs) compared to the intermediate-to-late induction of apoptotic genes (4 to 8 hrs) post-irradiation. Case (ii) considers the alternative mechanism whereby late-responsive genes require p53 to recruit specific transcriptional co-factors for their induction, as reviewed recently [67]; the time required to recruit and upregulate the specific co-factors all contribute to the delayed expressions of late p53-responsive genes. Case (iii) accounts for the scenario where a slow rate of transcription leads to corresponding slow protein expression kinetics, and case (iv) considers proteins that take a comparatively longer time to attain steady state upon induction of expression because of their slow degradation rates. Simulation results show that oscillatory p53 induces higher level of gene expressions in all four possible causes of slow dynamics of gene expression (Figure S4 of the Supplemental Data).

p53-transcriptional regulation of apoptosis

Here, we present a model that includes explicitly the downstream signaling pathways from p53 and AKT to apoptosis, and study the consequences of p53 oscillations on the regulation of caspase-dependent apoptosis. We focus on the mitochondrial pathway because excessive DNA damage induces apoptosis via the mitochondrial pathway. The role of the BCL2 protein family in regulating this intrinsic pathway is well established. Some members of this family are anti-apoptotic while others are referred to as pro-apoptotic because they bind and inhibit the action of the anti-apoptotic members. DNA damage increases the rate of p53-dependent transcription of pro-apoptotic members such as BAX [31] and BAD [34]. Examples of anti-apoptotic members are BCL-2 [68] and BCL-XL [40], [69]. Excessive downregulation of BCL-2 and BCL-XL triggers caspase-dependent apoptosis [69], [70]. Furthermore, the ratio of the levels of BAX to BCL-2 can be a determinant of apoptosis, as has been illustrated in both experimental and modeling studies [71], [72]. AKT enters the picture because it phosphorylates BAD, thereby inhibiting the binding of BAD with BCL-XL [35], [36], and consequently promoting cell survival. The p53-AKT Model is modified to include the aforementioned steps as shown in Figure 6 and Table 2 (the resulting network will be referred to as the Apoptotic Model). The following assumptions are made: (1) p53-dependent transcription of bax and bad is assumed to decay exponentially with the same timescale λ as the rate of DNA damage repair. This is because phosphorylation of p53, especially at Serine 46 by DNA damage signal transducers such as ATM and DNA-PK, is needed to express pro-apoptotic proteins (see Figure S1 of the Supplemental Data); (2) As p53 first induces the expression of cell-cycle arrest and DNA damage repair genes, followed by pro-apoptotic genes (see [66], [67], [82], [83]), a Hill-type function of time after irradiation is included as a factor in both v and v. This factor incorporates a time delay for the expressions of bax and bad, following the expressions of mdm2 and pten; (3) The rate equations and kinetic parameters values associated with the reactions of bax, BAX and BCL-2 are set to be identical to those of bad, BAD and BCL-XL, respectively. Note that Table S1 (in Supplemental Data) further summarizes and compares the kinetic parameter values used in the Model (Figure 1) and the Apoptotic Model (Figure 6).

Figure 6

Network model of the regulation of pro-apoptotic genes by p53 and AKT.

Table 2

Kinetic rate equations and parameters of the model depicted in Figure 6.

Differential Equations			Kinetic Parameters
d[bax]/dt	=	v17 - v18	k17, k24 = 0.006
d[BAX]/dt	=	v19 - v20 - v23	j17, j24 = 2
d[BCL-2]/dt	=	v21 - v22 - v23	km17, km24 = 400
d[bad]/dt	=	v24 - v25	k19, k26 = 0.04
d[BAD]/dt	=	v26 - v27 - v30 - v31+vm31	k21, k28 = 0.00048
d[BADp]/dt	=	v31 - vm31 - v32	k23, k30 = 600
d[BCL-XL]/dt	=	v28 - v29 - v30	k31 = 44
			j31 = 0.01
v17	=	{k17 * [p53]n3/(j17n3+[p53]n3)} * {t n3/(km17n3+tn3)} *exp(-λ*t)	km31 = 0.01
v18	=	dbax * [bax]	jm31 = 10
v19	=	k19 * [bax]	dbax, dbad = 0.001
v20	=	dBAX * [BAX]	dBAX, dBAD, dBADp = 0.00054
v21	=	k21	dBCL-2, dBCL-XL = 0.0036
v22	=	dBCL-2 * [BCL-2]	n3, n4 = 4
v23	=	k23 * [BAX] * [BCL-2]
v24	=	{k24 * [p53]n4/(j24n4+[p53]n4)} * {t n4/(km24n4+tn4)} *exp(-λ*t)
v25	=	dbad * [bad]
v26	=	k26 * [bad]
v27	=	dBAD * [BAD]
v28	=	k28
v29	=	dBCL-XL * [BCL-XL]
v30	=	k30 * [BAD] * [BCL-XL]
v31	=	k31 * [AKTa] * [BAD]/(j31+[BAD])
vm31	=	km31 * [BADp]/(jm31+[BADp])
v32	=	dBADp * [BADp]

The values of the kinetic parameters associated with the reactions of bax, BAX and BCL-2 are assumed to be identical to those of bad, BAD and BCL-XL, respectively. Time is in minutes; concentration is in µM.

In the computer simulations, the presence of DNA damage is reflected in the values of the parameters k and d. Furthermore, DNA damage is assumed to decay exponentially with timescale λ - this assumption is based on the observation that DNA damage repair reduces the number of DSBs exponentially after exposure to IR doses ranging from 0.2 to 80 Gy [73]. Therefore, both k and d are assumed to decay exponentially to their respective basal levels according to the following expressions: k and d.

Representative values of λ are used in the simulations: λ = 0, 0.0005, and 0.001 (corresponding, respectively, to no repair, slow, and moderate repair rates – as inferred from experiments [73]). For all simulations, the initial state of the Apoptotic Model is set to its steady state at ρ = 0 Gy. In accordance with experimental observations [69]–[72], protein levels of the anti-apoptotic BCL-2 and BCL-XL are used as markers for the progress of apoptosis. For the purpose of exploring the qualitative behavior of the Apoptotic Model, and since there are no reported quantitative measurements, we arbitrarily assumed that a cell commits to apoptosis upon the decrease of either BCL-2 or BCL-XL to less than the chosen apoptotic threshold of 1% of its steady state level at ρ = 0 Gy. The number of p53 pulses needed to deplete BCL-2 below the apoptotic threshold level is determined. As shown in Figure 7, fewer p53 pulses are required as ρ is increased, and in range of ρ ∼8 Gy to ∼13 Gy the number of p53 pulses needed are more or less constant. The sharp decrease in the number of p53 pulses for ρ >14 Gy is due to the switch to the upper branch of p53 steady states. Interestingly, for fixed ρ, more p53 pulses are needed as the repair rate is increased.

Figure 7

Number of p53 pulses (rounded up to the nearest half pulse, vertical axis) required to deplete BCL-2 below the defined apoptotic threshold level (see text) as function of IR intensity, ρ, at various DNA damage repair rates λ.

Note that the number of p53 pulses becomes zero after the early-switching point.

As shown in Figure 7, a wide range of ρ generates 3 to 4 pulses of p53. To investigate whether this small number of pulses is sufficient to induce higher levels of target gene expression, protein levels of a p53 target gene (Xp) induced by oscillatory and non-oscillatory p53 are compared at 1200 min (this time period corresponds to 3 to 4 p53 pulses). Indeed, our simulations show that 3 to 4 pulses of p53 are sufficient to induce higher levels of expression of both early and late responsive target genes for the entire range of ρ where limit cycles are manifested (see Figure S5 of the Supplemental Data). These results hold for all time points.

Discussion

Our earlier modeling work [41] demonstrated that the p53-AKT network (Figure 1A) has the potential to exhibit switching behavior between pro-survival (high AKT, low p53) and pro-death (high p53, low AKT) states. Originating from the mutual antagonism between p53 and AKT, the switching dynamics is associated with transitions between two stable steady states (bistability) that coexist under the same parameters [41]. Within the p53-AKT network is a negative feedback loop between p53 and MDM2 (steps 5 and 6 in Figure 1A). Kinetic models focusing on this p53-MDM2 loop have been proposed [1], [3], [6]–[9] to explain the experimentally observed oscillations of p53 and MDM2 [1]–[6]; however, the physiological importance of these oscillations remains unclear. The present work is based on the hypothesis that the physiological significance of these p53 oscillations could be found in their role in regulating the switching behavior of the p53-AKT network between pro-survival and pro-death states.

We have shown in this paper that the p53-AKT network (the Model shown in Figure 1B), in addition to its potential for bistable behavior, is capable of generating sustained oscillations in p53 and MDM2, and that these oscillations are predicted to occur only within a range of IR intensities (ρ) and only in the low-p53 (pro-survival) branch of steady states (Figure 2). Using model parameters that are biologically plausible, the Model predicts oscillation periods of 3.5 to 5.2 hrs – well within the reported range of experimental values [3]. We also found that the oscillation frequencies generally increase with increasing ρ (Figure S2 of Supplemental Data), in accord with experimental observations [3]. In addition, the Model reproduces the experimentally measured time separation (or delay) between peaks of p53 and MDM2 oscillations (1.2–1.9 hrs), and, furthermore, predicts that these time delays are insensitive to ρ (Figure S2 of the Supplemental Data).

Zhang, Brazhnik and Tyson (2007) [8] have suggested models in which the p53-MDM2 negative feedback loop is coupled with a positive loop; in particular, one of these models exhibits multiple steady states as well as large-amplitude sustained oscillations around the high-p53 steady states outside the multiple steady-state regime. These oscillations emerge out of a homoclinic bifurcation [8]. In contrast, our p53-AKT model exhibits sustained oscillations around the low-p53 steady state in both single-steady state and multiple-steady state ranges of parameters (Figure 2); furthermore, these oscillations arise out of supercritical Hopf bifurcation. Two key experimental observations [3] support the Hopf bifurcation scenario: first, observed p53 oscillations show more variability in amplitudes than in periods, and, second, the initial p53 pulse is larger than subsequent pulses, indicating that the initial response to DNA damage is not sluggish (which would be the case if it were a homoclinic bifurcation [74]). Simulations using our model do show that the initial response to DNA damage agrees with the second experimental observation (see also Figure S6 of the Supplemental Data); the Model generates p53 pulses with larger amplitudes, longer periods and larger time delays between peaks of MDM2 and p53 compared to those generated by a stand-alone p53-MDM2 oscillator (Figure S7 of the Supplemental Data). In addition, the Model coupled with downstream pathways to apoptosis predict that there is a range of IR intensities where the number of p53 pulses acts as the digital count for the induction of apoptosis (Figure 7).

A subtle prediction of the Model comes from the computer-simulated experiments of exposing cells to increasing radiation intensity – subtle because it is not readily evident from the steady state bifurcation diagrams (Figures 2 and 3). In these simulations, the initial condition of the system before IR exposure is always identical to that of the case where ρ = 0 (let this initial condition be symbolized by which is a vector in 8-dimensional concentration space). For values of the parameter ρ in between those corresponding to the knees of the steady-state curve in Figure 3, the unstable middle steady state acts as a boundary between the top steady state and oscillatory state. Because the system is multi-dimensional, the boundary that separates the basins of attractions of the top steady state and oscillatory state (for ρ<16, see Figure 3) is a multidimensional surface – let us call this surface . The position of depends on ρ, and as ρ is increased, the position of relative to moves from one that is located in the basin of attraction of the oscillatory state to one on the other side of that is within the basin of attraction of the top steady state; this switch occurs at ρ*, the early-switching point (so-called because ρ* is less than the ρ corresponding to the right knee of the steady state curve). A cartoon that illustrates the aforementioned ideas is shown in Figure 8. A scattering of initial conditions around leads to various probabilities of switching from the oscillatory state to the top steady state; for example, as shown in Figure 4, about half of all sampled initial conditions around are attracted to the top steady state at ρ = ρ* = 13.8 Gy (Figure 4) which is significantly lower than the value of 20.8 Gy corresponding to the right-knee of the steady-state curve (Figure 2). These simulations suggest that because of the instability of the lower steady states where oscillations exist (Figures 2A and 3), the fluctuations in p53 and MDM2 increase the probability of crossing – causing the abovementioned early switch which can be interpreted as a decrease in tolerance to death-inducing IR.

Figure 8

A cartoon of the p53-MDM2a phase portrait for the cases ρ<ρ* and ρ>ρ*, where ρ* is the early-switching point.

(A) For ρ<ρ*, the initial state corresponding to ρ = 0 Gy (indicated by x0) is located within the basin of attraction of the limit cycle. (B) At ρ>ρ*, x0 is now located in the basin of attraction of the high-p53 stable node (•). Note: the actual phase portrait is high-dimensional; the two-dimensional portraits above are given merely to illustrate the basic idea of the origin of the early-switching point.

Another important result presented in this paper is the consequence of p53 oscillations on the level of expression of p53-target genes. It is predicted that at a fixed ρ where oscillations exist, p53-target genes are expressed at higher levels compared to the case where there are no oscillations (Figure 5B). Furthermore, this increase in expression of target genes is more pronounced as the amplitude of p53 oscillations becomes larger (Figure 5C). In the context of the p53-AKT network (Figure 1A), the model predicts an increase in PTEN expression (thereby promoting the antagonism against the survival factor AKT) and an increase in BAX and BAD expression (thereby promoting apoptosis). Thus, IR exposure that induces p53 oscillations is predicted to sensitize cells towards apoptotic death. The issue of how p53 oscillations may influence the determination of cell fate (i.e., between cell cycle arrest/DNA damage repair versus apoptosis) is not considered in the present work.

Reported differential sensitivity of mice organs to IR-induced apoptosis [75] suggests that IR sensitivity is closely linked with p53 oscillations and the proliferation status of cells. The early-switch phenomenon and the increased expression of p53-target genes due to p53 oscillations, as discussed in the preceding paragraph, may be offered as an explanation for this IR-sensitivity. In terms of the proliferation status of cells, highly IR-sensitive organs such as spleen and small intestine – in which p53 oscillations have been observed [4] – are generally composed of proliferating cells, whereas non-proliferating cells in organs such as the brain are generally IR-insensitive. Moreover, in vitro studies of p53 oscillations have so far been carried out using proliferating cells such as breast cancer epithelial [1], [3], [5], [6] and fibroblast cells [2]. The fact that our Model predicts that the oscillations occur around the low-p53 steady states is consistent with these observed association between proliferating cells and p53 oscillations.

We would like to offer ways to test the predictions of our model in the laboratory. Several single cell time-lapse microscopy experiments (e.g., Ref. [3]) could be performed to validate key hypotheses generated from the models by taking advantage of the evidence that not every cell (of the same type) would manifest p53-MDM2 oscillations when irradiated with the same IR intensity. For instance, only 40% of MCF-7 cells showed oscillations upon 10 Gy of gamma radiation [3]. Strictly speaking, a cell could only be in one of the three states: low-p53 (non-oscillatory cells), high-p53 and oscillatory p53 (oscillatory cells). The key prediction that p53 oscillations induce higher level of target gene expression could be tested by semi-quantifying the expression level of a (transfected) luciferase reporter gene that possesses a p53 promoter sequence upon irradiation. Oscillatory cells are predicted to express higher intensity of fluorescence than non-oscillatory cells (of the same type); comparisons should be made among cells that expressed similar mean level of p53. Also, an interesting experiment would be to test whether higher p53-dependent expression of cell cycle and DNA repair genes in oscillatory cells could lead to faster cell cycle arrest and repair damage DNA than non-oscillatory cells upon irradiation. In contrast, the other key prediction that p53 oscillations lower the IR intensity level at which the system switches to high-p53 state is relatively trickier to perform experimentally. It involves the determination of cumulative IR dose that lead to a high-p53 state in each oscillatory and non-oscillatory cells by increasing the IR dose gradually. Oscillatory cells are predicted to switch to high-p53 state over a wider range of cumulative IR dose with lower median than non-oscillatory cells.

Lastly, for the future development of our model, we would like to point out that besides PTEN, the insulin growth factor-binding protein 3 (IGFBP3) connects p53 to AKT. Upon DNA damage, IGFBP3 is upregulated by both p53-dependent and independent transcription [76], [77]. IGFBP3 binds and sequesters IGFs (insulin growth factors) away from IGFRs (insulin growth factor receptors), and thereby inhibits AKT activation; active IGF-bound IGFRs induce the downstream activation of the PI3K/AKT survival pathway [78]. Surprisingly, through unknown mechanisms, IGFBP3 could also sensitize cells to the phosphorylation of AKT by IGFs, which leads to AKT activation [79], [80]. Thus, further experimental studies are needed to resolve the conflicting relationship between p53, IGFBP3 and AKT.

Materials and Methods

To determine the steady states of the Model, the left-hand sides of the ODEs in Table 1 are all set to zero and the corresponding systems of nonlinear algebraic equations were solved numerically using Maple (version 7.0). The steady states as functions of certain parameters are referred to as steady-state bifurcation diagrams. The local stability of the steady states is determined using standard linear stability analysis, involving the determination of the eigenvalues of associated Jacobian matrices. To obtain the time-series trajectories of the species concentrations in the model, the ODEs are integrated using a modified Rosenbrock formula of order 2 that is implemented in the MATLAB (The MathWorks, Natick, MA) platform (version 6.5, Release 13).

(0.13 MB DOC)

Click here for additional data file.

(0.05 MB DOC)

Click here for additional data file.

(0.03 MB DOC)

Click here for additional data file.

(0.09 MB DOC)

Click here for additional data file.

(0.08 MB DOC)

Click here for additional data file.

(0.08 MB DOC)

Click here for additional data file.

(0.08 MB DOC)

Click here for additional data file.

(0.04 MB DOC)

Click here for additional data file.