Dynamic network and epigenetic landscape model of a regulatory core underlying spontaneous immortalization and epithelial carcinogenesis

Dynamic network and epigenetic landscape model of a regulatory core underlying spontaneous immortalization and epithelial carcinogenesis

Méndez-López LF1,2,†    Davila-Velderrain J2,†    Enríquez-Olguín C3
Domínguez-Hüttinger E1,2
   Martinez-Garcia JC3*    Alvarez-Buylla ER1,2,*
Abstract

Tumorigenic transformation of human epithelial cells in vitro has been described experimentally as the potential result of a process known as spontaneous immortalization. In this process a generic series of cell–state transitions occur in which normal epithelial cells acquire a senescent state, later surpassed to attain first a mesenchymal state and then a final mesenchymal stem–like phenotype, with potential tumorigenic behavior. In this paper we integrate published data on the molecular components and interactions that have been described as key regulators of such cell states and transitions. A large network, that is provided, is constructed and then reduced with the aim of recovering a minimal regulatory core incorporating the necessary and sufficient restrictions to recover the observed cell states and their generic progression patterns in epithelial–mesenchymal transition. Data is formalized into logical regulatory rules that govern the dynamics of each of the network’s components as a function of the states of its regulators. The proposed core gene regulatory network attains only three steady–state gene expression configurations that correspond to the profiles characteristic of normal epithelial, senescent, and mesenchymal stem–like cells. Interestingly, epigenetic landscape analyses of the uncovered network shows that it also recovers the generic time–ordered transitions documented during tumorigenic transformation in vitro of epithelial cells. The latter strongly correlate with the patterns observed during the progressive pathological description of epithelial carcinogenesis in vivo.

1 Instituto de Ecología, Universidad Nacional Autónoma de México, Cd. Universitaria, México, D.F. 04510, México
2 Centro de Ciencias de la Complejidad (C3), Universidad Nacional Autónoma de México, Cd. Universitaria, México, D.F. 04510, México
3 Departamento de Control Automático, Instituto Politécnico Nacional, A. P. 14-740, 07300 México, DF, México

Corresponding authors: juancarlos_martinez-garcia@conciliencia.org, eabuylla@gmail.com
These authors contributed equally to this work

Introduction

Nearly 84% of cancers diagnosed in human adults are carcinomas (i.e., cancer of epithelial origin), and their emergence is strongly associated with both an underlying chronic inflammatory process and with aging [84]. The precise role and the contribution of these two processes to the origin, progression, and detected clinic behavior of epithelial cancers remains elusive, however. The current general assumption is that aging and inflammation increase the chance of accumulating somatic mutations, and this genetic instability ultimately leads to carcinoma. However, this view does not offer a logical or mechanistic explanation for well–documented observations. For example: (1) cancer cells show morphological and transcriptional convergences despite their diverse origin, (2) carcinogenesis recapitulates embryonic processes, (3) cancer behavior can be acquired in the absence of mutations through trans– or dedifferentiation, and (4) cancer cells can be “normalized” by several experimental means [85, 86, 87, 88]. Moreover, it is well–known that different carcinomas share the same cellular processes and histological stages or progression patterns, as well as robust associations with lifestyle factors [89]. These empirical observations suggests that, in analogy to normal development, the human genome is associated with an underlying robust mechanism restricting cell states and temporal progression patterns that are characteristic of epithelial carcinogenesis. In accordance with this view, other researchers have previously proposed that cancer can be considered a developmental disease [90, 91].

In systems biology it is common to understand both cell differentiation and development in terms of dynamical systems theory. In this framework, the genome of a cell is directly mapped into a global and multi–stable gene regulatory network (GRN) whose dynamics yields several (quasi)stationary and stable distinct phenotypic cellular states [92, 93, 94, 95, 96, 97]. That is, the same genome robustly generates multiple discrete cellular phenotypes through developmental dynamics [95, 98, 99]. These stable phenotypic states are called attractors and correspond to configurations of gene or protein activation states that underlie the cellular fates or phenotypes – i.e., which thus constitute biological observables. Therefore, developmental processes – cellular differentiation events in particular – are formalized in temporal terms as attractor’s (i.e., cell states) transitions. Here we adopt such approach to study the cell states attained and the time–ordered transitions observed during the tumorigenic transformation of epithelial cells cultured in vitro that surpass a senescent state; a process known as spontaneous immortalization.

Experimental findings in molecular and cell biology of cancer research have revealed that it is possible to recover cells with cancer–like phenotypes through some specific cellular transitions. This has been shown particularly in carcinomas [86, 100, 101, 102]. By a cellular transition we refer to a differentiation event in which a certain cell acquires a discretely different cellular phenotype. For example, the process called epithelial–mesenchymal transition (EMT) comprises a stereotypical cell state transition in which epithelial cells exposed, for example, to cytokines, are induced to undergo a discrete phenotypic change acquiring a mesenchymal phenotype [102, 100]. Interestingly, through inflammation–induced EMT epithelial cells surpass senescence, and undergo spontaneous immortalization. Cells that emerge from this process manifest mesenchymal stem–like properties and are capable of developing cancer in murine models [86, 101]. Furthermore, these cells are difficult to distinguish phenotypically and in terms of the transcription factors that they express from either the so–called cancer stem cells (also known as tumor initiating cells) or from embryonic stem cells [103, 104].

In the present work, we hypothesize that a generic series of cell state transitions widely observed and robustly induced by inflammation in cell cultures during spontaneous immortalization naturally result from the self–organized behavior emerging from an underlying intracellular GRN. During this process, normal epithelial cells first acquire a senescent state, to finally attain a mesenchymal stem–like cellular state with a potential tumorigenic behavior. We speculate that tissue–level conditions associated with a bad prognosis, such as a pro–inflammatory milieu, may increase the rate of occurrence of these same transitions in vivo promoting as a result the emergence and progression of epithelial cancer.

In an attempt to provide mechanistic insights into the regulation of the aforementioned observed cell–fates specification, as well as the time–ordered cell–state transitions, we propose here a cellular level GRN model that integrates the available experimental data concerning the main molecular components and interactions related to the emergence and progression of carcinomas. We propose a large GRN of 41 nodes that integrates cellular processes thoroughly studied experimentally, but which have not been integrated before into a single GRN. Specifically, the large GRN model includes key molecular regulators that: (1) characterize the cellular phenotypes of epithelial, mesenchymal, and senescent cells; (2) are involved in the induction of the cellular processes of replicative senescence, cellular inflammation, and EMT; and (3) characterize the phenotypic changes undergone by cells emerging from these processes (i.e. mesenchymal stem–like cells). To obtain a minimal regulatory core for further dynamical analyses we formally reduced the large GRN. We show that the proposed regulatory core module displays an orchestrating robust behavior akin to that seen in other developmental regulatory modules previously characterized with similar formal approaches (see, for example [92, 93, 105, 106]). Specifically, by proposing logical functions grounded on experimental data for this regulatory core module and by analyzing its behavior following conventional Boolean GRN dynamical approaches, we show that the uncovered minimal GRN converges only to three attractors. The uncovered states correspond to the expected gene expression configurations that have been observed for normal epithelial, senescent and stem–like mesenchymal cellular fates. Additionally, we also explore the GRN Epigenetic Landscape using a stochastic version of the model (following: [107, 108]) in order to address if the proposed GRN also restricts or underlies the generic temporal sequence with which cell states occur in cell cultures and which correlate with observed patterns of cell–type enrichment during pathological descriptions of carcinoma progression.

Results

Gene Regulatory Network Construction

Following a bottom–up and an expert knowledge approach we propose a set of cellular dynamical processes ubiquitous to epithelial carcinogenesis, namely: replicative cellular senescence, inflammation, and epithelial–mesenchymal transition (EMT). The cellular phenotypes epithelial, senescent, and mesenchymal cell–types – as well as a mesenchymal embryonic–like state; have been largely characterized as biological observables involved in such processes. We provide further definitions of these – and associated – phenotypes and processes in our complementary Text S1. We take this information as a methodological basis to integrate a generic dynamical network model of epithelial carcinogenesis. As a first step in network integration, based on an extensive literature search (see Methods and Text S1), we assembled a set of transcription factors (TFs) and additional molecules involved in the establishment and regulation of these cellular states and processes. Subsequently, we manually retrieved documented regulatory interactions among the molecules, considering only those supported by experimental evidence. For a detailed description of the published information for each interaction proposed see Text S1. The constructed large GRN is shown in Figure 1 (see Methods). TFs are represented in graphical terms by squares and the rest of the molecules by circles. The identified large network consists of 41 nodes and 97 interactions; it includes 12 TFs which can be considered as key regulators of the processes under consideration. Colors indicate the association that each node hold with specific cellular phenotypes or processes being considered: epithelial (green), mesenchymal (orange), inflammation (red), senescence and DNA damage (blue), cell–cycle (purple), and polycomb complex (yellow).

The Proposed Network is Enriched with Cancer Pathways

In order to provide additional partial support for the association of the bio-molecular set of regulatory interactions that we have manually curated based on published data with the processes under consideration, as well as with carcinoma, we performed a network–based gene set enrichment analysis (GSEA) (see Methods). Among the 13 pathways or processes reported as significant when taking the KEGG database as a reference, 9 (69%) correspond to cancer pathways, namely: Bladder cancer, Chronic myeloid leukemia, Pancreatic cancer, Glioma, Non–small cell lung cancer, Melanoma, Small cell lung cancer, Prostate cancer, and Thyroid cancer – note that 6 (66.6%) of these correspond to carcinomas. On the other hand, when taking the GO Biological Process database as reference, among the significant results we found: replicative senescence, cellular senescence, cell aging, activation of NF–B–inducing kinase activity, determination of adult life span, epithelial cell differentiation, and positive regulation of NF–B transcription factor activity (see Table 1). Using network topological gene set analysis (see Methods) we found that, in addition to pathway enrichment, the topological signature of the molecules in the proposed network also shows a topological signature that is similar to the one shown by reference cancer pathways included in the KEGG database (see Figure S1). These results provide partial support for the proposed molecular players: given the current state of knowledge according to annotated databases, the set of molecules manually included in the proposed large network seems to be representative of the cellular phenotypes and processes considered as prior biological knowledge in our model. In addition, the molecular components included in the proposed large network are tightly associated with reference pathways of epithelial cancers.

A Core Regulatory Network Module Underlying Spontaneous Immortalization

