The evolution of p53 network behavior

We study the evolution of the p53 core regulation network across the taxonomic span of humans to protozoans and nematodes. We introduce a new model for the core regulation network in mammalian cells, and conduct a formal analysis of the different network configurations that emerge in the evolutionary path to complexity. Solving the high dimensional equations associated with this model is typically challenging, and we develop a novel algorithm to overcome this problem. A key technical tool used is the representation of the distinct pathways in the core regulation networks as “modules”, such that the behavior of the composite of two or more modules can be inferred from the characteristics of each of the individual modules. Apart from simplifying the complexity of the algorithm, this modular representation also allows us to qualitatively compare the distinct types of switching behaviors each network can exhibit. This then allows us to demonstrate how our model for the core regulation network in mammalian cells matches experimentally observed phenomena, and contrast this with the plausible behaviors admitted by the network configurations in putative primordial organisms. We show that the complexity of the p53 core regulation network in vertebrates permits a range of behaviors that can bring about distinct cell fate decisions not possible in the putative primordial organisms. Significance Statement The p53 protein has been protecting organisms from tumors for a billion years. We study the link between the evolution of the p53 network structure and its corresponding tumor suppression strategies. We compare the dynamical behaviors in putative primordial organisms with simple networks with the vertebrate network that contains multiple feedback loops. We show that the vertebrate network, but not the ancestral network, can both repair moderate damage and induce apoptosis if too much damage accumulates, balancing the risk of cancer with the cost of too much cell death. Moreover, the complexity of the vertebrate network allows for adaptation, for example to increase p53 network sensitivity, which is consistent with recent research on large mammals.

tations that occur when the cellular DNA is damaged and not repaired to its original state.When the mutations interfere with the normal regulation of cell division, then cell proliferation could cause the tumor to become cancerous [1].While there is abundant evidence that exposure to chemical carcinogens and radiation has had an impact on cancer rates [2], there is also evidence that cancer is a general feature of multi-cellular organisms that has presented a long-standing evolutionary challenge [3].
A key component of the DNA damage response mechanism is the T P 53 gene which expresses the p53 tumor suppressor protein, the so-called "guardian of the genome" [4].p53 and its ancestors have been protecting metazoans from mutations arising from DNA damage for over one billion years [5].Recent research has shown that p53 plays this role by first sensing DNA damage, and then mediating various downstream processes in response.p53 can promote cell survival by upregulating genes that bring about cell-cycle arrest or DNA repair.However, p53 can also promote cell inactivation or death by bringing about permanent cell-cycle arrest or cell death [6][7][8][9].How p53 levels and dynamics determine different cell fates in response to DNA damage is a crucial component of an organism's tumor suppression mechanism [10][11][12][13].It is therefore no surprise that a mutation in p53 is implicated in over 50% of human cancers [14].
One of the oldest known mechanisms cells employ in response to DNA damage is programmed cell death, also known as apoptosis [15].Cells with DNA that is damaged above some "apoptotic threshold" level generally display sustained high levels of p53 [13,16,17], inducing the expression of proelimination genes [12,[16][17][18]] that bring about apoptosis, thus limiting the risk that a cell lineage will become cancerous.The underlying behavior can be understood in the context of Statistical Decision Theory by imagining that cells attempt to decide on the truthfulness of the statistical hypothesis: "The damaged DNA can be repaired to its original form".
It is clear that setting a low apoptotic threshold minimizes the risk of cancer since a cell that experiences potentially cancer-causing DNA damage would go to apoptosis.A feature of such networks is that they abet Type 1 errors.Such errors happen when the statistical null hypothesis is true (i.e., DNA can be repaired to its original form), but the cell is killed anyway.While not ideal, Type 1 errors would still be preferable to Type 2 errors, where a cell carrying a cancer-causing mutation could be allowed to live.A single Type 2 error could lead to lethal cancer, and this is in fact known to be one of the "hallmarks" of cancerous clones [2,19].That said, overly responsive networks with a low apoptotic threshold may result in developmental defects, reduced tissue growth, or high metabolic costs to replace the dead cells.A major task for p53 and the DNA damage response mechanism is therefore to effectively suppress the onset of tumors by selecting different cell fates that would increase whole organism fitness.This is achieved by trading-off between minimizing the number of Type 1 errors that take place, while also trying to suppress all Type 2 errors.
While p53 is a vital component of the damage response network, it acts as part of a network of interacting genes that

Significance Statement
The p53 protein has been protecting organisms from tumors for a billion years.We study the link between the evolution of the p53 network structure and its corresponding tumor suppression strategies.We compare the dynamical behaviors in putative primordial organisms with simple networks with the vertebrate network that contains multiple feedback loops.We show that the vertebrate network, but not the ancestral network, can both repair moderate damage and induce apoptosis if too much damage accumulates, balancing the risk of cancer with the cost of too much cell death.Moreover, the complexity of the vertebrate network allows for adaptation, for example to increase p53 network sensitivity, which is consistent with recent research on large mammals.

