Modeling and controlling the two-phase dynamics of the p53 network: a Boolean network approach

Although much empirical evidence has demonstrated that p53 plays a key role in tumor suppression, the dynamics and function of the regulatory network centered on p53 have not yet been fully understood. Here, we develop a Boolean network model to reproduce the two-phase dynamics of the p53 network in response to DNA damage. In particular, we map the fates of cells into two types of Boolean attractors, and we find that the apoptosis attractor does not exist for minor DNA damage, reflecting that the cell is reparable. As the amount of DNA damage increases, the basin of the repair attractor shrinks, accompanied by the rising of the apoptosis attractor and the expansion of its basin, indicating that the cell becomes more irreparable with more DNA damage. For severe DNA damage, the repair attractor vanishes, and the apoptosis attractor dominates the state space, accounting for the exclusive fate of death. Based on the Boolean network model, we explore the significance of links, in terms of the sensitivity of the two-phase dynamics, to perturbing the weights of links and removing them. We find that the links are either critical or ordinary, rather than redundant. This implies that the p53 network is irreducible, but tolerant of small mutations at some ordinary links, and this can be interpreted with evolutionary theory. We further devised practical control schemes for steering the system into the apoptosis attractor in the presence of DNA damage by pinning the state of a single node or perturbing the weight of a single link. Our approach offers insights into understanding and controlling the p53 network, which is of paramount importance for medical treatment and genetic engineering.


Introduction
It has been recognized that p53 acts as the 'guardian of the genome' by inducing cell-cycle arrest, senescence, and apoptosis in response to oncogene activation, DNA damage, and other stress signals [1,2]. This is validated by the loss of p53 function in many human tumors as a result of mutations in p53 or inactivation of the p53 signal transduction pathway [3][4][5][6]. The regulatory network centered on p53 includes many positive and negative feedback loops, and p53 is a hub point in several signal transduction pathways [7]. The DNA binding cooperativity of p53 modulates the decision between cell-cycle arrest and apoptosis in response to DNA damage by triggering two-phase dynamics [8]. Therefore, when a cell suffers from DNA damage, the cell cycle will be suspended, and DNA repair machines will be engaged. In response to DNA damage, the cells will display either G1 arrest or G2 arrest triggered by p53 at the early stage. Direct transfection of wild-type p53 into either cells with wild-type p53 or cells with mutant p53 will cause G1 arrest [9][10][11][12]. During the cycle arrest, the injured cell repairs itself. However, when the DNA damage is irreparable, p53 will trigger apoptosis by releasing caspase-3, inducing the destruction of the injured cell [13][14][15][16]. In these biological processes, p53 remains activated in the sensing pathway after DNA damage [17].
Quite recently, continuous dynamical models have been proposed for exploring the dynamics of the p53 network and uncovering the p53-induced mechanism of deciding cell fate [8,[18][19][20][21][22][23]. Despite successes in reproducing the two-phase dynamics, we continue to lack a comprehensive understanding of the dynamics and biological function of the p53 network. The difficulty stems from the complexity of a dozen nonlinear differential equations associated with several adjustable parameters in the continuous model. It is computationally exhaustive to estimate the parameter values and difficult to fully uncover the nonlinear dynamics of the complex dynamical system, such as the stability pertaining to parameter values and the controllability of the collective dynamics [24,25]. Therefore, we exploit a Boolean network model as a simplified alternative approach to modeling the p53 network with respect to its twophase dynamics in response to DNA damage [26][27][28][29][30][31]43]. The much smaller parameter space of the Boolean network, compared to the continuous model, allows us to exhaustively test all possible configurations of the interaction strengths of links and the degradation rate of proteins to identify an optimal model for producing the two-phase dynamics [43]. In this model, the repair phase and the apoptosis phase are captured by two different Boolean attractors accompanied by their attraction basins in the state space, where the basin size of each attractor represents the probability of entering the attractor from an arbitrary initial state. Further, the optimal model is valuable for measuring the significance of links in the network by perturbing the weights of links or by removing them. Based on their influences on the two-phase dynamics, we find that the links can be either critical or ordinary rather than redundant, indicating that the p53 network is vulnerable to attack or large mutations due to the absence of redundant links, but it has a certain degree of tolerance to small mutations at ordinary links. The paradox between vulnerability and robustness can be explained by evolutionary theory. The Boolean network model also facilitates the formulation of a practical control method to guide cell fate by pinning the state of a single gene or perturbing the weight of a single link. Two critical nodes, Wild type p53-induced phosphatase 1 (Wip1) and Murine double minute 2 (Mdm2), are identified by implementing the control scheme to all nodes and links separately. The fact that the two critical nodes are subject to two negative feedback loops accounts for their effectiveness in control. Our approach offers a comprehensive understanding of the dynamics of the p53 network and may provide theoretical clues for the treatment of p53-based tumors.