We performed a knowledge–based network reduction of the large GRN in Figure 1 in order to derive a smaller, core GRN module for which both a topology and architecture with fully defined logical functions could be established, and which could also be analyzed as a dynamical system (see Methods). In addition, such regulatory core should comprise the necessary and sufficient set of nodes and interactions that integrate the processes involved in the large network and that could explain, at least in part, the restricted set of the cell–states and time–ordered transitions among them during spontaneous immortalization and epithelial cancer emergence/progression. We were able to define a set of molecular species whose regulatory hierarchy, activity, and expression define the identity of the phenotypes of epithelial, mesenchymal, and senescent cells. We also converged to, and included, main regulators of replicative cellular senescence, inflammation–induced EMT, and determinants of an induced mesenchymal stem–like phenotype. Hence, after reduction we obtained a core GRN consisting of only 9 nodes: ESE–2, Snai2, NF–B, E2F, p53, p16, Rb, Cyclin, and Telomerase. Figure 1b shows the proposed core regulatory module (colored nodes) in the context of the larger proposed network. For details on how these 9 nodes were selected over the rest of the nodes see Text S1. In what follows we present a brief description of the nodes included in the reduced GRN, as well as some of the key molecular mechanism encoded in the regulatory logic. Although many of the nodes that are included in this regulatory core module have been thoroughly studied experimentally and in terms of their involvement in different types of cancer, the architecture and topology of the proposed regulatory core module is novel.

ESE–2

represents the activity of the TFs ESE–1, ESE–2, and ESE–3 (also known as ESX, E74–like factor 5, and EHF; respectively) – for a table with synonyms Table E1 in supplementary file. These proteins belong to the subgroup ESE (i.e. epithelium–specific) of the TF family ETS. ESE–2 promotes its own expression and the expression of the other ESE TFs [109, 110, 111]. On the other hand, ESE–2 represses Snai2 – one of the main EMT promoting TFs – expression by direct interaction with its promoter region [112].

p16

represents the activity of the INK4b–ARF–INK4a locus, which encodes for the proteins p16 and p14. Cellular senescence is molecularly characterized by the expression of the proteins p16 and p53 [113]. p16 indirectly inhibits E2F by inhibiting cyclins CDK 2,4 and 6, which in turn inhibit Rb [114, 115]. On the other hand, the INK4b–ARF–INK4a – and thus p16 – is regulated by the activity of Polycomb–group proteins by means of promoter hypermethylation [116].

p53

represents the protein with the same name. The shortage of telomeric DNA seems to be recognized as DNA damage promoting the activation of p53. In senescence, the activity of p16 and p53 over Rb, E2F and Cyclins invariably arrests the cell–cycle in the phase G1/G2 [117, 118].

Rb

represents the cell–cycle regulator with the same name. Rb prevents cycle progression by forming a complex with the TF E2F [119].

E2F

represents the TF with the same name. E2F regulates critical genes for adequate cell–cycle progression.

Cyclin

represents the activity of the complex Cyclin–dependent kinases (CDKs) known to inactivate Rb by phosphorylation. The latter, in turn, promotes the activity of E2F and cell–cycle progression [120].

NF-B

represents cellular inflammation by the activity of the TF NF–B. Accordingly, with this node we also represent the effect of the cytokines transforming growth factor–beta (TGF–), interleukin–6 (IL–6), and tumor necrosis factor alpha (TNF–). These three factors converge in the activation of NF–B by phosphorylating the inhibitor IB [121, 122].

TELasa

represents the enzyme telomerase. This enzyme is responsible for the de novo synthesis of telomeres. Most human cell–types do not express telomerase; however, it is expressed on immortalized epithelial cells, and it is thought to be responsible for telomere extension in tumors [123].

Snai2

this node includes the activity of the main TFs known to be directly associated with EMT regulation, namely: Snai2 (Slug), Snail, Twist1, Twist2, ZEB1, ZEB2, and FOXC2. These TFs repress (induce) the expression of genes specific to epithelial (mesenchymal) cells [124, 125]. It has been proposed that there is a regulatory hierarchy driving EMT in which Snail activates Snai2, Twist, Zeb, and FOXC2. The latter, in turn, regulates Snail and Snai2 in a positive manner [126, 124, 127, 128]. Regardless of a hierarchical interpretation, it is well–documented that these TFs maintain the mesenchymal phenotype in a coordinated fashion, showing co–expression patterns and regulatory crosstalk [127, 128]. It has been suggested that among these TFs, Snai2 may be the strongest suppressor of the epithelial phenotype [129]. However, we decided to represent the collective regulatory activity of the mesenchymal TFs using Snai2 based on the recent experimental demonstration of an antagonistic relation between Snai2 and ESE–2. Specifically, in vitro and in vivo studies showed that ESE–2 regulates the transcription of Snai2 [112].

According to our model reduction methodology, literature search, careful manual curation, and network–based enrichment analysis; we propose that the derived core GRN module (see Figure 2) includes a molecular set which is both necessary and sufficient to specify the identity of the aforementioned cellular phenotypes and to represent the main intracellular regulatory events driving spontaneous immortalization in a robust manner. We test our proposal by building and analyzing a mechanistic GRN dynamical model (see below).

Recovered Attractors of the Core GRN Module Correspond to Configurations that Characterize Expected Cellular Phenotypes

Based on the experimental data concerning the expression patterns of the genes incorporated in the proposed core GRN model in Figure 2 we assembled a table with a Boolean format of the state configurations expected to be recovered with the proposed GRN dynamical model. We refer to this configurations as the “expected attractors” – these correspond to the empirically observed genetic configurations. Furthermore, we integrated and formalized the experimental data concerning the interactions among the GRN nodes using Boolean logical functions that will rule the Boolean GRN dynamics and comprise the architecture of the proposed GRN. The set of formulated rules underlying the regulatory events is shown in Text S1 – each logical rule is presented both as a logical preposition and as a truth table. Using the set of nodes and their corresponding logical rules we completely define a mechanistic dynamical GRN model [130]. The exhaustive computer–based simulation analysis of this model (see Methods) recovered three fixed–point attractors. Interestingly, the recovered attractors showed perfect correspondence with the expected attractors representing cellular phenotypes (see Table 2). The three recovered attractors correspond to the expected epithelial, senescent, and mesenchymal stem–like phenotypes :

The normal epithelial cell phenotype

is represented by the attractor with ESE–2, E2F, Cyclin and NF–B activity. ESE–2 is an epithelial–specific TF which regulates a large number of genes specific to epithelial cells [131, 132]. NF–B shows ubiquitous expression through the different types of human cells; however, it is also positively regulated by TFs of the ESE family (i.e. ESE–2) [133]. Moreover, under inflammatory conditions the activity of NF–B is enhanced [134, 135]. On the other hand, E2F and Cyclin represent core regulators of cell–cycle entrance, and thus specify proliferative capability [136, 137].

The senescent cell phenotype

is represented by the attractor with ESE–2, Rb, p16, p53, and NF–B activity. Its biological counterpart would be an epithelial cell induced to replicative senescence, given (1) that it is expected to repress E2F [131]; and (2) that Rb, p16, p53, and NF–B are the molecular biomarkers of cellular senescence [138].

Messenchymal Stem-like phenotype

In the model proposed here, the attractor whose configuration shows Snai2, Cyclin, NF–B, and Telomerase activity – and inactivity of ESE–2, p16, Rb, p53, and E2F – would correspond to a mesenchymal stem–like phenotype with tumorigenic potential (see discussion below).

The correspondence between the recovered attractors and the expected cellular phenotypes strongly suggests that the proposed nine–node core GRN indeed constitutes a regulatory module that is robust to initial conditions and that comprises a set of necessary and sufficient components and interactions to restrict the system to converge to the cellular phenotypes observed during spontaneous immortalization.

Validation of the Uncovered Core Regulatory Module: Loss and gain–of–function Mutant and Robustness Analyses

In order to validate the Boolean GRN dynamics we tested if the same GRN module is able to recover observed attractors in loss and gain of function mutants. We simulated such mutants analogous to experimental observations reported in the literature. Specifically, we simulated loss– and gain–of–function mutations of ESE–2, Snai2, and p16 that have been reported in the literature. When simulating ESE–2 gain of function (by setting the expression state for this node permanently to “1” in the simulations), the GRN model recovers three attractors corresponding to three different phenotypes which have been experimentally described and are associated with ESE–2 over–expression: an epithelial senescent cell [139], a normal epithelial cell [112], and a metastable state with proliferative phenotype [140]. In the case of ESE–2 loss–of–function (simulated by setting the expression state of this node to “0” permanently), the model recovers an attractor corresponding to a mesenchymal phenotype, which is also consistent with observations [112]. For Snai2, gain–of–function simulation recovers one attractor corresponding to mesenchymal stem–like phenotype, which is consistent with observations from ectopic over–expression experiments of mesenchymal TFs [101, 141, 142]. Snai2 loss–of–function simulation, on the other hand, recovered two attractors corresponding to normal and senescent epithelial phenotypes, which is also consistent with observations [112, 143]. Finally, gain–of–function simulation of p16 recovered two attractors; one associated with a mesenchymal stem–like but incompletely senescent (due to the lack of p53) phenotype; the other corresponding to an epithelial senescent phenotype. The first prediction is consistent with the status of immortal and apoptosis–resistant shown by mesenchymal stem–like cells, as well as with the capability of mesenchymal TFs to abrogate senescence [144]. The second attractor is consistent with the potential for replicative senescence of epithelial cells. p16 loss–of–function simulation recovers two attractors corresponding to an epithelial cell and a mesenchymal stem–like cell. This prediction is consistent with the observed biological conditions for both phenotypes, where p16 is commonly repressed by polycomb proteins [145]. The recovered attractors in mutant conditions are shown in Figure S2 in supplementary file.