D R A F T
determine cell fate.Extensive genetic studies have uncovered a complex network of hundreds of genes [20] that cooperate to carry out tumor suppression as soon as damage is registered in the genome [21].There are a multitude of upstream sensors that detect DNA damage, the transducers that relay the information to p53 and the downstream actuators that are regulated by p53 to bring about different cell fates [22].In addition, the dynamics of p53 is known to be governed by its interaction with a set of core regulation network proteins including M DM 2, P T EN and ARF through multiple feedback pathways [23].Mutations in some of the respective genes have even been implicated in human cancers, although not to the extent of p53 [24].
It is therefore interesting that while homologs of T P 53 are inferred to have been present in the earliest metazoan ancestors [5], the same cannot necessarily be said about M DM 2, which is another core regulation protein in the mammalian p53 pathway.Based on these results, we introduce the underlying network configurations that are present in extant vertebrates and are inferred to have been present in the ancestors of these vertebrates.We also introduce a novel dynamical model that mathematically describes these networks.Our focus is on whether these alternative network structures would allow for qualitatively distinct forms of p53 regulation, which in turn would lead to qualitatively different strategies for dealing with the risk of cancer.Clearly, an individual's fitness will be low if it is at a high risk of developing fatal cancer.However, a policy whereby p53 and its core regulation proteins bring about apoptosis of too many cells at the slightest sign of DNA damage would also incur high metabolic costs.
A core question in the evolution of tumor suppression mechanisms therefore, is how the relevant genetic networks can optimally control the response in order to balance the risk of cancer with the risk of false responses to signals of DNA damage.To answer this question, we modularly study how the topology of alternative p53 core regulation networks integrate DNA damage information and determine the range of dynamic responses to DNA damage.Using a novel algorithm to solve for the equilibrium points of the different network models, we explore the qualitatively different types of switching behavior that the various network configurations can admit.This provides an insight into how the different core regulation genes play a functional role in determining the p53 response to DNA damage, and further reveals a novel mechanism by which ARF could play the role of an apoptosis enhancer [25].
We validate our model for the p53 network in vertebrates by using it to explain many recent experimental results in a way that previous models are unable to do.Specifically, our model predicts that for a moderate amount of damage, the p53 level will pulse, with the number of pulses proportional to the amount of damage.When the damage becomes too high, the p53 level elevates monotonically to a high level and remains there.This matches experiments on p53 in human breast cancer cells [10,26,27], and also more recent experiments on human U2-OS cells [13].An important prediction made by our model in this regard is that the monotonic switch of p53 to a high value in response to a large amount of DNA damage is a result of a property of the core regulation network model and not due to temporal effects, ensuring that the mechanism is robust.Our model also predicts that the upstream trans-duction kinases are responsible for activating p53 pulses when damage is present and that pulses in kinase levels are highly coupled with p53, matching recent experimental observations [10,28].
This work also elucidates how the p53 responses exhibited by different network structures relate to the potential costs and fitness benefits of different organisms.We show that the range of behaviors permitted by the vertebrate network is crucial in minimizing the frequency of Type 1 errors while eliminating Type 2 errors in human cells.We also show that the complexity of the network found in vertebrates allows for other distinct strategies for tumor suppression, as has been observed in elephants [29].The alternative network configurations in primordial organisms however, cannot express the range of behaviors that are observed in and permitted by vertebrate cells.

Taxonomic representation of regulatory genes and network configurations
Among clinical cancers, mutations in p53 are among the most commonly detected [4], and p53 is widely recognized as playing a key role in tumor suppression in multicellular organisms [7,23].The p53 family of genes is at least one billion years old [5], appearing in organisms as diverse as protists, cnidarians, and mammals.The tumor supression role of the p53 ancestral gene appears to have been preserved over a large time span; its function in cnidarians is to protect the germline gametes from DNA damage, and this function persists in arthropods, annelids, molluscs, and vertebrates [5].
In humans, p53 is part of a sophisticated network of hundreds of genes [20] that cooperate to mediate cell fate decisions such as the initiation of cell-cycle arrest, DNA repair, senescence or apoptosis, as soon as DNA damage is detected [21].This network includes a core regulation network [23] consisting of the proteins M DM 2 [30][31][32], P T EN , [33,34], and ARF [35][36][37] among a host of other upstream, downstream and intermediate species involved in sensing, transduction and regulation [21][22][23].
Recent studies have uncovered multiple aspects of how these core regulation proteins interact.It is well known that p53 transcriptionally activates the M DM 2 gene.M DM 2 in turn antagonizes p53 forming a negative feedback loop around p53 [30,31,38].p53 also activates the P T EN gene, and P T EN proceeds to down-regulate M DM 2 through a series of interactions [34,39,40], which forms a positive feedback loop around p53. ARF is known to cause the translocation and eventual degradation of M DM 2 [35,36,41], and it has been recently discovered that ARF can also mediate the degradation of M DM 2 [37], leading to a positive feedback loop around M DM 2. The network configuration of these four genes will henceforth be referred to as the full network Nv, and is illustrated in Fig. 1.Much of our understanding of the core regulation network interaction and dynamics in response to DNA damage comes from studies of human and mouse cell lines [42], and this network is most likely common to at least mammals.As we will see in this paper, this network configuration is crucial in explaining experimental observations about p53 behavior in human cells in response to DNA damage.

D R A F T
of vertebrates.Recent work suggests that M DM 2 and its paralog M DM 4 arose from a gene duplication event and became the primary negative regulators of the p53 family of genes [43].While invertebrate lineages contain a M DM 2 homolog, it bears little functional resemblance to the vertebrate M DM 2 and there is no direct evidence that it plays a role in regulating p53 dynamics.Examining the Network Nv in Fig. 1 reveals that the absence of this M DM 2 would prevent P T EN and ARF from affecting the dynamics of p53, because both P T EN and ARF only regulate p53 through feedback connections involving M DM 2. This means that in the absence of M DM 2, the response of p53 to DNA damage would only depend on the direct regulation of p53 by upstream inputs.This suggests that the dynamics of the network may depend critically on the feedback loops that M DM 2 mediates and motivated us to explore the possible intermediate network configurations that might have emerged in the evolution from the putative primordial gene network to that observed in modern vertebrates.
Evolutionary history of the p53 network.Based on the taxonomic representation of p53 and M DM 2, we hypothesize that the p53 function as a tumor suppressor evolved prior to the feedback loops mediated by M DM 2. This hypothesized ancestral network is labeled Ni and is shown in Fig. 1.Following from the work of Momand et al [43], transitions to the network observed in modern mammals must have proceeded through the intermediate networks shown in Fig. 1.Two possible evolutionary paths are possible from the ancestral network Ni to the mammalian network Nv.In both evolutionary paths, the p53−M DM 2 negative feedback interaction evolves first since this interaction is central to the regulation of p53 by either P T EN or ARF , leading to the network N1 shown in Fig. 1.In the first path, the interactions through which p53 inhibits M DM 2 through P T EN evolves first (leading to the network N2 as seen in Fig. 1), followed by the mutually inhibitive feedback loop between M DM 2 and ARF , leading to the network Nv.In the other path, the mutually inhibitive feedback loop between M DM 2 and ARF evolves first, leading to the network N3.This is followed by the interactions through which p53 inhibits M DM 2 through P T EN , leading again to the network Nv.
Our primary emphasis is on providing insight into the behaviors permitted by these different network configurations from a dynamical systems standpoint.As such, we do not make claims that the paths to complexity studied in this paper are the only possible routes that evolution might have evolved.Rather, we focus on exploring the different qualitative behaviors permitted by the intermediate network configurations and distinguish these from those observed in the mammalian network.

