Reverse engineering of drug induced DNA damage response signalling pathway reveals dual outcomes of ATM kinase inhibition

1 Inroduction

For the survival and normal functioning of a cell and ultimately of the whole organism, conservation and protection of the native DNA sequence and structure is necessary. However, our bodies are constantly exposed to diverse types of genotoxic insults from various sources that can damage cellular DNA and threaten the survival [1]. In order to maintain genomic integrity, eukaryotic cells have developed intricate DNA damage sensing and repair mechanism that combats the diverse sources of DNA damage and ensures survival [2]. The signalling responses that are triggered following DNA damage from endogenous or exogenous sources give rise to a pathway called DNA Damage Response (DDR) pathway [3]. The DDR pathway has evolved to be a complex, yet sensitive, highly integrated and interconnected pathway. Critical features of this highly robust repair mechanism are the ability to specifically recognize different types of DNA lesions and efficiently remove them [4]. Additionally, the different repair responses that are generated and the resulting network signalling triggered within the cell not only repairs the DNA lesion, but is also tightly linked with the cellular machinery that governs cell-fate decision e.g. cell cycle arrest to promote survival or apoptosis for programmed cell death [5]. Depending upon the scale and type of DNA damage, different repair responses can be activated with different outcomes for the cell (Fig. 1).

Central to DDR are the Ataxia telangiectasia mutated gene that codes for Ataxia Telangiectasia Mutated protein (ATM) and Ataxia Telangiectasia Mutated Rad3 related gene coding for ATR. ATM and ATR are Serine/Threonine kinases, functioning as core components of the DDR signalling pathway [6]. These kinases belong to phosphatidylinositol-3 kinase like kinase (PIKK) super family of large proteins having phosphatidylinositol-3/4 kinase (PI3K/PI4K) catalytic domain. ATM and ATR work as sensors of DNA damage and their kinase activities are rapidly induced following different types of DNA damage. Activation of these kinases requires autophosphorylation and this could be used as a biomarker for DDR activation [7, 8]. Following activation, these initiate signalling cascades by rapidly activating numerous downstream signal transducers and effector proteins comprising of kinases, phosphatases, transcriptional factors and chaperons. These effectors in turn modulate the cellular defence system, protein trafficking [9], mitochondrial regulation and metabolism [10], antioxidant pathway, cell cycle regulation, cell growth, recruitment of DNA repair enzymes and may trigger apoptosis if the DNA damage is beyond the repair capacity of the cell [11]. Pathways downstream of these kinases show extensive cross-talk, overlap and combination of specific and redundant activities making prediction of signal to response relationship difficult. Thus ATM function is implicated in number of responses resulting from DNA damage and its signalling decisions could be based on scale of damage, substrate availability and their degree of phospho induction [14]. This mandates a tightly regulated mechanism of interaction of ATM and ATR with their downstream effectors.

The targets of ATM- and ATR -mediated phosphory-lation in response to DNA damage are numerous. One of the crucially important targets of ATM/ATR is p53, a master transcription factor, functioning in the regulation of the cell cycle, induction of programmed cell death, regulation of cellular metabolism and control of aging [15]. p53 may be activated by the presence of damage in DNA or when the cell is deficient in nucleoside triphosphate precursors for DNA synthesis. Activating modifications in p53 results in decreasing its affinity towards the ubiquitin ligase MDM2 that in unstressed cells tags p53 for degradation via the ubiquitin-dependent pathways [16, 17]. ATM acts upstream of p53, phosphorylating it upon DNA damage, thus activating the p53-related pathways of cell cycle arrest and/or apoptosis. p53 may also be induced in an ATM-independent manner in cases when the damage is too severe and/or too extensive requiring the activation of various other signalling molecules [18].

Usually, moderate accumulation of p53 beyond certain threshold is associated with cell cycle arrests. This could be G1/S checkpoint arrest involving CHK2 and the retinoblastoma protein (pRB1) [22] or G2/M checkpoint arrest via ATM-CHK2 and ATR-CHK1-mediated signalling [23]. Suprathreshold levels of p53 are associated with induction of apoptosis via pRb-E2Fl, CHK2, p73 and the pro-apoptotic proteins of the BCL-2 family, e.g. BAX [26]. Among the other immediate phosphorylation targets of the ATM/ATR system are proteins that are directly involved in DNA repair, such as BRCA1, RAD17, SMC1, and CHK1 [31, 32]. Transient S-phase cell cycle arrest may be induced through NBS1, BRCA1, FancD2 and SMC1 [33], usually following ionising radiation.

Induction of DDR is a complex process involving activation of pathways that may be extensively mod- ulated by various factors e.g. type and scale of damage, the phase of the cell cycle, transcription status of the genomic region/s that have sustained damage, the differentiation status of the cell, the overall level of genotoxic (e.g. oxidative) stress, the presence of signs of neoplastic transformation etc. Nuclear proteins of the non-histone type, such as the high-mobility group (HMG) proteins may also actively modulate the recog- nition of lesions in DNA and the recruitment of the repair machinery at the damaged sites [20, 37]. The levels of HMG proteins may play a crucial role in maintaining the physiological minimum of unrepaired damage in DNA as many tumours are characterised by up-regulation of HMGA. However, inherently low levels of HMGA (e.g. resulting from polymorphisms in the respective gene/s) are associated with accumulation of oxidative lesions that ultimately may lead to DNA damage, aging and/or cancer [40, 41]. In addition, in terminally differentiated cells (e.g. adult neurons, cardiomyocytes and others) and sometimes, in cells that are only partially differentiated (e.g. monocytes) DNA repair may be selectively suppressed [42].

There may also be individual variation in the efficiency of response to genotoxic stress and this may significantly vary even between clinically healthy people. Several methodologies exist for assessment of capacity for DNA repair, providing a measure for overall repair capacity [46] or a composite measure of the DNA repair status according to individual biomarkers [50].

Over the last two decades the systems biology community has witnessed rapid advances in system identi- fication methods development including the following: linear, loglinear and nonlinear differential equations; artificial neural networks [53]; genetic algorithm [54]; evolutionary optimization with data collocation [55]; interval analysis [56]; alternating regression [57]; parameter estimation for noisy metabolic profiles using newton-flow analysis [58]; simulated annealing [59]; ant colony optimization algorithm for parameter estimation and network inference [60]; substitution of slopes for differentials; dynamic-flux estimation [61]; eigen-vector optimization [62]; transposive and repressive regression method [63]; and so on. For instance, the following presented ODE-based methods that use dynamic data and steady-state measurements to capture and identify complex systems: Sorribas and Cascante et al. [66]; Irvine [67]; Savageau et al. [68]; Tominaga and Okamoto [69]; etc. Other identification methods may involve techniques based on Boolean network (BN) [70], Linear programming [71, 70], and least-squares minimisation criterion [68, 72].

In deterministic modelling methods that conceptualise outcomes of causality in a system are sought to predict historical and new outcomes based on scientific and theoretical evidences. In formulating an automated deterministic model, the model structure and parameter search space domain may be adapted to recalibration based on known principles, theories, or known formalism. Temporal dynamics of complex processes involved in biological pathways such as the DDR pathway may be captured and profiled by measuring system properties over time. An assumption-free reverse engineering of only small quantities of these time series measurements taken at regular intervals may deterministically construct dynamic models of such complex systems [63, 73]. This requires a novel network inference algorithms that is capable of demonstrating complete recovery of a wide range of networks of interactions from experimental time series data [63] as illustrated in Figure 2.

In identifying some of the challenges associated with modelling approaches over the last decade, we learn that traditional modelling approaches may help explain some of the functions of the key components involved in biological systems. However, it may be necessary to avoid making any a priori assumptions about any underlying mechanisms involved within the systems.

Our proposed reverse engineering framework is developed specifically to address all the need to target and optimally utilise limited (extremely small e.g. 3 time points) time series data. We found that almost all existing inference algorithms require large quantities of time series data to work and none often guarantees data consistency and accurate reproducibility of experimental time series data under such data limitation. Voit et al. [74] cautioned on common mistakes and potential hindrances to effective identification of biological system identification using in vivo time series data. These include algorithmic difficulties of nonlinear regression analysis, validity and consequences of incorrect a priori assumptions made in model design.