It is important to note that, given that the uncovered regulatory module uncovered here is the result of a model reduction methodology where we permissively chose to represent multiple molecular species by the activity of some of the nodes, a direct interpretation of mutant simulations is not straightforward. Consequently, care should be taken when interpreting the results of the simulations or making predictions of mutant phenotypes yet to be experimentally tested and further explored in the context of the larger GRN in Figure 1, which is the focus of an ongoing study. With this in mind, instead of simulating additional mutant conditions, we further validated the dynamical GRN model by testing its robustness to perturbations of the logical rules. Specifically, we tested the robustness of the predicted attractors by generating a large set of perturbed networks (e.g, 10,000), calculating their respective attractors, and then counting the occurrences of the original attractors within the perturbed set. We generated each perturbed network by choosing a function of the network at random and flipping a single bit in this function [146]. We performed four complementary in silico based experiments following this general robustness analysis. First, we estimated the fraction of occurrences of the three original attractors (i.e., their robustness). Then, we repeated the experiment three times, but each time estimating the robustness of each individual attractor. For these four experiments we estimated a robustness (i.e., fraction of times) of 0.7439, 0.905, 0.923, and 0.902, respectively. Hence, out of 10,000 random networks generated by in silico perturbations to the logical rules, a major fraction recovered the original attractors; as it is expected for a developmental (core) regulatory module that is robust both to transient (initial) and genetic perturbations [93]. This result supports the view that the core GRN uncovered here is indeed a regulatory network module driving developmental dynamics. It also constitutes a mechanistic explanation (for definitions, see [130]) to the generic cell phenotypes observed during spontaneous immortalization in vitro and which correlate with the cellular description of carcinoma progression in vivo (see below).

Attractor Time–Ordered Transitions: Epigenetic Landscape of the Uncovered GRN Core Module

During the tumorigenic transformation of epithelial cells in culture, a generic time–ordered series of cell state transitions is observed and robustly induced by inflammation [86, 101]. Normal epithelial cell become senescent cells, which afterwards overcome this latter state acquiring a final mesenchymal stem–like phenotype. Interestingly, during the progressive pathological description of epithelial carcinomas in vivo the temporal pattern with which each of these different cell phenotypes enriches the tissue seems to be tightly ordered and is also generic to all types of such cancers irrespective of the tissue where they first appear. In order to test if the uncovered GRN core module not only underlies and restricts the types of cell phenotypes (attractors) but also their time–ordered transitions, following [108] we explored its associated Epigenetic Landscape (EL) by implementing a discrete stochastic model as an extension to the Boolean network model [95] (see Methods). By means of computer–based simulations we performed two analyses in order to uncover functional and structural constraints in attractor transitions. (1) We explored the temporal sequence of attractor attainment, and (2) we calculated the consistent global ordering of all the given attractors. Specifically, following [107], we found that the most probable temporal order of attractor attainment for a cell (population) initially on epithelial state is:

see Fig 4a. On the other hand, following [147] we defined a consistent global ordering of the uncovered attractors based on their relative stability (see Methods). Relative stability calculations are based on the mean first passage time (MFPT) between pairs of attractors. These, in turn, epitomize barrier heights in the EL by approximating a measure for the ease of specific transitions. Similar to the previous analysis, the uncovered global ordering of attractors is Epithelial Senescent Mesenchymal stem–like (Fig 4b). This corresponds to the only order in which the system can visit the three attractors following a positive net transition rate. These results indicate that, when considering only intracellular regulatory constraints alone, the uncovered GRN core module structures the epigenetic landscape in a way that a specific flow across the landscape is preferentially and robustly followed. We anticipate that observed transition rates in vivo are likely to depend on tissue–level processes and/or additional GRN components underlying epithelial cell sub–differentiation, that have not been considered here. These latter restrictions will be modeled in future contributions building up on the framework that has been put forward here.

Discussion

Multicellularity by definition implies a one–to–many genotype–phenotype map. The genome of a multicellular individual possesses the intrinsic potentiality to implement a developmental process by which all its different cell–types and tissue structures are ultimately established. In the last decades, a quite coherent theory to explain the development of multicellular organisms as the result of the orchestrating role of GRNs has been developed [92, 94, 95]. The main conclusion is that observable cell states emerge from the self–consistent multistable regulatory logic dictated by genome structure and obeyed by (mainly) transcription factors (TFs) resulting in stable, steady–states of gene expression. Cancer development and progression is also a phenomenon intrinsic to multicellular organisms. Furthermore, similar to normal development, cancer is robustly established as evidenced by its directionality and phenotypic convergence [85]. Is cancer somehow orchestrated by GRN dynamics as well? Several hypothesis have been presented in this direction such as the cancer attractor theory [91, 85], and the endogenous molecular cellular network hypothesis [148, 149]. In this contribution we also follow the viewpoint of an intrinsic regulatory network, but we focus on a specific developmental process at the cellular level: the robust cell state transitions observed during the tumorigenic transformation of human epithelial cells in culture induced by inflammation and resulting from surpassing a senescent state through EMT – i.e., tumorigenic transformation due to spontaneous immortalization. We propose that a mechanistic understanding of this process is an important first necessary step to unravel key cellular processes which might be occurring in vivo, where its rate of occurrence is likely to be regulated by tissue–level and systemic conditions directly linked with lifestyle choices, as well as additional regulatory interactions underlying epithelial cell sub–differentiation.

A Generic Molecular Regulatory Network

The predominant strategy in the molecular study of cancer and cellular tumorigenic transformation has been to focus on pathways and associated mutations. Aware that signaling pathways are actually embedded in complex regulatory networks here we assembled from curated literature a GRN comprising the main molecular regulators involved in key cellular processes ubiquitous to carcinogenesis following a bottom–up approach (see results). Subsequently, we followed a mechanistic approach to address the question of whether we assembled a set of necessary and sufficient molecular players and interactions to recover the cellular phenotypes and processes documented during the spontaneous immortalization of human epithelial cells in culture: we proposed, analyzed and validated an experimentally grounded core GRN dynamical model.

Small developmental regulatory modules have been shown to successfully include the necessary and sufficient set of components and interactions for explaining, as manifestations of intrinsic structural and functional constraints imposed by these GRNs, the dynamics of complex processes such as stem cell differentiation [150], cell–fate decision [151] and similar cellular processes during plant morphogenesis [92, 93, 107, 105]. We hypothesized that a similar core developmental module can be formulated in an attempt to explain the cell–fates observed during spontaneous immortalization of human epithelial cells in vitro resulting in a potentially tumorigenic state. In order to show this, we first reduced the proposed larger network into a regulatory core module, by eliminating transitory pathways within the network and by including compounded nodes while maintaining the core network structure and without affecting the dynamical output during each reduction step (for details, see Methods). We obtained a small set of main molecular players (Fig 2). We extracted from available literature the expression profiles of the generally observable cell states of interest in terms of this minimal set of molecules (see Table 2). Given our main hypothesis, we tested if the reduced molecular set and their regulatory logic formalized as a Boolean GRN model were able to recover the biologically observable expression profiles as stationary and stable network configurations (i.e., attractors). Interestingly, we found that the core GRN model only converges to the observed gene expression profiles in wild–type (see Table 2) and some mutant backgrounds (see results). This result strongly suggest that we have successfully included the key regulators and interactions at play during the establishment of cell states observed during the tumorigenic transformation of human epithelial cells resulting from spontaneous immortalization.

It is noteworthy that our model does not include any hypothetical interaction or component, a common practice in GRN modeling [93, 105, 151]. Our GRN model exclusively integrates available published experimental data; indeed, it was a surprising result that the observed dynamical behavior emerged naturally under such conditions. This suggests that despite incomplete information, there is enough molecular data to uncover important restrictions underlying cell behavior during transitions relevant to epithelial carcinogenesis. Consequently, we consider that the networks reported herein (both the large and the core GRNs) may serve as bona fide base models useful to integrate novel discoveries, as well as components underlying epithelial cellular sub–differentiation, while following a bottom–up approach in cancer network systems biology.

Attractor Time–Ordered Transitions

Discrete GRN models can be used to integrate regulatory mechanisms that not only recapitulate the observed gene expression patterns, but that also reproduce the observed developmental time–ordering of cell phenotypes. This can be done by considering stochasticity in the model in order to explore [95, 106, 108] and/or characterize [147] the associated EL. Importantly, by exploring noise–induced transitions we do not assume that noise alone is the driving force of the transitions, instead, we exploit noise as a tool to explore the GRN–based version of Waddington’s EL and to indirectly characterize its structure. Specifically, by calculating the relative stability of the attractors (see Methods) we approximate the in–between attractor barrier heights in the landscape. Furthermore, measures of relative stability can also be exploited to calculate net transition rates measuring the ease of specific inter–attractor transitions and to uncover the predominant developmental route across the epigenetic landscape [152]: ordered transitions sharing positive net transition rates will be preferentially followed. Our results show that such a developmental route follows the time–order of cellular phenotypic states epithelialsenescentmesenchymal stem–like (potentially tumorigenic). In other words, the constraints imposed by the GRN structure the associated EL in such a way that an epithelial cell in culture as a “ball” would naturally roll following such a path, in agreement with the observed spontaneous immortalization process.

Even in the case of the simple model presented here, it is interesting that of the many possible cell states and developmental routes, the core GRN network is canalized to the few steady–states and the developmental time–ordering consistent with the molecular characterization of cell phenotypes observed during spontaneous immortalization and correlating with carcinoma progression in vivo (see below). This suggests that specific progressive alterations or particular “abnormal” signaling mechanisms are not necessarily required for a cell to reach a potentially tumorigenic state. Additionally, robustness analysis performed on the same network showed that the recovered attractors are also robust to permanent alterations of the regulatory logic.

From Abstract Network Attractors and Dynamics to Biological Insight

We are aware of the high degree of simplification involved in the model proposed herein. Accordingly, we do not attempt to present it as a source of accurate predictions for either the occurrence or the future behavior of a phenomena as complex as carcinogenesis. Instead, we formulate the model in an attempt to provide some intuition into otherwise highly complicated processes, and to illuminate increasing body of confounding descriptions. Simple mechanistic models like the one presented here sacrifice detail and accuracy in exchange for understanding [153, 130]. What biological insights can be gained by the uncovered GRN dynamical model? Our simple GRN model strongly suggests that the generic series of cell state transitions widely observed and robustly induced by inflammation in cell culture from normal epithelial to immortalized senescent cells, and from this latter state to a final mesenchymal stem–like phenotype in the process defined as spontaneous immortalization naturally result from the self–organized behavior emerging from an underlying GRN novel architecture and topology.

Importantly, cells that emerge from spontaneous immortalization induced by cytokines display mesenchymal stem like phenotype and tumorigenic behavior – i.e., repress proteins p16 and p53, surpass senescence, and re–express telomerase [101]. Phenotypically, these cells are difficult to distinguish from the so–called cancer stem cells, tumor initiating cells or embryonic stem cells [103, 104]; are resistant to apoptosis; and have the ability to migrate and generate metastasis and form secondary tumors – all lethal traits characterizing cancer cells [86]. We, thus, speculate that tissue–level conditions associated with a bad prognosis, such as a pro–inflammatory milieu, may increase the rate of occurrence of these same transitions in vivo promoting as a result the development and progression of epithelial cancer. We substantiate this view by noting several independent empirical observations. (1) Histological diagnosis of carcinoma are generally preceded by a lesion called hyperplasia; senescent cells are abundant in hyperplasias and scarce in carcinomas [154]. (2) During chronological aging senescent cells increase in number within both normal tissues and hyperplasias. (3) Senescence is associated with the promotion of carcinogenesis by contributing with the loss of tissue architecture and promoting an inflammatory milieu [155]. (3) Overcoming the senescent barrier is fundamental in tumor progression [156, 157]. (4) The EMT process constitutes a well–characterized mean to overcome senescence under an inflammatory environment([158]).

