A multiscale model of the role of microenvironmental factors in cell segregation and heterogeneity in breast cancer development

We analyzed a quantitative multiscale model that describes the epigenetic dynamics during the growth and evolution of an avascular tumor. A gene regulatory network (GRN) formed by a set of ten genes that are believed to play an important role in breast cancer development was kinetically coupled to the microenvironmental agents: glucose, estrogens, and oxygen. The dynamics of spontaneous mutations was described by a Yule-Furry master equation whose solution represents the probability that a given cell in the tissue undergoes a certain number of mutations at a given time. We assumed that the mutation rate is modified by a spatial gradient of nutrients. The tumor mass was simulated by means of cellular automata supplemented with a set of reaction diffusion equations that described the transport of microenvironmental agents. By analyzing the epigenetic state space described by the GRN dynamics, we found three attractors that were identified with cellular epigenetic states: normal, precancer and cancer. For two-dimensional (2D) and three-dimensional (3D) tumors we calculated the spatial distribution of the following quantities: (i) number of mutations, (ii) mutation of each gene and, (iii) phenotypes. Using estrogen as the principal microenvironmental agent that regulates cell proliferation process, we obtained tumor shapes for different values of estrogen consumption and supply rates. It was found that he majority of mutations occurred in cells that were located close to the 2D tumor perimeter or close to the 3D tumor surface. Also, it was found that the occurrence of different phenotypes in the tumor are controlled by estrogen concentration levels since they can change the individual cell threshold and gene expression levels. All results were consistently observed for 2D and 3D tumors.

Table 1.Network Boolean rules for breast cancer development.
Gene Boolean rules TP53(x 1 ) x 2 ∧ x 4 ∧ ¬ x 5 ∧ ¬ x 6 ∧ x 8 ¬ x 9 ATR(x 2 ) x 4 ATM(x 3 ) x 1 BRAC1(x 4 ) x 2 ∧ ¬ x 8 ∧ ¬ x 10 HER2(x 5 ) x 5 MDM2(x 6 ) x 1 ∧ ¬ x 2 ∧ x 3 ∧ x 8 CHEK1(x 7 ) x 9 AKT1(x 8 ) x 2 ∧ x 4 ∧ ¬ x 8 P21(x 9 ) ¬ x 1 ∧ x 5 CDK2(x 10 ) x 9 ∧ ¬ x 10 represent, proteins, RNA, morphogenes, etc.The links that join the nodes represent positive or negative regulations.The time evolution of the state of node x i , is determined by the Boolean function Ψ i that depends on the states of its regulators x k , (with k = i), located at other nodes in the network as follows The Boolean function Ψ i is discrete and follows the algebraic rules of Boole's axiomatics.The stationary states of the dynamical mapping, called attractors, are determined by the condition x i (t) = x i (t + ∆t), for any ∆t > 0 that lead to activation or inhibition patterns defining a homeostatic state of the system.In the present context these attractor states represent normal, precancer or cancer cellular phenotypes.
An exhaustive exploration of the attractors of the breast cancer GRN was carried out using the R BoolNet package (version 4.0.2).The Knock-Out (KO) and Over-Expression (OE) functions of each node were represented by a 0 (KO) or 1 (OE) value.The action of the OE and KO of each node was also explored.The tables that show all the identified system attractors are presented in Figs 1 and S2.Many repeating attractors were found which were classified in groups and arranged in a single table that is presented in Fig 1 .These results are important because they facilitate the analysis of the GRN Boolean interactions and dynamics.Network attractors for breast cancer (see in Figs 1 and 2) can be summarized as follows: (a) The first attractor of Wild-Type character is identical for genes TP53-KO, ATR-KO, ATM-KO, BRCA-KO, HER2-KO, MDM2-KO, CHEK1-KO, AKT1-KO, P21-KO and CDK2-KO; (b) The second attractor of Wild-Type character is identical for genes TP53-KO, ATR-OE, ATM-KO, BRCA1-OE, HER2-KO, MDM2-KO, MDM2-OE, CHEK1-OE, AKT1-KO, P21-KO and CDK2-KO; (c) The third attractor of Wild-Type character is identical for genes TP53-KO, ATM-KO, HER2-KO, MDM2-KO, CHEK1-KO, AKT-KO, P21-KO and CDK2-KO; (d) The fourth, and last attractor of Wild-Type character is identical for genes TP53-KO, ATR-KO, ATM-KO, BRCA1-KO, HER2-OE, MDM2-KO, CHEK1-KO, CHEK1-OE, AKT1-OE and P21-OE; Finally, (e) the remaining attractors represent some other specific functions of some genes and are presented in the tables shown in Fig 2 .The complete list of attractors is shown in the tables 4-42 at the end of this Supplementary Material.In each table, the R BoolNet package provided each attractor frequency of success.As indicated at the tables bottom, an attractor can show up with a given frequency, because of this, in some cases there is a higher probability that certain attractors appear over others; this is the reason for which the dynamics network shows only bistable states.These results are used as a guide to built up and  Breast cancer attractors (Attr) considering Over-Expression (OE) and Knock-Out (KO) of specific genes.The upper table lists all the GRN attractors when gene TP53 is Over-Expressed (TP53-OE).On the other hand, the lower table lists the GRN attractors when genes ATM, MDM2 and AKT1 are either OE or KO. the tumor microenvironment change gene expression levels by activating or inhibiting genes which affect cell lineage [6].To model the effect of estrogen and oxygen in the the GRN dynamics we have used the gene interaction kinetics introduced previously in [2].