p53 network
We established the topology of the signaling network centered on the tumor suppressor p53, based on previous experimental findings of transcriptional pathways and the results of continuous dynamical models, as shown in figure 1. When DNA is damaged, ataxia telangiectasia mutated (ATM) will be activated by phosphorylation to ATM* [32]. Then, p53 will also be activated to p53* by phosphorylation [33], which is regulated by ATM*. p53* phosphorylated only at Ser-15 and Ser-20 (p53a) induces cell-cycle arrest, whereas p53* further phosphorylated at Ser-46 (p53k) can induce apoptosis [34]. p53* induces the expression of p53dependent damage-inducible nuclear protein 1 (p53D) and Mdm2 [35][36][37][38]. Subsequently, Mdm2 promotes the degradation of p53* [39][40][41] and p53D regulates the conversion from p53a The signaling network centered on p53*. The network consists of 13 nodes, 22 links, and 6 self-loops, where arrow links represent promotion, arrows with bar heads represent inhibition, and yellow arrows with bar heads represent self-degradation. The thicknesses of links and self-loops represent their weights. Thick links correspond to weight 2, and thin links correspond to weight 1. Self-loops are weighted as follows: to p53k. Wip1, which is triggered by p53a, helps to convert p53k to p53a. Meanwhile, Wip1 has a substantial inhibitory effect on ATM* by dephosphorylating ATM* to ATM [42]. The p53* signaling network has 13 nodes (proteins), 6 negative self-loops [43] representing self-degradation, and 22 links among nodes, including 7 inhibition links and 15 promotion links. Arrow links represent promotion, links with bar-heads represent inhibition, and self-loops colored yellow represent self-degradation. The links between p53* and p53 a and between p53* and p53 k indicate that p53* is composed of p53 a and p53k. Further, three feedback loops, p53-Mdm2, ATM-p53-Wip1, and p53-PTEN-Akt-Mdm2, are included in the network due to their significant roles in activating p53 [44][45][46].

