Cooperative RNA degradation stabilizes intermediate epithelial-mesenchymal states and supports a phenotypic continuum.

Multiple intermediate epithelial-mesenchymal transition (EMT) states reflecting hybrid epithelial and mesenchymal phenotypes were observed in physiological and pathological conditions. Previous theoretical models explaining multiple EMT states rely on regulatory loops involving transcriptional feedback, which produce three or four attractors. This is incompatible with the observed continuum-like EMT spectrum. Here, we used mass-action-based models to describe post-transcriptional regulations, finding that cooperative RNA degradation via multiple microRNA binding sites can generate four-attractor systems without transcriptional feedback. Furthermore, the newly identified intermediates-enabling circuits are common in the EMT regulatory network, and they can synergize with transcriptional feedback to support phenotypic continuum. Finally, our model predicted a role of miR-101 in multistate EMT, and we identified evidence from single-cell RNA-sequencing data that support the prediction. Our work reveals a previously unknown role of cooperative RNA degradation and microRNAs in EMT, providing a framework that can bridge the gap between mechanistic models and single-cell experiments.

Introduction

Epithelial-mesenchymal transition (EMT) is a cell state change required for embryogenesis, postnatal development, and some diseases’ progression including metastasis (Thiery et al., 2009; Watanabe et al., 2014). During EMT, epithelial cells lose their apical-basal polarity and gain the ability to migrate. The transition is not a binary switch: intermediate cellular phenotypes between epithelial (E) and mesenchymal (M) states have been found in development, mammalian cell lines, fibrosis, and tumors (Aiello et al., 2018; Grande et al., 2015; Hong et al., 2015; Pastushenko et al., 2018; Zhang et al., 2014). These intermediate EMT states may possess characteristics of both E and M phenotypes, which can be important for cell- or tissue-level functions such as collective cell migration (Campbell et al., 2019). Recent single-cell transcriptomic studies using individual cell lines provide additional evidence for the existence of multiple intermediate EMT states (Cook and Vanderhyden, 2020, 2022; Deshmukh et al., 2021; Pastushenko et al., 2018). To investigate the mechanistic basis of the intermediate states, theoretical and experimental approaches were used to demonstrate that interconnected transcriptional feedback loops can support intermediate EMT states (Hong et al., 2015; Lu et al., 2013b; Steinway et al., 2015; Tian et al., 2013). These existing mechanistic models can produce only one or two intermediate EMT states with a set of biochemical rate constants representing one cell type. Contrary to these models, it was recently shown that accurate description of EMT dynamics in a single mammary epithelial cell line requires several more intermediate states (Goetz et al., 2020). This finding is consistent with transcriptomic data showing a continuum-like EMT spectrum which possibly contains many stable intermediate cell states (Cook and Vanderhyden, 2020, 2022; Deshmukh et al., 2021; Pastushenko et al., 2018). These observations indicate that an improved theoretical foundation of the intermediate EMT states is needed to bridge the gap between models and experimental data. The improved models will in turn provide new insights into the developmental and pathological processes governed by epithelial-mesenchymal cell plasticity (Aiello et al., 2018; Pastushenko et al., 2018).

It was shown that EMT involves extensive post-transcriptional regulation by microRNA: more than one hundred miRNAs were found to be significantly associated with EMT (Cursons et al., 2018; Hussen et al., 2021). Both epithelial-associated microRNAs and mesenchymal-associated microRNAs have been identified, inhibiting mesenchymal genes (M genes) and epithelial genes (E genes), respectively. Mechanistically, this microRNAs-mediated inhibition often acts through mRNA degradation and translation repression upon microRNA-mRNA binding (Bartel, 2018). It was recently shown that multiple epithelial-associated microRNAs can inhibit in a strongly cooperative manner (Cursons et al., 2018), further supporting the roles of microRNAs in EMT. While existing models describing intermediate EMT states often include microRNA regulations, these models still require multiple transcriptional feedback loops to account for intermediate states (Hong et al., 2015; Lu et al., 2013b; Tian et al., 2013). It is therefore unclear whether intermediate EMT states can arise from simpler, more common gene regulatory networks. Recently, Li et al. used mass-action kinetics, a detailed description of biochemical reactions, to model microRNAs-mRNA reactions, including regulated degradation, with up to three microRNAs binding sites (Li et al., 2021). The models and subsequent experiments showed that these common microRNAs networks can generate two cell states (Li et al., 2021). This study and the extensive involvement of microRNAs in EMT suggest the possibility that reaction networks involving more than three microRNA-binding sites can be a hitherto unknown mechanism for intermediate EMT states.

In this work, we use models describing elementary RNA reaction networks to show that cooperative RNA degradation can generate intermediate EMT states in the absence of transcriptional feedback. We use bioinformatic approaches to demonstrate that the structures of gene regulatory networks allowing multiple intermediate states are widespread in the EMT system. Furthermore, we demonstrate that the transcriptional and post-transcriptional mechanisms can be combined to support larger numbers of intermediate EMT states in both modular and emergent manners. Finally, we use a large EMT model to show that cooperative RNA degradation can facilitate the formation of a phenotypic continuum. This model predicts the role of miR-101 in multistate EMT, a prediction we validate with recent single-cell RNA-sequencing data. Our work reveals a previously unknown mechanism for stabilizing intermediate EMT states and provides a new theoretical framework for understanding the perplexing EMT spectra in development and disease progression.

