Population Dynamics of Epithelial-Mesenchymal Heterogeneity in Cancer Cells

Phenotypic heterogeneity is a hallmark of aggressive cancer behaviour and a clinical challenge. Despite much characterisation of this heterogeneity at a multi-omics level in many cancers, we have a limited understanding of how this heterogeneity emerges spontaneously in an isogenic cell population. Some longitudinal observations of dynamics in epithelial-mesenchymal heterogeneity, a canonical example of phenotypic heterogeneity, have offered us opportunities to quantify the rates of phenotypic switching that may drive such heterogeneity. Here, we offer a mathematical modeling framework that explains the salient features of population dynamics noted in PMC42-LA cells: (a) predominance of EpCAMhigh subpopulation, (b) re-establishment of parental distributions from the EpCAMhigh and EpCAMlow subpopulations, and (c) enhanced heterogeneity in clonal populations established from individual cells. Our framework proposes that fluctuations or noise in content duplication and partitioning of SNAIL—an EMT-inducing transcription factor—during cell division can explain spontaneous phenotypic switching and consequent dynamic heterogeneity in PMC42-LA cells observed experimentally at both single-cell and bulk level analysis. Together, we propose that asymmetric cell division can be a potential mechanism for phenotypic heterogeneity.


Introduction
Intra-tumor heterogeneity is a major roadblock that thwarts multiple therapeutic approaches in the clinic [1]. It has earlier been largely thought of as existing at a genomic level, i.e., co-existence of many sub-clonal populations. Single-cell genomic analysis has helped construct the lineage trees mirroring clonal evolution [2]. However, recent preclinical (in silico, in vitro, in vivo) and clinical observations have emphasised that besides genetic heterogeneity, tumors exhibit substantial non-genetic heterogeneity as well, often referred to as phenotypic heterogeneity [3][4][5]. Non-genetic heterogeneity can facilitate 'bet hedging' in a cancer cell population, thus enhancing its fitness under stressed conditions (immune attack, targeted therapy, etc.) and enabling the survival of subpopulations that can eventually drive tumor relapse and/or metastasis [6,7]. Therefore, identifying the mechanisms underlying non-genetic heterogeneity is of fundamental importance.
A canonical example of intra-tumor phenotypic heterogeneity is along the epithelialmesenchymal axis. Epithelial-Mesenchymal Transition (EMT) and its reverse Mesenchymal-Epithelial Transition (MET) were initially considered as binary processes, but recent investigations across carcinomas, especially those at a single-cell level, have demonstrated that cancer cells can display many hybrid epithelial/mesenchymal (E/M) phenotypes

Asymmetric Distribution of Molecular Content on Cell Division
Following the method proposed earlier [24], we consider fluctuations in the levels of cellular content during its inheritance by daughter cells on cell division. These fluctuations arise due to both imperfect duplication during the cell cycle and later asymmetric partitioning to the daughter cells. We propose these fluctuations to be proportional to the amount of the molecular content available in the dividing parent cell itself. Considering SN AIL par 0 denotes SNAIL level in a cell right after its division. Now, during the cell cycle SNAIL content will approximately double, so that right before the next cell division we can write: where, η 1 is a stochastic scaling factor that determines the fluctuation due to imperfect molecule duplication.
Next, when a parent cell partitions its molecular content to two daughter cells during cell division, SNAIL levels in each daughter can be specified as: where, η 2 another random scaling factor determining the fluctuation in SNAIL levels due to partitioning at the time of cell division. We consider stochastic scaling factors η 1 and η 2 to be two independent normally distributed random variables with zero means and η 1 and η 2 as standard deviations, i.e., η 1 = η 1 N 1 (0, 1) and η 2 = η 2 N 2 (0, 1) where, N i (0, 1), i = 1, 2 represents a standard normal random variable. Hereafter, η 1 and η 2 are referred to as scaling factors for noise in SNAIL molecules' duplications & partitioning, respectively. Thus, Equations (2) and (3) can be rewritten as: Equations (4) and (5) are used to assign SNAIL values to the daughter cells when a cell division happens. Further, the same equations were used when stochastics effects were included in the other players of the EMT network (ZEB, mZEB and miR200) at the time of cell division.

