Modeling amplified p53 responses under DNA-PK inhibition in DNA damage response

During DNA double strand breaks (DSBs) repair, coordinated activation of phosphatidylinositol 3-kinase (PI3K)-like kinases can activate p53 signaling pathway. Recent findings have identified novel interplays among these kinases demonstrating amplified first p53 pulses under DNA-PK inhibition. However, no theoretical model has been developed to characterize such dynamics. In current work, we modeled the prolonged p53 pulses with DNA-PK inhibitor. We could identify a dose-dependent increase in the first pulse amplitude and width. Meanwhile, weakened DNA-PK mediated ATM inhibition was insufficient to reproduce such dynamic behavior. Moreover, the information flow was shifted predominantly to the first pulse under DNA-PK inhibition. Furthermore, the amplified p53 responses were relatively robust. Taken together, our model can faithfully replicate amplified p53 responses under DNA-PK inhibition and provide insights into cell fate decision by manipulating p53 dynamics.


Model construction
The model consists of 11 species and 65 kinetic parameters. The schematic diagram is shown in Supplementary Figure 1A and the model equations are listed in Supplementary Table 2. For p53, MDM2 and WIP1, we incorporated both mRNA (italics) and protein. Total levels of class-IV phosphoinositide 3-kinase (PI3K)-related kinases (PIKK) named Ataxiatelangiectasia-mutated (ATM), ataxia telangiectasia and Rad3-related (ATR) and DNA-dependent protein kinase (DNA-PK) were set to be constant during the simulation [1]. Eq.1 and Eq.2 describe the production and degradation of p53 and mdm2 mRNA. The degradation rate of p53 (k 2 ) was estimated from [2]. The degradation rate constant for mdm2 (k 4 ) was adapted from [3]. Eq.2 further incorporated P53 induced mdm2 expression. It has been shown that P53 is present at mdm2 promoter even under unstressed conditions and mdm2 promoter is not tightly associated with nucleosomes [4,5]. According to 'anti-repression' model, P53 is constitutively active but repressed by MDM2 [4]. Therefore, free p53 may transactivate a subset of genes such as MDM2. Activated p53 by posttranslational modifications may further transcriptionally upregulate downstream targets such as mdm2 and wip1 [6]. Therefore, Eq.3 describes wip1 turnover. The first two terms in Eq.4 show translation and degradation for P53. The combined P53 degradation rates (basal and MDM2-induced) were adjusted to be consistent with experiments [7]. ATM, ATR and DNA-PKcs can all phosphorylate P53 at Ser15 and activate P53 [4,6]. For simplicity, DNA-PKcs was abbreviated as DNA-PK throughout the text. Previous reports also showed that WIP1 can dephosphorylate activated ATM, and phosphorylated forms of both P53 and MDM2 [8]. Eq.5 describes the conversion from P53 to activated forms and deactivation by WIP1. We also described MDM2 mediated degradation of activated p53 as previously suggested [9]. Eq. 6 shows mdm2 translation and MDM2 turnover. ATM could phosphorylate MDM2 at Ser395 and induce MDM2 self-degradation [10,11]. Therefore, we assumed that MDM2 can perform self-degradation once it is phosphorylated in Eq. 6 and Eq. 7. Eq. 8 simply describes wip1 translation and WIP1 turnover. The first term in Eq. 9 represents DNA damage induced ATM activation. The second term shows ATM trans-autoactivation. WIP1 mediated ATM dephosphorylation was incorporated as described previously [12]. Eq. 10 also covers DNA damage induced ATR activation. It has been demonstrated that ATR can promote its own phosphorylation and activation via either autophosphorylation [13] or CtIP mediated resection [14]. We here described these effects simply by an autophosphorylation and activation term in Eq. 10. Previous experiments show that DNA-PK may undergo ATM mediated phosphorylation or autophosphorylation at several clusters [15] although recent findings discriminate these effects in further details [16]. Therefore, the second and third terms in Eq. 11 characterize ATM catalyzed DNA-PK phosphorylation and autophosphorylation, respectively for simplicity purpose. Basal deactivation rates are incorporated for PIKKs.
Noticeably, we assumed that DNA double strand breaks (DSBs) can directly activate PI3K-like kinases and subsequently induce p53 activation by phosphorylation (Supplementary Figure 1A and  ⋅ − + ) since some reports implied that the ATM kinase activity is inhibited by DNA-PK [17,18]. Earlier events in DSB repair such as MRE11-RAD50-NBS1 (MRN) and mediator of DNA damage checkpoint protein 1 (MDC1) recruitment, phosphorylation of γH2AX foci, binding of the KU70 and KU80 heterodimers to sites of DNA breaks as well as detailed downstream repair events were not incorporated in our model for simplicity [19]. It has been well documented that ATM relays DNA damage signals through Chk2 while ATR and DNA-PK can coordinate activation of downstream effectors via Chk1 [20,21]. We did not consider Chk1 and Chk2 in our current model for simplicity owing to their linear signaling transduction properties. Inhibition of DNA-PK catalytic subunit (DNA-PKcs) by specific inhibitors was simplified as DNA-PK inhibition or DNA-PKi. The numbers in Supplementary Figure 1A denote the parameters listed in Supplementary Table 2. Notably, numbers corresponding to basal deactivation rate of activated ATM, ATR and DNA-PK (31, 34 and 38) were not shown in the diagram. We further assumed that the inhibitors can completely block the activity of associated kinases. For treatment with ATM inhibitor, five parameters (k 12 , k 21 , k 28 , k 29 and k 36 ) were simultaneously set to be zero. For ATR inhibitor, three parameters (k 12 , k 32 and k 33 ) were restricted to be zero. For DNA-PK inhibitors, four parameters k 17 , k 35 , k 36 and k 37 were taken to be zero. Transcriptional and translational delays are adapted from Ma et al. [9]. The kinetic equations are listed in Supplementary Table 1  Given relatively low number of repair proteins, the DNA repair process was implemented through stochastic simulation. In DSB repair module, each locus in which a DSB is created can be in one of three states corresponding to intact DSB (state 1), DSB in complex with repair proteins (state 2), and (correctly or incorrectly) fixed DSB (state 3) (Supplementary Figure 1). At time step k, the number of DSBs in corresponding states (i.e. 1, 2 or 3) are represented by D(k), C(k), and F(k), respectively. By using subscripts '1' and '2' to differentiate simple DSBs (fast kinetics) and complex DSBs (slow kinetics), we have D(k) = D 1 (k). + D 2 (k), C(k) = C 1 (k) + C 2 (k) and F(k) = F 1 (k) + F 2 (k). The total number of repair proteins (RP) is assumed to be 20 as previously indicated [9]. At any time, variable repair proteins are free to bind to DSBs. To implement step size control, we chose a relatively small step size (Δt = 0.2, which is smaller than the shortest τ in τ-leap method in our simulation). That means, even the smallest 'τ' can cover more than one Δt in τ-leap simulation (Supplementary Figure 1). We updated the DNA repair module by consecutive Δt (i.e. n•Δt, where n=τ mod Δt, 'mod' denotes modulus after division). To retain compatibility with τ-leap size, the last step size within the current 'τ' was set to be τ-n•Δt. We simply modified the fixation rate of double strand breaks (k fix1 ) as: k fix1 = k fix1' • (ATM * +ATR * +[DNA-PK] * )/3 as activated ATM, ATR and DNA-PK can participate in the DSB repair process by various ways. The Monte Carlo algorithm for the evolution of DSBs during time [0, tfinal] is as follows: 1. Set the initial conditions. The initial total number of DSBs (DSB T ) was generated by a Poisson distribution. Set t=0. We further set D 1 (0) =0.7•DSB T , D 2 (0) =0.3•DSB T , C 1 (0) = C 2 (0) =0 and F 1 (0) = F 2 (0) = 0. Given that initially all repair proteins are free, the simulation starts with total RP=20. There is no DSB until intrinsic DSBs are encountered. If a new DSB occurs, we assume that it is in state 1 (i.e. intact DSB).