Results

Cooperative RNA degradation generates intermediate EMT states in the absence of transcriptional feedback

Based on recent theories and experiments showing post-transcriptional mechanisms for bistability (Li et al., 2021), we hypothesized that intermediate EMT states can arise without transcriptional feedback. To test the hypothesis, we first considered models with mass-action kinetics describing interactions between microRNAs and mRNAs as well as their synthesis and degradation. It was previously proved that with arbitrary positive rate constants, an mRNA with one microRNA-binding site (the MMI1 Model) can have only one stable steady state and that an mRNA with two microRNA-binding sites (the MMI2 Model) can have at most two stable steady states (Li et al., 2021). Systematic search with biologically relevant parameter sets showed that, like the MMI2 Model, an mRNA with three microRNA-binding sites (the MMI3 Model) can have at most two stable steady states (Li et al., 2021).

We therefore built a model containing an mRNA with four microRNA-binding sites: the MMI4 Model. In the first version of the MMI4 Model, the four microRNA-binding sites are bound by one type of microRNA. The two RNAs can form four types of complexes through complementarity-based binding, i.e. 1: complex where represents the number of microRNA molecules in each complex (Figure 1A, Table S1). The two RNAs are allowed to be degraded independently with distinct rate constants in the complexes. Details of all models can be found in the supplemental information. We randomly sampled parameter values within biologically plausible ranges (see STAR methods). Out of sampled parameter sets, we found 1,732 that generated three stable steady states (i.e. three attractors, i.e. tristability, Figures 1B and 1C) and none that generated four attractors (tetrastability). The third attractor at which the concentration of an EMT marker gene is neither minimal (E state) nor maximal (M state) for the system represents an intermediate state. The tristable parameter sets, i.e. those permitting an intermediate state, are generally robust to changes of at least 10% in any single relative RNA degradation rate (Figure S1). Their degradation rate constants, while individually plausible, are however collectively unusual. For example, the most extreme degradation rates are often found in unsaturated complexes (Table S2), which would require the mRNA to be strongly stabilized by the binding of additional microRNA molecules. While the structure of the model is applicable to EMT (all 7 matching target genes shown in Figure 1D and Figure S6A), the dependence on unusual parameter progression makes the single-microRNA MMI4 Model unlikely to drive multistability in vivo.

Figure 1

Multistability in the MMI4 Models

(A) Reactions in the one-microRNA MMI4 Model, which has four binding sites for the same microRNA on one mRNA. Wide rectangles, mRNA; squares, microRNA; horizontal arrows, transcription; colored arrows, RNA degradation; curved arrows, binding/unbinding.

(B) Example tristable systems from the one-microRNA MMI4 Model. Points represent attractors in the space of free mRNA vs. free microRNA concentration. Attractors of the same system/parameterization are joined by lines of the same color. AU, arbitrary units.

(C) Bifurcation diagrams showing the steady states of free mRNA (left) and free microRNA (right) as a function of the mRNA transcription rate from the brown system in B. Each steady state is colored the same in both plots. Dashed lines, unstable steady states. Parameter values in Table S2.

(D) Other EMT-related examples of the one-microRNA MMI4 Model.

(E) Reactions in the two-microRNA MMI4 Model, which has two binding sites for each of two microRNAs on one mRNA.

(F) Example tetrastable systems from the two-microRNA MMI4 Model in the space of free mRNA vs. free microRNA 1 (left) or free microRNA 2 (right). Order of attractors connected by each line is the same in both subplots.

(G) Bifurcation diagrams showing the steady states of free mRNA (left), free microRNA 1 (middle), and free microRNA 2 (right) as a function of the mRNA transcription rate from the green system in F. Parameter values in Table S4.

(H) Left: Scatterplot of functional cooperativities in mRNA degradation () and microRNA 1 degradation () rates due to second microRNA 1 binding in 5,000 3- or 4-attractor systems. Marginal distributions are Gaussian kernel density estimates. Right: Comparison of multistable cooperativity distributions for both microRNAs to distribution of 50,000 randomly sampled parameter sets.

(I) Other examples of the two-microRNA MMI4 Model in EMT genes.

See also Tables S1–S4 and Figures S1–S6.

We considered that allowing the four binding sites to be bound by different microRNAs allows the model to represent more post-transcriptional EMT circuits with potentially more realistic parameter sets. We therefore considered a modified version of the MMI4 Model which describes two microRNAs each with two binding sites on the mRNA (Figure 1D and Table S3). This structure indeed allowed the MMI4 Model to produce three- and four-attractor systems with more biologically plausible parameter progressions (Table S4) and low dependence on any individual degradation rate (Figure S3). These systems contain intermediate EMT states that may be related to the widely observed phenotypes (Figures 1F and 1G). To test the importance of functional cooperativity in stabilizing intermediate states, we randomly sampled parameter sets until 5,000 systems with at least three attractors had been obtained. Almost all, 4,973 of the 5,000, involved cooperative degradation of RNAs. The cooperativity is generated by the enhanced degradation rate constant of the mRNA in 1:2 complex with a microRNA compared to the 1:1 complex or by reduced degradation rate constant of the microRNA in 1:2 complex compared to the 1:1 complex (Figures 1H and S5). Cooperative mRNA degradation is observed experimentally (Grimson et al., 2007); functionally cooperative microRNA stabilization might arise through steric blocking of microRNA-degrading factors (Kai and Pasquinelli, 2010). Three instances of the two-microRNA MMI4 Model are shown in Figure 1I; all 45 EMT genes matching the model are shown in Figure S6B.