Continuous model of GRN
-See Eqs (1) in the main text-.In this approach, the values of the state variables and parameters that control the GRN dynamics -shown in Fig 1 in the main text-vary continuously within a certain range as a result of the estrogen and oxygen concentrations gradients that lead to diverse genetic states.Here we assume that the active and inactive part of a gene state are important for the GRN dynamics, then we can write, X a + X h = X 0 where X a represents the gene active part and X h is the gene inactive part, and X 0 represents the gene state.Bearing in mind that that the gene state remains unchanged and is constant during the GRN evolution, then we can describe the gene state dynamics in terms of relative variables (X a + X h )/X 0 = 1.By writing x a = X a /X 0 or equivalently x h = X h /X 0 we can express either x a = 1 − x h or x h = 1 − x a in the GRN dynamics.This approach has been used to analyze the multistability of the dynamics of some genetic circuits [7][8][9].Therefore, the temporal evolution of the state x i of gene i is determined by its network kinetic interactions described by the following system of coupled ordinary differential equations: The interaction rate constants between gene states and microenvironment are represented by the parameters α i .The action of the microenvironmental agents on gene activation is represented by y i while the activation-inhibition threshold parameter is represented by η i and the self-degradation constant is µ i .The self-regulatory positive and negative feedback gene interactions are described by the following equation, where λ i is the activation constant and δ i is the activation threshold.The genetic strength is represented by the Hill function H: where β j is the activation-inhibition parameters, ν j is the action rate at which gene j affects the state of gene i, the parameters k j represent the sigmoidal threshold function, γ j is the Hill coefficient which controls the steepness of the sigmoidal function and represents the cooperativeness of the regulatory transcription factor binding to genes.One should note that the Hill's function rewritten as in the second part of Eqs (4) suggests that the cooperativeness rate is regulated by the parameter k j .Finally, the value of γ j depends on the number of neighbor genes that contribute to the system cooperativeness to change the state of gene j.The system of differential Eqs (2) can be rewritten in terms of the dimensionless variables: xi = νi ki x i and δi = ki νi δ i , where ki = k 1/γi i .Then, the dimensionless equation is written as + j β j xi Ĥ(x j ; 1; 1; γ j ), with 1 ≤ i ≤ 10, November 13, 2023 5/39 where Ĥ(x j ; 1; 1; γ j ) = xγj j /(1 + xγj j ).This approach suggests that the stability of the system depends of fewer parameters.From this analysis we observe that there is a set of 52 parameters that are basic, but are not fully independent, that determine the dynamics system.However, this set of parameters can be further reduced if one considers the relevant time scales as: the gene interactions (β i , λ i ) and degradation of proteins (µ i ), the coexistence between the active and inactive parts of the genes (κ i ), and the unitary action rate at which a given gene affects others (ν i = 1).This leaves us with six free parameters to determine the multistability of the system which can be studied by following the standard bifurcation analysis.We also performed an extensive bifurcation analysis for the temporal scales represented by the parameters (α i ) and the threshold activation values (η i ) that are crucial to obtain the cellular phenotypes.In the following sections, we build up a multistability diagram of the dynamics by estimating the values of these parameters.

