{'fraction_non_alphanumeric': 0.03409434786994385, 'fraction_numerical': 0.011835513264516997, 'mean_word_length': 5.321360255047821, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deciphering the dynamics of Epithelial-Mesenchymal Transition and Cancer Stem Cells in tumor progression

Purpose of review: The epithelial-Mesenchymal Transition (EMT) and the generation of Cancer Stem Cells (CSC) are two fundamental aspects contributing to tumor growth, acquisition of resistance to therapy, formation of metastases, and tumor relapse. Recent experimental data identifying the circuits regulating EMT and CSCs have driven the development of computational models capturing the dynamics of these circuits and consequently various aspects of tumor progression. Recent findings: We review the contribution made by these models in a) recapitulating experimentally observed behavior, b) making experimentally testable predictions, and c) driving emerging notions in the field, including the emphasis on the aggressive potential of hybrid epithelial/mesenchymal (E/M) phenotype(s). We discuss dynamical and statistical models at intracellular and population level relating to dynamics of EMT and CSCs, and those focusing on interconnections between these two processes. Summary: These models highlight the insights gained via mathematical modeling approaches, and emphasizes that the connections between hybrid E/M phenotype(s) and stemness can be explained by analyzing underlying regulatory circuits. Such experimentally curated models have the potential of serving as platforms for better therapeutic design strategies.