EMT is regulated by a large gene regulatory network, in which many genes are targeted by multiple microRNAs (Hong et al., 2015; Huang et al., 2017; Lu et al., 2013b; Tian et al., 2013). Furthermore, the MMI2 Model can generate bistable systems (Li et al., 2021). We therefore asked whether connecting two MMI2 modules (the Chained-MMI2 Model, Figure 2A and Table S5) can enable intermediate EMT states. While this single connection can represent an EMT transcription factor regulating another, this network still does not contain any transcriptional feedback. We found that the Chained-MMI2 Model was able to generate intermediate EMT states in terms of the target gene expression (Figure 2B). Considering only direct regulations between EMT genes found in the TRRUST2 or OmniPath databases, we found 8 instances of the Chained-MMI2 Model in the EMT network (Methods, Figure S9A). If indirect regulations of up to 5 steps (i.e. 4 intermediate genes) are allowed, up to 171 combinations of EMT genes match the model.

Figure 2

Tetrastability from transcriptionally connected MMI2 targets

(A) Reactions in the Chained-MMI2 Model, in which the protein product of one MMI2 target gene transcriptionally regulates another MMI2 target gene. Either microRNA site on each mRNA is allowed to be bound first. Hollow arrows, translation; dashed pointed arrow, transcriptional activation; stars, proteins.

(B) Example tetrastable systems from the Chained-MMI2 Model in the space of Protein 2 vs. Protein 1 (left) or free microRNA 2 (right). Each system has four different expression levels of Protein 2, which are monotonically anticorrelated to free microRNA 2 levels.

(C) Reactions in the Co-targeting-MMI2 Model, in which two MMI2 target genes encode proteins that both transcriptionally regulate a third downstream gene. Either microRNA site on each mRNA is allowed to be bound first. Dashed blunt arrow, transcriptional repression.

(D) Example tetrastable systems from the Co-targeting-MMI2 Model in the space of the downstream mRNA expression vs. Protein 1 (left) or Protein 2 (right). Each system has four different expression levels of the downstream gene.

See also Tables S5–S8 and Figures S7–S9.

We next considered another scenario where one EMT gene is regulated by two EMT transcription factors (Batlle et al., 2000; Vannier et al., 2013), each involved in an MMI2 module (the Co-targeting-MMI2 Model, Figure 2C and Table S7). With this model, intermediate EMT states were observed in terms of the expression of the final target EMT gene (Figure 2D), arising in a combinatorial manner from the multiple upstream bistable systems. The 18 direct-regulation instances of this model in the EMT network are shown in Figure S9B. If allowing indirect regulations of up to 5 steps, up to 1,312 sets of EMT genes match the Co-targeting-MMI2 Model.

In summary, we showed that a model with a total of four microRNA-binding sites on either one or two EMT genes can generate intermediate EMT states. In each of the four versions of the model, cooperative RNA degradation stabilizes these intermediate states. Intermediates-enabling topologies are common in the EMT system. In fact, 55 of 423 classically defined EMT genes are involved in direct-regulation versions of these topologies, so many seemingly unidirectional interactions between well-known EMT genes could potentially contribute to the multistate dynamics associated with feedback (see STAR methods).

Modular and emergent synergies between post-transcriptional and transcriptional networks in generating EMT states

It is well known that transcriptional positive feedback loops can generate bistability (Gardner et al., 2000). This type of feedback loop is common in the core EMT regulatory network (Hong et al., 2015). Furthermore, hybrid transcriptional and post-transcriptional mechanisms can generate tristability according to previous models (Lu et al., 2013a, 2013b). We next asked whether combining a module consisting of a transcriptional feedback loop (transcriptional module) with a post-transcriptionally driven, intermediates-enabling module (post-transcriptional module) can generate even more intermediate EMT states (Jiménez et al., 2017). To test this, we first selected a bistability-enabling parameter set for a representative transcriptional module containing two mutually activating transcription factors (Figure 3A Module 1, Table S9). We next selected a parameter set for a tetrastability-enabling post-transcriptional module, the Chained-MMI2 Model (Figures 2A and 3A Module 2, Table S5). Without altering the values of parameters unique to each module, we considered new values for the parameters that conflict, namely the transcription and translation rates of the upstream gene in the Chained-MMI2 Model. We found that there existed values between each pair of original models’ values that allowed the addition of at least one intermediate state to the existing states generated by the post-transcriptional module (Figure 3B). Even in 6-attractor systems, only the dependence of mRNA 0 transcription on Protein 1 is a very sensitive parameter for multistability (Figure S10). It is also remarkable that the addition of intermediate states was achieved without altering the biochemical rate constants within each module, suggesting feasibility of this phenotypic change through evolution.

