Boolean Network Model for Cancer Pathways: Predicting Carcinogenesis and Targeted Therapy Outcomes

Abstract

A Boolean dynamical system integrating the main signaling pathways involved in cancer is constructed based on the currently known protein-protein interaction network. This system exhibits stationary protein activation patterns – attractors – dependent on the cell's microenvironment. These dynamical attractors were determined through simulations and their stabilities against mutations were tested. In a higher hierarchical level, it was possible to group the network attractors into distinct cell phenotypes and determine driver mutations that promote phenotypic transitions. We find that driver nodes are not necessarily central in the network topology, but at least they are direct regulators of central components towards which converge or through which crosstalk distinct cancer signaling pathways. The predicted drivers are in agreement with those pointed out by diverse census of cancer genes recently performed for several human cancers. Furthermore, our results demonstrate that cell phenotypes can evolve towards full malignancy through distinct sequences of accumulated mutations. In particular, the network model supports routes of carcinogenesis known for some tumor types. Finally, the Boolean network model is employed to evaluate the outcome of molecularly targeted cancer therapies. The major find is that monotherapies were additive in their effects and that the association of targeted drugs is necessary for cancer eradication.

About the Authors:

Herman F. Fumiã

Affiliation: Departamento de Física, Universidade Federal de Viçosa, Viçosa, Minas Gerais, Brazil

Marcelo L. Martins

* E-mail: [email protected]

Affiliation: National Institute of Science and Technology for Complex Systems, Rio de Janeiro, Brazil

Introduction

Cancer is a genetic disease derived, with few exceptions, from mutations on single somatic cells that disregard the normal control of proliferation, invade adjacent normal tissues, and give rise to secondary tumors (metastasis) on sites different from its primary origin [1]. In the human population, cancer refers to more than forms of a disease that can develop in almost every tissue in the body [2]. Today, cancer replaced heart disease as the leading cause of death among the United States citizens younger than years [3] and will probably become the leading one in some other parts of the world within a few years [4]. Altough each cancer type has unique features, all these diverse tumors evolve according to a universal scheme of progression [5] which involves genetic and epigenetic events as well as an intricate network of interactions among cancer cells and their host microenvironment (stromal cells and extracellular matrix).

The tumor growth is intrinsically multiscale in nature. It involves phenomena occurring over a variety of spatial scales ranging from tissue (for instance, invasion and angiogenesis) to molecular length scales (for example, mutations and gene silencing), while the timescales vary from seconds for signaling to years for tumor doubling times. Moreover, all those processes are strongly coupled. Indeed, an oncogene activation may confer a proliferative advantage to a given cell, promoting its clonal expansion and the depletion of the nutrient and oxygen supply which, in turn, affect the growth of cell clones. To survive in a hypoxic (low level of oxygen) environment, the transformed cells may acquire new traits such as resistance to apoptosis by a tumor suppressor gene inactivation or activated synthesis of growth factors that stimulate angiogenesis. Thus information flows not only from the finer to coarser scales, but between any pair of scales [6].

Despite the extensive information on the genetic and molecular basis of cancer currently available, the integration of this information into the physiological environment of the functioning cell and tissue remains a major challenge. Due to the complexity and redundancy of tumor survival and growth pathways, increasing resistance and tumor progression still is the rule for patients with advanced cancers. Better diagnostic and effective anticancer therapies demand a fundamentally systemic understanding of the disease, starting from the molecular level. There, complexity emerges from the large number of interacting proteins and the cross-talk between diverse cell signaling pathways. This vast network of complexity, characterized by multiple feedback loops, will not be understood by merely describing all its component pathways. An integrative approach considering simplified cell-wide models of protein interactions dependent on external environmental signals and accumulated genetic alterations is demanded. In addition, modeling protein interaction networks in cancer is essential to construct the “microscopic” (molecular) level in multiscale models of tumor growth [6].

Concerning the modeling of genetic interactions, Boolean networks is a promising framework [7]. Instead of providing quantitatively precise dynamical trajectories taken by complex networks, this class of discrete systems, with binary states, qualitatively predict the sequences of states accessed by these networks along their temporal evolution. This is especially true for signaling and regulatory circuits where activation and inhibition are the basic responses. Furthermore, Boolean models use much less parameters, such as biochemical reaction rates or binding affinities, often hard to measure, than do traditional differential equations approach. Successful applications of Boolean networks in biology include the reproduction of the yeast S. Pombe's cell-cycle [8], the mammalian cell cycle [9], the course of cell differentiation in early embryogenesis [10], the signaling mechanisms underlying T cell activation [11], and the behavior of the apoptotic pathway [12], [13].

In this study, we construct a Boolean network model integrating the main signaling pathways involved in cancer. These pathways and the network interconnecting them are discussed in the next subsection. Established the network topology, its dynamics is defined and validated in subsequent subsections. Then, the dynamical attractors, their stability to mutations, and the network's response to targeted “attacks” are reported in the following section. These results are interpreted in terms of the mutational events leading to carcinogenesis and cancer cell response to molecularly targeted agents. Finally, we conclude with a discussion where these biological interpretation is emphasized and confronted with experimental data on cancer genome and oncogenesis.

The Main Cancer Pathways

