A mathematical model of P53 gene regulatory networks under radiotherapy.

P53, a vital anticancer gene, controls the transcription and translation of a series of genes, and implement the cell cycle arrest and cell apoptosis by regulating their complicated signal pathways. Under radiotherapy, cell can trigger internal self-defense mechanisms in fighting against genome stresses induced by acute ion radiation (IR). To simulate the investigating of cellular responding acute IR at single cell level further, we propose a model of P53 gene regulatory networks under radiotherapy. Our model can successfully implement the kinetics of double strand breaks (DSBs) generating and their repair, ataxia telangiectasia mutated (ATM) activation, as well as P53-MDM2 feedback regulating. By comparing simulations under different IR dose, we can try to find the optimal strategy in controlling of IR dose and therapy time, and provide some theoretical analysis to obtain much better outcome of radiotherapy further.

Introduction

As one of the main tumor therapies, radiotherapy acts through the induction of DSBs to DNA, triggering the cellular self-defense mechanisms to induce apoptosis of cancerous cells via programmed apoptosis (Perez and Brady, 1998, Magné et al., 2006). Under continuous effect of acute IR, DSBs occur within the cell, and trigger the binding of the sensor protein kinase ATM to the nascent DNA ends, consequently increase the kinase activity of ATM (Li et al., 2001). As an important transcription factor within nuclear, P53 can be activated by DNA damage transferring through ATM activation. By regulating downstream genes and their signal pathways, active P53 further induce cell cycle arrest to repair DNA damage and cell apoptosis to eliminate the abnormal cells with genome damage or deregulated proliferation. Especially, these functional properties may modulate the activity of certain anticancer agents (Ritter et al., 2002, Pauklin et al., 2005).

The combined approaches of systems analysis, control theory, and computer science can stimulate new approaches to simulate the investigating of the cellular responding acute perturbations from outside (Dfichting et al., 1996). These approaches provide a good link between the diverging areas of biomedicine and mathematics (Chou, 2004, Chou and Zhou, 1982). Recently, some models have been proposed to simulate the kinetics of tumor therapy at cellular level (Ding et al., 2006). Especially, several models have been proposed to explain the damped oscillations of P53 in cell populations (Qi et al., in press-a, Qi et al., in press-b, Ma et al., 2005, Bar-Or et al., 2000). However, under radiotherapy, the cellular response mechanisms in fighting against genome stresses and the complicated regulations among vital genes need to be further addressed.