Figure 3

Synergies between transcriptional and post-transcriptional multistability

(A) The model resulting from combining a transcriptional mutual activation module (genes 0 and 1, purple shading) with a Chained-MMI2 module (genes 1 and 2, orange shading).

(B) Example five- and six-attractor systems produced by combining bistable transcriptional mutual activation parameter sets with tetrastable Chained-MMI2 parameter sets. All concentrations are in arbitrary units.

(C) The model resulting from adding transcriptional repression of the microRNA to the MMI2 Model.

(D) Bifurcation diagrams showing the steady states of protein (left) and free microRNA (right) levels with respect to the mRNA transcription rate. Tristability can emerge from the addition of the transcriptional repression to MMI2.

(E) The addition of transcriptional repressions to the two-microRNA MMI4 Model. The mRNA-microRNA complexes are hidden for compactness.

(F) Example five-attractor systems emerging from the addition of transcriptional repressions to the two-microRNA MMI4 Model.

See also Tables S9–S15 and Figures S10–S13.

Modular combination of transcriptional and post-transcriptional mechanisms for generating multistable systems is applicable to EMT, but there might also exist “emergent” synergy between transcriptional and post-transcriptional regulations: a transcriptional module may not be multistable by itself, but when it is combined with a multistable, post-transcriptional module, additional attractor(s) can arise. Indeed, we obtained an intermediate state by extending the bistability-enabling MMI2 module with a single transcriptional repression of the microRNA by the gene product (Figure 3D). Mechanistically, the emergence of this intermediate cell state was due to the emergent feedback loop between the microRNA and mRNA, consisting of both transcriptional and post-transcriptional regulations (Figure 3C and Table S12). It was previously shown that this type of hybrid feedback system is common in biology (Minchington et al., 2020). Importantly, the well-known Zeb1-miR200 feedback loop contains this network topology (Bracken et al., 2008). Our results indicate that the previously known tristability-enabling structure of the Zeb1-miR200 feedback loop and Zeb1 self-regulation is not the minimal topology for generating an EMT intermediate state (Lu et al., 2013b; Tian et al., 2013). Likewise, we found that adding transcriptional repression to the two-microRNA MMI4 Model permits the appearance of a fifth attractor (Figures 3E, 3F, and Tables S14, and S15).

In summary, our results suggest that the post-transcriptional and transcriptional networks commonly observed in the EMT system can be combined in both modular and emergent fashions to generate additional intermediate phenotypes.

A multiple-feedback-mechanism EMT network generates a 7-state EMT continuum

How does the RNA degradation-based mechanism for multistability contribute to the formation of intermediate EMT states in networks larger than simple motifs? To address this question, we first selected an EMT model published recently: Subbalakshmi et al. showed that a network containing three transcription factors and thirteen transcriptional regulations can generate a three- or four-attractor EMT system with biologically meaningful parameters (Subbalakshmi et al., 2022). While there are numerous existing EMT models for tetrastable systems, we chose the Subbalakshmi et al. model because only one microRNA was considered in the model, so it serves as a good basal model for us to test the effect of cooperative RNA degradation by adding more biologically relevant microRNA-binding sites sequentially to the system.

As the base model contains inhibition of ZEB1 and SNAI2 by miR-200 without explicit modeling of multiple complexes, we added an additional binding site on ZEB1 to the model (Figure 4A and Table S16) and selected a parameter set representing cooperative RNA degradation. As expected, the model with two binding sites of miR-200 produced a system with five attractors (Figures 4B and 4C), i.e. an additional intermediate EMT state compared to the previously published results. We next considered another microRNA, miR-101, that regulates ZEB1 (Guo et al., 2014) and is repressed by slug and snail (Huang et al., 2017; Zheng et al., 2015) (Figure 4D and Table S18). Note that these additional regulations do not introduce any additional transcriptional feedback loops. The inclusion of this microRNA allowed yet two additional intermediate EMT states, resulting in a seven-attractor EMT system (Figure 4E blue).

Figure 4

Combining transcriptional and post-transcriptional regulation produces many attractors and a continuum in epithelial-mesenchymal space

(A) Structure of the Subbalakshmi et al. model extended with two explicitly modeled miR-200 binding sites on ZEB1. mRNA-microRNA complexes are hidden for compactness.

(B) Example 5-attractor systems from the model in A shown in the gene expression space of Zeb1 protein vs. E-cadherin protein, Slug protein, or free miR-200.

(C) Quasi-potential diagram showing the stochastic gene expression landscape of the blue system in B under multiplicative noise in the space of E-cadherin protein vs. Zeb1 protein. Deeper regions and brighter color correspond to more likely gene expression states. Spheres, deterministic attractors.

(D) Structure of the Subbalakshmi et al. model further extended with miR-101 targeting one site on ZEB1, transcriptionally repressed by Slug and Snail.

(E) 7-attractor (blue), 6-attractor (orange), and 5-attractor (green) systems from the model in D. All concentrations are in arbitrary units.

(F) Quasi-potential (QP) diagram of the 7-attractor blue system in E.

(G) Quasi-potential diagram of the 5-attractor green system in E. Not every parameter set of the extended two-microRNA model exhibits the broad distribution found in F.

