Cancer: More than a geneticist’s Pandora’s box

Despite identical genetic constitution, a cancer cell population can exhibit phenotypic variations termed as non-genetic/non-mutational heterogeneity. Such heterogeneity – a ubiquitous nature of biological systems – has been implicated in metastasis, therapy resistance and tumour relapse. Here, we review the evidence for existence, sources and implications of non-genetic heterogeneity in multiple cancer types. Stochasticity/noise in transcription, protein conformation and/or external microenvironment can underlie such heterogeneity. Moreover, the existence of multiple possible cell states (phenotypes) as a consequence of the emergent dynamics of gene regulatory networks may enable reversible cell-state transitions (phenotypic plasticity) that can facilitate adaptive drug resistance and higher metastatic fitness. Finally, we highlight how computational and mathematical models can drive a better understanding of non-genetic heterogeneity and how a systems-level approach integrating mathematical modeling and in vitro/in vivo experiments can map the diverse phenotypic repertoire and identify therapeutic vulnerabilities of an otherwise clonal cell population.


Introduction
Cancer is thought to originate from an individual normal cell that gains genetic mutation(s) offering it growth advantage over other cells. A clonal population of cells can arise from this 'first' cancer cell; this population has identical genetic composition (Nowell 1976;Greaves and Maley 2012). As time progresses, some of these cells can gain additional mutations, thus leading to sub-clones, some of which can be more fit as compared to others and thus undergo natural selection to become the predominant sub-clone(s). Recent studies across clonal populations in cancer cells as well as other biological contexts (such as microorganisms; Davidson and Surette 2008) have proposed that phenotypic variations exist among genetically identical cells (Brock et al. 2009). These phenotypic variations are referred to as non-genetic heterogeneity (NGH), highlighting their non-genetic/non-mutational origin. Such non-genetic mechanisms can include a combination of various processes such as stochasticity or noise in gene expression (Balázsi et al. 2011), asymmetry in distribution of molecules during cell division (Huh and Paulsson 2011), variability in epigenetic status of cells (Bell and Gilan 2020), promiscuity in protein-protein interaction networks due to disorder in protein structures (Lin et al. 2019) and ability of cells to exhibit multiple phenotypes (multistability) as a result of feedback loops embedded in a gene regulatory network (Evans and Zhang 2020;Hari et al. 2020).
Consequently, cells in a clonal population can have variability in their protein levels, chromatin accessibility status, metabolites, etc., which can eventually manifest as different phenotypes. Non-genetic heterogeneity has been shown to offer survival advantages to cells in fluctuating environments such as drug treatment, hypoxia and nutrient deprivation, by facilitating 'bet-hedging', an evolutionary strategy to achieve fitness across environmental conditions (Veening et al. 2008a, b;Pisco and Huang 2015;van Boxtel et al. 2017;Bell and Gilan 2020;Sahoo et al. 2021a). However, unlike genetic heterogeneity that has been extensively investigated for decades (McGranahan and Swanton 2017), understanding the underlying mechanisms and implications of non-genetic heterogeneity in cancer research is still in its infancy (Marine et al. 2020;Barzgar Barough et al. 2021;Lewis and Kats 2021;Shlyakhtina et al. 2021).
The 'one gene-one enzyme' hypothesis postulated in 1941 (Beadle and Tatum 1941) was among the first to propose a one-to-one mapping between the genotype and phenotype. However, as we now realize, a phenotype is the outcome of extensive cross-talk among many biological factors that form regulatory networks at various length scales and time scales, revealing pleiotropy (Tyler et al. 2009) in the genotype-phenotype map. Regulatory networks in biological systems consist of multiple feedback loops which can produce more than one phenotype depending upon cell-intrinsic (e.g., rates of transcription, translation and protein degradation) and cell-extrinsic (e.g., temperature, oxygen and nutrient availability) factors (Brandman and Meyer 2008). Thus, despite having identical genetic composition, clonal cancer cells can display differences in phenotypic properties such as growth rates (Vega et al. 2004), tumour initiation capabilities (Mani et al. 2008;Pasani et al. 2021) and the ability to evade therapeutic attacks (Sharma et al. 2010;Paek et al. 2016;Shaffer et al. 2017). Intriguingly, cells can switch back and forth among these different phenotypes, driving reversible changes which may or may not be inherited, unlike mutational effects which are irreversible in nature and 'hard-wired' to be passed on to progeny.
In this review, we first discuss the tacitly assumed essentiality and sufficiency of mutations for some aspects of cancer progression and highlight some observations that cannot be completely explained by the somatic mutation theory alone. Next, we provide evidence of non-genetic heterogeneity observed in clonal cancer cell populations and offer mechanistic details for some of its sources. Finally, we discuss the implications of such heterogeneity in various hallmarks of cancer. We also highlight the contribution made by mathematical models to understand non-genetic heterogeneity and suggest that a systems-level understanding integrating theory and experimental observations may provide a better understanding of nongenetic heterogeneity and how this knowledge may be leveraged to restrict aggressive clinical outcomes.

