A Mathematical Model for the Detection Mechanism of DNA Double-Strand Breaks Depending on Autophosphorylation of ATM

Background After IR stress, DNA double-strand breaks (DSBs) occur and repair proteins (RPs) bind to them, generating DSB-RP complexes (DSBCs), which results in repaired DSBs (RDSBs). In recent experimental studies, it is suggested that the ATM proteins detect these DNA lesions depending on the autophosphorylation of ATM which exists as a dimer before phosphorylation. Interestingly, the ATM proteins can work as a sensor for a small number of DSBs (approximately 18 DSBs in a cell after exposure to IR). Thus the ATM proteins amplify the small input signals based on the phosphorylation of the ATM dimer proteins. The true DSB-detection mechanism depending on ATM autophosphorylation has yet to be clarified. Methodology/Principal Findings We propose a mathematical model for the detection mechanism of DSBs by ATM. Our model includes both a DSB-repair mechanism and an ATM-phosphorylation mechanism. We model the former mechanism as a stochastic process, and obtain theoretical mean values of DSBs and DSBCs. In the latter mechanism, it is known that ATM autophosphorylates itself, and we find that the autophosphorylation induces bifurcation of the phosphorylated ATM (ATM*). The bifurcation diagram depends on the total concentration of ATM, which makes three types of steady state diagrams of ATM*: monostable, reversible bistable, and irreversible bistable. Bistability exists depending on the Hill coefficient in the equation of ATM autophosphorylation, and it emerges as the total concentration of ATM increases. Combining these two mechanisms, we find that ATM* exhibits switch-like behaviour in the presence of bistability, and the detection time after DNA damage decreases when the total concentration of ATM increases. Conclusions/Significance This work provides a mathematical model that explains the DSB-detection mechanism depending on ATM autophosphorylation. These results indicate that positive auto-regulation works both as a sensor and amplifier of small input signals.


