Skip to main content

Gene regulatory network underlying the immortalization of epithelial cells



Tumorigenic transformation of human epithelial cells in vitro has been described experimentally as the potential result of spontaneous immortalization. This process is characterized by a series of cell–state transitions, in which normal epithelial cells acquire first a senescent state which is later surpassed to attain a mesenchymal stem–like phenotype with a potentially tumorigenic behavior. In this paper we aim to provide a system–level mechanistic explanation to the emergence of these cell types, and to the time–ordered transition patterns that are common to neoplasias of epithelial origin. To this end, we first integrate published functional and well–curated molecular data of the components and interactions that have been found to be involved in such cell states and transitions into a network of 41 molecular components. We then reduce this initial network by removing simple mediators (i.e., linear pathways), and formalize the resulting regulatory core into logical rules that govern the dynamics of each of the network components as a function of the states of its regulators.


Computational dynamic analysis shows that our proposed Gene Regulatory Network model recovers exactly three attractors, each of them defined by a specific gene expression profile that corresponds to the epithelial, senescent, and mesenchymal stem–like cellular phenotypes, respectively. We show that although a mesenchymal stem–like state can be attained even under unperturbed physiological conditions, the likelihood of converging to this state is increased when pro–inflammatory conditions are simulated, providing a systems–level mechanistic explanation for the carcinogenic role of chronic inflammatory conditions observed in the clinic. We also found that the regulatory core yields an epigenetic landscape that restricts temporal patterns of progression between the steady states, such that recovered patterns resemble the time–ordered transitions observed during the spontaneous immortalization of epithelial cells, both in vivo and in vitro.


Our study strongly suggests that the in vitro tumorigenic transformation of epithelial cells, which strongly correlates with the patterns observed during the pathological progression of epithelial carcinogenesis in vivo, emerges from underlying regulatory networks involved in epithelial trans–differentiation during development.


Nearly 84% of cancers diagnosed in human adults are carcinomas (i.e., cancer of epithelial origin). Although epithelial carcinogenesis has been strongly associated with a chronic inflammatory process and aging [1], the precise role of these two processes to the origin and progression of carcinomas remains elusive. The current general assumption is that aging and inflammation increase the chance of accumulating somatic mutations, and that this genetic instability constitutes the cause of carcinoma. But this view does not explain several well–described experimental and clinical observations. For instance, cancer behavior can be acquired in the absence of mutations through trans–or de–differentiation and is characterized for recapitulating embryonic processes. Cancer cells can be “normalized” by several experimental means and commonly show morphological and transcriptional convergence despite their diverse origin and mutations [24]. In addition, different carcinomas share similar time–ordered patterns of progression, as well as clear associations with lifestyle factors in many cases [5]. These empirical observations suggest that, in analogy to normal development, epithelial carcinogenesis is a consequence of conserved or generic system–level mechanisms that restrict the malignant phenotypes that can be attained, as well as the temporal patterns of progression that describe the transitions between them. In accordance with this latter view, it has been proposed that cancer can be considered a developmental disease [6].

In agreement with this developmental view of cancer, numerous experimental findings in molecular and cell biology of cancer research have revealed that it is possible to recover cells with cancer–like phenotypes through the induction of de–differentiation events. This has been shown particularly in carcinomas [3, 7], since inflammatory cytokines induce a de–differentiation event of epithelial cells denominated Epithelial–Mesenchymal Transition (EMT) in which cells acquire a mesenchymal stem–like phenotype with a tumorigenic potential able to develop cancer in mice [3].

In systems biology it is common to study cell differentiation processes that underlie development and pathogenesis from the point of view of dynamical systems theory. In this framework, the information encoded by the genome, in addition to epigenetic mechanisms, can be mapped to a gene regulatory network (GRN), that shows multiple stable steady states, each of which corresponds to a particular phenotypic cellular state. Further, the GRN also underlies the epigenetic landscape (EL), that restricts the time–ordered patterns of transition between the phenotypes [813]. Thus, the same genome and GRN can robustly generate multiple discrete cellular phenotypes through developmental dynamics [11, 14, 15]. These stable phenotypic states are called attractors and correspond to configurations of gene or protein activation states that underlie the cellular fates or phenotypes. Therefore, dynamic developmental processes – particularly, cellular differentiation and morphogenetic patterning – can be formalized in temporal terms as transitions between attractors (i.e., cell states). Here, we adopt such approach to study how tumorigenic transformation due to spontaneous immortalization via EMT emerges from the regulatory interactions between different molecular players that are known to contribute to the tumorigenic transformation of epithelial cells. We hypothesize that: (1) a generic series of cell state transitions describing the phenotypic transformation of epithelial cells first to senescent cells and finally to mesenchymal stem–like cells, that is widely observed in epithelial cell cultures during spontaneous immortalization, naturally result from the self–organized behavior emerging from an underlying intracellular GRN; and (2) that pro–inflammatory tissue–level conditions, which are associated with a bad prognosis, increase the likelihood of this tumorigenic process, promoting the emergence and progression of epithelial cancer.

To test our hypothesis, we propose here a cellular–level GRN that, for the first time, integrates those molecular components and their interactions that have been experimentally shown to play an important role during the emergence and progression of carcinomas. It includes the key molecular regulators that: (1) characterize the cellular phenotypes of epithelial, mesenchymal stem–like, and senescent cells; (2) are involved in the induction of the cellular processes of replicative senescence, cellular inflammation, and epithelial–mesenchymal transition (EMT); and (3) are involved in the phenotypic changes undergone by cells emerging from these processes (i.e., mesenchymal stem–like cells). We then obtained a reduced regulatory core for further dynamical analyses by removing linear cascades while maintaining the feedback loops. 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 modeling approaches (see, for example [8, 9, 16, 17]). Specifically, by proposing logical functions grounded on the available experimental data for this regulatory core module, and by analyzing its behavior following conventional Boolean GRN dynamical approaches, we show that the uncovered GRN converges to exactly the three attractors corresponding to the expected gene expression configurations characterizing the epithelial, senescent, and mesenchymal stem-like phenotypes. Additionally, using a stochastic version of the model to explore the GRN EL (following the methodology proposed in [13]), we found that the proposed GRN also explains the commonly observed temporal sequence by which epithelial cells acquire the potentially tumorigenic mesenchymal stem–like phenotype.

Our results suggest that the proposed core GRN incorporates a set of necessary and sufficient components and interactions to explain the emergence of gene configurations characteristic of epithelial, senescent and mesenchymal cells, as well as the time–ordered sequence of cellular transformations observed in the spontaneous immortalization process that, in turn, underlies the tumorigenic transformation of epithelial cells.


Gene regulatory network underlying spontaneous immortalization

Following a bottom–up approach, we performed an extensive literature search to gather the most relevant experimental functional molecular data describing the cellular–level processes involved in epithelial carcinogenesis, namely: replicative cellular senescence, inflammation, and EMT (see Additional file 1: Table S4). We found a set of 41 molecular players (12 TFs and 29 signaling molecules) which are involved in epithelial or mesenchymal cell differentiation, cellular inflammation, senescence, DNA damage, cell cycle, or in epigenetic silencing; as well as 97 regulatory interactions between them. For the first time, we integrated this previously scattered experimental information into the GRN represented in Fig. 1 a. To further support that the set of regulatory interactions that we manually curated based on published data are indeed representative of the cellular–level processes underlying epithelial carcinogenesis, we performed a network–based Gene Set Enrichment Analysis (GSEA) of the GRN, using both the KEGG and the GO Biological Process databases as reference. We found that among the 12 pathways or processes reported as significant when taking the KEGG database as a reference, 10 (≈ 83%) correspond to the cancer types bladder cancer, chronic myeloid leukemia, non-small cell lung cancer, glioma, melanoma, pancreatic cancer, prostate cancer, small cell lung cancer, thyroid cancer, from which 6 (66.6%) correspond to carcinomas. When taking the GO Biological Process database as reference, we found that the molecules considered in our regulatory network are significantly enriched for several of the biological processes known to play important roles during spontaneous immortalization of epithelial cells, namely replicative senescence, cellular senescence, cell aging, positive regulation of ephitelial to mesenchymal transition determination of adult life span, among others (Table 1). Additionally to these GSEA, we performed a Network-based topological gene set enrichment analysis (described in Methods) and found that, in addition to the enrichment of the pathways and processes described above, the molecules in the proposed network show also a topological signature that strongly resembles the structure of the cancer pathways included in the KEGG database (see Additional file 2: Figure S1). The complete enrichment results are included in Additional file 3: Supplementary Tables. These results provide further support for the relevance of the proposed molecular players, manually assembled from an exhaustive literature search, and for the novel regulatory module proposed here.