The Boolean dynamics
In a standard Boolean network, each node has two states, S i = 1 and S i = 0, representing the active and inactive states, respectively, of a protein or gene. The state of the entire network is determined by all nodal states S i (i.e., = S S S S S { , , ,..., } n 1 2 3 , where n is the number of nodes in the network). Because the state of each node is binary, a network with n nodes has 2 n types of overall states. In general, a certain number of Boolean attractors exist in the state space, determined by the Boolean rules. There are two possible types of attractors: fixed-point and limited-cycle. For a fixed-point attractor, there is a single state in the attractor, and for a limitedcycle attrator, there are a certain number of states that form a finite directed cycle. Once the Boolean system evolves into an attractor from an arbitrary state, the system will be trapped in the attractor and cannot escape. Another important concept in Boolean dynamics is the basin of an attractor, defined as the set of initial states that flow into the attractor. Each of the initial states is often denoted by a vertex in the state space, so that the basin of the attractor can be characterized by a set of vertices connecting to an attractor.
To better mimic the real p53 network, we slightly modified the setting of the standard Boolean networks, as described above. Specifically, we assume that the state of ATM*, S * ATM , ranges from 1 to 10, reflecting different degrees of the concentrations of ATM* that depend on the amount of DNA damage. Because some empirical evidence indicates that a great deal of p53 is phosphorylated and associated with large amounts of ATM, we also set the states S S * , p53 p53a , and S p53k of p53*, p53a, and p53k in the range between 0 and 10. Moreover, because distribution involves PIP and Akt, the states S Akt , S * Akt , S PIP2 , and S PIP3 are set in the range 0 to 2. Using the biochemical reactions discovered in experiments, we have the following restriction equations: where T ATM , T PIP , and T Akt denote the constant total amount and are set at 10, 2, and 2, respectively. The states of the other nodes are binary, either 0 or 1.
Next we introduce the initial state space of our modified Boolean network, which is imperative for quantifying the basins of the attractors. Due to the restriction (equation (1)) of some nodes and the initial existence of some inactive nodes, their node states are not totally free in the initial state. In other words, the initial state space is a subspace of the entire state space composed of all possible combinations of nodal states. Specifically, the initial states of some nodes are entangled (i.e., for the pairs (Akt, Akt*) and (PIP2, PIP3), there are only three combinations of initial states: (0, 2), (1, 1) and (2, 0)). In addition, because the concentrations of p53*, p53a, p53k, and p53D at the beginning of DNA damage are negligibly small, we set their initial states as 0. Other than these nodes, the states of Wip1, Mdm2, and PTEN are initially free and can be independently set to be either 0 or 1. Moreover, we set the initial state of ATM* as a constant to reflect the amount of DNA damage, with upper and lower bounds of 10 and 0, respectively, during evolution. The initial state of ATM depends on that of ATM*, according to equation (1). Thus, the initial state space of the modified Boolean network contains 72 states.
Starting from a given initial state, the Boolean network evolves according to the following Boolean rule: for an inhibition link, and = a 0 ij for the absence of link from j to i. σ i represents the regulation capability of nodes. With respect to some experimental evidence and without loss of generality, we set σ for nodes with high state values and upper bounds that are relatively small (i.e., σ = 0.3 i for p53*, p53a, p53k, ATM*, and ATM). For the other nodes, we set σ i to be 1. S i max stands for the upper bound of the nodal state of i.
For the nodes that are not involved in any distribution processes as described in equation (1) i i including Wip1, p53D, PTEN, Mdm2, and p53*. For the other nodes that are involved in the distribution processes, is determined by intermediate states, the total amount and distributions, as follows: where⌊⌋ · represents the rounding operation because of the integer values of all nodal states, and T i is the total amount of S i and its counterpart, S , i c in equation (1). For example, if i refers to p53a, iʼs counterpart is p53k, and vice versa, and T i is , according to equation (1). For ATM, S ATM and S * ATM are counterparts and = = T T 10 i ATM . A similar distribution rule is also applied to PIP and Akt.
After all of the nodal states are updated according to the Boolean rule (equation (2)) and the distribution rule (equation (4)), we subsequently consider the degradation of the states of some nodes without inhibition links pointing at them. The set of nodes includes ATM, Wip1, Akt*, p53D, PTEN, and PIP2. The degradation of their states is captured by weighted, negative self-loops of the nodes, as shown in figure 1. The degradation is implemented as follows. If a node i with a self-loop is activated at t and remains activated for τ time steps, at the step τ + 1 after t, S i becomes − S 1 i . Thus, the degradation process is approximated in this manner. It is imperative to implement degradation to the nodes that are not negatively regulated by the other nodes; otherwise, insofar as they are activated, they will not return to silence forever, which is clearly different from real situations.
Prior to presenting our results, we need to define the repair attractor and apoptosis attractor and their attraction basins. Although our Boolean network model is not standard, the fact that there are no stochastic effects ensures the existence of attractors in the form of either fixed points or limited cycles. Insofar as the systemʼs state enters an attractor along an arbitrary trajectory, the state will be trapped in the attractor, and the system will be unchanged or exhibit periodic behavior. Based on this a priori knowledge, we explore the systemʼs dynamics by relying on the attractors. In particular, an attractor is determined to be a repair attractor if all of the vertex-states in the attractor are subject to the repair state. Note that here, the vertex-state refers to the state composed of the states of 13 nodes in the state space, and that a vertex in an attractor corresponds to a combined state in the state space. A vertex in an attractor is said to be subject to the repair state if > S S p53a p53k in the vertex. The definition of the repair state can be explained in terms of the biological functions of p53a and p53k. p53 a plays the key role in triggering cell-cycle arrest that is necessary for repairing DNA damage. In contrast, p53 k is the key to causing the destruction of an injured but irreparable cell. Thus, > S S p53a p53k indicates that the cell is likely to be repaired, rather than killed. Analogously, in an apoptosis attractor, all vertex-states are subject to the apoptosis state, in which < S S p53a p53k . In fact, a third type of attractor may exist, in which repair and apoptosis vertices are mixed. We call that attractor a mixture attractor.
The attraction basin associated with an attractor is defined as the set of initial states that eventually flow into the attractor. In other words, for each vertex in the set, there exists a trajectory to the attractor. The size of an attraction basin is defined by the number of vertices in the basin. We denote the sizes of repair basins and apoptosis basins as N R and N A , respectively. The sizes of repair and apoptosis basins characterize the respective probabilities of entering the two attractors. As an extreme example, if N R equals 72, which is the number of vertices (states) in the entire state space, the DNA damage is totally reparable. In contrast, if N A equals 72, death is the only fate. Therefore, N R and N A can be used to quantify the fate of a cell in response to different amounts of DNA damage.