Parameters estimation
The analysis of the dynamics described by Eqs (2) was carried out by setting the temporal scale in terms of the parameters α i since it is the main temporal scale that defines the interactions between relative gene states and microenvironment.For simplicity, we assumed that the dynamics of all GRN genes occurs at the same time scale, that is, α i = α 0 .Because of this, we set α 0 ± 1, where the positive value corresponds to a gene activation and the negative value corresponds to an gene inhibition.The local concentrations of estrogen and oxygen, are represented by y i , and vary in the interval [0, 1].One should recall that the transport of these microenvironmental agents occurs by diffusion and it is described by the set of reaction-diffusion Eqs (1) in the main text.
The microenvironment action thresholds, η i , are chosen following some phenotypic plasticity rule and take values between 0 and 1.In this case, the values of η i depend on the metabolic and genetic activity in the neighborhood of the cell and are measured with a cooperative enzymatic function.The relation of the intrinsic and the extrinsic plasticities with the microenvironment action thresholds, η i is given by Eq (14) in the main text.To obtain the gene inhibition parameters values, µ i , we first considered that the average lifetime of a protein and its transcription time is 20 minutes [10,11] and that the average time of protein degradation is about 3 to 6 minutes [12,13].Thus, taking the proportion of these times and the values of α i , it was found that the self-degradation rates, µ i , should take values between 0.15 and 0.3.We chose the intermediate value µ i = 0.2 for all genes.
By assuming that the activation and inhibition of the GRN genes drive the system behavior, it is reasonable to consider that the values of the activation-inhibition parameters are of the same magnitude as the self-degradation rates.Thus, assuming a balance between these quantities and the degradation, all the kinetic parameters take the value β i = 0.2, except for some of them that are related to the most outstanding multistability control circuits in the GRN [14][15][16][17][18], in which case the parameter β i may change.On the other hand, to guarantee that the negative or positive self-regulation preserves the state of the gene, the value of the λ i should double the value of β i .Taking into account that only two genes HER2 and CDK2 present self-regulation, one can set λ 5 = λ 10 = 0.4 and the remaining are set to zero.For breast cancer, it has been found that HER2 tends to remain active and overexpressed during cancer development [19][20][21][22], while CDK2 is partially degraded [2,23].It does not degrade completely because it must keep the cell cycle initiation pathway active [23,24].Therefore, we chose the values δ 5 = 1 and δ 10 = 0.1.This choice of these values yields overexpression of HER2 and keeps CDK2 in a latent state.
Taking into account that the increase in estrogen receptor failure may be associated with a cooperative phenomenon [20,22] the self regulation of gene HER2 was multiplied November 13, 2023 6/39 by a Hill function.In this case, the Hill function models the gene affinity with itself and the possible binding sites under the cooperative action among other genes.Thus, we chose the values β5 = 1, ν5 = 0.5 and γ5 = 1.In the case of gene CDK2 no evidence has been found that there exist a self-cooperative mechanism that inhibits to itself.Because of this, the Hill function associated with CDK2 self regulation is H(x 10 ; β10 ; ν10 ; γ10 ) = 1.In all other cases, the capacities associated with the Hill functions are the same and take the value k j = 0.5.This latter value guarantees that both, the active and the inhibit states of the gene have the same probability of being achieved.The Hill's functions vary in the range 0 ≤ H(x j ; 1; ν j ; γ j ) ≤ 1 when the activation-inhibition parameters β j are normalized by the values of k j as can be seen in the Eq (4).This allows us to analyze the dynamics without writing the equations in terms of dimensionless variables.Finally, the cooperativeness is quantified with the Hill coefficients by counting the number of input interactions that each gene has in the GRN (see Fig 1 in main text).Thus, we chose the values: For many genes in the GRN it is poorly understood how a variation of the action rate constant ν i values modify their state.However, it is known that for some circuits, such as the TP53-MDM2-P21 the values of the action rate constants are crucial for the GRN behavior.Because of this, the action rate constants were chosen as ν j = 1 for most of the cases.Given the relevance of genes TP53, MDM2 and P21 in the appearance of cancer cells and apoptosis processes the values of these parameters were ν 1 = 1.5, ν 6 = 1.2 and ν 9 = 0.5.[8,[14][15][16][17][18]24].These values take into consideration the fact that the TP53 agent, which is considered the guardian of the genome, has a greater influence than the other genes.Moreover, it has been shown [8,14,24] that the TP53-MDM2 circuit influences the network equilibrium states and cell fate, where the affinity of breaking the activation-inhibition of MDM2 over TP53 is estimated to be 80%.In addition, these kinetic values guarantee that in the TP53-P21 circuit the apoptosis process is not completed until the CDK2-P21 pathway is activated, as has been shown in [10,24].The choice of ν 9 = 0.5 limits P21 processes while the value ν 1 = 1.5 speeds up processes related to TP53.Some studies have shown that genes can undergo changes in the average transcription factors of up to two orders of magnitude due to the action of some microenvironmental agents [24,25].In particular, it has been demonstrated that the transcription levels of MDM2 could be increased by up to two orders of magnitude [24].This suggests that the TP53 circuit should speed up its transcription to maintain network stability.Furthermore, it is known that variations in the transcription levels of MDM2 should compensate to maintain CDK2 activity [8,24], since it is essential for the control of the other parts of the network.Thus, the MDM2 and CDK2 parameter values were chosen to be β 6 = 2.0 and β 10 = 2.0, one order of magnitude larger than the initial values, β i = 0.2, indicated above.(See Table 2 for a summary of parameter values).
It is known that an excess of estrogen in the microenvironment increases the possibility of cancer development, since it does not only modify the cell cycle but also alters some signaling pathways [19][20][21][22].Therefore, the magnitude of parameter values that drive the HER2 gene should be the largest in the network, so that, we set β 5 = 20.On the other hand, the HER2-TP53-P21 circuit participates in DNA repair, cell cycle arrest, and apoptosis [2].Thus, it is expected that the HER2 repression towards TP53 and the activation of HER2-P21 must be compensated.Because of this, the TP53 activation parameter was set to about half of that corresponding to HER2 gene, that is, β 1 = 9, and the activation of gene P21 was set at β 9 = 3.These values are sufficient for the activation of P21, without affecting the expression of the gene TP53.A summary of the values of all parameters involved in Eqs ( 2) is presented in Table 2.