We must point out, however, that transition rates during spontaneous immortalization, if occurring in vivo, may be regulated by tissue–level, self–organizational processes not considered in our cellular level model. For example, the likelihood of spontaneous immortalization in vivo may be increased by extracellular perturbations that inevitably occur during aging; mainly, by inflammation and tissue remodeling resulting from an increased population of senescent cells. The cellular level network models reported here are, nevertheless, a valuable building block for more detailed multi–level models integrating further sources of tissue–level constraints such as cell cycle progression, cell–cell interactions, differential proliferation rates, and mechanical forces.

Summarizing, in this contribution we propose an experimentally grounded GRN model for spontaneous immortalization. We report one large GRN model (41 nodes) and one core GRN developmental module (9 nodes), both useful and necessary for further integration of signaling and mechanical processes in multi–level, more detailed modeling efforts. We explore by analyzing the dynamical behavior of the latter if the uncovered GRN topology and architecture underlies the gene expression configurations that characterize normal epithelial, senescent, and mesenchymal stem–like cell–fates well documented during tumorigenic transformation in vitro and which correlate with those observed in the progressive pathological description of epithelial carcinogenesis in vivo. Overall, our results suggest that tumorigenic transformation in vitro due to spontaneous immortalization can be understood and modeled at a cellular level generically as a developmental system undergoing cell–state transitions resulting from the structural and functional constraints imposed, in part, by the interactions included in the proposed GRN. They also suggest that similar transitions may be occurring in vivo and might be relevant for carcinoma development and progression. This view is consistent with the robustness, generic patterns, and directionality observed during the development of human cancers derived from epithelial tissues. Particularly, based on our results, we hypothesize that replicative senescence and chronic inflammation are likely to increase the occurrence of spontaneous immortalization in vivo promoting the development of epithelial carcinogenesis. Testing such hypothesis awaits the development of multi–level models taking the ones presented here as building blocks, which is the subject of ongoing investigation.

Materials and Methods

Literature Search

A total of 159 references, considering both references in extended view material (see Text S1) and main text, were carefully and manually reviewed in order to first define a minimal set of cellular phenotypes and processes (for definitions, see Text S1) which enable a generic representation of epithelial carcinogenesis on the basis of cell state transition events. Subsequently, a set of associated, experimentally described molecular regulators was extracted from the literature, including their regulatory interactions.

Network Assembly

The network (see Fig. 1) was assembled manually by adding nodes (genes/proteins) and edges (activating or inhibitory interactions) describing direct mechanisms reported in the available literature to have an influence on both the specification of the cellular phenotypes and the development of the cellular process defined in (Text S1). The initial network was created based on experimentally grounded knowledge from 159 references (including reviews and research papers) and consists of 41 nodes and 97 edges. The literature included data known before 2014. Support for each of the proposed interactions is listed in Text S1.

Network–based Gene Set Enrichment Analysis

The bioinformatics tools EnrichNet [159] and TopoGSA [160] were used to perform network–based gene set enrichment analysis and topology–based gene set analysis, respectively. Briefly, EnrichNet maps the input gene set into a molecular interaction network and calculates distances between the genes and pathways/processes in a reference database. TopoGSA also maps the input gene set into a network, and then it computes its topological statistics and compares it against the topology of pathways/processes in a reference database. Here a connected human interactome graph extracted from the STRING database and the KEGG and GO Biological Process databases were used as reference molecular interaction network and databases. Both analyses were performed using the Cytoscape plugin Jepettp [161].

Network Reduction

In order to extract a representative core regulatory model from the initial network and to obtain a more computationally tractable one, which reasonably unfolds the regulatory pathways, a reduction methodology was followed based on certain simplifying assumptions – supported by previous results in molecular biology studies – and on mathematical results from dynamical systems and graph theory. Here we briefly describe the main steps. The step–by–step reduction process is included in Text S1.

Simplifying assumptions:
  • ESE-2 groups activities of ESE-1, ESE-3, EGF, Her-2/neu.

  • Snai2 groups activities of Snail, Twist (Twist, in turn, groups activities of Twist1 and Twist2), Zeb and FOXC2.

  • p16 groups p14 and NF–B node groups the inflammatory response activated by growth factors, mitogens and cytokines.

Reduction process:

(1) Simple mediator nodes (i.e., those nodes with in–degree and out–degree of one) were removed iteratively. (2) Nodes with in–degree of one and out–degree greater than one were removed iteratively. These steps (1 and 2) does not alter the attractors of the Boolean network under the asynchronous update, as mathematically proved in [162]. (3) Redundant interactions of selected nodes (based on biological arguments) resulting in self–regulation were included in single nodes/interactions (for details, see Text S1). (4) Selected nodes (based on biological knowledge again) with in–degree greater than one and out–degree of one were removed. The final steps (3 and 4) are supported by the mathematical analysis made in [163] in which the authors prove that the methodology preserves relevant topological and dynamical properties.

It is noteworthy that fixed point attractors are time–independent, so they are the same in both synchronous and asynchronous update methods. Complex attractors (in which the system oscillates among a set of states), on the other hand, depend on the update method. Consistently, the update method used in the model is irrelevant for the obtained results. This last assertion is valid because the model shows only fixed point attractors, which means, under the mathematically proved reduction methods applied, that the large network describes a qualitative long time behavior conserved in the reduced one. Besides, the methodology applied in order to obtain the reduced network enables the analysis of a resulting regulatory graph which is biologically meaningful and dynamically consistent with the network constructed with available molecular biology experimental data.

The final reduced network is shown in Figure 2. We refer to this network and its corresponding logical rules as the core regulatory module.

Dynamical Gene Regulatory Network Model

A Boolean network models a dynamical system assuming both discrete time and discrete state variables. This is expressed formally with the mapping:

(1)

where the set of functions are logical prepositions (or truth tables) expressing the relationship between the genes that share regulatory interactions with the gene , and where the state variables can take the discrete values or indicating whether the gene is expressed or not at a certain time , respectively. An experimentally grounded Boolean GRN model is then completely specified by the set of genes proposed to be involved in the process of interest and the associated set of logical functions derived from experimental data [106]. The set of logical functions for the core regulatory module used in this study is included in Text S1 – both as logical prepositions and truth tables. The dynamical analysis of the Boolean network model was conducted using the package BoolNet [146] within the R statistical programming environment (www.R-project.org).

Epigenetic Landscape Exploration

Including Stochasticity

In order to extend the Boolean Network into a discrete stochastic model and then study the properties of its associated EL, the so–called stochasticity in nodes (SIN) model was implemented following [107, 106, 108]. In this model, a constant probability of error is introduced for the deterministic Boolean functions. In other words, at each time step, each gene “disobeys” its Boolean function with probability . Formally:

(2)

The probability that the value of the now random variable is determined or not by its associated logical function is or , respectively.

Attractor Transition Probability Estimation

An attractor transition probability matrix with components:

(3)

representing the probability that an attractor is reached from an attractor , was estimated by numerical simulation following [107]. Specifically, for each network state in the state space () a stochastic one–step transition was simulated a large number of times (). The probability of transition from an attractor to an attractor was then estimated as the frequency of times the states belonging to the basin of the attractor were stochastically mapped into a state within the basin of the attractor .

Following the discrete time Markov chains (DTMCs) [164] theoretical framework, the estimated transition probability matrix was integrated into a dynamic equation for the probability distribution:

(4)

where is the probability distribution over the attractors at time , and is the transition probability matrix. This equation was iterated to simulate the temporal evolution of the probability distribution over the attractors starting from a specific initial probability distribution.

Attractor Relative Stability and Global Ordering Analyses

In addition to the calculation of the most probable temporal cell–fate pattern (see [107]), a discrete stochastic GRN model enables the study of the ease for transitioning from one attractor to another [152]. Specifically, a transition barrier in the EL epitomizes the ease for transitioning from one attractor to another. The ease of transitions, in turn, offers a notion of relative stability. It has recently been proposed that the GRN has a consistent global ordering of all cell attractors and intermediate transient states which can be uncovered by measuring the relative stabilities of all the attractors of a Boolean GRN [152, 147]. Here, the relative stabilities of the cell states were defined based on the mean first passage time (MFPT). Specifically, a relative stability matrix was calculated which reflects the transition barrier between any two states based on the MFPT. Here, in all cases, the MFPT was estimated numerically. Using the transition probabilities among attractors, a large number sample paths of a finite Markov chain were simulated. The MFPT from attractor to attractor corresponds to the averaged value of the number of steps taken to visit attractor for the fist time, given that the entire probability mass was initially localized at the attractor . The average is taken over the realizations. Following [152], based on the MFPT values a net transition rate between attractor and can be defined as follows:

(5)

This quantity effectively measures the ease of transition as a net probability flow. For all the calculation involving stochasticity, the robustness of the results was assessed by taking three different values for the probability of error . Stability of the results was assessed by manually changing the number of simulated samples until results become stable.

The consistent global ordering of all attractors uncovered with the core GRN was defined based on the formula proposed in [147]. Briefly, the consistent global ordering of the attractors is given by the attractor permutation in which al transitory net transition rates from an initial attractor to a final attractor are positive. This is schematically represented in Figure 4b. Calculated transition probability, MFPT, and net transition rate matrices are included in Text S2. R source code implementing all the calculations and analyses is available upon request.

Authors’ contributions

ERAB and JMG coordinated the study and with the other authors established the overall logic and core questions to be addressed. All the authors conceived and planned the modeling approaches. FML recovered the information from the literature to establish the model and provided expert knowledge in cancer biology. JDV established many of the specific analyses to be done, and programmed and ran all the modeling and analyses. FML and CEO formalized experimental data into regulatory logic. ERAB, JMG, JDV and FML participated in the interpretation of the results and analyses. JDV wrote most of the paper with help from ERAB and JMG and input from FML. All authors proofread the final version of the ms submitted.