Fig. 1
figure 1

Gene Regulatory Network underlying spontaneous immortalization. a Gene regulatory network for epithelial carcinogenesis. b A Core Regulatory Network Module Underlying Spontaneous Immortalization and Epithelial–Mesenchymal Transition. c Predicted gene expression profiles characterizing the epithelial, senescent and mesenchymal stem–like cells. d Cellular inflammation increases the size of the basin of attraction of the mesenchymal stem–like phenotype

Table 1 Most significant pathways and processes in the GRN (Fig. 1 a) shown by a network–based gene set enrichment analysis

In conclusion, based on these analyses and the current state of knowledge according to annotated databases, the set of molecules manually included in the proposed large network seem 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 network are tightly associated with reference pathways of epithelial cancers. Additionally, robustness analysis performed on the network showed that the recovered attractors are also robust to permanent alterations of the regulatory logic (see Additional file 4: Text).

To study the dynamic and steady–state properties of the network underlying epithelial carcinogenesis, we reduced the initial network (Fig. 1 a) into a more mathematically and computationally tractable regulatory core network by collapsing the linear pathways while retaining the feedback loops. This reduced network, represented in Fig. 1 b, retains only the main players and regulatory interactions that underlie epithelial and mesenchymal differentiation, cell cycle progression, senescence, and cellular inflammation.

Epithelial differentiation is represented by the tumor suppressor Epithelium-Specific TF ESE–2, which acts as a master regulator of this process by triggering the expression of epithelium-specific genes while repressing mesenchymal markers, such as Snail [18]. ESE–2 also promotes its own expression and the expression of the other ESE TFs [19], and it is hence considered as the main representative of the TF family [18]. The differentiation of epithelial cells into a mesenchymal stem–like phenotype is orchestrated by the TF Snai2, which further induces its own expression via the activation of Twist-Zeb-FOXC2 [20]. The progression of the cell cycle is controlled by Rb, E2F, Cyclin, and Telomerase (here TELase). While both E2F [21] and cyclins [22] are required for cell cycle progression, the tumor suppressor Rb acts as an inhibitor of this process by forming an inactive heteromer with E2F [23]. The enzyme TELase is responsible for the de novo synthesis of telomeres, a process that allows cells to surpass the cell cycle checkpoints and become immortalized [24]. Indeed, high levels of TELase are typical of tumor initiating cells that become resistant to therapy [3, 25]. The senescence of epithelial cells is mainly controlled by the two tumor supressors p53 and p16 [26]. They both induce replicative senescence by reducing the activity of cyclins [27] and by inducing an Rb-mediated inhibition of E2F [28]. Cellular inflammation is characterized by the increased activity of the TF NF– κB. Many (micro)environmental stimuli, including pathogens, cytokines, interleukins, and antigens; trigger the expression of NF– κB, resulting in an immune response characterized by increased levels of cytokines and enzymes such as phospholipase A2, cyclooxygenase, and lipoxygenase [29].

The above described molecular players (ESE–2, Snai2, Rb, E2F, cyclin, TELase, p53, p16, and NF– κB) are tightly interconnected by the following regulatory interactions: NF– κ B positively controls its own expression via the induction of pro-inflammatory cytokines [30, 31]. It also promotes the epithelial and mesenchymal stem–like phenotype by positive interactions with ESE–2 [32] and Snai2 [33], respectively; and it induces the cell cycle process by increasing cyclin [34] and inhibiting p53 expression [35]. ESE–2 forms an auto–activation feedback loop, inhibits Snai2 [18], induces cyclin, inhibits TELase [18], and activates of NF– κB [36]. Snai2 induces its own expression [37], increases the activity of TELase [38], and decrease the expression of the epithelial–specific TF ESE–2 [18]. It also represses the cell cycle process by reducing the transcription of cyclin [39] and E2F [40], represses cellular senescence by decreasing the expression of p16 [41] and p53 [42], and induces cellular inflammation by increasing the expression of NF– κB [43]. p16 promotes its own activity forming a positive feedback loop [44], stabilizes p53 via the inactivation of its inhibitor MDM2 [45], and inhibits the cell cycle progression by the induction of Rb [46] and the inhibition of cyclin [23, 28]. These molecular processes are known to induce the senescent phenotype. p16 also contributes to the cellular inflammation observed in senescent cells by inducing the activation of NF– κB [47]. p53 inhibits E2F [48]. Rb inhibits E2F [23]. E2F induces cyclin [49] and inhibits the senescence marker p16 [50]. Cyclin inhibits Rb activity [51] and stimulates E2F transcription [52]. Finally, TELase inhibits the senescence markers p16 [53] and p53 [54].

To analyze its long term and dynamic properties, we then formalized the above described regulatory interactions (Fig. 1 b), by translating the nodes and their corresponding logical interaction rules into a mechanistic Boolean dynamical GRN model [55]. The corresponding logical rules and truth tables can be found in the Additional file 5: Logic rules and truth tables.

The epithelial, senescent, and mesenchymal stem–like phenotypes are attractors of the GRN

It has been experimentally shown that during the process of spontaneous immortalization epithelial cells acquire first a senescent and finally a mesenchymal stem–like cell phenotype [56], each of which can be described in terms of their gene expression profiles. Here, we performed long-term (technical term) simulations of our Boolean model of the GRN with the aim to recover these three characteristic expression patterns. We found that indeed our GRN (Fig. 1 b) converges to three attractors, each corresponding, respectively, to the epithelial, senescent, and mesenchymal stem–like phenotypic markers (Fig. 1 c).

Epithelial cells are characterized by the high expression of the TF ESE–2, which acts as a master regulator. Being part of a constantly renewing tissue, epithelial cells show an increased expression of the cell cycle inducers Cyclin and E2F and a decreased expression of the cell cycle inhibitor Rb [57] and of the senescence markers p53 and p16 [58]. Epithelial cells do not express TELase, since the activity of this enzyme is inhibited in response to induction of differentiation events [59]. Being a constitutively expressed gene, epithelial cells show also high levels of the TF NF– κB [60]. In accordance to this empirical activity profile, approximated as a Boolean vector, the first attractor of our GRN shows an absence of activity for p53, p16, Rb, TELase and Snai2; and presence of NF– κB, Cyclin, E2F and ESE–2 (Fig. 1 c, left).

Senescent epithelial cells conserve the high levels of the epithelial markers ESE–2 [61] and the low expression of the mesenchymal marker Snai2 [18]; but, in contrast to normal epithelial cells, these cells show an increased expression of the two tumor suppressor proteins p16 and p53 [62], both of which contribute to the cell cycle arrest by inhibiting the cell cycle inducers Cyclins and E2F [27, 28, 63] and by inducing the cell cycle suppressor Rb [64]. As in homeostatically cycling epithelial cells, TELase is inactive in this cell type [59]. Also in this case, NF– κB expression is constitutively active. These features are recovered in our second attractor, which shows activity of NF– κB, p53, p16, Rb and ESE2, and absence of TELase, Snai2, cyclin, and E2F (Fig. 1 c, center). Also in this case, our in silico analysis predicts activity of NF– κB in epithelial senescent cells.