Multistability
To analize the multistability of the dynamical system described by Eqs (2), we used a set of different algorithms described in [7,9,[26][27][28].They involve the analysis of two dimensional sections of the phase space as well as a numerical continuity approach by varying the values of the system parameters to find out the system states and their stability, that is, the basins of attraction in phase space .To this end we used the R-script Grind.R (https://tbb.bio.uu.nl/rdb/grind.html& https://tbb.bio.uu.nl/rdb/grindR/grind.R version 7-10-20).This script uses a powerful R routines package to analyze the solutions of a system of ordinary differential equations (ODEs).The results of this intensive and extensive analysis yielded all attractors as well as wild-type, KO and OE gene states that were found with the Boolean model.In Figs 3 and 4, are shown the evolution of a normal cell genetic states for different values of the microenvironmental threshold parameter η i .It was found that for η i < 0.5, with i = 1, 2, ...10, the genetic states are conserved, while for η i > 0.5 there appear basins of attraction that are identified with precancer and cancer phenotypes.We also carried out a detailed analysis to estimate the values of the other parameters that led to a precancer phenotype.We found that for each one of the gene states the activation threshold value, should be in the range 0 < δ i < 0.8.These results are in complete agreement with most of the thresholds activation values found in [2], where a stochastic model was analyzed and the thresholds were expressed as percentages.
On the other hand, it was found that for each one of the states, gene inhibition was obtained for a threshold value η i = 0, whereas gene activation was obtained for η i > 0.8.In Table 3, are identified the GRN Boolean attractors with those obtained with the continuous model.Gene activation/inhibition threshold parameter values are also indicated in Table 3.In the first column is listed each GRN gene whose expression leads to one of the three phenotypes that were identified as attractors in the dynamics of the set of Eqs (2), namely, (i) normal, (ii) precancer, and (iii) cancer.To analyze the multistability properties of the genetic states, variations and perturbations of the parameters were carried out.It was found that the stability depends on the activity of genes TP53(x 1 )-MDM2(x 6 ), TP53(x 1 )-P21(x 9 ) and ATR(x 2 )-BRCA1(x 4 ), that undergo positive feedback in the GRN.Stability changes occurred only in some cases as will be described in the following paragraphs.
Fig 5 shows a two dimensional (2D) projection of the phase space trajectories onto the plane x 1 -x 6 of the TP53(x 1 )-MDM2(x 6 ) positive feedback equilibrium states.To obtain this projection, in Fig 5A we have considered variations of the coordinates x 1 and x 6 and have kept constant the values of the remaining variables x i .This projection shows the points where the nullclines intersect identifying the steady states of the system obtained for parameters values fixed at y i = 0.2 and η i = 0.5 with i = 1, 6.The stability was determined by calculating the Jacobian and the corresponding eigenvalues.The squares represent unstable states while the circles correspond to stable states.Fig 5B shows the equilibrium points in a 2D projection of the phase space trajectories onto the plane x 1 -x 6 , and have kept constant the values of the remaining variables x i , and chose a set of parameter values that led to a normal phenotype.Similarly, Figs 5C and 5D show the equilibrium points in a 2D projection of the phase space trajectories onto the plane x 1 -x 6 , that are identified with precancer and cancer phenotypes, respectively.The values of microenvironmental concentrations and threshold values were held at y i = 0.2 and η i = 0.5, for i = 1, ...10.
Figs 6A-6D show the projections of the phase space trajectories onto the plane x 1 -x 6 to demonstrate how the stability of a normal cell genotype changes for different values of the concentrations of microenvironmental agents and the activation threshold.In Fig 6A the equilibrium states correspond to a normal genotype.For the concentration values of the microenvironmental agents η i = 0.25 with i = 1, ...10, the unstable and stable states that occurred at x 1 ≈ 0.25 and x 1 ≈ 0.35, disappeared, while the unstable state at x 1 = 1, remains invariant as shown in Fig 6B .On the contrary, when the concentration values of the microenvironmental agents double their value (y i = 0.4, with i = 1, ...10), the unstable state that occurred at low values of x 1 disappeared while the stable and unstable states at x 1 ≈ 0.35, and at x 1 = 1, remain invariant as shown in Fig 6C .As the activation threshold is increased up to a value η i = 0.75, the distance between the first unstable state and the stable state decreases, however, the unstable state at x 1 = 1 remains invariant again.These results demonstrate how the multistability of the system depends on the values of the parameters associated to the genetic and epigenetic characters of the system.
Figs  ) circuit.The analysis of the circuit stability led to similar findings as those described above.The bifurcation diagrams show one stable and one unstable states represented by two separate lines.This suggests that the TP53-MDM2 circuit contributes more to the GRN stability than the TP53-P21 circuit.Finally, in Fig 17, the ATR(x 2 )-BRCA1(x 4 ) circuit is analyzed.In this case there is only one stable state which represents a normal phenotype.Therefore, it can be concluded that regardless of the parameter values the stability of the circuit ATR(x 2 )-BRCA1(x 4 ) is robust.November 13, 2023       Attractor states a and c can be identified with normal phenotypes, whereas attractor states b, d, e, f, g, h, and i can be identified with a precancer phenotype.This precancer identification is based on the fact that some genes are either over-expressed or inhibited permanently in the Boolean GRN representation (see Figs 1 and 2).Finally, the state j is the only attractor that can be identified with a cancer phenotype.
The microarrray shown in Fig 19 demonstrates the effect of increasing estrogen threshold values on the expression of all GRN genes.The initial threshold values η oi , with i = 1, ...10 were increased in the following percentages: 10%, 25%, 50%, 75% and 100%, with respect to the values η oi = (0.3, 0.45, 0.15, 0.45, 0.15, 0.45, 0.3, 0.15, 0.3, 0.5).It is observed that the gradual increase of η oi changes de proportion of phenotypic states suggesting that the GRN is more adaptable to increasing values of η oi .The microarray presented in Fig 20 shows the effect of increasing the estrogen initial threshold value η o5 in gene HER2 as indicated in the figure caption.By comparing these microarrays it is found that the HERB2 initial threshold value η o5 favors the a precancer phenotype.By going from top to the bottom in the microarrays we note that the increase in η o5 reduces the proportion of precancer phenotypes and favors the cancer phenotype.These results are consistent with findings reported in [29] where it was found that the OE of HER2 leads to cancer phenotype and accelerates tumorigenesis as well as the occurrence of different Positive feedback equilibrium states for the TP53(x 1 )-P21(x 9 ) circuit.(A) 2D projection onto the plane x 1 − x 9 obtained by maintaining the other variables at a zero constant value.The red and blue lines represent the trajectories x 1 and x 9 , respectively.The unstable state is represented by the square and the stable state is represented by the circle.(B) Equilibrium points for the variables x 1 and x 9 when the other gene actions take constant values that led to a normal phenotype.(C) Equilibrium states for a precancer phenotype and (D) Equilibrium states for a cancer phenotype.In the four panels the values of microenvironmental concentrations and thresholds were held at y i = 0.2 and η i = 0.5, with i = 1, ...10.phenotypes.In addition, this analysis is in complete agreement with previous findings on the role of estrogen [19][20][21][22]30] and give support to the model and results presented in this paper.