In a seminal paper [14], Hanahan and Weinberg proposed a logical framework for understand the diversity and complexity of cancer. The key notion is that along the multistep process of tumor pathogenesis, normall cells need to acquire six biological traits in order to ultimately become malignant. These are the capabilities of “sustaining proliferative signaling, evading growth suppressors, resisting cell death, enabling replicative immortality, inducing angiogenesis, and activating invasion and metastasis” [15]. The acquisition of such traits is ensured by genome instability. All genes are potentially subjected to mutations. Thereby, there are several alternative ways to achieve the same cell phenotypic transformation. But there are many fewer signaling pathways controlling cell response than genes. Rather than individual genes, it seems more appropriate focusing on pathways that have a role in many tumors [16]. They include those involving receptor tyrosine kinase (RTKs), phosphatidylinosital 3-kinase (PI3K)/AKT, WNT/-Catenin, transforming growth factor- (TGF-)/Smads, retinoblastoma protein (Rb), hypoxia-inducible transcription factor (HIF-1), p53 and ataxia-telangiectasia mutated (ATM)/ataxia-telangiectasia and Rad3-related (ATR) protein kinases. These major pathways regulating cell death and proliferation share some genes and exhibit a substantial cross-talk among them. Some pathways were detailed studied in isolation [12], [17], [18], but in order to integrate them a natural organizing principle is represent these pathways as a network. Given the complexity and lacunas present in its structure, an operational alternative is work with a simplified model network.

The constructed, simplified signaling network is illustrated in Figure 1. It contains nodes and edges. These nodes represent a significant subset of proteins involved in cancer and the network edges produce numerous parallel pathways and alternative routes through which transformed cells can sustain aberrant gene expression patterns, survive and develop further malignancy. The network has input nodes for applying distinct environmental stimuli and stresses to the cell such as hypoxia, carcinogens, nutrients depletion, proliferative and growth suppressive signalings.

[Figure omitted. See PDF.]

Figure 1. Simplified cancer network.

Links correspond to interactions between proteins and to each node is associated a threshold function, eq. 1. Activating interactions are represented by arrows and inhibiting interactions by lines with a bar. The input nodes are shown in red.

https://doi.org/10.1371/journal.pone.0069008.g001

The topological structure of this cancer network was characterized in terms of its shortest path length, clustering coefficient and connectivity or degree distributions [19]. These quantities were compared with their average counterparts for random networks with the same number of nodes (see Table 1). It can be noticed that the cancer network has a much higher clustering coefficient than random networks. Major features of complex networks are their connectivity distributions. Directed networks are characterized by in-degree and out-degree distributions. is the probability that a node in the network has inputs or is “regulated” by other vertices. In turn, is the probability that a node “regulates” other vertices. In Figure 2 these distributions for the cancer network are shown.

[Figure omitted. See PDF.]

Figure 2. Connectivity distributions for the cancer network.

(a) In-degree and (b) out-degree distributions. The insets suggest exponential and power law distributions for the number of nodes regulating and regulated by a given node, respectively.

https://doi.org/10.1371/journal.pone.0069008.g002

[Figure omitted. See PDF.]

Table 1. Topological properties of the cancer network and their average counterpats for an ensemble of random networks.

https://doi.org/10.1371/journal.pone.0069008.t001

A Boolean Dynamics for Cancer Pathways

Each protein i , a node in the network, is represented by a binary state , . When the protein is functionally active. On the contrary, when the protein is inactive. The network state at a given time is specified by its protein activity pattern . Each protein interacts with other input proteins. These inputs are all nodes in the network from which a directed link is sent towards the protein , including eventually itself. So, is the in-degree of the node . In turn, each link can be either a activating or inactivating interaction. Activation or inhibition can be the result of distinct biochemical mechanisms such as transcriptional regulation, phosphorylation, enzymatic or biding interactions.

The dynamics of the network proceeds in discrete time steps through the simultaneous (parallel) update of the states of its nodes according to the rule.(1)

Here, is the interaction strength from input on protein . An activation interaction is positive and an inhibitory one is negative. The threshold function is the unitary step function ( if but if ). Finally, is the activation threshold of protein . Thus, every protein evaluates the present stimulus from all its inputs. If the overall stimulus it receives at time overcomes its activation threshold, the protein activates, or stays active if it was already active; otherwise, it turns inactive or stays inactive. Text S1 lists all the update rules of the model (see Text S1). As emphasized in reference [8], a dramatic simplification inserted in this evolution rule consists in neglect any difference in the time scales of the biochemical interactions involved.

Results

The Basic Phenotypes

Since the state space of a Boolean network with nodes contains different configurations, its deterministic dynamics, viewed as a flow in this state space, ultimately will drive the system towards attractors. Such attractors are particular subsets of states, either a fixed point, i. e, a single network configuration, or a limit cycle of period , comprised by states sequentially visited by the network dynamics. These attractors correspond to specific protein activation patterns and can be interpreted as distinct cell phenotypes.

Thirty two million initial states, associated to all environmental conditions, flowed into attractors ( fixed points and limit cycles). All attractors are listed in the Text S1. The state space is hierarchically organized. At a higher level, it is partitioned into subsets of states by distinct environmental conditions. At a lower level, every of these subsets is subdivided into basins of attraction associated to distinct attractors. None of these attractors can be reached starting from initially distinct environmental conditions. Thus, the repertoire of cell behaviors (attractors) is determined univocally by the cell microenvironment.

Although distinct, these attractors can be classified in groups characterized by specific phenotypes, some of them comprised of very similar elements, i. e, having small Hamming distances among them (Text S1). The phenotypes were defined taking into account the states of a small subset of nodes, instead of all proteins on the network. Considering the effects of mutations, reported on the next subsection, they include the following basic cell phenotypes: apoptotic, characterized by active caspases; glycolytic, with H1f1 activated under normoxia; immortalized, in which hTert is active; migratory, associated to inactivate E-cadherin; mutator, corresponding to inactive Atm/Atr proteins in the presence of DNA damage; proliferative, in which cyclins are activated along the cell cycle in the correct sequence; and quiescent, with cyclins inactive or activated in a wrong sequence. In terms of such phenotypes, the network response to the diversity of microenvironments is highly constrained, as illustrated in Figure 3. Among the attractors, correspond to apoptotic, to proliferative and to quiescent phenotypes, which attract , , and of tested initial states, respectively.