Here, we seek to create data consistent, self-reconfigurable dynamic models [63, 73] of the DDR pathway using only semi quantitative experimental data and then analyse such constructed predictive models to generate in silico topological maps in order to generate new insights into the topological protein networks underlying the measured signalling profiles following two different magnitudes of DNA damage. Systems analysis and study of these interaction maps are then used to identify key signalling alterations in the DDR network that may explain or interpret the measured and observed DNA damage dependent consequences of ATM inhibition and signal to response inter-relationship. Novel network inference (reverse engineering) algorithm is used to identify unknown signalling alterations that result from both ATM inhibition and without. The reverse engineering strategy introduced is applicable to any complex system that requires signalling targets to be identified. We demonstrate an effective system identification procedure for uncovering potential signalling targets for sensitising genotoxic outcomes in cancer cells.

2 Materials and Methods

The need to develop a mathematical and theoretical framework for creating multiple dynamic (deterministic, reconfigurable, self-organising, unsupervised, and automatically constructible) models for optimally utilising experimental time series data has resulted in a novel reverse engineering strategy that often guarantees the production of highly consistent system identification solutions. In practice, this may involve an automated inference or extraction of unsupervised models of time series data that are useful for making accurate simulation of historical data and predictions about other unknown system states and network of interactions.

In order to characterise time dependent DDR signalling profiles induced following genotoxic challenge and illustrate how such activation takes cell fate decision, we have used immortalised human keratinocyte cell line (HaCat) model and generated quantitative time series data of kinetics of DDR substrate phospho inductions following different treatment regimes. Cells were treated with a lower (0.1µM) and higher (0.4µM) concentration of the widely used radiomimetic drug, Doxorubicin (Dox) for treatments, with and without 10µM of ATM kinase inhibitor, KU55933 (Ku). The two concentrations of genotoxic agent were chosen in order to delineate the corresponding signalling dynamics at a lower and higher degree of DNA damage and characterise the signalling alterations upon a repairable DNA damage and irreparable state (apoptosis). The time series experiments under different conditions (Table 1) were performed targeting all major proteins that respond to DNA damage (Table 2) and carried out their semi quantitative analysis. Relative measurement of protein inductions were performed by quantitative Enzyme Linked Immunosorbent Assay (ELISA) while parallel cytotoxicity assays resulting from these treatments were performed using Neutral Red (NR) uptake assay in order to determine the consequences of the activated signalling profiles on cell health and survival.

The acquired time series data are analysed and modelled using one of the network inference algoithms introduced in [63, 65]. These reverse engineering techniques are reintroduced in subsections 2.2.3 and 2.2.4.

These algorithms only require time series measurements with at least three time points to work and often converge to a unique network of interactions (solution) as the number of time points increases. Since the algorithms only work on time series data with regular time intervals, we have applied them on subsets of the above-mentioned experimental data, i.e. time points 0hr, 2hr and 4hr; time points 0hr, 4hr and 8hr; time points 0hr, 4hr, 8hr and 12hr; time points 0hr, 4hr, 8hr, 12hr and 16hr; time points 0hr, 4hr, 8hr, 12hr, 16hr and 20hr; and time points 0hr, 4hr, 8hr, 12hr, 16hr, 20hr and 24hr, gradually and progressively. The temporal models constructed and topological networks of interactions generated for each of these 6 subsets have been viewed and analysed to inform and guide intermediate biological protocols and experiments for validation and verification purposes. Such an optimal utilisation of semi quantitative experimental data combined with simultaneous generations of topological maps is being used to direct and advance our investigative study into the role of ATM in DDR pathway.

In this section, we describe the biological experiments that were performed to produce the time series data, fundamental deterministic modelling approach that has been developed to process such data, an- other alternative recast technique and modelling strategy for identifying such complex systems, and systems analysis/modelling challenges encountered during the process.

2.1 Experimental methods

2.1.1 Cell lines, culture conditions and treatments:

Immortalised human keratinocytes (HaCaT) were maintained in DMEM without phenol red, supplemented with 10% foetal bovine serum (FBS), 2 mM glutamine, 1 mM sodium pyruvate, 100 µg/ml streptomycin and 100U/ml penicillin in an atmosphere of 5% CO₂. 1mM stock of Doxorubicin (Sigma-Aldrich) was prepared in sterile H₂O, filtered and stored in the dark at 4°C. Final concentration of treatments was established on the day of experiment by diluting the stock in media. 10mM stock of ATM kinase inhibitor, KU55933 (Calbiochem) was prepared in DMSO (Fisher Bioreagents), aliquoted and stored at -80°C in dark. Final concentration of 10µM was achieved by diluting the stock in media on the day of experiment.

2.1.2 Enzyme Linked Immunosorbent Assay (ELISA)

ELISA was performed for cells grown in opaque flat bottom 96-well tissue culture plates using Amplex Red detection system (Invitrogen, UK). Either media only or media containing cells at a density of 2.5 x 103 were seeded per well in quadruplicates and allowed to grow for 24 hours (hr). Following this, cells were either left untreated or treated with drugs for different time points as indicated. At the end of treatments, cells were washed thrice with PBS, fixed with 150µL of 3.5% Paraformaldehyde (3.5gm Paraformaldehyde, 100µL 0.1M NaOH in 100mL PBS, heated to 70°C to dissolve, cooled and pH adjusted to 7.2 using concentrated HCl) for 30 minutes (min). Cells were washed with ice cold PBS and permeabilized using 0.5% triton X-100 made in PBS, washed with ice cold PBS and quenched using quenching buffer (50mM NH₄Cl solution made with H₂0) and washed thrice with ice cold PBS. Cells were next blocked using blocking buffer (2.5gm Bovine serum albumin, 5mL Foetal bovine serum (FBS), 5µL Triton X-100 in 50ml PBS) for 1 h at room temperature, washed with PBS and either incubated with blocking buffer alone (no primary antibody controls) or blocking buffer containing primary antibodies as indicated, for 2 hours (hr) at room temperature. Following incubation, cells were washed twice with washing buffer (0.01% Tween20 made with PBS) and incubated for 1 h with Horseradish peroxidise conjugated secondary antibodies diluted 1 in 1000 in washing buffer. Meanwhile, Amplex Red detection reagent was prepared by diluting 50µL of 10mM Amplex Red (previously made in DMSO), 500µl of 20mM H₂O₂ (Sigma-Aldrich) in 4.45 mL of the provided 1x reaction buffer to achieve final concentration of 100µM Amplex Red reagent/2mM H₂O₂. 100µl of the reagent was added per well and incubated in incubator for 30 min. Readings were taken in a 96-well fluorescence multi plate reader (Modulus template, Turner Biosystem) and fluorescence intensities in each well were measured using excitation/emission of 530/590 nm. Values are means of the fluorescence acquired from quadruplets of treatments and normalized to the mean value of untreated control (UT) for that protein and expressed as 1. The data is based on three independent experiments.

2.1.3 Cell cytotoxicity assay

To assess the cytotoxic effects of drug treatments, Neutral- Red (NR) uptake assay was performed. For this, either media alone or media containing cells at a density of 2.5 x 10³ were seeded per well in quadruplicates and allowed to grow for 24 hr. Following this, cells were either left untreated or treated with drugs for different time points as indicated. NR reagent was prepared by adding 0.33g of neutral red powder (Sigma-Aldrich) in 100ml of distilled sterilised H₂O, filter sterilised and stored in dark at room temperature. On the day of assay towards the end of drug treatments, NR reagent was diluted 1 in 100 in media (to obtain a final concentration of 33µg/mL NR), cells washed thrice with warm PBS and 125µl of the media containing neutral red reagent was added per well. Plates were incubated at 37°C for 5 h. Following this, cells were gently washed thrice with PBS. To release the up taken NR, cells were incubated with 100µl of NR desorb (1% glacial acetic acid and 50% ethanol made with distilled H₂O), shaken on a rotary shaker at 100 rpm protected from light, followed by letting the plate stand for 5 min. Readings were taken in 96-well multi plate reader (Modulus template, Turner Biosystem) and absorbance in each well was measured at 540 nm.Values are means of the absorbance acquired from quadruplet of treatments and normalized to the mean value of untreated control (UT) and expressed as % cell death with the formula 100 - [(Abs₅₄₀ treated sample/Abs₅₄₀ untreated sample) x 100]. The data is based on three independent experiments.

2.1.4 Immunocytochemistry