Modular models of the p53 pathways
The diagrams shown in Fig. 1 are conceptual in that, while each block is labeled with the name of a single protein, it typically represents several interacting chemical species.For example, the block labeled M DM 2 includes the dynamics of the M DM 2 protein in its active an inactive forms, as well as its corresponding mRNA.We should thus think of each block as a "module" rather than an individual protein.
To predict how the network behavior changes as modules are added, the different modules in Fig. 1 must be associated with precise mathematical models, so as to exhibit two key properties: dynamic modularity and parametric modularity.Dynamic modularity ensures that the properties of each module do not change upon interconnection with other modules, and is therefore an essential property to ensure that the independent analyses of different modules can be combined to make precise statements about the behavior of the entire network, and specifically, how the behavior changes as modules are removed or inserted into the network.This can enable the analyses of large networks and also the design of novel networks [44][45][46].Parametric modularity, which has not been explicitly mentioned as often in the literature, implies that the network parameters that are associated with a given module do not appear in any other module.This property ensures that each module can be analyzed independently to understand how various parameters or inputs to that module affect the outputs from the module.Fig. 2 shows a mathematically precise "block-diagram" representation of the full network configuration from Fig. 1.Each module, labeled from M1 to M4, is characterized by a system of ordinary differential equations (ODEs) depicted in Fig. S1 that describe the behavior of the species in that module.The solid arrows between the modules represent the signals that flow between them.For example, the arrow originating at M2 and ending at M1 represents the mechanism by which active M DM 2 represses p53.This signal is said to be an "output" of M2 and an "input" to M1.A more comprehensive block diagram representation of this network and the corresponding block diagram representations of the other network configurations are shown in Figs.S2 and S3 of the Supporting Information, respectively.
The species associated with each module are as follows.Module M1 is associated only with the dynamics of the p53 protein.Module M2 consists of the dynamics of M DM 2 in both its active (M DM 2a) and inactive (M DM 2) forms, along with its mRNA (mRN Am).Module M3 consists of the dynamics of P T EN and its mRNA (mRN Ap), along with the intermediate species P IP 2 and its phosphorylated form P IP 3, and AKT and its active phosphorylated form AKT a. Finally, Module M4 contains the dynamics of the ARF protein.
The inputs to the full network are D1, which represents kinases that activate p53 and cause the degradation of M DM 2 in response to DNA damage [40], and D2, which represents transducers like P ARP 1 and E2F 1 which activate ARF in response to DNA damage [25,47].A more detailed description of the model and its parameters is provided in the Supporting Information.