See also Tables S16–S19, Figures S14, and S15.

To test the effect of cooperative RNA degradation with transcriptional noise, we performed stochastic simulations of the one-microRNA Subbalakshmi et al. network (Figures 4A–4C) and the extended model with two microRNAs (Figures 4D–4G), applying the same level of multiplicative noise to all RNAs, complexes, and proteins. Based on the steady-state distributions of molecules from at least 480 simulations each representing a cell from a population, we constructed the quasi-steady-state landscape to visualize the multi-attractor systems under the influence of noise. As expected, cooperative RNA degradation via multiple binding sites gave rise to gene expression states near the additional attractors (Figure 4F yellow). Notably, it also resulted in broader distributions of gene expression further from the attractors (Figure 4F orange, compare to Figure 4C). That is, a deterministic attractor—which may be highly sensitive to several aspects of Zeb1 regulation (Figure S15)—is not required for a region of gene expression space to be populated in the presence of noise. While the wider distribution of epithelial and mesenchymal marker genes is not simply an artifact of simulating the more complex model (Figure 4G), the new complexes involving the additional microRNA do generally create more axes in concentration space. Fluctuations in gene expression space can then be amplified in impact by functional cooperativity in degradation rates, pushing the system toward a different state.

Our new model (Figure 4D) proposed a role of miR-101 and its stepwise changes of expression level in EMT progression. To test this prediction, we next analyzed two single-cell RNA-sequencing (scRNA-seq) datasets from TGF-β-driven EMT. Transcriptome data for 3,568 A549 cells from a time course experiment (Cook and Vanderhyden, 2020) and for 11,220 MCF10A cells from a dose-dependent, nearly steady-state experiment (Panchy et al., 2022), were visualized using Uniform Mani-fold Approximation and Projection (UMAP). Consistent with earlier analyses (Cook and Vanderhyden, 2020; Panchy et al., 2022), cells were distributed in continuous space for each dataset (Figure 5A), suggesting an EMT continuum. Furthermore, median expression of ZEB1 positively correlated with the progression of EMT, whereas that of CDH1 had the opposite trend (Figures 5B and 5C) (the negative correlation of CDH1 in A549 cells was not significant due to low detection rate). Interestingly, the mean target gene expressions for both miR-200b/c and miR-101 were positively correlated with the progression of EMT, suggesting stepwise changes of their expression levels (Figures 5D and 5E). These effects were not due to the time/dose-dependent changes of overall gene expression patterns (see STAR methods), and they are consistent with the predictions from our model (Figure 4E).

Figure 5

EMT continuums and progression of microRNA activities validated by scRNA-seq data

Two scRNA-seq datasets (GSE147405 and GSE213753) for A549 and MCF10A respectively were visualized with UMAP and expression profiles of individual genes.

(A) UMAP plots show distributions of 3568 A549 cells and 11,220 MCF10A cells. Colors indicate treatment conditions.

(B and C) Boxplots show medians and interquartile ranges of ZEB1 and CDH1 expression under each condition. ∗ indicates the exclusion of cells with no expression. r denotes the Pearson correlation coefficient for quantities on both axes.

(D and E) Boxplots show medians and interquartile ranges of miR-200c and miR-101 activities inferred from their target gene expressions (see STAR methods).

In summary, in a larger EMT network containing both transcriptional and post-transcriptional regulations, cooperative RNA degradation via multiple microRNA-binding sites gave rise to additional attractors and a broader distribution of gene expression, reflecting the EMT continuum observed in recent single-cell transcriptome data.

Discussion

Intermediate, or hybrid, EMT phenotypes have been widely observed in several biological contexts. In this work, we have used mathematical models to demonstrate a new mechanism for generating intermediate EMT states based on first principles of gene regulation. This finding can serve as a step toward the reconciliation of the observed EMT continuum with transcriptomic studies and the three or four discrete EMT states captured by previous models. While the observed EMT continuum may be alternatively explained by large subpopulations of cells en route to different attractors, recent work using cell-state transition models to explain experimental data showed that models with only a few states cannot describe the time course EMT data accurately (Goetz et al., 2020). The newly identified post-transcriptional mechanism for generating intermediate states provides a foundation for the additional EMT states necessary to explain expression data at the gene regulation level. A recent study showed that similar post-transcriptional reactions can generate oscillations on slow timescales in addition to multistability (Nordick et al., 2022), which suggests another possibility that the EMT continuum may be supported by a combination of point attractors and cyclic attractors.

Cooperativity of microRNA-binding sites has been widely observed (Buck et al., 2014; Cursons et al., 2018; Lai et al., 2019). In particular, Cursons et al. showed high cooperativity between multiple microRNAs in controlling EMT (Cursons et al., 2018). In this work, we modeled the cooperativity of microRNA-binding sites in the form of synergistic mRNA degradation, which was also experimentally observed, and synergistic microRNA stabilization. There are other forms of biologically plausible synergy between binding sites, such as cooperative binding affinities, cooperative translational inhibition, and competition for one microRNA’s binding to 3′UTRs of multiple mRNAs (e.g. miR-200 in the last model of this study) (Briskin et al., 2020; Grimson et al., 2007). We expect that some of these mechanisms can be used to support stable intermediate cell phenotypes. Future work is necessary to systematically compare the functions of these molecular mechanisms and to identify their existence in specific biological contexts. Nonetheless, the prevalence of the multi-site microRNA interactions with individual and groups of mRNAs in EMT and other systems suggests that these regulatory networks have nontrivial emergent functions.