To simulate the investigations of the cellular self-defense mechanisms in response to DNA damage further, we propose a model of P53 gene regulatory networks under radiotherapy. By using differential equations, we implement the modules of DSBs generation and repair, ATM activation, and P53-MDM2 feedback regulation, as well as the dynamic interactions among these modules versus continuous radiation time. In our simulations, the dynamic regulations among DSB-protein complexes (DSBCs), ATM, P53, and MDM2, as well as the kinetics of toxins eliminating triggered by active P53 (P53#) are presented under continuous effect of acute IR, here we use the superscript # to represent the activate state. Meanwhile, the time threshold of ATM activation, and dynamic oscillations between P53 and MDM2, as well as the vital role of P53# in toxins eliminating are analyzed by comparing different simulation results.

Using differential equations and graphic approaches to study various dynamical and kinetic processes of biological systems can provide useful insights, as indicated by many previous studies on a series of important biological topics, such as enzyme-catalyzed reactions (Zhou and Deng, 1984, Chou and Zhou, 1982), hepatitis B viral infections (Xiao et al., 2005), and visual analysis of SARS-CoV (Wang et al., 2005). In this study, we are to use differential equations and graphic methods to investigate the dynamic and kinetic processes of P53 gene regulatory networks under radiotherapy.

Method

Model overview

Under radiotherapy, DNA damage trigger the cellular self-defense mechanisms to respond abnormal signals, and further eliminate the cancerous cells (Dewey et al., 1995). Based on the latest studies on biomedicine (Oren, 2003, Ritter et al., 2002) and former models (Qi et al., in press-a, Qi et al., in press-b, Ma et al., 2005, Pauklin et al., 2005), we propose a model of P53 gene regulatory networks under radiotherapy by using existing methods. Due to the complicated interactions among the related genes and their signal pathways under continuous IR (Dewey et al., 1995, Vogelstein et al., 2000), our model is a simplification of real biosystem in response to continuous perturbations on the basis of some assumptions. The integrated model scheme is shown in Fig. 1 .

Fig. 1

The integrated model scheme of P53 gene regulatory networks under radiotherapy. It includes the modules of DSBs generation and their repair, ATM activation, as well as P53-MDM2 feedback loop. As acute IR is applied, resulting DSBs are binded with repair proteins into DSBCs. ATM is activated by the transferring of DSBCs from damage sites, and subsequently trigger the regulatory mechanisms of P53-MDM2 feedback loop, controlling downstream genes and signal pathways to induce cell cycle arrest and cell apoptosis further.

The first part is the module of DSBs generating and their repair process. As acute IR dose is applied into a cell, DSBs stochastically occur in accordance with the initial strength of IR dose, and subsequently bind with repair proteins into DSBCs at damage sites (Rothkamm et al., 2003, Rothkamm and Lobrich, 2003, Kühne et al., 2004). The second part is ATM activation module. Followed the transferring of DSBCs signal from damage sites, ATM is activated rapidly with the cooperation of DSBCs and positive self-feedback of active ATM (ATM#). The module of P53-MDM2 feedback loop is the kernel part of P53 gene regulatory networks, which is sensitive to DNA damage transferring through ATM activation. P53 is activated rapidly followed the increase of ATM# within the cell (Bar-Or et al., 2000; Banin et al., 1998) and subsequently regulate its downstream genes and signal pathways to induce cell cycle arrest and cell apoptosis, controlling the process of DNA damage repair, and abnormal cells apoptosis (Weller, 1998, Ritter et al., 2002). The implementation of the three modules is described in the following paragraphs in detail.

DSBs generation and repair module

During the period of radiotherapy, DSBs occur and trigger two major repair mechanisms in eukaryotic cells: homologous recombination (HR) and nonhomologous end joining (NHEJ) (Rothkamm et al., 2003, Rapp and Greulich, 2004). DSB repair is a first-order process if break ends associated with the same DSB are rejoined, and a second-order process if the break ends associated with two different DSBs (Daboussi et al., 2002). Meanwhile about 60–80% of DSBs are rejoined quickly, whereas the remaining 20–40% of DSBs rejoined more slowly, with the precise relative contributions of slow and fast processes depending on the cell types (Daboussi et al., 2002, Budman and Chu, 2005).

To account for such dynamics under continuous IR, in first part of our model, we propose a module of DSBs generation and repair process based on previous studies (Ma et al., 2005). As shown in Fig. 2 , this module contains both a fast and a slow pathways with distinct reaction rates, each of which is composed of a reversible binding of DSB lesions with repair proteins into DSBCs, as well as an irreversible process from the DSBCs to fixed DSBs (Daboussi et al., 2002, Rapp and Greulich, 2004).

Fig. 2

The module of DSBs generating and their repair process, including a fast and a slow repair pathway. Resulting DSBs induced by acute IR can be in one of four states: intact DSB, DBSC, Fr, and Fw, subscripts ‘1’ and ‘2’ refer to the fast kinetics and the slow one, respectively.

Although DSBCs signal will subside whether the DSB was correctly repaired or not (Ma et al., 2005), the misrepair of DSBs have profound consequence for the subsequent viability of the cell, especially affect the outcome of radiotherapy (Dewey et al., 1995, Magné et al., 2006). Therefore, we obviously distinguish between correct repair part (F r) and misrepair part (F w) of DSBs, and further deal F w as a part of toxins remaining within the cell, which can be eliminated indirectly by the regulatory mechanisms of P53 gene regulatory networks during and after radiotherapy (Ritter et al., 2002, Pauklin et al., 2005).

Some experimental studies suggest that the quantity of the resulting DSBs within different IR dose domains obey a Poisson random distribution (Li et al., 2001). Given the different numbers of DSBs present at different IR doses domains, we assume that the stochastic number of the resulting DSBs per time scale is proportional to the number generated by a Poisson random function during the period of acute radiation (Ma et al., 2005). The DSBs generation process is formulated as follows:where [DT] is the total resulting DSBs in both fast and slow repair processes, k the parameter to set the number of DSBs per time scale, and a ir is the parameter to set the number of DSBs per IR dose.

In our model, we deal that 70% of the initial DSBs are fixed by fast repair kinetics, moreover, the adequate repair proteins are available around the damage sites. In each damage locus, DSB can be in one of four states including intact DSB, DSBC, F r, and F w. At each time step k, the concentrations of DSBs in four states are formulated as follows: where D(k), C(k), F r(k) and F w(k) represent the total concentrations of DSBs, DSBCs, F r, and F w in both fast and slow repair processes, subscripts ‘1’ and ‘2’ refer to the fast and the slow kinetics, respectively. The others formulations used in DSBs repair process are listed as follows: where [D], [C], and [F w] represent the respective concentrations of DSB, DSBC, and F w; k dc, k cd, and k fw are the transition rates among above three states of DNA damage; k dc and k cross represent the repair rates in the first-order and the second-order process (Daboussi et al., 2002).

ATM activation module

In our model, we deal that DSBCs are the main signal transferring from DSBs to the downstream P53-MDM2 feedback regulatory loop through ATM activation (Ma et al., 2005). As a DNA damage detector, ATM is extremely sensitive to DSBCs transferring from damage sites (Tribius et al., 2001, Pauklin et al., 2005). Under continuous effect of acute IR, the increasing DSBCs prompt the phosphorylation of inactive ATM monomers, as a result, the concentration of ATM# increase rapidly with the cooperation of DSBCs and the positive self-feedback from ATM# for intermolecular autophosphorylation (Pauklin et al., 2005; Cuddihy and Bristow, 2004). ATM activation act an important role in prompting the regulatory mechanisms of P53 gene networks further (Oren, 2003; Kohn and Pommier, 2005).

Shown in Fig. 3 is the module scheme of ATM activation under radiotherapy, which is built on the basis of the previous models (Qi et al., in press-a, Qi et al., in press-b, Ma et al., 2005). ATM monomer is activated under the cooperation of DSBC and self-feedback of ATM# (Banin et al., 1998, Tribius et al., 2001, Pauklin et al., 2005, Ma et al., 2005). In this module, the rate of ATM activation is assumed to be a function of the amount of DSBCs and ATM#, especially, the total concentration of ATM including ATM dimer, ATM monomer and active ATM is a constant within the cell (Ma et al., 2005). The related formulations used in the reversible kinetics of ATM activation are listed as follows: where [ATMd], [ATM] and [ATM#] represent the concentrations of ATM dimer, ATM monomer, and active ATM monomer, respectively; k undim the rate of ATM undimerization, and k dim the rate of ATM dimerization; k ar the rate of ATM monomer inactivation, and k af is the rate of ATM monomer activation. In addition, f is the function of ATM activation, the term a 1 C implies the fact that DSBs somehow activate ATM molecules at a distance, a 2[ATM#] indicates the mechanism of autophosphorylation of ATM, and a 3 C[ATM#] represents the interaction between the DSBCs and ATM# (Bakkenist and Kastan, 2003, Ma et al., 2005).

Fig. 3

The module scheme of ATM activation with the cooperation of DSBCs signal and self-feedback of ATM#, under continuous effect of acute IR, ATM is activated rapidly from ATM monomers triggered by DSBCs transferring.

P53-MDM2 feedback loop module

The module scheme of P53-MDM2 feedback loop is plotted in Fig. 4 . This negative feedback loop formed by P53 and MDM2 is the core part of P53 gene regulatory networks in response to genome stresses (Vogelstein et al., 2000, Ritter et al., 2002, Kohn and Pommier, 2005). Followed damage signal transferring, ATM# prompt the phosphorylation of P53, and further elevate the transcriptional activities of P53# (Vogelstein et al., 2000, Bar-Or et al., 2000, Ritter et al., 2002). As a result, MDM2, a P53-specific ligase and antagonist of P53 (Ma et al., 2005), is prompted by the increase of P53#. Reversely, P53 is depressed by the increase of MDM2 protein (MDM2p). This negative feedback loop can produce oscillations in response to the sufficiently strong IR dose (Ma et al., 2005).

Fig. 4

The module of P53-MDM2 feedback loop. P53 is translated from P53mRNA and phosphorylated by ATM#. P53# increase the transcribing of MDM2, whereas, MDM2 protein promote a fast degradation of P53 and a slow degradation of P53#. Meanwhile, ATM# stimulate the degradation of MDM2p, and increase the regulatory activation of P53# indirectly. Especially, the toxins are eliminated directly by P53# according to the assumptions.

To account for a decreased binding affinity between inactive P53 and P53#, we assume that MDM2p-induced degradation of inactive P53 is faster than that of P53#. Moreover, only P53# can prompt the activation of the target genes through their complicated regulation pathways, causing an overall effect of increased transcriptional activity (Pauklin et al., 2005, Ma et al., 2005). Especially, we deal that F w are the part of toxins accumulated within the cell, which can be eliminated directly by regulatory mechanisms of P53#. The main differential equations used in this module are as follows:

where [P53R], [P53P], [P53#], [MDM2R], [MDM2P], and [Toxins] represent the concentrations of P53 mRNA, P53 protein, active P53, MDM2 mRNA, MDM2 protein, and misrepair part of DSBs, respectively. SP53, and SMDM2 are the basal induction rates of P53 mRNA, and MDM2 mRNA. The other parameters are presented in Table A1, Table A2, Table A3 of Appendix A.

Table A1

The parameters used in the module of DSBs generation and repair process

Parameters	Description	Constant
kt	The rate of DSBs generation per time scale	0.01
air	The number of DSBs generation per IR dose	35
a1	The percentage of DSs processed by fast repair	0.70
a2	The percentage of DSs processed by slow repair	0.30
kdc1	The rate of DSBs transition to DSBCs in fast repair process	2
kdc2	The rate of DSBs transition to DSBCs in slow repair process	0.2
kcd1	The rate of DSBCs transition to DSBs in fast repair process	0.5
kcd2	The rate of DSBCs transition to DSBs in slow repair process	0.05
kfw1	The rate of DSBCs transition to Fw in fast repair process	0.001
kfw2	The rate of DSBCs transition to Fw in slow repair process	0.0001
kcross	The rate of DSB binary mismatch in second-order repair process	0.001

Table A2

The parameters used in the module of ATM activation process

Parameters	Description	Constant
kdim	ATM dimerization rate	8
kundim	ATM undimerization rate	1
Kaf	ATM phosphorylation rate	1
Kar	ATM dephosphorylation rate	3
a1	Scale of the activation function of ATM phosphorylation	1
a2	Scale of the activation function of ATM phosphorylation	0.08
a3	Scale of the activation function of ATM phosphorylation	0.8

Table A3

The parameters used in the module of P53-MDM2 feedback regulatory loop

Parameters	Description	Constant
SP53	Basal induction rate of P53 mRNA	0.001
drp	Degradation rate of P53 mRNA	0.02
krp	Translation rate of P53 mRNA	0.12
k#p	Dephosphorylation rate of P53#	0.2
kap#	ATM#-dependent phosphorylation rate of P53	0.6
kmp	MDM2p-dependent degradation rate of P53	0.1
km#	MDM2p-dependent degradation rate of P53#	0.02
dpp	Basal degradation rate of P53	0.02
dp#	Basal degradation rate of P53#	0.008
SMDM2	Basal induction rate of MDM2 mRNA	0.002
k#m	P53#-dependent MDM2 transcription rate	0.03
kmrp	Translation rate of MDM2 mRNA	0.02
dmr	Degradation rate of MDM2 mRNA	0.01
dmp	Basal degradation rate of MDM2p	0.003
kmat	ATM#-dependent degradation rate of MDM2p	0.01
ktfw	Toxins translation rate of misrepair of DSBs (Fw)	0.01
kp	Michaelis constant of ATM#-dependent P53# phosphorylation	1.0
k	Michaelis constant of P53#-dependent MDM2 mRNA transcription	1.0
kd	Threshold concentration for MDM2p-dependent P53 degradation	0.03
n	Hill coefficient of MDM2 transcription rate	4
kat	Threshold concentration for ATM#-dependent MDM2p degradation	1.60
kd#	Threshold concentration for MDM2p-dependent P53# degradation	0.32

Results

To ensure the rightness of the simulation results, we consider the fact that the valid parameter sets should obey the following rules (Qi et al., in press-a, Qi et al., in press-b, Bar-Or et al., 2000): (1) the model must contain oscillations. This is important as there has been experimental evidence that oscillations occur between P53 and MDM2 after cell stress; (2) the mechanism used to mathematically describe the degradation of P53 by MDM2 is accurate only for low concentrations of P53; (3) the concentration of P53# is much higher than that of inactive P53 after the system reaching an equilibrium. Based on above rules and the existing parameters used in the previous models (Qi et al., in press-a, Qi et al., in press-b; Ma et al., 2005; Bar-Or et al., 2000), we obtained the kinetics of cellular responding radiotherapy through simulation platform in MATLAB7.0. The detailed parameter sets used in our model can be found in Table A1, Table A2, Table A3 of Appendix A.

Dynamic process of DSBs generation and repair

During the simulation process, we applied 15 Gy IR to generate DSBs fraction. In order to agree with the experimental results that the measured 30–40 DSBs per Gy IR occurred in the single cell (Budman and Chu, 2005, Ma et al., 2005), the stochastic number of resulting DSBs were generated by using a Poisson random function with a mean of 35x as continuous IR dose of x Gy was applied. It means that 35  kx DSBs are generated on average within a single cell, x is IR dose, and k is a changeable parameter according to different time scale.

The total resulting DSBs (DT) are split into the simple breaks (D 1) in the fast repair process and complex ones (D 2) in the slow repair process (Lan, Ma et al., 2005). Initially, D 1(0) = 0.7DT, D 2(0) = 0.3DT and C 1(0) =  C 2(0) =  F r(0) =  F w(0) = 0. During DNA damage repair process, transitions occur among DSBs states with respective rates. In addition, there are adequate repair proteins available around damage sites. Shown in Fig. 5a is a stochastic trace of the resulting DSBs versus continuous radiation time.

Fig. 5

The kinetics of DSBs generating and their repair process under constant 15 Gy IR: (a) the stochastic trace of DSBs generation vs. radiation time; (b) the relationships among intact DSBs, DSBCs and Fw. DSBCs increase dramatically, intact DSBs and Fw are accumulate with respective rates.

Compared with the Monte Carlo methods (Ma et al., 2005), our simulation method is more suitable for the conditions that many resulting DSBs and adequate repair proteins exist in the cell. Meanwhile, the simulation results can describe the stochastic kinetics of DSBs generation under constant effects of the acute IR. Shown in Fig. 5b are the dynamic traces of resulting DSBs, DSBCs, and intact DSBs in response to 15 Gy IR in both fast and slow repair kinetics. Due to the adequate repair proteins available around damage sites and the irreversible step of the final repair process, DSBCs increase dramatically against radiotherapy time. Meanwhile, the remaining of the intact DSBs increase slowly due to the DSBs increasing and the limited ability of cellular repair process, and F w are accumulated within the cell without the elimination mechanisms of P53 gene regulatory networks.

The kinetics of ATM activation

ATM activation module was established to describe the switch-like dynamics of the ATM activation in response to DSBCs increasing. As a sensor of genome damage, ATM activation is important to relay DNA damage signal, and further trigger the regulatory functions of P53-MDM2 feedback loop (Banin et al., 1998, Bakkenist and Kastan, 2003, Pauklin et al., 2005). To simulate the mechanism of ATM activation under constant effect of 15 Gy IR, we illustrate a switch-like behavior of ATM plotted in Fig. 6 .

Fig. 6

The kinetics of ATM activation with the cooperation of DSBCs and self-feedback of ATM# under IR = 15 Gy. (a) The kinetics of ATM dimer, ATM monomer, and ATM# under the increase of DSBCs, ATM# switch to “on” state after about 7 min against the decreasing ATM dimer and ATM monomer. (b) The kinetics of ATM activation within 50 min, the concentration of ATM# reach saturation state, and keep constant concentration after about 17 min.

Under the cooperative functions of DSBCs and the positive self-feedback of ATM#, as shown in Fig. 6a, the concentration of ATM# increases fast against the decreasing of ATM dimer and ATM monomer, both of which reach the cross point behind the IR by about 7 min, simultaneously, ATM# switching to “on” state. Due to the increasing DSBCs and the constant quantity of ATM, Fig. 6b shows that the concentration of ATM# reaches the saturation state after about 17 min, reaching a dynamic equilibrium versus radiation time. The step-like traces in Fig. 6 suggest that ATM module can produce an on-off switching signal to the P53-MDM2 feedback loop.

The regulatory mechanisms of P53-MDM2 feedback loop

P53-MDM2 feedback loop is a vital part in regulating the downstream genes and their signal pathways to induce cell cycle arrest and cell apoptosis in fighting against genome stresses indirectly (Bar-Or et al., 2000, Ritter et al., 2002). In response to the input signal of ATM#, P53-MDM2 module generates one or more oscillations.

Shown in Fig. 7a are the dynamic traces of P53# and MDM2 protein under continuous application of 15 Gy IR from time 0. Upon the activation by ATM# and decreased degradation by MDM2, the total amount of P53 proteins increases quickly. Due to the P53-dependent induction of MDM2 transcription, the increase of MDM2 proteins is sufficiently large to lower the P53 level, which in turn reduces the amount of the MDM2 proteins. The oscillation pulses shown in Fig. 7a have a period of 400 min, and the phase difference between P53 and MDM2 is about 100 min. Moreover, the first pulse is slightly higher than the second, quite similar with the experimental observations (Lahav et al., 2004) and the simulation results (Ma et al., 2005). Meantime, based on the assumption in this module that toxins including F w can be eliminated by the regulation effect of P53# directly, we plot the different kinetics of F w elimination with and without the regulatory functions of P53#. As shown in Fig. 7b, with the eliminating effect of P53# the toxins keep low level with some oscillations against the toxins increasing dramatically without P53#.

Fig. 7

The kinetics of P53-MDM2 feedback loop in response acute 15 Gy IR. (a) The oscillating traces between P53# and MDM2p with the same frequency and different swing. (b) The different kinetics of Fw elimination with and without the regulating effect of P53#. (c) The oscillating kinetics of P53# and MDM2 within 1500 min, both of which reaching a dynamic equilibrium vs. continuous therapy time. (d) The oscillating kinetics of Fw elimination within 1500 min, the quantities of Fw keep low level due to the elimination functions of P53#.

Meanwhile, the oscillations between P53# and MDM2 shown in Fig. 7c get much weaker in response to constant regulation of ATM#, and then reach a dynamic equilibrium versus radiation time. It suggests that the self-defense mechanisms of P53 gene regulatory networks have been “on” state in response to the genome stress induced by acute IR dose. As a result, the kinetics of F w elimination shown in Fig. 7d suggest that toxins within the cell keep low level due to the eliminating mechanisms of P53 gene regulatory networks under radiotherapy.

Discussion and conclusion

Based on the previous models and existing methods, a mathematical model of P53 gene regulatory networks under radiotherapy was proposed. By using a set of differential equations, combined with the Poisson random function, the stochastic kinetics of DSBs generating was implemented under constant effect of acute IR dose. Meanwhile, we assume that there are adequate repair proteins available around damage sites. Therefore, the increasing DSBC signal trigger the switch-like ATM activation, and further prompt the regulatory functions of P53 gene networks in response to genome stresses. In addition, we deal that F w is a part of toxins accumulated within the cell during radiotherapy, and further assume that the toxins can be eliminated directly by the regulatory functions of P53 gene regulatory networks.

In our simulations, it is demonstrated according to our model that ATM exhibits a strong sensitivity and switch-like behavior in response to the number of DSBs, fully consistent with the necessary outcome to transfer the stress signal to the P53-MDM2 feedback loop and arrest the cell cycle in response to the acute IR (Qi et al., in press-a, Qi et al., in press-b, Bar-Or et al., 2000). Once the IR dose is sufficiently large, the P53-MDM2 feedback loop will produce oscillations. Especially, the number and amplitude of the oscillations are different according to the IR dose domains and radiation time (Qi et al., in press-a, Qi et al., in press-b, Bar-Or et al., 2000). Also, we can present the importance of P53# in response to the genome stresses under radiotherapy, and simulate the internal self-defense mechanisms to eliminate toxins further.

Due to the more complicated mechanisms of real biosystem under radiotherapy, it is necessary to further verify and improve our model based on the latest studies and biomedicine experiences. Our mode, although simplified, does provide a framework for the theoretical analysis of the mechanisms of P53 gene regulatory networks under radiotherapy.