Introduction
Recently, research suggests that biological functions depend on specific small components called network motifs [1]. In these motifs, positive and negative feedbacks are very important for bistability or oscillatory behaviours, respectively. For example, positive feedback in mitogen-activated protein kinase (MAPK) cascade produces bistability of phosphorylated MAPK which contributes to an all-ornone cell fate switch [2,3], and the production of self-sustained biochemical ''memories'' of transient stimuli [4]. A synthetic regulatory network of a mutually inhibitory double negative feedback loop in Escherichia coli also provides bistability, and a simple theory that predicts the conditions necessary for bistability has been suggested [5]. Also, stochasticity of gene expression in a single cell has been recently observed [6]. These stochastic single-molecule events determine a cell's phenotype depending on positive feedback [7,8,9]. However, understandings of functions for these positive feedbacks are limited.
Generally, there are several factors which damage DNA in cells. Signal-transduction pathways are rapidly activated after exposure to DNA-damaging agents. ATM, the gene that is mutated in the human disease ataxia-telangiectasia (AT), is important for activating signalling pathways in mammalian cells following exposure to ionizing radiation (IR) or oxidative stress where DNA double strand breaks (DSBs) are generated [10]. For example, hydrogen peroxide (H 2 O 2 ), one type of oxidative stress, is a normal metabolite in the cell whose steady-state concentration is in the range 10 28 -10 29 M [11,12], and is one of the products to protect the mammalian host from the invasion of bacillus [13]. However, if it is not properly controlled, it can cause severe damage to a cell. In the presence of Fe 2+ , H 2 O 2 can generate free radical (OH N ). Also, the Haber-Weiss reaction can form OH N in an interaction between O .{ 2 and H 2 O 2 in the presence of Fe 2+ or Fe 3+ [12]. These oxygen free radicals and H 2 O 2 are spoken of as reactive oxygen species (ROS). DNA strand breaks are due to free radicals in these reactions. Depending on its concentration, H 2 O 2 induces two types of DNA lesions: DNA single-and double-strand breaks. DNA single-strand breaks (SSBs) are dominant under H 2 O 2 stress, but these lesions are efficiently repaired and do not appear to mediate the cytotoxic response [14]. On the other hand, DNA double-strand breaks (DSBs) seldom occur under H 2 O 2 stress, but they are toxic and potentially induce apoptosis. In addition, after IR stress, DNA double-strand breaks are generated and ATM proteins are phosphorylated. The true mechanism of this process has not been understood yet. However, one of the targets of ATM phosphorylation is suggested to be the Nbs1 (nibrin) protein, which associates with the conserved DSB repair factors Mre11 and Rad50 [15]. The phosphorylated ATM phosphorylates itself, and it is suggested that ATM is autophosphorylated within 15 min after exposure to 0.5 Gy IR, which induces only 18 DNA breaks, approximately [16,17]. However, it has not been known how ATM detects a small number of DSBs and activates signalling cascades. In this paper, we propose a mathematical model of the ATM phosphorylation process after a stochastic generation of a small number of DSBs regardless of the source of damage.
In our model, we assume that DSBs are generated under DNA damage, repair proteins (RPs) bind to them and become DSB-RP complexes (DSBCs), and DSBs are repaired. The numbers of DSBs and DSBCs are very small and these processes are stochastic. We will see that we can calculate theoretical values for the mean numbers of DSBs and DSBCs by using the number of repair proteins and rate constants of repair processes. The produced DSBs and DSBCs phosphorylate ATM which autophosphorylates itself. We will find that DSBs are not successfully repaired and the number of DSBs increases when the number of repair proteins is small, but when sufficient repair proteins exist, the number of DSBs is suppressed to low levels. Also, we will find that autophosphorylation of ATM induces bifurcation of the phosphorylated ATM (ATM*). Depending on the total concentration of ATM, the fixed points of ATM* will have three types of steady state diagrams: monostable, reversible bistable, and irreversible bistable diagrams. Of these steady-state diagrams, bistability emerges when the total concentration of ATM increases, and the concentration of ATM* exhibits switch-like behaviour in the presence of such bistabilities. Furthermore, we will see that the time to detection after the DNA damage decreases when the total concentration of ATM increases. Figure 1 shows a diagram of our model. DSBs are induced by some stress signal, and repair proteins bind to DSBs, generating DSBCs. The repaired complex produces a repaired DSB. We model these processes as stochastic processes. Details are shown in the next section. The generated DSBs and DSBCs are detected by ATM, and it phosphorylates itself to amplify the stress signal. In this paper, we mainly focus on the repair processes of DSBs and the detection mechanism of ATM. The negative feedback from p53 is also treated in the Discussion section.

Calculation methods for DSB-production model
In this paper, we assume the following schemes of a stochastic production mechanism of DSBs and their repair processes: where DSB denotes DNA double strand breaks, RP denotes repair proteins, DSBC denotes DSB and repair protein complexes, and RDSB denotes the repaired DSB. The constants, c 1 , c z 2 , c { 2 , c 3 , represent the stochastic rate constants [18] (or reaction parameters [19]). Also, the associated rate laws (or hazard functions) are h i X ,c i ð Þ, where i is a reaction type and X~X DSB ,X RP ,X DSBC ,X RDSB ð Þ is the current state (the number of molecules (or sites) of each reaction species) of the system. These chemical reactions occur stochastically, thus the fluctuations of the number of molecules which are produced in these reactions are stochastic processes [18]. For example, the production of DSBs is a zeroth-order reaction, and the hazard of the reaction is The repair process of the DSBs is a second-order reaction, and the combined hazard of the reaction is The failed and succeeded repair processes are first-order reactions, and we respectively denote the combined hazards of each reaction as The above equations allow us to calculate X i for each molecular type i by using the Gillespie algorithm (we use Gillespie's direct method [20]). For example, a time course of X DSB zX DSBC is shown in

Comparison of theoretical and simulation results
In this section, we compare simulation and theoretical results of mean values of DSBs and DSBCs. Methods for calculating the theoretical mean values of DSBs and DSBCs are shown in the Materials and Methods section. The precise mechanism of DNA damage processes induced by stress signals has not been clarified yet, and we cannot estimate the stochastic rate constants. In this paper, we assume that the stochastic rate constants c z 2 , c { 2 , and c 3 are not affected by stress signals, and their values are defined in Table 1. In addition, we assume the parameter c 1 is proportional to the strength of stress signals: For example, when we define the parameters c 1~1 0, and X max RP~1 00, we can estimate the mean numbers of DSBs and DSBCs (X X DSB andX X DSBC ). Then the time (t DSB C ð Þ ) until the X DSB C ð Þ approaches the mean aŝ t DSB *0:0525: ð10Þ Figure 3 shows the comparison between simulation and theoretical results. The E X i t ð Þ ½ denotes the ensemble average of X i at time t. In Figure 3A, we can see that E X DSB t ð Þ ½ converges on the steady state, but its value is a little different from the theoretical value. The time until the ensemble approaches the steady state (,2.0 min) is larger than the theoretical value (t DSB ). In   Figure 3B, E X DSBC t ð Þ ½ also converges to the steady state, and the value is approximately the same as the theoretical value. However, the time until the ensemble approaches the steady state is approximately 10.0 min. This is larger than the theoretical value (t DSBC ). Also, we compare the theoretical and simulation values of the mean numbers of DSBs and DSBCs as a function of the maximum number of repair proteins (X max RP ) at time t = 100 [min]. For a small number of repair proteins (30ƒX max RP ƒ100), theoretical values can explain the simulation values well. However, as the maximum number of repair proteins increases, the difference between the theoretical and simulation results increases for both E X DSB ½ and E X DSBC ½ .
Small differences in the maximum number of repair proteins induce a large difference in DSB generation The time course of the model depends largely on the maximum number of repair proteins, X max RP . Figure 4 shows the time course of X DSB and X DSBC with c 1~1 0. If X max RP wX X DSBC~2 0 ð Þ, the number of DSBs is very small, and the number of DSBCs approachesX X DSBC~2 0. However, the number of DSBCs existing at the same time should be smaller than the maximum number of repair proteins, X max RP . Thus, if the maximum number of repair proteins is smaller than the mean value of X DSBC , the production rate of DSBs becomes higher than the rate of the DSB-repair processes, and the free repair proteins are given out to repair DSBs; as a result, the number of DSBs increases. Figure 4 shows the time courses of X DSB C ð Þ with several maximum numbers of repair proteins. When X max RP~2 0, the number of DSBCs is approximately 20, and it cannot exceed 20. Then the number of DSBs gradually increases to nearly 100. On the other hand, when X max RP is larger than 30, X DSB does not continuously increase. As the number of repair proteins X max RP increases, the number of DSBs gradually decreases (Figures 4B, C, and D). Therefore, a time course of the ensemble average of X DSB zX DSBC gradually increases for a small number of repair proteins X max RP~2 0, but it converges to a constant value when the number of repair proteins is large ( Figures 4E and F). In the following sections, we assume the number of repair proteins is sufficient, and the expected value of X DSB zX DSBC is constant.

Autophosphorylation of ATM induces bifurcation
The ATM module is modeled as follows: where the total concentration of ATM satisfies , the numbers of DSBs and DSBCs satisfy X DSB zX DSBC ð Þc Stress ½ , and other rate constants are as described in Table 1. The right-first term denotes autophosphorylation of ATM by ATM* (details are shown in the Materials and Methods section). The true value of the Hill coefficient n A is unknown, but the steady state features of ATM* depend on this value. We assume that this value is n A~2 (details of the reasons are described in the Materials and Methods section), which can induce bistability of ATM* as we show in Figure 5. The term indicates that ATM is phosphorylated by DSBs and DSBCs. Here X i and [ATM] have different dimensions, and therefore k DSB includes a role to adapt the dimensions. If we plot the fixed points above the (j ATM ,k DSB ) plane, we get the cusp catastrophe surface [21] shown in Figure 5A. This figure shows the steady state concentrations of ATM* and the stability of them. For some parameter regions, there are three fixed points (one is unstable, but the other two are stable). We call these regions bistable regions. We show the parameter regions of three fixed points in Figure 5B. In this figure, we show two bifurcation curves and they meet at a point j ATM ,k DSB ð Þ 0:825,0:0029 ð Þ , which we call a cusp point. Next, we fix the value of a parameter j ATM or k DSB , and plot the fixed points as a function of j ATM or k DSB . In Figure 5C, we show the fixed points as a function of k DSB with constant j ATM values. As we can estimate from Figure 5B, there exists bistability when the parameter j ATM satisfies the condition 0ƒj ATM ƒ0:825. Specifically, we select four values j ATM~0 :0,0:65,0:78, and 0.85. When j ATM~0 , there is a bistable region, but it dominates a very narrow range of k DSB . When j ATM~0 :65, or 0.78, the steady state concentrations of ATM* have bistability. One is irreversible (0.65), and the other is reversible (0.78). For the other value of the parameters, j ATM~0 :85, it is monostable as we can expect from Figure 5B. Also, in Figure 5D, we show the fixed points of ATM* as a function of j ATM with constant k DSB values. As Figure 5B suggests, bistability exists in the condition that 0ƒk DSB ƒ0:0029. When k DSB~0 :0 or 0:002, the steady state concentration of ATM* has bistability, but for the other two parameter values it has monostability as we can estimate from Figure 5B. Therefore, the ATM model generates bistability when it contains an autophosphorylation mechanism with the Hill coefficient of n A~2 mentioned above.

The total concentration of ATM determines the bifurcation diagram
Here we consider the effect of the total concentration of ATM on the bifurcation diagram. The bistable regions with the rate constants k DSB and j ATM are shown in Figure 6A. When the total concentration of ATM increases, the bistable region expands. For a certain rate constant pair (k DSB~0 :0023 and j ATM~0 :78), we plot the fixed points as a function of X DSB zX DSBC ( Figure 6B). The steady state concentration of phosphorylated ATM exhibits irreversible bistability when the total concentration of ATM satisfies ATM tot ½ 3:0 (mM), reversible bistability when ATM tot ½ 2:0 (mM), or monostability when ATM tot ½ 1:0 (mM). Therefore, the total concentration of ATM determines the bifurcation diagram. Figure 6C is the bistable region for X DSB zX DSBC and ATM tot ½ . When ATM tot ½ is higher than 1.89 and lower than 2.19, the model exhibits reversible bistability. When ATM tot ½ is higher than 2.19, the model exhibits irreversible bistability. In the presence of noise, characteristics of the time-dependent concentration of ATM for the three cases become different. Here the term noise means the repair-process noise for the DSB-production model in equation (1).

Significance of ATM bifurcation as a DSB sensor
Here we consider the stochastic case where DSBs are generated stochastically, which means X DSB zX DSBC ð Þ is a stochastic process. To connect these stochastic processes and the deterministic ATM model, we used a hybrid simulation algorithm (see Materials and Methods and reference [22]). We define the parameters of DSBs and DSBCs as in Table 1. In addition, the parameters of the ATM model are k ATM1~kATM2~1 :0. Furthermore, the initial condition of [ATM*] is 0. In this case, the concentration of phosphorylated ATM becomes a stochastic process. Figure 7 shows examples of the time courses of the concentration of phosphorylated ATM with several parameters of ATM tot ½ . In the deterministic case, we calculate the steady state [ATM*] for the constant X DSB zX DSBC ð Þ , but in the stochastic case X DSB and X DSBC become stochastic processes. Therefore the concentration of [ATM*] fluctuates because of the fluctuations of DSBs and DSBCs. When the total concentration of ATM is small ( ATM tot ½ 1:89 (mM)), the concentration of ATM* is suppressed to a low value. When ATM tot ½ 2:0 (mM), transitions of the concentration of ATM* between low and high values occur. When ATM tot ½ 2:19 (mM), a transition of the concentration of ATM* from low to high values occurs. These three cases reflect the bifurcation diagrams of deterministic cases. For example, when ATM tot ½ is smaller than 1.89, the fixed points of ATM* are monostable and ATM* fluctuates around a low concentration. However, when 1:89ƒ ATM tot ½ ƒ2:19, the fixed points of ATM* are reversible and bistable. In this case, once the numbers of DSBs and DSBCs cross over the bifurcation point, ATM* jumps to the higher concentration state. This high concentration state can return to the lower state because of the reversibility of ATM*. When ATM tot ½ w2:19, once ATM* jumps to the higher state, ATM* cannot return to the lower state because of the irreversibility of ATM*.
Here we define the transition threshold T ATMÃ as the value that [ATM*] passes through if transition occurs. When the concentration of ATM* passes the threshold, we assume the DSBs and DSBCs are detected by the ATM module. In Figure 8, we show the detection rate as a function of time. This figure shows that if the total concentration of ATM increases, the detection rate increases, and the response becomes quick. This trend can also be explained from the bifurcation diagrams of the deterministic model. When the total concentration of ATM is between 1.89 and 2.19, ATM* becomes reversible and bistable. As the total concentration increases, the bifurcation point of X DSB zX DSBC decreases as shown in Figure 6B. Therefore, when the total concentration of ATM is high, ATM* jumps to the higher concentration state even when the numbers of DSBs and DSBCs are small. Thus the detection rate increases.

Discussion
In this work, we modeled the stochastic repair processes of DSBs and a detection mechanism which is based on the autophosphorylation of ATM. In our first model, we could simulate time-dependent fluctuations of the numbers of DSBs and DSBCs, and we proposed theoretical mean values of DSBs and DSBCs. Depending on the species and strength of stress signals, DSBs may rarely occur. Even in this case, it is suggested that ATM can detect DSBs [16]. We propose that ATM can detect a small number of DSBs by using an ATM autophosphorylation mechanism, which induces bifurcation of ATM*. In the presence of bistability, ATM* exhibited switch-like behavior. Also, we suggested that the total concentration of ATM determines the bifurcation diagrams, and as the total concentration of ATM increases, the detection rate also increases. Therefore, we conclude that the positive auto-regulation works as a sensor of small fluctuating DSBs and amplifier for the detected signals.

A theoretical method for determining stochastic rate constants
In experiments, dynamics of repair processes of DSBs are still unknown. However, based on the experimental result in which ATM can detect about 18 DSBs in a cell, we defined the parameter values of c 1 , c z 2 , c { 2 , and c 3 which induces approximately 20 DSBs plus DSBCs. Then we considered how ATM detects these small numbers of DSBs and DSBCs. To define these four parameters, we described theoretical mean values of DSBs and DSBCs by using the four parameters (equations (19) and (20) in Materials and Methods). Interestingly, the mean number of DSBCs (X X DSBC ) only depends on the production rate of DSBs (c 1 ) and the success rate of DSBCs (c 3 ). On the other hand, the mean number of DSBs (X X DSB ) complexly depends on all rates and the maximum number of repair proteins (X max RP ) which is inversely proportional toX X DSB . We compared these theoretical results and simulation results. As expected,X X DSB gradually decreases as X max    (Figure 3). However, the difference between theoretical and simulation results forX X DSB andX X DSBC increases as X max RP increases. More specifically, the theoreticalX X DSB converged to 0, but the simulation result converged to 0.5. This may be because of the fact that in the simulation results X DSB and X DSBC are small and they are stochastic processes depending on hazard functions h m X ,c m À Á , which define the reaction rates. Therefore, mean-field approximation of theoretical results cannot evaluate the true mean numbers of DSBs and DSBCs. For example, X DSB usually becomes 0 or 1 when X max RP is large. In this situation, if X DSB~0 , the hazard function h 2 z X ,c 2 z ð Þ0, and repair processes are independent of the repair proteins. Therefore, DSBs can be generated regardless of the number of repair proteins. If X DSB~1 , the hazard function h 2 z X ,c 2 z ð Þ~c 2 z X RP ð Þbecomes very large, and a generated DSB is immediately repaired. In addition, when X max RP is very large, each process occurs at approximately the same rate. Thus the mean value of DSBs does not converge to 0 but some other value between 0 and 1.

RP
Many experiments should be done to confirm our results. We have to observe the numbers of DSBs and DSBCs as functions of time in a single cell. The number of repair proteins characterizes repair processes of DSBs and also should be observed. Even when these observations are successful, it is difficult to estimate the stochastic rate constants for repair processes. For example, when we get the time until the number of DSBCs approaches the mean t DSBC , we can estimate the stochastic rate constant c 3 *1=t DSBC (see Materials and Methods). Then we can also estimate the stochastic rate constant c 1 *X X DSBC t DSBC . The other stochastic rate constants c + 2 cannot easily be estimated by using t DSB and X X DSB . However, if we assume that c z 2 &c { 2 , we can estimate the stochastic rate constants as where the maximum number of repair proteins should be X max RP *100 because of the discrepancy between the theoretical and simulation results. Details are further discussed in the Materials and Methods section.

Biological interpretations for the emergence and consequences of bistability of ATM*
Here we discuss biological interpretations of our results. First, we consider the significance of the ATM autophosphorylation mechanism. The autophosphorylation of ATM increases the active ATM proteins, which again increases the active ATM proteins. Thus it works as a positive feedback loop of ATM*. This mechanism is described by a nonlinear Hill equation as we show in the Materials and Methods section. In this equation, a Hill coefficient n A is an important factor for nonlinearity. For example, as the Hill coefficient n A increases, the nonlinearity of the equation also increases. For our model, we assume that inactive ATM proteins exist as dimers and two ATM* proteins autophosphorylate the dimer which is indicated by Bakkenist et al. [16]. In this case, we can generally define the Hill coefficient as n A~2 , which means that two ATM* proteins bind to an ATM dimer protein and phosphorylate it. Then bistability emerges in this system for appropriate values of parameters. In Figure 5, we consider conditions for parameters j ATM and k DSB where bistability arises. Thus, our results indicate an physical importance of the ATM dimer proteins, in that the dimer leads to the high value of the Hill coefficient which emerges bistability of the steady state concentration of ATM*.
Biological importance of a positive feedback loop is considered by Xiong et al. [4]. They provide experimental evidence that in a physiological process of cell fate induction, Xenopus oocyte maturation, a bistable signalling system of a MAPK cascade converts a transient stimulus into a reliable, self-sustaining, effectively irreversible pattern of protein activities. Also, in DNA damage response processes, our model suggests that bistability, which is generated from autophosphorylation, amplifies transient damage signals (DSBs and DSBCs), and sustains the active state in the presence of DSB and DSBC noise. Other experiments indicate that bistability arises when there is a positive feedback loop in gene regulation networks [5,23,24]. In general, it is suggested that hysteresis, which is caused by bistability, leads to robustness against noise rather than an ultra-sensitive response in which the system is monostable [25]. Therefore, biological meanings of bistability which have two states (on and off states) will be (i) amplifying the input signal by the mechanism of hysteresis, and (ii) maintaining an on or off state in the presence of noise.
Other experimental results suggest that positive feedback of ATM contributes to genomic stability [26]. In the experiments, authors indicate that the initial DNA damage signal induces ATM activation and recruitment, and results in early H2AX phosphorylation immediately adjacent to DSBs. Phosphorylated H2AX then binds to MDC1 (mediator of DNA damage checkpoint protein 1) and causes additional activation of ATM. ATM then phosphorylates its substrates, resulting in checkpoint activation and DNA repair [26]. In the experiments, authors concluded that positive feedback of ATM amplifies the damage signal and it is vital in controlling proper DNA damage response and maintaining genomic stability [26]. Assumptions of our model are different from those of the above experiments in that positive feedback is caused by the autophosphorylation of ATM. However, our model also suggested that positive feedback of ATM amplifies DSB and DSBC damage signals and maintains a state in the presence of fluctuation of DSBs and DSBCs. Therefore, our model provides theoretical bases of the genomic stability which is indicated in the previous experimental results.

Effects of a negative feedback for the phosphorylated ATM concentration
Recent studies suggest that the p53/Mdm2 negative feedback loop generates oscillation of p53 and Mdm2 [27,28]. These pulses are initiated by DNA damage and the signaling kinase, ATM. Batchelor et al. suggest that the negative feedback between p53 and ATM, via Wip1, is essential for maintaining the uniform shape of p53 pulses [29]. We consider how this negative feedback affects ATM's dynamical behaviour. For simplicity, we only consider the case that the fixed points of ATM* have irreversible bistability in which the detection time is short. In the absence of negative feedback from p53, once ATM* is activated to the higher state, the concentration of ATM* is sustained with the same value because of the irreversibility. This phenomenon means that ATM works as a memory module of DNA damage. Figure 9A shows time courses of the concentration of ATM* without negative feedback. In this case, even when we add some transient stress (0#t#100 [min]), the concentration of ATM* sustains high values. Therefore, once the DNA damage is detected, the concentration of ATM* maintains high values without further DNA damage. In the presence of negative feedback from p53, the concentration of ATM* exhibits instability, and it oscillates under the constant stress. Figure 9B shows time courses of ATM* with a negative feedback loop. The concentration of ATM* is oscillating under some constant stress. However, it suddenly decreases after some transient stress (0#t#100 [min]). Detailed models are described in the Materials and Methods section. These results indicate that when DNA damage occurs, ATM detects DSBs as quickly as possible, and the concentration of ATM* becomes high. This high concentration state sustains until p53 is activated to suppress ATM phosphorylation. In addition, the concentration of p53 oscillates, which results from reactivation from ATM when DNA damage persistently exists.

Further challenges for the initial DNA damage response model
Experiments are needed to confirm whether ATM* has bistability which is caused by the autophosphorylation of ATM. As we show above, bistability collapses when the negative feedback from p53* exists. Therefore the phosphorylation of p53 should be blocked when we identify bistability. In addition, there are many signaling pathways from ATM [10,30]. The phosphorylation sites of p53 vary depending on the DNA damage agents, and it has not been known whether p53 is always phosphorylated and generates negative feedback to ATM for all stress signals. Quantitative data of the concentrations of p53 and ATM under several stress signals have not been observed enough. Further observations might provide new insight into DNA damage responses and signaling processes after damage. In addition, other experiments suggest that after exposure to H 2 O 2 , the p53 and ERK proteins are phosphorylated, which induce apoptosis or survival, respectively [31]. Also, the decision of the two exclusive fates is stochastically determined in independent cells. In this paper, our model only addresses the initial responses to DNA damage, but we may expand this model and clarify the mechanism of stochastic and exclusive decision of apoptosis.

Theoretical description of the mean numbers of DSBs and DSBCs
We can estimate the production rates of molecules per unit time.
Here we denote such rates as v RDSB~c3 X DSBC : ð14Þ The number of molecules for X i is time dependent (X t i ), but some of them are assumed to be in the mean convergence,X X i , meaning lim t??
If there are enough repair proteins, we can assume that the numbers of DSBs, DSBCs, and free RP are in the mean convergence, and we describe them asX X DSB ,X X DSBC , andX X RP . In addition, when we assume that the maximum number of repair proteins (X max RP ) is fixed and larger thanX X DSBC , we can estimatê X X RP~X max RP {X X DSBC . In this situation, the mean production rateŝ v v DSB andv v DSBC approach their steady states and become 0, and the mean ratev v RDSB approaches c 3X X DSBC . Here we estimate mean values of DSBs and DSBCs. In the steady states, X DSB and X DSBC satisfy When we add equation (16) to equation (17), we have the mean number of DSBCs by using the stochastic rate constants: Then we substitute equations (18) and (19) to (16), and we can describe the mean number of molecules for DSBs by using the stochastic rate constants and the maximum number of repair proteins X max RP :X Because the number of DSBs which are generated per unit time is c 1 , it takes at leastX X DSB C ð Þ c 1 to generate enough DSB(C)s. Thus the time until the number of DSBCs approaches the mean, X X DSBC , is calculated by and that for DSBs is calculated by Steady state analysis Here we show the calculation method of steady states and their stability [21]. The steady state concentration of ATM* is calculated by equation (11). Here we denote it as: where x denotes the concentration of ATM*. The steady states of x satisfy dx/dt = 0 and they are solutions of When we set a solution of f(x) = 0 as x~x x, the stable steady states satisfy We used Mathematica (version 5.2) to solve equation (24) and estimate the stability of solutions. The source code to find steady state concentrations of ATM* and the stability of them is included in Text S1.

The Gillespie algorithm
In this section, we introduce a method for calculating stochastic repair processes of DSBs. The Gillespie algorithm is one of the famous simulation methods [18,20]. It is clear that the time course of the state of the reaction system (the number of molecules of each type) can be regarded as a continuous Markov process with a discrete state space, because of the fact that the reaction hazards depend only on the current state of the system. Here we show a method for stochastic simulation of the time-evolution of the system. In a reaction system which contains m reactions, the hazard of reaction R i obeys h i X ,c i ð Þ, so the total hazard for all reactions is In this reaction system, each chemical reaction occurs following a Poisson process, and therefore time to the next reaction (dt) obeys an exponential distribution with parameter h 0 X ,c ð Þ, thus dt obeys Also, the probability that the ith reaction occurs after this time interval dt is proportional to h i X ,c i ð Þ, independent of the time to the next event. Therefore, the reaction type will be i with probability h i X ,c i ð Þ=h 0 X ,c ð Þ. Using the time to the next reaction and the reaction type, the state of the system can be updated, and simulation can continue. This simulation procedure was first proposed by Gillespie and is known as the ''Gillespie algorithm'' (or ''Gillespie's direct method'') [20]. The concrete procedure of this algorithm is as follows: Gillespie's direct method 4. Calculate the time to the next event dt which follows the exponential distribution. 5. Set t = t+dt. 6. Select a reaction type, i, based on probabilities h i X ,c i ð Þ=h 0 X ,c ð Þ for each reaction i. 7. Update X according to reaction i. 8. Output X and t. 9. If tvT max , return to step 2.
The source code of this algorithm is included in Text S1.

Hybrid simulation algorithm
The stochastic behavior of DSBs can be calculated by the Gillespie algorithm [19,20]. In equation (11), the numbers of DSBs and DSBCs, X DSB and X DSBC , are random variables. Here, we have to connect the continuous equations and the stochastic ones to solve the ATM sensor equation (equation (11)). In other words, the time until the next reaction occurs in the discrete (stochastic) regime, dt, is not constant, and it does not match the time step of the numerical algorithm in the continuous regime, dt. To settle the problem, we used a hybrid simulation method [18,22] in which some processes are simulated discretely while other processes are handled in a continuous manner by differential equations. Kiehl et al.'s method treats the effects of round trip conversions between discrete and continuous variables [22]. In our model, the numbers of DSBs and DSBCs affect the phosphorylation of ATM, but ATM proteins do not affect repair processes of DSBs, and thus the discrete reactions do not depend on the continuous reactions. That is to say, there are no time-varying hazard functions which are affected by ATM. Thus the hybrid simulation algorithm for our model can be simplified to the following procedure. In this algorithm, the time step of the numerical algorithm in the continuous regime is fixed as dt. At first, we determine initial conditions of both discrete and continuous values, and calculate reaction hazards of discrete reactions h i X ,c i ð Þ ð Þat time t. From these hazards, we can calculate the time to the next reaction (dt), and the reaction type i of the discrete regime by using Gillespie's direct method [20]. After these preparations, we compare the step sizes of discrete and continuous regimes, and select a smaller step size to update the simulation time as t: = t+min{dt,dt}. For both cases, we update the continuous variables to the values for the new time by using a numerical solution algorithm for ordinary differential equations (Euler, Runge-Kutta, etc.). Only when some discrete reaction has occurred, we update the discrete variable according to the reaction type i. The concrete procedure of this algorithm is as follows: 1. Set time t = 0, and determine initial conditions of X j t~0 ð Þ (j = DSB, DSBC, RP, RDSB) and [ATM*](t = 0).

Calculate the discrete reaction hazards
Þat time t. 3. From the hazards of the discrete reactions, select a discrete time step size dt and type of reaction i by using Gillespie's direct method. 4. If dt.dt (no discrete reaction has taken place), set t = t+dt, and update the continuous variables to the values appropriate for this new time. 5. If dt#dt (some discrete reaction has occurred), set t = t+dt, update the continuous values to those appropriate for this new time, and update the discrete variables according to the reaction type which is selected in step 3. 6. If t is less than the simulation duration, return to step 2.
The source code of this algorithm is included in Text S1.
Hill equation for the autophosphorylation mechanism of ATM units as ''molecules''. In previous experiments, it is suggested that ATM is autophosphorylated in 15 min after exposure to 0.5 Gy IR, which causes only 18 DNA double-strand breaks in a cell [16]. In addition, equation (20) indicates that if the maximum number of repair proteins is large, the mean number of DSBs becomes small. This result is also supported by our simulation results in Figure 4D. Therefore, we estimate that the main products in DNA damage processes are DSBCs which phosphorylate ATM. Based on this assumption, we estimate the theoretical time until DSBCs reach the steady state t DSBC v15 (min) because the time until ATM is autophosphorylated is within 15 min. Simulation results indicate that when t DSBC~2 (min), the mean number of DSBCs becomes the steady state in 15 min ( Figure 3B). In this case, the stochastic rate constant c 3 becomes c 3~1 =t DSBC~0 :5 (see equation (21)). Also we estimate the mean number of DSBCs asX X DSBC *20, and define the stochastic rate constants as c 1~c3X X DSBC~1 0 (see equation (19)). When the maximum number of repair proteins is large, bothX X DSB and t DSB become small, and the stochastic rate constants c z 2 and c { 2 have small effects onX X DSB and t DSB . Here we assume that the association rate of repair proteins (c z 2 ) is larger than the dissociation rate (c { 2 ), as the previous work suggested [28]. Thus we simply define these parameters as c z 2~0 :25 and c { 2~0 :025. As we showed in the previous section, the Hill coefficient n A denotes the number of molecules which bind to inactive ATM. Experimental results indicate that inactive ATM exists as a dimer [16], and thus we can assume two phosphorylated ATMs bind to an inactive ATM dimer. Therefore we estimate the Hill coefficient as n A~2 in our model. In addition, the rate constants k ATM1 and k ATM2 are estimated to be 0.1,10 in Ma et al.'s work [28]. Their model successfully explained the previous experimental results [28,32]. Thus we simply selected two values k ATM1~1 and k ATM2~1 which are between 0.1 and 10.
After the above preparations, we finally estimated parameters j ATM , k DSB , and ATM tot . To begin with, we fixed the total concentration of ATM as ATM tot ½ 2:0 (mM), and calculated the steady state concentrations of ATM* as functions of j ATM and k DSB ( Figure 5). Then we found the region of j ATM and k DSB where bistability of ATM* exists. We selected the rate constants as j ATM~0 :78 and k DSB~0 :0023 from bistable regions, because we predict that bistability plays an important role in DNA damage detection. Next, we studied dependencies of the total concentration of ATM on bistability ( Figure 6). This figure suggests that if bistability of ATM* exists, the total concentration of ATM needs to be larger than 1.89 mM.
In addition, when we substitute equations (44), (45), and (47) into equation (41), c { 2 can be approximated as In equation (47), we have to estimate a value of X max RP . However, the discrepancy between theoretical and simulation results makes it difficult to estimate the value. In particular, even if there are enough repair proteins, the mean number of DSBs from a simulation converges toX X DSB *0:5, which is not 0 as theoretical results predict. In equation (41), X max RP *100 when the mean number of DSBs satisfiesX X DSB *0:5.

Supporting Information
Text S1 This file includes 4 source codes of C program and 1 source code of Mathematica to calculate our models, which use Gillespie algorithm or hybrid simulation algorithm (for C codes), and steady state analysis (for a Mathematica code). Details are shown in the README text. Found at: doi:10.1371/journal.pone.0005131.s001 (0.01 MB ZIP)