Results
Given the topology of the Boolean network centered on p53, we aim to determine the configuration of link weights (interaction strengths) and the weights of self-loops (degradation rates) based on the two-phase dynamics mapped into the attractors and attraction basins in the modified Boolean dynamics. This is implemented by exhaustively testing all possible configurations, which is feasible because of the simplicity of our modified Boolean network model, compared to continuous dynamical models. The ascertained weighted Boolean network allows us to identify the significance of each link and classify the links via attacks and perturbations of the links. Based on the importance of each link, we are able to devise practical control schemes by pinning the state of a critical node or perturbing the weight of a critical link to steer the system to the desired attractor pertaining to the desired final state.

Two-phase dynamics of p53
First, we attempt to determine the configuration of the weights of the links and self-loops in order to reproduce the two-phase dynamics in response to DNA damage. Specifically, for small amounts of DNA damage induced by radiation, the damage is believed to be repaired by the p53 network, corresponding to entering the repair attractor. In contrast, if the DNA damage is too severe to be fixed, the system enters the apoptosis attractor, accounting for the eventual death of the cell. Moreover, as the amount of DNA damage increases, the probability of being successfully repaired should decrease, accompanied by an increase in the death rate. Taken together, these phenomena offer criteria for identifying an optimal configuration by exhaustively enumerating all possible configurations. However, it is generally computationally prohibitive to achieve this goal, because of the continuous characteristics of the weights of links and self-loops. This problem can be solved by coarse-graining the continuous values. In other words, we assume that the weights of links or self-loops are integers, which we will show are effective for recovering the two-phase dynamics and for modeling the function of the p53 network. This approximation leads to a considerable reduction of the solution space, allowing us to find useful configurations associated with empirically observed biological functions.
Specifically, for a given configuration, we adjust the initial level of ATM* that captures the amount of DNA damage to explore the dependence of the two types of attraction basins on DNA damage. By examining all possible configurations with respect to the basins, we eventually found a configuration that meets the criteria, as shown in figure 1. We can see that the basin area of the apoptosis attractor is a monotonic increase function of the level of ATM*, as shown in figure 2. In particular, for minor DNA damage (e.g., ATM* = 1), there is no apoptosis attractor; however, for severe DNA damage (e.g., ATM* = 10), the apoptosis basin occupies the entire state space, indicating that the cell is irreparable and death is the only possible fate. , demonstrating that the apoptosis attractor is absent and the cell is reparable. For severe damage (e.g., ATM* = 10), we find = N 72 A , indicating that the initial state space is fully occupied by the apoptosis basin and that death is the only fate. Moreover, N A is a monotonic increase function of ATM*, consistent with empirical findings and common sense. Figure 3 shows the two attractors and their basins for different ATM* levels. We find that both attractors exist in the form of limited cycles. For ATM* = 1, there is a repair attractor and a mixture attractor, but no apoptosis attractor. At ATM* = 2, an apoptosis attractor arises, with a single initial state entering the attractor. The scenario of ATM* = 3 is similar to that of ATM* = 2. At ATM* = 4, there are two repair attractors, and the basin of the apoptosis attractor contains only two initial states. At ATM* = 6, the basin of the apoptosis attractor becomes much larger than that of ATM* < 6, indicating a higher probability of eliminating the cell. As However, because the system is deterministic without any stochastic effect, once entering an attractor, the trajectory will be trapped and never escape. Due to the complexity of the state space, we omit the trajectory from initial states (vertices) to attractors, and linking vertices to their attractors, directly with each link pointing at a vertex in the cycle. ATM* continuously increases from 7, the expansion of the apoptosis basin is sustained, along with the shrinking of the repair basin and the reduction of the number of repair attractors. At ATM* = 9, repair attractors either disappear or are replaced by a mixture attractor. At ATM* = 10, there is only a single apoptosis attractor, resulting in the death of the irreparable cell. All of these findings pertaining to the two types of attractors for different ATM* levels are consistent with empirical observations, validating the configuration of p53 networks we have found.
The validated p53 network allows us to explore the importance of links and to steer the system to enter the desired phase by pinning some critical nodes or perturbing the weights of some critical links.