Mesenchymal stem–like cells express the mesenchymal marker Snai2 [18], which acts as one of the key players during the EMT by orchestrating the repression of epithelium–specific genes [65] such as ESE–2 [18] and by inducing the expression of mesenchymal markers [63]. In contrast to the senescent cells, these mesenchymal stem-like cells have a strong proliferative potential, shown by a decreased expression of p53, p16, Rb; and by a high expression of cyclin [66]. Cell cycle progression is further promoted in these potentially tumor initiating cells by the high levels of TELase [3, 25]. This cell type also shows a constitutive expression of NF– κB [60, 67]. Our third attractor recovers this genetic configuration, showing inactivity of p53, p16, Rb, E2F and ESE–2; and activity of NF– κB, TELase, Snai2, and Cyclin (Fig. 1 c, right).

Cells that have acquired this latter mesenchymal stem–like phenotype have a strong tumorigenic potential, since they show most of the hallmarks of cancer: The low expression of ESE–2 accompanied by the over–activation of Snai2 enable cells to sustain proliferative signals and to evade growth suppressors by undergoing a de–differentiation process [68]. Further, the constitutive activity of Snai2 is associated to an motile and invasive phenotype, and with the avoidance of immune destruction [69] and deregulation of cellular energetics [70]. The de–activation of the senescence markers p16 and p53 and of the tumor suppressor Rb, as well as the over–expression of TELase, confers genome instability that is associated to a mutation–prone phenotype [71] and enable cells to acquire replicative immortality and to resist cell death. Moreover, high levels NF– κB suggest chronic and tumor-promoting cellular inflammation [72].

The correspondence between the recovered attractors in the GRN model simulation and the experimentally observed gene or protein configurations in the studied cellular phenotypes strongly suggests that the proposed core GRN indeed constitutes a regulatory module that comprises a set of components and interactions able to restrict the system to converge to the cellular states observed during spontaneous immortalization. We conclude that the derived core GRN module (Fig. 1 b) includes a set of sufficient and necessary molecular players and interactions that specify epithelial, senescent and stem–like mesenchymal cells.

Cellular inflammation accelerates the spontaneous immortalization of epithelial cells

Inflammation has been recognized as one of the key drivers of carcinogenesis, partly due to its implication in the EMT [79]. Indeed, several pro-inflammatory cytokines such as TFG- β [80] and IL-6 [81], some of which are produced by senescent cells [82], have shown to be strong inducers of EMT. Cells that are exposed to such a pro-inflammatory micro-environment show an over-activation of NF– κB [31, 47], which induces EMT by triggering the expression of mesenchymal markers including Snail [20] and silences the expression of p16 and p53 [83]. To assess whether our model reproduces this increased propensity of de–differentiation into a mesenchymal stem–like phenotype in the presence of cellular inflammation, we calculated the relative size of the basins of attraction of the epithelial, senescent, and mesenchymal stem–like phenotypes (Fig. 1 c), with and without the assumption of a constant over–activation of the NF– κB node in the GRN model. We found that cellular inflammation increases the size of the mesenchymal stem–like attractor from 56.25 to 75% while decreasing the region of convergence of the epithelial (from 17.97 to 6.25%) and of the senescent (from 25.78 to 18.75%) phenotypes (Fig. 1 d). These results are in accordance with the experimental results stated above, in which cellular inflammation is recognized as an important driver of spontaneous immortalization by promoting the induction of de–differentiation of epithelial cells into a mesenchymal stem–like phenotype.

The proposed GRN reproduces the characteristic phenotypes of 6 different mutant conditions

To further validate our model, we tested if the GRN (Fig. 1 b) is able to reproduce the phenotypic configurations (in the form of attractors) observed in 6 different mutant conditions. Specifically, we simulated loss– and gain–of–function conditions of ESE–2, Snai2, and p16 (by setting the expression state for the corresponding node constitutively to “1” or “0”, respectively), and compared the resulting attractors to the corresponding gene expression profiles reported in the literature. When simulating ESE–2 loss–of–function the model recovers a single attractor, equivalent to a mesenchymal stem–like phenotype with increased Snai2 expression as experimentally observed [18] (Fig. 2 a). Simulations of the ESE–2 gain of function mutant results in three attractors that are consistent with the experimental description of the phenotypes resulting from ESE–2 over–expression: an epithelial senescent cell [61], a normal epithelial cell [18], and a EMT intermediate state presenting simultaneous features of both epithelial and mesenchymal states with proliferative phenotype [73] (Fig. 2 b). Simulating Snai2 loss–of–function gives rise to two attractors corresponding to normal and senescent epithelial phenotypes, both consistent with experimental observations [18] (Fig. 2 c). Simulating Snai2 gain–of–function mutation results in a single attractor corresponding to mesenchymal stem–like phenotype, which is also consistent with experimental observations derived from ectopic over–expression of mesenchymal TFs [74] (Fig. 2 d). These results are also consistent with the TGF- β-dependent induction of EMT, which occurs via the activation of Snail by downstream components of the TGF- β signaling pathway [75], and which have been successfully reproduced in a recent model of TGF- β-driven EMT [76]. p16 loss–of–function simulations recover two attractors corresponding to an epithelial and a mesenchymal stem–like cell, which is also consistent with experimental observations [77] (Fig. 2 e). 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 immortalization and apoptosis–resistance shown by mesenchymal stem–like cells, as well as with the capability of mesenchymal TFs to abrogate senescence [78]. The second attractor is consistent with the potential for replicative senescence of epithelial cells [56] (Fig. 2 f). In conclusion, we found that our model simulations are consistent with six mutant conditions reported in the literature, providing further validation for our proposed GRN.

Fig. 2
figure 2

Predicted attractors of loss – and gain of function mutants of the GRN for ESE2 (a, b), Snai2 (c, d) and p16 (e, f). Percent (%) represent the size of the corresponding basin of attraction

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 can be robustly induced [3]. Normal epithelial cells become first senescent cells, a state which they afterwards overcome to acquire a 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, we performed two independent EL: (1) We explored the temporal sequence of attractor attainment, and (2) calculated the consistent global ordering of all the given attractors. To do so, we followed [13, 84] and explored the EL associated to the GRN by implementing a discrete stochastic model as an extension to the deterministic Boolean network model [11] (detailed in Methods). The results of our analysis regarding the temporal sequence of attractor attainment (following the methodology proposed in [13]) show that the most probable temporal order of attractor attainment for a population of cells that initially show and epithelial phenotype correspond to the expected transition from an epithelial to a senescent to a mesenchymal stem–like phenotype (Fig. 3 a). The estimated transition probability matrices are given in Additional file 4: Table S1.

Fig. 3
figure 3

Temporal sequence and global order of cell–fate attainment pattern under the stochastic Boolean GRN model during epithelial carcinogenesis. a Maximum probability p of attaining each attractor, as a function of time (in iteration steps). The most probable sequence of cell attainment is: epithelial → senescent → mesenchymal stem-like. The error probability used in this case was ξ=0.05. 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. c Schematic representation of the time–ordered phenotype transitions along the epigenetic landscape, showing the in–between attractor barrier highs in the landscape

