Coupling the modules of EMT and stemness: A tunable 'stemness window' model.

Metastasis of carcinoma involves migration of tumor cells to distant organs and initiate secondary tumors. Migration requires a complete or partial Epithelial-to-Mesenchymal Transition (EMT), and tumor-initiation requires cells possessing stemness. Epithelial cells (E) undergoing a complete EMT to become mesenchymal (M) have been suggested to be more likely to possess stemness. However, recent studies suggest that stemness can also be associated with cells undergoing a partial EMT (hybrid E/M phenotype). Therefore, the correlation between EMT and stemness remains elusive. Here, using a theoretical framework that couples the core EMT and stemness modules (miR-200/ZEB and LIN28/let-7), we demonstrate that the positioning of ‘stemness window’ on the ‘EMT axis’ need not be universal; rather it can be fine-tuned. Particularly, we present OVOL as an example of a modulating factor that, due to its coupling with miR-200/ZEB/LIN28/let-7 circuit, fine-tunes the EMT-stemness interplay. Coupling OVOL can inhibit the stemness likelihood of M and elevate that of the hybrid E/M (partial EMT) phenotype, thereby pulling the ‘stemness window’ away from the M end of ‘EMT axis’. Our results unify various apparently contradictory experimental findings regarding the interconnection between EMT and stemness, corroborate the emerging notion that partial EMT associates with stemness, and offer new testable predictions.


REGULATORY CIRCUIT DETAILS EMT regulatory network
EMT regulatory network includes two highly interconnected modules -miR-34/SNAIL and miR-200/ZEB, both of which are mutual inhibitory circuits [1][2][3]. miR-34/SNAIL has been shown to act as a noise-buffering integrator, while miR-200/ZEB acts as a three-way switch allowing the co-existence of three phenotypes -epithelial (E -high miR-200, low ZEB), hybrid epithelial/ mesenchymal (E/M -intermediate miR-200, intermediate ZEB), and mesenchymal (M -low miR-200, high ZEB), and acts as the decision-making module during EMT. Coupling miR-34/SNAIL to miR-200/ZEB does not change the qualitative behavior of EMT decision-making [4,5]; hence we treat SNAIL as the external input signal to the miR-200/ZEB circuit.
The miR-200 family contains two subgroups based on the seed sequences -one subgroup comprises miR-141 and miR-200a, and the other includes miR-200b, miR-200c, and miR-429. The transcription factor ZEB family includes two isoforms -ZEB1 and ZEB2. The 3′ UTR region of ZEB1 has a total of eight conserved binding sites for miR-200 -three for the first subgroup and five for the second subgroup. The 3′ UTR region of ZEB2 has nine conserved binding sites for miR-200, specifically three for the first subgroup and six for the second subgroup [6]. Experimental evidence suggests that the expression of miR-200c standalone can induce MET and restore the expression of E-cadherin [7]. Therefore, here, in our miR-200/ZEB module, we considered six binding sites on the 3′ UTR region of ZEB for the binding of miR-200 and three binding sites on the promoter region of miR-200 for the binding of ZEB. In addition, ZEB can also activate its own transcription through SMAD complexes [8] and we assume two binding sites in the promoter region of ZEB for its self-activation. SNAIL is a well-known EMT inducer, which can transcriptionally inhibit miR-200 and activate ZEB [2,9]. Both these regulations have been assumed to happen through two binding sites [4].

Stemness regulatory network
Stemness is governed by a mutually inhibitory loop between RNA-binding factor family LIN28 and microRNA family let-7. Being RNA-binding factors, LIN28A and LIN28B, the two members of LIN28 family, can repress the biogenesis of let-7 by blocking its processing to mature miRNAs. LIN28B inhibits the processing of primary let-7 transcripts by the microprocessor Drosha, and LIN28A recruits a TUTase (Zcchc11/ TUT4) to block the processing of let-7 precursors by Dicer in the cytoplasm [10,11]. Being a microRNA, let-7 can inhibit the translation process of LIN28 [10]. Further, let-7 can promote its own processing and LIN28 can promote its own translation by binding to seven consensus sites on its own mRNA [12,13].
LIN28/let-7 circuit also behaves as a three-way switch -with the three states being U (Up -high LIN28, low let-7), D (Down -low LIN28, high let-7) and D/U (Down/Up -medium LIN28, medium let-7) [14]. LIN28 promotes the translation of pluripotency factor OCT4 [15]. Both very high and very low levels of OCT4 lead to differentiation, but its intermediate levels are strongly associated with gaining stemness [16][17][18][19]. Therefore, cells with intermediate levels of LIN28 (or D/U state) are most likely to lie in the 'stemness window' or correspond with the maximum likelihood of gaining stemness.