[Figure omitted. See PDF.]

Figure 3. Network's responses to distinct environmental conditions.

Three phenotypes (apoptotic, proliferative and quiescent) are generated in response to all distinct environmental conditions. Here, a microenvironment is specified by the binary sequence of values associated to input nodes (carcinogens, growth factors, nutrient supply, growth suppressors, hypoxia). For instance, the microenvironment (11000) corresponding to a carcinogenic and mitogenic background leads the cell to either an apoptotic (in 99.8% of the initial states) or a quiescent phenotype (rarely, ). In our network, carcinogens elicit DNA damage and TNF- is the suppressive growth signal.

https://doi.org/10.1371/journal.pone.0069008.g003

In particular, for some environmental conditions the behaviors emergent from the whole network were compared with those of a normal cell. Under normoxia and adequate nutrient supply, the network reaches from all compatible initial states a fixed point interpreted as a quiescent phenotype. If, in addition to normoxia and nutrient abundance, growth factor signaling is receipt, the network always evolves to a limit cycle associated to a proliferative phenotype. So, as expected, normal cells are totally dependent for their proliferation upon mitogenic signals [20]. Furthermore, under hypoxia, adequate nutrient supply and absence of growth factors, the network is attracted from all compatible initial states to a fixed point corresponding to an apoptotic phenotype. Consistently, programmed cell death is the typical response of somatic cells to many forms of stress such as hypoxia and nutrient deprivation [20].

Under a hypoxic environment with nutrients and growth factors, the network exhibits bistability. It is either attracted to the quiescent phenotype (in of the initial states) or to the apoptotic phenotype (in of the initial states). Accordingly, it is known that hypoxia-dependent activation of HIF1 inhibits Myc, leading to cell cycle arrest [21]. Also, HIF1 can bind to and stabilize p53, resulting in apoptosis or growth arrest. Bistability is also observed if DNA damage is introduced in a scenario of normoxia, adequate nutrient supply, and mitogenic signaling: around of the compatible initial states are attracted to the apoptotic phenotype, while a very small fraction () of them reach the proliferative phenotype. Again, it is widely known that the cell cycle is easily interrupted and apoptosis triggered by DNA damage in normal cells. But a proliferative response, although rare, endows altered somatic cells with a proliferative capacity. In our model, this proliferative response is associated to two distinct limit cycles. In one of them the anti-apoptotic signals – Bcl2, Bcl-Xl, and Mdm2 – are consistently active whereas pro-apoptotic signals – Bad, Bax, p53 – are inactive. In the other, Bad and Bax are inactive but p53 and anti-apoptotic signals oscillate in such a way that whenever p53 is activated, the same occurs with Bcl2, Bcl-Xl, and Mdm2. Then, the result is that caspases, the effectors of apoptosis, are always inactive along these limit cycles.

Summarizing, the whole network generates responses coherent with those observed in a normal cell under different somatic environments, indicating the fundamental correctness of the model.

Mutational Events and Carcinogenesis

We additionally checked the robustness (stability) of attractors to mutations in network nodes and/or links. This is a central feature because incipient cancer cells need to acquire hallmark traits to ultimately become malignant [15] and genome instability underlies these acquisitions. Once a mutation was introduced, the node DNA damage is permanently turned on, activating the Atm/Atr pathway. We focused on the attractors associated to two environmental conditions, namely, adequate nutrient supply and either normoxia or hypoxia, frequently present in early carcinogenesis.

Under normoxia and adequate nutrient supply, it was found that mutations in proteins transform the formerly quiescent, normal cell into a proliferating one. These proteins, as well as the nature of the driver mutations and their efficacy are listed in Table 2. In turn, under hypoxia and adequate nutrient supply, mutations in proteins enable the transformed cell to evade apoptosis formerly induced by hypoxia (Table 3). The protein Nf-B is common to Tables 2 and 3, hence it can enable a transformed cell to simultaneously acquire proliferative capacity and evading apoptosis.

[Figure omitted. See PDF.]

Table 2. Driver mutations under normoxia.

https://doi.org/10.1371/journal.pone.0069008.t002

[Figure omitted. See PDF.]

Table 3. Driver mutations under hypoxia.

https://doi.org/10.1371/journal.pone.0069008.t003

We also investigated the effect of defective DNA integrity sensors that impair a cell to detect the occurrence of mutations. Now, the node DNA damage is permanently turned off and does not activates the Atm/Atr pathway. Under normoxia and adequate nutrient supply, the number of mutated proteins that transform a quiescent, defective cell into a proliferating one increases to . Yet, the number of mutated proteins that confer to a hypoxic, defective cell the capacity to evade apoptosis also increases to . These proteins include Akt, Bcl2, Egfr, Nf-B, p53, Pi3k, Pten, Ras, and Wnt (see Tables S3 and S4). In the absence of an intact DNA damage repair pathway, in which Atm and Atr play central roles, our results indicate that network attractors become more prone to structural changes or, in biological terms, exhibit increased genomic instability.

Finally, we investigated if nodes whose mutations can confer hallmark capabilities to the transformed cell have special status in network topology. For all them their betweenness centrality were determined [19]. Three groups have been observed. The first one, comprised of eight nodes (Akt, Hif1, hTert, Ikk, mTor, Myc, Nf-B, and p53), has more than twice the network average centrality . Further, the group average connectivity is and seven of its elements are highly connected (). These results indicate the centrality of the nodes in this group for which converge or through which crosstalk distinct signaling pathways. The second group, containing nodes with , includes Mdm2 and Pdk1. In addition, the group average connectivity is and half of them exhibits intermediate connectivities (). At last, the third group is characterized by small () and connectivities (). Among its elements, Bcl2,Bcl-x, Egfr, Gli, Nf1, Phd, Pi3k, Pkc, Pten, Ras, Snail, Vhl, and Wnt, nine have small connectivities . Even though these nodes are not topologically central, almost all are nearest neighbors (direct regulators) of central nodes from the first group. Hence, they assume major roles in network dynamics.