In agreement with these results, following [85] we found a consistent global temporal ordering of the uncovered attractors. This analysis is based on calculating the relative stability of the three different attractors, which is done by computing the Mean First Passage Time (MFPT, shown in Additional file 4: Table S2) 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. 3 b). This path corresponds to the only order in which the system can visit the three attractors following a positive net transition rate (Additional file 4: Table S3). These results indicate that, when considering 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 (Fig. 3 c). We anticipate that observed transition rates in vivo are likely to be altered or favored by 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 single cell and multi–level models building up on the framework that has been put forward here.


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 coherent theory to explain the development of multicellular organisms as the result of the orchestrating role of GRNs has been developed [8, 10, 11]. The main conclusion is that observable cell states emerge from the self–consistent multi–stable regulatory logic dictated by genome structure and obeyed by (mainly) 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 [2]. Might also cancer be orchestrated by GRN dynamics? Several hypothesis have been presented in this direction such as the cancer attractor theory [2, 6], and the endogenous molecular cellular network hypothesis [86, 87]. In this contribution we adhere to the viewpoint of an intrinsic regulatory network, but 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 cellular inflammation and resulting from surpassing a senescent state through EMT – i.e., tumorigenic transformation due to spontaneous immortalization via EMT [63, 88]. We propose that a mechanistic understanding of this process is an important first 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 GRN underlying the epithelial, senescent and mesenchymal stem–like phenotypes

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 curated literature into a GRN comprising the main molecular regulators involved in key cellular processes ubiquitous to carcinogenesis, following a bottom–up approach. Subsequently, we followed a mechanistic approach to address the question of whether we assembled a set of sufficient and necessary molecular players and interactions able to recover the cellular phenotypes and processes documented during the spontaneous immortalization of human epithelial cells in culture. Thus, in this contribution, we 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 [89], cell–fate decision [90] and similar cellular processes during plant morphogenesis [8, 9, 13, 16]. 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 (Fig. 1 a) into a regulatory core module, consisting of a small set of key molecular players and the regulatory interactions between them (Fig. 1 b). Although the components of the core module have been previously shown to be involved in EMT, the proposed architecture and topology of this core module had not been proposed before. We extracted from available literature the expression profiles of the generally observable cell states of interest in terms of this minimal set of molecules, and tested whether the reduced GRN, formalized into logical rules, is able to recover the biologically observable expression profiles as stationary and stable network configurations (i.e., attractors). Indeed, we found that the core GRN model converges exactly to the three observed gene expression profiles that correspond to the three phenotypes that wild type cells acquire during the process of spontaneous immortalization (Fig. 1 c). Our model analysis also shows that cellular inflammation increases the size of the basin of attraction of the potentially tumorigenic mesenchymal stem–like phenotypic attractor (Fig. 1 d), which is in agreement with the strong tumorigenic effect of inflammation that has been consistently reported in the experimental and clinical literature. Further, our in silico simulations of mutant conditions are also congruent with the corresponding expression profiles (Fig. 2). These results strongly suggest that we have successfully included the key regulators and interactions at play during the establishment of cell steady–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 [9, 16, 90]. Our GRN model exclusively integrates available published functional 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 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.

Time–ordered transitions of the phenotypic attractors

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 [11, 17, 84] and/or characterize [85] 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 multidimensional 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 path across the epigenetic landscape [91]: ordered transitions sharing positive net transition rates will be preferentially followed. Our results show that such a developmental path follows the time–order of cellular phenotypic states from epithelial, to senescent and finally mesenchymal stem–like cells with a tumorigenic potential (Fig. 3). In other words, the constraints imposed by the GRN structure the associated EL in such a way that an epithelial cell in culture “ball” would naturally roll following such a path, in agreement with the observed spontaneous immortalization process (Fig. 3 c).

It is interesting that 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 which correlate with carcinoma progression in vivo. This suggests that specific progressive alterations or particular “abnormal” signaling mechanisms are not necessarily required for a cell to reach a potentially tumorigenic state.


In this contribution we propose an experimentally grounded GRN model for spontaneous immortalization. For the first time, we integrated a wealth of empirical evidence into a GRN (41 nodes, Fig. 1 a) which we reduced to a core GRN developmental module (9 nodes, Fig. 1 b) that converges to the three phenotypes observed during the spontaneous immortalization of epithelial cells (Fig. 1 c). Simulations of cellular inflammation lead to an increase in the size of the basin of attraction of the potentially tumorigenic mesenchymal stem–like phenotype (Fig. 1 d). Our model also reproduces several experimentally reported mutant conditions (Fig. 2), as well as the time–ordered phenotypic transitions undergone by epithelial cells during the process of spontaneous immortalization (Fig. 3). The proposed GRN constitutes thus a integrative, coherent, and experimentally validated framework which can be used in the future for the integration and analysis of additional signaling and mechanical processes that affect, and are affected by, the oncogenic transformation of the epithelial tissue. To this end, first we will explicitly incorporate into the current model EMT–inducing signaling pathways, (triggered by pro-inflammatory cytokines), increased TFG- β concentrations, or changes in the mechanical properties of the surrounding environment. Analysing how these pathways affect the phenotypic outcome, and how different nutritional or pharmacological conditions modulate these trans-differentiation processes (using for example the methodologies proposed in [92, 93]) can aid the design of better prevention and treatment or even oncogenic transformation reversion strategies. To explore this complex interplay between the transcriptional and signaling events driving the phenotypic transitions in individual cells, and the changing properties of the epithelial tissue, we will develop multi-scale tissue level models. Such a framework will enable the systematic assessment of the role on the oncogenic process of several tissue–level constraints such as cell cycle progression, cell–cell interactions, differential proliferation rates, as well as physical fields (e.g., mechanical forces).

Our proposed integrative, quantitative, and dynamical framework contributes to the understanding of the mechanisms underlying the onset and progression of epithelial cancer, which is necessary for devising new and more effective preventive strategies that halt or slow the progression of these diseases.


Construction of the network underlying epithelial carcinogenesis

A total of 159 references, considering both references in the main text and in the Supplementary Information [see Additional file 1: Table S4], were carefully and manually reviewed in order to determine a minimal set of cellular phenotypes and processes which enable a generic representation of spontaneous immortalization (also reported in epithelial carcinogenesis) on the basis of cell state transition events.

Following a bottom–up and an expert knowledge approach we propose a set of cellular dynamical processes ubiquitous to spontaneous immortalization (also reported in epithelial carcinogenesis), namely: replicative cellular senescence EMT driven by inflammation. The cellular phenotypes, epithelial, senescent, and mesenchymal stem–like produced by EMT induction have been largely characterized as biological observables involved in such processes. We take this information as a methodological basis to integrate a generic dynamical network model of spontaneous immortalization. As a first step in network integration, we assembled a set of TFs and additional molecular species involved in the establishment and regulation of these cellular states and processes. Subsequently, we manually retrieved documented regulatory interactions among the molecular species, considering only those supported by experimental evidence. The resulting nodes and regulatory interactions were then assembled manually into the network represented in Fig. 1 a with nodes representing genes and proteins and edges representing activating or inhibitory interactions.

Network–based gene set enrichment analysis

The bioinformatics tools EnrichNet [94] and TopoGSA [95] 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 first maps the input gene set into a reference network, to then compute its topological statistics, to finally compare these against the topology of pathways/processes in a reference database. Here a connected human interactome graph extracted from the STRING database, the KEGG and the GO Biological Process databases were used as reference molecular interaction network and database, respectively. Both analysis were performed using the Cytoscape plugin Jepetto [96].

Network reduction

The regulatory core underlying spontaneous immortalization (Fig. 1 b) was obtained from the large network (Fig. 1 a) by iteratively reducing all the simple mediator (i.e. those with in-degree and out-degree of one) and source (those with no regulators) nodes. This reduction process has been shown to conserve the attractors of the original Boolean network under an asynchronous update method [97]). During the reduction process, we kept the main molecular players controlling cellular senescence (p53 and p16), Cell Cycle (TELase, E2F, Rb and cyclin), differentiation into epithelial cells (ESE-2), differentiation into mesenchymal cells (Snai2), and cellular inflammation (NF- κB), since different activity configurations of these molecular players define the transcriptional identity of the three cellular phenotypes we sought to reproduce with the proposed model.

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:

$$ x_{i}(t+1) = F_{i}\left(x_{1}(t),x_{2}(t),\ldots,x_{k}(t)\right), $$

where the set of functions F i are logical propositions (or truth tables) expressing the relationship between the genes that share regulatory interactions with the gene i, and where the state variables x i (t) can take the discrete values 1 or 0 indicating whether the gene i is expressed or not at a certain discrete–time t, 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 [17]. The set of logical functions for the core regulatory module used in this study is given in the Additional file 5. The dynamical analysis of the Boolean network model was conducted using the package BoolNet [98] within the R statistical programming environment (

Size of the basins of attraction

The size of the basins of attraction of the GRN corresponds to the percentage of initial Boolean network configurations converging to that attractor, and were calculated for our GRN under wt (Fig. 1 d) and under different mutant conditions (Fig. 2) by exhaustively exploring all the 29=512 possible initial configurations of the GRN.

Epigenetic Landscape characterization

Stochastic version of the Boolean GRN

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 model was implemented following [13, 17, 84]. 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:

$$\begin{array}{@{}rcl@{}} & P_{x_{i}(t+1)}\left[F_{i}\left(\mathbf{x}_{reg_{i}}(t)\right)\right] = 1- \xi, \end{array} $$
$$\begin{array}{@{}rcl@{}} & P_{x_{i}(t+1)}\left[1 - F_{i}\left(\mathbf{x}_{reg_{i}}(t)\right)\right] = \xi. \end{array} $$

The probability that the value of the now random variable x i (t+1) is determined or not by its associated logical function \(\phantom {\dot {i}\!}F_{i}(\mathbf {x}_{reg_{i}}(t))\) is 1−ξ or ξ, respectively.

Attractor Transition Probability Estimation

An attractor transition probability matrix Π with components:

$$ \pi_{ij} = P\left(A_{t+1}=j|A_{t}=i\right), $$

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

Following the Discrete Time Markov Chains [99] theoretical framework, the estimated transition probability matrix was integrated into a dynamic equation for the probability distribution:

$$ P_{A}(t+1) = \Pi P_{A}(t), $$

where P A (t) is the probability distribution over the attractors at time t, and Π is the transition probability matrix (Additional file 4: Table S1). This equation was then iterated to simulate the temporal evolution of the probability distribution over the attractors starting from a specific initial probability distribution (Fig. 3 a).

Attractor relative stability and global ordering analyses

In addition to the calculation of the most probable temporal cell–fate pattern [13], a discrete stochastic GRN model enables the study of the ease for transitioning from one attractor to another [91]. 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 stability of all the attractors of a Boolean GRN [85, 91]. Here, the relative stability of each one of the cell states were defined based on the MFPT. Specifically, a relative stability matrix M was calculated which reflects the transition barrier between any two states based on the MFPT. Here, the MFPT was numerically estimated. Using the transition probabilities among attractors, a large number sample paths of a finite Markov chain were simulated. The MFPT from attractor i to attractor j corresponds to the averaged value of the number of steps taken to visit attractor j for the fist time, given that the entire probability mass was initially localized at the attractor i. The average is taken over the realizations. Following [91], based on the MFPT values a net transition rate between attractor i and j can be defined as follows:

$$ d_{i,j} = \frac{1}{MFPT_{i,j}} - \frac{1}{MFPT_{j,i}}. $$

This quantity effectively measures the ease of transition as a net probability flow, and is shown for our GRN in Additional file 4: Table S2. For all the calculation involving stochasticity, the robustness of the results was assessed by taking three different values for the probability of error (0.01,0.05,0.1). 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 [85]. Briefly, the consistent global ordering of the attractors is given by the attractor permutation in which all transitory net transition rates from an initial attractor to a final attractor are positive. This is schematically represented in Fig. 3 b. Calculated transition probabilities, MFPT, and net transition rate matrices are included in Additional file 4: Table S1, S2, and S3, respectively. R source code implementing all the calculations and analyses is available upon request.



Epigenetic landscape


Epithelial-Mesenchymal transition


Gene set enrichment analysis


Gene regulatory network


Mean first passage time


  1. Anand P, Kunnumakara AB, Sundaram C, Harikumar KB, Tharakan ST, Lai OS, Sung B, Aggarwal BB. Cancer is a preventable disease that requires major lifestyle changes. Pharm Res. 2008; 25(9):2097–116.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Huang S. On the intrinsic inevitability of cancer: from foetal to fatal attraction. In: Seminars in cancer biology. vol. 21, No. 3. Academic Press: 2011. p. 183–99.

  3. Mani SA, Guo W, Liao MJ, Eaton EN, Ayyanan A, Zhou AY, Brooks M, Reinhard F, Zhang CC, Shipitsin M, et al.The epithelial-mesenchymal transition generates cells with properties of stem cells. Cell. 2008; 133(4):704–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Ben-Porath I, Thomson MW, Carey VJ, Ge R, Bell GW, Regev A, Weinberg RA. An embryonic stem cell–like gene expression signature in poorly differentiated aggressive human tumors. Nat Genet. 2008; 40(5):499–507.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Kelloff GJ, Sigman CC. Assessing intraepithelial neoplasia and drug safety in cancer-preventive drug development. Nat Rev Cancer. 2007; 7(7):508–18.

    Article  CAS  PubMed  Google Scholar 

  6. 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. No. 7. Academic Press: 2009. p. 869–76.

  7. Li CW, Xia W, Huo L, Lim SO, Wu Y, Hsu JL, Chao CH, Yamaguchi H, Yang NK, Ding Q, et al.Epithelial–mesenchymal transition induced by TNFa requires NF-kB–mediated transcriptional upregulation of Twist1. Cancer Res. 2012; 72(5):1290–300.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Mendoza L, Alvarez-Buylla ER. Dynamics of the genetic regulatory network for Arabidopsis thaliana flower morphogenesis. J Theor Biol. 1998; 193(2):307–19.

    Article  CAS  PubMed  Google Scholar 

  9. 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. Plant Cell Online. 2004; 16(11):2923–939.

    Article  CAS  Google Scholar 

  10. Huang S, Kauffman SA. ComplexGRN complex GeneComplex GRN Regulatory Networks-from Structure to Biological Observables: Cell Fate DeterminationGene regulation, cell fate determination. In: Encyclopedia of complexity and systems science. New York: Springer: 2009. p. 1180–1213.

    Google Scholar 

  11. Alvarez-Buylla ER, Azpeitia E, Barrio R, Ben itez M, Padilla-Longoria P. From ABC genes to regulatory networks, epigenetic landscapes and flower morphogenesis: making biological sense of theoretical approaches. Semin Cell Dev Biol. 2010; 21(1):108–17.

    Article  CAS  PubMed  Google Scholar 

  12. Huang S. Reprogramming cell fates: reconciling rarity with robustness. Bioessays. 2009; 31(5):546–60.

    Article  CAS  PubMed  Google Scholar 

  13. Álvarez-Buylla ER, Chaos Á, Aldana M, Benítez M, Cortes-Poza Y, Espinosa-Soto C, Hartasánchez Da, Lotto RB, Malkin D, Escalera Santos GJ, Padilla-Longoria P. Floral morphogenesis: Stochastic explorations of a gene network epigenetic landscape. PLoS ONE. 2008; 3(11). doi:10.1371/journal.pone.0003626.

  14. Huang S. The molecular and mathematical basis of Waddington’s epigenetic landscape: a framework for post-Darwinian biology?Bioessays. 2012; 34(2):149–57.

    Article  CAS  PubMed  Google Scholar 

  15. Davila-Velderrain J, Alvarez-Buylla ER. Bridging genotype and phenotype. Frontiers in Ecology, Evolution and Complexity. 2014;144–54. CopIt ArXives.

  16. 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 Syst Biol. 2010; 4(1):134.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Azpeitia E, Davila-Velderrain J, Villarreal C, Alvarez-Buylla ER. Gene regulatory network models for floral organ determination. Flower Development: Methods and Protocols. 2014;441–69.

  18. Chakrabarti R, Hwang J, Blanco MA, Wei Y, Lukačišin M, Romano RA, Smalley K, Liu S, Yang Q, Ibrahim T, et al.Elf5 inhibits the epithelial–mesenchymal transition in mammary gland development and breast cancer metastasis by transcriptionally repressing Snail2. Nat Cell Biol. 2012; 14(11):1212–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Zhou J, Ng AY, Tymms MJ, Jermiin LS, Seth AK, Thomas RS, Kola I. 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–732.

    Article  CAS  PubMed  Google Scholar 

  20. Dave N, Guaita-Esteruelas S, Gutarra S, Frias À, Beltran M, Peiró S, de Herreros AG. Functional cooperation between Snail1 and twist in the regulation of ZEB1 expression during epithelial to mesenchymal transition. J Biol Chem. 2011; 286(14):12024–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Iaquinta PJ, Lees JA. Life and death decisions by the E2F transcription factors. Curr Opin Cell Biol. 2007; 19(6):649–57. doi:10.1016/

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Sherr CJ, Roberts JM. Living with or without cyclins and cyclin-dependent kinases. Genes Dev. 2004; 18(22):2699–711. doi:10.1101/gad.1256504.

    Article  CAS  PubMed  Google Scholar 

  23. 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–61.

    Article  CAS  PubMed  Google Scholar 

  24. Harley CB, Futcher AB, Greider CW. Telomeres shorten during ageing of human fibroblasts. Nature. 1990; 345(6274):458–60.

    Article  CAS  PubMed  Google Scholar 

  25. 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):2888.

    Article  Google Scholar 

  26. Zglinicki TV, Saretzki G, Ladhoff J, Adda F, Jackson SP. Human cell senescence as a DNA damage response. Mech Ageing Dev. 2005; 126:111–7. doi:10.1016/j.mad.2004.09.034.

    Article  Google Scholar 

  27. 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. Mol Cell Biol. 1999; 19(3):1981–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Li J, Poi MJ, Tsai MD. Regulatory mechanisms of tumor suppressor P16INK4A and their relevance to cancer. Biochemistry. 2011; 50(25):5566–582.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Hoesel B, Schmid JA. The complexity of NF- κB signaling in inflammation and cancer. Mol Cancer. 2013; 12:86. doi:10.1186/1476-4598-12-86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Khalaf H, Jass J, Olsson PE. Differential cytokine regulation by NF-kappaB and AP-1 in Jurkat T-cells. BMC Immunol. 2010; 11:26. doi:10.1186/1471-2172-11-26.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Freudlsperger C, Bian Y, Contag Wise S, Burnett J, Coupar J, Yang X, Chen Z, Van Waes C. TGF- β and NF- κB signal pathway cross-talk is mediated through TAK1 and SMAD7 in a subset of head and neck cancers. Oncogene. 2013; 32(12):1549–59. doi:10.1038/onc.2012.171.

    Article  CAS  PubMed  Google Scholar 

  32. Wu J, Duan R, Cao H, Field D, Newnham CM, Koehler DR, Zamel N, Pritchard MA, Hertzog P, Post M, Tanswell AK, Hu J. Regulation of epithelium-specific Ets-like factors ESE-1 and ESE-3 in airway epithelial cells: potential roles in airway inflammation. Cell Res. 2008; 18(6):649–63. doi:10.1038/cr.2008.57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Wang Y, Yue B, Yu X, Wang Z, Wang M. SLUG is activated by nuclear factor kappa B and confers human alveolar epithelial A549 cells resistance to tumor necrosis factor-alpha-induced apoptosis. World J Surg Oncol. 2013; 11(1):1. doi:10.1186/1477-7819-11-12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Guttridge DC, Albanese C, Reuther JY, Pestell RG, Baldwin AS. NF- κ B controls cell growth and differentiation through transcriptional regulation of Cyclin D1. Mol Cell Biol. 1999; 19(8):5785–99.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Gurova KV, Hill JE, Guo C, Prokvolit A, Burdelya LG, Samoylova E, Khodyakova AV, Ganapathi R, Ganapathi M, Tararova ND, Bosykh D, Lvovskiy D, Webb TR, Stark GR, Gudkov AV. Small molecules that reactivate p53 in renal cell carcinoma reveal a NF-kappaB-dependent mechanism of p53 suppression in tumors. Proc Natl Acad Sci U S A. 2005; 102(48):17448–53. doi:10.1073/pnas.0508888102.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Longoni N, Sarti M, Albino D, Civenni G, Malek A, Ortelli E, Pinton S, Mello-Grand M, Ostano P, D’Ambrosio G, Sessa F, Garcia-Escudero R, Thalmann GN, Chiorino G, Catapano CV, Carbone GM. ETS transcription factor ESE1/ELF3 orchestrates a positive feedback loop that constitutively activates NF- κB and drives prostate cancer progression. Cancer Res. 2013; 73(14):4533–547. doi:10.1158/0008-5472.CAN-12-4537.

    Article  CAS  PubMed  Google Scholar 

  37. De Herreros AG, Peiró S, Nassour M, Savagner P. Snail family regulation and epithelial mesenchymal transitions in breast cancer progression. J Mammary Gland Biol Neoplasia. 2010; 15(2):135–47. doi:10.1007/s10911-010-9179-8.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Tsai CC, Chen YJ, Yew TL, Chen LL, Wang JY, Chiu CH, Hung SC. Hypoxia inhibits senescence and maintains mesenchymal stem cell properties through down-regulation of E2A-p21 by HIF-TWIST. Blood. 2011; 117(2):459–69. doi:10.1182/blood-2010-05-287508.

    Article  CAS  PubMed  Google Scholar 

  39. Vega S, Morales AV, Ocaña OH, Valdés F, Fabregat I, Nieto MA. Snail blocks the cell cycle and confers resistance to cell death. Genes Dev. 2004; 18(10):1131–43. doi:10.1101/gad.294104.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Siletz A, Schnabel M, Kniazeva E, Schumacher AJ, Shin S, Jeruss JS, Shea LD. Dynamic transcription factor networks in epithelial-mesenchymal transition in breast cancer models. PLoS ONE. 2013; 8(4). doi:10.1371/journal.pone.0057180.

  41. Yang MH, Hsu DS-S, Wang HW, Wang HJ, Lan HY, Yang WH, Huang CH, Kao SY, Tzeng CH, Tai SK, Chang SY, Lee OK-S, Wu KJ. Bmi1 is essential in Twist1-induced epithelial-mesenchymal transition. Nat Cell Biol. 2010; 12(10):982–92. doi:10.1038/ncb2099.

    Article  PubMed  Google Scholar 

  42. Vichalkovski A, Gresko E, Hess D, Restuccia DF, Hemmings BA. PKB/AKT phosphorylation of the transcription factor Twist-1 at Ser42 inhibits p53 activity in response to DNA damage. Oncogene. 2010; 29(24):3554–565. doi:10.1038/onc.2010.115.

    Article  CAS  PubMed  Google Scholar 

  43. Wu K, Bonavida B. The activated NF-kappaB-Snail-RKIP circuitry in cancer regulates both the metastatic cascade and resistance to apoptosis by cytotoxic drugs. Crit Rev Immunol. 2009; 29(3):241–54.

    Article  CAS  PubMed  Google Scholar 

  44. Overhoff MG, Garbe JC, Koh J, Stampfer MR, Beach DH, Bishop CL. Cellular senescence mediated by p16INK4A-coupled miRNA pathways. Nucleic Acids Res. 2014;2003; 42(3):1606–18. doi:10.1093/nar/gkt1096.

    Article  CAS  PubMed  Google Scholar 

  45. Lowe SW, Sherr CJ. Tumor suppression by Ink4a-Arf: progress and puzzles. Curr Opin Genet Dev. 2003; 13(1):77–83.

    Article  CAS  PubMed  Google Scholar 

  46. Villacañas Ó, Pérez JJ, Rubio-Martínez J. Structural analysis of the inhibition of Cdk4 and Cdk6 by p16INK4a through molecular dynamics simulations. J Biomol Struct Dyn. 2002; 20(3):347–58.

    Article  PubMed  Google Scholar 

  47. Chien Y, Scuoppo C, Wang X, Fang X, Balgley B, Bolden JE, Premsrirut P, Luo W, Chicas A, Lee CS, Kogan SC, Lowe SW. Control of the senescence-associated secretory phenotype by NF- κB promotes senescence and enhances chemosensitivity. Genes Dev. 2011; 25(20):2125–136. doi:10.1101/gad.17276711.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Kim H, You S, Kim IJ, Foster LK, Farris J, Ambady S, Ponce de León Fa, Foster DN. Alterations in p53 and E2F-1 function common to immortalized chicken embryo fibroblasts. Oncogene. 2001; 20(21):2671–82. doi:10.1038/sj.onc.1204378.

    Article  CAS  PubMed  Google Scholar 

  49. Lee RJ, Albanese C, Fu M, D’Amico M, Lin B, Watanabe G, Haines GK, Siegel PM, Hung MC, Yarden Y, Horowitz JM, Muller WJ, Pestell RG. Cyclin D1 is required for transformation by activated Neu and is induced through an E2F-dependent signaling pathway. Mol Cell Biol. 2000; 20(2):672–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Bracken AP, Kleine-Kohlbrecher D, Dietrich N, Pasini D, Gargiulo G, Beekman C, Theilgaard-Mönch K, Minucci S, Porse BT, Marine JC, et al.The Polycomb group proteins bind throughout the INK4A-ARF locus and are disassociated in senescent cells. Genes Dev. 2007; 21(5):525–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Kato J, Matsushime H, Hiebert SW, Ewen ME, Sherr CJ. Direct binding of cyclin D to the retinoblastoma gene product (pRb) and pRb phosphorylation by the cyclin D-dependent kinase CDK4. Genes Dev. 1993; 7(3):331–42.

    Article  CAS  PubMed  Google Scholar 

  52. Duronio RJ, O’Farrell PH, Xie JE, Brook A, Dyson N. The transcription factor E2F is required for S phase during Drosophila embryogenesis. Genes Dev. 1995; 9(12):1445–55.

    Article  CAS  PubMed  Google Scholar 

  53. Suo G, Sadarangani A, Tang W, Cowan BD, Wang JYJ. Telomerase expression abrogates rapamycin-induced irreversible growth arrest of uterine fibroid smooth muscle cells. Reprod Sci. 2014; 21(9):1161–70. doi:10.1177/1933719114532839.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Herbig U, Jobling WA, Chen BPC, Chen DJ, Sedivy JM. Telomere shortening triggers senescence of human cells through a pathway involving ATM, p53, and p21(CIP1), but not p16(INK4a). Mol Cell. 2004; 14(4):501–13.

    Article  CAS  PubMed  Google Scholar 

  55. Davila-Velderrain J, Martinez-Garcia JC, Alvarez-Buylla ER. Descriptive vs. mechanistic network models in plant development in the post-genomic era. Plant Funct Genomics Methods Protoc. 2015;455–79.

  56. Romanov SR, Kozakiewicz BK, Holst CR, Stampfer MR, Haupt LM, Tlsty TD. Normal human mammary epithelial cells spontaneously escape senescence and acquire genomic changes. Nature. 2001; 409(6820):633–7. doi:10.1038/35054579.

    Article  CAS  PubMed  Google Scholar 

  57. Creixell P, Schoof EM, Erler JT, Linding R. Navigating cancer network attractors for tumor-specific therapy. Nat Biotechnol. 2012; 30(9):842–8.

    Article  CAS  PubMed  Google Scholar 

  58. Ressler S, Bartkova J, Niederegger H, Bartek J, Scharffetter-Kochanek K, Jansen-Dürr P, Wlaschek M. p16INK4A is a robust in vivo biomarker of cellular aging in human skin. Aging Cell. 2006; 5(5):379–89. doi:10.1111/j.1474-9726.2006.00231.x.

    Article  CAS  PubMed  Google Scholar 

  59. Sharma HW, Sokoloski JA, Perez JR, Maltese JY, Sartorelli AC, Stein CA, Nichols G, Khaled Z, Telang NT, Narayanan R. Differentiation of immortal cells inhibits telomerase activity. Proc Natl Acad Sci U S A. 1995; 92(26):12343–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Oeckinghaus A, Ghosh S. The NF- B Family of Transcription Factors and Its Regulation. Cold Spring Harb Perspect Biol. 2009; 1(4):000034–000034. doi:10.1101/cshperspect.a000034.

    Article  Google Scholar 

  61. Fujikawa M, Katagiri T, Tugores A, Nakamura Y, Ishikawa F. ESE-3, an Ets family transcription factor, is up-regulated in cellular senescence. Cancer Sci. 2007; 98(9):1468–75.

    Article  CAS  PubMed  Google Scholar 

  62. Vernier M, Bourdeau V, Gaumont-Leclerc MF, Moiseeva O, Bégin V, Saad F, Mes-Masson AM, Ferbeyre G. Regulation of E2Fs and senescence by PML nuclear bodies. Genes Dev. 2011; 25(1):41–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Zeisberg M, Neilson EG, et al.Biomarkers for epithelial-mesenchymal transitions. J Clin Invest. 2009; 119(6):1429–37.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. 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.

    Article  CAS  Google Scholar 

  65. Hajra KM, Chen DYS, Fearon ER. The SLUG zinc-finger protein represses E-cadherin in breast cancer. Cancer Res. 2002; 62(6):1613–8.

    CAS  PubMed  Google Scholar 

  66. Roy HK, Iversen P, Hart J, Liu Y, Koetsier JL, Kim Y, Kunte DP, Madugula M, Backman V, Wali RK. Down-regulation of SNAIL suppresses MIN mouse tumorigenesis: modulation of apoptosis, proliferation, and fractal dimension. Mol Cancer Ther. 2004; 3(9):1159–65. doi:3/9/1159.

    CAS  PubMed  Google Scholar 

  67. Wu Y, Zhou BP. TNF-alpha/NF-kappaB/Snail pathway in cancer cell migration and invasion. Br J Cancer. 2010; 102(4):639–44. doi:10.1038/sj.bjc.6605530.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Devarajan E, Song Y-h, Krishnappa S, Alt E. Epithelial—mesenchymal transition in breast cancer lines is mediated through PDGF-D released by tissue-resident stem cells. Int J Cancer. 2012; 131:1023–31. doi:10.1002/ijc.26493.

    Article  CAS  PubMed  Google Scholar 

  69. Rasmusson I. Immune modulation by mesenchymal stem cells. Exp Cell Res. 2006; 312(12):2169–179. doi:10.1016/j.yexcr.2006.03.019.

    Article  CAS  PubMed  Google Scholar 

  70. Xing F, Saidou J, Watabe K. Cancer associated fibroblasts (CAFs) in tumor microenvironment. Front Biosci. 2010; 15:166–79. doi:10.1097/MPG.0b013e3181a15ae8.Screening. NIHMS150003

    Article  CAS  Google Scholar 

  71. Hackett JA, Feldser DM, Greider CW. Telomere dysfunction increases mutation rate and genomic instability. Cell. 2001; 106(3):275–86. doi:10.1016/S0092-8674(01)00457-3.

    Article  CAS  PubMed  Google Scholar 

  72. Dolcet X, Llobet D, Pallares J, Matias-Guiu X. NF-kB in development and progression of human cancer. Virchows Arch. 2005; 446(5):475–82. doi:10.1007/s00428-005-1264-9.

    Article  CAS  PubMed  Google Scholar 

  73. Lee JM, Dedhar S, Kalluri R, Thompson EW. The epithelial–mesenchymal transition: new insights in signaling, development, and disease. J Cell Biol. 2006; 172(7):973–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Battula VL, Evans KW, Hollier BG, Shi Y, Marini FC, Ayyanan A, Wang R-y, Brisken C, Guerra R, Andreeff M, et al.Epithelial-Mesenchymal Transition-Derived Cells Exhibit Multilineage Differentiation Potential Similar to Mesenchymal Stem Cells. Stem Cells. 2010; 28(8):1435–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Moustakas A, Heldin CH. Mechanisms of TGF β -Induced Epithelial — Mesenchymal Transition. J Clin Med. 2016; 5:1–34.

    Article  Google Scholar 

  76. Steinway SN, Zañudo JGT, Ding W, Rountree CB, Feith D, Loughran TP, Albert R. Network modeling of TGFb signaling in hepatocellular carcinoma epithelial-to-mesenchymal transition reveals joint sonic hedgehog and Wnt pathway activation. Cancer Res. 2014; 21:5963–977.

    Article  Google Scholar 

  77. Kim WY, Sharpless NE. The regulation of INK4 ARF in cancer and aging. Cell. 2006; 127(2):265–75.

    Article  CAS  PubMed  Google Scholar 

  78. Weinberg RA. Twisted epithelial–mesenchymal transition blocks senescence. Nat Cell Biol. 2008; 10(9):1021–3.

    Article  CAS  PubMed  Google Scholar 

  79. Smit MA, Peeper DS. Epithelial-mesenchymal transition and senescence: two cancer-related processes are crossing paths. Aging (Albany NY). 2010; 2(10):735.

    Article  CAS  Google Scholar 

  80. Xu J, Lamouille S, Derynck R. TGF-beta-induced epithelial to mesenchymal transition. Cell Res. 2009; 19(2):156–72. doi:10.1038/cr.2009.5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Sullivan NJ, Sasser AK, Axel A, Vesuna F, Raman V, Ramirez N, Oberyszyn TM, Hall BM. Interleukin-6 induces an epithelial–mesenchymal transition phenotype in human breast cancer cells. Oncogene. 2009; 28(33):2940–947.

    Article  CAS  PubMed  Google Scholar 

  82. Coppé J-P, Patil CK, Rodier F, Sun Y, Muñoz DP, Goldstein J, Nelson PS, Desprez PY, Campisi J. Senescence-associated secretory phenotypes reveal cell-nonautonomous functions of oncogenic RAS and the p53 tumor suppressor. PLoS Biol. 2008; 6(12):301.

    Article  Google Scholar 

  83. Beausejour CM. A. Krtolica,. "Reversal of human cellular senescence: roles of the p53 and p16 pathways. EMBO J. 2003; 22(16):4212–222.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Davila-Velderrain J, Martinez-García JC, Alvarez-Buylla ER. Modeling the epigenetic attractors landscape: towards a post-genomic mechanistic understanding of development. Front Genet. 2015; 6:160.

    Article  PubMed  PubMed Central  Google Scholar 

  85. 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. Biosystems. 2016; 142:15–24.

    Article  PubMed  Google Scholar 

  86. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  87. Zhu X, Yuan R, Hood L, Ao P. Endogenous molecular-cellular hierarchical modeling of prostate carcinogenesis uncovers robust structure. Prog Biophys Mol Biol. 2015; 117(1):30–42.

    Article  CAS  PubMed  Google Scholar 

  88. Venkov CD, Link AJ, Jennings JL, Plieth D, Inoue T, Nagai K, Xu C, Dimitrova YN, Rauscher FJ, Neilson EG. A proximal activator of transcription in epithelial-mesenchymal transition. J Clin Invest. 2007; 117(2):482–91. doi:10.1172/JCI29544.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  89. Li C, Wang J. Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths. PLoS Comput Biol. 2013; 9(8):1003165.

    Article  Google Scholar 

  90. Zhou JX, Brusch L, Huang S. Predicting pancreas cell fate decisions and reprogramming with a hierarchical multi-attractor model. PLoS ONE. 2011; 6(3):14752.

    Article  Google Scholar 

  91. Zhou JX, Qiu X, d’Hérouël AF, Huang S. Discrete gene network models for understanding multicellularity and cell reprogramming: from network structure to attractor landscape. Computational Systems Biology. 2014;241–76.

  92. Steinway SN, Zañudo JGT, Michel PJ, Feith D, Loughran TP, Albert R. Combinatorial interventions inhibit TGF β-driven epithelial-to-mesenchymal transition and support hybrid cellular phenotypes. npj Syst Biol Appl. 2015; 1:15014.

    Article  Google Scholar 

  93. Davila-Velderrain J, Villarreal C, Alvarez-Buylla ER. Reshaping the epigenetic landscape during early flower development: induction of attractor transitions by relative differences in gene decay rates. BMC Syst Biol. 2015; 9(1):20–28.

    Article  PubMed  PubMed Central  Google Scholar 

  94. Glaab E, Baudot A, Krasnogor N, Schneider R, Valencia A. EnrichNet: network-based gene set enrichment analysis. Bioinformatics. 2012; 28(18):451–7.

    Article  Google Scholar 

  95. Glaab E, Baudot A, Krasnogor N, Valencia A. TopoGSA: network topological gene set analysis. Bioinformatics. 2010; 26(9):1271–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  96. 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–30.

    Article  CAS  PubMed  Google Scholar 

  97. Saadatpour A, Albert R, Reluga TC. A reduction method for Boolean network models proven to conserve attractors. SIAM J Appl Dyn Syst. 2013; 12(4):1997–2011.

    Article  Google Scholar 

  98. Müssel C, Hopfensitz M, Kestler HA. BoolNet—an R package for generation, reconstruction and analysis of Boolean networks. Bioinformatics. 2010; 26(10):1378–80.

    Article  PubMed  Google Scholar 

  99. Allen LJ. An introduction to stochastic processes with applications to biology. CRC Press. 2010.