The role of each p53 network configuration in determining switching behavior
In this section, we discuss the different dynamic behaviors that are permitted by each p53 network configuration we identified on the pathways to complexity.A specific type of behavior that we study is the ability of the p53 concentration to behave like a "switch" in response to DNA damage.
Multistability and biological switches.Biological networks commonly use switch-like behavior to turn a continuous input signal (usually an extrinsic stimuli) into an "all-or-nothing" output signal response (usually species or compound concen-   trations).A common mechanism employed by these networks to bring about this behavior is bistability.A bistable network has two sets of stable steady-states, in which the state of the network transitions from a neighborhood of one set of steadystates (which we call A) to the other (which we call B) when an input signal grows above a certain threshold (which we call T1).A unique property of these networks is that once the transition occurs, even if the input later falls below T1, the state may still remain in a neighborhood of B. Reversion to a neighborhood of A is generally observed only when the input signal falls below another threshold T0, where T0 < T1.This property is known as hysteresis, and refers to the ability of a bistable system to "remember" that the input stimulus was above T1 long after that stimulus is removed, until it falls below T0 [48].The lac operon network in E. Coli, the Wee1-Cdc2 network, the N F κβ response in an apoptosis network and the cell cycle oscillator in Xenopus laevis are among the many naturally occurring biological networks whose responses have been modeled as bistable switches [48][49][50].In some cases, the lower threshold T0 does not exist, which makes the transition from A to B irreversible.It is also plausible to have more than two sets of stable steady-states, which is known as multi-stability.
Monostable networks on the other hand have states that slide along a continuum of steady-states [49] in response to input signaling.Such networks can only exhibit switch-like behavior if the system is ultrasenstive; that is, if the inputoutput relationship of the network is described by a sharp sigmoidal curve [48].When the input signal crosses the exponential phase of the sigmoid, the output response of the network transitions from a low to a high state.In fact, the dynamics of transcription factor target genes as a function of transcription factor concentration is typically modeled as a sigmoidal "Hill function", which can be viewed as an ultrasensitive switch.However, such switches are memoryless: once the input signal is removed or is reduced, the system immediately returns to its original state [48].
The different p53 network configurations can display various switching mechanisms with respect to the two damage sensing inputs [D1] and [D2].Our modular study revealed that as [D1] is varied, the network configurations Ni and N1 operate as monostable switches, in the sense that when an increase in [D1] makes the p53 level rise, a subsequent decrease in [D1] to its original value causes the concentration of p53 to revert to its initial level.The steady-state concentration of p53 as a function of [D1] is shown in Fig. 3.The main difference between the two networks Ni and N1 is that negative feedback of p53 with M DM 2 typically makes the p53 steady-state in N1 to be lower than that of Ni.
Networks with a positive feedback loop admit switching with "memory", allowing for sustained high levels of p53 expression once damage crosses a threshold.With model parameters similar to those from [40], the network N2 operates as a bistable switch with respect to the input [D1].From this study, we can plot the bifurcation diagram of N3 with respect to the parameter [D1] in Fig. 4(a).The solid blue line signifies a continuum of stable steady-state concentrations of p53 as a function of [D1], while the red dashed lines represent unstable equilibrium points.It can be seen that when [D1] is between the two thresholds T0 and T1, there are two stable  To illustrate the behavior of the bistable network, we run a simulation as shown in the red curve of Fig. 4(c).We observe the manifestation of the "memory" property of bistable networks, where after rising to a point in set of equilibrium points labeled as "B", the p53 level returns to a point in the set "A" after [D1] crosses the threshold T0 from above.Fig. 4(c) further contrasts this to the behavior of the monostable networks Ni (blue curve), where the switch has no memory.
A change in the parameters corresponding to the rate of initiation of p53 and ubiquitination of M DM 2 by D1 causes the threshold T0 from Fig. 4(a) to become less than 0, as shown in Fig. 5(a).In this case, the network exhibits an irreversible switch from the set of points in Region "A" to those in Region "B".This is illustrated in Figs.5(b)-(c).The bifurcation diagram and simulations are normalized by the p53 level at [D1] = 0.2 for clarity.
The network N3 can also behave as a bistable switch with respect to [D1] provided that [D2] also increases along with [D1].This is a reasonable assumption, since both [D1] and [D2] will increase with DNA damage.Qualitatively, this network admits the same behavior as N2.The bifurcation diagram for this network is provided in Fig. S4 of the Supporting Information.
Unlike in the monostable networks, bistability allows for the p53 levels to remain high for a sustained period of time once the damage sensor level has crossed threshold T1, even if the damage sensor level subsequently falls below this level.This can induce p53 to up-regulate downstream pathways that bring about cell-cycle arrest and DNA repair, as well as those that bring about apoptosis or senescence.Since apoptotic pathways are typically slower than repair pathways [51], it is more likely that apoptosis will be the cell fate either if the damage sensor level remains high for sufficiently long, or if the network operates as an irreversible bistable switch, since a sustained high levels of p53 expression is known to bring about apoptosis [12,16,18].
The full network Nv can admit both a reversible and an irreversible switch, allowing for distinct p53 responses to different levels of DNA damage.The addition of the ARF module M4 to N2, or the P T EN module M3 to N3, can bring about tri-stability in the p53 levels with respect to [D1].That is, there are three sets of stable steady-states between which the p53 level can toggle, which we label "A", "B" and "C".The bifurcation diagram illustrating the steady-states of p53 in this network Nv is shown in Fig. 6(a).From this diagram, we observe that the original switching behavior which was a property of the bistable networks is retained; when [D1] crosses T1 the p53 level switches to a high value, and this switch is "turned-off" only when [D1] falls below T0.In addition to this behavior, the tri-stable network admits a third threshold T2, that can be seen in Fig. 6(a).Once [D1] crosses T2, p53 transitions irreversibly to a high value; no matter how [D1] changes following this transition, p53 will not return to its nominal level.This behavior is illustrated in Figs.6(b)-(c).The bifurcation diagram and simulations are normalized by the p53 level at [D1] = T0 for clarity.
While the bistable networks can operate as either reversible or irreversible switches, the full network is able to simultaneously admit both switching behaviors.From the perspective of our network, the level of damage-sensor [D1] crossing T1 is likely associated with a scenario that enough damage has been triggered to initiate a p53 repair response, but not to bring about a terminal cell fate.On the other hand, [D1] crossing  T2 likely means that too much damage has been sustained, causing an irreversible switch in the p53 level to a high value and eventually bringing about either senescence or apoptosis.It turns out that the level of D2 also plays an important role in shaping the p53 response.We recall that D2 represents transducers like P ARP 1 and E2F 1 which activate ARF in response to DNA damage.As [D2] increases, the rate of production of ARF is increases.This in turn causes the middle branch of stable equilibrium points (labeled 'B') to "shrink", as illustrated in Fig. 6(a).A lower threshold T2 would mean that the cell has a higher chance of going to senescence or apoptosis, since the irreversible switch in p53 levels is triggered for lower levels of [D1].This behavior illustrates how ARF activation can trigger apoptosis, as has been observed experimentally [25].The dynamical behaviors permitted by the different networks are summarized in Fig. 1, which also shows the possible evolutionary paths from putative primordial organisms to extant vertebrates.

Simple models of DNA damage transduction and repair explain experimentally observed results
We now show how our tri-stable model for the p53 core regulation network in vertebrates Nv, coupled with simple models of DNA damage transduction and repair, explains many experimentally observed results.We model the activity of the transduction kinases such as AT M and CHK1 using a single variable D1, similar to what was done in [10].We model the dynamics of D1 as In this equation, the variable #DSB represents the number of DNA double-strand breaks (DSBs), which activate D1 through the function g1(#DSB).D1 then goes on to activate p53 and cause the accelerated degradation of M DM 2 [32,40], per the model in Fig. 2. The feedback term g2([p53]) represents negative feedback between p53 and the transduction kinases [10,28], and the rate constant γ represents the p53-independent degradation of the kinases.We assume that the dynamics of Eq. ( 1) is slow compared to the dynamics of the p53 core regulation network, which is a key element in explaining the experimentally observed behavior.