For immunocytochemistry, exponentially growing HaCat cells were seeded at a density of 5 x10⁴ cells in complete media onto poly-L lysine (Sigma-Aldrich) coated cover slips placed in a 12-well tissue culture plates. Next day, following relevant treatments, cells were washed three times with ice cold PBS, fixed in 3.5% paraformaldehyde prepared as mentioned before and cells incubated at room temperature for 30 min. Following this, cells were gently washed twice with 1 ml of PBS, permeabilized with 0.3% triton X-100 for 10 min, and following three washed with PBS, blocked with blocking solution (1% goat serum, 1% bovine serum albumin in PBS containing 0.05% Triton X-100) for 30 min. Cells were then incubated with relevant primary antibodies (Table 3) diluted in blocking solution for 1h, washed three times with 0.1% Triton X-100 in PBS for 5 min, and then incubated with either Alexa Fluor 488 conjugated anti-rabbit or anti mouse antibodies (green fluorescence) or Alexa fluor 568 (red fluorescence) conjugated antibodies (Invitrogen) for 30 min. Only secondary antibody incubations (no primary antibody control) was also performed (data not shown). After subsequent washing three times with the 0.1% Triton X-100 in PBS for 5 min, cover slips with cells were mounted on slide using 4’,6-Diamidino-2-Phenylindole, Dihydrochloride (DAPI) containing mounting reagent (Vectashield, Vector Laboratories) to provide for nuclear reference. Cells were imaged under relevant filters with a Leica DMiRe2 electronic microscope.

2.1.5 Imaging and analysis

Fluorescence images of immunocytochemistry were collected under relevant excitation and emission filters de- pending on the fluorotype under Leica DMiRe2 electronic microscope equipped with iXonEM +897 EMCCD camera (ANDOR Technologies Ltd). Images were visualised using multi dimensional microscopy software Andor Module iQ Core. Colocalization assays were performed and determined with software integral features supplied by Andor IQ core software features. Data were generally expressed as mean for individual sets of experiments.

2.1.6 Coomassie staining of cells

Coomassie staining of cells following different treatments were performed in order to normalise the protein induction data of ELISA on cell density. At the end of treatments, cells in 96 well plates were washed with PBS and fixed using 30% methanol and 10% acetic acid made in distilled H₂O for 5 min. The fixing solution was removed and 100µl of coomassie stain (0.05% coomassie brilliant blue made with fixing solution and filtered) was added and plates incubated for 20 min. Following this, 50µl of 0.1M NaOH was added per well to release the dye and cells incubated for 30 min. Readings were taken in 96-well multiplate reader (Modulus template, Turner Biosystem) at 620-650nm.

2.2 Reverse engineering methods

2.2.1 Deterministic mathematical modelling

Reverse engineering may be viewed as a network inference challenge that seeks to obtain from data accurate approximation of underlying networks of interactions or optimal estimation of model parameters. Hence effective methods for constructing models that are data consistent and self-reconfigurable are required. The demand for data consistency is necessary in automated model construction in satisfying minimum requirements. Such constructed models are expected to be consistent in their capability to simulate historical time series data and predict new outcomes. However, system identification challenges can be extremely difficult. We have developed novel and efficient inference algorithms that are capable of dealing with contemporary modelling challenges including systems identification and parameter estimation [63, 65]. Our ultimate goal is to formulate a theoretical framework for automating the construction of dynamic models and in silico topological maps that describe and represent acquired time series data sets. Such assumption-free predictive models are powerful tools that may be used to determine unknown evolutionary dynamics and produce sufficient explanation to experimental data.

2.2.2 The recommended reverse engineering method

Our data-driven modelling approach is based on the principle that a set of fundamental inference methods that satisfactorily demonstrate successful recovery of networks of interactions from artificially simulated ex- periments and assessment tests may be useful for effective reverse engineering and mathematical modelling of real experimental data [63, 65] (Fig. 2). Further development and assessment of those methods have identified and suggested an optimal method that is being used in this reverse engineering study. We have successfully identified a unique modelling strategy that can be used for analysing and modelling time series data without making a priori assumptions or preconceived knowledge about the underlying network of interactions of the target system. This method is particularly useful for inferring a network of interactions that is able to explain the experimental data supplied [63]. In further study [65], an alternative power-law based recast transformation method was integrated into the modelling framework to establish a secondary implementation of a multimodel integrated solution (Fig. 3) [63, 65]. The inference methods we have developed [63, 65], illustrated in Figure 3 and summarised in subsections 2.2.3 and 2.2.4, are used to construct jacobian and power-law based models from data.

Figure 3 illustrates the proposed dynamic modelling strategy used in this study which highlights three basic types of output results: heatmap; plotted graphs; and an in-silico topological maps of network of interactions. This schematic diagram (Fig. 3) depicts the two mathematical methods of analytically processing time series data of DDR pathway to automatically construct either a jacobians or power-law models of the real biological systems approximated by real time series measurements. These two methods are based on the algorithms developed and presented in [63, 65].

2.2.3 Theoretical understanding of inverse problems using jacobian-based implementation

In a system of ordinary differential equations an inverse problem may be defined in the form

with the solution given as

where X represents system states (e.g. from one state, i.e. X_before, to another state, i.e. X_after. Both X_after and X_before) are known vectors of same length, Ẋ can be calculated and J, the unknown jacobian (transition matrix) must be inferred from data. The unknown transition matrix J may be formulated in the matrix form

(1)

where is a known state vector, i.e. the t^th vector of the given time series X. According to Idowu and Bown [63], the inverse problem may be reexpressed as the multi-state reformulation

redefined as

assuming that the time series (state vector) measurements are measured at regular time intervals of t_c. This new equation is completely solvable using the transpositive regression method (TRM) [63]. The exponential matrix E:

must first be derived before calculating the parameters of J.

2.2.4 Theoretical understanding of inverse problems using power-law method

Altenatively, a power-law based (e.g. Half-system , S-system, or the Generalised Mass Action (GMA)) model of a dynamical system may be used in solving an inverse problem. For instance, the Half-system representation is in the form:

where i = 1 ... n; m is the number of dependent variables; g_ij (kinetic orders) quantify the overall net efect of each X_j on the production (and/or degradation) of X_i ; α_i are called rate constants. In matrix form, this model representation may be defined as

(2)

which is equivalent to the logarithmic equations

(3)

such that

(4)

where the unknown parameters are the matrix collection M.

Hence conceptualising equation 4 (variant or redenition of the X(_after) = M * X(_before) form) as though the unknown matrix M is to be calculated, this inverse problem is completely solvable using the TRM method.

2.2.5 Jacobian-to-Half-system model integration

In [65] Idowu and Bown presented new methods that establish the relation between parameters of jacobian model and those of Half-system, otherwise known as kinetic orders based on the fundamental relation

(5)

where X_j and X_k are selectable entries in the given time series data and the partial derivatives ∂Ẋ_i and ∂X_k are interrelated entries of the inferred jacobian matrix which are consistent with the actual data. Ultimately the following matrix-based equation can be derived [65]:

(6)

from which the parameters (kinetic orders) of the Half-system model,

is factorable in matrix product factors

(7)

As a consequence of the above factorisation data-consistent jacobian relation may be expressed also as the relation

(8)

3 Results

3.1 Biogical systems: results

3.1.1 Time dependent treatment of 100nM Doxorubicin demonstrated cytoprotective function of ATM signalling:

In order to characterise the mechanism and dynamics of DDR signalling following double stranded DNA damage, time series treatments were performed with commonly used chemotherapeutic drug, Dox and the resulting DDR pathway activation was studied. For this, HaCat cells seeded in 96 well plates for 24 h were either left untreated, treated with a low concentration of 100nM of the radiomimetic drug, Dox alone, or with the a co-treatment of 10µM ATM kinase inhibitor, KU. In addition, cells grown in a separate 96 well plates were treated in the same manner and subjected to cytotoxicity analysis via NR-uptake assay. This was done to allow for a correlative study of the changing dynamics of DDR signalling following treatments with the associated consequences for cell viability. Doxorubicin treatment in this section of experiments was set to 100nM because this is a physiologically relevant concentration currently being employed in chemotherapy [26]. Immunofluorescent labelling of total ATM revealed that in HaCat cells, ATM is mostly localised in the nucleus along with some cytoplasmic localisation (Fig. 4). The localisation of ATM in the cytoplasm was previously identified to be in Golgi apparatus [9]. Treatment with 100nM Dox resulted in no major change in levels of total ATM. However, the appearance of nuclear foci following time dependent Dox treatment was evident. Appearance of these damage induced ATM foci illustrated that DDR pathway is active in HaCat cell line (Fig. 4a). KU treatment caused inhibition of the foci demonstrating ATM dependence for the foci formation (Fig. 4b). Immunofluorescent labelling of phosphorylated ATM on serine 1981 (pATM) on the other hand showed time dependent nuclear induction following Dox treatment (Fig. 5a) while KU treatment caused disruption of the observed Dox induced nuclear increase (Fig. 5b). This illustrated the requirement of kinase activity of ATM for its own phosphorylation, a feature of ATM that is well studied elsewhere [7].

