{'fraction_non_alphanumeric': 0.036237897648686034, 'fraction_numerical': 0.019271553711387737, 'mean_word_length': 5.266214069045212, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Systems-level network modeling deciphers the master regulators of phenotypic plasticity and heterogeneity in melanoma

Summary Phenotypic (i.e. non-genetic) heterogeneity in melanoma drives dedifferentiation, recalcitrance to targeted therapy and immunotherapy, and consequent tumor relapse and metastasis. Various markers or regulators associated with distinct phenotypes in melanoma have been identified, but, how does a network of interactions among these regulators give rise to multiple “attractor” states and phenotypic switching remains elusive. Here, we inferred a network of transcription factors (TFs) that act as master regulators for gene signatures of diverse cell-states in melanoma. Dynamical simulations of this network predicted how this network can settle into different “attractors” (TF expression patterns), suggesting that TF network dynamics drives the emergence of phenotypic heterogeneity. These simulations can recapitulate major phenotypes observed in melanoma and explain de-differentiation trajectory observed upon BRAF inhibition. Our systems-level modeling framework offers a platform to understand trajectories of phenotypic transitions in the landscape of a regulatory TF network and identify novel therapeutic strategies targeting melanoma plasticity.


INTRODUCTION
Melanoma is a highly aggressive cancer arising from melanocytes, the pigment-producing cells of the body. Metastatic melanoma is the deadliest skin malignancy with an abysmal five-year survival rate of 22.5% (Rebecca and Herlyn, 2020). Over 50% of melanomas were found to be driven by BRAF V600E mutation, leading to promising developments in targeted therapy (i.e., BRAF inhibitors -vemurafenib, dabrafenib) over the past decade. Although specificity in genetic mutation and the targeted nature of the treatment provide short term resolution, cases of long-term relapse and recurrence remain prevalent (Savoia et al., 2020). The involvement of non-genetic factors in conferring these abilities to tumors was identified based on drug-resistant subpopulations present in genetically identical cells (Konieczkowski et al., 2014;Hartman et al., 2020). Most of the patients rapidly acquire resistance to these drugs and show signs of relapse within 6-8 months of treatment (Sun et al., 2014;Rebecca and Herlyn, 2020) which limits the cure for metastatic melanoma. Thus, understanding how these non-genetic factors help cells overcome targeted treatment strategies is crucial to develop better therapeutic outcomes.
A major hurdle in treating melanoma is the high degree of phenotypic plasticity and heterogeneity. Melanoma comprises cell subpopulations with diverse phenotypes -'proliferative' and 'invasive', reminiscent of heterogeneity along the spectrum of epithelial-mesenchymal transition (EMT) seen in carcinomas (Jolly and Celia-Terrassa, 2019). These subpopulations in melanoma show varying drug sensitivities and can interconvert among one another in vitro and in vivo during metastasis (Hoek et al., 2008;Verfaillie et al., 2015), thus giving rise to 'moving targets' that is difficult for a 'static' therapy to eliminate. Such reversible cell-state transitions indicate that cellular reprogramming in melanoma is largely driven at transcriptomic and/or epigenetic levels, instead of accumulation of specific DNA mutations. Although individual molecules associated with these phenotypes have been identified, a systems-level investigation of the nonlinear dynamics of cell-state transitions in melanoma remains largely uncharacterized.
MITF, master regulator of melanocytes, is among the most well-studied molecules in melanoma. MITF high cells mark a proliferative phenotype and express melanocytic genes such as MLANA and TYR along with Ki67, whereas MITF low cells denote an invasive one and express WNT5A and DKK1. Tumors formed by either proliferative or invasive melanoma cell lines contained both MITF high and MITF low cells, indicating bidirectional phenotypic switching between these two (Hoek et al., 2008). An intermediate melanocytic/ invasive state expressing NFATC2, SOX6, and ETV4 has also been recently reported . Further, in response to BRAF inhibitors, some melanoma cells can exhibit a neural crest stem celllike (NCSC) state, expressing markers such as NGFR, and proliferating relatively slowly (Fallahi-Sichani et al., 2017;Su et al., 2017Su et al., , 2019Rambow et al., 2018;Tsoi et al., 2018). Some of these cells following this dedifferentiation trajectory can further progress to an undifferentiated invasive (also called 'mesenchymal-like') state, indicating that cells may acquire an NGFR+ state that lies en route transitions between proliferative and invasive phenotypes, similar to recent reports about hybrid epithelial/mesenchymal phenotypes (Pastushenko et al., 2018) in EMT. Thus, a gradient of MITF activity (''revised MITF rheostat model'') can define phenotypic heterogeneity in melanoma: the above mentioned four phenotypes (melanocytic/ proliferative, intermediate/transitory, NCSC, undifferentiated/invasive) and the two relatively less-characterized 'hyper-differentiated' and 'starved' ones ( Figure 1A) (Rambow et al., 2019).
Besides MITF, other reported regulators of phenotypic heterogeneity and/or therapy resistance in melanoma include AXL, c-JUN, BRN2, and PAX3. MITF low /AXL high cells can confer resistance to multiple targeted drugs (Mü ller et al., 2014). Similarly, MITF low /c-JUN high cells in mouse and human tumors are associated with increased myeloid cell infiltration and inflammation-driven dedifferentiation (Riesenberg et al., 2015). The association of MITF with BRN2 and PAX3 is more complicated: overexpression of BRN2 has been seen to both increase and decrease MITF activity in different cell lines, and a reconciliatory explanation has been postulated involving PAX3 that upon BRAF inhibition, MITF is not activated by PAX3 instead by BRN2 which happens to be a weak activator . This PAX3/BRN2 rheostat model-based proposed regulation of MITF is also consistent with the identified role of BRN2 in migration of neural crest cells (Fane et al., 2019). Furthermore, ZEB1 and ZEB2 can regulate MITF and their relative levels can associate with a proliferative-invasive switch, as identified by in vivo lineage tracing (Vandamme et al., 2020). Put together, these associations and/or regulatory interactions offer important mechanistic insights into the roles of these individual molecules in phenotypic plasticity in melanoma. However, these interactions have not been investigated from a dynamical systems perspective to understand how the emergent properties of a regulatory network involving these players can give rise to multiple ''attractor'' states in melanoma to explain the diverse set of observations of phenotypic switching.
Here, we integrate correlative metrics among different genes with dynamic modeling of underlying gene regulatory network to explain the emergence of distinct phenotypes in melanoma that can switch among themselves as noted experimentally. Our model simulations can recapitulate the four phenotypes seen in melanoma, identify transcriptional regulators of distinct cell states and phenotypic heterogeneity in melanoma, as well as explain de-differentiation trajectory observed upon BRAF inhibition. This modeling framework offers a platform to predict the diverse dynamic trajectories of phenotypic switching in melanoma and to identify novel therapeutic strategies.

Two distinct teams of molecules underlying phenotypic heterogeneity
First, we identified a set of genes that regulate and/or correspond to different phenotypes reported in melanoma, based on existing literature. MITF, MLANA, TYR, TRPM1, and other melanocytic genes are upregulated in proliferative cells, whereas AXL and WNT5A are affiliated with an invasive phenotype (Weeraratna et al., 2002;Hoek et al., 2008;Mü ller et al., 2014). Further, EMT-inducing transcription factors ZEB1 and ZEB2 can play opposite roles in melanoma, with ZEB2 enabling a more proliferative state while ZEB1 promoting invasion and metastasis (Vandamme et al., 2020). Other molecules such as PAX3, SOX10, and c-JUN regulate phenotypic plasticity too (Riesenberg et al., 2015;Smith et al., 2019). Thus, we calculated pairwise correlations between every two genes, across multiple datasets of melanoma cell lines and human samples.
The correlation matrices revealed an intriguing trend: the invasive phenotype genes ZEB1, JUN, WNT5A, and AXL formed one team of players which were all positively correlated with one another (upper red triangle in Figures 1B and S1A). Similarly, genes corresponding to proliferative phenotype -PAX3, ZEB2, SOX10, TYR, MLANA, TRPM1, and MITF -were positively correlated among one another, forming another such team (lower red triangle in Figures 1B and S1A). The players of these two teams were found to be negatively correlated with one another (purple rectangle in Figures 1B and S1A). These results suggest that the proliferative and invasive phenotype genes tend to form two 'teams' of mutually opposing players iScience Article that can govern the emergence of phenotypic plasticity and heterogeneity in melanoma. Such 'teams' have been witnessed between EMT-inducing and EMT-inhibiting factors (Jia et al., 2020b), as well as those between master regulators of 'classic' and 'variant' subtypes in small cell lung cancer (Chauhan et al., 2021), elucidating a potential common design principle for regulatory networks involved in cancer cell plasticity.

Presence of well-defined proliferative and invasive phenotypes
To examine whether these two 'teams' correspond to distinct phenotypes in melanoma in different datasets, we used K-means clustering on top 3,000 genes with the highest variance (Tsoi et al., 2018  iScience Article clustering is used to then identify the optimal number of clusters in these datasets, based on AIC and BIC scores calculated for K = 1 to K = 15. The highest peak difference in these scores was noticed from K = 1 to K = 2, indicating that two clusters are optimal ( Figure S1B).
These two clusters segregate well when projected on their first two principal components ( Figure 2A). We next checked the enrichment of reported proliferative and invasive signatures in these two clusters across datasets using gene set enrichment analysis (GSEA) (Hoek et al., 2008;Verfaillie et al., 2015). The two clusters obtained for each dataset show enrichment of either the proliferative or invasive signatures reported earlier ( Figures 2B, S2A and S2B), thus highlighting that these two clusters can be mapped onto proliferative and invasive phenotypes as initially identified.
The presence of these two phenotypic clusters was further confirmed in melanoma samples at bulk (GSE112509) and single-cell (GSE81383) levels ( Figures S2C and S2D). These samples also showed the emergence of two 'teams' of players earlier seen for melanoma cell lines ( Figure S2E). Put together, these results underscore that these two distinct phenotypes: proliferative and invasive can be seen in diverse publicly available melanoma datasets -at bulk and single-cell RNA-seq levels, as well as cell lines and primary/metastatic samples.

Identifying a master regulator network for phenotypic heterogeneity in melanoma
To identify the genes regulating diverse phenotypes in melanoma, we used an unbiased approach (Figure 3A) by first classifying all genes into modules based on correlation among them through WGCNA (Langfelder and Horvath, 2008). These modules comprise co-expressed genes; some of these modules may correspond to specific phenotypes (Udyavar et al., 2017). To calculate the overall expression of each of these modules in a given melanoma tissue sample in GSE4843, we computed corresponding eigengene values (first principal component of all genes within a module). Module eigengene values for salmon (henceforth called invasive module) and yellow module (henceforth called proliferative module) were significantly different across the proliferative and invasive samples (Bonferroni adjusted p value <0.001) ( Figure S3A, i). In proliferative samples, eigengene values of the proliferative module were upregulated and those of the invasive module were downregulated ( Figure 3B, i). Conversely, in invasive samples, eigenvalues of the proliferative module were downregulated and those of the invasive module were upregulated ( Figure 3B, ii). Thus, these two modules can be considered to correspond to two distinct phenotypes which largely exhibit mutually exclusive gene expression patterns ( Figure S3A, ii).
Next, we sought to decipher the master regulators of these two modules of interest, using two orthogonal and unbiased approaches. First, we use EnrichR (Chen et al., 2013) which scans publicly available databases to identify potential regulators of genes comprising our modules of interest. Second, we infer a backbone regulatory network using ARACNE (Margolin et al., 2006) on GSE4843 for a set of transcriptional regulators. Nodes in this network were selected based on Fishers' exact test to identify the master regulators corresponding to each module (Lambert et al., 2018). Finally, we found common regulators identified by these two complementary approaches. The master regulators identified formed two mutually inhibiting groups or 'teams', similar to the interactions seen in known regulators of phenotypic heterogeneity, across different in vitro and tumor sample datasets ( Figures 3C, S3B). Interestingly, upon repeating the same process with two other datasets (CCLE and GSE10916), the candidate master regulators identified did not generate as strong ''teams''-like behavior in the remaining datasets, as compared to regulators identified from GSE4843 ( Figure S3C), thus, justifying our choice for selecting this dataset.
The network obtained via ARACNE is undirected in nature. Thus, we arrived at a master regulatory network for melanoma by collating information on directed links in the network ( Figure 4A) through manual curation of the above-mentioned regulators by using public databases (ChEA, ENCODE) and examining relevant experimental literature (Table S1).

Emergent network dynamics recapitulate phenotypic heterogeneity in melanoma
Next, we simulated the dynamics of the 17-node regulatory network identified ( Figure 4A) using an ensemble of kinetic parameter sets to capture the cell-to-cell variability in gene expression levels and/or biochemical reaction rates, through the RACIPE (Random Circuit Perturbation) method (Huang et al., 2017). The steady state values of expression levels obtained via RACIPE (further referred to as the simulated dataset) for each gene captures the possible expression levels of these genes within a cell.
Hoek Invasive   iScience Article dataset had optimal clustering for K = 2, based on silhouette width estimation for K-means clustering (K = 2 to 10) ( Figure S4A, i). The two clusters were spatially separated along the first principal component (PC1) (Figure 4B, i) and exhibit differential expression of several genes in the network ( Figure 4B, ii). To identify relevant genes involved in this classification, we calculated the contribution of each gene to PC1 ( Figure 4B, iii) and the extent of representation of each gene by PC1, based on their squared cosines ( Figure S4A, ii). Based on both of these metrics, MITF, FOS, KLF4, NFIC, TBX3, NR2F1, SMAD3, ETV5, JUN, TFE3, and NR3C1 were considered to be important for segregation of the two clusters. The distributions for these genes exhibit bimodality ( Figures 4C and 4D), further supporting their differential expression in two clusters. Several poorly contributing genes (AHR, SMAD4, FOXF1, MAFB, STAT5A, and TFAP2A - Figure 4B, iii) do not exhibit significant bimodality based on Hartigan's dip test and Bimodality Coefficient (Hartigan and Hartigan, 1985;Pfister et al., 2013) (Figures 4C and S4B) in the simulated data. Therefore, these genes are less likely to be able to distinguish between the two clusters. Many of the genes important for segregation of above mentioned two clusters (TBX3, JUN, KLF4, NR2F1, and SMAD3) exhibited bimodality in single-cell transcriptomics for a cohort of single cells cultured from biopsies of three different patients (GSE81383; Figures S4C and S4D).
To confirm whether these two clusters correspond to the two phenotypes (proliferative, invasive), we compared the simulated dataset with multiple experimental datasets for the expression of 11 genes relevant to segregation along PC1 (MITF, FOS, KLF4, NFIC, TBX3, NR2F1, SMAD3, ETV5, JUN, TFE3, and NR3C1). Proliferative and invasive samples (as identified earlier in Figure 2A by GSEA of respective signatures) exhibit differential expression of these 11 genes, as predicted by the simulated data ( Figures 5A and  S5A). The samples are also spatially separated along the first two principal components for these 11 genes, across diverse in vitro and tumor sample datasets ( Figure 5B). In addition, correlation plots for these datasets and the simulated dataset show similar patterns, where these 11 genes seem to split into forming two teams of 8 genes and 3 genes ( Figures 5C and S5B), similar to our earlier observations ( Figure 1B). One of these teams includes MITF, and the other includes JUN, thus reminiscent of antagonism between them and correspondingly opposite phenotypes associated with these molecules (Riesenberg et al., 2015). Together, this analysis confirmed that dynamic simulations for the master TF regulatory network we inferred can recapitulate the existence of the proliferative and invasive phenotypes.

Relative proximity and transition paths among four co-existing phenotypes in melanoma
Because the mean silhouette scores for K = 3 and K = 4 suggested their ability to offer efficient clustering as well ( Figure S4A, i), and previous reports pinpoint four possible phenotypes in melanoma (Tsoi et al., 2018), we looked at the possibility of four clusters corresponding to the four phenotypes. In the simulated data, these four clusters form distinct groups along the first two principal components ( Figure 6A, i). Interestingly, the four clusters appear to be further classification of the two clusters previously observed ( Figures  6A, ii and4B, ii), similar to the classification observed in experimental data (Tsoi et al., 2018). To identify genes that distinguish between the proliferative and invasive sub-clusters, we used Linear Discriminant Analysis (LDA). LDA identifies a linear combination of genes that maximizes separability of clusters, thus, the loading scores of each gene were analyzed to identify genes that contribute the most. A cut-off of 0.6 was set to classify a gene as a discriminant. MITF was identified as a discriminant between the two proliferative sub-clusters ( Figure S5C, i) and TBX3, NR2F1, KLF4, and FOXF1 as discriminants for the two invasive sub-clusters ( Figure S5C, ii). To check if these four clusters map onto the four phenotypes observed, we examined the expression of these determinant genes in experimental datasets. To classify a given dataset into these four phenotypes, we used K-means clustering on the basis of previously defined markers of the 4 phenotype -MLANA for the melanocytic phenotype, ETV4 for transitory phenotype, AXL for undifferentiated phenotype and NGFR for NCSC ( Figure S6A). Differential expression of discriminant genes was observed at single cell resolution in melanoma cell lines (GSE134432) as well as in tissue culture samples from tumours (GSE4843) ( Figure 6B). This analysis suggests that the four phenotypes observed in melanoma can exist at an individual cell level as well, thus facilitating heterogeneity at a bulk level.
Interestingly, we observed that the proximity between clusters for the simulated dataset reflected features of the dedifferentiation trajectory followed during drug treatment (Su et al., 2017;Tsoi et al., 2018). Based on the projection of expression data on the PC axes ( Figure 6A, i), the melanocytic cluster is found to be closest to the transitory cluster, and the undifferentiated invasive cluster is closest to the NCSC one, suggesting that during their transition to undifferentiated phenotype, melanocytic cells pass through intermediate phenotypes that are transcriptionally similar. To confirm this prediction about relative positioning of iScience Article the phenotypes in the experimental datasets, we measured the shortest Euclidean distance between the clusters, based on expression levels of master regulators ( Figure S6B). Validating the predictions made by the model, the intermediate phenotypes (transitory and NCSC) are found to be closest to the melanocytic and undifferentiated phenotypes which explains the path followed in the dedifferentiation trajectory (from melanocytic to transitory to NCSC to undifferentiated). Furthering our confidence in cell-switching trajectories, a pseudo-time trajectory analysis for de-differentiation (based on the 17 genes in our network) seen in single-cell RNA-seq data demonstrated that cells pass through the transitory phenotype, then achieve the NCSC phenotype and finally reach the undifferentiated phenotype ( Figure 6C). Thus, the cell-fate landscape formed by the emergent dynamics of the master TF network can not only explain the existence of four phenotypes earlier reported in melanoma, but also identify most likely transition paths or trajectories that cells follow while transitioning among these states.
We next asked how unique is the outcome of the regulatory network inferred here in enabling these 4 states. Thus, we generated 100 random variations ('mutated networks') of the original ('wild-type'/WT) network shown in Figure 4A, by randomly shuffling links among various edges. For each such mutated network, we calculated team strength between invasive and proliferative modules (Chauhan et al., 2021), as well as quantified how different its phenotypic distribution was from that of the WT network. The team strength for these random networks is relatively weak as compared to that seen in the WT network; moreover, the weaker the team strength for a given network, the further the phenotypic distribution is from that of the WT one ( Figure S7A, i-ii). This analysis highlights that the phenotypic heterogeneity in melanoma emerges from dynamics of underlying network inferred here.
We also examined what kinetic parameters contribute the most to multistability and consequently phenotypic plasticity and heterogeneity. One proposed way of reducing multistability of a system is to decrease its non-linearity, which can be done by decreasing the Hill coefficients (Gardner et al., 2000). Upon setting the maximum limit of Hill coefficients to be chosen by RACIPE as 3, as compared to the default limit of 6, we observed a decrease in frequency of multistability seen for the WT network ( Figure S7B). We also used an unbiased method of identifying important parameters that segregate scenarios giving rise to monostability and multistability. We identified the top three parameters that give rise to multi-stability to be inhibitory constants for molecules across the two teams ( Figure S7C, i, FOS in the proliferative team, and NR3C1 and AHR in the invasive team team). This further reinstates the importance of mutual antagonism between the two teams in giving rise to multiple phenotypes. Interestingly, removal of these edges, significantly decreases the number of positive feedback loops having less than 5 ( Figure S7C, ii) and 6 ( Figure S7C, iii) edges, which has been characterized to play an important role in facilitating multistability in large networks (Hari et al., 2020). Overall, the ''teams'' structure is crucial to enable 4 phenotypes observed in our network.

Network-based modeling explains population dynamics during BRAF inhibition
Resistance to MAPK/BRAF inhibitors is accompanied by an increase in invasive cell population (Hoek et al., 2008;Tsoi et al., 2018;Su et al., 2019Su et al., , 2020. We wanted to check if our network could explain such phenotypic switching in the presence of drug treatment. One of the downstream targets of MAPK and BRAF is MITF. Hence, to represent drug treatment, we simulated the network under MITF knockdown (KD) conditions using RACIPE. Interestingly, we observed a shift in the density of steady-state solutions (i.e., phenotypic composition) from proliferative towards invasive phenotypes ( Figure S8A). This shift is caused by a significant decrease in the frequency of cells belonging to melanocytic clusters ( Figure 7A). Although there is no noticeable change in proportion of undifferentiated phenotype, we see a significant increase in the frequencies of the transitory and NCSC cell phenotypes. This change is also accompanied by an increase in the coefficient of variance for the distributions corresponding to the proliferative and invasive phenotypes, suggesting higher phenotypic heterogeneity in response to the drug, and also an increased dominance of these intermediate phenotypes in the overall cell population composition ( Figure 7B). Similar trends have also been observed experimentally, where short term MAPK inhibition leads to a large increase in Starved Melanoma Cell (SMC, a possible form of the transitory phenotype) and NCSCs (Rambow et al., 2019).
To identify whether knocking down any of the other genes in our master TF network can also drive phenotypic switching, we performed in silico knockdowns of each gene. Our results showed that knockdown of none of the other genes was as potent as that of MITF in terms of inducing phenotypic switching, suggesting that this effect is unique to MITF knock down ( Figure S8B). Such a drastic effect observed upon

DISCUSSION
Phenotypic heterogeneity in melanoma has been associated with drug sensitivity, immunotherapy resistance and dedifferentiation of melanoma cells (Reinhardt et al., 2017;Su et al., 2017;Tsoi et al., 2018).
Here, we identify a transcription factor (TF) regulatory network whose emergent dynamics can settle into multiple ''attractors'' states, each corresponding to a specific phenotype. The concept of ''attractors'' is embedded in the mathematical realization of the Waddington's epigenetic landscape, a powerful approach to understand the dynamics of cellular differentiation in development and increasingly in cancer too (Zhou and Huang, 2011;Li et al., 2016;Su et al., 2017). Cells in different ''attractors'' may interconvert among one another driven by biological noise, for instance, stochasticity in gene expression which can result in phenotypic transitions, as seen in various cancer cell populations (Gupta et al., 2011;Shaffer et al., 2017Shaffer et al., , 2020. Thus, a tumor may contain cells belonging to different ''attractors'' in various timedependent and space-dependent proportions (Yeo and Guan, 2017;Wouters et al., 2020). Quantifying corresponding cell-state transitions and their landscapes is thus essential for a better understanding of diverse adaptive strategies cancer cells can take Agozzino et al., 2020).
The four phenotypes seen here for melanoma is a refined classification of two broader states: proliferative and invasive, similar to the classification system used previously (Tsoi et al., 2018). Thus, our model proposes the possibility of two ''macro-states'' corresponding to the proliferative and invasive phenotypes, that comprise ''micro-states'' corresponding to the sub-clusters ( Figure 7C). These microstates can further be constituted of an ensemble of microstates; for instance, the hyper-differentiated cells may be thought of as a microstate for the melanocytic phenotype, and starved melanoma cell (SMC) for the intermediate phenotype (Rambow et al., 2019). Transitions among microstates have lower barriers and are postulated to occur relative easily as compared to transition between macro-states, which might need a stronger ''driving force'' such as drug treatment (i.e. BRAF inhibition) or genomic silencing. A similar model has been proposed in the context of EMT (Goetz et al., 2020), where such ''micro-states'' may facilitate cellular adaptation.
Several parallels can be drawn between phenotypic switching in melanoma and EMT in epithelial cancers. First, the invasive phenotype has upregulation of several mesenchymal genes and EMT-regulators such as ZEB1, WNT5A, CDH2 (Su et al., 2017;Vandamme et al., 2020), as well as that of matrix metalloproteases (MMPs) that can change the extracellular matrix (ECM) stiffness (Arozarena and Wellbrock, 2017), as typically observed in EMT (Wei et al., 2015;Deng et al., 2020). However, ZEB1 and ZEB2, both EMT-inducing factors, promote opposite phenotypes in melanoma (Denecker et al., 2014). In addition, MITF directly regulates ECM and focal adhesion pathways, which have been associated with inducing phenotypic switching (Dilshat et al., 2021;Spoerri et al., 2021). Also, EMT-ECM bidirectional crosstalk was recently implicated in phenotypic switching (Deng et al., 2021), endorsing a melanoma-ECM-EMT axis governing phenotypic heterogeneity. Other EMT regulators such as TWIST and AP-1 (Feldker et al., 2020) are also reported to drive phenotypic plasticity in melanoma (Caramel et al., 2013;Verfaillie et al., 2015). Second, heterogeneity along the EMT axis seen in circulating tumor cells (CTCs) (Yu et al., 2013;Bocci et al., 2021) is recapitulated in melanoma CTCs too (Aya-Bonilla et al., 2020). Third, cooperation among phenotypes with varying EMT status has been witnessed during metastasis and tumor formation (Jolly and Celia-Terrassa, 2019), reminiscent of recent observations in melanoma (Rowling et al., 2020). Further investigations into the dynamics of phenotypic transitions can elucidate whether melanoma cells exhibit hysteresis (Karacosta et al., 2019) and/ or spontaneous/stochastic state switching  that characterize EMT.
Our study characterizes phenotypes based on dynamics of the TF regulatory network. However, previous studies have also identified translational, epigenetic and metabolic features of distinct phenotypes, suggesting multiple interconnected levels of regulation of phenotypic plasticity in melanoma (Bettum et al., 2015;Hugo et al., 2015;Falletta et al., 2017;Lionetti et al., 2020). Variations in chromatin accessibility, DNA and histone modifications can be mapped to transcriptional heterogeneity seen in melanoma (Hugo et al., 2015;Verfaillie et al., 2015). Two regulatory networks corresponding to proliferative and invasive phenotypes have been identified, with SOX10 and MITF acting as key regulators of the proliferative network and TEAD and AP-1 as regulators of the invasive one. iScience Article epigenetic signatures, the proliferative and invasive phenotypes have been characterized to have distinct metabolic profiles too. A switch from proliferative to invasive phenotypes is supported by metabolic reprogramming from oxidative phosphorylation (OXPHOS) to glycolysis mediated by MITF target gene, PGC1a which promotes OXPHOS. Lowering MITF and PGC1a can shift cells toward glycolysis (Vazquez et al., 2013). A recent theoretical model predicted the possibility of four metabolic profiles -OXPHOShigh/glycolysis-low, OXPHOS-low/glycolysis-high, OXPHOS-low/glycolysis-low, and OXPHOS high/ glycolysis-high (Jia et al., 2020a). The first two profiles can be mapped onto the proliferative and invasive phenotypes, based on previous experimental evidence. The ability of intermediate phenotype to capture features of proliferative and invasive cells raises the possibility of its ability to exhibit both OXPHOS and glycolysis. The low/low phenotype been identified to map onto an idling population of cells, with no net growth. NCSCs, considered to be slow-growing and to maintain minimal residual disease (Liguoro et al., 2020), may possibly represent the 'idling' population of melanoma cells (Paudel et al., 2018) exhibiting both low glycolysis and OXPHOS.
The physiological relevance of phenotypic switching is observed in drug resistance in melanoma. The more invasive phenotypes, NCSC and undifferentiated cells, are observed to emerge along differing time scales upon drug exposure. While shorter duration of BRAF and RAF/MEK inhibition gives rise to NCSC phenotype (NGFR+ cells), prolonged exposure facilitates undifferentiated cells (Fallahi-Sichani et al., 2017;Su et al., 2017). How does adaptation at different time scales contribute to emergence and maintenance of drug-tolerant persisters (Shaffer et al., 2017;Ahmed and Haass, 2018;Schuh et al., 2020), and the role of reversible cellular reprogramming (Roesch et al., 2010;Karki et al., 2021) in enabling such persisters requires further careful investigation. Understanding the dynamics and trajectories of phenotypic switching in these scenarios can help guide better combinatorial (Boshuizen et al., 2018;Luo et al., 2018) and/or sequential (Goldman et al., 2015) therapies that can target vulnerabilities of diverse phenotypes. Our simulations were able to recapitulate the phenomenon of phenotypic switching observed upon BRAF inhibition, thus offering our mechanism-based model as a possible platform to identify potent interventions.
Overall, we have identified a master TF regulatory network that explains the (co-) existence of diverse phenotypes observed in melanoma cell lines and tumors, and the phenomenon of phenotypic switching, which is crucial in melanoma development, progression and metastasis. To the best of our knowledge, ours is the first non-linear dynamical analysis of melanoma cell states as well as a meta-analysis across datasets to decode the principles of phenotypic heterogeneity and plasticity. Our model depicting different ''attractor'' states also provides a framework to develop future therapeutic strategies, by identifying the key drivers of distinct drug-resistant sub-populations. For instance, a recent CRISPR-screening based study has already confirmed the role of two molecules -AHR and SMAD3 -in conferring drug resistance to invasive phenotypes (Gautron et al., 2021). Thus, by targeting the master regulators of given phenotype(s) that we can identify through systems-level integrative experimental and computational analysis, as done for small cell lung cancer (Udyavar et al., 2017;Wooten et al., 2019), we can "push'' the cells towards more drug-sensitive phenotypes, improving the outcomes of conventional therapy options in use.

Limitations of the study
Our data-driven analysis identifies a transcriptional regulatory network governing the existence of multiple phenotypes in melanoma. This network includes some well-reported master regulators such as MITF, AHR, and SMAD3, but also includes additional players such as KLF4 and NR2F1 which need further investigation from an experimental standpoint. In addition, other layers of regulation (epigenetic, metabolic, and translational) are not yet integrated into our framework.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

AUTHOR CONTRIBUTIONS
MKJ designed and supervised research. MP performed research and analyzed the data. Both authors contributed to writing the manuscript.

DECLARATION OF INTERESTS
The authors declare no conflict of interest.
iScience Article iScience Article from the cluster means. The algorithm calculates cluster centroids and iteratively repeats the assignment of datapoints to the clusters, till no further changes occur. The kmeans function from the stats package was used to perform k-means clustering on the datasets.

Selection of optimal number of clusters
Akaike Information Criterion (AIC), Bayesian Information criterion (BIC) and Silhouette scores were used to determine optimal number of clusters to classify the data into. The metrics are calculated as follows: AIC = À 2logðlikelihoodÞ + 2K BIC = À 2logðlikelihoodÞ + KlogðNÞ where K is the number of parameters estimated by each model, and N is the sample size for a given model. For K-means clustering, twice-negative log likelihood is estimated as the total within cluster sum of squares. To determine the optimal number of clusters, the model with the largest dip in AIC and BIC score was selected. Plots were generated for the difference in score of the n th and (n-1) th model.

Identification of two phenotypes
a) K-means clustering: It was used to generate two clusters of cells based on the top 3,000 genes with the highest variance.
b) GSEA: Each cluster was checked for enrichment of Hoek proliferative, Hoek invasive (Hoek et al., 2008), Verfaillie proliferative and Verfaillie invasive gene sets (Verfaillie et al., 2015) based on a ranked gene list calculated using the signal-to-noise ratio for all genes in the expression matrix. Gene Set Enrichment Analysis (GSEA) software (Subramanian et al., 2005) was used for creating the ranked gene list and calculating enrichment scores.

PCA
Principal component analysis (PCA) was used to visualize multidimensional gene expression and simulation data. factomineR and factoextra packages were used to calculate and visualize the datapoints on the principal components. To determine the correlation between variables and the representation of variables by the principal components, a correlation circle with squared cosines was plotted.

LDA
Linear discriminant analysis (LDA) is a dimensionality-reduction method that provides linear combinations of variables (discriminant functions) to maximize the separability and differences between multiple classes of data. To identify variables that contribute the most to the separation of two classes of data, we plotted the weights or loading scores of each variable in the first discriminant function. A cut off value of 0.6 was set to select genes that contribute to the separation of the subclusters. The package MASS (Venables and Ripley, 2002) was used to calculate the discriminant functions.

Identification of gene modules for phenotypes
To generate the network underlying the two phenotypes, GSE4843 (Hoek et al., 2006) was subjected to Weighted Gene Correlation Network Analysis (WGCNA) (Langfelder and Horvath, 2008). The topological overlap matrix was constructed based on the Pearson correlation matrix for the gene expression data. A soft threshold power of 4 was identified as the least power that generates a scale free network fit with R 2 > 0.9. To identify relevant co-expression modules, unsupervised average-linkage hierarchical clustering was used and adaptively pruned at height 0.995. Minimum cluster size was set to 100 to ensure a sufficient size of modules. Merging of modules was repeated at dissimilarity threshold of 0.25 to ensure that similar modules were merged (Tremblay et al., 2019) (Table S1).
Modules of interest were identified based on extent of significant differential expression of genes across the two phenotypes based on Bonferroni adjusted p value (<0.001). Module eigengene values and heatmaps for expression levels of all genes in each module were generated.

Identification of master regulator network
To identify the master regulators for the differentially expressed genes obtained from WGCNA, we used geWorkbench ( (Margolin et al., 2006). A p value of 10 À7 was set to determine the mutual information threshold and the software was run for 100 bootstraps with data processing inequality set to 0. Fisher's exact test was used to identify master regulators from a list of candidate master regulators (Lambert et al., 2018). Only those transcription factors (TFs) enriched for in the WGCNA modules with p value < 0.05 were considered for further analysis. This list was cross validated against CHEA (Lachmann et al., 2010), ENCODE  and ARCHS4 (Lachmann et al., 2018) databases by using EnrichR (Chen et al., 2013) to identify potential TFs regulating each module. Only those TFs identified by both analyses (ARACNE and EnrichR) were considered as master regulators (Table S2).
Interactions between the master regulators was established by manual curation of literature and publicly available databases (Table S3). Nodes having no incoming or outgoing edges were removed, to finally arrive at a 17 gene network.

RACIPE
RAndom CIrcuit PErturbation (RACIPE) (Huang et al., 2017) was used to generate an ensemble of ordinary differential equation (ODE) models. Each model represents a collection of modified Hills equations for each gene, with randomized kinetics parameters sampled from user-defined ranges. Each model is solved for multiple initial conditions, till steady state solutions are obtained. An ensemble of steady state solutions for multiple models is used to represent heterogeneity within a population of cells. For our analysis, we used 10,000 ODE models, each with 100 initial conditions. Additional details on parameters have been given in Table S4.

Trajectory analysis
To estimate transcriptional similarity between clusters (i.e. phenotypes) we calculated the Euclidean distance between cluster centers. Euclidean distance (D) between two points is given as the square root of sum of squares of distances and is calculated as follows: D = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X n i = 1 ðx 2 i À x 1 i Þ 2 q where n represents number of genes, x 1i represents value of i th gene for point 1 and x 2i represents value of i th gene for point 2.

Pseudo-time analysis
For identifying the pseudotime trajectory of cells in GSE134432, we used the monocle package. The network genes (17 master regulators) were used for ordering and dimension reduction. The method used for dimension reduction was 'DDRTree' which also allows ordering of samples.

Edge removal
For edge removal analysis after LDA ( Figure S7), the edges that were removed one at a time were: (i) Inhibition of FOS by NR3C1, (ii) Inhibition of FOS by NFIC and (iii) Inhibition of AHR by FOS. For the control case, we considered the original network ( Figure 4A), and two random edges that contribute to teams structure (one activating within a team and the second, inhibiting across teams) -activation of KLF4 by JUN and Inhibition of TFE3 by MITF, removed one at a time.

Modelling BRAF/MAPK inhibition
To model the effect of BRAF/MAPK inhibitors on the four phenotypes, we simulated the network using RA-CIPE with and without MITF knock down (experiment and control cases, respectively), because MITF is a direct target of MAPK and BRAF. To cluster the samples, k-means clustering was used. Both the simulated datasets were projected onto the principal components for the control network and the density distribution of points along PC1 was measured. To determine the Coefficient of Variance (CV), we fit a Gaussian curve to each peak using package mixtools. The curve was used to identify the mean (m) and standard deviation (s) and calculate CV as follows: CV = s/m. To assess the extent of bimodality within genes, we used 2 metric: Hartigan's dip test and the Bimodality coefficient (BC) (Hartigan and Hartigan, 1985;Pfister et al., 2013

QUANTIFICATION AND STATISTICAL ANALYSIS
Statistical analysis was performed using R version 3.6.3, unless mentioned otherwise. Details of statistical analysis, definitions for significance and abbreviations can be found in figure legends.