Fig 2 .
Fig 2.Breast cancer attractors (Attr) considering Over-Expression (OE) and Knock-Out (KO) of specific genes.The upper table lists all the GRN attractors when gene TP53 is Over-Expressed (TP53-OE).On the other hand, the lower table lists the GRN attractors when genes ATM, MDM2 and AKT1 are either OE or KO.
Epigenetic mechanisms are essential for normal development and maintenance of tissue specific gene expression patterns.Disruption of these mechanisms may lead to altered cell functions.Because of this, epigenetic mechanisms play an important role in breast cancer development.It has been found that oxygen and estrogens distributed among November 13, 2023 4/39

Fig 3 .
Fig 3. Evolution of genetic states for different values of the microenvironmental agent threshold values (A) η i = 0 and (B) η i = 0.25, with i = 1, 2, ...10.It starts with a genetic state that corresponds to a normal phenotype and ends with the same state.

Fig 4 .
Fig 4. Genetic states evolution for different microenvironmental agent threshold values.(A) η i = 0.5 and (B) η i = 0.75 with i = 1, 2, ...10.It begins with a genetic state that corresponds to a normal phenotype and ends with precancer or cancer phenotype.We have considered that gene Over-Expression(OE) remains at value equal to one while Knock-Out (KO) stays at zero value.The sudden change in the genetic states behavior is originated from the OE and KO threshold values.
Fig 9  shows the stability of the T P 53(x 1 ) − M DM 2(x 6 ) circuit bifurcation diagram as a function of kinetic rate constants ν 1 and ν 6 for a normal phenotype.It was found that the stability range is short as indicated by the red lines.However, in the instability parameter region there may appear molecular processes such as inhibition of apoptosis or DNA damage that can change the cell phenotype.Fig10shows the ranges of Hill coefficients values γ i ≥ 2 for which the cooperative effect contributes to the system stability.Finally, as a preamble to future work, in Fig11the stability of the system is visualized as a function of the parameters α i , with i = 1, 6.These analysis show that asynchronous interactions, α i = ±1, may modify the stability of the system.In are shown the results of a similar analysis for the TP53(x 1 )-P21(x 9

Fig 5 .
Fig 5. Projection of the phase space trajectories onto the x 1 − x 6 plane for the TP53(x 1 )-MDM2(x 6 ) positive feedback equilibrium states for microenvironmental concentration values, y i = 0.2 and thresholds values η i = 0.5, with i = 1, ...10.The red and blue lines represent the nullclines of trajectories x 1 and x 6 , respectively.Unstable states are represented by squares and stable states by the circles.(A) Trajectories x 1 − x 6 obtained by maintaining the other variables at a zero constant value.(B) Phase space trajectories onto the plane x 1 − x 6 and have kept constant the values of the remaining variables x i , and chose a set of parameter values that led to a normal phenotype.(C) As in Fig (B) the trajectories correspond to a precancer phenotype, (D) Same as in previous panels but for a cancer phenotype.

Fig 6 .
Fig 6.Phase diagrams x 6 versus x 1 for the TP53(x 1 )-MDM2(x 6 ) circuit for a normal phenotype for different combination values of microenvironmental agent y i and activation threshold parameter η i .(A) y i = 0.2 and η i = 0.5, (B) y i = 0.1 and η i = 0.5, (C) y i = 0.4 and η i = 0.75, and (D) y i = 0.2 and η i = 0.75 with i = 1, ...10, in the four panels.Unstable states are represented with squares while the stable with circles.The red line represent the nullcline for variable x 1 and the blue line the nullcline for x 6 .

Fig 7 .
Fig 7. Multistability diagrams for the genetic states of the TP53(x 1 )-MDM2(x 6 ) circuit for a normal phenotype.Panels (A) and (B) show the genetic states x 1 and x 6 as a function of microenvironmental agent y 1 .Panels (C) and (D) show the genetic states x 1 and x 6 as a function of activation threshold parameter η 1 .The blue lines indicate unstable states while the red lines represent stable states.

Fig 8 .
Fig 8. Multistability diagrams for the genetic states of the TP53(x 1 )-MDM2(x 6 ) circuit for a normal phenotype, as a function of the β 1 parameter, panels (A) and (B), and β 6 parameter, panels (C) and (D).The blue lines indicate unstable states while the red lines represent stable states.

Fig 9 .
Fig 9. Multistability diagrams for the genetic states of the TP53(x 1 )-MDM2(x 6 ) circuit for a normal phenotype, as a function of kinetic rate constants ν 1 , panels (A) and (B) and ν 6 , panels (C) and (D).The blue lines indicate unstable states while the red lines represent stable states.

Fig 11 .
Fig 11.Multistability diagrams for the genetic states of the TP53(x 1 )-MDM2(x 6 ) circuit for a normal phenotype.The blue lines indicate unstable states while the red lines represent stable states.(A) x 1 versus the interaction rate α 1 .(B) x 6 versus α 1 .(C) x 1 versus α 6 and (D) x 6 versus α 6 .

Fig 12 .
Fig 12.  Positive feedback equilibrium states for the TP53(x 1 )-P21(x 9 ) circuit.(A) 2D projection onto the plane x 1 − x 9 obtained by maintaining the other variables at a zero constant value.The red and blue lines represent the trajectories x 1 and x 9 , respectively.The unstable state is represented by the square and the stable state is represented by the circle.(B) Equilibrium points for the variables x 1 and x 9 when the other gene actions take constant values that led to a normal phenotype.(C) Equilibrium states for a precancer phenotype and (D) Equilibrium states for a cancer phenotype.In the four panels the values of microenvironmental concentrations and thresholds were held at y i = 0.2 and η i = 0.5, with i = 1, ...10.

Fig 13 .
Fig 13.Multistability diagrams for the genetic states of the TP53(x 1 )-P21(x 9 ) circuit for a normal phenotype.The blue lines indicate unstable states while the red lines represent stable states.(A) x 1 versus microenvironmental agent y 1 .(B) x 9 versus microenvironmental agent y 1 .(C) x 1 versus activation threshold η 1 and (D) x 9 versus activation threshold η 1 .

Fig 17 .
Fig 17.Positive feedback equilibrium states for ATR(x 2 )-BRCA1(x 4 ) circuit.(A) 2D projection onto the plane x 2 − x 4 obtained by maintaining the other variables x i at a zero constant value.The red and blue lines represent the x 2 and x 4 trajectories, respectively.The only stable state is represented by the circle.(B) Equilibrium points for the variables x 2 − x 4 when the other gene actions take nonzero constant values and led to a normal phenotype.(C) Equilibrium states of a precancer phenotype and (D) Equilibrium states of a cancer phenotype.In the four panels the values of microenvironmental concentrations and thresholds were held at y i = 0.2 and η i = 0.5, with i = 1, ...10.

Fig 18 .
Fig 18. Microarrays showing gene expression and cell phenotypes.Normal phenotypes are represented by the states a and c.The states b, d, e, f, g, h, and i are identified with attractors in which some genes are Over-Expressed or Knocked-Out, suggesting a precancer phenotype.The state j is identified with the attractor that represents cancer phenotypes.The dotted line indicates an overlapping between states h and i.The initial threshold parameters for the set of ten genes were: η oi = (0.3, 0.45, 0.15, 0.45, 0.15, 0.45, 0.3, 0.15, 0.3, 0.5).

Fig 19 .Fig 20 .
Fig 19.Microarrays showing the effects of increasing estrogen threshold values in the expression of all ten GRN genes.From top to bottom, the increase was as follows: 10%, 25%, 50%, 75% and 100%, in comparison with the values assigned to obtain the microarrays shown in Fig 18.The parameters and color bars are the same as in Fig 18

Table 2 .
Gene expression parameter values for the continuous model.

Table 3 .
Summary of gene states for Boolean and states are indicated with blue lines.They were obtained for parameter values that led to a normal cell phenotype.In what follows it is shown that multistability maintains under small perturbations.Fig 7 shows four multistability diagrams x 1 − y 1 , x 1 − η 1 , x 6 − y 1 , and x 6 − η 1 , for the genetic states of the TP53(x1)-MDM2(x6) circuit with parameters that led to a normal phenotype.Fig 8 also shows four multistability diagrams x 1 − β 1 , x 1 − β 6 , x 6 − β 1 , and x 6 − β 1 for the genetic states of the TP53(x1)-MDM2(x6) circuit with parameters that led to a normal phenotype.As in previous figures the blue lines represent unstable states while the red lines correspond to stable states.These states correspond to the circles (stable state) and squares (unstable state), respectively, at the intersection of the nullclines shown in Fig 5B.

Table 4 .
Wild-Type.Attractors with one state.

Table 5 .
Wild-Type.Attractors with two states.

Table 10 .
ATR KO.Attractors with one state.

Table 11 .
ATR KO.Attractors with two states.

Table 12 .
ATR OE.Attractors with one state.

Table 13 .
ATR OE.Attractors with two states.

Table 14 .
ATM KO.Attractors with one state.

Table 15 .
ATM KO.Attractors with two states.

Table 16 .
ATM OE.Attractors with one state.

Table 17 .
ATM OE.Attractors with two states.