Multi-stability in cellular differentiation enabled by a network of three mutually repressing master regulators

Network topology and interaction logic determine states it supports

In this review paper we summarize a recent progress on the problem of describing range of dynamics supported by a network. We show that there is natural connection between network models consisting of collections of multivalued monotone boolean functions and ordinary differential equations models. We show how to construct such collections and use them to answer questions about prevalence of cellular phenotypes that correspond to equilibria of network models.

Assessing biological network dynamics: comparing numerical simulations with analytical decomposition of parameter space

Mathematical modeling of the emergent dynamics of gene regulatory networks (GRN) faces a double challenge of (a) dependence of model dynamics on parameters, and (b) lack of reliable experimentally determined parameters. In this paper we compare two complementary approaches for describing GRN dynamics across unknown parameters: (1) parameter sampling and resulting ensemble statistics used by RACIPE (RAndom CIrcuit PErturbation), and (2) use of rigorous analysis of combinatorial approximation of the ODE models by DSGRN (Dynamic Signatures Generated by Regulatory Networks). We find a very good agreement between RACIPE simulation and DSGRN predictions for four different 2- and 3-node networks typically observed in cellular decision making. This observation is remarkable since the DSGRN approach assumes that the Hill coefficients of the models are very high while RACIPE assumes the values in the range 1-6. Thus DSGRN parameter domains, explicitly defined by inequalities between systems parameters, are highly predictive of ODE model dynamics within a biologically reasonable range of parameters.

Theoretical investigation of functional responses of bio-molecular assembly networks

Cooperative protein-protein and protein-DNA interactions form programmable complex assemblies, often performing non-linear gene regulatory operations involved in signal transductions and cell fate determination. The apparent structure of those complex assemblies is very similar, but their functional response strongly depends on the topology of the protein-DNA interaction networks. Here, we demonstrate how the coordinated self-assembly creates gene regulatory network motifs that corroborate the existence of a precise functional response at the molecular level using thermodynamic and dynamic analyses. Our theoretical and Monte Carlo simulations show that a complex network of interactions can form a decision-making loop, such as feedback and feed-forward circuits, only by a few molecular mechanisms. We characterize each possible network of interactions by systematic variations of free energy parameters associated with the binding among biomolecules and DNA looping. We also find that the higher-order networks exhibit alternative steady states from the stochastic dynamics of each network. We capture this signature by calculating stochastic potentials and attributing their multi-stability features. We validate our findings against the Gal promoter system in yeast cells. Overall, we show that the network topology is vital in phenotype diversity in regulatory circuits.

Universal structures for adaptation in biochemical reaction networks

At the molecular level, the evolution of life is driven by the generation and diversification of adaptation mechanisms. A universal description of adaptation-capable chemical reaction network (CRN) structures has remained elusive until now, since currently-known criteria for adaptation apply only to a tiny subset of possible CRNs. Here we identify the definitive structural requirements that characterize all adaptation-capable collections of interacting molecules, however large or complex. We show that these network structures implement a form of integral control in which multiple independent integrals can collaborate to confer the capacity for adaptation on specific molecules. Using an algebraic algorithm informed by these findings, we demonstrate the existence of embedded integrals in a variety of biologically important CRNs that have eluded previous methods, and for which adaptation has been observed experimentally. This definitive picture of biological adaptation at the level of intermolecular interactions represents a blueprint for adaptation-capable signaling networks across all domains of life, and for the design of synthetic biosystems.

Spatiotemporal Patterning Enabled by Gene Regulatory Networks

Spatiotemporal pattern formation plays a key role in various biological phenomena including embryogenesis and neural network formation. Though the reaction-diffusion systems enabling pattern formation have been studied phenomenologically, the biomolecular mechanisms behind these processes have not been modeled in detail. Here, we study the emergence of spatiotemporal patterns due to simple, synthetic and commonly observed two- and three-node gene regulatory network motifs coupled with their molecular diffusion in one- and two-dimensional space. We investigate the patterns formed due to the coupling of inherent multistable and oscillatory behavior of the toggle switch, toggle switch with double self-activation, toggle triad, and repressilator with the effect of spatial diffusion of these molecules. We probe multiple parameter regimes corresponding to different regions of stability (monostable, multistable, oscillatory) and assess the impact of varying diffusion coefficients. This analysis offers valuable insights into the design principles of pattern formation facilitated by these network motifs, and it suggests the mechanistic underpinnings of biological pattern formation.