Dynamics of Core EMT Regulatory Network
The dynamics of a core regulatory circuit involving interaction in canonical epithelial (miR200) and mesenchymal (mRNA ZEB and ZEB protein) markers was modelled to explain EMT and MET, based on SNAIL levels [22]. miR200 and ZEB (mRNA and protein) mutually repress each other, and SNAIL suppresses miR200 levels and activates ZEB at the mRNA level. The steady state response of the circuit was analysed for a relevant biological parameter set, which gave a bifurcation diagram showing distinct possible stable ranges of ZEB and miR200 based on SNAIL levels as shown in Figure 1D. The ordinary differential equations describing the regulation dynamics and model parameters have been described in Tables 1 and 2 [22]. The system's ODEs are listed below: [·] represents the concentration of a molecular species within a cell. H is the shifted Hill function.
The functions Y µ , Y m , and L describe the post-transcriptional regulation of mRNA activity by micro-RNAs and have been described in [22].
Here, µ is the concentration of the micro-RNA and n is the number of micro-RNA binding sites on the mRNA. For the inhibition of ZEB1 mRNA by miR-200, n = 6 and µ 0 = µ 0 200 . The values of all kinetic parameters are listed in Tables 1 and 2.     Each cell in the system is represented by a set of four variables that hold the levels of miR200, mRNA ZEB, ZEB protein and SNAIL protein for that cell. Random SNAIL values are sampled from a log-normal distribution with median 200 × 10 3 and coefficient of variance 1 and all possible stable states corresponding to that SNAIL value are used to initialize the cells' variables. For example, if sampled SNAIL = 200 K and as for this value all three phenotype-E, E/M and M-are stable. So, three cells are initialised with steady state values of all variables corresponding to each phenotype. Initialization of a cell from a phenotypic state is stopped when its required count in the population is achieved.

Avg. Birth and Death Rate of Cells
In the cell population, the division rate of cells of a particular phenotype follows the logistic equation shown below: And the death rate of cells of a particular phenotype is as follows: : average death rate of an individual cell N phen : total cells of a phenotype N tot : total cells in the population K: Carrying capacity of the system

Population Growth
The population growth is simulated using Gillespie's Stochastic Simulation (SSA) algorithm [25], where six events-three division and three death events for each phenotype are considered. The propensity of occurrence of an event is determined by its average rate as described above. The SSA algorithm tells what the next event will be and at what time point. Now, if the next event is a division of a cell of the E phenotype and will occur at t 1 time, then a cell is uniformly sampled from the pool of cells of that phenotype, and its molecular levels are updated using ODEs for the time gap (t 1 -t 0 ), where t 0 is the time point of last most recent event. Then, a new cell is initialised in the population with molecular levels same as that of the parent E cell, but with perturbed SNAIL levels. Similarly, the parent cell SNAIL levels are perturbed to account for the second daughter cell on division. For a cell death of a phenotype, a cell is uniformly sampled from the pool of cells of that phenotype, and it is erased from the population. Molecular levels of all the other undivided/unaffected cells are updated using ODEs for the time gap (t 1 -t 0 ).

Cell Doubling Quantification
Images for PMC42-LA cells were captured on PhaseFocus LiveCyte Image Scanner (Phase Focus, Sheffield, UK) with 10× magnification; individual images were captured every 11 min for a span of 48 h. Imaging selected regions of interest (ROI) were 750 × 750 µm. Sixty individual selected cells were randomly selected and then manually tracked from cytokinesis of a cell to two daughter cells to next cytokinesis to determine the exact cell doubling time.