Significance of links
All of the links can be classified into three categories according to their significance in maintaining the biological dynamics of the p53 network: (i) critical links, if perturbing their weights induces changes in the attraction basins in the state space of either ATM* = 1 or ATM* = 10; (ii) ordinary links, if their removal induces changes in the attraction basins for ATM* = 1 or ATM* = 10, but perturbing their weights has no influence; and (iii) redundant links, if neither perturbing their weights nor removing them leads to any change in the attraction basins for ATM* = 1 or ATM* = 10. In other words, the significance of a link can be measured by the sensitivity of the two-phase dynamics to its manipulation. Critical links are named 'critical' because small perturbations to their weights affect the two-phase dynamics and ruin the biological function of the p53 network. In contrast, the dynamics are unaffected by either perturbing or removing redundant links, hence the name 'redundant'. Ordinary links are intermediate between critical and redundant links.
In our weighted p53 network, the weights of promotion and inhibition links are subject to low and high levels, represented by 1, 2, and −1, −2, respectively. The perturbation of a link weight refers to switching the level of the link weight, for instance, from the original weight of 1 to 2, or from the original weight of −2 to −1. Interestingly, by implementing weight perturbation and link removal, we find 14 critical links and 8 ordinary links, but no redundant links, as shown in figure 4. This finding suggests that the p53 network is vulnerable to attacks or large mutations due to the absence of redundant links. However, we should keep in mind that the network model used in the analysis is a simplified model, and some redundant pathways have already been ignored in the modeling. Even for this simplified network model, it is tolerant of small mutations occurring at ordinary links. In the sense of the evolutionary theory, small beneficial mutations are acceptable and accumulable, but large mutations are likely to lead to death. Another striking finding is that all inhibition links are critical, which may have implications for biological experiments.

Controlling the two-phase dynamics
Based on the weighted p53 networks, one important task is to devise a practical control strategy to steer the biological dynamics, with a particular interest in reaching apoptosis, which is important for blocking the formation of tumors resulting from the survival of irreparable cells. We attempt to pin the state of a single node or perturb the weight of a single link to drive a system with DNA damage to the apoptosis phase. We have tested each of the nodes and links and identified two effective nodes, Mdm2 and Wip1, and one effective link, Mdm2→p53. As shown in figure 5, by pinning the state of Mdm2 to zero, the apoptosis basin expands remarkably for different amounts of DNA damage (from ATM* = 1 to ATM* = 10), compared to the original two-phase dynamics in the absence of pinning control. The expansion of the apoptosis basin indicates a higher probability of entering the attractor and triggering the destruction of the cell. Perturbing the weight of Mdm2→p53 subject to Mdm2 is also feasible for achieving the apoptosis phase for different levels of ATM* (e.g., Mdm2→p53= −1). In particular, when we block the interaction (i.e., Mdm2→p53 = 0), the repair attractor vanishes and apoptosis becomes the only fate.
Controlling Wip1 is more effective than controlling Mdm2, as shown in figure 5. There are two manners of controlling Wip1: pinning the state of Wip1 to zero or adjusting the degradation rate of Wip1. We can see that the apoptosis attractor becomes the exclusive attractor by pinning , regardless of the amounts of DNA damage quantified by ATM*. Adjusting the weight of the self-loop of Wip1 to be 1 or 2 gives rise to similar results to the pinning control of Wip1ʼs state with the full occupation of the initial state space by the apoptosis basin. These results demonstrate that self-destruction of a cell can be ensured merely by controlling Wip1, which is the most critical node for triggering the apoptosis process, according to our findings.
In fact, we have tested all nodes and links, but we found no valid nodes and links other than those described above. This may be attributed to the fact that Mdm2 and Wip1 are the key nodes subject to two negative feedback loops, ATM*-p53*-Wip1 and p53*-PTEN-Akt-Mdm2, in the p53 network.