Coupling of the two core networks
These two core networks are coupled in two ways: miR-200 inhibits LIN28 by binding to its mRNA, and let-7 inhibits HMGA2 that activates ZEB, hence indirectly repressing ZEB [20][21][22].

MODEL FORMULATION
In the combined circuit of miR-200/ZEB/LIN28/let-7, there are multiple modes of regulation -transcriptional, translational, miRNA-mediated regulation, and regulation of miRNA processing. For all modes except miRNAmediated regulation with a large number of binding sites of miRNA on its target mRNA, we use shifted Hill functions H S (X, λ) defined as H S (X, λ) = H − (X) + λH + (X), a weighted sum of positive and negative Hill functions. Parameter λ in the shifted Hill functions is the weight factor that represents fold-change in production rate from its basal level, due to the binding of regulatory factor. For activation, λ > 1 and shifted Hill function is represented by H S+ ; for repression, 0 < λ < 1 and shifted Hill function is represented by H S− ; and for no change, λ = 1. λ X, Y denotes the effect of X on Y. Transcriptional regulation has been canonically modeled using Hill function [23], and inhibition of let-7 by LIN28 has been shown experimentally to behave like an inhibitory Hill function [24].
The terms representing miRNA-mediated translational effects capture both actions of miRNA -mRNA degradation and mRNA sequestration so as to prevent translation. Since a microRNA is usually 22nt long, and the length of the seed sequence it recognizes on a mRNA is only 7~8 nt, here we assume that the bindings of different microRNAs are independent. We assume the binding rate between microRNA and mRNA is r + and

SUPPLEMENTARY INFORMATION
www.impactjournals.com/oncotarget/ Oncotarget, Supplementary Materials 2015 the unbinding rate is r − , and that these binding/unbinding processes between microRNA and mRNA are much faster as compared to production and degradation of proteins, therefore miRNA and mRNA are always in equilibrium, i.e.
where γ mi is the individual degradation rate of mRNA bound to i molecules of miRNA. Similarly, the term denoting the degradation rate of microRNA due to binding to many mRNAs is mY μ (μ) = a n i=0 iγ μi C i n [m i ] = m a n i=0 iγ μi C i n M i n (μ) , where γ μi is the individual degradation rate for a microRNA molecule when i molecules of microRNA are bound to a mRNA. For small number of binding sites (~2) of microRNA on its target mRNA, the net silencing effects of microRNA can be approximated as a Hill function [4]. Therefore, for the cases when let-7 inhibits LIN28 and ZEB, and when miR-200 inhibits LIN28, shifted Hill functions are used to represent the effects of microRNAs.
On including OVOL, two more components -OVOL mRNA(m_ O ) and OVOL protein (O) are also incorporated.
The dynamics of miR-200 (μ 200 ) can be described by the following equation: The dynamics of ZEB mRNA (m Z ) and ZEB protein (Z) are described by the following equations: where g m Z and g Z are the innate production rates of ZEB mRNA and ZEB protein respectively, and k m Z and k Z are their respective innate degradation rates. H S+ 1 Z, λ Z, m Z 2 denotes transcriptional self-activation of ZEB, and H S+ 1 S, λ S, m Z 2 denotes transcriptional activation of ZEB by SNAIL. Y m 1 μ 200 2 represents the degradation of ZEB mRNA due to forming mRNA-miRNA complexes with miR-200, L 1 μ 200 2 denotes the translational inhibition of ZEB by miR-200, and H S− 1 u l , α u l , m Z 2 represents the inhibition from let-7.
When including OVOL, equation (1) for dynamics of miR-200 and equation (2) for ZEB mRNA are changed accordingly as follows: where Dynamics of OVOL mRNA (m O ) and OVOL protein (O) are described by these equations: where g m O and g O are innate production rates of OVOL mRNA and protein respectively, and k m O and k O are their respective degradation rates.
The dynamics of let-7 (μ l ) is described by the following equation: where g μ l and k μ l are the innate production and degradation rates of let-7 respectively, H S+ 1 μ l , λ μ l , μ l 2 represents the self-activation of let-7, H S− 1 L, λ L, μ l 2 represents inhibition by protein LIN28 and H S+ (N, λ N, μ l ) denotes activation by external signal (protein NF-kB (represented by N)).
Dynamics of LIN28 mRNA (m l ) and protein (L) are described by the following equations:

PARAMETERS ESTIMATION
The innate degradation rates of miRNAs, mRNAs, and proteins are selected based on their experimentally measured half-lives. Typical half-life of mammalian proteins is around 10 hours [25], so we chose 0.1/hour as the degradation rate for protein ZEB, OVOL and LIN28. The half-life of mRNA is around several hours [26], so we chose 0.5/ hour as the degradation rate for the mRNAs of ZEB, OVOL and LIN28. Generally, microRNA is more stable than mRNA [27,28], therefore we chose 0.05/hour as the degradation rate for miR-200 and let-7.
The expression levels of various molecules were estimated based on the typical concentrations and ratios in eukaryotic cells. The typical length scale for dimensions of a mammalian cell is around 10μm and typical protein concentration is in the range 10 nM~1 μM [29]. Hence, there are around 6 × 10 23 × 10 −6 /10 −3 × (10 × 10 −6 ) 3 or about one million molecules. Because the typical ratio of protein/mRNA is around 2800 [30], the level of mRNA should be in range of thousands of molecules. In addition, the number of microRNA is around ten thousand molecules [31]. Based on these typical levels, we chose production rates accordingly. The typical translation rate for one gene is around 140 proteins/(mRNA*hour) [30], here we chose 100 proteins/(mRNA*hour) for ZEB and 200 proteins/(mRNA*hour) for LIN28. The coupling strength parameters -α1 and α2 -vary from 0 (no coupling) to 1 (very strong coupling). All parameters used in the model have been detailed in Tables S1-S5.

PARAMETERS SENSITIVITY ANALYSIS
To test the sensitivity of our predictions based on the parameters we listed above, we conduct parameter sensitivity analysis by varying each parameter at one time. All the parameters -production rates, degradation rates, thresholds and weight factors in the shifted Hill functions are varied by ± 20%. The number of binding sites for the different interactions have been kept fixed because most of them are directly determined from experimental data.
First, we conduct parameter sensitivity for EMT regulatory circuit -miR-200/ZEB. For some values of the driving signal SNAIL, this circuit typically has three states -epithelial (E), hybrid epithelial/mesenchymal (E/M), and mesenchymal (M). Here, we measure how changes in the parameters affect the range of SNAIL levels for which the hybrid E/M phenotype exists (either alone in or combination with other phenotypes). Not surprisingly, the absolute values of SNAIL levels for the existence of E/M phenotype changes with every parameter change ( Figure S1), but here we focus on the range of SNAIL levels (i.e. the difference in the minimum and the maximum values of SNAIL) for which the cells can attain any of the three phenotypes, and not the absolute levels of SNAIL. Comparing with the control case (case 0, denoting the parameters listed above in the Tables S1-S5), the changes in the production rates and degradation rates of all species don't affect the range of SNAIL levels for which the hybrid state exists very much, as shown in Figure  S1A. However, the range of SNAIL levels for which the hybrid E/M state exists decreases when threshold of the shifted Hill function representing the inhibition by ZEB on miR-200 decreases (case G2), threshold of the shifted Hill function representing the self-activation of ZEB increases (case H1), the weight factor in the shifted Hill function representing the self-activation of ZEB increases (case N1) and weight factor in the shifted Hill function representing the activation on ZEB by SNAIL increases (case O1). All these cases represent an increase in effective ZEB levels that can inhibit miR-200 and hence drive EMT more strongly. In contrast, the range of SNAIL levels for which the hybrid E/M state exists increases when the threshold in the shifted Hill function representing the inhibition of miR-200 by ZEB increases (case G1), the threshold in the shifted Hill function representing ZEB self-activation decreases (case H2), the weight factor in the shifted Hill function representing ZEB self-activation decreases (case N2) and the weight factor in the shifted Hill function representing the activation of ZEB by SNAIL decreases (case O2). During all these cases, the effective ZEB levels that can inhibit miR-200 decreases, which result in a higher effective miR-200 levels and contribute to a larger range for which E/M exists ( Figure S1B). These results suggest that changes in the parameters which decrease the levels of miR-200 decrease the overall range of SNAIL levels for existence of E/M phenotype, or in other words, the model is sensitive towards low values of miR-200.
Similar sensitivity analysis was conducted for stemness regulatory circuit -LIN28/ let-7 in our previous manuscript [14], suggesting that changes in parameters which cause a decrease in the levels of let-7 result in the decrease of range of NF-kB levels for which the intermediate D/U state (medium LIN28, medium let-7) exists.
As the final step, we also varied the 'stemness window' parameters that are based on experimental data suggesting that the intermediate levels of OCT4 correspond to the maximum likelihood of being 'stem-like', and that both strong overexpression and strong down-regulation can force the cell to differentiate [16][17][18] [14]. OCT4 is an output of the coupled circuit, and translational activation of OCT4 by LIN28 is represented by a positive Hill function , where L I represents LIN28 levels, and L O c I 0 represents threshold levels of LIN28 for OCT4 activation = 480 * 10 3 molecules. The value of this Hill function lies between 0 and 1. Here, based on experimental data, we consider a 'guessed range' for the levels of OCT4 pertaining to 'stemness window' on this scale of (0-1), i.e. (0.25-0.65). If OCT4 values pertaining to any phenotype lies in this regime, that phenotype is within the 'stemness window' or has a heightened likelihood of gaining stemness. The LIN28 threshold levels (= 480 * 10 3 molecules) are chosen based on the overall scale of LIN28 values (~1000 * 10 3 molecules). Thus, the range of corresponding LIN28 levels for the 'stemness window' is [277 * 10 3 , 655 * 10 3 ] molecules. This range was varied by 20% to understand how the definition of the 'stemness window' affects our results, and we did not find any significant change in the phase diagrams obtained, for instance, results pertaining to inhibiting OVOL as shown in Figure 6, thereby suggesting the robustness of the model to our definition of 'stemness window'.
Collectively, the parameter sensitivity analysis suggests that as long as the levels of let-7 or miR-200 are not too low, our model is quite robust to parameter changes. This observation highlights an important facet -extremely low levels of let-7 and/or miR-200 will anyway render one or both coupling links irrelevant ( Figure S2) -inhibition from miR-200 to LIN28 and/ or inhibition from let-7 to ZEB would be more or less independent of α1 and/or α2 respectively. Therefore, in the parameter regime where coupling links can potentially largely affect the dynamics of the two circuits, our model is quite robust.