ELISA for pATM following the above treatments showed consistent results. pATM induction was evident at the first time point of 2 hr treatment with a sustained induction for up to the final time point of 24 h of treatment (Fig. 6a). This showed that a lower scale of DNA damage with 100nM Dox was still enough to induce pATM induction. However, the kinetics of induction were slow with < two fold maximum induction at 24 hr time point. Co-treatment of ATM kinase inhibitor, KU, with Dox resulted in disruption of pATM induction for up to 12 hr of treatment. Furthering the treatment beyond 12 h resumed pATM induction in a time dependent manner (Fig. 6b). Hence, it appeared that while KU caused earlier disruption of pATM levels, at later time points, the induction slightly resumed resulting in overall delayed kinetics of pATM activation. A possible explanation of this is the fact that KU is a reversible inhibitor of ATM and may only form a transient effective inhibitory complex and lose its affinity during long term treatments in the presence of DNA damage. Furthermore, KU was previously shown to rapidly upregulate total ATM levels resulting in a transient induction of pATM which may account for that later induction.

Cell cytotoxicity analysis following same treatments revealed that inhibition of ATM during 100nM Dox treatment sensitised cells to the genotoxic agent at almost all the time points tested (Fig. 7). For example, at 24 hr time point, there was only 30% cell death in Dox treated cells as compared to 45% with the inclusion of KU. This confirmed the role of ATM signalling in cytoprotection and survival at this particular extent of DNA damage. This is consistent with the typical protective function of ATM whereby it is known to exert cell cycle arrest and promote DNA repair in order to ensure cell survival at a repairable scale of DNA damage. This is through its phosphorylation of P53 which exerts CIP/KIP mediated cell cycle arrest by sequestering E2F1 activities [23], as well as its activation of checkpoint kinases Chk1 and Chk2, which phosphorylate and inhibit phosphatase activities of Cdc25 [75] stopping the cell cycle progression. While ATM is also implicated in apoptosis induction [26], it could be inferred from fig. 7 that at 100nM Dox, its activation resulted in cytoprotection.

Unsurprisingly, pATR S-428 did not show any major change through all the time points of Dox treatment (Fig. 6a). This can be explained by the fact that ATR activity is mostly implicated in DNA damage following UV exposure or replication folk stalling [76] and that a lower magnitude of double stranded DNA damage alone is not sufficient to induce pATR. Interestingly however, with the addition of ATM inhibitor, pATR showed an induction post 12 h of Dox and KU co-treatment and continuing further until the 24 hr time point (Fig. 6a). This later activation of ATR following double stranded DNA damage is consistent with a previous report which indicated that double stranded DNA damage may still activate ATR but with delayed kinetics as compared to ATM activation [77]. pATR induction due to sustained DNA damage in the absence of ATM activity could either be as a result of emergency cellular response where ATR activation may represent a compensation for the loss of ATM cytoprotective signalling, or could be to rather initialise apoptosis. Among the two possibilities, the latter seems likely because cell cytotoxicity analysis revealed correlation between pATR induction and cell death during KU + Dox treatments pointing out to a pathway involving KU ┤ ATM = pATR → Apoptosis. As such, the role of ATR in apoptosis induction has already been reported following treatments with chemotherapeutic agents [78, 79].

Phosphorylated forms of checkpoint kinases, Chk1 and Chk2, the main effectors for ATM and ATR signalling [80] showed similar behaviour following the milder genotoxic challenge with inductions following 12 hr of treatment (Fig. 6a). pChk2 demonstrated an overall greater induction as compared to pChk1, fitting well with pChk2 being the primary effector molecule of pATM following double stranded DNA damage [81]. The fact that both pChk1 and pChk2 showed induction supported the earlier hypothesis of ATM mediating checkpoint kinase dependent cell cycle arrest. The late kinetics of activation of checkpoint kinases possibly pointed towards G2/M arrest of the cell cycle rather than induction of apoptosis through ATM → pATM → pChk1/pChk2 ┤ Cdc25 = G2/M cell cycle pathway. This assumption is based on the fact that KU treatment inhibited ATM dependent checkpoint kinase induction, which was accompanied by greater degree of cell death whereas the activation of these checkpoints in ATM functional state correlated with reduced cell death (Fig. 7). Contrastingly, examination of phosphorylated form of P53 on Ser-15, another substrate of ATM, revealed time dependent induction that continued steadily for 20 h of 100nM Dox treatment (Fig. 6b). ATM dependent pP53 induction following mild DNA damage has been reported to result in G1/S cell cycle arrest [82]. Hence, this pP53 induction can contribute in causing cell cycle arrest and induce DNA repair signalling possibly via pP53 → CIP/KIP ┤ cyclin/Cdk = Cell cycle arrest , in addition to ATM → Chk1/2 ┤ Cdc25 pathway.

pP53 induction was disrupted with the addition of KU during Dox treatments for 12 hr. This reconfirmed ATM dependence for the observed Dox induced pP53 increase and supported the above mentioned pathway. Treatments longer than 12 hr resulted in ATM independent increase in pP53 (Fig. 6b). Time dependent and ATM independent increase in pP53 tightly correlated with a surge in cell death in the cytotoxicity assay probably demonstrating P53 mediated and ATM independent apoptosis post 12 hr (Fig. 7). The likely underlying signalling pathway could involve DNA DAMAGE + KU = ATR → pP53 → Apoptosis. Another likely kinase responsible for P53 induction could be DNA-PK [83]. Next, levels of E2F1 transcription factor, implicated in both cell cycle progression [84] as well as apoptosis [85] were examined. It was found that lower Dox treatment rapidly resulted in reduction of E2F1 levels reaching its lowest at 4 hr of treatment and demonstrating a single oscillatory behaviour (Fig. 6b). Hence, at lower scale of damage, pATM activation that resulted in pP53 induction was also accompanied by reduction in E2F1 levels. It is already known that activation of P53 results in Rb induced block of E2F1 activity causing G1/S arrest [86]. However, the actual ATM dependent reduction of E2F1 protein levels could also be explained by the fact DNA damage induced ATM was reported to phosphorylate and destablise MDM2 protein [87] which not only results in the observed P53 stabilisation, but may also leads the loss of E2F1 stability [88]. Taken together, the earlier assumption of pATM induction at earlier time points of 100nM Dox causing a G1/S arrest is supported these observations, while the involvement of checkpoint kinases also suggested G2/M checkpoint arrest.

Further support of pATM induced G1/S arrest via DNA damage → ATM → P53 → Rb ┤ E2F1 = Cell cycle arrest and ATM ┤ MDM2 → E2F1 pathway is provided by the fact that addition of KU with Dox treatment resulted in not just loss of this repressive effect on E2F1, but also led to a time dependent induction with a correlation with cell death. Since this induction was ATM independent, other mechanisms may have been involved to account for it. Indeed E2F1 activation following Dox treatment has also been reported in pATM and pChk2 independent mechanisms [89, 90].

The γ-H2AX levels, a marker for double stranded breaks, as expected, were induced starting at first time point examined and showed a continuous time dependent increase following 100nM Dox treatment. Treat- ment with ATM inhibitor resulted in total disruption of this induction at all the time points tested. This indicated total ATM dependence for γ-H2AX formation after doxorubicin treatment. Similarly, BRCA1, an important downstream substrate of ATM also showed ATM dependence for its induction (Fig. 6b). While BRCA1was reported to induce apoptosis in both P53 dependent and independent mechanisms [91], under the conditions tested, and it terms of ATM signalling, BRCA1 exhibited a complete dependence on ATM for its induction (Fig. 6b). This is consistent with a previous report showing ATM dependence for phosphorylation of BRCA1 at serine 1524 [92]. Since BRCA1 is a major component of ATM induced Homologous recombination repair, this further supports the earlier assumption of ATMs role in cytoprotection in mild DNA damage. Taken together, the data generated with a lower DNA damage of 100nM Doxorubicin demon- strated that ATM kinase functioned as a DNA damage sensor causing G1/S or G2/M arrest mediated by checkpoint kinases and P53 and DNA repair by BRCA1 induction rather than apoptosis, indicating its role in cytoprotection. Pathway analysis also revealed differential kinetics of phospho induction of ATM/ATR substrates as well as distinct sensitivities to ATM inhibitor.

3.1.2 Treatment with higher concentration of 400nM Doxorubicin demonstrated time dependent dual roles of ATM signalling:

Once the kinetic parameters of DDR activation at relatively low level of DNA damage were determined, we next set out to perform treatments and examine the activation kinetics at a higher scale of DNA damage. To achieve this, the concentration of Dox was increased by 4x to 400nM. This was done to determine whether increasing DNA damage would still induce ATM signalling towards cytoprotection as observed with a lower concentration of Dox or alternative switch it towards apoptosis and hence present a different underlying DDR profile. We also wanted to examine whether increased DNA damage would cause the resulting DDR activation more ATM independent.