Dominance of Epithelial Cells in the Population over Time Irrespective of Initial Distribution
Here, we have developed a population dynamics framework to explain the emergence of epithelial-mesenchymal heterogeneity in a given population, and the contribution of spontaneous state switching in enabling this heterogeneity ( Figure 1A). Specifically, we consider phenotypic switching to occur during cell division ( Figure 1B), where two factors can contribute to a daughter cell having a phenotype different than its parent cell: noise or fluctuations in (i) content duplication and that in (ii) partitioning of biomolecules, particularly in SNAIL (depicted by f (SN AIL par , η 1 , η 2 )) ( Figure 1C). (For more information on formalism used to include content duplication and partitioning noise, please refer to Methods Section 2.1). Each cell contains the ZEB1/miR-200 feedback loop driven by SNAIL, and the levels of these molecules define the state of each cell. Depending on the levels of SNAIL, cells may acquire a phenotype among all the stable ones, as shown in the bifurcation diagram ( Figure 1D) [22]. At SNAIL = 150 K molecules, all cells can adopt only an epithelial state (lower blue curve in Figure 1D); at SNAIL = 200 K molecules, a cell can acquire any of the three states-epithelial (E; lower blue curve), mesenchymal (M; top blue curve) or hybrid E/M (middle blue curve), while at SNAIL = 250 K molecules, all cells adopt a mesenchymal state (top blue curve in Figure 1D). Thus, asymmetry in content duplication and/or partitioning of SNAIL levels can alter the SNAIL values sufficiently enough to allow a phenotypic switch. For instance, if one daughter cell has SNAIL = 250 K for a parent cell with SNAIL = 200 K, then the daughter cell will be mesenchymal in nature irrespective of the phenotype of the parent cell (E, hybrid E/M or M). We have also implemented in silico passaging to mimic the experimental protocol for conducting these experiments, where 10% of the cell population is passaged maintaining the distribution of cells in different phenotypes when the entire population reaches 80% of its carrying capacity ( Figure 1E). Further, the division rate of each subpopulation of cells is considered to follow the logistic growth rate, whereas the death rate is directly proportional to the subpopulation size ( Figure 1F). Together, these factors are incorporated in a population dynamics model including cell division which may be accompanied by a phenotypic switch (please see Methods Sections 2.2 and 2.3 for more details about the population dynamics model).
Using this framework, we investigated how the population distribution emerged over time as we started with different initial fractions, and whether we can recapitulate the dominance of epithelial (EpCAM high ) subpopulation over a mesenchymal (EpCAM low ) one as seen experimentally (Figure 2A) [20]. We first experimentally quantified the doubling time of PMC42-LA cells to be 22.67 ± 2.77 h ( In these simulations, we considered an average doubling time of the population to be 20 h, η 1 (scaling factor for noise due to SNAIL molecules' duplication) = 0.2 and η 2 (scaling factor for noise due to SNAIL molecules' partitioning) = 0.1, and tracked the population distribution as a function of time. This choice of values represents typical coefficient of variance values in protein levels reported in human H1299 lung carcinoma cells [26]. We observed that all the initial fractions converged to an epithelial dominant population over a period of 16 weeks ( Figure 2B), i.e., greater than or equal to 80% population being E (EpCAM high ). Particularly, for a mixed initial fraction (E:E/M:M = 0.4:0.2:0.4), most of the hybrid E/M cells switch phenotype either to E or M within 2 weeks of time, after which population is mostly comprised of E and M cells, which eventually converges to an epithelial dominant one. Concomitantly, there is also a shift in the distribution of SNAIL levels, such that the range of SNAIL values observed tend to correspond to an epithelial phenotype as well by week 16 and 32, as compared to week 0 ( Figure 2C), thereby explaining a gain in the epithelial-dominated subpopulation as seen experimen-tally, i.e., the EpCAM high subpopulation constituting the majority of the PMC42-LA cell line [20]. Further experiments revealed that this population distribution can be reproduced by the FACS-sorted EpCAM high and EpCAM low subpopulations when cultured individually, thus reminiscent of our simulations showing the asymptotic dominance of epithelial subpopulation irrespective of initial phenotypic distributions.
Over longer simulation times in our model, the dominance of the epithelial fraction grew even stronger ( Figure S1A). Thus, while our model encapsulates the dominance of EpCAM high subpopulation in PMC42-LA cells, it cannot accurately reproduce the experimentally observed 80:20 EpCAM high : EpCAM low ratio. This lacuna indicates the role of various important factors (both cell-autonomous and non-cell-autonomous: chromatin status and cellular communication, for instance, respectively) which can influence stochastic fluctuations during cell division induced spontaneous switching, thus altering this ratio. Nonetheless, our simple phenomenological model can reproduce salient features of population dynamics reported in the PMC42-LA cell line [20].
To assess how fluctuations in other players of the EMT network (miR-200, ZEB) influence population dynamics, we introduced stochasticity in their content duplication and partitioning during cell division, rather than just in SNAIL levels. The population dynamics for this scenario, using the parameters and initial fractions described above, gave qualitatively similar results of epithelial dominance over mesenchymal, though the time taken to gain this dominance was shorter in this case (Figure S1B-C), indicating that additional noise can accelerate the system dynamics. These results suggest that accounting for asymmetry in the levels of SNAIL is sufficient to capture the qualitative population dynamics for PMC42-LA cells.
The dominance of the epithelial phenotype in a population over time irrespective of initial fractions of phenotypes points towards the possibility that hybrid E/M and mesenchymal cells switch more frequently to epithelial as compared to epithelial states switching to hybrid E/M and mesenchymal. This one-sided higher switching rate can be explained by the multiplicative nature of noise considered (lower levels of SNAIL in epithelial cells invoke further fewer fluctuations during division) and have been quantified later.

Time to Attain Epithelial Dominance Depends on Initial Fractions, Doubling Times of Phenotypes and the Extent of Stochastic Fluctuations during Cell Division
Upon simulating the population dynamics starting with purely E and purely M phenotypes, we noticed differences in the mean epithelial (E) fraction at week 16 ( Figure 3A). When starting with the purely M phenotype, it took 16 weeks to arrive at an epithelialdominant population as compared to starting with the hybrid E/M (8-12 weeks) or fully E ones. This trend raised the possibility that while initial phenotypic distribution may not alter the steady state itself, it can change the time taken to arrive at it. For the scenario starting with purely hybrid E/M cells, within 4 weeks the population structure had approximately 60-70% epithelial cells, thus its dynamics post the 4 week timepoint is understandably similar to that seen for E:E/M:M = 0.7:0.1:0.2 scenario (compare middle panel in Figure 3A with the top left panel in Figure 2B).
Besides initial phenotypic distribution, another factor that can impact the population dynamics is average doubling time (DT). So far, we considered the same DT for E, E/M, and M phenotypes. However, experimental evidence suggests the slowing down of proliferation rate of cells undergoing EMT [20,27]. Thus, we considered the case of increased DT during EMT, by keeping the average DT of E/M and M phenotypes as 1.5, 2, and 2.5 times more than that of the E phenotype. We observed that the population maintained its epithelial dominance, and converged to a stable phenotypic distribution faster than in the case when all cells doubled at the same rate ( Figure 3B), irrespective of the initial condition. This trend can be explained by a higher resilience of E cells to switch to a hybrid E/M or M phenotype during cell division, now coupled with their higher proliferation rate, thereby offering the epithelial subpopulation an additional advantage to amplify their population fraction. Importantly, this trend was already seen at the DT of hybrid E/M and M cells being 1.5 times that of the E cells, hence indicating that a 50% increase in doubling time for cells undergoing EMT may be sufficient in influencing the population structure. The slight initial decrease in the epithelial fraction noticed for the purely E case ( Figure 3A,B) can be explained by appreciating that SNAIL levels for the initial cell population were sampled from a log-normal distribution, whose median was centred on the SNAIL level where all phenotypes were stable (tristable region in Figure 1D); therefore, the cells were highly susceptible to undergo phenotypic switching within the first few cell divisions. Despite this initial dip, an epithelial dominant population emerged eventually. Next, we investigated how the extent of stochastic fluctuations in SNAIL molecules being duplicated and partitioned (η 1 and η 2 respectively) influenced phenotypic distribution over time. When we varied η 2 , while maintaining the values of η 1 = 0.2 and average population DT = 20 h, we noticed that for the mesenchymal dominated initial fraction (E:E/M:M = 0:0.2:0.8), the fraction of epithelial cells was higher for a higher η 2 value for the same time point (Figure 3C, left). However, not much observable effect on this fraction was noticed when starting with an epithelial dominated population ( Figure 3C, right). Similar observations were made when we varied η 1 instead of varying η 2 ( Figure 3D). Thus, amplifying fluctuations in either duplication or partitioning of SNAIL molecules seemed to enhance the chance of phenotypic switching for an M cell much more than for an E cell. When we accounted for heterogeneity in average DT along with increasing fluctuations in SNAIL levels during cell division, the fast proliferating and relatively stable E cells grew much faster than the slow proliferating and more plastic E/M and M cells, enriching for epithelial cells (Figure S2A-D).
Finally, we analysed how the population dynamics were altered when average DT was increased for all phenotypes, given the experimentally observed average DT can often depend on the confluency of cells in a petri dish. Thus, for this, we simulated population dynamics keeping the average DT of all three phenotypes as 30 h, which led to overall slower dynamics ( Figure S2E). Instead of plotting against absolute time units, we also took the number of cell cycles as the x-axis, whose one unit is the population's average DT. This helped to compare the overall changes in phenotypic fractions between DT of 20 and 30 h scenarios. We found that given an equal number of cell cycles, the changes in the E fractions were similar ( Figure S2F). These observations help us to conclude that even if all cells, on average, divided slower, the population growth and phenotypic switching trajectory would be similar to when cells divided faster when normalised with average DT.

Phenotypic Switching Probability and Rate in Cell Division Events Depends on the Cells' Location on E-M Axis
After characterising the population dynamics at various time points as a function of different model parameters, we wanted to better understand it from a cell division perspective. In our framework, a cell can undergo one of the three division types: (1) symmetric division-when both daughter cells have the same phenotype as the parent cell, (2) asymmetric division-when one daughter has a phenotype different than the parent, and (3) divergent division-when both daughters have a phenotype different than the parent ( Figure 4A). To quantify the probability of cells undergoing one of the three division types, we analysed certain cells occupying possible stable phenotypes spread across the SNAIL ranges ( Figure 4B). Iterating cell division events at a given SNAIL value, we tracked the phenotypes of daughter cells, at specific η 1 and η 2 values, thus calculating different division probabilities over an ensemble of iterations. At η 1 = 0.2, η 2 = 0.1, in SNAIL levels regions where either E or M phenotypes were the only stable state (monostable regions in bifurcation diagram- Figure 4B; SNAIL = 100 K, 300 K), more than 90% events were of symmetric division ( Figure 4C). However, as the SNAIL levels corresponded to a multistable region (SNAIL values = 150 K, 250 K), there was an increasing tendency to undergo asymmetric division, which was higher for M than for E cells. In different bi-stable regions (SNAIL = 189 K for {E, M} and SNAIL = 219 K for {E/M, M}), with further increasing probability of asymmetric division, the divergent division also became more prominent and was the most dominant division for hybrid E/M cells ( Figure 4C). This trend explains the sudden drop in hybrid E/M fraction of population to very low levels within two weeks, when starting from a hybrid dominant population ( Figure 3A). Further, in the tri-stable region (SNAIL = 206 K for {E, E/M, M}), the probabilities for both divergent and asymmetric division were increased ( Figure 4C). Put together, the probability of phenotypic switching at cell division is the highest in the tri-stable region (intermediate SNAIL levels~200 K) and decreases for cells as their corresponding SNAIL levels either increase or decrease. Next, we quantified these probabilities for varying η 1 and η 2 values. While η 1 = 0 resulted in either asymmetric or symmetric division of E and M cells across SNAIL levels (i.e., preventing divergent division), η 2 = 0 leads to only symmetric and divergent divisions for these two phenotypes (i.e., preventing asymmetric division) ( Figure S3A,B). Additionally, higher values of η 1 and η 2 amplify the chances of divergent and asymmetric division, respectively, across SNAIL ranges (increasing η 1 in Figure S3A,C,E and increasing η 2 in Figure S3B,D,F). Thus, η 1 and η 2 -the factors that represent noise during cell division-can alter the probabilities of undergoing symmetric, asymmetric and divergent division types for a cell with a SNAIL level ( Figure S3).
We observed that cells with SNAIL levels well away from the multi-stable phenotypic regions have mostly undergone symmetric division. However, when we started with such a homogenous or largely homogeneous population and tracked the phenotypes of daughter cells over multiple cell cycles, we noted phenotypic switching in which at least one daughter cell took a different phenotype (Figure 3). Thus, we quantified how many cell divisions it took for a cell to give rise to one of the cells in its progeny with a different phenotype than its own. We observed the progeny up to 12 generations/cell cycles. We saw that the cells with SNAIL levels in a multi-stable region switch phenotype within one or two cell cycles ( Figure 4D). We also noticed a skew between the resilience of E and M cells to phenotypic switching in their mono-stable regions, i.e., E cells required more cell cycles to give rise to a non-similar progeny cell than the M cells did (compare the behavior seen at SNAIL = 300 K and SNAIL = 250 K with that at SNAIL = 100 K and SNAIL = 150 K in Figure 4D). This difference may underlie the phenomenon of E cells dominating over E/M and M cells in the population over time. However, this skew vanished in the bi-stable and tri-stable regions, where all three phenotypes became equally susceptible to undergo asymmetric switching within a few generations/cell cycles (SNAIL = 189 K, 206 K, 219 K in Figure 4D).
We also examined how η 1 and η 2 values varied the number of cell cycles over which progeny diversification was observed. An increase in either η 1 and η 2 caused faster switching for all phenotypes of cells, i.e., a smaller number of cell divisions was required, on average, for a cell to give rise to a different phenotype (increasing η 1 in Figure S4A-C and increasing η 2 in Figure S4D-F). However, η 2 contributed more as compared to η 1 (compare Figure S4A-C with Figure S4D-F; quantified in Figures S5 and S6).

Heterogeneity in E Fraction at Initial Time Points among Single Cell Clones
So far, we have focused on population dynamics when starting with an initial cell population; however, heterogeneity has also been observed experimentally in single-cell clones established from cell lines [20,21]. For instance, in distinct single-cell clones established from PMC42-LA, varying distributions of EpCAM high : EpCAM low subpopulations were seen after the initial two passages ( Figure 5A) [20]. We interrogated whether our model could reproduce this heterogeneous behavior.
We performed a population dynamics simulation starting with a single E, E/M and M cell, maintaining η 1 = 0.2, η 2 = 0.1, and average DT as 20 h. We observed heterogeneity in the E fraction at week 4 when multiple such single cell simulation runs were performed ( Figure 5B). Thus, as a proof of principle, our model could recapitulate the experimentally observed heterogeneity in the EpCAM low fraction. In these simulations for single-cell clones, at the week 4 time point, the heterogeneity in the fraction of E cells was the highest when the seeding cell was mesenchymal (M). Among the M clones, the highest E fraction noticed was close to the highest E fraction noticed for single-cell clones established from the E or E/M initial phenotype. However, in M clones, we observed instances where the E fraction was as low as 27% ( Figure 5B).
To examine this heterogeneous behaviour of individual clones more closely, we probed the levels of SNAIL in the seeding (individual) cell for each of these clones. This led us to identify the range of SNAIL levels in the individual cells that were all 'cultured' in silico to develop a clone ( Figure 5C). From this range, we identified representative SNAIL values and independently established single-cell clones from them. Interestingly, the single-cell clones showed heterogeneity in the E fraction at 4 weeks, despite being seeded with the same SNAIL level ( Figure 5D), reminiscent of stochastic effects at lower (cell) numbers. Additionally, as expected, the average E fraction decreases as seeding SNAIL levels are increased (compare the average E fraction at SNAIL = 600 K vs. that at SNAIL = 50 K in Figure 5D). Another feature we noticed is that the clones established from cells with the same initial phenotype (M), but different initial SNAIL levels had varying E fractions at week 4 (compare the average E fraction at SNAIL = 206 K vs. that at SNAIL = 600 K in Figure 5D). This extent of diversity is less when the initial cell belongs to an epithelial or a hybrid E/M phenotype. This difference between the extent of variability noticed can possibly explain why we see more heterogeneity in the E fraction when starting from an initially mesenchymal cell as compared to an initially epithelial or hybrid E/M one ( Figure 5B-D). Finally, when we continued the single-cell (clonal) simulations for a longer duration, we observed a decrease in heterogeneity in the E fractions with time ( Figure 5E). Further, the E fraction for all single cell clones increased overall. This difference in short-term vs. long-term behavior can be possibly rationalised by our population dynamics simulations earlier showing predominance of epithelial phenotypes irrespective of initial phenotypic distributions ( Figures 2B and 3A) if we consider the clonal distribution noticed at week 4 as the initial condition for simulations being continued until week 16 or later.

Discussion
Understanding the molecular mechanisms underlying epithelial-mesenchymal plasticity and heterogeneity can contribute to better therapeutic strategies [28]. These mechanisms can be context-specific with varying degrees of contribution to genetic and/or non-genetic heterogeneity. Epigenetic alterations, for instance, can govern the rate of bidirectional switching among the phenotypes, enabling reversible or irreversible EMT, as well as driving resistance to undergo EMT [29][30][31]. Cell-cell communication through autocrine and/or paracrine signaling with other tumor cells as well as stromal cells can also shape the E-M phenotypic heterogeneity patterns in a population [19,[32][33][34]]. Another contributing factor can be differential activation of many signaling pathways implicated in EMT, thus generating a varied phenotypic repertoire of states in the multidimensional EMT landscape [35][36][37][38][39]. Here, we highlight one other possible reason driving E-M heterogeneity-phenotypic switching due to asymmetric cell division driven by noise in the processes of content duplication and in the partitioning of biomolecules. We investigate the influence of such fluctuations on levels of SNAIL-a driver of miR-200/ZEB feedback loop-during cell division in determining the phenotypic distribution of population, but our framework is applicable to investigate the population dynamics emerging from stochastic partitioning of molecules involved in other multi-stable EMT networks [40,41] as well.
Asymmetric cell division is an evolutionarily conserved mechanism used by prokaryotes, as well as eukaryotes to generate cell-to-cell heterogeneity, and mediate cell-fate decisions [42,43]. This phenomenon has been observed in multiple cancers, particularly relating to cancer stem cell (CSC) differentiation [44][45][46]. Glioma CSCs can differentiate into either glial or neural cells by asymmetric partitioning of CD144 and Numb levels during cell division [45]. Similarly, colon CSCs can have asymmetric partitioning of ALDH1 and CK20 levels, leading to cell differentiation [44]. The frequency of asymmetric division depends on micro-environmental conditions such as growth factors in the case of glioma CSCs, and on regulatory molecules such as Notch and miR-34a in colon CSCs. Consequent heterogeneity in molecular content among cells can lead to different signaling responses, as measured by distinct IFN-γ and oncostatin M response in fibroblasts [47], and can lead to the presence of rare drug-resistant cells, such as those seen in melanoma [48]. Not only biomolecules (RNAs, proteins) but also entire organelles, such as mitochondria and membrane lipids can be asymmetrically partitioned, with implications in cell division rates and/or stemness traits [49,50]. In the EMT literature, to the best of our knowledge, the asymmetric partitioning of cell-fate determinant transcription factors (here, SNAIL) has not yet been reported experimentally, but the frequency of the modes of cell division (symmetric division, asymmetric division, divergent division) has been shown to vary with EMT induction [51]. Although our modeling framework does not yet specifically incorporate molecular mechanisms regulating this phenomenon [42], our results suggest that one possible consequence of fluctuations during cell division can be phenotypic switching and heterogeneity among subpopulations. Recent reports in glioblastoma have demonstrated that asymmetric enrichment of EGFR and p75NTR in a daughter cell during cell division conferred enhanced resistance to the standard-of-care therapies, such as radiation and temozolomide [52]. While we do not yet know about differences, if any, in the drug resistance features of EpCAM high and EpCAM low sub-populations in PMC42-LA cells, the varied drug-resistance features seen in single-cell clones established from PMC42-LA [20] can be a putative outcome of underlying asymmetric cell division. Approximately 10-30% of cells undergoing TGFβ-driven EMT were seen to exhibit asymmetric cell division, as traced by NUMB distribution in daughter cells [51], but whether this asymmetry led to phenotypic switching was not tracked per se. Therefore, our model suggests that blocking cell division can be a possible way to restrict phenotypic plasticity and/or heterogeneity. Preliminary experimental observations made recently support this prediction by our model [51].
Our model can recapitulate the observations for the PMC42-LA system, an epithelialdominant subline. However, what mechanisms may explain phenotypic heterogeneity in a mesenchymal-dominant population, such as EM3 or M clone from SUM149 cell line [21], remains to be investigated further. One factor that can alter the model outcomes is the way noise during cell division is incorporated. We have considered multiplicative noise (fluctuations in SNAIL proportional to its levels); however, our previous effort encapsulating additive noise (constant magnitude of fluctuations in SNAIL, irrespective of its levels) could explain the spontaneous phenotypic switching observations in the prostate cancer PKV cell line [24]. Whether cancer cells exhibit additive or multiplicative noise during cell division remains unknown experimentally. Further, this noise and/or its consequences can be influenced by mutually dependent factors, such as chromatin status and diffusible cytokines [53]. These factors have not yet been explicitly incorporated into our framework.
Despite the above-mentioned limitations, our model recapitulates various observations for the PMC42-LA system: (a) stable dominance of the EpCAM high subpopulation, (b) repopulation of parental distributions starting with only one subpopulation, and (c) enhanced heterogeneity in EpCAM high : EpCAM low ratio of cells in single-cell derived clones. We predict that these single-cell derived clones converge to EpCAM high dominant distribution in longer time-scales, a prediction that remains to be experimentally verified. Thus, we demonstrate that during cell division, stochasticity in content duplication and partitioning of molecules involved in EMT can lead to spontaneous state switching, and therefore, generate non-genetic heterogeneity.
Future efforts are directed towards integrating continuous stochastic fluctuations in EMT drivers with asymmetric cell division which happens at discrete time-steps [54]. Addressing these questions will involve mathematical models that can decode the emergent dynamics at multiple scales-regulatory levels (transcriptional, epigenetic), length (intracellular, non-cell-autonomous effects by cytokines) and time (cell division, chromatin remodelling, stochastic gene expression).

Conclusions
Proportional noise/fluctuations in duplication and distribution of the parent cell SNAIL's level during cell division can possibly explain the experimentally observed population-level dynamics of epithelial-mesenchymal heterogeneity in PMC42-LA cells. Future endeavors would involve incorporating the regulatory mechanism of asymmetric division and analysing the contributions of various stochastic and deterministic (regulatory) processes to it. Another direction would be to consider an integration of asymmetric cell division with other phenotype stabilising mechanisms, such as epigenetics and cellcell communication.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biom12030348/s1; Figure S1: Longer time simulations and similarities between effects of fluctuation in SNAIL and all players in core EMT network during cell division; Figure S2: Temporal changes in E fraction for combinations of avg. doubling time (DT) ratio, η 1 , and η 2 values; and higher cell doubling time; Figure S3: Phenotypic switching probability for various scaling factors (η 1 and η 2 ) values across SNAIL levels; Figure S4: Average cell cycles/generations required for first asymmetric switching for various scaling factors (η 1 and η 2 ) values across SNAIL levels; Figure S5: Statistical analysis of differences in minimum cell cycles for required asymmetric division by a cell of given phenotype and SNAIL level on varying η 1 and keeping η 2 fixed (0.1); Figure S6: Statistical analysis of differences in minimum cell cycles required for asymmetric division for a cell of given phenotype and SNAIL level on keeping η 1 fixed (0.2) and varying η 2 .
Author Contributions: P.J. performed simulations, S.B. performed experiments, E.W.T. and M.K.J. conceived and supervised research. All authors contributed to data analysis and writing of the manuscript. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.