Acknowledgments

This work was supported by grants CONACYT 240180, 180380, 167705, 152649 and UNAM-DGAPA-PAPIIT: IN203113, IN 203214, IN203814, UC Mexus ECO-IE415. J.D.V acknowledges the support of CONACYT and the Centre for Genomic Regulation (CRG), Barcelona, Spain; while spending a research visit in the lab of Stephan Ossowski. This article constitutes a partial fulfillment of the graduate program Doctorado en Ciencias Biomédicas of the Universidad Nacional Autónoma de México, UNAM in which J.D.V. developed this project. J.D.V receives a PhD scholarship from CONACYT. The authors acknowledge logistical and administrative help of Diana Romo.

References

  • [1] Anand P, Kunnumakara AB, Sundaram C, Harikumar KB, Tharakan ST, Lai OS, et al. Cancer is a preventable disease that requires major lifestyle changes. Pharmaceutical research. 2008;25(9):2097–2116.
  • [2] Huang S. On the intrinsic inevitability of cancer: from foetal to fatal attraction. In: Seminars in cancer biology. vol. 21. Elsevier; 2011. p. 183–199.
  • [3] Mani SA, Guo W, Liao MJ, Eaton EN, Ayyanan A, Zhou AY, et al. The epithelial-mesenchymal transition generates cells with properties of stem cells. Cell. 2008;133(4):704–715.
  • [4] Huang S. Non-genetic heterogeneity of cells in development: more than just noise. Development. 2009;136(23):3853–3862.
  • [5] Ben-Porath I, Thomson MW, Carey VJ, Ge R, Bell GW, Regev A, et al. An embryonic stem cell–like gene expression signature in poorly differentiated aggressive human tumors. Nature genetics. 2008;40(5):499–507.
  • [6] Kelloff GJ, Sigman CC. Assessing intraepithelial neoplasia and drug safety in cancer-preventive drug development. Nature Reviews Cancer. 2007;7(7):508–518.
  • [7] Virchow RLK. Cellular pathology. John Churchill; 1860.
  • [8] Huang S, Ernberg I, Kauffman S. Cancer attractors: a systems view of tumors from a gene network dynamics and developmental perspective. In: Seminars in cell & developmental biology. vol. 20. Elsevier; 2009. p. 869–876.
  • [9] Mendoza L, Alvarez-Buylla ER. Dynamics of the genetic regulatory network for Arabidopsis thaliana flower morphogenesis. J Theor Biol. 1998;193(2):307–319.
  • [10] Espinosa-Soto C, Padilla-Longoria P, Alvarez-Buylla ER. A gene regulatory network model for cell-fate determination during Arabidopsis thaliana flower development that is robust and recovers experimental gene expression profiles. The Plant Cell Online. 2004;16(11):2923–2939.
  • [11] Huang S, Kauffman S. Complex gene regulatory networks-from structure to biological observables: cell fate determination. Encyclopedia of Complexity and Systems Science Meyers RA, editors Springer. 2009;p. 1180–1293.
  • [12] Alvarez-Buylla ER, Azpeitia E, Barrio R, Benítez M, Padilla-Longoria P. From ABC genes to regulatory networks, epigenetic landscapes and flower morphogenesis: making biological sense of theoretical approaches. Seminars in cell & developmental biology. 2010;21(1):108–117.
  • [13] Huang S. Reprogramming cell fates: reconciling rarity with robustness. Bioessays. 2009;31(5):546–560.
  • [14] Kaneko K. Characterization of stem cells and cancer cells on the basis of gene expression profile stability, plasticity, and robustness. Bioessays. 2011;33(6):403–413.
  • [15] Huang S. The molecular and mathematical basis of Waddington’s epigenetic landscape: A framework for post-Darwinian biology? Bioessays. 2012;34(2):149–157.
  • [16] Davila-Velderrain J, Alvarez-Buylla ER. Bridging the Genotype and the Phenotype: Towards An Epigenetic Landscape Approach to Evolutionary Systems Biology. bioRxiv. 2014;.
  • [17] Xu J, Lamouille S, Derynck R. TGF--induced epithelial to mesenchymal transition. Cell research. 2009;19(2):156–172.
  • [18] Battula VL, Evans KW, Hollier BG, Shi Y, Marini FC, Ayyanan A, et al. Epithelial-Mesenchymal Transition-Derived Cells Exhibit Multilineage Differentiation Potential Similar to Mesenchymal Stem Cells. Stem Cells. 2010;28(8):1435–1445.
  • [19] Li CW, Xia W, Huo L, Lim SO, Wu Y, Hsu JL, et al. Epithelial–mesenchymal transition induced by TNF- requires NF-B–mediated transcriptional upregulation of Twist1. Cancer research. 2012;72(5):1290–1300.
  • [20] Morel AP, Lièvre M, Thomas C, Hinkal G, Ansieau S, Puisieux A. Generation of breast cancer stem cells through epithelial-mesenchymal transition. PloS one. 2008;3(8):e2888.
  • [21] Neph S, Stergachis AB, Reynolds A, Sandstrom R, Borenstein E, Stamatoyannopoulos JA. Circuitry and dynamics of human transcription factor regulatory networks. Cell. 2012;150(6):1274–1286.
  • [22] Azpeitia E, Benítez M, Vega I, Villarreal C, Alvarez-Buylla ER. Single-cell and coupled GRN models of cell patterning in the Arabidopsis thaliana root stem cell niche. BMC systems biology. 2010;4(1):134.
  • [23] Azpeitia E, Davila-Velderrain J, Villarreal C, Alvarez-Buylla ER. Gene regulatory network models for floral organ determination. In: Flower Development. Springer; 2014. p. 441–469.
  • [24] Álvarez-Buylla ER, Chaos Á, Aldana M, Benítez M, Cortes-Poza Y, Espinosa-Soto C, et al. Floral morphogenesis: stochastic explorations of a gene network epigenetic landscape. Plos one. 2008;3(11):e3626.
  • [25] Davila-Velderrain J, Martínez-García J, Alvarez-Buylla ER. Modeling the Epigenetic Attractors Landscape: Towards a Post-Genomic Mechanistic Understanding of Development. Name: Frontiers in Genetics. 2015;6:160.
  • [26] Zhou J, Ng A, Tymms MJ, Jermiin LS, Seth AK, Thomas RS, et al. A novel transcription factor, ELF5, belongs to the ELF subfamily of ETS genes and maps to human chromosome 11p13-15, a region subject to LOH and rearrangement in human carcinoma cell lines. Oncogene. 1998;17(21):2719–2732.
  • [27] Ma XJ, Salunga R, Tuggle JT, Gaudet J, Enright E, McQuary P, et al. Gene expression profiles of human breast cancer progression. Proceedings of the National Academy of Sciences. 2003;100(10):5974–5979.
  • [28] Escamilla-Hernandez R, Chakrabarti R, Romano RA, Smalley K, Zhu Q, Lai W, et al. Genome-wide search identifies Ccnd2 as a direct transcriptional target of Elf5 in mouse mammary gland. BMC molecular biology. 2010;11(1):68.
  • [29] Chakrabarti R, Hwang J, Blanco MA, Wei Y, Lukačišin M, Romano RA, et al. Elf5 inhibits the epithelial–mesenchymal transition in mammary gland development and breast cancer metastasis by transcriptionally repressing Snail2. Nature cell biology. 2012;14(11):1212–1222.
  • [30] Vernier M, Bourdeau V, Gaumont-Leclerc MF, Moiseeva O, Bégin V, Saad F, et al. Regulation of E2Fs and senescence by PML nuclear bodies. Genes & development. 2011;25(1):41–50.
  • [31] McConnell BB, Gregory FJ, Stott FJ, Hara E, Peters G. Induced expression of p16 INK4a inhibits both CDK4-and CDK2-associated kinase activity by reassortment of cyclin-CDK-inhibitor complexes. Molecular and cellular biology. 1999;19(3):1981–1989.
  • [32] Villacañas Ó, Pérez JJ, Rubio-Martínez J. Structural analysis of the inhibition of Cdk4 and Cdk6 by p16INK4a through molecular dynamics simulations. Journal of Biomolecular Structure and Dynamics. 2002;20(3):347–358.
  • [33] Bracken AP, Kleine-Kohlbrecher D, Dietrich N, Pasini D, Gargiulo G, Beekman C, et al. The Polycomb group proteins bind throughout the INK4A-ARF locus and are disassociated in senescent cells. Genes & development. 2007;21(5):525–530.
  • [34] Fang L, Igarashi M, Leung J, Sugrue MM, Lee SW, Aaronson SA. p21Waf1/Cip1/Sdi1 induces permanent growth arrest with markers of replicative senescence in human tumor cells lacking functional p53. Oncogene. 1999;18(18):2789–2797.
  • [35] Mao Z, Ke Z, Gorbunova V, Seluanov A. Replicatively senescent cells are arrested in G1 and G2 phases. Aging (Albany NY). 2012;4(6):431.
  • [36] Chellappan SP, Hiebert S, Mudryj M, Horowitz JM, Nevins JR. The E2F transcription factor is a cellular target for the RB protein. Cell. 1991;65(6):1053–1061.
  • [37] Byeon IJL, Li J, Ericson K, Selby TL, Tevelev A, Kim HJ, et al. Tumor Suppressor p16INK4A: Determination of Solution Structure and Analyses of Its Interaction with Cyclin-Dependent Kinase 4. Molecular cell. 1998;1(3):421–431.
  • [38] Beauséjour CM, Krtolica A, Galimi F, Narita M, Lowe SW, Yaswen P, et al. Reversal of human cellular senescence: roles of the p53 and p16 pathways. The EMBO journal. 2003;22(16):4212–4222.
  • [39] Freudlsperger C, Bian Y, Wise SC, Burnett J, Coupar J, Yang X, et al. TGF- and NF-B signal pathway cross-talk is mediated through TAK1 and SMAD7 in a subset of head and neck cancers. Oncogene. 2012;32(12):1549–1559.
  • [40] Harley C, Futcher A, Greider C. Telomeres shorten during ageing of human fibroblasts. Nature. 1990;345(6274):458–460.
  • [41] Mani SA, Yang J, Brooks M, Schwaninger G, Zhou A, Miura N, et al. Mesenchyme Forkhead 1 (FOXC2) plays a key role in metastasis and is associated with aggressive basal-like breast cancers. Proceedings of the National Academy of Sciences. 2007;104(24):10069–10074.
  • [42] Zeisberg M, Neilson EG, et al. Biomarkers for epithelial-mesenchymal transitions. The Journal of clinical investigation. 2009;119(6):1429–1437.
  • [43] Bolós V, Peinado H, Pérez-Moreno MA, Fraga MF, Esteller M, Cano A. The transcription factor Slug represses E-cadherin expression and induces epithelial to mesenchymal transitions: a comparison with Snail and E47 repressors. Journal of cell science. 2003;116(3):499–511.
  • [44] Dave N, Guaita-Esteruelas S, Gutarra S, Frias À, Beltran M, Peiró S, et al. Functional cooperation between Snail1 and twist in the regulation of ZEB1 expression during epithelial to mesenchymal transition. Journal of Biological Chemistry. 2011;286(14):12024–12032.
  • [45] Casas E, Kim J, Bendesky A, Ohno-Machado L, Wolfe CJ, Yang J. Snail2 is an essential mediator of Twist1-induced epithelial mesenchymal transition and metastasis. Cancer research. 2011;71(1):245–254.
  • [46] Hajra KM, Chen DY, Fearon ER. The SLUG zinc-finger protein represses E-cadherin in breast cancer. Cancer research. 2002;62(6):1613–1618.
  • [47] Davila-Velderrain J, Martinez-Garcia J, Alvarez-Buylla E. Descriptive vs. Mechanistic Network Models in Plant Development in the Post-Genomic Era. Plant Functional Genomics: Methods and Protocols. 2015;p. 455–479.
  • [48] Siegel R, Naishadham D, Jemal A. Cancer statistics, 2013. CA: a cancer journal for clinicians. 2013;63(1):11–30.
  • [49] Jemal A, Bray F, Center MM, Ferlay J, Ward E, Forman D. Global cancer statistics. CA: a cancer journal for clinicians. 2011;61(2):69–90.
  • [50] Stratton MR, Campbell PJ, Futreal PA. The cancer genome. Nature. 2009;458(7239):719–724.
  • [51] Hudson TJ, Anderson W, Aretz A, Barker AD, Bell C, Bernabé RR, et al. International network of cancer genome projects. Nature. 2010;464(7291):993–998.
  • [52] Weinstein JN, Collisson EA, Mills GB, Shaw KRM, Ozenberger BA, Ellrott K, et al. The cancer genome atlas pan-cancer analysis project. Nature genetics. 2013;45(10):1113–1120.
  • [53] Yaffe MB. The scientific drunk and the lamppost: massive sequencing efforts in cancer discovery and treatment. Science signaling. 2013;6(269):pe13.
  • [54] Creixell P, Schoof EM, Erler JT, Linding R. Navigating cancer network attractors for tumor-specific therapy. Nature biotechnology. 2012;30(9):842–848.
  • [55] DePinho RA, Polyak K. Cancer chromosomes in crisis. Nature genetics. 2004;36(9):932–934.
  • [56] Fujikawa M, Katagiri T, Tugores A, Nakamura Y, Ishikawa F. ESE-3, an Ets family transcription factor, is up-regulated in cellular senescence. Cancer science. 2007;98(9):1468–1475.
  • [57] Lee JM, Dedhar S, Kalluri R, Thompson EW. The epithelial–mesenchymal transition: new insights in signaling, development, and disease. The Journal of cell biology. 2006;172(7):973–981.
  • [58] Cano A, Pérez-Moreno MA, Rodrigo I, Locascio A, Blanco MJ, del Barrio MG, et al. The transcription factor snail controls epithelial–mesenchymal transitions by repressing E-cadherin expression. Nature cell biology. 2000;2(2):76–83.
  • [59] Sun Y, Song GD, Sun N, Chen JQ, Yang SS. Slug overexpression induces stemness and promotes hepatocellular carcinoma cell invasion and metastasis. Oncology Letters. 2014;7(6):1936–1940.
  • [60] Liu Y, El-Naggar S, Darling DS, Higashi Y, Dean DC. Zeb1 links epithelial-mesenchymal transition and cellular senescence. Development. 2008;135(3):579–588.
  • [61] Weinberg RA. Twisted epithelial–mesenchymal transition blocks senescence. Nature cell biology. 2008;10(9):1021–1023.
  • [62] Kim WY, Sharpless NE. The Regulation of¡ i¿ INK4¡/i¿/¡ i¿ ARF¡/i¿ in Cancer and Aging. Cell. 2006;127(2):265–275.
  • [63] Müssel C, Hopfensitz M, Kestler HA. BoolNet—an R package for generation, reconstruction and analysis of Boolean networks. Bioinformatics. 2010;26(10):1378–1380.
  • [64] Zhou JX, Samal A, d’Hèrouël AF, Price ND, Huang S. Relative Stability of Network States in Boolean Network Models of Gene Regulation in Development. arXiv preprint arXiv:14076117. 2014;.
  • [65] Wang G, Zhu X, Gu J, Ao P. Quantitative implementation of the endogenous molecular–cellular network hypothesis in hepatocellular carcinoma. Interface focus. 2014;4(3):20130064.
  • [66] Zhu X, Yuan R, Hood L, Ao P. Endogenous molecular-cellular hierarchical modeling of prostate carcinogenesis uncovers robust structure. Progress in biophysics and molecular biology. 2015;.
  • [67] Li C, Wang J. Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths. PLoS computational biology. 2013;9(8):e1003165.
  • [68] Zhou JX, Brusch L, Huang S. Predicting pancreas cell fate decisions and reprogramming with a hierarchical multi-attractor model. PloS one. 2011;6(3):e14752.
  • [69] Zhou JX, Qiu X, d’Herouel AF, Huang S. Discrete Gene Network Models for Understanding Multicellularity and Cell Reprogramming: From Network Structure to Attractor Landscapes Landscape. In: Computational Systems Biology Second Edition Elsevier. 2014;p. 241–276.
  • [70] Lander AD. The edges of understanding. BMC biology. 2010;8(1):40.
  • [71] Chen Z, Trotman LC, Shaffer D, Lin HK, Dotan ZA, Niki M, et al. Crucial role of p53-dependent cellular senescence in suppression of Pten-deficient tumorigenesis. Nature. 2005;436(7051):725–730.
  • [72] Campisi J. Cellular senescence: putting the paradoxes in perspective. Current opinion in genetics & development. 2011;21(1):107–112.
  • [73] Narita M, Lowe SW. Senescence comes of age. Nature medicine. 2005;11(9):920–922.
  • [74] Yildiz G, Arslan-Ergul A, Bagislar S, Konu O, Yuzugullu H, Gursoy-Yuzugullu O, et al. Genome-wide transcriptional reorganization associated with senescence-to-immortality switch during human hepatocellular carcinogenesis. PloS one. 2013;8(5):e64016.
  • [75] Smit MA, Peeper DS. Epithelial-mesenchymal transition and senescence: two cancer-related processes are crossing paths. Aging (Albany NY). 2010;2(10):735.
  • [76] Glaab E, Baudot A, Krasnogor N, Schneider R, Valencia A. EnrichNet: network-based gene set enrichment analysis. Bioinformatics. 2012;28(18):i451–i457.
  • [77] Glaab E, Baudot A, Krasnogor N, Valencia A. TopoGSA: network topological gene set analysis. Bioinformatics. 2010;26(9):1271–1272.
  • [78] Winterhalter C, Widera P, Krasnogor N. JEPETTO: a Cytoscape plugin for gene set enrichment and topological analysis based on interaction networks. Bioinformatics. 2014;30(7):1029–1030.
  • [79] Saadatpour A, Albert R, Reluga TC. A reduction method for Boolean network models proven to conserve attractors. SIAM Journal on Applied Dynamical Systems. 2013;12(4):1997–2011.
  • [80] Naldi A, Remy E, Thieffry D, Chaouiya C. Dynamically consistent reduction of logical regulatory graphs. Theoretical Computer Science. 2011;412(21):2207–2218.
  • [81] Allen LJ. An introduction to stochastic processes with applications to biology. CRC Press; 2010.
  • [82] Li J, Poi MJ, Tsai MD. Regulatory mechanisms of tumor suppressor P16INK4A and their relevance to cancer. Biochemistry. 2011;50(25):5566–5582.
  • [83] Yamakoshi K, Takahashi A, Hirota F, Nakayama R, Ishimaru N, Kubo Y, et al. Real-time in vivo imaging of p16Ink4a reveals cross talk with p53. The Journal of cell biology. 2009;186(3):393–407.