Firstly, cytotoxicity assay of cells treated with either 400nM Dox alone or with the addition of KU revealed greater cell death in both cases as compared to that at lower Dox concentration (Fig. 8). Interestingly however, time dependent co-treatment of 400nM Dox and KU resulted in enhanced sensitisation and cell death only at earlier time points of 2, 4 and 8 hr of treatments as compared to 400nM Dox alone. On the contrary, at 12, 16, 20 and 24 hr of co-treatments, there was desensitisation with diminished cell death upon the addition of KU as compared to treatments with Dox only (Fig. 8).

From this, it could be inferred that at earlier time points of Dox treatment where DNA damage repair could still be manageable, ATM kinase might mainly perform a cytoprotective function by signalling in cell cycle arrest and DNA repair. Whereas, at a greater and unmanageable extent of DNA damage i.e. beyond 8 hr of 400nM Dox treatment, its function might have switched towards apoptosis, as there was higher cell death in functional ATM state and comparatively lower when it was inhibited. This context dependent functional switch in ATM activity must also implicate an altered and distinct underlying signalling dynamics bringing about the observed re-adjustment in cell fate.

Treatment of HaCat cells with 400nM Dox resulted in a rapid pATM induction reaching >3 fold at the first time point of 2 hr and maintaining time dependent induction up to 16 hr of treatment (>5 fold increase, Fig. 9a). Co-treatment with KU caused disruption of the Dox induced upregulation of pATM as expected with profound effects on cell survival as well. For example, 2 hr treatment of 400nM Dox resulted in 15% cell death while addition of KU increased cell death to 25%. During 12-24 hr of Dox and KU treatment, time points where cells showed lower cell death as compared to Dox only, pATM induction remained inhibited as opposed to very high inductions in Dox only (Fig. 9a), supporting the earlier assumptions i.e desensitisation is caused due to inhibition of apoptotic function of ATM during extensive DNA damage.

Interestingly, at a higher damage of 400nM Dox, pATR also showed induction, starting at 2 hr and maximum at 16 hr of treatment. The kinetics of this induction, however, were slower than those of pATM consistent with a previous report showing slower kinetics of ATR activation following DSBs [77]. Furthermore, co-treatment of KU + 400nM Dox resulted in a rapid and a high induction of pATR starting at 2 hr post treatment and maximum following 12 hr of treatment (Fig. 9a). While Jazayeri A et al., [93] reported a functional requirement of ATM for ATR function in double stranded DNA damage, this result suggested ATM independent induction of pATR greater than that observed with milder Dox + KU treatment. Analy- sis of this induction with the results obtained from cytotoxicity assay indicated correlation of this transient pATR induction during Dox + KU treatment with cell death. For example, with the addition of ATM inhibitor during 400nM Dox, pATR was up-regulated at 2, 4 and 8 hr time points where cells showed higher cell death, whereas, it was rapidly down regulated at 12, 16 and 20 hr, time points where cells showed an overall lower cell death (compare figs. 8 and 9). Hence one candidate for inducing higher apoptosis in ATM inhibited states until 8 hr of treatment could be ATR activity.

In terms of checkpoint kinase activities of Chk1 and Chk2, while lower Dox treatment revealed comparable Dox induced activation and KU induced disruption of these kinases, at a higher scale of DNA damage with 400nM Dox, these kinases showed both significantly different inductions following Dox, as well as sensitivities toward ATM inhibition (Figs. 8 and 9). Overall, pChk2 showed earlier induction following Dox treatment as compared to pChk1, starting at 2 hr time point and continuing up to 16 hr of treatment and staying constant following that. pChk1, on the other hand, showed a delayed induction but a very high single spike at 16 hr of treatment and reducing at further time points of 20 and 24 hr.

Addition of ATM inhibitor along with 400nM Dox resulted in total disruption of pChk2 levels at all the time points, confirming that the dependency of Chk2 on ATM for its Threonine 68 phosphorylation was unaltered even at higher scale of DNA damage during which, greater cell death was seen than any other treatment in this study. These results further demonstrate that pChk2 induction, which was totally dependent on ATM activity, was implicated in effecting the cytoprotective function of ATM at relatively lower scale of DNA damage, whereas promoted apoptotic signalling at a higher scale of DNA damage.

This result is also supported by previously well-explained ATM dependent role of Chk2 in apoptotic induction [94, 95]. On the other hand, pChk1 induction at higher concentration of Dox became less independent on ATM activity unlike pChk2 or at low concentration of Dox. pChk1 showed ATM independent induction at a higher Dox concentration starting at 2 hr of treatment and remaining > 6 fold high from 8-16 hr of treatment (Fig. 9a). Interestingly, these kinetics mirrored those observed for pATR in 400nM Dox + KU treatments and furthermore, were tightly correlated with cell death. For example, at 2, 4 and 8 hr of 400nM Dox + KU, there was induction of pATR and pChk1 accompanied by higher cell death. By 16, 20 and 24 hr of treatment, both pATR and pChk1 levels reduced with a corresponding reduction in comparative cell deaths. Since pChk1 is a direct substrate of pATR, these results point towards an ATM independent and ATR dependent apoptotic induction via Dox + KU = ATR →pChk1 →pP53 →Apoptosis.

Another candidate for pChk1 activation could be DNA-PK, which has been shown to have a functional interaction with it [96] and thus pChk1 induction at higher scale of damage in ATM inhibited state may demonstrate an overlapping function of DNA-PK which may still ensure DDR activation [97]. Overall, both in ATM active and inhibited states, pChk1 induction fitted well with the degree of cell death in such treatments. Hence, while lower scale of DNA damage, pChk1 was classified as effecting the cytoprotective function of ATM, at higher scale of DNA damage, it appeared to promote cytotoxicity both in ATM de- pendent and independent manner. Owing to a rapid induction of pATM and pChk2 following 400nM Dox treatment, it was expected to induce pP53 S-15 levels. Figure 9b shows rapid and very high pP53 induction around 9 fold at the first time point of 2 hr of treatment with time dependent increase and maximum at 24 hr time point (> 20 fold induction). Hence, so far, the higher degree of cell death with single agent Dox treatment is proposed to have an underlying DDR pathway involving ATM →Chk2/Chk1 →P53 →Apoptosis. Addition of ATM inhibitor with 400nM Dox disrupted pP53 induction only showing a slight transient rise at 12 hr time point. Nevertheless, pP53 levels in Dox + KU treated cells always remained below those in Dox treated cells, unlike what was observed in (Fig. 6b) where at 24 hr time point, KU addition to 100nM Dox resulted in greater pP53 induction than at any time point of Dox treatment alone. E2F1 is already mentioned earlier to have a prominent role in ATM mediated apoptosis following DNA damage [95, 98]. Indeed, E2F1 was found to be induced in a time dependent manner and showed greater than 12 fold induction at 20 hr post 400nM Dox treatment (Fig. 9b). Interestingly, while a lower scale damage of 100nM Dox resulted in ATM dependent reduction in E2F1 levels (proposed earlier via the pathways ATM →P53 →RB ┤E2F1 →Cell cycle arrest and ATM ┤ MDM2 — E2F1, Fig. 6b), which were abrogated after ATM inhibition, at a higher scale damage of 400nM Dox, E2F1 on the contrary showed ATM dependent induction. This clearly demonstrated dual fates of E2F1 during lower and higher scale damage, being degraded at lower scale damage to prevent cell cycle progression in order to arrest and repair cells, hence ensuring cell survival, while activated at higher scale damage, to promote ATM dependent apoptosis. The proposed pathway for this apoptotic mechanism could be ATM →(directly or via pChk2) →E2F1→p300/ARF →P53 or P73 →Apoptosis. These results reemphasized the importance of ATM function in performing context dependent multiple cellular roles, which in this instance was DNA damage dependent differential regulation of E2F1.

γ-H2AX S-139 showed a greater scale of time dependent induction during 400nM treatments as compared to 100nM demonstrating differences in scale of DNA damage incurred (compare Fig. 6b and 9b). These levels were completely inhibited with the addition of KU showing ATM dependence at most of the time points tested. However, at 20 and 24 hr of Dox and KU treatment, there was slight upregulation of γ-H2AX. pBRCA1 S-1524 showed a time dependent induction in its levels following treatment with 400nM Dox while disruption of this induction during ATM inhibition. This again demonstrated a sole requirement of ATM function during these DNA damage dependent inductions. While BRCA1 phosphorylation has been reported to occur in an overlapping manner by other DDR kinases [92], the antibody used in the current research detects S-1524 which was demonstrated to get completely disrupted after ATM inhibition.