Colorectal carcinogenesis.

Here, we investigate if cancer cells need only a few driver mutations (those that change phenotypic traits) to deal with all environmental constraints and advance towards a fully malignant phenotype. As a paradigm, the carcinogenesis of the colorectal cancer was considered. So, the first mutations introduced in the network were Apc deletion and Tcf interactions with their targets strengthened by a factor . These mutations lead to a structural instability of the network “phase portrait”. Now, there are attractors, apoptotic, proliferatives, and quiescent, which attract , and of tested initial states, respectively. The number of proliferative attractors and the sizes of their basins increased at the expenses of the quiescent attractors. However, no anti-apoptotic advantage was observed. As expected, the network response was environmental dependent. Under normoxia and nutrient availability, these mutations lead to a proliferative phenotype. However, this proliferative advantage is lost under hypoxia or genotoxicity (DNA damage), when apoptotic phenotypes are observed.

Next in the sequence, a new mutation – Ras constitutively activated – was implemented. As a result the network dynamics exhibits attractors, apoptotic, proliferatives, and quiescent, which attract essentially the same fractions of tested initial states (, and , respectively) as before (Apc and Tcf mutated). However, this additional mutation conferred to the network a small chance to overcome the hypoxic barrier, sustaining quiescent phenotypes for of tested initial states under a hypoxic, but nutrient rich environment. Further, constitutive Ras activation also can lead to a proliferative phenotype in of initial states in a normoxic, nutrient rich, but genotoxic environment.

In the sequence, Smad4 was deleted. This mutation increases to the number of attractors, apoptotic, proliferative, and quiescent. The network acquires proliferative phenotypes for all initial states under normoxia and adequate nutrient supply, even receiving inhibitory growth signaling provided by a constitutively active Tgf-. It is worthy to mention that the network exhibits only apoptotic and quiescent attractors if Tgf- is constitutively active and Smad4 is functional (undeleted). So, Smad4 deletion in accumulation with the previous mutations endows the transformed network with the capacity to evade suppressive growth signals.

The next mutations were Pten deletion and doubling of Akt interaction strengths. In consequence, the number of attractors further increases to , apoptotic, proliferative, and quiescent, which attract , and of tested initial states, respectively. Again, different microenvironments elicit distinct responses. Under normoxia and adequate nutrient supply the network always exhibit aggressive (proliferative, glycolitic and immortalized) phenotypes. But if hypoxia replaces normoxia, in addition to proliferative, glycolitic and immortalized phenotypes which attract of the initial states, there are quiescent attractors toward which of initial states converge. Adding growth suppressors or DNA damage to the former microenvironment can at most lead to quiescence. For instance, in normoxic, nutrient rich and genotoxic microenvironment, of initial states are driven to proliferative, glycolytic and immortalized attractors, whereas of them are driven to quiescent attractors. Therefore, since hypoxia or functional DNA damage sensors can lead to quiescent phenotypes, some constraints persist impairing tumor growth.

The last mutation was p53 deletion. Its result is decrease to the number of attractors, apoptotic and proliferative, both attracting of the initial states. Indeed, apoptosis for of the initial states is the minimum value possible because in our network active TNF- leads to p53-independent activation of caspases. Nevertheless, the main result is that the network always exhibits proliferative, glycolytic and immortalized phenotypes in microenvironments with adequate nutrient supply, hypoxic or normoxic, even genotoxic, which activate DNA damage sensors, and under growth suppressor signaling. Almost all barriers to tumor growth were overcome after this sequence of few mutations.

In summary, as shown in Figure 4, our simulations reveal that each driver mutation in the canonical route for the colorectal cancer [22] contributes to increase either the proliferative capacity or the resistance to apoptosis of the transformed cell. In particular, although Smad4 is mutated in only of colorectal cancers, its mutation in concert with the others in the classical colorectal carcinogenesis model generates more aggressive tumor cells. Indeed, their associated proliferative phenotypes attract of the initial states against only in the absence of the Smad4 mutation. Further, the model indicates that other mutations outside this classical route of colorectal carcinogenesis also leads to proliferative and apoptotic resistant phenotypes. These are the cases, for instance, of Pten, or p53, or Atm, or Fadd, or Chk deletions after Apc and Ras mutations. Alternatively, the constitutive activation of Pi3k, or Akt, or Bcl2, or Mdm2 again after Apc and Ras mutations decreases apoptosis and increases proliferation.

[Figure omitted. See PDF.]

Figure 4. Network response to driver mutations in colorectal carcinogenesis.

Fraction of initial states evolving into apoptotic, proliferative or quiescent attractors (phenotypes) for all environmental conditions after the sequential accumulation of each driver mutation in colorectal cancer.

https://doi.org/10.1371/journal.pone.0069008.g004

The Outcomes of Targeted Therapies