Modelling of indirect cell-cell interaction networks mediated by IFNγ/IL-4 cytokine involved in atopic dermatitis

Atopic dermatitis (AD) is an immune-driven inflammatory skin disease that is known to have a significantly high life-time prevalence in the human population. T-helper (Th) immune cells play a key role in the pathogenesis of AD which is marked by defects in the skin barrier function along with a significant increase in the population of either Th1 or Th2 sub-types of Th cells. The progression of AD from the acute to chronic phase is still poorly understood, and here we explore the mechanism of this transition through the study of a mathematical model for indirect cell-cell interactions among Th and skin cells via the secreted cytokines IFNγ and IL-4, both known to have therapeutic potential. An increase in the level of cytokine IFN γ can catalyse the transition of AD from an acute to a chronic stage, while an increase in the level of cytokine IL-4 has the reverse effect. In our model, the transition of AD from the acute to chronic stage and vice versa can be abrupt (switch-like) with hysteresis: this bistable behaviour can potentially be used to keep AD in the acute phase since therapy based on suppression of IFNγ can retard the transition to the chronic phase.

What Makes a Functional Gene Regulatory Network? A Circuit Motif Analysis

One of the key questions in systems biology is to understand the roles of gene regulatory circuits in determining cellular states and their functions. In previous studies, some researchers have inferred large gene networks from genome wide genomics/transcriptomics data using the top-down approach, while others have modeled core gene circuits of small sizes using the bottom-up approach. Despite many existing systems biology studies, there is still no general rule on what sizes of gene networks and what types of circuit motifs a system would need to achieve robust biological functions. Here, we adopt a gene circuit motif analysis to discover four-node circuits responsible for multiplicity (rich in dynamical behavior), flexibility (versatile to alter gene expression), or both. We identify the most reoccurring two-node circuit motifs and the co-occurring motif pairs. Furthermore, we investigate the contributing factors of multiplicity and flexibility for large gene networks of different types and sizes. We find that gene networks of intermediate sizes tend to have combined high levels of multiplicity and flexibility. Our study will contribute to a better understanding of the dynamical mechanisms of gene regulatory circuits and provide insights into rational designs of robust gene circuits in synthetic and systems biology.

Functional Resilience of Mutually Repressing Motifs Embedded in Larger Networks

Elucidating the design principles of regulatory networks driving cellular decision-making has important implications for understanding cell differentiation and guiding the design of synthetic circuits. Mutually repressing feedback loops between 'master regulators' of cell fates can exhibit multistable dynamics enabling "single-positive" phenotypes: (high A, low B) and (low A, high B) for a toggle switch, and (high A, low B, low C), (low A, high B, low C) and (low A, low B, high C) for a toggle triad. However, the dynamics of these two motifs have been interrogated in isolation in silico, but in vitro and in vivo, they often operate while embedded in larger regulatory networks. Here, we embed these motifs in complex larger networks of varying sizes and connectivity to identify hallmarks under which these motifs maintain their canonical dynamical behavior. We show that an increased number of incoming edges onto a motif leads to a decay in their canonical stand-alone behaviors. We also show that this decay can be exacerbated by adding self-inhibition but not self-activation loops on the 'master regulators'. These observations offer insights into the design principles of biological networks containing these motifs and can help devise optimal strategies for the integration of these motifs into larger synthetic networks.

Rewiring the altered tryptophan metabolism as a novel therapeutic strategy in inflammatory bowel diseases

OBJECTIVE

The extent to which tryptophan (Trp) metabolism alterations explain or influence the outcome of inflammatory bowel diseases (IBDs) is still unclear. However, several Trp metabolism end-products are essential to intestinal homeostasis. Here, we investigated the role of metabolites from the kynurenine pathway.

DESIGN

Targeted quantitative metabolomics was performed in two large human IBD cohorts (1069 patients with IBD). Dextran sodium sulphate-induced colitis experiments in mice were used to evaluate effects of identified metabolites. In vitro, ex vivo and in vivo experiments were used to decipher mechanisms involved. Effects on energy metabolism were evaluated by different methods including Single Cell mEtabolism by profiling Translation inHibition.

RESULTS