In this study, we showed similar and modularizable performances of transcriptional and post-transcriptional mechanisms in generating intermediate EMT states. These two mechanisms are different in their cellular locations. While it may be beneficial for cells to combine both nuclear (transcriptional) and cytosolic (post-transcriptional) machineries to achieve the desired goal of stabilizing intermediate states, the post-transcriptional mechanism may be advantageous in terms of avoiding some sources of noise. This is because transcription is subject to significant noise levels due to the low numbers of DNA coding for the regulatory products, whereas post-transcriptional mechanisms involve larger numbers of molecules (Kataruka et al., 2022), which reduce intrinsic noise. Therefore, we expect that the proposed RNA-centric mechanism for stabilizing intermediate EMT states can be an efficient strategy for cells to adopt hybrid phenotypes with cytosolic reactions without the need for transcriptional regulatory systems.

The motivation to build our models derived from inconsistencies between existing EMT models that predict a paucity of EMT intermediate states, and experimental single-cell transcriptomic data that have been interpreted to support a wealth of states in a phenotypic continuum (Cook and Vanderhyden, 2020; Deshmukh et al., 2021; Panchy et al., 2022). It should be noted that some continuum-like expression distributions revealed by scRNA-seq may be due to technical noise, whose effects on EMT progression warrant future investigation. Nonetheless, it is plausible that some of the intermediate states should be favored, perhaps in relation with microenvironmental factors such as nutrient availability and cytokines. The models we present here can guide experimentation designed to validate the role of microRNAs to stabilize a constrained number of intermediate cell phenotypes both in physiologic and pathologic systems beyond EMT. For example, in small-cell lung cancer, a phenotypic continuum spanning neuroendocrine (NE) and non-neuroendocrine (non-NE) cell states was recently described based on archetype analysis (AA) of experimental data (Groves et al., 2022; Hausser and Alon, 2020). NE to non-NE transitions bear many similarities to EMT, particularly as it pertains increased metastatic properties of non-NE cells, and the similarity of transcriptional signatures (Ireland et al., 2020). Similar to the EMT phenotypic continuum, AA identified a phenotypic continuum for SCLC NE and non-NE subtypes, which gene expression enrichment links to cellular task they are optimized for, such as secretion, proliferation, or motility. In this continuum, cells may be specialists at one task or suboptimal generalists at one or more tasks. Thus, the intermediate generalist phenotypes arise from task trade-offs, so that they can perform suboptimally to be in tune with the microenvironment that they experience. The continuum then becomes dominated by Pareto optimality. Interrogating the role of microRNAs in NE to non-NE transitions intermediates will provide mechanistic underpinnings for plasticity and high metastatic propensity of SCLC tumors. Our modeling approach is a valuable starting point for streamlining these experiments. In summary, future studies may reveal that post-transcriptional mechanisms are widely used by mammalian cells for generating intermediate states, both for EMT and other differentiation systems. Our mathematical models will aid in designing experiments to test this possibility.

Limitations of the study

While it is expected that several gene regulatory mechanisms can generate multistable systems (Angeli et al., 2004; Li et al., 2017; Markevich et al., 2004), the predicted co-existing cell states from cooperative RNA degradation models were not validated experimentally in this study, and they warrant future work. The observed EMT continuum at the transcriptome level (Cook and Vanderhyden, 2020; Deshmukh et al., 2021; Panchy et al., 2022) does not exclude the possibility that only a few cell states emerge post-translationally to govern cellular phenotypes in a more discrete manner. Future work is needed to quantitatively characterize the continuity of EMT spectrum in terms of functional phenotypes.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Tian Hong (hongtian@utk.edu).

Materials availability

This study did not generate new unique reagents.

Method details

Model construction

Each model is described by a system of ordinary differential equations (ODEs). Dynamics of molecular species are described by two types of function. Mass-action kinetics is used to model elementary reactions, i.e. the basal processes of RNA transcription/maturation, constitutive RNA decay, mRNA-microRNA binding and unbinding, regulated decay of each RNA member of each complex, translation, and protein decay. If applicable to the model, Hill functions are used to describe regulated transcription rates as in HiLoop models (Nordick and Hong, 2021):for repression, orfor activation, where is the concentration of the regulator protein, is the threshold of regulation, is the cooperativity of regulation, and is the proportion of maximum transcription rate controlled by this regulator. Multiple regulatory influences are combined by adding their Hill functions.

The models based on the Subbalakshmi et al. network follow the form used there and by Lu et al. (Lu et al., 2013b; Subbalakshmi et al., 2022)whereand is the fold-change in target expression resulting from overexpression of the regulator. Multiple regulations are combined by multiplying their functions.

In all models, each type of RNA within each type of complex is assigned a regulated degradation factor, for mRNAs or for microRNAs, representing the change in its degradation rate in the complex relative to its free state (Nordick et al., 2022).