The rationale of targeted therapy is inhibit critical, functional nodes in the oncogenic network to elicit the cessation of the tumorigenic state through apoptosis, necrosis, senescence, or differentiation [23]. We performed a survey of nodes in our Boolean model whose inhibition or activation (reintroduction of wild-type proteins) either increase the basins of attraction of apoptotic and quiescent phenotypes or decrease those associated to proliferative phenotypes. Specifically, as a model for fully developed colorectal cancer cells, a network carrying mutations in Apc, Ras, Smad4, Pten, and p53, was considered. For this network, only two phenotypes were observed under all microenvironmental conditions and from every initial state tested: apoptotic () and proliferative (). All nodes were individually targeted. Our simulations reveal that the inhibition of Pdk1, Akt, E2F, Cyclins D and E, and Mdm2 enhances quiescent phenotypes ( for Pdk1 and Akt, and for E2f, cyclins and Mdm2) and also impair the proliferative ones ( for Pdk1 and Akt, and for E2f, cyclins and Mdm2). The same is observed with the activation of p53, Rb, Cdh1, p21, p27, and Smad4 (quiescence increases to and proliferation vanishes for p53, Rb, Cdh1, p21 and p27, whereas for Smad4 both fractions changes to ). In turn, the activation of FADD and Caspases enhances apoptosis to of chance. Paradoxically, the model predicts that either inhibition or activation of Cyclins A and B as well as Cdc20 enhances quiescence and impair proliferation ( and , respectively).

Concerning the combinatorial application of drugs currently under clinical trial [23], the effects of several targeted inhibitors on our colorectal cancer cell model were investigated under three environmental conditions, namely, normoxia and adequate nutrient supply, hypoxia and adequate nutrient supply, and normoxia, adequate nutrient supply and carcinogenic stress. For these conditions our colorectal cancer cell model is always attracted to a proliferative, immortalized and glycolytic phenotype. According our simulations the combined use of mTor and hTert inhibitors reverses two malignant hallmarks: glycolysis and immortalization; the association of mTor and Cdks inhibitors prevents glycolysis and cell cycling; and the combination of Cdks and hTert inhibitors prevents proliferation and immortalization. So, since in our model monotherapies with mTor, Cdks, and hTert inhibitors prevents glycolysis, cell division, and immortalization, respectively, the therapeutic outcome seems to be additive or “linear” (i. e., obeys a superposition principle). In turn, the association of p53 reactivators and Bcl-2/Bcl-Xl inhibitors leads to cell cycle arrest under hypoxia and adequate nutrient supply or normoxia, adequate nutrient supply and carcinogenic stress, but has no effect under normoxic and nutrient rich microenvironments. Another association involving p53 reactivators and Vegf inhibitors blocks cell cycling in of the initial states tested in these three environmental conditions. In our model, the action of Vegf inhibitors is simulated by forcing hypoxia and nutrient depletion on the input nodes.

Finally, we tested mono- and two-drugs therapies using mTor, Vegf, Ras, and Pkc inhibitors in cells at distinct stages of malignancy. Again, our colorectal cancer cell model was considered. Our results reveal that a full malignant cell (network carrying mutations in Apc, Ras, Smad4, Pten, and p53) at most lost its glycolytic phenotype if mTor is inhibited by a therapeutic agent, but its limitless proliferative capacity remains unchanged. In contrast, a cell in the previous carcinogenic stage (network with mutations in Apc, Ras, Smad4, and Pten) can have its chance of proliferation decreased under both monotherapy and two-drugs association. Indeed, the fraction of proliferative phenotypes decreases from to about , , and using Vegf, mTor Vegf, and Pkc Vegf inhibitors, respectively. Also, the therapeutic outcome depends on the external environment. For instance, combined mTor and Pkc inhibition reduces the proliferation chance from about to under normoxia, adequate nutrient supply and carcinogenic stress, but this change is significantly smaller under hypoxia and adequate nutrient supply, from about to . So, these results demonstrate that each treatment distinctly affects cells in different grades of malignancy and eventually clones will emerge, rendering the therapy ineffective.

Discussion

We constructed a Boolean dynamical system integrating the main cancer signaling pathways in a simplified network. The dynamics of this network is controlled by attractors associated to apoptotic, proliferative and quiescent phenotypes that qualitatively reproduce the behaviors of a normal cell under diverse microenvironmental conditions. Indeed, the network response is highly constrained with , , and of the initial states attracted to apoptotic, proliferative and quiescent phenotypes, respectively. So, under persistent stress, apoptosis or cell cycle arrest are the rule. Further, cell proliferation is tightly regulated, occurring almost only in a normoxic environment and in the presence of growth signaling. As observed in our model, GF signaling significantly increases the stability of the surviving (proliferative and quiescent) phenotypes while inhibits apoptosis. This result is consistent with the findings of Mai and Lieu [13] that, using a Boolean network integrating both the intrinsic and extrinsic pro-apoptotic pathways with pro-survival GF signaling, demonstrated that apoptosis can be induced either easily or difficultly depending on the balance between the strengths of pro-apoptotic and pro-surviving signals.

Our simulational results demonstrate that perturbations in some network nodes elicit phenotypic transitions. We interpreted them as driver mutations and can represent either the constitutive activation or inactivation of a node or yet an increase in the interaction strengths of a node with its targets. Under normoxia and adequate nutrient supply, we found that mutations in Egfr, Gli, Nf1, Nf-B, Pi3k, Pkc, Pten, Ras, and Wnt transform the formerly quiescent, normal cell into a proliferating one. The resultant clonal expansion often leads to hypoxia. Additional mutations in Akt, Bcl2, Bcl-Xl, Ikk, Nf-B, p53 and Snail enable the transformed cell to evade apoptosis formerly induced by hypoxia. These driver mutations predict by our model are included among the approximately of genes in the human genome causally implicated in tumor progression by diverse census of cancer genes recently performed [24], [25], [26]. The predicted drivers clusters on certain signaling pathways as, for instance, in the classical Mapk/Erk (Egfr, Nf1 and Ras), Pi3k (Pi3k, Pkc, Pten, Akt), p53 and Wnt signaling pathways. Also, sequencing data reveal that some of them are significantly mutated in cancers: Pi3k, Pten, and Akt in breast cancer [26], [27]; Ras and p53 in either breast and colorectal cancers [26], [28]; p53 and Nf1 in ovarian carcinoma [29]; p53 and Pten in small-cell lung cancer [30]; and Egfr, p53, Nf1, and Pi3k in human glioblastoma [31]. In turn, Gli is target of translocation or amplification in either breast and colorectal cancers [28], Bcl2 is amplified in squamous cell lung cancers [32] and the Wnt signaling pathway is altered in of all colon and rectal cancer [33]. The other predicted drivers realize a main result obtained from our model: driver nodes are not necessarily central in the network topology, but at least they are direct regulators of central components towards which converge or through which crosstalk distinct cancer signaling pathways.