In mice and humans, intestinal inflammation severity negatively correlates with the amount of xanthurenic (XANA) and kynurenic (KYNA) acids. Supplementation with XANA or KYNA decreases colitis severity through effects on intestinal epithelial cells and T cells, involving Aryl hydrocarbon Receptor (AhR) activation and the rewiring of cellular energy metabolism. Furthermore, direct modulation of the endogenous tryptophan metabolism, using the recombinant enzyme aminoadipate aminotransferase (AADAT), responsible for the generation of XANA and KYNA, was protective in rodent colitis models.

CONCLUSION

Our study identified a new mechanism linking Trp metabolism to intestinal inflammation and IBD. Bringing back XANA and KYNA has protective effects involving AhR and the rewiring of the energy metabolism in intestinal epithelial cells and CD4⁺ T cells. This study paves the way for new therapeutic strategies aiming at pharmacologically correcting its alterations in IBD by manipulating the endogenous metabolic pathway with AADAT.

Dynamics of hepatocyte-cholangiocyte cell-fate decisions during liver development and regeneration

The immense regenerative potential of the liver is attributed to the ability of its two key cell types – hepatocytes and cholangiocytes – to trans-differentiate to one another either directly or through intermediate progenitor states. However, the dynamic features of decision-making between these cell-fates during liver development and regeneration remains elusive. Here, we identify a core gene regulatory network comprising c/EBPα, TGFBR2, and SOX9 which is multistable in nature, enabling three distinct cell states – hepatocytes, cholangiocytes, and liver progenitor cells (hepatoblasts/oval cells) – and stochastic switching among them. Predicted expression signature for these three states are validated through multiple bulk and single-cell transcriptomic datasets collected across developmental stages and injury-induced liver repair. This network can also explain the experimentally observed spatial organization of phenotypes in liver parenchyma and predict strategies for efficient cellular reprogramming. Our analysis elucidates how the emergent dynamics of underlying regulatory networks drive diverse cell-fate decisions in liver development and regeneration.

•

Identified minimal regulatory network to model liver development and regeneration

•

Changes in phenotypic landscapes by in-silico perturbations of regulatory networks

•

Ability to explain physiological spatial patterning of liver cell types

•

Decoded strategies for efficient reprogramming among liver cell phenotypes

Lack of Correlation between Landscape Geometry and Transition Rates

Biological cells can exist in a variety of distinct phenotypes, determined by the steady-state solutions of genetic networks governing their cell fate. A popular way of representing these states relies on the creation of landscape related to the relative occupation of these states. It is often assumed that this landscape offers direct information regarding the state-to-state transition rates, suggesting that these are related to barrier heights separating landscape minima. Here, we study a toggle triad network exhibiting multistability and directly demonstrate the lack of any direct correlation between properties of the landscape and corresponding transition rates.

Emergent dynamics of a three-node regulatory network explain phenotypic switching and heterogeneity: a case study of Th1/Th2/Th17 cell differentiation

Naïve helper (CD4+) T-cells can differentiate into distinct functional subsets including Th1, Th2, and Th17 phenotypes. Each of these phenotypes has a 'master regulator' - T-bet (Th1), GATA3 (Th2) and RORγT (Th17) - that inhibits the other two master regulators. Such mutual repression among them at a transcriptional level can enable multistability, giving rise to six experimentally observed phenotypes - Th1, Th2, Th17, hybrid Th/Th2, hybrid Th2/Th17 and hybrid Th1/Th17. However, the dynamics of switching among these phenotypes, particularly in the case of epigenetic influence, remains unclear. Here, through mathematical modeling, we investigated the coupled transcription-epigenetic dynamics in a three-node mutually repressing network to elucidate how epigenetic changes mediated by any 'master regulator' can influence the transition rates among different cellular phenotypes. We show that the degree of plasticity exhibited by one phenotype depends on relative strength and duration of mutual epigenetic repression mediated among the master regulators in a three-node network. Further, our model predictions can offer putative mechanisms underlying relatively higher plasticity of Th17 phenotype as observed in vitro and in vivo. Together, our modeling framework characterizes phenotypic plasticity and heterogeneity as an outcome of emergent dynamics of a three-node regulatory network, such as the one mediated by T-bet/GATA3/RORγT.

A robust method for designing multistable systems by embedding bistable subsystems