Discussion
Reproducing the two-phase dynamics of p53. We have developed a Boolean network model to recover the two-phase dynamics of the p53 network in response to DNA damage [47]. The topology of the network is established based on both previous experiments on biological pathways and continuous dynamical models in the literature. Our main contribution lies in enumerating all possible interaction weights and degradation rates of the nodes to identify an optimal configuration of the Boolean network with regard to the two-phase dynamics. Despite the simplicity of our Boolean network model in contrast to the continuous dynamical models, the two-phase dynamics as an emergent phenomenon in a self-organized complex system can be well captured by our model. More importantly, the Boolean network model allows us to explore the dynamics of the p53 network are a comprehensive manner in terms of Boolean attractors and their attraction basins in a finite state space. Subsequently, our model enables the identification of significant nodes and links in the network based on the influence of their changes to the repair and apoptosis attractors, which also facilitates the design of practical control schemes to guide the system toward a desired fate. Our model reveals the relations between the p53 network topology and the corresponding physiologic function, including cell self-repair and apoptosis. In 2012, Bashan et al found that there is a specific network structure underlying and characterizing each physiological state, which implies that network topology and function interplay robustly [48,49].
Classification of links. Each link in the p53 network has been classified into one of three categories-critical, ordinary, and redundant-according to its importance to the emergence of two-phase dynamics. We have found 14 critical links among the 22 total links. Perturbing or removing the critical links leads to violation of the two-phase dynamics. The non critical links are all ordinary links. The removal of the ordinary links ruins the two-phase dynamics, but perturbing their weights has no affect. Thus, the p53 network has no redundant links and is irreducible and vulnerable to external attacks on its links. However, the network is resilient to small mutations at ordinary links, without much effect on the two-phase dynamics. From the perspective of evolutionary theory, the p53 network has typical characteristics, in the sense that it permits small mutations that can be accumulated during evolution, but suffers from large mutations that usually lead to death. Recently, Bartsch and Ivanov investigated the network of physiologic interactions between the brain, cardiac, and respiratory systems across different sleep stages, and found that there was a strong relationship between network connectivity, patterns in the strength of network links, and physiologic function [50].
Controlling the two-phase dynamics. Based on our Boolean model of the p53 network, we have devised a practical control method that can be implemented at a single node or link to force the system into the apoptosis attractor in the presence of DNA damage. In particular, we have identified two important nodes, Wip1 and Mdm2, that are the most effective for this and perturbing the weight of link Mdm2→p53* to be −1 or its weight to be 0, and adjusting the negative self-feedback (weight of self-loop) of Wip1 to be 1 or 2. Each type of control scheme is independent of the others (i.e., the apoptosis phase can be achieved by controlling a single node or a single link, which is feasible from a practical perspective).
control. We can carry out this control by either pinning the state of Wip1 or Mdm2 at a fixed value or by perturbing the weight of the link, Mdm2→p53. In fact, successfully achieving apoptosis by controlling Wip1 and Mdm2 is consistent with some empirical evidence about the relationship between these two nodes and apoptosis. For instance, previous experiments have demonstrated that the absence of Wip1 or Mdm2 causes apoptosis in the system. The negative regulator, Wip1, has an important role in inhibiting ATM, and the inhibition of ATM may lead to dissociation of ATM from damage sites, thereby allowing repair proteins to access the damage [46,51,52]. Thus, the absence of Wip1 will hinder the repair process, inducing cell death. Meanwhile, the absence of Mdm2 induces the malfunction of its transcriptional factor, p53*. Thus, the apoptosis states may be ascribed to the mutation of p53, which has been reported in the literature [53][54][55].
In summary, our Boolean network model for the p53 network has several advantages and attractive characteristics, including simplicity compared with continuous dynamical models, successful reproduction of the two-phase dynamics in response to DNA damage, the classification of links based on their significance with evolutionary meanings, and the design of practical control strategies with empirical support. Although continuous dynamical models, as proposed in [8,18], can describe more details of system, a considerable number of parameters will lead a huge phase space, which is impossible to scan completely. However, it is easy to grasp the nature of a system described by a Boolean model. On the other hand, the limitation of our model lies in the fact that the state space grows exponentially as the size of the network increases. It is impractical to apply our model to large networks. Taken together, we believe our approach offers new insights on understanding and controlling the dynamics of a p53 network by means of Boolean dynamics, and it also has implications for cancer treatment and genetic engineering.