Tables

KEGG – Pathway or Process XD–score q-value Overlap/Size
Bladder cancer 1.1447 0 12/38
Chronic myeloid leukemia 0.86866 0 17/69
p53 signaling pathway 0.78477 0 14/62
Pancreatic cancer 0.68155 0 14/70
Glioma 0.68155 0 12/60
Non–small cell lung cancer 0.66586 0 10/51
Melanoma 0.65574 0 12/62
Small cell lung cancer 0.56447 0 14/82
Prostate cancer 0.54821 0 14/84
Cell cycle 0.54821 0 20/120
Cytosolic DNA–sensing pathway 0.48155 0.00001 6/40
Thyroid cancer 0.36155 0.00784 3/25
NOD-like receptor signaling pathway 0.35612 0.00001 7/59
GO Biological Process XD-score q–value Overlap/Size
replicative senescence 3.13328 0 8/10
cellular senescence 0.73328 0.02244 2/10
cell aging 0.43328 0.00608 3/24
activation of NF–B–inducing kinase activity 0.43328 0.04656 2/16
determination of adult lifespan 0.33328 0.40382 1/10
epithelial cell differentiation 0.32721 0.13188 2/33
positive regulation of NF–B transcription factor activity 0.30109 0 8/87
Table 1: Significant pathways and processes according to network–based gene set enrichment analysis
Cellular Phenotype Recovered Attractor (Active) “Expected Attractors” References
Epithelial Ese–2, NF-B, E2F, Cyclin Ese–2, NF–B, Cell Cycle(+) [112]
Senescent p16, p53, Ese–2, NF–B, Rb p16, p53, NF–B, Cell Cycle(-) [165, 166, 125]
Mesenchymal stem-like Snai2, Telomerase, NF–B, Cyclin Snai2, Telomerase, NF–B, Cell Cycle(+) [125, 112]
Table 2: Predicted and Observed Attractors

Figure Legends