STABLE STATES OF THE COUPLED EMT AND STEMNESS CIRCUITS AT DIFFERENT SNAIL LEVELS
Inhibition of miR-200 on LIN28 and inhibition of let-7 on ZEB couple the two circuits -miR-200/ZEB and LIN28/let-7 together. We explore how different initial conditions, i.e. different levels of SNAIL at no coupling affect the stable steady states obtained at a given set of coupling parameters for the combined circuit (Figure 3-5).
Starting from monostable phase {E} as enabled by SNAIL levels= 180 * 10 3 molecules, α1 (strength of inhibition of LIN28 by miR-200) affects the coupled dynamics much more than α2 (strength of inhibition of ZEB by let-7) because levels of miR-200 are quite high initially and hence any changes in the set of stable states change depends much more on α1 than α2 ( Figure S2B). Conversely, starting from monostable phase {M} as enabled by SNAIL = 240 * 10 3 molecules, α1 (strength of inhibition of LIN28 by miR-200) affects the coupled dynamics of EMT and stemness much less than α2 (strength of inhibition of ZEB by let-7) because levels of miR-200 are quite low initially and therefore the set of stable states change depends much more largely on α2, specially at lower values of α1 ( Figure S2A).

STABLE STEADY STATES OF THE COUPLED CIRCUITS AT DIFFERENT VALUES OF (α1, α2)
In Figure 2, we illustrate that two phases that have the same number of stable states but are obtained