Our Boolean model for cancer pathways allowed us to explore carcinogenesis at the molecular level. Carcinogenesis is an evolutionary process driven by a sequential acquisition of stepwise, somatic-cell mutations with concomitant subclonal selection [34]. Cells in a nascent tumor are continuously facing environmental stresses often in response to inflammation. Inflammation elicits cytokine-induced hyperplasia, genotoxicity, and ROS-induced cell death. In hyperplasic epithelia, abnormal cell growth leads first to hypoxia and selection for a glycolytic phenotype, resulting in increased acidity and nutrient and growth factor deprivation. Severe chronic hypoxia can select for apoptosis resistance or mutated p53, further promoting the accumulation of mutations. Low pH can generate DNA damage and glucose deprivation strongly reinforces the selection for activated oncogenes. The spatial heterogeneity of microenvironments within a primary tumor selects for cell migration, which may produce invasion and metastasis. So, as was demonstrated by our simulations of the colorectal carcinogenesis, in order to evolve, cancer cells perform sequential and random searches for phenotypic solutions to overcome the barriers imposed by their altering environment. Ultimately, each driver mutation in a carcinogenic route contributes to increase either the proliferative capacity or the resistance to apoptosis of the transformed cell, while several “passenger” mutations can accumulate along the process. Our Boolean model supports distinct carcinogenic routes characterized by specific sets of critical mutations embedded within varied spectra of passenger mutations. Thus, as observed in genomic analysis of diverse cancers, different combinations of mutated genes exist and most samples of a given tumor type differ from all others [26].

The development of effective targeted therapies copes with a major challenge. From the experimental view point it is very hard to determine which proteins have critical roles in tumor progression to be targeted pharmacologically. This is where mathematical models integrating diverse cancer pathways can be helpfull. Using our model, we find that the inhibition of some nodes (Pdk1, Akt, E2F, Cyclins D and E, and Mdm2) or the activation of others (p53, Rb, Cdh1, p21, p27, and Smad4) enhances quiescent phenotypes and also impairs the proliferative ones. However, all the monotherapies tested were ineffective to simultaneously reverse all the malignant hallmarks and seems to be additive or “linear” in their effects. Thus, our simulations indicate that the association of targeted drugs is necessary for cancer control or eradication. Further, we find that each treatment distinctly affects cells in different grades of malignancy or subject to distinct microenvironmental conditions. Although a combined targeted therapy can eradicate a set of cancer cells with certain phenotypes, eventually clones will emerge, rendering the therapy ineffective. So, these results provide additional support to the view that a combinatorial series of drugs applied concurrently to block cancer pathways and alter the tumor environment is needed to eradicate all of the cancer cells in a tumor [23]. However, in contrast to the present model focusing an isolated cell, tumors are highly heterogeneous systems on larger length scales patterned in patches (local microenvironments) that continuously change in space and time. Thus, a single drug association can have a limited effectiveness or to be efficacious demands a large number of therapeutic agents which imposes unacceptable clinical risks.

Finally, we briefly comment on the limitations and prospects of our modeling approach. In its present form, the cancer network and its Boolean dynamics are able to reproduce some biologically relevant features of carcinogenesis. It will certainly be true that further improvements on the topology of the cancer network will lead to better results and predictions. Indeed, adding new nodes, links and signaling pathways to the network can generate new convergences and redundancies essential to achieve more accurate predictions concerning cancer gene candidates and effective targets for therapies. Despite these rather straightforward secondary advances, the primary limitation of the model lies in its own single-scale nature constrained to the molecular level of a cell. Ultimately, cell-cell interactions in a changing environment determine which of the possible transformed cell clones will be generated and selected along cancer progression and, in consequence, the efficacy and safety of molecularly targeted therapies. Hence, Boolean network models for cancer pathways might be suitable to describe the microscopic scale in multiscale models of cancer growth and therapy [6]. Currently, we are pursuing this goal by extending our previous multiscale models considered in references [35] and [36].

Materials and Methods

The network integrating the main signaling pathways involved in cancer was constructed based on the current literature and protein-protein interactions map reported on KEGG database [1], [16], [37]. Specifically, subgraphs of the PI3K-AKT, mTOR, MAPK, HIF1, TGF-, WNT, NF-B, TNF, cell cycle, p53, and apoptosis KEGG pathways were included in our network model. Detailed descriptions of these pathways can be found in several comprehensive reviews [38], [39], [40], [18], [41], [42], [43], [44], [45], [46]. Its input nodes, i. e., those which are not regulated by other nodes and therefore associated to environmental cues, represent growth factors, nutrient and oxygen supplies, carcinogens (or mutagens), and apoptotic signal (Tnf-). This simplified network focuses attention on those gene products and signaling pathways that are implicated more universally in cancer and which seems to be sufficient to generate a simple model still capable of explaining some general features of carcinogenesis.