D R A F T
We model the dynamics of #DSB as a simple deterministic death process, where µ#DSB [p53] is the rate of repair, that is known to be regulated by p53 [52][53][54].Since damage is induced in the cell on the order of minutes [55], while the response can last for many hours or even days [26,27], the number of DSBs induced by damaging agents is introduced as the initial condition of Eq. ( 2).D2 represents the activity of proteins such as P ARP 1 [47] and E2F 1 [25], which are known to play a role in activating ARF in response to DNA damage.We model [D2] as a simple monotonically increasing function of the number of DNA DSBs [D2] = g3(#DSB), [3] with the rate of production of ARF governed by [D2].The specific parameters and functions from Eqs. ( 1)-( 3) are given in the Supplemental Information.
We couple these models with the full network Nv and report on the resulting simulations.In Fig. 7(a), we show that the p53 levels remain low for a damage resulting in 5 DNA DSBs, which is a low level of DNA damage.In this case, the DNA damage is repaired before a p53 pulse is initiated.For damage resulting in 10 DSBs, a single p53 pulse is observed as shown in Fig. 7(b).This matches experimental results from [26], which show that 0.3Gγ of damage is sufficient to initiate a pulse, where 1Gγ of damage results in between 30-35 DSBs [27,56,57].In Fig. 7(c), our simulation shows that DNA damage that brings about 150 DSBs causes the initiation of multiple p53 pulses for over 60 hours until the damage is repaired, matching experimental results that suggest that the number of pulses increases as the number of DSBs increases [4,26,27], and that these pulses could last for days after damage is induced [27].All pulses in our simulation have a period of between 4-7 hours, matching these experimental results in human breast cancer cells.
Early modeling work suggested that these pulses were primarily a result of the negative feedback interaction between p53 and M DM 2, along with some positive feedback or delays [55,[58][59][60][61].However, more recent experimental results, in which pulses in p53 are observed to be highly coupled with pulses in the upstream transduction kinases, suggest that p53 pulses are brought about by the interaction of p53 with these kinases [10,28].Another important observation from [10] is that transient activation of the kinases results in a pulse of p53, suggesting that there is an excitable mechanism controlling p53 pulses [10,28].It was noted that while a single negative feedback loop between p53 and these kinases (with delays) would be sufficient to generate sustained oscillations, other p53 interactions would be necessary to explain the observed excitability.
In our model, the interaction between D1 and the p53 network Nv can be described as an excitable mechanism known as a relaxation oscillator [62].The slow negative feedback between D1 and p53 described in Eq. ( 1) operates over the faster dynamics of the p53 core regulation network.The pulses in p53 are observed because of the resulting cyclic transfers between the p53 equilibrium points in Regions 'A' and 'B' [Fig.6(a)].The pulses in D1 levels are also seen to be highly coupled with the pulses in p53 levels.In this way, our simula-tion results strongly matches both experimental observations, and hypotheses about the underlying mechanism that brings about the observed behavior without recourse to additional mechanisms or interactions.
In Fig. 7(d), our simulation shows the p53 levels monotonically elevating to a high level for a large number of DNA DSBs, and remain high without pulsing.Even with the progression of the DNA repair process, the p53 levels remain high, due to the irreversible switch permitted by the tri-stable model.This behavior also matches recent experimental observations in human U2-OS cells, in which the p53 pulsing is observed for moderate levels of DNA damage, but very high levels of DNA damage trigger a strong monotonic elevation of the p53 level that does not return to a low level after rising [13].
It is worth noting that the experimental results are from different cells.Our goal is to show that the experimentally observed behavior is possible, and not to identify specific parameter values for the different cell types.

Discussion
Some of the earlier p53 modeling work was carried out under the assumption that the transduction kinase levels were proportional to the amount of damage, and that pulses were primarily observed due to the negative feedback between p53 and M DM 2 [58,59].Other modeling work, while considering the positive feedback loop, assumed that the upstream kinase concentrations were proportional to damage and caused autonomous pulsing in p53 when the kinase levels were sufficiently high [55,60,61,63].In contrast, our model takes into account recent experimental results which suggest that upstream transduction kinases are responsible for initiating p53 pulses, and that the pulses of these kinases are highly coupled with pulses of p53 [10,28].Moreover, these previous models primarily explain plausible mechanisms for the observation of oscillations in p53, and do not postulate the irreversible switch in the concentration of p53.
There are other models of the p53 core regulation network in humans that display both the pulsing and sustained high levels of p53 [40,63,64].However, our model matches recent experimental observations that demonstrate pulsing in p53 for a moderate amount of damage, and a strong monotonic elevation in p53 for high levels of damage [13], contrasting to these previous models in which p53 pulses before switching to a high value [40,63,64].
To our knowledge, our model of p53 core regulation in vertebrates is the first where the irreversible switch of p53 to a high level is an inherent property of the tri-stable network structure, as opposed to a consequence of temporal behavior.It is also worth noting that all the observed behavior can be explained precisely using the dynamical properties of the network ODE model, and do not necessitate the artificial introduction of external signals, switches or clamps.
A novel behavior that our model predicts is the role played by ARF in bringing about an irreversible cell fate.It has long been postulated that ARF enhances apoptosis by sequestering away M DM 2, thereby causing an increase in p53 levels.Although models to capture this behavior have been proposed [35,65], the details of the underlying dynamics are still largely unknown [66].A recent study has further pointed to a novel mechanism through which M DM 2 can inhibit ARF [37].In our model, this mutually inhibitory feedback between active  Our model predicts that an increased rate of production of ARF is an important factor in bringing about the irreversible switch of p53 to a high level.In this way, our model proposes a novel and mathematically sound reasoning for how ARF might behave as an apoptosis enhancer.
From an evolutionary standpoint, we can understand the selection pressures on the tumor suppressor network as balancing a trade-off between minimizing the number of Type 1 errors, while eradicating Type 2 errors.Type 1 errors occur when the statistical null hypothesis that DNA can be repaired to its original state is true, but the cells are killed anyway, while Type 2 errors occur when the null hypothesis is false, but the cells are not killed.While Type 1 errors could have metabolic costs in that too many cells are killed, Type 2 errors are potentially fatal.
From this perspective, the functional role of the p53 core regulation network in humans is clear.For low amounts of damage, the p53 level remains low.This is associated with the fact that the damage is deemed too low for a tumor to develop, and hence a significant p53 response is not initiated.
For moderate amounts of damage, p53 exhibits pulsatile behavior, with the number of pulses proportional to the amount of damage.At this point, the damage triggers the repair mechanisms by inducing pulses in p53 until all the damage has been repaired.This is consistent with experimental results which show that majority of cells that exhibit pulsed p53 are able to grow and divide after recovery from DNA damage [12].For a sufficiently large amount of damage, the p53 level switches monotonically to a high level and remains high.In cells that exhibit sustained high levels of p53 signaling, most cells go into senescence or apoptosis [12,13,18].This switch to fixed high p53 levels removes the possibility that a reduction in the damage signal could allow un-repaired cells to proliferate.In this way, the cell is able to minimize the number of Type 1 errors by repairing DNA damage for a range of moderate levels of damage, but also has a clear threshold to bring about an irreversible switch to apoptosis when the amount of damage sustained is too high.
If every cell had an equal probability of developing cancer per unit time across all organisms, then large long-lived mammals would have a very high risk of developing cancer as

D R A F T
compared to smaller short-lived organisms.However, there is a lack of correlation between the body size of an organism and its cancer risk [67], and this is known as "Peto's Paradox".There are many hypotheses to resolve Peto's Paradox, as outlined in [67].Two of these hypotheses, having a redundancy of tumor suppressor genes and a more sensitive apoptotic process, were found to be true in a recent study on elephants [29].It was found that elephant cells have 20 copies of the T P 53 gene as compared to the single copy observed in human cells.Elephant cells were also found to carry out p53-dependent apoptosis at a much higher rate than human cells for the same amount of DNA damage [29].
To test the effect of increasing the number of copies of the T P 53 gene, we took the full vertebrate network Nv described above and increased the basal and kinase-dependent rates of production of p53, while keeping all other parameters constant.We found that this would cause the level of threshold T2, shown in Fig. 6(a), to decrease.This change results in p53 making an irreversible switch to a high level and hence bring about apoptosis for a smaller number of DSBs, which would explain the experimental results observed in elephants [29].A key point to note is that the increased apototic rates result from the temporal dynamics of the network, rather than from a simple increase in the basal concentration of p53.Our model also predicts that increasing the ARF sensitivity to DNA damage would also have a similar effect, and it would be interesting to study this gene in other large long-lived mammals.
The networks N2 and N3 were found to admit both reversible and irreversible bistable switching behavior.It is worth noting that these networks are not capable of exhibiting tristability regardless of the values of the model parameters.
For the reversible switch shown in Fig. 4(a), the p53 dynamics of these networks would be qualitatively similar to that of Nv in the presence of low and moderate amounts of DNA damage when coupled with the transducer and repair dynamics we introduced.When the DNA damage is very high, the p53 level monotonically elevates to a high level.However, this switch is reversible if the DNA damage is repaired before the cell goes into senescence or apoptosis.This is illustrated for N2 in Fig. S5 of the Supporting Information.However, since the cell would have experienced a large amount of damage, the DNA repair may not have been successful or may have introduced genome destabilizing mutations.Hence, while the reversible bistable switches can minimize Type 1 errors by repairing moderate levels of damage, the elimination of Type 2 errors in the bistable networks requires the fast initiation of apoptotic pathways once the p53 level becomes high.
In the case of the irreversible switch shown in Fig. 5(a), p53 switches irreversibly to a high level once [D1] crosses threshold T1.However, in this case p53 is not able to admit pulsatile behavior, increasing the chance of the cell going to senescence or apoptosis.Hence, this policy would be very effective in eliminating Type 2 errors, but could also increase the frequency of Type 1 errors.In essence, the bistable networks N2 and N3 are either able to admit p53 pulses or an irreversible switch of p53 to a high level triggering apoptosis, but not both phases of behavior as the tristable network is able to.These networks would probably have been more likely to behave as irreversible switches, since the elimination of Type 2 errors is crucial in the cell's ability to avoid tumors.The network configurations Ni and N1 will only admit a single continuum of stable equilibrium points [68], and hence can only exhibit ultrasensitive switching behavior.In the absence of more complex mechanisms, the p53 levels are incapable of pulsing.Moreover, when the p53 level is high, it is more sensitive to changes in D1 levels than the bistable switch.In these organisms, it is therefore likely that the exponential phase of switching happens for reasonably low amounts of damage, and that the pro-apoptotic pathways were very sensitive to rising levels of p53, to ensure that Type 2 errors do not occur.This hypothesis is consistent with an earlier proposal on the operation of the p53 family of genes in early organisms [5].
In conclusion, our model shows that one adaptive value of the p53 core regulation network in vertebrates is its ability to balance the risk of cancer with the cost of too much cell death, and this behavior is clearly observed in experiments on human cells.The complexity of the p53 network in vertebrates also permits other organisms to operate at different points on this trade-off curve.On the other hand, inferred ancestral organisms with alternative p53 network configurations admit a narrower range of behaviors, which in turn limits the range of the possible points which they can lie on this trade-off curve.

Materials and Methods
The computation of the bifurcation diagrams in Figs.4(a), 5(a) and 6(a) requires the computation of equilibrium points for the system of ODEs in the corresponding networks.In practice, this reduces to the problem of computing the solutions to multivariate polynomial equations, which is known to be an NP-complete problem [69].While simple equations can often be solved using commercially available numerical solvers, more complicated systems such as the ones that arise in our analysis, exhibit multiple equations with many possible solutions that are close to each other.In this case, solvers or numerical continuation software for the computation of bifurcation diagrams often miss solutions that might exist, which was the case when attempting to evaluate the equilibrium points of the tri-stable and bistable networks.
To overcome this problem, we developed an algorithm to solve for the equilibrium points of a network as a function of the inputs to the network.This algorithm was used to solve for the p53 equilibrium points as a function of the concentrations of the transducer inputs D 1 and D 2 , to produce the plots in Figs.4(a), 5(a) and 6(a).
A key aspect in the operation of the algorithm is the modular decomposition of the p53 network into the modules shown in Fig. 2.
The procedu]re that the algorithm employs involves gridding the equilibrium values of the input and output signals for each module to form a table known as the Steady-State Behavior (SSB) of a module [70].Typically, this table is obtained by gridding the equilibrium values of the outputs for a range of constant inputs, and this special case of the SSB has been termed the Static Characteristic Function (SCF) [71,72] of the module.It is worth noting that these tables could be multi-valued, as there could be multiple output equilibria for a given input as is evidenced in the bistable and tri-stable networks.
Once the SSB of each module has been determined, the algorithm proceeds to systematically eliminate the latent variables in each of these tables, which are typically the input and output signals from each module that are interconnected to other modules in the network.This procedure is outlined in detail in the Supporting Information, and is crucial to obtaining and understanding the bifurcation behavior that each network configuration can admit.u 23 and u 24 .u 23 " rD 1 s represents the concentration of upstream transduction kinases like AT M , which mediate the down-regulation of active M DM 2 [5][6][7].u 24 represents the ARF -mediated translocation and eventual degradation of active M DM 2 by ARF [8][9][10], which results in the interconnection y 41 " u 24 , where y 41 is the only output from Module M 4 given by the expression The ODEs describing the dynamics of the module M 3 are depicted in Fig. 1(c).In this model, the mRNA of P T EN , represented by mRN A p , is transcribed at some basal rate p 5 , with a degradation rate constant d 4 .The transcription of mRN A p is also governed by the input u 31 , which describes the transcriptional activation of P T EN by p53 [4,6,11] resulting in the interconnection where y 12 is an output from Module M 1 given by the expression mRN A p is translated to P T EN with rate constant p 6 , and degraded with rate constant d 5 .P T EN mediates the de-phosphorylation of P IP 3 into P IP 2. The total concentration of P IP 2 and P IP 3 is given by the constant tot 1 .The rate of phosphorylation of P IP 2 is given by the term b 6 ¨f6 ptot 1 ŕP IP 3sq, while the rate of dephosphorylation of P IP 3 is given by b 7 ¨rP T EN s ¨f7 prP IP 3sq.P IP 3 then mediates the phosphorylation of AKT to AKT a (active AKT ).The total concentration of AKT and AKT a is given by the constant tot 2 .The rate of phosphorylation of AKT is given by the term b 9 ¨rP IP 3s f 9 ptot 2 ´rAKT asq, while the rate of dephosphorylation of AKT a is given by b 8 ¨f8 7prAKT asq [6,12].
The ODEs describing the dynamics of module M 4 are depicted in Fig. 1(d).In this model, p 7 is the basal rate of production of ARF , with u 41 " rD 2 s the extrinsic rate of production that represents the activation of ARF in response to DNA damage through transducers like P ARP 1 [13] and E2F 1 [14] .d 6 is the rate constant for the degradation of ARF .ARF is also degraded by an extrinsic process mediated by u 42 , which represents a novel recently-discovered mechanism through which M DM 2 can mediate the degradation of ARF [15].This results in the interconnection where y 22 is an output of Module M 2 given by the expression 9 rp53s " p 1 `h1 pu 11 q ´d1 ¨rp53s ´u12 ¨f1 rp53s Since there are multiple known M DM 2 ´ARF binding sites [15], we assume that these two interactions occur at separate binding sites, which results in the degradation rates of M DM 2a given by b 11 rARF s f 11 rM DM 2as, and that of ARF given by b 10 rM DM 2as f 10 rARF s, taking on distinct values.
In Fig. 2 of the main text, we showed a simplified block diagram representation of the full network N v consisting of the interconnections of the modules shown in S1.In Fig. S2, we show a more detailed block diagram representation of N v .In addition, Fig. S3 shows the block diagram representations of the remaining different network configurations that we study in the main text.

A.2 Model parameters
The functions f i pxq, i P t1, 2, ¨¨¨, 10, 11u that appear in the models in Figure S1 are given by and the functions h i pxq, i P t1, 2u, which represent the kinase activities, are given by h i pxq " c i x n i x n i `kn i i .