Although multistability is an important dynamic property of a wide range of complex systems, it is still a challenge to develop mathematical models for realising high order multistability using realistic regulatory mechanisms. To address this issue, we propose a robust method to develop multistable mathematical models by embedding bistable models together. Using the GATA1-GATA2-PU.1 module in hematopoiesis as the test system, we first develop a tristable model based on two bistable models without any high cooperative coefficients, and then modify the tristable model based on experimentally determined mechanisms. The modified model successfully realises four stable steady states and accurately reflects a recent experimental observation showing four transcriptional states. In addition, we develop a stochastic model, and stochastic simulations successfully realise the experimental observations in single cells. These results suggest that the proposed method is a general approach to develop mathematical models for realising multistability and heterogeneity in complex systems.

Modular, robust and extendible multicellular circuit design in yeast

Division of labor between cells is ubiquitous in biology but the use of multi-cellular consortia for engineering applications is only beginning to be explored. A significant advantage of multi-cellular circuits is their potential to be modular with respect to composition but this claim has not yet been extensively tested using experiments and quantitative modeling. Here, we construct a library of 24 yeast strains capable of sending, receiving or responding to three molecular signals, characterize them experimentally and build quantitative models of their input-output relationships. We then compose these strains into two- and three-strain cascades as well as a four-strain bistable switch and show that experimentally measured consortia dynamics can be predicted from the models of the constituent parts. To further explore the achievable range of behaviors, we perform a fully automated computational search over all two-, three- and four-strain consortia to identify combinations that realize target behaviors including logic gates, band-pass filters and time pulses. Strain combinations that are predicted to map onto a target behavior are further computationally optimized and then experimentally tested. Experiments closely track computational predictions. The high reliability of these model descriptions further strengthens the feasibility and highlights the potential for distributed computing in synthetic biology.

Effects of immunosuppressants on T-cell dynamics: Understanding from a generic coarse-grained immune network model

Long-term immunosuppressive therapy is a drug regimen often used to lower aggressive immune responses in various chronic inflammatory diseases. However, such long-term therapy leading to immune suppression may trigger other adverse reactions in the immune system. The rising concern regarding the optimal dose and duration of such treatment has motivated us to understand non-classical immunomodulatory responses induced by various immunosuppressive steroid and secosteroid drugs such as glucocorticoid and vitamin D supplements. The immunomodulatory actions of such immunosuppressants (that govern the adaptive immune response) are often mediated through their characteristic control over CD4+ T-cells involving pro- and antiinflammatory T-cells. Several early studies attempted to decode temporal and dose-dependent behaviors of such pro- and anti-inflammatory T-cells using the chemical dynamics approach. We first summarize these early works. Then, we develop a minimal coarse-grained kinetic network model to capture the commonality in their immunomodulatory functions. This generic model successfully reproduces the characteristic dynamical features, including the clinical latency period in long-term T-cell dynamics. The temporal behavior of T-cells is found to be sensitive to specific rate parameters and doses of immunosuppressants. The steady-state analysis reflects the transition from an early classified weakly regulated (autoimmune-prone) immune state to a strongly regulated state (immunocompromised state), separated by an intervening state of moderate/balanced regulation. An optimal dose and duration are essential in rescuing balanced immune regulation. This review elucidates how developing a simple generic coarse-grained immune network model may provide immense information that helps diagnose inefficacy in adaptive immune function before and after administering immunosuppressants such as glucocorticoid or vitamin D.

KLF4 Induces Mesenchymal-Epithelial Transition (MET) by Suppressing Multiple EMT-Inducing Transcription Factors

Epithelial-Mesenchymal Plasticity (EMP) refers to reversible dynamic processes where cells can transition from epithelial to mesenchymal (EMT) or from mesenchymal to epithelial (MET) phenotypes. Both these processes are modulated by multiple transcription factors acting in concert. While EMT-inducing transcription factors (TFs)-TWIST1/2, ZEB1/2, SNAIL1/2/3, GSC, and FOXC2-are well-characterized, the MET-inducing TFs are relatively poorly understood (OVOL1/2 and GRHL1/2). Here, using mechanism-based mathematical modeling, we show that transcription factor KLF4 can delay the onset of EMT by suppressing multiple EMT-TFs. Our simulations suggest that KLF4 overexpression can promote a phenotypic shift toward a more epithelial state, an observation suggested by the negative correlation of KLF4 with EMT-TFs and with transcriptomic-based EMT scoring metrics in cancer cell lines. We also show that the influence of KLF4 in modulating the EMT dynamics can be strengthened by its ability to inhibit cell-state transitions at the epigenetic level. Thus, KLF4 can inhibit EMT through multiple parallel paths and can act as a putative MET-TF. KLF4 associates with the patient survival metrics across multiple cancers in a context-specific manner, highlighting the complex association of EMP with patient survival.