Summarizing the binomial τ-leap stochastic simulation
For a system with N molecular species, the vector X=(X 1 (t),…, X N (t)) denotes the number of molecules X i at time t, 1≤i≤N (N = 11 in our model). Each species is involved in M chemical reaction channels (M = 36 in our model) represented by R j , j=1, 2, …, M with a propensity function a j . Then, a j dt denotes the probability that one reaction R j will happen in the time interval [t, t + dt] and the v ij is the stoichiometric coefficient of species S i in reaction R j .
In binomial τ-leap method, each of the M reactions is allowed to fire in a given order in [t, t + τ]. The maximum firing number k max (j) is determined by the limiting substrates, i.e. the species which would be completely depleted if the reaction went to completion. Identification of the limiting reactants is the critical step for implementation of the algorithm.
During the stochastic simulations, the parameters should be rescaled by a factor Ω. For kinetic parameters describing zero-order reaction (e.g. with unit μM• min -1 ), corresponding parameters in stochastic simulations should be multiplied by Ω. For second order reaction parameters (e.g. with unit μM -1 • min -1 ), the corresponding parameters should be divided by Ω. The first order parameters are unchanged (See Supplementary Table 2 for parameters). All reactants are multiplied by Ω. Ω = 6×10 4 in the model.
Briefly, the binomial τ-leap algorithm is represented as follows: (1) Obtaining the stoichiometric coefficients v ij , the initial number of all species X(0) and the rate constants.
(3) Repeat step 4-6 until the final time t final has been reached.
(4) Obtaining the propensity function a j (X(t)) using the current population X(t).
(5) Compute τ and update t = t + τ. τ is defined as follows: Where f denotes a coarse-grained factor and is set to be 10 4 as suggested [23].
(6) For j = 1 to M reactions, (a) Computing k max (j) = min(v ij <0)(int(X i ' /|v ij |)), where X' denotes the updated X before execution of any reaction step and int( ) is the largest number of integer function.
(b) Compute p = a j τ/ k max (j) , sample k j from the binomial distribution as defined below: BD denotes binomial distribution. If a j τ > k max (j) , we set The equations are formulated with ordinary differential equations and delay differential equations. The italics denote mRNA species.