A.2.1 Tristable network N h
The parameters for this network are given in Table S1.Most parameters are consistent with the values and ranges from [6], except for the parameters relating to the ARF -M DM 2 interactions and ARF reactions, that did not appear in [6].The parameters relating to the entry of damage into the core regulation network through the kinases (p 1 , p 4 ,c 1 ,c 2 ,k 1 and k 2 ) are chosen somewhat arbitrarily as was done in [6].Moreover, while [6] chose to use µM as their unit of measurement, we used arbitrary unites (A.U.) in line with what is done in experimental studies [16,17].Since the amplitude of responses are known to vary significantly between different cells [18], the ratios between the parameters play a far more crucial role in bringing about the dynamical behavior as opposed to the parameters themselves.The steady-state solution to the network with these parameters results in the bifurcation diagram 6(a), which has been normalized with respect to the p53 level at T 0 .

A.2.2 Bistable networks N 2 and N 3
The study of N 2 was performed with the same parameters as shown in Table S1, with the external input u 24 set to 0 (although the results could be generalized to any value of u 24 ).These parameters  In the main text, we showed simulations of the full network in response to varying levels of damage.We stated that the reversible bistable network is able to admit both pulsing for moderate damage and a monotonic switch in p53 to a high level for high damage, but that this switch was reversible.We demonstrated this through a simulation of the network N 2 with an initial #DSB " 450, and illustrate our result in Fig. S5.