Coexistence of Hopf-born rotation and heteroclinic cycling in a time-delayed three-gene auto-regulated and mutually-repressed core genetic regulation network

In this work, we study the behavior of a time-delayed mutually repressive auto-activating three-gene system. Delays are introduced to account for the location difference between DNA transcription that leads to production of messenger RNA and its translation that result in protein synthesis. We study the dynamics of the system using numerical simulations, computational bifurcation analysis and mathematical analysis. We find Hopf bifurcations leading to stable and unstable rotation in the system, and we study the rotational behavior as a function of cyclic mutual repression parameter asymmetry between each gene pair in the network. We focus on how rotation co-exists with a stable heteroclinic flow linking the three saddles in the system. We find that this coexistence allows for a transition between two markedly different types of rotation leading to strikingly different phenotypes. One type of rotation belongs to Hopf-induced rotation while the other type, belongs to heteroclinic cycling between three saddle nodes in the system. We discuss the evolutionary and biological implications of our findings.

Immune phase transition under steroid treatment

The steroid hormone glucocorticoid (GC) is a well-known immunosuppressant that controls T-cell-mediated adaptive immune response. In this work, we have developed a minimal kinetic network model of T-cell regulation connecting relevant experimental and clinical studies to quantitatively understand the long-term effects of GC on pro-inflammatory T-cell (T_{pro}) and anti-inflammatory T-cell (T_{anti}) dynamics. Due to the antagonistic relation between these two types of T cells, their long-term steady-state population ratio helps us to characterize three classified immune regulations: (i) weak ([T_{pro}]>[T_{anti}]), (ii) strong ([T_{pro}]<[T_{anti}]), and (iii) moderate ([T_{pro}]∼[T_{anti}]), holding the characteristic bistability. In addition to the differences in their long-term steady-state outcome, each immune regulation shows distinct dynamical phases. In the presteady state, a characteristic intermediate stationary phase is observed to develop only in the moderate regulation regime. In the medicinal field, the resting time in this stationary phase is distinguished as a clinical latent period. GC dose-dependent steady-state analysis shows an optimal level of GC to drive a phase transition from the weak or autoimmune prone to the moderate regulation regime. Subsequently, the presteady state clinical latent period tends to diverge near that optimal GC level where [T_{pro}]:[T_{anti}] is highly balanced. The GC-optimized elongated stationary phase explains the rationale behind the requirement of long-term immune diagnostics, especially when long-term GC-based chemotherapeutics and other immunosuppressive drugs are administrated. Moreover, our study reveals GC sensitivity of clinical latent period, which might serve as an early warning signal in diagnosing different immune phases and determining immune phasewise steroid treatment.

Operating principles of circular toggle polygons

Decoding the dynamics of cellular decision-making and cell differentiation is a central question in cell and developmental biology. A common network motif involved in many cell-fate decisions is a mutually inhibitory feedback loop between two self-activating 'master regulators' A and B, also called as toggle switch. Typically, it can allow for three stable states-(high A, low B), (low A, high B) and (medium A, medium B). A toggle triad-three mutually repressing regulators A, B and C, i.e. three toggle switches arranged circularly (between A and B, between B and C, and between A and C)-can allow for six stable states: three 'single positive' and three 'double positive' ones. However, the operating principles of larger toggle polygons, i.e. toggle switches arranged circularly to form a polygon, remain unclear. Here, we simulate using both discrete and continuous methods the dynamics of different sized toggle polygons. We observed a pattern in their steady state frequency depending on whether the polygon was an even or odd numbered one. The even-numbered toggle polygons result in two dominant states with consecutive components of the network expressing alternating high and low levels. The odd-numbered toggle polygons, on the other hand, enable more number of states, usually twice the number of components with the states that follow 'circular permutation' patterns in their composition. Incorporating self-activations preserved these trends while increasing the frequency of multistability in the corresponding network. Our results offer insights into design principles of circular arrangement of regulatory units involved in cell-fate decision making, and can offer design strategies for synthesizing genetic circuits.