Binding and unbinding in RNA complexes is always modeled explicitly to ensure that each microRNA molecule acts on only one mRNA at once, though Subbalakshmi et al.’s allowance of incomplete translation repression by microRNA binding (Huntzinger and Izaurralde, 2011) is followed. For all mRNA-microRNA binding events in all models, the microRNA dissociation rate and each association constant .

For mRNAs with multiple binding sites for a type of microRNA, any binding/unbinding order is allowed, but complexes that differ only in the positions of sites bound by their microRNAs are grouped together into one variable, which represents the total concentration of all forms of the complex. This is accomplished by multiplying the single-microRNA association rate by the number of remaining free sites and multiplying the dissociation or microRNA degradation rate by the number of occupied sites.

The reactions, rate functions, and parameters in each model are listed in the supplementary information. Antimony model code for use with the Tellurium modeling package (Choi et al., 2018) can be found in the associated GitHub repository. All variables and parameters are dimensionless quantities (Nordick et al., 2022).

For simulation, the set of reactions is converted to a system of ODEs, one for each RNA, protein, or complex. A species’ rate of change is the sum of the rates of the reactions that produce or consume it, weighted by the net change in the species’ amount caused by the reaction. As an example, the ODEs for the one-microRNA MMI4 Model (derived from Table S1) arewhere is the concentration of free mRNA, is the concentration of free microRNA, and is the total concentration of complexes with microRNA molecules bound to an mRNA molecule except and .

One reaction can contribute terms to multiple species’ rate-of-change equations. The regulated mRNA decay reactions, for example, which occur at rate , destroy one unit of complex but produce units of free microRNA , so they contribute a term to and a term to .

Testing model parameterizations for multiple attractors

Each parameterized system was simulated starting from at least 100 initial conditions which span at least five powers of 5 in a log-uniform manner, drawn from a Sobol quasirandom sequence (Bratley and Fox, 1988). The initial concentration of each free RNA species was determined by one dimension of the Sobol hypercube, protein initial concentrations were set to their coding mRNA’s initial concentration times the translation rate, and other species’ initial concentrations were zero. Deterministic simulation with Tellurium (Choi et al., 2018) proceeded for at least 100 time units up to a maximum of 1,000 time units until the system reached a steady state as determined by the 2-norm of the derivatives vector falling below . Simulation endpoint concentration vectors and were considered equivalent attractors if their difference’s 2-norm was less than times the number of species , or if .

Numerical bifurcation analysis

Bifurcation diagrams were created with the AUTO2000 plugin for Tellurium 2.2.0 (Choi et al., 2018). For bifurcations with respect to transcription rate, continuation proceeded backwards from a high monostable signal value. For bifurcations with respect to degradation rate, continuation proceeded forward from a low degradation rate and was restarted from a different steady state when needed to follow all attractors.

Stochastic simulation and quasi-potential landscape

We performed stochastic simulations for the modified 3-TF model with various microRNA binding sites by adding an independent multiplicative noise term to each ODE. Divergence arising from negative concentrations was reduced by applying the Zero-Reaction remedy (Chen and Cao, 2019). Starting with a population of an equal number of cells at each deterministic attractor, we solved the stochastic ODEs at a noise intensity of 0.2 for 200 time units using DifferentialEquations.jl (Rackauckas and Nie, 2017).

To visualize the epithelial-mesenchymal gene expression space in which the stochastic system fluctuates in the long term, the concentrations of Zeb1 and E-cadherin protein were extracted for each simulated cell at intervals of 5 time units starting at time 150. This pair of molecules were chosen because they were very widely studied and are considered core effectors/regulators of EMT (Celià-Terrassa et al., 2018; Sanchez-Tillo et al., 2010; Wellner et al., 2009). The concentrations of protein-coding genes are readily inferred from RNA-sequencing data. By contrast, microRNA concentrations are more difficult to obtain. To avoid distortions from temporarily negative concentrations, timepoints with either component less than were not used. Quasi-potential landscapes (Li and Wang, 2014) were rendered with potential , where is the probability density function computed by scikit-learn’s 2-dimensional Gaussian kernel density estimate (Pedregosa et al., 2011) of bandwidth 0.14 on the base-2 logarithm of E-cadherin protein and the base-10 logarithm of Zeb1 protein, as Zeb1’s expression was more variable under the parameter values tested.

Single-cell RNA-sequencing data origin

We analyzed two datasets for validating some predictions from our final model that includes miR-101. Single-cell RNA-sequencing data for TGF-β induced EMT in A549 lung epithelial cell line (Cook and Vanderhyden, 2020) (GSE147405), and in MCF10A breast epithelial cell line (Panchy et al., 2022) (GSE213753) were obtained from NCBI Gene Expression Omnibus (accessions in key resources table). The A549 cells were treated with a single dose of TGF-β and their data have time labels of 0-7 days. The MCF10A cells were treated with multiple doses (0-200 p.m. chosen for this work) of TGF-β for 14 days.

Quantification and statistical analysis

Parameter sampling for multistable systems