Figure 1. Gene regulatory network for epithelial carcinogenesis. Nodes represent genes, and arrows (bars) represent experimentally characterized activation (arrow-heads) or repression (flat-heads) interactions. Genes corresponding to TFs are represented by squares and the rest by circles. (a) Colors indicate association with specific phenotypes and processes: epithelial (green), mesenchymal (orange), inflammation (red), senescence and DNA damage (blue), cell–cycle (purple), and polycomb complex (yellow). (b) Core gene regulatory module in the context of the global network. Colored nodes represent the final set of molecules obtained after the network reduction methodology was applied (see Methods) and which were included in the core GRN model.

Figure 2. Core gene regulatory network module for epithelial carcinogenesis Nodes represent either single or subsets of genes (see Results); arrows-heads represent activations and flat–heads repression interactions. Five of the nodes are involved in the specification of the cellular phenotypes: Epithelial (Ese–2), Senescent (p16, p53), and Mesenchymal stem–like (Snai2, TELasa). Three nodes are tightly associated with cell–cycle regulation (Rb, E2F, Cyclin), while node NF–B represents cellular inflammation.

Figure 3. The core gene regulatory module in the context of the Hallmarks of Cancer approach. The antagonistic activity state ESE–2 (-) and Snai2 (+) enable cells to sustain proliferative signals and evade growth suppressors by undergoing a dedifferentiation process. The state p16(-), Rb(-), p53(-), and TELasa (+) enable cell to acquire replicative immortality, resist cell death, as well as present genome instability and a mutation–prone phenotype by surpassing cellular senescence. High levels of cytokines and NF–B(+) expose cells to tumor promoting inflammation. The constitutive activity of Snai2(+) epitomizes the intrinsic phenotypic features of the cells emerging from the process of inflammation–induced EMT: activating invasion, avoiding immune destruction, and deregulating cellular energetics.

Figure 4. Temporal sequence and global order of cell–fate attainment pattern under the stochastic Boolean GRN model during epithelial carcinogenesis. (a) Maximum probability of attaining each attractor, as a function of time (in iteration steps). Vertical lines mark the time when maximal probability of each attractor occurs. The most probable sequence of cell attainment is: epithelial(E) senescent(S) mesenchymal(cancer–like)(M). The value of the error probability used in this case was . The same patterns were obtained with the 3 different error probabilities tested (data not shown). (b) Schematic representation of the possible transitions between pairs of attractors. Arrows indicate the directionality of the transitions. Above each arrow a sign () or () indicates whether the calculated net transition rate between the corresponding attractors is positive or negative. Red arrows represent the globally consistent ordering for the 3 attractors: the order of the attractors in which all individual transition has a positive net rate, resulting in a global probability flow across the EL.