Topological signatures in regulatory network enable phenotypic heterogeneity in small cell lung cancer

Phenotypic (non-genetic) heterogeneity has significant implications for the development and evolution of organs, organisms, and populations. Recent observations in multiple cancers have unraveled the role of phenotypic heterogeneity in driving metastasis and therapy recalcitrance. However, the origins of such phenotypic heterogeneity are poorly understood in most cancers. Here, we investigate a regulatory network underlying phenotypic heterogeneity in small cell lung cancer, a devastating disease with no molecular targeted therapy. Discrete and continuous dynamical simulations of this network reveal its multistable behavior that can explain co-existence of four experimentally observed phenotypes. Analysis of the network topology uncovers that multistability emerges from two teams of players that mutually inhibit each other, but members of a team activate one another, forming a 'toggle switch' between the two teams. Deciphering these topological signatures in cancer-related regulatory networks can unravel their 'latent' design principles and offer a rational approach to characterize phenotypic heterogeneity in a tumor.

Multi-Stability and Consequent Phenotypic Plasticity in AMPK-Akt Double Negative Feedback Loop in Cancer Cells

Adaptation and survival of cancer cells to various stress and growth factor conditions is crucial for successful metastasis. A double-negative feedback loop between two serine/threonine kinases AMPK (AMP-activated protein kinase) and Akt can regulate the adaptation of breast cancer cells to matrix-deprivation stress. This feedback loop can significantly generate two phenotypes or cell states: matrix detachment-triggered pAMPK^high/ pAkt^low state, and matrix (re)attachment-triggered pAkt^high/ pAMPK^low state. However, whether these two cell states can exhibit phenotypic plasticity and heterogeneity in a given cell population, i.e., whether they can co-exist and undergo spontaneous switching to generate the other subpopulation, remains unclear. Here, we develop a mechanism-based mathematical model that captures the set of experimentally reported interactions among AMPK and Akt. Our simulations suggest that the AMPK-Akt feedback loop can give rise to two co-existing phenotypes (pAkt^high/ pAMPK^low and pAMPK^high/pAkt^low) in specific parameter regimes. Next, to test the model predictions, we segregated these two subpopulations in MDA-MB-231 cells and observed that each of them was capable of switching to another in adherent conditions. Finally, the predicted trends are supported by clinical data analysis of The Cancer Genome Atlas (TCGA) breast cancer and pan-cancer cohorts that revealed negatively correlated pAMPK and pAkt protein levels. Overall, our integrated computational-experimental approach unravels that AMPK-Akt feedback loop can generate multi-stability and drive phenotypic switching and heterogeneity in a cancer cell population.

Hybrid E/M Phenotype(s) and Stemness: A Mechanistic Connection Embedded in Network Topology

Metastasis remains an unsolved clinical challenge. Two crucial features of metastasizing cancer cells are (a) their ability to dynamically move along the epithelial-hybrid-mesenchymal spectrum and (b) their tumor initiation potential or stemness. With increasing functional characterization of hybrid epithelial/mesenchymal (E/M) phenotypes along the spectrum, recent in vitro and in vivo studies have suggested an increasing association of hybrid E/M phenotypes with stemness. However, the mechanistic underpinnings enabling this association remain unclear. Here, we develop a mechanism-based mathematical modeling framework that interrogates the emergent nonlinear dynamics of the coupled network modules regulating E/M plasticity (miR-200/ZEB) and stemness (LIN28/let-7). Simulating the dynamics of this coupled network across a large ensemble of parameter sets, we observe that hybrid E/M phenotype(s) are more likely to acquire stemness relative to "pure" epithelial or mesenchymal states. We also integrate multiple "phenotypic stability factors" (PSFs) that have been shown to stabilize hybrid E/M phenotypes both in silico and in vitro-such as OVOL1/2, GRHL2, and NRF2-with this network, and demonstrate that the enrichment of hybrid E/M phenotype(s) with stemness is largely conserved in the presence of these PSFs. Thus, our results offer mechanistic insights into recent experimental observations of hybrid E/M phenotype(s) that are essential for tumor initiation and highlight how this feature is embedded in the underlying topology of interconnected EMT (Epithelial-Mesenchymal Transition) and stemness networks.