Introduction
Metastasis and tumor relapse are insuperable clinical challenges that claim most cancer-related deaths [1]. The metastatic cascade has extremely high rates of attrition, because of the multi-step and challenging sequence of events leading to a secondary tumor-cancer cells detaching from their home organ, circulating in the bloodstream, lodging and colonizing in a foreign environment, all while escaping recognition and attack by the immune system and other clinical interventions.
As a first step in the metastatic cascade, cancer cells in a solid tumor tissue often undergo an Epithelial-to-Mesenchymal Transition (EMT) characterized by a loss of cell-cell adhesion and acquisition of migratory and invasive traits [2]. Afterwards, these disseminated cells travel through the bloodstream and colonize a distant organ, giving rise to macrometastases [2,3]. EMT in cancer is not necessarily a cell-autonomous and binary process. Cells can attain intermediate, or hybrid, epithelial/mesenchymal (E/M) states and/or involve their neighbors to form more aggressive clusters of circulating tumor cells (CTCs), the main drivers of metastases [4][5][6]. EMT is regulated at multiple layers including transcriptional, translational and post-translational levels. Although tissue-specific factors and micro-environment can play an important role, some common traits of EMT include transcriptional repression of E-Cadherin that mediates cell-cell junctions and adhesion, and activation of one or more EMT-inducing transcription factors (EMT-TFs) such as SNAI1, SNAI2, ZEB1, ZEB1, TWIST1 that induce cell scattering, motility, and invasion [2].
For colonizing a secondary tumor site, cells need the ability to give rise to different cell types that constitute a tumor -a trait typical of cancer stem cells (CSC). Cells with such stem-like properties are also typically resistant to various clinical treatments, and are often implicated in tumor relapse. The conventional so-called 'CSC hypothesis' envisions a small fraction of CSCs that can self-renew (symmetric division) or generate differentiated cells (asymmetric division) [7,8]. This process implies a hierarchical lineage of tumor cells that irreversibly lose their stem properties, which resembles the stem cell hierarchy in normal tissues [9]. However, recent studies have emphasized that stemness can be a dynamic cell state such that be acquired or lost by cells [10][11][12][13]; in other words, some tumor cells can dedifferentiate and re-gain stemness via epigenetic and/or environmental factors such as abnormal cancer metabolism and EMT.
The connection between EMT and CSC was first postulated by Brabletz et al. [14] in 2005 as 'migrating cancer stem cells' by suggesting that the concepts of EMT and CSCs, considered independent of one another, were not sufficient to explain various traits of cancer progression. Afterwards, experimental evidence accumulated suggesting that stemness can be gained during EMT [2,5,[15][16][17][18][19]. Recent experiments have shown that cells in intermediate E/M states bear a higher metastatic potential as compared to the cells that have undergone a complete EMT. This behavior, together with enhanced drug-resistance traits of the hybrid E/M phenotype(s) [20,21], underlies the clinical implications of hybrid E/M phenotypes [5,22,23].
Recent studies have made significant progress in identifying the molecular networks regulating EMT, CSCs, and their interconnections [24]. These networks are formidably complex, and capable of giving rise to emergent non-linear behavior. Identification of these networks has driven a surge in deciphering their underlying principles from a dynamical systems perspective, and has thus led to the development of many computational models to capture the dynamics of these transitions, models that may deal with intracellular and intercellular circuits, or may stick to a population level description. Here, we review both of these classes of models. We review a set of models that attempt to clarify the possible set of states accessible to EMT and their possible relevance to tumor progression and metastasis. We review a set of models that consider the population structure of the tumor and its implications for drug resistance. Finally, we discuss models that aim at gaining a comprehensive understanding of the connection between these two crucial axes of cancer.

Mathematical models of EMT
Computational models developed for EMT can be categorized broadly into two classes: mechanism-based models, and data-based models. While the former adopts a 'bottom-up' approach and focus on elucidating the properties of molecular networks identified experimentally, the latter adopt a 'top-down' approach starting with high-dimensional data and aim to reverse engineer the networks and/or trace the trajectories of these transitions using statistical methods.

Decoding the dynamics of cellular transitions: Mechanism-based models of EMT
The first set of mechanism-based models for EMT regulation -developed independently by two groups -focused on a small set of nodes, and captured the dynamical features emerging from the interconnections among those nodes (Fig. 1A, left). These models included the EMT-suppressing microRNA families miR-34, miR-200 and the families of EMT-TFs ZEB and SNAIL [25,26]. Both models predicted that this network can be tristable, and could give rise to a hybrid epithelial/ mesenchymal (E/M) phenotype in addition to epithelial and mesenchymal phenotypes (Fig. 1A, right) [25,26]. These models suggested that more than one EMT phenotype can be accessible to a cell due to underlying multistability, hence giving rise to sub-populations of epithelial, hybrid E/M and mesenchymal cells in a genetically identical population. This phenomenon was observed and later characterized in detail in multiple cancer cell lines [5,[27][28][29]. However due to different modeling approaches, these models differed on the dynamics of attaining this hybrid E/M phenotype. Experimental support for both these models has been observed [27,30], highlighting the heterogeneity and multiplicity of hybrid E/M phenotype(s) present in different cell lines.
Further follow-up work has identified several intracellular phenotypic stability factors (PSFs) that can stabilize a hybrid E/M phenotype, including OVOL2, GRHL2, Np63α and NRF2 [31][32][33][34][35]. Their role as PSFs have also been validated experimentally in vitro and in vivo [32,[34][35][36]. Moreover, higher levels of these PSFs correlated with a worse patient survival, emphasizing the clinical implications of hybrid E/M phenotypes [4,22]. Among those, NRF2 has been recently proposed to be maximally expressed, and therefore being a hallmark, in these hybrid E/M phenotypes [35]. Different energy landscape approaches have been also developed for the aforementioned EMT circuit [37] as well as for larger gene regulatory circuits [38] to compute the transition rates between EMT states. These aim to predict the relative abundance of different phenotypes in an isogenic population.
EMT can be also induced by biochemical signals coming from neighboring cells. Boareto et al. [39] elucidated the connection between EMT and the Notch signaling pathway, a cell-cell contactbased, evolutionary conserved signaling mechanism implicated in different hallmarks of cancer including angiogenesis and therapy resistance. The model predicted that Notch-Jagged signaling, but not Notch-Delta signaling, among cells can foster the formation of clusters of hybrid E/M cells by promoting a similar hybrid E/M phenotype in neighboring cells [39]. This prediction was further supported by gene expression analysis of CTC clusters exhibiting higher levels of Jagged, compared to single CTCs in inflammatory breast cancer [33]. Therefore, a hybrid E/M phenotype can be stabilized not just by intracellular PSFs directly coupled to EMT core circuit, but also by cell-cell signaling. As another example, Bocci et al. [40] predicted that Numb -an inhibitor of Notch signaling -can also stabilize a hybrid E/M phenotype; this prediction was validated experimentally in multiple independent studies [40,41].
As the network grows in size (such as going from Fig 1A to Fig. 1B), identifying parameters becomes more and more challenging. Thus, computational models that have focused on larger networks have typically been simulated using Boolean modeling approaches, where the state of gene expression is either on (active) or off (inactive). Cohen et al [42] developed a Boolean network to evaluate the combinatorial effect of different mutations on EMT and metastatic potential using transcriptome data from TGF--induced EMT. Similarly, Steinway et al [43] constructed a circuit for TGF--induced EMT using data from hepatocellular carcinoma (HCC). Their model predicts the activation of several pathways during EMT such as Sonic Hedgehog and Wnt. Following up, the authors further showed that certain perturbations could give rise to one or more hybrid E/M states, and identified key circuitry components as possible targets to inhibit TGF--driven EMT [44]. Font-Clos et al [45] recently constructed a genetic network that describes both EMT and its reverse, MET. An energy landscape approach showed two main attractors, or stable states, corresponding to epithelial and mesenchymal phenotypes, and multiple local minima, or relatively less stable states, corresponding to hybrid E/M phenotypes (Fig. 1C). They further mapped RNA-seq data from both lung adenocarcinoma and embryonic differentiation during EMT/MET and compared it to the predicted phenotype profiles, to validate the existence of multiple different intermediate E/M states.
In an attempt to combine the advantages of both continuous small-scale models and Boolean large-scale models, Huang et al [46] devised a Random Circuit Perturbation algorithm (RACIPE) where the expression levels of genes are continuous, but the parameters for all regulatory links are randomly chosen within a biologically relevant range. RACIPE generates an ensemble of mathematical models, each with different sets of parameters, and identify the robust dynamical states emerging from a given network topology. Applying RACIPE to an EMT circuit composed of 9 microRNAs and 13 TFs (Fig. 1B) highlighted two different hybrid E/M states [46] that could be stabilized further by stochasticity or noise [47].

Figure 1. Mathematical models that characterize the landscape of cellular plasticity mediated via EMT and CSCs. (A)
Left: a gene regulatory circuit for EMT proposed by Hong et al. [34]. Right: a bifurcation diagram of ZEB mRNA as a function of EMT-TF SNAIL adapted from Lu et al. [25] [65]. (E) Left: a core gene regulatory circuit including regulation of EMT via the miR-200/ZEB axis and stemness via the LIN28/let-7 axis proposed by Jolly et al. [12]. Right: varying the coupling strength between the EMT circuit and the stemness circuit shifts the 'stemness window' along the EMT axis in the model of Jolly et al. [12]. (F) The EMT score of different cancer stem cell lines adapted from Bocci et al. [13] shows the spread of CSC properties along the EMT spectrum.  [49]. These regulators correlate with a worse clinical outcome and their knockdown can prevent TGFβ-driven EMT [49].
Two different groups developed methods to analyze gene expression data of a certain cell line or tumor and calculate an 'EMT score' predicting the positioning along the EMT spectrum. Tan [51]. Another, more rigorous analysis was proposed by Leggett et al. [52] that relies on single cell microscopy to classify cells as epithelial or mesenchymal with high precision during TGF-β driven EMT [52].
As the connections between molecular and morphological traits of EMT continue to be explored in detail [53], a synergistic crosstalk among the computational models described above and their integration with experimental data can provide novel and crucial insights into EMT dynamics.

Modeling the CSC fraction during tumor progression
An important direction where mathematical approaches have offered significant insights into the CSC dynamics and its relationship with tumor progression is a set of population dynamics models that aim at understanding the temporal dynamics and the mechanisms regulating the CSC fraction, the fraction of cells with stem cell properties in a tumor [54,55]. Dhawan et al. [56] considered two compartments, the CSC-like cells and the non CSC-like cells, to elucidate the increased plasticity observed in human mammary epithelial cells under hypoxia. In their model, individual cells can both differentiate (from CSC to non-CSC) and dedifferentiate (from non-CSC to CSC). Integrating their model with gene expression analysis, the authors showed that hypoxia generates a shift toward a more stem-like population and increases EMT features [56].
The role of cell dedifferentiation to a stem-like state has been also investigated by Jilkine et al. [57] via a hybrid model that describes deterministically the development of the differentiated cell population but considers stochastic accumulation of mutations to better describe the small CSC population. The authors concluded that dedifferentiation to a stem-like state can speed up tumor progression by enlarging the CSC population [57].
Among the different possible mechanism for dedifferentiation, metabolic reprogramming is especially frequent in the context of cancer [58]. Liu et al. [59] devised a probabilistic framework to specifically investigate metabolic reprogramming that converts somatic cells into pluripotent stem cells. Zhou et al. [60] developed a population model of tumor growth that integrates the differential growth rate of CSCs and differentiated cells as well as transitions among cell phenotypes. In particular, phenotype switching contribute to maintaining an equilibrium ratio of cell sub-populations [60]. Extending this idea, Wang et al. [61] proposed a population model that combines hierarchical organization (irreversible loss of stem traits upon differentiation) and stochastic switching (stemness can be gained by switching to a stem-like state). In this model, CSCs can self-renew (symmetric division in two CSCs), differentiate (symmetric division into two nonstem cells) and asymmetrically divide into a CSC and a daughter differentiated cell. Additionally, differentiated cells can proliferate (symmetric division) but also switch to the progenitor CSC state [61]. The combination of hierarchical and stochastic processes can reproduce the CSC/differentiated cell fraction observed in a human colon cancer cell population [61]. A similar idea has been proposed by Zhou et al. [62] as well, showing that phenotypic plasticity allowing back and forth transitions between stem-like and non-stem states is crucial to establish an equilibrium cell fraction of CSCs [62].
A different approach to model CSC-driven tumor progression was proposed by Poleszczuk et al. [63]. They proposed an agent-based model where CSCs can gain migratory traits based on acquiring stochastic mutations. Such approach enables to simulate the spatiotemporal dynamics of the cancer cell population, and investigates cell heterogeneity that arises during tumor development due to mutations. In this model, CSC can divide symmetrically or asymmetrically and also have a migration potential that translates into discrete movements on a two-dimensional lattice [63], reminiscent of the idea of 'migrating cancer stem cell' proposed by Brabletz [14].

CSC, tumor progression and therapy: from modeling to the clinic
A subset of recent models have focused their attention toward identify optimal therapy schedules for the cancer cell [64]. In this context, CSCs are considered as important target as they can be therapy-resistant and therefore drive tumor progression and relapse. For instance, the model developed by Dhawan et al. [56] discussed earlier can be generalized for the context of understanding drug tolerance by introducing one or more additional cell sub-populations exhibiting resistance to different treatments [65]. Specifically, the authors show via integrating in vivo experiments and mathematical modeling that chemotherapy can change the rates of conversions among different phenotypes and promote a chemotherapy-tolerant state (Fig. 1D) [65].
A more data-driven approach aims to relate the CSC population with tumor progression and response to therapy. For instance, Werner et al. [66] proposed a computational method to quantify the fraction of tumor-initiating cells (i.e. CSCs) by analyzing the tumor's macroscopic growth rate as a function of time. This patient-specific method can be applicable to many types of tumors and provides an estimate of the CSC fraction that can be used in a clinical setting to rationalize the optimal therapy [66]. Zhou et al. [67] applied a statistical approach to compute the transition rate between CSC and differentiated cells in colon cancer cells and showed phenotypic plasticity with back and forth transitions [67]. Furthermore, Yu et al. [68] gathered the differential response of CSCs and differentiated cells to radiotherapy for different tumor types including glioblastoma, lung, prostate, and breast cancer, and fitted this tumor-specific information in a stochastic mathematical model to explain the different inter-tumor responses to radiation therapy [68].
Not all models of cancer cell-therapy interplay need to employ a population dynamics approach. Instead, Chen et al. [69] used an energy landscape approach to investigate the transitions of breast cancer cells between phenotypes which are sensitive, hypersensitive or independent to hormone therapy regulated by the ER signaling network. The authors implemented different treatment strategies including sequential treatment (multiple drugs) and intermittent treatment (alternation of treatment and 'holiday' periods) [69]. The effects of continuous vs. intermittent treatments was also explored in the context of prostate cancer, where a small-scale model predicted that cells could oscillate between a therapy-sensitive and a therapy-recalcitrant phenotype [70]. The authors further modeled different hormonal treatments for prostate cancer that were predicted to synchronize oscillations among different cells, thus restricting population heterogeneity [71].
In general, we feel that there is much room for progress in constructing mechanistic models for CSC-driven tumor progression and the emergence of drug-resistant phenotypes, because most models related to CSCs have focused on identifying causes underlying varying fractions of tumor population that can behave as CSCs [72]. In this direction, Nazari et al. [73] recently proposed a mathematical model for the role of inflammatory cytokines in mediating CSC-driven tumor growth. This model couples the ligand-receptor interaction at the molecular scale with CSC selfrenewal and proliferation at the cellular level [73], and could reproduce the decreased tumor volume in mouse models with knockdown of IL-6, an inflammatory cytokine that can act as driver of tumor growth.

Towards an integrated understanding of EMT and CSC
Aside from separate models for EMT and CSC dynamics as discussed above, multiple computational models have investigated the connection between EMT and CSC. Turner et al. [74] interrogated the connection between EMT and CSC through a phenomenological population-level model with two possible scenarios of EMT-mediated enrichment of the CSC population. First, cells can dedifferentiate back to a CSC state while undergoing EMT. Secondly, EMT increases the probability of symmetric, self-renewal division of cells that are already stem-like [74]. The authors used the model to fit the experimental data on CSC fraction and mammosphere expansion, indicating that both processes may play an important role in supporting cancer progression [74]. Later, Gupta et al [75] showed that breast cancer cells can exist in different sub-populations with varying functional attributes: luminal, basal and stem-like. They demonstrated that the overall population, when perturbed, re-establishes an equilibrium fraction of the three cell phenotypes. This robustness could be explained by a population model where cells can undergo stochastic phenotypic transitions between the three different states [75]. Moreover, the stem-like cell line SUM 149 has been shown to exhibit the traits of a hybrid E/M phenotype, hence suggesting a possible correlation between a partial EMT and stemness [33].
A different type of example, focusing on spatially resolving the population structure, is the multiscale model proposed by Sfakianakis et al. [76] that couples CSC and EMT to describe the extracellular matrix invasion by tumor cells. This phenomenological model considers EMT driven by growth factors (individual cell scale) as a binary switch between an epithelial-like and a mesenchymal-like phenotype, and the population dynamics and growth of the tumor mass (multicell scale). Note however that the models discussed so far proposed mechanisms for CSC-driven tumor progression and maintenance of the CSC fraction, but did not provide a molecular rationale for the acquisition of CSC traits.
Li and Wang [38] reconstructed a core gene regulatory circuit with relevant players determining CSC properties such as miR-145 and OCT4, and core regulators of EMT -miR-200 and ZEB. The authors applied an energy landscape approach to predict the co-existence of multiple cellular phenotypes. In their model, a cell either assumes a 'normal' state or a 'cancer' state, both of which could or could not exhibit stem-like traits, thus a total of four possible cell phenotypes [38]. In this framework, p53 represents a degree of cancerization and ZEB represents a degree of stemness. In other words, the predicted 'normal stem cell' and 'CSC' states highly express ZEB, hence implicitly suggesting that stemness is gained along with EMT [38].
Finally, the models developed by Jolly and colleagues explicitly proposed a mechanism-based rationale to elucidate the connection between EMT and CSC: the stemness circuit comprising LIN28, let-7, and OCT4 is connected to the EMT circuit already discussed by Lu et al [25] (Fig. 1E, left). The CSC phenotype was defined as a state with intermediate levels of OCT4 that have been shown to correlate with stem-like traits [77,78]. These models proposed that a CSC phenotype is highly correlated with a hybrid E/M phenotype [79], but intracellular factors such as OVOL [12] or cell-cell communication via Notch signaling [13] could move the predicted 'stemness window' toward the epithelial or mesenchymal end of the EMT spectrum (Fig. 1E, right). Experimental evidence for this dynamic 'stemness window' concept was provided by Bocci et al [13] by computing the 'EMT metric score' [23], of different human CSC lines using publicly available datasets and showing that CSC traits can be scattered along the EMT spectrum based on contextspecific activation of signaling pathways, therefore resulting in epithelial, hybrid E/M and mesenchymal CSC (Fig. 1F) [13]. Furthermore, this model proposed a strong overlap of a hybrid E/M phenotype, CSC properties and enhanced Notch-Jagged signaling [13], a pathway implicated in both drug resistance and in clusters of CTCs, the key drivers of metastasis [33]. Thus, this work opens a window to understanding spatial arrangement of cells with stem-like traits, as seen experimentally [80], since receptor-ligand interactions via Notch that can generate different types of spatial patterns including lateral inhibition (i.e. neighbor cells have opposite phenotypes) or lateral induction (i.e. neighbor cells have similar phenotypes) [81,82].

Conclusion
EMT and CSC represent two crucial biological axes that bolster tumor progression, metastasis and tumor relapse [2,4]. Although the molecular details of multiple steps of tumor development continue to be identified, it is largely accepted that EMT often plays a crucial role in regulating epigenetic, morphological and functional cell properties during tumor progression and metastasis formation [2,4]. Similarly, it is well accepted that the acquisition of stem-like properties by cancer cells potentiate tumor maturation and resist various treatments, driving tumor relapse. Only recently, we have been gaining insights into how, when and where these two dynamic processes can influence one another (Fig. 2). In this context, mathematical modeling has proven itself as a potent tool to interpret existing data and formulate new predictions that can be tested experimentally. In the context of EMT, mechanism-based computational models have suggested that intermediate cell states enabling hybrid phenotypes with mixed epithelial (E) and mesenchymal (M) characters, as opposed to a binary E-M switch scenario [2,4]. Novel in vivo and in vitro analysis recently highlighted the existence of such hybrid states able to co-express E and M markers and possessing mixed morphological traits of cell-adhesion and motility [18,19], and have highlighted their enhanced metastatic potential [83]. The next crucial steps will include a more comprehensive attempt to integrate data-based models, mechanism-based models, and time-course and singlecell experimental data, to formulate a more quantitative characterization of these malignant hybrid E/M state(s).

Figure 2. Interconnections between EMT and CSC axes. The spectrum of EMT states can range from epithelial (E) to mesenchymal (M) phenotypes, with a variable number of partial E/M states (x-axis). Both cancer stem cells and differentiated cancer cells can assume different states depending on genetic and epigenetic factors (y-axis). Horizontal and vertical transitions represent
In the context of CSC, one set of models consider the dynamics of a CSC population employing the tools of population dynamics and agent-based modeling. Such class of models could provide predictions about CSC fraction or population dynamics under perturbations, hence potentially providing strategies for containing CSC-driven tumor progression. Additionally, coupling mathematical modeling with clinical data of therapy response enables predictive tools that can shed light on the CSC-therapy interplay. Such models can provide information on, for instance, adaptive response, differential drug sensitivity, or phenotypic plasticity in a cancer cell population.
Recent experimental observations have led to a class of models that can offer insights into coupling between EMT in cancer cells with the acquisition of stem-like properties. A first set of models relates phenomenologically the acquisition of stem traits with the EMT process, hence explaining how CSC-EMT interaction can support tumor progression and maintain a certain fraction of different cancer cell phenotypes. Moreover, a second class of models investigated the coupling between EMT and CSC at the level of gene regulatory networks, showing a correlation between the cell phenotypes enabled by an EMT regulatory circuit and CSC. A common feature across these models is envisioning the acquisition of stemness as a dynamical process correlated to EMT [12,13,38], which can be reversed by interactions with the tumor microenvironment and perturbations of a cancer cell population [13]. Crucially, recent mathematical modeling and experiments suggested a correlation among hybrid E/M states and stem cell properties [5,12,13,18,19]. However, CSC traits are not exclusively observed in intermediate states, and the crosstalk between tumor, micro-environment and therapies is likely to play a major role in modulating the plasticity properties of cancer cells (Fig. 2), as shown by recent experiments highlighting subsets of CSCs in multiple cancer types [20,80].
Considered together, these computational models developed for EMT, CSC, or their interplay have contributed not only in deciphering the mechanisms underlying specific experimental observations, but also have driven the next set of experiments by generating testable predictions. Such bidirectional crosstalk can significantly accelerate our goal of understanding and consequently targeting these processes for therapeutic benefit.