B Damage sensing and repair model
The functions for Eqs. ( 1)-(3) of the main text are given by with the parameters given in Table S5.S5: Parameters altered for core regulation network models Network N1 and Ni.
As was noted in the main text, the dynamics of D 1 are assumed to be slow compared to that of the core regulation network, leading to the observed behavior.We simplified our simulations by linearizing g 2 in the following way:

C Algorithm to compute network equilibrium points
In this section, we describe the algorithm used to compute the equilibrium points of our networks and obtain the bifurcation diagrams.

C.1 The Steady-State Behavior of a module
The procedure that the algorithm employs involves gridding the equilibrium values of the input and output signals for each module to form a table known as the Steady-State Behavior (SSB) of a module [19].Typically, this table is obtained by computing the equilibrium values of the outputs for a range of constant inputs, and this special case of the SSB has been termed the Static Characteristic Function (SCF) [20,21]

C.2 A modular approach to computing the equilibrium point of a network
To solve for the equilibrium value of rp53s of the network N v as a function of rD 1 s and rD 2 s, we need to systematically eliminate the input and output signals flowing between the modules, also known as the latent variables in the network [19].We demonstrate this by focusing on the interconnection between Modules M 1 and M 2 as is shown in Fig. S2, where we see the associations y 11 " u 21 and y 21 " u 12 .To eliminate these latent variables, we identify the rows i and j in Tables S6 and S7 respectively such that rD 1 s i " rD 1 s j y i 11 " u j 21 y j 21 " u i 12 , leading to the Table S10.Mathematically, the construction of Table S10 S6 and Table S7 into the u12 vs y11 (green line) and y21 vs u21 (other lines) planes respectively, for different values of u24 and u22 with rD1s " 0.4.The gray intersection points correspond to the equilibrium values of u12 and u21 for the given u24, u22 and rD1s.
The next step is to identify the rows m and k from Tables S10 and S8 respectively such that The elimination of these latent variables results in Table S11, whose construction can be mathematically represented by the SSB B 123 " prD 1 s , rp53s , u 24 , y 22 q P R 4 ě0 | Dpu 22 , u 31 q P R 2 ě0 s.t.prD 1 s , u 31 , rp53s , u 22 , u 24 , y 22 q P B 12 , pu 31 , u 22 q P B 3 ( . Computationally, the construction of Table S11 is achieved by finding the intersection points between each of the lines mapping from u 22 to y 12 from Table S10, and each of the lines mapping from y 31 to u 31 from Table S8, for all combination of other variables, as outlined in Fig. S7.S10 and Table S8 into the u22 vs y12 (black line) and y31 vs u31 (other lines) planes respectively, for different values of rD1s with u24 " 0.0015.The gray intersection points correspond to the equilibrium values of u31 and u22 for the given u24 and rD1s.The bifurcation can clearly be seen as rD1s increases, with rD1s " 0.1 and rD1s " 0.5 corresponding to a single equilibrium point, while rD1s " 0.3 corresponds to three equilibrium points.Finally, we identify the rows n and l from Tables S11 and S9 respectively such that The elimination of these latent variables results in Table S12, which is the combined steady-state behavior of the entire network N v .Mathematically, Table S12 corresponds to the set B 1234 " prD 1 s , rD 2 s , rp53sq P R 3 ě0 | Dpu 24 , u 42 q P R 2 ě0 s.t.prD 1 s , rp53s , u 24 , u 42 q P B 123 , prD 2 s , u 24 , u 42 q P B 4 ( , which was precisely what we wanted to compute.Computationally, the construction of Table S11 is achieved by finding the intersection points between each of the lines mapping from u 24 to y 22 from Table S11, and each of the lines mapping from y 41 to u 41 from Table S8, for all combination of rD 1 s and rD 2 s, as outlined in Fig. S8. 9 rp53s " ¨¨M

Fig. 3 .
Fig. 3. Steady-state concentration of p53 as a function of [D1] shows a continuum of steady-states in both Ni and N1.Due to the negative feedback on p53, N1 admits lower p53 steady-state values than Ni.The p53 levels are normalized by their value when [D1] = 0.

Fig. 4 .
Fig. 4. (a) Bifurcation diagram of N2, illustrating the reversible bistable switch.Assuming the p53 level starts at a point in "A", increasing [D1] beyond T1 results in the p53 level making a switch to a point in "B".For the p53 level to revert to a point in "A", the [D1] needs to fall below T0.(b) The level of D1 in our simulation.(c) The corresponding p53 response contrasting the behaviors of the monostable network Ni to the bistable network N2 as [D1] is varied.While the bistable switch has memory, the monostable switch does not.

Fig. 6 .
Fig. 6.(a) Bifurcation diagram of the network Nv, illustrating the tristable switch.The switching behavior of the bistable network is retained.However, there is an additional threshold T2 which when [D1] crosses, results in an irreversible switch of the p53 level to the set of points "C".(b) Simulated D1 levels.[D2] is kept fixed.(c) The resulting p53 behavior shows that the tri-stable network can behave both as a reversible and an irreversible switch for different values of maximal [D1].The green line represents the reversible switch, while the pink line represents the irreversible switch.

Figure S3 :
Figure S3: Modular block-diagram representations of the various p53 pathway configurations found in the evolutionary data.The specific ODEs describing each module are shown in Fig. S1.(a) Network Ni (b) Network N1 (c) Network N2 (d) Network N3 Modular block-diagram representations of full network Nv.

Table S3 :
Parameters altered for the core regulation network model in the Network N3

Table S8 :
of the module.It is worth noting that these tables could be Table of equilibrium values for M3.The superscript k denotes the k th row of the table.

Table S9 :
Table of equilibrium values for M4.The superscript l denotes the k th row of the table.

Table S10 :
corresponds to computing the SSB set B 12 , where B 12 " prD 1 s , y 12 , rp53s , u 22 , u 24 , y 22 q P R 6 ě0 | Dpu 21 , u 12 q P R 2 ě0 s.t.prD 1 s , u 12 , u 21 , y 12 , rp53sq P B 1 , prD 1 s , u 22 , u 24 , u 12 , y 22 , u 21 q P B 2 ( .Computationally, the construction of Table S10 is achieved by finding the intersection points between each of the lines mapping from u 12 to y 11 from Table S6, and each of the lines mapping from u 21 to y 21 from Table S7, for all combination of other variables, as outlined in Fig. S6.rD 1 s rp53s u 22 u 24 y 22 y 12 Table of equilibrium values when modules M1 and M2 are combined.

Table S11 :
Table of equilibrium values when modules M1, M2 and M3 are combined.