The likelihood of the single-microRNA MMI4 Model to generate multiple attractors was tested by sampling all regulated degradation factors independently from a log-uniform distribution on (a range consistent with previous experiments and models (Eichhorn et al., 2014; Nordick et al., 2022)), the microRNA association constant from a log-uniform distribution on , and the microRNA transcription rate from a log-uniform distribution on . The conditions for the two-microRNA MMI4 Model to produce multiple attractors were investigated by sampling from the same parameter regions independently for each microRNA, discarding systems with at most two attractors. Kernel density estimates were computed by SciPy (Virtanen et al., 2020) with standard settings.

When simply searching for example parameter sets in other models rather than characterizing the parameter space, subregions that provided computational efficiency were selected empirically from plausible regions similar to previous work (Nordick et al., 2022); all sampling distributions and parameter sets of interest for each model are listed in supplementary tables.

Parameter sensitivity analysis

To determine how much a parameter must change to abolish high multistability, the system was repeatedly re-simulated as described above. First, the original parameter value was doubled or halved up to twelve times until the perturbed parameter set had fewer attractors than the original. To find the precise limit, a binary search was then conducted between the last many-attractor and first fewer-attractor power of 2. The binary search was terminated when the base-2 logarithms of the narrowed parameter value endpoints differed by less than 0.0001.

To address derailment of the binary search from numerical imprecision in attractor detection, values 1% above and below each endpoint were re-simulated for verification. If a value inside the many-attractor range was found to lack an attractor, or if a value outside the range was found to still have many attractors, the limit search process was repeated with more stringent integrator settings. The initial time was increased to 500 units, the maximum time was increased to 4000 units, and the initial timestep was decreased to 0.25. Parameter limits that still failed verification, i.e. could not be determined to within 1%, were discarded.

Enumeration of instances of network topologies

We examined a previously curated list of core EMT genes (Tan et al., 2014). Among them, 232 were classified as E genes and 191 as M genes. In addition, we used a list of 133 microRNAs, each with experimental evidence supporting its role in EMT (Cursons et al., 2018). With these lists and the TargetScan program for microRNA binding site prediction (Agarwal et al., 2015), we first identified 46 EMT genes regulated by a total of four binding sites of one or two EMT microRNAs, i.e. the MMI4 Models. Including those, we also identified 93 EMT genes matching the MMI2 Model.

For the topologies that combine transcriptional and posttranscription regulation, the TRRUST2 network (Han et al., 2018) and transcriptional subgraph of the OmniPath network (Türei et al., 2016) as preprocessed for HiLoop (Nordick and Hong, 2021) were added together except for regulations given opposite signs by the two networks. The existence of a regulatory path from each EMT MMI2 gene to every other, up to the specified path length of 1 or 5, was tested in the combined transcriptional network using NetworkX (Hagberg et al., 2008). Each ordered pair with such a directed path was considered an instance of the Chained-MMI2 topology.

The list of EMT MMI2 genes that, directly or indirectly via the combined network, transcriptionally regulates each EMT gene was similarly obtained. If an unordered pair of EMT MMI2 genes regulated the same EMT target gene without common dependencies on any intermediate genes, the partially ordered triple of regulators and target was considered an instance of the Co-targeting-MMI2 topology. The condition of no common intermediates avoids counting subnetworks in which one MMI2 gene regulates the target through the other MMI2 gene, an arrangement which may have different dynamics, and non-minimal subnetworks in which the target is downstream of an already co-targeted gene.

Note that these distinct instances are only based on different combinations of EMT genes. Considering combinations of EMT-associated microRNA binding sites would yield much larger numbers of instances.

Single-cell RNA-sequencing data analysis

The Uniform Mani-fold Approximation and Projection (UMAP) was used to examine the overall distributions of cells for both datasets. Z-scores were computed for each gene across each dataset. Representative EMT genes were used to visualize the progression of EMT via their expression across time points and doses. For genes that were poorly detected (ZEB1 and CDH1 in A549 cells, and ZEB1 in MCF10A cells), only cells with non-zero expression were used for this visualization. To infer microRNA expression levels, we computed the mean expression levels of the genes potentially targeted by miR-200b/c, and miR-101 respectively. We obtained lists of predicted target genes from TargetScan (Agarwal et al., 2015). To avoid the trivial case in which the mean target gene expression was merely driven by ZEB1 (regulated by the microRNAs in our model) and its downstream effector gene VIM, we excluded these two genes from this inference analysis. After the exclusion, miR-101 and miR-200b/c had 954 and 1194 predicted target genes respectively. We normalized the mean target gene expression by subtracting the overall mean expression for each time point or dose from the target gene expression.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Single-cell RNA-seq data for TGF-β induced EMT in A549 lung epithelial cell line	(Cook and Vanderhyden, 2020)	GEO: GSE147405
Single-cell RNA-seq data for TGF-β induced EMT in MCF10A breast epithelial cell line	(Panchy et al., 2022)	GEO: GSE213753
Software and algorithms
Multistable model parameterization search scripts	This paper	https://github.com/BenNordick/MMI4
Tellurium 2.2.0	(Choi et al., 2018)	https://tellurium.analogmachine.org/
SciPy 1.1.0	(Virtanen et al., 2020)	https://scipy.org/
NetworkX 2.5.1	(Hagberg et al., 2008)	https://networkx.org/
Scikit-learn 0.20.0	(Pedregosa et al., 2011)	https://scikit-learn.org/
Python 3.8.5	Python Software Foundation	https://python.org/