OVOL PULLS THE DYNAMIC 'STEMNESS WINDOW' TOWARDS E END OF EMT AXIS
In Figure 4, we analyze how OVOL affects the likelihood of the phenotypes E, M and E/M in gaining stemness, when cells are in monostable {M} phase enabled by the SNAIL levels = 340 * 10 3 molecules at no coupling between the circuits (α1 = α2 = 0). Here, we show the case when SNAIL levels = 360 * 10 3 molecules at no coupling (α1 = α2 = 0). Under very strong coupling, in the absence of OVOL, only the mesenchymal (M) phenotype can gain stemness ( Figure S4A), but in presence of OVOL, only hybrid (E/M) phenotype can gain stemness under very strong coupling ( Figure S4B). These results resonate with those presented in the main text that OVOL can aid the hybrid E/M phenotype to gain stemness and preclude M phenotype from gaining stemness.
We further explore the effect of OVOL in mediating the positioning of the 'stemness window' on the 'EMT axis'. Coupling OVOL to miR-200/ZEB enables a special feature not seen in the miR-200/ZEB circuit without OVOL -a monostable {E/M} phase [32]. Starting from this phase initially as enabled by SNAIL levels = 300 * 10 3 molecules, we add an external activation signal on OVOL to explore how the effect of over-expression of OVOL on EMT-stemness interplay. Overexpressing OVOL enables the E phenotype to gain stemness for almost the entire range of (α1, α2). Further, the hybrid E/M phenotype can also gain stemness at low α2 (weak inhibition of ZEB by let-7) ( Figure S5). This result is in stark contrast to the phase diagram obtained when OVOL was inhibited by an external signal ( Figure  6) where for low α2 (weak inhibition of ZEB by let-7), both E/M and M phenotypes can gain stemness, as well as the epithelial phenotype cannot gain stemness unless accompanied by the E/M phenotype for any values of (α1, α2).
Overall, OVOL can pull the 'stemness window' towards E end of the axis, and are consistent with previous experimental as well as theoretical observations [32,33] that the overexpression of OVOL induces a Mesenchymal to Epithelial Transition (MET), thereby enabling a much higher likelihood of epithelial phenotype to gain stemness. However, in the context of cancer, only rarely can cell lines be classified as purely epithelial; they are usually EMT-prone due to factors such as genomic plasticity, hence, OVOL is down-regulated during tumor progression as shown earlier [33]. Consequently, E phenotype might not be likely to gain stemness in the context of cancer (type III EMT) as in the case of wound healing and embryonic development (types I and II EMT). Figure S1: Parameter sensitivity analysis for the EMT regulatory circuit miR-200/ZEB driven by input signal SNAIL. Alphanumeric codes on x-axis represent the cases with different changed parameters by ± 20%. Number 1 after the alphanumeric code represents the increase of the parameter by 20% and number 2 after the alphanumeric code represents the decrease of the parameter by 20%. Case 0 represents the result in this work, which is the control case. A. represents the cases with ± 20% changes in the production and degradation rates. A1 and A2 represent the increase and decrease in the degradation rate of miR-200 (denoted by k μ 200 ) by 20% respectively. B1 and B2 represent the increase and decrease in the degradation rate of ZEB mRNA by 20% (denoted by k m z ) respectively. C1 and C2 represent the increase and decrease in the degradation rate of ZEB protein (denoted by k z ) by 20% respectively. D1 and D2 represent respectively the increase and decrease in the production rate of miR-200 (denoted by g μ 200 ) by 20%. E1 and E2 represent the increase and decrease in production rate of ZEB mRNA (denoted by g m z ) by 20% respectively. F1 and F2 represent the increase and decrease in the production rate of ZEB protein (g Z ) by 20% respectively. B. represents the cases with ± 20% changes in the thresholds and weight factors in the shifted Hill functions. G1 and G2 represent respectively the increase and decrease in the threshold levels for the inhibition of miR-200 by ZEB (denoted by Z 0 μ 200 ) by 20%. H1 and H2 represent respectively the increase and decrease in threshold levels of the shifted Hill function representing ZEB self-activation (denoted by Z 0 m Z ) by 20%. I1 and I2 represent the 20% increase and decrease in the threshold levels of SNAIL inhibition on miR-200 (denoted by S 0 μ 200 ) respectively. J1 and J2 represent the increase and decrease in the threshold levels of SNAIL activation on ZEB (denoted by S 0 m Z ) by 20%. K1 and K2 represent the increase and decrease in the threshold levels of miR-200 (denoted by μ 0 200 ) by 20%. L1 and L2 represent the increase and decrease in the weight factor of ZEB inhibition on miR-200 (denoted by λ Z, μ 200 ) by 20%. M1 and M2 represent the respective increase and decrease in the weight factor of SNAIL inhibition on miR-200 (denoted by λ S, μ 200 ) by 20%. N1 and N2 represent the respective increase and decrease in the weight factor of ZEB self-activation (denoted by λ Z, m Z ) by 20%. O1 and O2 denote respective increase and decrease in weight factor of SNAIL activation on ZEB (denoted by λ S, m Z ) by 20%. Dotted rectangles highlight the cases when the range of SNAIL values for which the cells can attain E/M phenotype vary the most, as compared to the control (case 0).