Are mutations necessary for cancer progression?
For many years, DNA mutations have been considered to be the main cause of cancer initiation wherein a mutation in an oncogene or tumour suppressor gene can drive abnormal cell growth. Such mutations are implicitly assumed to be necessary and sufficient for cancer progression, leading to the most popular and widely accepted concept in the field of cancer biology -somatic mutation theory (SMT) -which posits a mutation in a single somatic cell as the first step of cancer (Sonnenschein et al. 2014). This mutation in a cell is considered to be sufficient to perturb its cell cycle regulation, leading to uncontrolled cell proliferation. However, many observations in cancer biology do not appear, at least prima facie, in consensus with SMT, such as spontaneous regression of childhood neuroblastoma without any cytotoxic therapy (Haas et al. 1988), ependymomas in children (MacK et al. 2014), reprogramming of cancerous cells to normal cells when implanted into normal microenvironments (Mintz and Illmensee 1975;McCullough et al. 1997;Kasemeier-Kulesa et al. 2008;Bussard et al. 2010), stromal induction of carcinogenesis in epithelial cells (Barcellos-Hoff and Ravani 2000;Maffini et al. 2004Maffini et al. , 2005Barclay et al. 2005), and the cycling of hepatic cells from cancerous to normal (dormancy) by the dialling up or down of the wild-type MYC oncogene (Shachaf et al. 2004;Shachaf and Felsher 2005).
Thus, an alternative to SMT, although yet not that popular, after SMT has been proposed recently -tissue organization field theory (TOFT) (Soto and Sonnenschein 2011). TOFT considers cancer to be a tissuebased disease (instead of a cell-based disease as considered by SMT) akin to development gone awry. According to TOFT, cancer arises as a result of simultaneous occurrence of two steps: a disturbed interaction between parenchyma and stroma resulting in altered tissue organization and a weaker inhibitory control exerted by tissues over cell proliferation (Sonnenschein and Soto 2015). TOFT is further supported by observations that primary tumours can metastasize to only few specific organs, indicating that the stroma of an organ plays an important role in determining whether or not cancer cells attached to it can start metastatic growth in it (Tarin 2011). This concept is endorsed by recent observations suggesting that cancer cells can carry their own stroma ('soil' as per Paget's 'seed and soil' hypothesis; Paget 1889) for successful metastasis (Duda et al. 2010). While the discussion about similarities, differences and complementarity between SMT and TOFT is beyond the scope of this review, the dynamics of cancer metastasis offers an intriguing scenario to dissect the role of mutations in cancer progression.
Despite tremendous advances in high-throughput single-cell sequencing in the last decade, no unique mutational signature has yet been identified for metastasis (Celià-Terrassa and Kang 2016; Welch and Hurst 2019). However, non-genetic variability due to network architecture between metastasis suppressors and metastasis drivers has been shown to generate a subpopulation of pro-metastatic cells without gaining any additional genetic changes (Lee et al. 2014). This variability further leads to cellular or phenotypic plasticity -the ability of cells to reversibly switch their phenotypes, and has been proposed as a hallmark of metastasis (Bhatia et al. 2020;Biswas and De 2021;Sacchetti et al. 2021). One of the key axes of phenotypic plasticity during metastasis is epithelial-mesenchymal plasticity (EMP), where cells can transition among a range of phenotypes spread over a spectrum ranging from more epithelial (usually more adhesive and less invasive) to more mesenchymal (usually more invasive and having reduced cell-cell adhesion) (Gupta et al. 2019;Pastushenko and Blanpain 2019). EMP can have transcriptional (Yang et al. 2004;Cieply et al. 2012;Ocaña et al. 2012;Roca et al. 2013;Subbalakshmi et al. 2020) and/or epigenetic (Ruscetti et al. 2016;Jia et al. 2019;Nihan et al. 2019;Serresi et al. 2021) regulatory control. Thus, the dynamics of complex interconnected networks at various levels of regulation can dictate the propensity of a cell to slide along the 'EMP axis'. Moreover, EMP can influence other axes of plasticity such as stemness/tumour initiation potential (Celià-Terrassa et al. 2012;Jolly et al. 2015), resistance to cell death caused by anchorage independence (anoikis) (Huang et al. 2013), resistance to various targeted therapies and immunotherapy across cancers (Chouaib et al. 2014;Tripathi et al. 2016;Dongre et al. 2017;Sahoo et al. 2021a;Shafran et al. 2021), increasing cancer cell fitness during the metastatic cascade. Thus, EMP is a canonical example of phenotypic plasticity and consequent non-genetic heterogeneity implicated in successful metastasis.
Metastasis is a highly inefficient process (Luzzi et al. 1998;Cameron et al. 2000) during which the local microenvironment of disseminated cells is quite dynamic, and cells need to adapt rapidly to survive the bottlenecks they face. The timescale of obtaining the 'right' mutation that can tunnel cells through that bottleneck is over multiple cell divisions, thus being largely inconsequential to the probability of cells surviving that bottleneck. Moreover, while a newly acquired mutation may enhance the survival likelihood of a circulating tumour cell for a given bottleneck, it can compromise the ability of the cell to adapt to additional bottlenecks due to irreversible changes in its phenotypic repertoire. Put together, fast and reversible adaptations at a phenotypic level (i.e., Shachaf and Felsher 2005, non-genetic) seem to be playing an instrumental role in enabling cancer metastasis as compared to slow and irreversible adaptations available at a genetic level. Hence, it is not surprising that while cells from various sub-clones of the primary tumour have been seen in circulating tumour cells and capable of metastasizing (Lyberopoulou et al. 2015;Simeonov et al. 2021), no unique mutational footprints have yet been deciphered, unlike other hallmarks of cancer, for which mutations in various oncogenes and/or tumour suppressor genes have been pinpointed (Hanahan and Weinberg 2011; Mantovani et al. 2019). Such plasticity during metastasis can lead to phenotypic (non-genetic) heterogeneity, as witnessed in circulating tumour cells (CTCs) across cancer types (Yu et al. 2013;Bocci et al. 2021). Consistent with the observed impact of plasticity on the evolvability of cellular traits (Fierst 2011), non-genetic heterogeneity has been identified to impact evolutionary dynamics in lung cancer beyond genetic heterogeneity as well (Sharma et al. 2019).
Besides metastasis, phenotypic plasticity and nonmutational heterogeneity have been implicated in the emergence of adaptive drug resistance (Boumahdi and de Sauvage 2020; Qin et al. 2020;Oren et al. 2021), particularly through drug-tolerant persister (DTPs) -a subpopulation of cells that can survive sustained therapeutic attack by entering a reversible slow-proliferation state (figure 1). DTPs adapt to environmental fluctuations through epigenomic, transcriptional and metabolic reprogramming events, and are capable of expanding into a colony (Shen et al. 2020b). Initially reported in lung cancer (Sharma et al. 2010 Oren et al. 2021;Rehman et al. 2021). Intriguingly, DTPs can act as a reservoir subpopulation through which genetically mutated cells can emerge to stabilize diverse drug resistance mechanisms at a longer timescale (Ramirez et al. 2016). Thus, genetic and non-genetic mechanisms can be thought to cooperate to allow cancer cell adaptation at different timescales during the emergence of drug 'resistance' (Salgia and Kulkarni 2018;Hayford et al. 2021).

Are mutations sufficient for cancer progression?
Investigations into mechanisms of cancer initiation have also questioned the sufficiency of genetic mutations in cancer cells. For example, transformation of a normal cell to a cancerous melanoma cell has been shown to be triggered by imbalance in physiological factors along with exposure to the environmental carcinogen ultraviolet B (UVB) (Berking et al. 2004). Using human skin grafting experiments in immune-compromised mice, it was found that exposing normal melanocytes to increased levels of fibroblast growth factor, stem cell factor and endothelin-3, along with exposure to UVB, could transform normal melanocytes to a cancerous melanoma within four weeks of treatment, while treatment with individual growth factor along with UVB had no effect (Berking et al. 2004). This study suggests that only an external carcinogen is not always sufficient to initiate cancer; instead, some internal physiological imbalance is crucial as a permissive key to trigger neoplastic transformation. Consistently, another study in UV-induced melanoma argues that the susceptibility or resistance of mice to develop cancer strongly depends upon the presence of variants in the modifier genes along with the pathogenic genetic mutation (Ferguson et al. 2019). Thus, apart from genetic mutation(s), the overall genetic make-up of an organism and/or perturbation in the local environment may contribute to the induction of carcinogenesis, offering possible reconciliation between SMT and TOFT.
A commonly asked question in cancer biology is, 'If mutations in cancer-associated genes are sufficient for neoplastic transformation, why do normal cells carrying similar somatic mutations not get transformed and develop cancer?' With advancements in DNA-sequencing technologies, it is now possible to detect lowfrequency somatic mutations in normal cells (Kennedy et al. 2019). The existence of somatic mosaicism (genetically distinct somatic cells harboured by an individual through DNA structural abnormalities, epigenetic changes and errors in chromosome partitioning) is well established (Youssoufian and Pyeritz 2002). Thus, genetic instability is not necessarily a unique property of cancer cells but inherent to all somatic cells, further emphasizing the role of altered microenvironments and/or other permissive cues in tumour initiation (Lichtenstein 2018). For example, Figure 1. Non-genetic heterogeneity in a cell population and its impact on therapeutic efficacy. Upon exposure to a given therapy, a majority of cancer cells die (fractional killing), but a few of them can survive either due to pre-existing mutations and/or pre-existing (non-genetic) heterogeneity or due to additional non-genetic adaptations such as a phenotypic switch. Application of another therapy may lead to further phenotypic diversification. deleterious, age-associated, somatic mutations in the TP53 gene, commonly associated with serous ovarian cancer, are also detected at a very low frequency in the peritoneal fluid from women without cancer (Krimmel et al. 2016). Similarly, cancer-associated somatic mutations arising due to clonal expansion of hematopoietic cells have been reported in healthy individuals (Genovese et al. 2014;Xie et al. 2014). Such mutations are also detected in solid tissues of healthy individuals such as skin, colon, endometrium, brain, etc. (Kennedy et al. 2019). Conversely, in ependymomas -common childhood brain tumoursno significant recurrent mutations were detected in a cohort (MacK et al. 2014).
Together, these studies provide a strong indication that the presence of genetic mutations may not always result in tumour initiation, and that other permissive cues within a cell and/or its local microenvironment may drive neoplastic transformation of normal cells and tumour progression as well.

Evidence of non-genetic heterogeneity in clonal population of cells
The evidence of non-genetic heterogeneity in clonal populations was first noted in microorganisms. For example, Dictyostelium discoideum shows a bimodal (two-peak) distribution of motility speed and calcium content within 15 minutes upon starvation. The observed differences also reverted within 15 minutes after restoration of the nutrient medium, indicating reversible state transitions (Goury-Sistla et al. 2012). Similarly, variations in the levels of intracellular cyclic AMP (cAMP) in Saccharomyces cerevisiae has been shown to be associated with heterogeneity in growth rate and in acute stress tolerance. Perturbing this heterogeneity can impact these functional traits: while increase in intracellular cAMP levels increases susceptibility to acute heat stress, PKA inhibition decreases it, suggesting that underlying population structures get altered in opposite directions by these perturbations (Li et al. 2018).
While cell-extrinsic perturbations can reveal nongenetic heterogeneity, it may also arise due to cellintrinsic variations in functioning of homeostasis in cell organelles. For example, variability in mitochondrial membrane potential among cells was correlated with that in proliferation rate, stress tolerance and resistance to therapy in yeast, identified using high-throughput automated microscopy (Dhar et al. 2019). Similar to unicellular organisms (Elowitz et al. 2002), eukaryotes also exhibit non-genetic cell-to-cell variability. For instance, the clonal population of mouse hematopoeitic cells showed a distribution of levels of the stem cell marker Sca-1. While the Sca-1 high subpopulation with increased expression of PU.1 preferentially differentiated to a myeloid lineage, the Sca-1 low subpopulation expressing higher levels of GATA1 showed greater preference towards erythroid differentiation, highlighting how stochasticity can impact lineage choice in mammalian progenitor cells (Chang et al. 2008).
More recently, a population of cancer cells with identical genetic background has been observed to show differential gene expression and protein levels and consequently functional readouts such as response to drugs and tumour initiation. . Single-cell expression variability is not unique to cancer cells; it has been witnessed among homogenous population of non-cancerous cells too, with important functional implications. For instance, highly variable genes in lung airway epithelial cells were enriched with collagen formation; those in dermal fibroblasts were found to be involved with keratinization, and those in lymphoblastoid cells were enriched with cytokine signaling (Osorio et al. 2020).
Multilineage differentiation programs operated in solid tissues have been proposed as potentially responsible for non-genetic heterogeneity observed in cancer cells. For example, six molecular subtypes of normal Fallopian tube epithelium (FTE; cells of origin of serous ovarian cancer (SOC)) were identified using transcriptomic analysis of 4000 normal FTEs, which was used to deconvolute non-genetic heterogeneity observed in high-grade SOC (Hu et al. 2020). Similarly, using single-cell PCR gene-expression profiling, non-genetic transcriptional variability observed in human colon cancer was demonstrated to be similar to different lineages of normal colon epithelium (Dalerba et al. 2011). Different single-cell transcriptomics or proteomics methods are offering unprecedented insights into elucidating patterns of heterogeneity in a homogenous cell population. For example, co-sequencing of microRNA-mRNA in individual cells using the half-cell genomics approach showed that variability of microRNA levels may drive non-genetic heterogeneity among cells . Similarly, two distinct cell populations within a melanoma tumour were observed to be characterized by variable expression levels of microphthalmia-associated transcription factor (MITF) (Tirosh et al. 2016;Rebecca and Herlyn 2020). Further, subpopulations with varying differential EphA cluster morphologies and intrinsic migration potential were observed in breast cancer cells using single-cell assays (Ravasio et al. 2020). Identification of such heterogeneity has revealed that multiple cancer subtypes may coexist within an individual tumour (Yeo and Guan 2017). Relative proportions of cells exhibiting distinct subtypes constituting a tumour are expected to be highly dynamic and under constant drug-induced evolutionary pressures.
An outstanding example of phenotypic plasticity in cancer is EMP, which includes transition of cells among epithelial (E), mesenchymal (M) and hybrid E/M states (Celià-Terrassa and Jolly 2020). Epithelialmesenchymal transition (EMT) and its reverse mesenchymal-epithelial transition (MET) -which together constitute EMP -are fundamental processes in development and wound-healing where they facilitate the movement of cells from one location to another (Nieto et al. 2016). This property of EMT-MET is exploited by and benefits cancer cells, where it not only confers cell motility (Pearson 2019)  At a cell morphological level, EGF-induced EMT in breast cancer cells can be classified into three distinct reversible morphological states and function in a dosedependent manner: cobble, spindle and circular (Devaraj and Bose 2019). Similarly, using a Z-cad dual-sensor system with an epithelial and a mesenchymal marker together, dynamic changes in breast cancer cells undergoing EMT or MET can be observed, which can help isolate the subpopulation displaying mesenchymal properties from a population consisting of predominantly epithelial-like cells (Toneff et al. 2016). Moreover, the percentage of cells in E, hybrid E/M and M states at various timepoints during EMT induction can be quantified using such a sensor and/or single-cell RNAseq data, highlighting patterns of heterogeneity (Jia et al. 2019;Cook and Vanderhyden 2020;Deshmukh et al. 2021). Reversible changes in the frequency of epithelial and mesenchymal cell states have been seen not only in vitro but also in the circulating tumour cell (CTC) composition of cancer patients with each cycle of response to therapy (Yu et al. 2013). Further, there may be heterogeneity in hybrid E/M states as well, as identified in primary skin and mammary tumours (Pastushenko et al. 2018)  A major reason enabling non-genetic heterogeneity in EMP has been multistability in underlying regulatory networks (Steinway et al. 2015;Font-Clos et al. 2018;Watanabe et al. 2019), which can often introduce asymmetry in the 'paths' taken by cells during EMT vs. those taken during MET. Indeed, transcriptional profiling of metastatic prostate cancer reveals that the expression profiles of cells at various timepoints are not just the reverse of those seen when cells underwent EMT (Stylianou et al. 2019). Similar patterns of hysteresis were seen in proteomics of lung cancer cells (Karacosta et al. 2019). Therefore, despite much investigation, we do not have a unique EMP signature that can be applied in a pan-cancer manner to identify whether cells are undergoing EMT or MET at a given timepoint and, by extension, a signature that can estimate the metastatic potential of cells.
Non-genetic heterogeneity has also been reported for a trait connected with EMT-Cancer Stem Cells (CSCs). CSCs were earlier thought to occupy the apex of cell differentiation (i.e., hierarchical model), but recent evidence has shown that CSCs and non-CSCs can interconvert among one another at both molecular and functional levels (Tang 2012;Thankamony et al. 2020 Another axis connected to EMT that has growing evidence of phenotypic plasticity and heterogeneity is that of drug resistance. Of course, resistance can emerge due to de novo existing mutations in a subpopulation of cells and/or additional mutations gained during therapeutic assault, but non-genetic factors can play a role in driving drug resistance too (Salgia and Kulkarni 2018;Marine et al. 2020;Rebecca and Herlyn 2020). For instance, a partial or full EMT can drive ER? breast cancer cells into resistance to tamoxifen and/or docetaxel; intriguingly, tamoxifen-resistant cells tend to be more mesenchymal than their sensitive counterparts, indicating a mutual causal connection between these axes (Prieto-Vila et al. 2019; Sahoo et al. 2021a). Similarly, non-genetic heterogeneity may provide a mechanism to adapt to drug treatment. By using quantitative proteomics and computational modeling, it was shown that immediately after exposure to RAF/MEK inhibitors such as vemurafenib, the BRAF V600E melanoma 'persister' cells show adaptive changes involving brief pulsatile reactivation of the MAPK pathway which can activate ERK signaling in neighbouring cells too. These pulses enable 'persister' cells to escape cell cycle arrest and sustain long-term resistance at a non-genetic level. This study provides mechanistic detail of the role of nongenetic heterogeneity in emergence of drug resistance in a genetically identical population (Gerosa et al. 2020). Additional analysis of drug-tolerant 'persisters' in melanoma has indicated how vemurafenib treatment can trigger cell-state transitions into a more undifferentiated phenotype which is therapeutically resilient (Su et al. 2017(Su et al. , 2019Pillai and Jolly 2021). Such transitions are often reversible, as seen for EMT (Tripathi et al. 2020), thus enabling resumption of growth upon drug removal. Also, these drug-tolerant 'persisters' can serve as a reservoir of cells some of which may acquire additional mutations at prolonged timescales, leading to 'stabilization' of the drug-resistant phenotype, as seen in lung cancer cells treated with gefitinib (Ramirez et al. 2016).
The concept of persisters was initially reported in bacterial populations which, when exposed to antibiotic drugs, show differential killing rates such that the majority of bacterial cells ( Reinforcing observations are reported in luminal breast cancer using single-cell RNA sequencing, where a subpopulation of 'pre-adapted' cells, with reduced levels of estrogen receptor (ERa) and increased properties of quiescence and migration, undergo transcriptional reprogramming upon drug treatment and gather copy number changes to gain long-term resistance to endocrine therapy (Hong et al. 2019). Another study investigating the mechanism of drug resistance in triple-negative breast cancer has shown that after treatment with doxorubicin combined with cyclophosphamide, the residual tumour cells maintained the sub-clonal architecture of an untreated tumour; however, their transcriptomic, proteomic and histological profiles were different from that of the untreated tumour profiles. Once the drug treatment was stopped, residual tumours gave rise to drug-sensitive tumours with similar expression profiles as that of the untreated tumour, indicating reversible chemotherapy tolerance (Echeverria et al. 2019). Similar analysis of trajectories of escapees from ALK inhibitor treatment was seen in NSCLC, where both genetic and non-genetic mechanisms contributed to gradual adaptation and gained resistance (Vander Velde et al. 2020).
Considered together, emerging evidence along multiple axes of cellular plasticity -EMT, CSCs and drug resistance -strongly endorses the role that non-genetic heterogeneity and consequent possible cell adaptation trajectories can play in facilitating tumour survival, therapy resistance and relapse.

Sources of non-genetic heterogeneity
What are the mechanisms that can result in non-genetic heterogeneity within a clonal population of cells as observed in microorganisms and cancer cells, as well as during embryonic development of a metazoan? Here, we briefly review some canonical sources of non-genetic heterogeneity in different scenarios, and their implications. Transcriptional processes are inherently stochastic, i.e., the production of mRNA and turnover rate of mRNA and protein can produce gene expression noise, which arises due to random binding of transcription factors to the gene promoter (Mcadams and Arkin 1997;Raser and O'Shea 2004). Non-continuous transcription may fluctuate the promoter between an 'ON' state (which results in a transcription burst, producing mRNA transcripts at a high rate) and an 'OFF' state (which pauses transcription) (Raj et al. 2006;Singh et al. 2012;Kumar et al. 2015;Friedrich et al. 2019). The amount of RNA produced from a particular gene depends upon the frequency, amplitude and duration of the transcription burst which may be specific for a particular cell depending upon extrinsic noise such as concentration and availability of general transcription factors (Elowitz et al.  (Sigal et al. 2006).

Epigenetic regulation and transcriptional noise
Transcriptional noise has been implicated as a source for cell-fate decision-making in yeast (Blake et al. 2006), bacteria (Süel et al. 2006) and many mammalian systems (Moris et al. 2016). More recently, it has been shown to influence cancer progression. In leukemia, depletion of acetyl-transferase KAT2A enhances transcriptional bursting and variability, depletes leukemic stem-like cells and delays disease progression (Domingues et al. 2020).

Conformational noise
Intrinsically disordered proteins (IDPs) are proteins, which, unlike many other proteins, lack a well-defined 3D structure either locally or throughout the protein and display a high degree of flexibility which can confer 'conformational noise' due to promiscuity in their interactions (Mahmoudabadi et al. 2013). This 'conformational noise' may manifest differently in different cells, depending on the fluctuations that such IDPs may go through in individual cells. Due to their high conformational flexibility, IDPs serve as a determinant of protein activity and can often be present as hub proteins in biochemical networks (Haynes et al. 2006;Patil et al. 2010). IDPs display faster binding/unbinding rates with their ligands and may undergo transition from disordered to ordered states upon binding with their partner or posttranslational modification, and thus may amplify promiscuity and bring stochasticity in interactions in biochemical networks (Chakrabortee et al. 2010;Bah et al. 2015;Jolly et al. 2018a;Lin et al. 2018). This promiscuity may lead to dynamic rewiring of proteinprotein interaction networks, thus potentially impinging on transcriptional and translational noise also. For example, many drivers of EMT such as ZEB1 and OVOL1/2  have been predicted to contain intrinsically disordered regions (Mooney et al. 2016), adding to the long list of oncogenes and/or tumour suppressor genes where IDP regions have been reported consistently (Lin et al. 2019). Depending on which proteins ZEB1 and/or OVOL1/2 stably interact with in a given time interval, the phenotypic outcome of cells in a population may be varied.
For example, an IDP associated with prostate cancer -Prostate Associated Gene 4 (PAGE4) -has been implicated in phenotypic heterogeneity (Kulkarni et al. 2017;Singh et al. 2021). PAGE4 is phosphorylated by two kinases, namely, Homeodomain Interacting Protein Kinase 1 (HIPK1) and CDC-Like Kinase 2 (CLK2). PAGE4 activates the Activator Protein-1 (AP-1). HIPK1-PAGE4 has a stronger affinity to AP-1 when compared with CLK2-PAGE4. Experimental studies and mathematical modeling have shown that as a result of differential phosphorylation of PAGE4 which leads to an altered AP-1/androgen receptor (AR) regulatory circuit, prostate cancer cells can have a spectrum of phenotypes with varying sensitivity to the standard-of-care androgen deprivation therapy (ADT). Furthermore, the fact that a majority of transcription factors are intrinsically disordered (Staby et al. 2017;Niklas et al. 2018;Tsafou et al. 2018;Zhang and Tjian 2018) lends further credence to the argument.

Tumour microenvironment
The tumour microenvironment shows a high degree of heterogeneity in terms of angiogenesis which modulates oxygen and nutrient availability, composition of stromal and immune cells, endothelial cell density and extracellular matrix composition (Quail and Joyce 2013;Saxena and Jolly 2019). Hypoxia, for instance, confers certain phenotypes on tumour cells present in hypoxic regions, such as stemness (Louie et al. 2010), EMT (Liu et al. 2017) and chemoresistance (Chen et al. 2015). Similarly, varying spatial localization of stromal cells and immune cells with respect to cancer cells can govern the autocrine/paracrine impact they can have on each other. For example, cancer-associated fibroblasts (CAFs) have been shown to expand the breast cancer stem cell population by secreting prostaglandin (PGE2) and IL-6 (Rudnick et al. 2011), to promote migration and invasion mediated by high expression of COX-2 in nasopharyngeal carcinoma (Zhu et al. 2020) and to facilitate neo-angiogenesis in hepatocellular cancer via placental growth factor (Liu et al. 2019). Further, not all fibroblasts are pro-tumour; a subset of them, as identified via single-cell analysis, can also support antitumour immunity (Hutton et al. 2021). Similarly, neutrophils may exert pro-or anti-tumour effects and these subpopulations may also interconvert among one another (Furumaya et al. 2020), reminiscent of macrophage phenotypic heterogeneity exhibited along the M1-M2 axis. Also, feedback loops formed between immune cells and/or cancer cells of varying phenotypes can allow for dynamic phenotypic composition in a tumour, thus influencing its prognosis (Li et al. 2019). Thus, both cell-autonomous (transcriptional, conformational noise and epigenetic changes) and non-cellautonomous (tumour-stroma interactions) can amplify non-genetic heterogeneity in cancer.

Mathematical models to understand non-genetic heterogeneity
In 1957, Waddington proposed an epigenetic landscape model to explain how a pluripotent stem cell can differentiate into multiple lineages represented as valleys (Waddington 1957) (figure 2A). From a dynamical systems theory perspective, these valleys represent 'attractors' in a high-dimensional landscape. These 'attractors' correspond to stable gene expression patterns defining a phenotype, and for a given gene regulatory network (GRN), these 'attractors' can be identified by simulating their emergent non-linear dynamics (Wang et al. 2011;Ferrell 2012;Jia et al. 2017). Many GRNs driving phenotypic plasticity are multistable in nature, i.e., they have multiple attractors.
Thus, cells having these GRNs can acquire more than one phenotype, and can also switch among them under the influence of biological perturbations or noise (Wooten and Quaranta 2017;Li and Balazsi 2018;Sahoo et al. 2020) ( figure 2B). Considering the large number of genes in eukaryotes and their web of complex interactions, there can be many possible 'attractors', but the conditions imposed by the underlying network topology can restrain the 'solution space' to enable acquisition of only a few attractors. These conditions are akin to those imposed by energy minimization principles during protein folding, such that only a limited number of protein configurations are achieved, starting from a given amino acid sequence. Thus, cells can be postulated to traverse in a landscape where each 'attractor' corresponds to stable phenotypes and has a specific basin of attraction (Agozzino et al.  ). Further, the presence of different subpopulations in multiple cancers is well established. For instance, using a Markov model of stochastic cell transition, it was shown that a subpopulation of cells returns to equilibrium phenotypic proportions with time, and thus breast cancer stem cells arise from non-cancer stem-like cells de novo, highlighting the role of stochasticity in enabling phenotypic heterogeneity in a clonal cell population (Gupta et al. 2011). Another study in small-cell lung cancer (SCLC) used Boolean modeling of the underlying GRN to show that the network dynamics can stabilize neuro-endocrine/epithelial (NE) or non-neuroendocrine/mesenchymal-like (ML) phenotypes which act as attractors. Additionally, they also found that when NE and ML cells were treated with cytotoxic drugs, these cells converged towards a hybrid state displaying surface markers of both NE and ML, possibly as a strategy to escape the cytotoxic effects of the treatment (Udyavar et al. 2017). Similarly, in melanoma, identification and simulation of an underlying GRN enabled four different attractors which mapped to the four distinct phenotypes reported experimentally -proliferative, neural crest-like, intermediate/transitory and undifferentiated (Rambow et al. 2018;Pillai and Jolly 2021). This computational analysis could also recapitulate the cell-state transition trajectory observed experimentally upon treatment with vemurafenib through single-cell analysis (Su et al. 2019), offering a platform to identify novel perturbations that can enrich or deplete certain phenotypes.
Mathematical models have helped construct the landscapes of cell-state transitions associated with nongenetic heterogeneity in cancer, such as those for EMT, CSCs, metabolic reprogramming and drug resistance, using both deterministic and stochastic approaches (such as Gillespie simulations; Gillespie 2007) (Font-Clos et al. 2018;Kang et al. 2019;Sarkar et al. 2019;Lang et al. 2021;Sahoo et al. 2021a). Particularly, in EMT, the concept of partial EMT, also referred as hybrid E/M phenotypes, has been largely championed by mathematical models decoding the emergent dynamics of highly inter-connected mutually inhibitory feedback loops involving miR-200/ZEB and miR-34/SNAIL (Lu et al. 2013;Tian et al. 2013). These models predicted that, contrary to previous assumptions, hybrid E/M states are not mere intermediates or 'metastable' states, but are, instead, stable phenotypes that cells can acquire. Further, such mechanistic mathematical models have made experimentally testable predictions about factors stabilizing hybrid E/M states. Validating those predictions in vitro led to identification of 'phenotypic stability factors' (PSFs) -GRHL2, OVOL1/2, NRF2, NUMB, NFATc, among others (Hong et al. 2015;Biswas et al. 2019;Bocci et al. 2019b;Pastushenko and Blanpain 2019;Subbalakshmi et al. 2020). Mathematical models have also elucidated the cell-state transition dynamics upon EMT induction, identifying multiple 'micro-states' and/or hybrid E/M phenotypes that cells acquire en route to EMT in a dose-and/or time-dependent manner (Zhang et al. 2014;Steinway et al. 2015;Font-Clos et al. 2018;Celià-Terrassa et al. 2018;Sha et al. 2020;Deshmukh et al. 2021). Predictions made by mathematical models for coupled EMT-stemness networks (Jolly et al. 2014) about the high tumour-initiating potential of hybrid E/M phenotypes have also been recently validated in vitro and in vivo (Bierie et al. 2017;Pastushenko et al. 2018;Kröger et al. 2019). For instance, the presence of hybrid E/M cells was associated with the worst survival in breast cancer patients and enriched for stem-like cells in different types of breast cancer cell lines with properties such as increased mammosphere formation and higher ALDH1 levels (Grosse-Wilde et al. 2015). Another insight gained by mathematical models of EMT has been that various positive feedback loops in a GRN drive plasticity among epithelial, mesenchymal or hybrid E/M phenotypes (Hari et al. 2020). Indeed, breast cancer cells with the miR-200/ZEB positive feedback loop perturbed via CRISPR had reduced metastatic potential in vivo (Celià-Terrassa et al. 2018), suggesting that mathematical models can not only elucidate the dynamical principles of non-genetic heterogeneity and cell-fate transitions in cancer, but also pinpoint specific therapeutic vulnerabilities to be tested.
The functional role of feedback loops and gene expression noise in enabling drug resistance was recently investigated using synthetic gene network circuits to deconvolute noise from the mean expression of the puromycin-resistance gene with inducible positive and negative feedback loops in Chinese hamster ovary cells. This study demonstrated that the greater noise emerging from the positive feedback loop increased drug resistance at higher concentrations of puromycin, but at lower drug concentrations, it delayed long-term adaptation. Further, a positive correlation between low noise as a result of negative feedback circuits and mutational adaptation driving stable drug resistance was observed (Farquhar et al. 2019). The cross-talk of various positive and negative feedback loops in a cell can therefore influence its sensitivity to various chemotherapeutic assaults, leading to fractional killing (Spencer et al. 2009;Paek et al. 2016;Miura et al. 2018;Guinn et al. 2020). Such non-genetic heterogeneity in a clonal cell population, upon the influence of drug-induced reprogramming, can lead to a rare and stably resistant subpopulation of cells (Shaffer et al. 2017). Population dynamics models capturing such reversible (non-genetic) and/or stable (genetic) resistance scenarios can suggest combinatorial and/or sequential therapeutic strategies to prevent or delay the emergence of tumour (re)growth (Gunnarssson et al. 2020;Cassidy et al. 2021;Sahoo et al. 2021a).

Conclusion
Non-genetic heterogeneity can confer fitness advantage to a clonal population of cancer cells during metastasis, acquisition of therapy resistance and tumour progression. Such heterogeneity can arise from various sources of biological noise within a cell as well as due to multistable dynamics of various underlying networks. Integrated and iterative mathematical-experimental approaches have been instrumental in identifying the sources and implications of non-genetic heterogeneity in cancer. Developing therapeutic strategies which can target the sources of such heterogeneity in isogenic cancer cells may result in higher efficacy in preventing metastasis and tumour progression. Attractors: attractors usually refer to discrete stable steady states of a dynamical system. In the context of cellular decision-making, attractors can be considered as 'valleys' seen in the metaphorical Waddington landscape. They denote the different 'phenotypes', or states, that the cell population can attain during differentiation. Switching among attractors is often called trans-differentiation or de-differentiation, depending on the attractors and trajectories involved in such switching. Noise: noise refers to stochastic effects in various cellular processes such as transcription and translation, which can drive cell-to-cell (phenotypic) heterogeneity. In cells, noise is usually categorized into intrinsic (noise due to stochastic effects that influence the mRNA and protein levels of a given gene in a cell) and extrinsic (noise due to factors that can influence the mRNA and protein levels of multiple genes simultaneously, including extracellular factors such as the microenvironment). Variability: differences between cells, individual organisms, or groups of organisms of a given species caused either by genetic differences (genotypic variability) or by the effect of environmental factors on expression of genetic players (phenotypic variability). It is often quantified using the Fano factor, defined as the ratio of variance to the mean.  Cycling hypoxia induces chemoresistance through the