Download references


The authors acknowledge logistical and administrative help of Diana Romo.


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. JDV receives a postdoctoral scholarship from the Centro de Ciencias de la Complejidad (C3)/Laboratorio Nacional del C3, UNAM. EDH receives a UNAM-DGAPA postdoctoral scholarship.

Availability of data and materials

The authors declare that all datasets are presented in the main manuscript and Additional files 4, 2, 3, 5 and 1.

Authors’ contributions

ERAB and JCMG 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 models and data analyses. FML and CEO formalized experimental data into regulatory logic. ERAB, JCMG, JDV, EDH and FML participated in the interpretation of the results and analyses. JDV wrote an important part of the paper with significant contributions from ERAB, EDH, FM and JCMG. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Author information

Authors and Affiliations


Corresponding authors

Correspondence to Juan Carlos Martínez-García or Elena R. Alvarez-Buylla.

Additional files

Additional file 1

Table S4. Empirical evidence supporting the construction of the large GRN. (PDF 451 kb)

Additional file 2

Figure S1. Network -based topological gene set enrichment analysis shows that the topological signature of the molecules in the proposed network is structurally close to the cancer pathways included in the KEGG database. (EPS 39 kb)

Additional file 3

Supplementary Tables. (XLXS 543 kb)

Additional file 4

Supplementary materials. (PDF 69 kb)

Additional file 5

Logic rules and truth tables. (PDF 177 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Méndez-López, L., Davila-Velderrain, J., Domínguez-Hüttinger, E. et al. Gene regulatory network underlying the immortalization of epithelial cells. BMC Syst Biol 11, 24 (2017).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: