Experimental design , validation and computational modeling uncover DNA damage sensing by DNA-PK and ATM †

Reliable and efficient detection of DNA damage constitutes a vital capability of human cells to maintain genome stability. Following DNA damage, the histone variant H2AX becomes rapidly phosphorylated by the DNA damage response kinases DNA-PKcs and ATM. H2AX phosphorylation plays a central role in signal amplification leading to chromatin remodeling and DNA repair initiation. The contribution of DNA-PKcs and ATM to H2AX phosphorylation is however puzzling. Although ATM is required, DNA-PKcs can substitute for it. Here we analyze the interplay between DNA-PKcs and ATM with a computational model derived by an iterative workflow: switching between experimental design, experiment and model analysis, we generated an extensive set of time-resolved data and identified a conclusive dynamic signaling model out of several alternatives. Our work shows that DNA-PKcs and ATM enforce a biphasic H2AX phosphorylation. DNA-PKcs can be associated to the initial, and ATM to the succeeding phosphorylation phase of H2AX resulting into a signal persistence detection function for reliable damage sensing. Further, our model predictions emphasize that DNA-PKcs inhibition significantly delays H2AX phosphorylation and associated DNA repair initiation.


Introduction
Cells are constantly affected by DNA damage, resulting from ionizing g-irradiation (IR), genotoxic or replication stress and reactive oxygen species.DNA damage, including single and double strand breaks (DSB), base modification, deletions or point mutations, seriously affects genome stability and cell integrity if not properly detected and repaired by the DNA damage response (DDR). 1 Upon DNA damage detection, higher order chromatin has to be made accessible by various modifications before DSB can be repaired. 2Among several DNA-damage associated histone modifications, phosphorylation of H2AX is widely accepted as an indicator of DSB.H2AX becomes rapidly phosphorylated at serine 139 (gH2AX) to generate foci at the DSB site. 3The assembly of chromatin remodeling complexes at the DSB site greatly depends on gH2AX and enables the accessibility of the damaged DNA to repair proteins. 4epending on the stimulus, gH2AX is induced by different members of the phosphoinositide 3-kinase like kinase (PIKK) family; ataxia telangiectasia mutated (ATM), ataxia telangiectasia and Rad3-related (ATR) and DNA-dependent protein kinase catalytic subunit (DNA-PK cs ).ATR phosphorylates H2AX upon replicative stress, 5 whereas ATM and DNA-PK cs are responsible for this phosphorylation upon DNA DSB, which are induced by IR. 6 ATM and DNA-PK cs have been studied on a qualitative basis focusing on their impact of repair pathway choice for rebuilding damaged DNA either via rapid (classical) non-homologous end joining cNHEJ and/or slow homologous recombination repair (HR) pathway. 7,8As for the pathway choice, the interplay between ATM and DNA-PK cs regarding IR-induced H2AX phosphorylation remains puzzling.Because although ATM is required, 9 DNA-PK cs can substitute for it. 10n this work we follow a model-based approach to analyze the contribution of DNA-PK cs and ATM to H2AX phosphorylation during the initial DNA damage sensing stage.Cucinotta et al. 11 have created a dynamic model solely focused on DNA-PK cs to predict dose and dose-rate effects on gH2AX dynamics.Very recently, a mechanistic model describing DNA damage complexity dependent sub-pathway choice in cNHEJ repair has been presented. 12Although several other mechanistic models of DNA-PK cs and cNHEJ repair exist, [13][14][15][16] mechanistic modeling of ATM dynamics in the context of DNA damage is rare. 17 computational model for ATM and DNA-PK cs interactions with regard to gH2AX activation integrating biochemical time course data is missing so far.We describe an iterative workflow to identify a predictive dynamic model involving ATM/DNA-PK cs mediated H2AX phosphorylation.Starting from several models, optimal experimental design (OED) was applied to optimize experiments for model identification.The identified model was used to analyze the dynamic contribution of ATM and DNA-PK cs to H2AX phosphorylation.

Model identification
2.1.1Defining network structures for cH2AX activation upon IR.The network structures (Fig. 1A) have been constructed based on meta-analysis 7,[17][18][19][20] focusing on the initial activation dynamic within the nucleus and the interplay between ATM and DNA-PK cs .DDR initiates with recognition of damaged DNA (DDNA1).Ku70/80 as a sensor for cNHEJ associates to the damage site (RC11) forming the DNA-PK complex (RC12). 21Then, the catalytic subunit of DNA-PK is either phosphorylated by active ATM or/and autophosphorylated at the T2609 cluster to initiate cNHEJ. 8The MRN complex (Mre11-Rad50-Nbs1), a sensor for the HR pathway, can also co-localize to the damage site to promote ATM autophosphorylation at Ser1981.
Failure of DNA repair via cNHEJ potentially allows HR proteins to access the damage site.This is modeled by splitting the initial DSB pool (DDNA) into DDNA1 and DDNA2, whereby DDNA2 is associated to HR and/or alternative non-homologous end joining (aNHEJ). 22Phosphorylation of H2AX can be achieved by active DNA-PK cs or active ATM.We generated four alternative models describing various interplays between ATM, DNA-PK cs and gH2AX (Fig. 1A).
2.1.2Experimental design for model calibration and identification.For model calibration purpose, an initial time course of H2AX phosphorylation in response to IR was studied in MDCK cells in a dose-dependent manner using 0.5, 1, 2, 5, 40 Gy.gH2AX levels increased with IR dose, while concurrently signal attenuation was delayed (Fig. 1B).These results agree with data from Burma et al. 9 From the competing network structures, we derived ordinary differential equation models and calibrated them (see Materials and methods).Simulations of the initial data set for all models are shown in Fig. 1C.Based on w 2 statistics, none of the models could be rejected at a significance level of a 0.05 = 0.05 (Table 1, OED 0).p-Values of Anderson-Darling (AD) residual statistics also indicated that all models seemed adequate for the initial data.
To discriminate between models, we subsequently designed (i) an IR double-pulse (Fig. 2A-D) and (ii) an IR double-pulse in combination with kinase inhibitors (Fig. 3).The IR doublepulse was parameterized with 2 design variables, namely interpulse time D 1 and second pulse dose D 2 , whereas the first pulse was fixed at 1 Gy.The objective was to maximize O = [T red hV ihSi] T .Herein T red is the reduced, modified T criterion to measure discriminative power, 23 whereas hV i, hSi represent mean model prediction variance and variance-entropy.The latter two criteria measure parameter information and distribution within the gH2AX signal (see Materials and methods).
For OED I, the optimal design D I * was chosen by trading off T red , hV i and hSi (Fig. 2B).Recalibration of all models to data from OED 0 and I, and additional inclusion of p53-P data (Fig. 2E) from titration experiments did not allow for model discrimination (all p-values 4 a 0.05 for both fit statistics; Table 1), but reduced prediction variances (Table 2).
Kinase inhibitors were employed for OED II to better dissect DNA-PK cs and ATM contributions.Titration of two highly specific inhibitors, namely Nu7441 and Ku55933 for DNA-PK cs and ATM, respectively, identified the optimal concentration for each.Further, we used the phosphorylation of p53 at S15 as a read-out to show the specificity of the inhibitors.Two successive pulses with different intensities (1 and 20 Gy) show in the immunoblot that the contribution of DNA-PK cs to this particular phosphorylation of p53 is marginal (Fig. 2E).This confirms earlier data. 24,25ED II was designed for three different inhibitor settings, namely Nu7441 and/or Ku55933.The estimated optimal design D II * potentially allowed for discrimination (Table 2, T red * c 1, Fig. 3A).The initial gH2AX peak showed a comparable reduction for both inhibitors.Phosphorylation of H2AX after the second pulse seemed to decay more rapidly for inhibited ATM compared to inhibited DNA-PK cs .Both inhibitors together showed synergistic effects on gH2AX (Fig. 3B).
According to the fit statistics of OED II (Table 1) only model A2 cannot be rejected in terms of w 2 .However, we find significant AD p-values for all four models, whereas models A2 and B2 have non-significant AD 3s p-values, which account only for residuals smaller than 3s.This behavior may be attributed to outliers in one of the experimental conditions (Fig. 1C and 3C) owing to experimental variations or deficits of the models in describing experimental conditions of OED 0, I, II.We selected model A2 as the final model for further analysis, since it was the only model with p-values of w 2 and AD 3s statistics exceeding a 0.05 for all 3 experimental runs.and ATM.To investigate the contribution of DNA-PK cs and ATM to H2AX phosphorylation, we analyze their times of maximal peak activity post irradiation.We simulated a single IR pulse from 1 mGy to 100 Gy (Fig. 4A-C).Active DNA-PK cs (DNA-PK cs -P) responds directly after irradiation within 2-10 minutes and shows fast signal attenuation.Response time of active ATM (ATM-P) in terms of maximal activity is delayed with respect to gH2AX and much more dose-dependent ranging from 10 to 56 minutes.These model predictions are in line with the literature: DNA-PK cs activation peaks at 10 minutes after IR treatment, whereas ATM has its peak activity at around 20 minutes. 26ccording to the model predictions, phosphorylation of H2AX is biphasic, following a dose-independent temporal activation order: the first activation phase of gH2AX right after stimulation is associated to DNA-PK cs , whereat the second phase is linked to ATM-P (Fig. 4A).The gH2AX signal decays on the scale of hours and correlates with ATM-P.This dynamics of fast initial and prolonged response is known from coherent feed forward loops, which serve as a signal persistence detector. 27At doses below 1 dGy peak level of gH2AX is dominated by DNA-PK cs , whereas ATM dominates above 1 dGy (Fig. 4B and C).For larger dose levels, ATM auto-phosphorylation results into a prolonged activation phase, with gH2AX peak activity shifted from 10 minutes at 10 Gy to 40 minutes at 100 Gy.
2.2.2 DNA-PK cs compensates inhibited ATM.Simulations of gH2AX dynamics with inhibited DNA-PKcs or/and ATM show that exclusive inhibition of ATM is nearly compensated by DNA-PK cs replacing the ATM associated activation phase of gH2AX by a prolonged DNA-PK cs associated phase (Fig. 5A left and B black vs. magenta).In contrast, DNA-PK cs inhibition results into loss of the DNA-PK cs associated activation phase.Owing to slower activation kinetics, ATM cannot compensate this delay (Fig. 5A left and B black vs. red).At doses where DNA-Pk cs dominates, gH2AX peak activity is delayed by roughly 45 minutes.Simulations of simultaneous inhibition of DNA-PK cs and ATM show a 3-to 10-fold reduction in gH2AX peak level, depending on IR dosage, whereas exclusive inhibition of either DNA-PK cs or ATM is not as much affecting peak activity of gH2AX (Fig. 5A right).For all inhibition scenarios, the biphasic phosphorylation kinetics of H2AX is lost.

Cell culture and treatment with c-irradiation
MDCK cells (ATCC CCL-34) were routinely cultured in RPMI-1640 supplemented with 10% fetal calf serum, glutamine and 100 U mL À1 penicillin and 100 mg mL À1 streptomycin, and incubated at 37 1C in a 5% CO 2 humidified incubator.The MDCK cells were seeded at a density of 2 Â 10 6 per 10 cm culture dish and cultured for 24 hours.The cells were irradiated with the Biobeam GM 2000 (Gamma-Service Medical GmbH, Germany) at a dose rate of 3.332 Gy min À1 using either single or double-pulse conditions.After a single pulse of 40, 5, 2, 1 or 0.5 Gy, the cells were harvested at 30, 90, 180, 300 and 720 minutes.The double-pulse consists of a single pulse of 1 Gy, followed 6 hours later by a second pulse of 20 Gy.The cells were harvested at 15, 35, 60, 160, 240, 370, 420 and 450 minutes.The inhibitors, Ku55933 (ATM, Tocris Bioscience, Germany) and Nu7441 (DNA-PK cs , Tocris Bioscience, Germany), used in the double-pulse setting, were added 30 minutes before first irradiation at a final concentration of 1 mM, either alone or together.The titration of the inhibitors were performed at 0,
The primary antibodies were incubated overnight in 5% skim milk in TBS-T at 4 1C.The membranes were washed thrice in TBS-T.The appropriate HRP-conjugated secondary antibody was added at a dilution of 1 : 5000 in 5% skim milk in TBS-T for 1 hour at room temperature, followed by three washes in TBS-T.The membranes were developed using a chemiluminescence substrate (Millipore, Germany).The respective bands were visualized using the ChemoCam Imager (Intas, Germany), followed by the estimation of the band intensities using ImageJ. 29ntibodies used in this work were as follows: LaminB2 (sc-133722) and HDAC1 (sc-7872) were obtained from Santa Cruz (USA, CA).gH2AX (ab26350) was from Abcam (UK).The secondary anti-rabbit-HRP or anti-mouse-HRP antibodies were from Jackson ImmunoResearch Laboratories Inc. (USA, PA).

Building of a dynamic signaling model network of cH2AX activation
Initially, 4 dynamic models in the form of ordinary differential equation systems were derived from the network structures in Fig. 1A and implemented in MATLAB using the solver CVODES. 30etails on the choice of kinetic rate laws are given in the ESI † Section S2.After the poor discrimination performance of OED I, we extended the models to contain also p53.The tumor suppressor p53 is an important effector protein during DDR.Phosphorylation of p53 at Ser15 by ATM promotes its release from MDM2 and results in p53 activation. 24,31Activation of p53 by DNA-PK cs has also been described. 32However, DNA-PK À/À MEFs show normal p53 activation. 25We did not find evidence for a DNA-PK cs contribution to the p53 phosphorylation (Fig. 2E), which agrees with earlier data. 33Therefore, we implemented the p53 activation as an ATM-dependent process only.As described in detail in the ESI, † 19 kinetic and 8 scaling parameters were estimated by maximizing the likelihood function, whereas the variance has been estimated from data replicates.Parameter estimation was performed for each model in an iterative manner, according to the 3 datasets, OED 0/I/II.Optimization of the likelihood function was performed iteratively, using a hybrid strategy.We combined a genetic algorithm ('ga' function from the global optimization toolbox of MATLAB), which was used to obtain a population of suitable starting solutions for a local, gradientbased optimization.Here we used an interior-point algorithm ('fmincon' function from the optimization toolbox of MATLAB).
Before analyzing DNA-PK cs , ATM and gH2AX dynamics with model A2, we performed an identifiability analysis based on the profile likelihood to assess the uniqueness of the model prediction and to also derive prediction uncertainty bands (see Fig. 4B and C).This analysis revealed that 8 kinetic parameters were not fully identifiable for the given optimization constraints, i.e. upper and lower bounds restricting the parameters to fall within 4 orders of magnitude.Six of these parameters were non-significant at the upper bound, whereas the other two were non-significant at the lower bound.One parameter was structurally non-identifiable.The non-identifiable parameters were not decisive for the question of kinase contribution to H2AX phosphorylation.More details on the identifiability analysis, parameter dependencies and impact on the prediction power are given in the ESI † in Section 2.

Experimental design criteria for model identification
Model identification is the process of comparing plausibility amongst models from a pool of competing computational models in the light of given experimental data.Plausibility is typically derived from some kind of lack-of-fit measure, for instance w 2 statistics.Experimental design for model identification aims at generating new experimental conditions and therefore data, to support this identification process in an optimal way using the models at hand.In the early phase of modeling a biochemical system with ODEs, parameters are typically very uncertain.Consequently, model predictions including design criteria are uncertain as well.Accounting for these uncertainties during design robustifies the optimal experiment against these uncertainties.In this work we use a multi criterion approach to identify optimal stimulus designs for model identification.We use three criteria that measure discriminative power, parameter information and its distribution along the time points of the model predictions for gH2AX.