Almost all activating interactions have a strength , whereas the inhibitory ones have . The exceptions are listed in Table 4 and totalizes three of the network links. Similarly, the majority of the activation thresholds has , but there are nodes among the proteins with different values. They are also listed in Table 4. The assignment of each threshold value was guided by the local network topology (the inputs of each node) and very general dynamical arguments or expected biological responses of normal cells (see comments in Table 4). For instance, some nodes are “regulated” only by inhibitory inputs. For them, in order to ensure their activation and the spreading of their signals on the network if those inputs are inactive. Other nodes are binary protein complexes whose activation thresholds were set to because their two constituents must be activated to form the complexes.

[Figure omitted. See PDF.]

Table 4. Interaction strengths and activation thresholds with special values.

https://doi.org/10.1371/journal.pone.0069008.t004

In our simulations, every environmental condition, each one corresponding to a fixed set of input nodes, were analyzed individually. Fixed the input nodes, initial states for the other nodes were randomly chosen according the following protocol. Firstly, distinct and equally spaced values for the probability of a node to be active, , , , , were used. Secondly, configurations for the non-input nodes were randomly generated for each value of .

In order to analyze the stability of the network attractors, two types of mutations were implemented: constitutive activation or inactivation of a chosen protein ( or , respectively), and protein overexpression (constitutive activation with concomitant increasing in its interaction strengths with its target proteins). In the simulations, all nodes were subjected to such mutations and their effects were observed for randomly chosen initial states.

The simulations were performed in gfortran using a microcomputer with a Intel i3-2310 M, 2.1 GHz dual core processor and 4 Gb RAM. The network was drawn with Cytoscape [47].

Supporting Information

[Figure omitted. See PDF.]

Figure S1.

Attractors observed under distinct microenvironmental conditions. These microenvironments are described by their binary codes listed at the left margin of the figure. Specifically, (normoxic and nutrient rich), (normoxic and plenty of nutrients and growth factors), (hypoxic and nutrient rich), (hypoxic, plenty of nutrients and growth factors), and (normoxic, plenty of nutrients and growth factors, and under genotoxic stress – carcinogens). These conditions frequently occur either during different stages of carcinogenesis or in certain regions within spatially heterogeneous solid tumors.

https://doi.org/10.1371/journal.pone.0069008.s001

(PDF)

Table S1.

Fixed points of the cancer network. Fixed points for all 32 possible environmental conditions.The basin size of a fixed point was estimated as the fraction of initial states driven to that attractor.

https://doi.org/10.1371/journal.pone.0069008.s002

(PDF)

Table S2.

Limit cycles of the cancer network. Limit cycles for all 32 possible environmental conditions. The basin size of a cycle was estimated as the fraction of initial states driven to that attractor.

https://doi.org/10.1371/journal.pone.0069008.s003

(PDF)

Table S3.

Driver mutations under normoxia. New driver mutations under normoxia and adequate nutrient supply in the context of defective DNA integrity sensors.

https://doi.org/10.1371/journal.pone.0069008.s004

(PDF)

Table S4.

Driver mutations under hypoxia. New driver mutations under hypoxia and adequate nutrient supply in the context of defective DNA integrity sensors.

https://doi.org/10.1371/journal.pone.0069008.s005

(PDF)

Text S1.

Supporting information text file.

https://doi.org/10.1371/journal.pone.0069008.s006

(PDF)

Acknowledgments

We would like to thanks our colleagues S. C. Ferreira and L. R. Paiva for helpful discussions and comments on the manuscript.

Author Contributions

Conceived and designed the experiments: HFF MLM. Performed the experiments: HFF. Analyzed the data: HFF MLM. Contributed reagents/materials/analysis tools: HFF MLM. Wrote the paper: MLM.

Citation: Fumiã HF, Martins ML (2013) Boolean Network Model for Cancer Pathways: Predicting Carcinogenesis and Targeted Therapy Outcomes. PLoS ONE 8(7): e69008. https://doi.org/10.1371/journal.pone.0069008

References

1. Weinberg R (2007) The biology of cancer. New York: Garland Science.

2. Kumar V, Abbas AK, Aster JC, Fausto N (2010) Robbins and Cotran pathologic basis of disease. 8th ed. Philadelphia: W.B. Saunders Co.

3. Siegel R, Ward E, Brawley O, Jermal A (2011) Cancer Statistics, 2011: the impact of eliminating socioeconomic and racial disparities on premature cancer deaths. CA Cancer J. Clin. 61: 212–236.

4. Murray CJL, Lopez AD (1997) Alternative projections of mortality and disability by cause, 1920–2020: Global burden of disease study. Lancet 349: 1448–1504.

5. Clark WH (1991) Tumour progression and the nature of cancer. Br. J. Cancer 64: 631–644.

6. Martins ML, Ferreira Jr SC, Vilela MJ (2009) Multiscale model for biological systems. Curr. Opin. Colloid Interface Sci. 15: 18–23 (2010)..

7. Kauffman SA (1993) The origins of order: Self-organization and selection in evolution. New York: Oxford University Press.

8. Davidich MI, Bornholdt S (2008) Boolean network model predicts cell cycle sequence of fission yeast. PLoS One 3(2): e1672.

9. Fauré A, Naldi A, Chaouiya C, Thieffry D (2006) Dynamical analysis of a generic Boolean model for the control of the mammalian cell cycle. Bioinformatics 22(14): e124.

10. Silva HS, Martins ML (2003) A cellular automata model for cell differentiation. Phys. A 322: 555–566.

11. Saez-Rodrigues J, Simeoni L, Lindquist JA, Hemenway R, Bommhardt U, et al. (2007) A logical model provides insights into T cell receptor signaling. PLoS Comput. Biol. 3(8): e163.