Taken together, these data indicated a role of ATM in modulating both cell survival and apoptosis and that this decision was based on the scale of DNA damage and involves re-orchestration of underlying DDR signalling pathway. Interestingly, we also found that ATM substrates not only show distinct sensitivities towards ATM kinase inhibition, but also that this sensitivity shows variation both in time dependent as well as context dependent manner. From these results, it could be concluded that at earlier time points or at lower scale of DNA damage, ATM function ensured cell survival by causing cell cycle arrest and DNA repair, while at later time points or higher scale of DNA damage, ATM mediated cell death as inhibition of its function caused higher survival, hence making it an effective strategy to sensitise cells only when combined with a particular state of DNA damage in these cell lines.

3.2 Reverse engineering of time series data: results

3.2.1 Data input

The network innference algorithms presented in [63, 65] and reintroduced in 2.2.3 and 2.2.4 require time series data inputs to be specified in regular time intervals. We first obtained small quantities of experimental time series data (i.e. with only 3 time points) obtained from biological experiments and modelled those stratified data to inform intermediate experimental design to produce additional data. We progressively increased the number of time points data acquired until we were able to obtain experimental data with time points 0-4-8-12-16-20-24hr and these were first preprocessed to produce heatmap representation of those data 10. We further modelled those data sets to produce in silico topological maps of networks of interactions for interpretation and further analysis. Time series data may be preprocessed to generate a heatmap table that represents evolutionary dynamics of individual components in colour-coded depictions (Fig. 10). Here, we applied the techniques described in subsections 2.2.3 and 2.2.4 to infer networks of interactions from the biological time series data acquired and presented in Figure Preprocessed Time Series Data.

3.2.2 Application of reverse engineering and mathematical modelling

We applied assumption-free reverse engineered methods to model time measurements of DDR signalling pathway (pATM, pATR, pP53, pChk1, pChk2, pBRCA, pE2F1 and pH2AX) in (Fig. 10) recorded at timepoints 0, 4, 8, 12, 16, 20, and 24 hr and constructed data consistent predictive models to inform biological protocols and experiments using the inference procedure presented in [63] and reintroduced in 2.2.3. The mathematical modelling results are presented in Figures 11 and 12 which produce the jacobian network model representations of the four (4) experimental data sets. Each jacobian result is obtained using either of the matrix-based methods presented in 2.2.3 and 2.2.4. We further ensured that the constructed network models are data consistent (i.e. are all capable of simulating historical data utilised) by using to independently simulate (predict) the dynamics of the systems from 0-24hr without providing any information about the systems. The results of the simulation and prediction assessment reveal that nearly all the predicted data points are accurate. Data plots of both the predicted data and original time series data reveal that the constructed models are generally consistent as evidenced in Figures 13, 14, 15, 16 that show the plots of both data. Predicted data that do not match exactly are often characteristic of time series data that contain duplicate rows or columns. Such unsimulable (unreproducible) data are regarded as containing linearly dependent data. The TRM [63] inference method always seeks to infer data consistent models and networks of interactions from data.

3.2.3 Description of in silico topological maps

Topological maps of the network models inferred from time series are used to represent possible the most basic direct or indirect interactions that may exist between the protein kinases of the DDR signalling pathway. Each topological map may be viewed as a meta-descriptive network that represents the net change in evolutionary signalling dynamics of each kinase of the DDR pathway within certain time intervals. Such extractable information are useful for system identification purposes. Also they may be used to inform new experimental design or conduct systems analysis of the complex biological systems the produce the data. The information generated from such representations often reveal new knowledge which are useful for the understanding of the biological pathways. As the derived topological maps are schematic diagrams that may hold information about regulated mechanisms of interaction of ATM and ATR with their downstream effectors, the edges (links) within the networks may be viewed as valid (direct or indirect) interactions (signals, signal cascade or signalling process) between the protein kinases. In these maps there are basically two fundamental links to be considered: blue lines; and red lines. The blue lines represent any positive interaction (e.g. gene activation, protein phosphorylation etc) while the red lines represent the opposite effects (e.g. inactivation, inhibition, dephosphorylation etc). In both of these representations we assume that both the source (initiator of a signal) and target (recipient of the signal) are identifiable. To avoid confusion, we suggest that the target proteins are always at the end of the arrows that represent the interactions, whether positive or negative.

3.2.4 Further reverse engineering of 0.4µM Dox+KU data

Since the derived 0.4µM Dox+KU network of interactions in Figure 12 does not fully reflect the oscillatory patterns observed in the predicted dynamics shown in Figure 16, it does seem that multiple time invariant jacobian representations may be more appropriate for characterising this data rather than a single repre- sentation. In order to infer biphasic (i.e. in two phases) time invariant networks of interactions from single time series data, two non overlapping subsets of the such a data set may be infer independently. Each inferred jacobian matrix represents the network of interactions of the system at each phase. To this aim we further divided the 0.4µM Dox+KU data in two subsets (i.e. 0-8hr and 12-24hr) and inferred a network of interactions and in silico topological map from each data segment. The results produced are presented in 17 and 18.

4 Analysis of the mathematical models of DDR

The computational models for DDR signalling pathway helped elucidate how the activities of different proteins in the complex signalling network following DNA damage are triggered with each effector having distinct dynamic activation profile. Furthermore, it was demonstrated how the different DDR substrates and effectors show variable sensitivity towards ATM kinase inhibition (some activated while some inhibited) as well as that for each substrate alone, this sensitivity showed context dependent alteration. The DDR model also aided in elucidating how the network topology correlates with cell cytotoxicity or protection and provided for sensitivity analysis. Understanding this mechanism is particularly important while designing inhibitory strategies targeting DDR proteins to achieve cell sensitisation. This is critical since a prerequisite of anti cancer strategy based on the above principle is the characterisation of relative contribution of target molecules in a particular node of signalling network in modulating cell fate decisions within a treatment regime. This alone requires integrative pathway analysis and understanding in a systematic manner.

The mechanism of the induced DDR signalling pathway was further dissected by manipulating the activity of ATM in order to examine the resulting alteration in DDR and help elucidate the role of ATM in DDR signalling following genotoxic challenge. This was done to clearly demonstrate and delineate the role of ATM between DNA repair or survival and cytotoxicity or apoptosis. Quantitative data of signalling dynamics of important nodes within the DDR pathway revealed complex signalling relationships of ATM with its downstream substrates that function in physiological processes including DNA repair, cell cycle arrest, cell cycle re-entry, proliferation and apoptosis. Based on the data of signalling kinetics and apoptotic cell death, a number of important biological queries were formulated. For example, to help infer better signal- response relationships, to perform sensitivity analysis of the signalling network based on manipulation of ATM signalling and to see whether different concentrations of the same chemotherapeutic agent would produce the same signalling profile or reveal unique signatures. Finally it was also aimed to decipher the altered topological network following ATM inhibition to identify whether the emergent new links within network topology forms the signalling basis for the altered response observed or whether it is a compensatory mechanism to overcome the lost ATM signal.

This model in the first instance revealed a series of complex dynamics and inter-relationships between the different signalling nodes and established overall correlations between their magnitudes and kinetics which were not obvious in the experimentally determined biological data. Systems analysis of time series data firstly revealed consistency with the biological data. Most importantly, following DNA damage at both scales of damage tested, an oscillatory pattern of induction was seen in all the DDR proteins (Figs. 13 and 14). Oscillation of components of intra-cellular signal transduction pathways is a property that denotes feedback loops in biological systems. These feedback loops cause fluctuation of the proteins locked in such reactions and are important feature in maintaining a dynamic equilibrium between different components within a signal transduction pathway. Such oscillations were previously found out to occur for DDR proteins in both biological experiments and computational modelling [99, 100]. Other studies have shown direct involvement of ATM and Chk2 kinase activities in such oscillations [101, 102].

Previous experimental and modelling approaches have also attempted to identify the physiological conse- quences of these protein oscillations on cell fate following DNA damage [103]. This is an important outcome of modelling that not only fits well with the known oscillatory pattern of induction for DDR proteins, but also advocates the importance of in-silico prediction of system dynamics that work on experimentally determined data which highlighted oscillatory pattern of induction not obvious in the biological data points alone. Furthermore, the oscillatory pattern of induction was found to be ATM dependent as inhibition of ATM activity abrogated oscillations (Figs. 15 and 16). Additionally, this abrogation was also seen during higher scale of DNA damage upon KU treatment. This illustrates the importance of oscillatory induction in different treatment settings as reported previously [100, 104]. In one study [101], a modelling approach was used that proposed the involvement of ATM kinase activity in downstream protein oscillations following DSBs. These results are in agreement with the proposed mechanism and point towards the requirement of ATM kinase function in maintaining oscillatory mode via recurrent initiation mechanism involving feedback loops.