Fig. 1
Fig. 1 Network structure and initial data (OED 0).(A) The network structures of four different models based on meta-analysis is shown as an interaction graph.Interactions are modeled via state transitions (arrows with squares), enzyme catalysis (lines with circles) and complex formation (joined lines).Stimulus and inhibitors have round-edge boxes.Abbreviations: IR ionizing irradiation; DDNA1 initial, damaged DNA; RC11 Ku70/80 to DDNA1 association; RC12 Ku70/80-DNA-PK cs complex; RDNA1/2 repaired DNA (cNHEJ/aNHEJ or HR); RC20 MRN complex to DDNA1 association; RC21|ATM MRN-ATM complex at damage site; RC22 RAD52 mediated repair complex; DDNA2, unsuccessful cNHEJ repair moved to aNHEJ/HR.Four mechanisms have been considered for branching (A1, A2, B1, B2).A and B refer to the location of the catalytic activity of ATM and indices 1 and 2 refer to the kinetic law used.For model index 1 branching to DDNA2 is catalyzed by the total amount of damaged DNA.Index 2 does not use the total amount of damaged DNA.(B) MDCK cells were irradiated with different doses and the insoluble nuclear extracts were analyzed by immunoblot.Lamin B2 or HDAC1 served as loading control.(C) Model simulation and quantified experimental data for OED 0 using the estimated band intensities of gH2AX.Data represent mean AE 2SD of 3-5 independent experiments.

Fig. 2
Fig. 2 Parameterization of the stimulus design, design criteria and respective immunoblots.(A) Parameterization of the stimulus design for OED I/II.(B) Design criteria predicted from the model simulations are plotted over the feasible design space.The optimal design point for OED I D I * and corresponding criteria O I * = [T red hV ihSi] T are indicated.(C) A representative immunoblot from an experiment based on D I * is shown.MDCK cells were irradiated as indicated and the insoluble nuclear extracts were analyzed by immunoblot.Lamin B2 served as loading control.(D) Corresponding model simulation describe the acquired data for gH2AX (model colors as in Fig. 1).Data represent mean AE 2SD of 3 independent experiments.(E) MDCK cells were irradiated as indicated.Inhibitors Ku55933 and Nu7441 were used at different concentrations and whole cell lysates were analyzed for p53-P and gH2AX.GAPDH served as loading control.Model simulation and quantified experimental data for p53-P are shown.Data of a single experiment.

Fig. 3 (Table 2 6 Â
Fig. 3 (A) D II * is obtained as in Fig. 2(B).(B) MDCK were incubated with 1 mM of the indicated inhibitor and irradiated as indicated.The insoluble nuclear extracts were analyzed by immunoblot.HDAC1 served as loading control.(C) The corresponding model simulations compare the acquired data for gH2AX before and after OED II (mean AE 2SD of 2-4 independent experiments, model colors as in Fig. 1).

Fig. 4
Fig. 4 Model predictions for the dynamic contribution of DNA-PK cs and ATM to gH2AX.(A) Simulated time courses of active DNA-PK cs and ATM and resulting biphasic gH2AX activity for IR pulses of different dose levels (1 mGy to 100 Gy).At larger dose, ATM shows a damped oscillation as a result of a positive feedback (autophosphorylation), which contributes to peak level of gH2AX at doses above 10 Gy. (B) Model prediction of the corresponding dose response in terms of time points at maximal activity of gH2AX, DNA-PK cs and ATM.Thin lines indicate 95% confidence regions of the model predictions, estimated from simulations along the profile likelihood.(C) Ratio of maximal DNA-PK cs -P to ATM-P.Thin lines indicate 95% confidence region of the model predictions, estimated as in (B).

Fig. 5
Fig. 5 Model predictions for the dynamic contribution of DNA-PK cs and ATM to gH2AX at inhibition for different dose levels (1 mGy to 100 Gy).(A) Peak time and peak level of gH2AX for indicated inhibitors (color code), (B) corresponding time courses.

Table 1
Fit statistics for initial (OED 0) and optimized experiments (OED I and II) Anderson-Darling p-values are indicated as AD.AD 3s indicates p-values of AD statistics where residuals larger than 3s have been excluded.The number of data points N data do not include the time point t = 0 [min].N y and N S indicate the number of estimated kinetic and scaling parameters 2.2 Model predictions 2.2.1 Biphasic control of H2AX phosphorylation by DNA-PK cs