References

  • [84] Anand P, Kunnumakara AB, Sundaram C, Harikumar KB, Tharakan ST, Lai OS, et al. Cancer is a preventable disease that requires major lifestyle changes. Pharmaceutical research. 2008;25(9):2097–2116.
  • [85] Huang S. On the intrinsic inevitability of cancer: from foetal to fatal attraction. In: Seminars in cancer biology. vol. 21. Elsevier; 2011. p. 183–199.
  • [86] Mani SA, Guo W, Liao MJ, Eaton EN, Ayyanan A, Zhou AY, et al. The epithelial-mesenchymal transition generates cells with properties of stem cells. Cell. 2008;133(4):704–715.
  • [87] Huang S. Non-genetic heterogeneity of cells in development: more than just noise. Development. 2009;136(23):3853–3862.
  • [88] Ben-Porath I, Thomson MW, Carey VJ, Ge R, Bell GW, Regev A, et al. An embryonic stem cell–like gene expression signature in poorly differentiated aggressive human tumors. Nature genetics. 2008;40(5):499–507.
  • [89] Kelloff GJ, Sigman CC. Assessing intraepithelial neoplasia and drug safety in cancer-preventive drug development. Nature Reviews Cancer. 2007;7(7):508–518.
  • [90] Virchow RLK. Cellular pathology. John Churchill; 1860.
  • [91] Huang S, Ernberg I, Kauffman S. Cancer attractors: a systems view of tumors from a gene network dynamics and developmental perspective. In: Seminars in cell & developmental biology. vol. 20. Elsevier; 2009. p. 869–876.
  • [92] Mendoza L, Alvarez-Buylla ER. Dynamics of the genetic regulatory network for Arabidopsis thaliana flower morphogenesis. J Theor Biol. 1998;193(2):307–319.
  • [93] Espinosa-Soto C, Padilla-Longoria P, Alvarez-Buylla ER. A gene regulatory network model for cell-fate determination during Arabidopsis thaliana flower development that is robust and recovers experimental gene expression profiles. The Plant Cell Online. 2004;16(11):2923–2939.
  • [94] Huang S, Kauffman S. Complex gene regulatory networks-from structure to biological observables: cell fate determination. Encyclopedia of Complexity and Systems Science Meyers RA, editors Springer. 2009;p. 1180–1293.
  • [95] Alvarez-Buylla ER, Azpeitia E, Barrio R, Benítez M, Padilla-Longoria P. From ABC genes to regulatory networks, epigenetic landscapes and flower morphogenesis: making biological sense of theoretical approaches. Seminars in cell & developmental biology. 2010;21(1):108–117.
  • [96] Huang S. Reprogramming cell fates: reconciling rarity with robustness. Bioessays. 2009;31(5):546–560.
  • [97] Kaneko K. Characterization of stem cells and cancer cells on the basis of gene expression profile stability, plasticity, and robustness. Bioessays. 2011;33(6):403–413.
  • [98] Huang S. The molecular and mathematical basis of Waddington’s epigenetic landscape: A framework for post-Darwinian biology? Bioessays. 2012;34(2):149–157.
  • [99] Davila-Velderrain J, Alvarez-Buylla ER. Bridging the Genotype and the Phenotype: Towards An Epigenetic Landscape Approach to Evolutionary Systems Biology. bioRxiv. 2014;.
  • [100] Xu J, Lamouille S, Derynck R. TGF--induced epithelial to mesenchymal transition. Cell research. 2009;19(2):156–172.
  • [101] Battula VL, Evans KW, Hollier BG, Shi Y, Marini FC, Ayyanan A, et al. Epithelial-Mesenchymal Transition-Derived Cells Exhibit Multilineage Differentiation Potential Similar to Mesenchymal Stem Cells. Stem Cells. 2010;28(8):1435–1445.
  • [102] Li CW, Xia W, Huo L, Lim SO, Wu Y, Hsu JL, et al. Epithelial–mesenchymal transition induced by TNF- requires NF-B–mediated transcriptional upregulation of Twist1. Cancer research. 2012;72(5):1290–1300.
  • [103] Morel AP, Lièvre M, Thomas C, Hinkal G, Ansieau S, Puisieux A. Generation of breast cancer stem cells through epithelial-mesenchymal transition. PloS one. 2008;3(8):e2888.
  • [104] Neph S, Stergachis AB, Reynolds A, Sandstrom R, Borenstein E, Stamatoyannopoulos JA. Circuitry and dynamics of human transcription factor regulatory networks. Cell. 2012;150(6):1274–1286.
  • [105] Azpeitia E, Benítez M, Vega I, Villarreal C, Alvarez-Buylla ER. Single-cell and coupled GRN models of cell patterning in the Arabidopsis thaliana root stem cell niche. BMC systems biology. 2010;4(1):134.
  • [106] Azpeitia E, Davila-Velderrain J, Villarreal C, Alvarez-Buylla ER. Gene regulatory network models for floral organ determination. In: Flower Development. Springer; 2014. p. 441–469.
  • [107] Álvarez-Buylla ER, Chaos Á, Aldana M, Benítez M, Cortes-Poza Y, Espinosa-Soto C, et al. Floral morphogenesis: stochastic explorations of a gene network epigenetic landscape. Plos one. 2008;3(11):e3626.
  • [108] Davila-Velderrain J, Martínez-García J, Alvarez-Buylla ER. Modeling the Epigenetic Attractors Landscape: Towards a Post-Genomic Mechanistic Understanding of Development. Name: Frontiers in Genetics. 2015;6:160.
  • [109] Zhou J, Ng A, Tymms MJ, Jermiin LS, Seth AK, Thomas RS, et al. A novel transcription factor, ELF5, belongs to the ELF subfamily of ETS genes and maps to human chromosome 11p13-15, a region subject to LOH and rearrangement in human carcinoma cell lines. Oncogene. 1998;17(21):2719–2732.
  • [110] Ma XJ, Salunga R, Tuggle JT, Gaudet J, Enright E, McQuary P, et al. Gene expression profiles of human breast cancer progression. Proceedings of the National Academy of Sciences. 2003;100(10):5974–5979.
  • [111] Escamilla-Hernandez R, Chakrabarti R, Romano RA, Smalley K, Zhu Q, Lai W, et al. Genome-wide search identifies Ccnd2 as a direct transcriptional target of Elf5 in mouse mammary gland. BMC molecular biology. 2010;11(1):68.
  • [112] Chakrabarti R, Hwang J, Blanco MA, Wei Y, Lukačišin M, Romano RA, et al. Elf5 inhibits the epithelial–mesenchymal transition in mammary gland development and breast cancer metastasis by transcriptionally repressing Snail2. Nature cell biology. 2012;14(11):1212–1222.
  • [113] Vernier M, Bourdeau V, Gaumont-Leclerc MF, Moiseeva O, Bégin V, Saad F, et al. Regulation of E2Fs and senescence by PML nuclear bodies. Genes & development. 2011;25(1):41–50.
  • [114] McConnell BB, Gregory FJ, Stott FJ, Hara E, Peters G. Induced expression of p16 INK4a inhibits both CDK4-and CDK2-associated kinase activity by reassortment of cyclin-CDK-inhibitor complexes. Molecular and cellular biology. 1999;19(3):1981–1989.
  • [115] Villacañas Ó, Pérez JJ, Rubio-Martínez J. Structural analysis of the inhibition of Cdk4 and Cdk6 by p16INK4a through molecular dynamics simulations. Journal of Biomolecular Structure and Dynamics. 2002;20(3):347–358.
  • [116] Bracken AP, Kleine-Kohlbrecher D, Dietrich N, Pasini D, Gargiulo G, Beekman C, et al. The Polycomb group proteins bind throughout the INK4A-ARF locus and are disassociated in senescent cells. Genes & development. 2007;21(5):525–530.
  • [117] Fang L, Igarashi M, Leung J, Sugrue MM, Lee SW, Aaronson SA. p21Waf1/Cip1/Sdi1 induces permanent growth arrest with markers of replicative senescence in human tumor cells lacking functional p53. Oncogene. 1999;18(18):2789–2797.
  • [118] Mao Z, Ke Z, Gorbunova V, Seluanov A. Replicatively senescent cells are arrested in G1 and G2 phases. Aging (Albany NY). 2012;4(6):431.
  • [119] Chellappan SP, Hiebert S, Mudryj M, Horowitz JM, Nevins JR. The E2F transcription factor is a cellular target for the RB protein. Cell. 1991;65(6):1053–1061.
  • [120] Byeon IJL, Li J, Ericson K, Selby TL, Tevelev A, Kim HJ, et al. Tumor Suppressor p16INK4A: Determination of Solution Structure and Analyses of Its Interaction with Cyclin-Dependent Kinase 4. Molecular cell. 1998;1(3):421–431.
  • [121] Beauséjour CM, Krtolica A, Galimi F, Narita M, Lowe SW, Yaswen P, et al. Reversal of human cellular senescence: roles of the p53 and p16 pathways. The EMBO journal. 2003;22(16):4212–4222.
  • [122] Freudlsperger C, Bian Y, Wise SC, Burnett J, Coupar J, Yang X, et al. TGF- and NF-B signal pathway cross-talk is mediated through TAK1 and SMAD7 in a subset of head and neck cancers. Oncogene. 2012;32(12):1549–1559.
  • [123] Harley C, Futcher A, Greider C. Telomeres shorten during ageing of human fibroblasts. Nature. 1990;345(6274):458–460.
  • [124] Mani SA, Yang J, Brooks M, Schwaninger G, Zhou A, Miura N, et al. Mesenchyme Forkhead 1 (FOXC2) plays a key role in metastasis and is associated with aggressive basal-like breast cancers. Proceedings of the National Academy of Sciences. 2007;104(24):10069–10074.
  • [125] Zeisberg M, Neilson EG, et al. Biomarkers for epithelial-mesenchymal transitions. The Journal of clinical investigation. 2009;119(6):1429–1437.
  • [126] Bolós V, Peinado H, Pérez-Moreno MA, Fraga MF, Esteller M, Cano A. The transcription factor Slug represses E-cadherin expression and induces epithelial to mesenchymal transitions: a comparison with Snail and E47 repressors. Journal of cell science. 2003;116(3):499–511.
  • [127] Dave N, Guaita-Esteruelas S, Gutarra S, Frias À, Beltran M, Peiró S, et al. Functional cooperation between Snail1 and twist in the regulation of ZEB1 expression during epithelial to mesenchymal transition. Journal of Biological Chemistry. 2011;286(14):12024–12032.
  • [128] Casas E, Kim J, Bendesky A, Ohno-Machado L, Wolfe CJ, Yang J. Snail2 is an essential mediator of Twist1-induced epithelial mesenchymal transition and metastasis. Cancer research. 2011;71(1):245–254.
  • [129] Hajra KM, Chen DY, Fearon ER. The SLUG zinc-finger protein represses E-cadherin in breast cancer. Cancer research. 2002;62(6):1613–1618.
  • [130] Davila-Velderrain J, Martinez-Garcia J, Alvarez-Buylla E. Descriptive vs. Mechanistic Network Models in Plant Development in the Post-Genomic Era. Plant Functional Genomics: Methods and Protocols. 2015;p. 455–479.
  • [131] Siegel R, Naishadham D, Jemal A. Cancer statistics, 2013. CA: a cancer journal for clinicians. 2013;63(1):11–30.
  • [132] Jemal A, Bray F, Center MM, Ferlay J, Ward E, Forman D. Global cancer statistics. CA: a cancer journal for clinicians. 2011;61(2):69–90.
  • [133] Stratton MR, Campbell PJ, Futreal PA. The cancer genome. Nature. 2009;458(7239):719–724.
  • [134] Hudson TJ, Anderson W, Aretz A, Barker AD, Bell C, Bernabé RR, et al. International network of cancer genome projects. Nature. 2010;464(7291):993–998.
  • [135] Weinstein JN, Collisson EA, Mills GB, Shaw KRM, Ozenberger BA, Ellrott K, et al. The cancer genome atlas pan-cancer analysis project. Nature genetics. 2013;45(10):1113–1120.
  • [136] Yaffe MB. The scientific drunk and the lamppost: massive sequencing efforts in cancer discovery and treatment. Science signaling. 2013;6(269):pe13.
  • [137] Creixell P, Schoof EM, Erler JT, Linding R. Navigating cancer network attractors for tumor-specific therapy. Nature biotechnology. 2012;30(9):842–848.
  • [138] DePinho RA, Polyak K. Cancer chromosomes in crisis. Nature genetics. 2004;36(9):932–934.
  • [139] Fujikawa M, Katagiri T, Tugores A, Nakamura Y, Ishikawa F. ESE-3, an Ets family transcription factor, is up-regulated in cellular senescence. Cancer science. 2007;98(9):1468–1475.
  • [140] Lee JM, Dedhar S, Kalluri R, Thompson EW. The epithelial–mesenchymal transition: new insights in signaling, development, and disease. The Journal of cell biology. 2006;172(7):973–981.
  • [141] Cano A, Pérez-Moreno MA, Rodrigo I, Locascio A, Blanco MJ, del Barrio MG, et al. The transcription factor snail controls epithelial–mesenchymal transitions by repressing E-cadherin expression. Nature cell biology. 2000;2(2):76–83.
  • [142] Sun Y, Song GD, Sun N, Chen JQ, Yang SS. Slug overexpression induces stemness and promotes hepatocellular carcinoma cell invasion and metastasis. Oncology Letters. 2014;7(6):1936–1940.
  • [143] Liu Y, El-Naggar S, Darling DS, Higashi Y, Dean DC. Zeb1 links epithelial-mesenchymal transition and cellular senescence. Development. 2008;135(3):579–588.
  • [144] Weinberg RA. Twisted epithelial–mesenchymal transition blocks senescence. Nature cell biology. 2008;10(9):1021–1023.
  • [145] Kim WY, Sharpless NE. The Regulation of¡ i¿ INK4¡/i¿/¡ i¿ ARF¡/i¿ in Cancer and Aging. Cell. 2006;127(2):265–275.
  • [146] Müssel C, Hopfensitz M, Kestler HA. BoolNet—an R package for generation, reconstruction and analysis of Boolean networks. Bioinformatics. 2010;26(10):1378–1380.
  • [147] Zhou JX, Samal A, d’Hèrouël AF, Price ND, Huang S. Relative Stability of Network States in Boolean Network Models of Gene Regulation in Development. arXiv preprint arXiv:14076117. 2014;.
  • [148] Wang G, Zhu X, Gu J, Ao P. Quantitative implementation of the endogenous molecular–cellular network hypothesis in hepatocellular carcinoma. Interface focus. 2014;4(3):20130064.
  • [149] Zhu X, Yuan R, Hood L, Ao P. Endogenous molecular-cellular hierarchical modeling of prostate carcinogenesis uncovers robust structure. Progress in biophysics and molecular biology. 2015;.
  • [150] Li C, Wang J. Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths. PLoS computational biology. 2013;9(8):e1003165.
  • [151] Zhou JX, Brusch L, Huang S. Predicting pancreas cell fate decisions and reprogramming with a hierarchical multi-attractor model. PloS one. 2011;6(3):e14752.
  • [152] Zhou JX, Qiu X, d’Herouel AF, Huang S. Discrete Gene Network Models for Understanding Multicellularity and Cell Reprogramming: From Network Structure to Attractor Landscapes Landscape. In: Computational Systems Biology Second Edition Elsevier. 2014;p. 241–276.
  • [153] Lander AD. The edges of understanding. BMC biology. 2010;8(1):40.
  • [154] Chen Z, Trotman LC, Shaffer D, Lin HK, Dotan ZA, Niki M, et al. Crucial role of p53-dependent cellular senescence in suppression of Pten-deficient tumorigenesis. Nature. 2005;436(7051):725–730.
  • [155] Campisi J. Cellular senescence: putting the paradoxes in perspective. Current opinion in genetics & development. 2011;21(1):107–112.
  • [156] Narita M, Lowe SW. Senescence comes of age. Nature medicine. 2005;11(9):920–922.
  • [157] Yildiz G, Arslan-Ergul A, Bagislar S, Konu O, Yuzugullu H, Gursoy-Yuzugullu O, et al. Genome-wide transcriptional reorganization associated with senescence-to-immortality switch during human hepatocellular carcinogenesis. PloS one. 2013;8(5):e64016.
  • [158] Smit MA, Peeper DS. Epithelial-mesenchymal transition and senescence: two cancer-related processes are crossing paths. Aging (Albany NY). 2010;2(10):735.
  • [159] Glaab E, Baudot A, Krasnogor N, Schneider R, Valencia A. EnrichNet: network-based gene set enrichment analysis. Bioinformatics. 2012;28(18):i451–i457.
  • [160] Glaab E, Baudot A, Krasnogor N, Valencia A. TopoGSA: network topological gene set analysis. Bioinformatics. 2010;26(9):1271–1272.
  • [161] Winterhalter C, Widera P, Krasnogor N. JEPETTO: a Cytoscape plugin for gene set enrichment and topological analysis based on interaction networks. Bioinformatics. 2014;30(7):1029–1030.
  • [162] Saadatpour A, Albert R, Reluga TC. A reduction method for Boolean network models proven to conserve attractors. SIAM Journal on Applied Dynamical Systems. 2013;12(4):1997–2011.
  • [163] Naldi A, Remy E, Thieffry D, Chaouiya C. Dynamically consistent reduction of logical regulatory graphs. Theoretical Computer Science. 2011;412(21):2207–2218.
  • [164] Allen LJ. An introduction to stochastic processes with applications to biology. CRC Press; 2010.
  • [165] Li J, Poi MJ, Tsai MD. Regulatory mechanisms of tumor suppressor P16INK4A and their relevance to cancer. Biochemistry. 2011;50(25):5566–5582.
  • [166] Yamakoshi K, Takahashi A, Hirota F, Nakayama R, Ishimaru N, Kubo Y, et al. Real-time in vivo imaging of p16Ink4a reveals cross talk with p53. The Journal of cell biology. 2009;186(3):393–407.
Figure 1: Gene regulatory network for epithelial carcinogenesis.
Figure 2: Core gene regulatory network module for epithelial carcinogenesis.
Figure 3: The core gene regulatory module in the context of the Hallmarks of Cancer approach.
Figure 4: Temporal sequence and global order of cell–fate attainment pattern under the stochastic Boolean GRN model during epithelial carcinogenesis.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
354603
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description