Figure 11 represents map of network topology, which shows that ATM activity is induced at lower scale of DNA damage. pP53 S-15 was shown to be the immediate downstream substrate of ATM to be induced by ATM activity and hence was positively regulated. Since the cytotoxicity assay showed lower cell death at this point as compared to cells with a blocked ATM function, this up regulated pP53 may exert CIP/KIP mediated cell cycle arrest and promote ATM dependent repair of the damaged DNA as mentioned earlier.

Importantly, the topological map in figure 11 shows that the induced pP53 S-15 levels caused downregulation of E2F1 at lower scale of DNA damage. This link supports the signalling conclusions drawn earlier and is in agreement with the experimental data. Additionally, there is a minor negative link between pATM and E2F1, which also supports pATM mediated E2F1 inhibition, which would lead to cell cycle arrest. Furthermore, another computationally determined indication of cell cycle arrest is provided by the induction of pChk1 through pATM activity as well as revealing no direct link between pATM and pChk2. Previous study has reported the role of pChk1 in cell cycle arrest following Dox treatments [104], and ATM mediated activation of Chk2 that promotes apoptosis [105]. Hence, the activation of ATM at lower scale may indi- cate its role in cytostatic and cytoprotection against genotoxicity to promote DNA repair. The experimental data revealed no major change in pATR induction at 0.1µM Dox. Interestingly, the in silico topological map shows pP53 mediated inhibition of pATR in response to lower scale of Dox treatment. P53 has already been previously reported to downregulate ATM expression [106] as well as pChk1 [107]. A possible mechanism of the p53 dependent pATR inhibition seen in signalling network could be P53 mediated transcriptional upregulation of Cyclin G that is known to recruit PP2A [108] which in turn has recently been reported to antagonize ATR activity [109, 110]. Another inhibitory link could be provided by WIP1, a positive downstream target of P53 [111] shown to reverse the ATR mediated DDR pathway [112]. Furthermore, the negative influence of pP53 on pChk2 and γ-H2AX in the topological map could also be explained via suppressive effects of WIP1 on pChk2 [113] and γ-H2AX [114]. Overall, these results are consistent with biological data where while both pATM and pP53 showed time dependent induction, pATR levels were mostly kept in check at all the time points tested. Altogether, these links indicate that at a lower scale of DNA damage, while pP53 levels are induced in ATM dependent manner, the cell may avoid apoptosis by keeping in check the activities of pChk2, pChk1 and pATR and exert cell cycle arrest and DNA repair. Activation of the DNA repair component of ATM signalling is also supported by positive links between pBRCA1 and pP53 as well as pBRCA1 and pATM. BRCA1 has been shown to act as a co-activator of P53 [115] and to direct P53 function towards cell cycle arrest and DNA repair by upregulating p21 and GADD45 [116].

Functional assumptions of pATM signalling from topological network map are further strengthened by the corresponding changes in the signalling network upon ATM inhibition (Fig. 11 lower panel). Most prominently, the activation of pP53 by pATM was diminished. The negative links between pP53 to pChk1, pChk2, E2F1 and γ-H2AX switched to positive influences. Mathematically derived topological maps uncovered the underlying signalling profiles associated with altered cytotoxicity responses seen in different treatment regimens. The altered network of signalling in the ATM active and inhibited states explained the associated physiological response and provided a mechanistic overview that aided in establishing signal to response relationship. Experimentally determined DDR dynamics showed a later induction of pATR and E2F1 levels in response to ATM inhibition following Dox treatment. An interesting observation was the emergence of a strong positive signal of E2F1 towards pATR upon ATM inhibition. This computationally determined important new signal in the network illustrates a possible mechanism of greater apoptosis seen in ATM inhibited states. Secondly, in the presumed role of ATM in cytoprotection at lower scale of damage where ATM-induced pP53 negatively influenced pATR levels, KU treatment resulted in the loss of this negative influence. Thus, pATR upregulation may result from two separate mechanisms. One is loss of suppressive effect of p53 and the other is upregulation via E2F1.

Altogether, signals following lower concentration of Dox treatment suggest that while functional ATM keeps E2F1 activity in check, probably causing cell cycle arrest and time for DNA repair, disruption of ATM kinase ceases E2F1 sequestration that may directly upregulate ATR while disruption of suppressive effect of pP53 towards pATR would further triggers apoptotic signalling. ATR activity has already been known to signal apoptosis in response to treatments with different chemotherapeutic agents [78, 79]. Hence the inhibitory state of ATR levels at lower DNA damage suggests cellular signalling preference towards DNA repair.

As mentioned before, KU treatment caused greater sensitivity to 0.4µM Dox only at initial time points of 2, 4 and 8 hr treatments while decreased sensitivity following 12, 16, 20 and 24 hr of treatment (Fig. 8) probably indicating a functional switch and biphasic ATM activity. The dual consequence of ATM kinase inhibition could be explained by previous reports illustrating ATM functions both in cell cycle arrest and DNA repair whereby it promotes cytoprotection 17, 22), as well as functions as a central component in triggering apoptotic cell death [27, 28, 30]. This indicates that its specific downstream signalling preference is context dependent, and an alteration to its substrate preference may occur as a function of time of treatment or dosage of drug.

If indeed such explanation holds for the observed dual outcomes of ATM inhibition, both the underlying DDR signalling kinetics as well as signalling connections within DDR network are expected to alter during this transition of ATM function. It is important to appreciate that the signalling alterations detected during ATM inhibition can either be causative of such altered phenotypic response, or be activated as an emergency response in an attempt to over come the lost ATM signal.

Since at lower scale of DNA damage, co-treatment with KU sensitised cells at all the time points of treatment indicating a sustained cytoprotective function of ATM, a single overall in silico network of topology was established. On the other hand, at higher concentration of Dox, KU treatment caused specific time dependent sensitisation and desensitisation, due to which two independent computationally determined topological networks were established between 0, 2, 4 and 8 hr of time point (cell sensitization following KU) (Fig. 17) and 12, 16, 20 and 24 hr time point (cell desensitisation with KU) (Fig. 18).

Analysis of the topological networks revealed that E2F1 was strongly upregulated by ATM, both between 0-8 hr and 12-24 hr of 0.4 µM Dox treatment. However, there was a stronger influence of ATM on E2F1 in 12-24 hr network as compared to 0-8 hr. Upregulation of E2F1 via ATM is well documented to signal P53 dependent and independent apoptosis. pChk1 was found to be positively up regulated by ATM both between 0-8 hr and 12-24 hr networks. However again, there was a stronger influence of ATM on pChk1 in the latter network. A Switch to greater Chk1 induction thus may trigger apoptosis. Another important signalling alteration between the two networks was observed between BRCA1 and E2F1. ATM was shown to greatly upregulate BRCA1 in the 0-8 hr network that was further shown to negatively regulate E2F1. Similarly, ATR was shown to upregulate BRCA1 in the early network. BRCA1 activation as mentioned may elicit cytoprotective pathways. In the 12-24 hr network on the contrary, this pATM → pBRCA1 signal diminished and pBRCA1 instead upregulated E2F1. Furthermore, ATR signalling switched to downregulate pBRCA1. These links overall may result in E2F1 accumulation in the cell, providing underlying signalling mechanism to the experimentally observed upregulation of E2F1 indicating a trigger of ATM dependent apoptosis.

E2F1, the downstream substrate of ATM was shown to be positively influenced in both the early and later network of 0.4µM Dox (as opposed to being down regulated during 0.1 µM Dox). Once activated, while E2F1 showed negative signalling towards ATR in the earlier topological network, in 12-24 hr network, this negative link was disrupted. This may represent a mechanism for ATR accumulation at higher scale of damage, which again might cause chk1 mediated apoptotic signal as mentioned earlier.

An interesting signalling alteration was seen between ATM and Chk1. ATM was shown to upregulate pChk1 both in the 0-8 hr (Fig. 17) and 12-24 hr (Fig. 18) networks. Strikingly, this upregulated pChk1 was shown to be locked in a positive feedback loop with BRCA1 (which was also directly activated by ATM) in earlier network, while disruption of this and emergence of a new positive feedback loop with E2F1 (which was then also directly activated by ATM) in the later 12-24 hr network. This important signalling alteration may be a hallmark of switch of ATM from a role in cytoprotection to a role in cytotoxicity mediated by Chk1.