12. Schlatter R, Schmich K, Vizcarra IA, Scheurich P, Sauter T, et al. (2009) ON/OFF and beyond - A Boolean model of apoptosis. PLoS Comput. Biol. 5(12): e1000595.

13. Mai Zhongxing, Liu Haiyan (2009) Boolean network-based analysis of the apoptosis network: irreversible apoptosis and stable surviving. J. Theor. Biol. 259: 760–769.

14. Hanahan D, Weinberg RA (2000) The hallmarks of cancer. Cell 100: 57–70.

15. Hanahan D, Weinberg RA (2011) Hallmarks of cancer: The next generation. Cell 144: 646–674.

16. Vogelstein B, Kinzel KW (2004) Cancer genes and the pathways they control. Nat. Med. 10(8): 789–799.

17. Oda K, Matsuoka Y, Funahashi A, Kitano H (2005) A comprehensive pathway map of epidermal growth factor receptor signaling. Mol. Syst. Biol. 1: 8–24.

18. Calzone L, Gelay A, Zinovyev A, Radvanyi F, Barillot E (2008) A comprehensive modular map of molecular interactions in RB/E2F pathway. Mol. Syst. Biol 4: 173.

19. Newman MEJ (2010) Networks: An Introduction. New York: Oxford University Press.

20. Evan GI, Vousden KH (2001) Proliferation, cell cycle, and apoptosis in cancer. Nature 411: 342–348.

21. Lendahl U, Lee KL, Yang H, Poellinger L (2009) Generating specificity and diversity in the transcriptional response to hypoxia. Nat. Rev. Genet. 10: 821–832.

22. Bunz F (2008) Principles of cancer genetics. Berlin: Springer.

23. Luo J, Solimini NL, Elledge S (2009) Principles of cancer therapy: oncogene and non-oncogene addiction. Cell 136: 823–837.

24. Santarius T, Shipley J, Brewer D, Stratton MR, Cooper CS (2010) A census of amplified and overexpressed human cancer genes. Nat. Rev. Cancer 10: 59–64.

25. Stratton MR, Campbell DJ, Futreal PA (2009) The cancer genome. Nature 458: 719–724.

26. Stephens PJ, Tarpey PS, Davies H, Loo PV, Greenman C, et al. (2012) The landscape of cancer genes and mutational processes in breast cancer. Nature 486: 400–404.

27. Banerji S, Cibulskis K, Rangel-Escareno C, Brown KK, Carter SL, et al. (2012) Sequence analysis of mutations and translocations across breast cancer subtypes. Nature 486: 405–409.

28. Sjöblom T, Jones S, Wood LD, Parsons DW, Lin J, et al. (2006) The consensus coding sequences of human breast and colorectal cancers. Science 314: 268–274.

29. The Cancer Genome Atlas Research Network (2011) Integrated genomic analysis of ovarian carcinoma. Nature 474: 609–615.

30. Pleasance ED, Stephens PJ, O'Meara S, McBride DJ, Meynert A, et al. (2010) A small-cell lung cancer genome with complex signatures of tobacco exposure. Nature 463: 184–190.

31. The Cancer Genome Atlas Research Network (2008) Comprehensive genomic characterization defines human glioblastoma genes and core pathways. Nature 455: 1061–1068.

32. The Cancer Genome Atlas Research Network (2008) Comprehensive genomic characterization of squamous cell lung cancers. Nature 489: 519–525.

33. The Cancer Genome Atlas Research Network (2012) Comprehensive molecular characterization of human colon and rectal cancer. Nature 467: 330–337.

34. Greaves M, Maley CC (2012) Clonal evolution in cancer. Nature 481: 306–313.

35. Ferreira Jr SC, Martins ML, Vilela MJ (2002) Reaction-diffusion model for the growth of avascular tumor. Phys.Rev.E 65: 021907.

36. Paiva LR, Binny C, Ferreira Jr SC, Martins ML (2009) A multiscale mathematical model for oncolytic virotherapy. Cancer Res. 69: 1205–1211.

37. KEGG Pathway database, Kyoto Encyclopedia of Genes and Genomes. Available at http://www.genome.jp/kegg/pathway.html.

38. Vanhaesebroeck B, Stephens L, Hawkins P (2012) PI3K signalling: the path to discovery and understanding. Nat. Rev. Mol. Cell Biol. 13(3): 195–203.

39. Keith B, Johnson RS, Celeste Simon M (2012) HIF1 and HIF2: sibling rivalry in hypoxic tumour growth and progression. Nat. Rev. Cancer 12(1): 9–22.

40. Yarden Y, Pines G (2012) The ERBB network: at last, cancer therapy meets systems biology. Nat. Rev. Cancer 12(8): 553–563.

41. Massague J (2012) TGF signalling in context. Nature Rev. Mol. Cell Biol. 13: 616–630.

42. Anastas JN, Moon RT (2013) WNT signalling pathways as therapeutic targets in cancer. Nat. Rev. Cancer 13: 11–26.

43. Cotter TG (2009) Apoptosis and cancer: the genesis of a research field. Nat. Rev. Cancer 9: 501–507.

44. Batchelor E, Loewer A, Lahav G (2009) The ups and downs of p53: understanding protein dynamics in single cells. Nat. Rev. Cancer 9: 371–377.

45. Balkwill F (2009) Tumour necrosis factor and cancer. Nat. Rev. Cancer 9: 361–371.

46. Karin M, Cao Y, Greten FR, Li ZW (2002) NF-B in cancer: from innocent bystander to major culprit. Nat. Rev. Cancer 2: 301–310.

47. Smoot M, Ono K, Ruscheinski J, Wang PL, Ideker T (2010) Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics 27(3): 431–432. Cytoscape is available for download at http://www.cytoscape.com.