In terms of ATM and Chk2 signalling, in earlier network, pATM upregulated pChk2 that further formed a similar positive feedback loop with BRCA1 while inhibiting ATR. However, in the later network, this positive loop was shown to be disrupted and pChk2 switched its signal to a positive influence on ATR that further inhibited BRCA1.

Extensive data analysis of computationally determined topological maps between time points 0-8 hr and 12-24 hr have thus successfully identified the above mentioned critically important signalling alterations.

The underlying signalling changes that the in silico modeling has identified not only explained the difference in cellular response to ATM inhibition between the two sets of treatments, but also invited for further experimental testing. For example, the altered link between ATR and BRCA1, feedback loops involving Chk1 and BRCA1 in the earlier network and pChk1 and E2F1 with the disruption of the first feedback loop in the later network, the switch of E2F1 and pChk2 influences on ATR from negative (0-8 hr network) to positives (in 12-24h network) demonstrate novel computationally derived insights that could be targeted to manipulate cellular response to treatment regimens.

5 Discussion

Decades of research efforts in attempts to understand the DDR pathway have added significant knowledge to our comprehension of the mechanisms involved in DNA damage and repair and induction of apoptosis [117]. Researchers have identified that cellular sensitivity to genotoxic agents in course of cancer therapy could be achieved by modulating the function of key proteins involved in DNA repair, cell cycle arrest and apoptotic mechanisms [118]. The role of ATM in the above cellular processes is of immense importance and thus has been the centre of attention in the past decade or so. ATM has been targeted by way of either inhibiting kinase activity [119, 120] its expression [121, 122], its trafficking [9] or overexpression of its dominant negative form [123] [For review, see 2]. However, due to the central role of ATM/ATR pathway in governing myriads of other molecular interactions, it is often difficult to devise inhibitory treatment regimens allowing the predictability of both sensitivity and selectivity.

A prerequisite of a molecularly targeted anticancer approach is a detailed understanding of the underlying abhorrent signalling network contributing to tumour development [124]. This involves changes in both the types and extents of molecular interactions governing key process in cellular homeostasis. It is now known that numerous key proteins have multiple functions in different pathways and that the downstream signalling choice of a particular protein depends on a number of variables and is context dependent. In this micro- environment, the effect of inhibiting ATM protein is not fully predictable in terms of signalling consequences and the associated cellular response, especially owing to the fact that ATM function can contribute in both cell cycle arrest to allow for DNA repair [24, 25], as well as in apoptosis [125, 30]. Prediction is further complicated by the finding that ATM may self-regulate its own protein levels [75, 117]. In this complex and unpredictable situation, the efficacy of potential ATM inhibitors would have to be assessed in terms of the effects it exerts, not only on ATM activity and its immediate substrates e.g. ATM → pATM → pP53 → DNA repair, but on a number of connected key proteins in the DNA repair, cell cycle and apoptotic pathways and the overall impact of such pathways on cell fate [2]. Therefore, it is increasingly becoming obvious that complexity in determining the absolute effects of kinase inhibition necessitates parallel theoretical and experimental considerations. This could be best addressed by employing a systems biology approach to the problem, involving mathematical modelling of biological processes [73, 118].

Over the years, tremendous progress has been in producing qualitative signalling data that provide for the explanation for different cellular events. In terms of therapeutic intervention of such pathways, a mere descriptive knowledge is not enough. The decision making property of these pathways may involve oscillations and concentration thresholds that require a quantitative approach, calculations and numeral analysis for data collection and interpretation [126]. Once such time series data is available for use in mathematical logic, they can not only shed light on the mechanism of functioning of these vital cellular pathways but could also be predictive in nature in terms of outcomes of interventions.

Following treatments both at high and lower scale of DNA damage, with and without ATM inhibition, attempts were made to correlate status of individually examined important nodes in the DDR and explain how they interact with each other to bring about the phenotypic changes in the cell leading to cell death and thus to establish relationships between signalling network to the cellular response. To do this, we first determined the phosphorylational events in DDR following different treatments. We then quantified their distinct kinetics through time courses and finally performed NR-uptake based cell cytotoxicity assays. This allowed parallel study of temporal induction of DDR phosphorylational signalling and their associated effects on cell health and survival. Since ATM inhibition was also used along with the drug treatments, the precise role of ATM in either DNA repair or cell death was also inferred. This enabled the interpretation of pathway activation on overall cell fate. Hence, experimentally we adopt a method of drug interventions, one which is based on varying levels of drug dosages to stratify experimental trials to provided rich quantitative data for information extraction purposes and knowledge discovery. The novel mathematical and computational modelling strategies assessed previously using artificially controlled in silico experiments [127] are now applied to real biological data. In this way, the same strategies that worked during in silico experimentation are now applied here in time series data mining to identify potential treatment regimes or inform the design of new therapeutic targets and alternatives in cancer treatments.

There is a need for collecting quantitative dynamic information of DDR signalling proteins at high temporal resolution [128, 129]. Mathematical model construction based on such information would not only provide novel insights into how thresholds, localisation and specific interactions of proteins within DDR signalling are triggered and regulated during the course of genotoxicity, but would also help establish maps of network interactions that would further aid in identifying spatio-temporally regulated critical links and their contributions in pathways responsible for generating a specific cellular response during a treatment regime. Furthermore, such a systems approach would elucidate how signalling links within a network of interaction influence overall cellular behaviour in a given type of treatment condition. Once such critical signalling links are established in cancer cells, these could be exploited to devise treatment portfolios to achieve targeted cellular sensitivity. In the past, several researchers have targeted a variety of dynamical processes within biological systems for their computational analysis ranging from studies into parallel reaction pathways [130], HBV infections [131] and in silico analysis of Sar-CoV [132]. In terms of DDR signalling pathway, most of systems biology studies have been undertaken to examine and model the oscillatory patterns of DDR protein induction following DNA damage in light of the known biological insights [133, 134, 135, 137] with some reports focussing on elucidating the role of such oscillations in determining cell fate i.e. DNA repair, cell cycle arrest or apoptosis [88, 100, 103, 138]. However, owing to an intricate nature of context dependent signalling networks, high degree of cross-talk and pathway choices and insufficient quantitative data of signalling dynamics, a clear inter-relationship among cellular signalling to cellular response to its ultimate consequence on cellular phenotype has not been fully determined.

The current work attempted to establish a link between signalling networks emerging from ATM that functions at the core of DDR pathway upon double stranded DNA damage, to its ultimate impact on cell survival as well as to determine a cellular fail-safe mechanism upon ATM inhibition and how the cell responds by re-arranging signalling links to adapt to such inhibition and maintain the same phenotypic consequence. Previous mathematical modelling attempts have shown that the design and construction of a deterministic mathematical model of the molecular interactions that underpin the DDR require new kind of time series quantitative data [12]. This data must be consistent, have high temporal resolution and spatial consideration and provide the kinetic parameters of larger number of different proteins involved in a single pathway. While qualitative data are easily available, quantitative data pertaining to key signalling molecules that would allow speedy calibrations and provide kinetic parameters for the construction of mathematical model is scarce.

Assuming that quantitative data provide a representation of the unknown dynamics of signal activation of divergent proteins involved in DDR pathway is available, it follows that such data may be analysed to generate new and useful hypotheses about the possible or potential biological signalling, if done carefully. Furthermore, we demonstrate that by using simple and practicable methods informative and data-consistent network model may be inferred and constructed from such experimental time series. The constructed predictive model is then used to populate the original data set with better and richer quantitative data set for studying and understanding system dynamics and behavioural patterns.

Successful computational modelling of molecular pathways triggered by multifunctional proteins like ATM can not only provide predictive behaviour of this kinase in a specific treatment regimen but also highlight key differences in pathway choices between normal and cancer cells within a single treatment regime. These predictions could then be tested at experimental level to examine the degree of accuracy of these models. Thus, tremendous assistance could be provided in researching potential drug treatment options for certain types of cancers based on in silico predictions. Following this, and if other potential ways are found for inhibition of the DDR pathway, a secondary calibration stage could be carried out along with the initial clinical trials. This would permit the model to be developed as a clinical application to be used in conjunction with the drug treatment regimes. Such an application would enable clinicians to determine the optimum timing and level of drug dose and radiation for use in the individual patient where the DDR pathway status has already been established by the measurement of key protein levels in that patient. In the terms of this approach the optimal outcome is taken to be the highest number of cancer cells damaged by the drug and/or radiation, which are directed toward apoptosis.