Modeling the Transitions between Collective and Solitary Migration Phenotypes in Cancer Metastasis

Cellular plasticity during cancer metastasis is a major clinical challenge. Two key cellular plasticity mechanisms —Epithelial-to-Mesenchymal Transition (EMT) and Mesenchymal-to-Amoeboid Transition (MAT) – have been carefully investigated individually, yet a comprehensive understanding of their interconnections remains elusive. Previously, we have modeled the dynamics of the core regulatory circuits for both EMT (miR-200/ZEB/miR-34/SNAIL) and MAT (Rac1/RhoA). We now extend our previous work to study the coupling between these two core circuits by considering the two microRNAs (miR-200 and miR-34) as external signals to the core MAT circuit. We show that this coupled circuit enables four different stable steady states (phenotypes) that correspond to hybrid epithelial/mesenchymal (E/M), mesenchymal (M), amoeboid (A) and hybrid amoeboid/mesenchymal (A/M) phenotypes. Our model recapitulates the metastasis-suppressing role of the microRNAs even in the presence of EMT-inducing signals like Hepatocyte Growth Factor (HGF). It also enables mapping the microRNA levels to the transitions among various cell migration phenotypes. Finally, it offers a mechanistic understanding for the observed phenotypic transitions among different cell migration phenotypes, specifically the Collective-to-Amoeboid Transition (CAT).

The effective (reduced) circuit, where μ 1 represents the inhibition on RhoA by either miR-34 or miR-200, and μ 2 represents similar inhibition on Rac1. A solid arrow denotes activation, and a solid bar denotes repression. A solid line represents transcriptional regulation, a dotted line represents indirect regulations on GTP loading or hydrolysis process via GEFs or GAPs, and a dashed line represents translational inhibition by microRNAs (c) Dynamical system characteristics of the effective circuit. The plot shows the nullclines and possible steady states corresponding to equation (1). When μ 1 = μ 2 = 100 molecules, the circuit can be tri-stable ( . Green solid circles denote the stable steady states, and green hollow circles denote the unstable steady states. Each stable state can be associated with a cell phenotype, as depicted by a cartoon beside them. circuits of EMT/MET and AMT/MAT. These two regulatory circuits are coupled through the microR-NAs, miR-200 and miR-34, that inhibit both RhoA as well as Rac1 [19][20][21][22][23][24] Fig. 1a. Therefore, the coupled circuit contains links that are both translational (through microRNAs) and post-translational (through GTPases) in nature. These interactions can potentially give rise to novel dynamical features that are not observed in circuits composed of only transcription factors, or those with both transcription factors and microRNAs 15,25 . This inherent complexity in the coupled system explains why the transitions between collective and solitary migration phenotypes have received limited theoretical attention.
Here, we follow and extend our previous theoretical framework to study small GTPase-based regulations in RhoA/Rac1 circuit by adding microRNA-mediated inhibition on the production of RhoA and Rac1. We show that this coupled circuit can introduce another stable steady state in addition to the three stable steady states (A, A/M, M) for the standalone Rac1/RhoA circuit 14 . We propose to associate this extra stable state with the hybrid epithelial/mesenchymal (E/M) phenotype, a phenotype that enables collective cell migration. We also found that the cells moving collectively (E/M) can switch to being in either an amoeboid (A) or a mesenchymal (M) phenotype directly, depending on the relative strength of inhibition of RhoA and Rac1 by the microRNAs. Our results further indicate that the high levels of microRNAs, miR-200 and miR-34, can restrict the transitions towards solitary migration phenotypes, even in the presence of EMT-stimulating signals, such as Hepatocyte Growth Factor (HGF). Our model is consistent with several experimental observations and provides the first step towards understanding how different levels of miR-200 and miR-34, and the Rho GTPases in a cell govern phenotypic plasticity during carcinoma metastasis under the influence of external signals in the tumor microenvironment.

Results
Mathematical Model of the Coupled Circuit. The EMT regulatory circuit regulates the Rac1/RhoA circuit via miR-34 and miR-200 Fig. 1a. Here, we use an effective (reduced) model ( Fig. 1b) to represent the association between these two circuits -miR-34/SNAIL/miR-200/ZEB (the core circuit for EMT/ MET) and RhoA/Rac1 (the core circuit for AMT/MAT) (Fig. 1a). As shown in Fig. 1b, the effect of EMT regulatory circuit is treated as two external signals to the Rac1/RhoA circuit, where (μ 1 ) represents the inhibition on RhoA by either of the two microRNAs (miR-34 or miR-200) while (μ 2 ) represents a similar inhibition on Rac1.
We have previously developed the theoretical framework for small GTPase-based regulatory circuit such as Rac1/RhoA circuit 14 , as well as that for microRNA-based regulatory circuit like miR-34/SNAIL and miR-200/ZEB 15 . Here, we use both these frameworks to study the dynamics of the coupled effective circuit. The deterministic dynamics of the circuit can be modeled by two rate equations as given below: where ⁎ R c represents the active Rac1 (Rac1-GTP), and ⁎ R h represents the active RhoA (RhoA-GTP). g R c is the production rate for Rac1, g R h and g R A h are basal and excitatory production rates for RhoA respectively.
are the corresponding degradation rates for Rac1 and RhoA. I R c and I R h represent two external signals that drive the GTP hydrolysis/loading process of the Rac1/RhoA circuit, such as Grb2 and Gab1 in c-MET signaling, whose effects have been analyzed in our previous work 14 . The Hill function ( ) represents the transcriptional self-activation of RhoA. B and J functions are defined to represent the total GTP loading and hydrolysis rates including the intrinsic ones for RhoA and Rac1 and also the activated ones by GEFs or GAPs (See Supplementary Information section 1) 14 . The inhibition effect of microRNAs is described by the function P([μ]) whose value ranges from 0 to 1 -0 denotes a strong silencing effect, and 1 denotes no effect 15 . This effective model (eq. (1)) was used for stability and bifurcation analyses while the detailed model shown in Supplementary Information section 1 was used for the dynamic simulations below. The model derivation and the parameter estimations are described in details in the Supplementary Information.
Now, considering the connections between Rac1/RhoA circuit and EMT regulatory circuit, i.e. in presence of low levels of microRNAs (miR-34 and miR-200), the system can display another stable state that has relatively low levels of both active RhoA and Rac1, as denoted by (LL, or (0, 0)) ( Fig. 1c). At certain levels of microRNAs, the system can also behave as a four-way switch such that all these steady states-LL, HL, HH, and LH -co-exist ( Supplementary Fig. S1). We hereby propose to associate this new (LL) state with the hybrid epithelial/mesenchymal (E/M) phenotype. This association is corroborated by experimental observations of the presence of balanced, moderate levels of active Rac1 and RhoA in cells of the hybrid E/M phenotype not only during cancer metastasis (type III EMT), but also during wound healing (type II EMT), and embryonic development (type I EMT) 26-28 . The RhoA/Rac1 Circuit Response to microRNA Signals. Next, as the first step towards analyzing the response of the Rac1/RhoA circuit to different levels of the external microRNA signals (μ 1 and μ 2 ), we consider μ 1 and μ 2 to be equal at all times; and that allows us to consider them together as a single external signal denoted by μ (Fig. 2). Different levels of the signal μ enable different sets of co-existing phenotypes and hence different multi-stable phases, each marked by a distinct color in the figure. At a high level of μ , cells have relatively lower activities of both RhoA and Rac1, and can thereby adopt only the LL state, i.e. a hybrid E/M phenotype (marked by the pink region). As the level of microR-NAs decreases, cells may undergo the transition to adopt a LH state (M phenotype) or a HL state (A phenotype) either deterministically or stochastically. With further decrease in the levels of μ , cells can attain relatively high levels of both active RhoA and active Rac1, thereby adopting the HH state (A/M phenotype) state. Our results are consistent with experiments showing that miR-34 and miR-200 can act as gatekeepers of metastasis and that their decrease leads to collective as well as individual migration 29,30 . Two-parameter Bifurcation Diagram. Next, we analyze the response of the circuit when the two microRNA signals (μ 1 and μ 2 ) can change independently (Fig. 3a). We considered the initial condition such that the cells can adopt any of the three individually migrating phenotypes -A, M, and A/M (marked by black dashed line in Fig. 2 where μ = 25,000 molecules). At different levels of these two Depending on how μ 1 and μ 2 change temporally, the cells follow different trajectories in this phase diagram and thus go through different phenotypic transitions as is illustrated in Fig. 3b-d. When both the microRNAs -one inhibiting RhoA and the one inhibiting Rac1 -are decreased at the same rate, we see the co-existence of all three solitary migration phenotypes-M, A/M and A phenotypes, i.e. plasticity of cell migration phenotypes is quite rich (Fig. 3b). However, if only one of the microRNAs decreases, cells finally adopt either A or M phenotype (Fig. 3c,d).
Temporal Dynamics. We demonstrate the temporal dynamics of a phenotypic transition of a cell from collective migration (E/M phenotype) to individual migration (A, M and A/M phenotypes). In the E/M phenotype (LL state), cells have high levels of microRNAs (μ 1 and μ 2 ). When these signals where the two bifurcation parameters are the two input signals μ 1 (y-axis) and μ 2 (x-axis). Each phase, marked by a distinct color (shown in bottom), corresponds to a different combination of coexisting phenotypes (Phase plane diagrams for each phase are showed in Supplementary Fig. S4). (b) One-parameter bifurcation diagram when both μ 1 and μ 2 change simultaneously at the same rate, allowing them to merge into one parameter (here it is referred to as μ ), as shown by the trajectory I in the two-parameter bifurcation phase diagram (Fig. 3a). (c) One-parameter bifurcation diagram for the circuit driven by varying levels of μ 1 for a fixed μ 2 = 4,000 molecules, corresponding to the trajectory II in the two-parameter bifurcation phase diagram (Fig. 3a). (d) One-parameter bifurcation diagram for the circuit driven by varying levels of μ 2 at a fixed μ 1 = 4,000 molecules, corresponding to the trajectory III in the two-parameter bifurcation phase diagram (Fig. 3a).
decrease, cells gradually gain the ability to migrate and can eventually start moving individually in one of the solitary migration phenotypes-A, M and A/M. Notably, in different biological contexts, the two microRNAs (the one inhibiting RhoA (μ 1 ) and the other inhibiting Rac1 (μ 2 )) may decrease at different rates, thus leading to a difference in the dynamics of phenotypic transition. As shown in Fig. 4, if the microRNA inhibiting RhoA (μ 1 ) decreases faster than the microRNA inhibiting Rac1 (μ 2 ) (Fig. 4a,b), the cells transit from the E/M phenotype to A phenotype, i.e. they undergo direct collective to amoeboid transition (CAT) as observed in fibrosarcoma cells 18 . Similar behavior is observed when both microRNAs decrease at same rate (Fig. 4c,d). In contrast, when (μ 2 ) decreases faster (Fig. 4e), E/M to M phenotype transition may occur (Fig. 4f), i.e. the cells undergo a complete EMT.
Effective Landscape. In order to understand the relative stability of the co-existing stable states (phenotypes), we construct the effective landscape by considering the biological noise that can be present due to fluctuations in gene expression. Noise can originate from many sources such as birth/death of species and binding/unbinding of proteins, and can induce switches among its stable steady states spontaneously without the action of any external signal. The effective potentials of the landscape are defined as the negative logarithm of the probability (P) of each state (− ln(P)) 31,32 . Therefore, the more common (or frequently adopted) steady state (or phenotype) would have a lower effective potential as shown by blue color. The larger the blue area surrounds a steady state, the more frequently the state or the phenotype is observed.
Here, we have simulated Langevin dynamics of our model by adding Gaussian white noise (see the method section for details) for three cases (Fig. 5a), each of which corresponds to a different level of μ . Without this inhibition signal (μ ), four steady states (HL, LH, HH, LL) can be detected. Notably, without the effect of noise, i.e. in the deterministic analysis, the system had three stable states (HL, LH, HH) ( Fig. 1), but noise can induce transitions to another steady stable state (LL). The states representing the A, E/M and M phenotypes are relatively more frequently observed than the state representing the state of the A/M phenotype, as indicated by the higher effective potentials for the HH basin (Fig. 5b). Also, we observed more transitions between the A and E/M or M and E/M, rather than the other possible When μ 1 decreases faster than μ 2 (a), the cell undergoes a transition from the LL to the HL state (b). When μ 1 and μ 2 decrease at the same rate (c), cell still switches from the LL to the HL state (d). When μ 1 decreases slower than μ 2 (e), the cell switches to the LH state from the LL state (f). Fig. S7a). However, with increasing levels of μ , the HH state disappears and the system could still stay in the HL, LH or LL states, but the probability of staying in the LL state increases while that for the LH and the HL states decreases ( Fig. 5c and Supplementary Fig. S7b). For even higher levels of μ , the transitions from the collective to the solitary migration is largely inhibited (Fig. 5d and Supplementary Information section 7) therefore indicating that microRNAs can potentially stabilize the E/M phenotype against the biological noise.

The switch responds to the input signals from c-Met pathway at different microRNA levels.
Until now, we have considered the effect of microRNAs on the AMT/MAT circuit in the absence of any EMT-inducing signal. Here, we analyze the response of the AMT/MAT switch or circuit in the presence of an EMT-inducing signal -the HGF/c-Met pathway. c-MET is a tyrosine kinase receptor that is encoded by an oncogene, and is activated when the signaling molecule HGF binds to it 33 . The downstream effectors of the c-MET/HGF pathway include Grb2 and Gab1 that can regulate cell migration phenotypes by activating RhoA or Rac1 34-38 respectively.
To consider the regulatory functions of all these signals combined together, we first calculated the two-parameter bifurcation diagrams (See Supplementary Information section 8) when the Rac1/RhoA regulatory circuit is driven by two of these signals -(a) μ and Grb2 ( Supplementary Fig. S8a), and (b) μ and Gab1 (Supplementary Fig. S8b). As expected, high levels of microRNAs suppress the stimulation of solitary migration by Grb2 and Gab1 signals and therefore retain the cells in the LL or hybrid E/M phenotype (Supplementary Fig. S8a). When the level of μ or microRNAs is further decreased, high Grb2 signal can induce the cell to undergo complete EMT to attain the M phenotype, while the regulatory function of Gab1 signal depends on the level of signal μ . Gab1 induces the cells to the M phenotype at the intermediate level of μ , but to A phenotype at the low level of μ ( Supplementary Fig. S8b). Therefore, the regulatory functions of Grb2 and Gab1 in inducing cell migration depend on the amount Next, we have calculated the two-parameter bifurcation diagrams of Rac1/RhoA circuit driven by Grb2 and Gab1 signals for various values of μ (Fig. 6a,b). For the high level of μ , the expression of Rac1 and RhoA is strongly inhibited, thus the plasticity of cell migration driven by Grb2 and Gab1 is limited only to {LL} and {LH} phases, indicating that only EMT/MET could occur (the 1st-2nd panels to the left, Fig. 6a). However, with a decrease in the levels of μ , cells can also adopt the HL state (A phenotype) (the 3rd panel to the left, Fig. 6a); and this phenotypic plasticity for cells only increases with the decrease in the signal μ . At extremely low levels of μ , the cells can adopt any of the three solitary migration phenotypes -A, M, or A/M -depending on the Gab1 and Grb2 levels (the 4th-5th panels to the left Fig. 6a). Again, these results are consistent with the experiments suggesting that the miR-200 and miR-34 inhibit cell migration and consequently metastasis, and their loss can lead the cells to display a spectrum of solitary migration phenotypes 29,30 . Also, we notice that this series of phase diagrams shows an asymmetry for the induction of different solitary migration phenotypes by Gab1 and Grb2 signals at different microRNA levels. Cells in the E/M phenotype may undergo a transition to the M phenotype (i.e. complete EMT) even at high levels of μ , but they can transit to the A or the A/M phenotypes only at low levels of μ . This asymmetry might underlie why different cell lines might prefer different migration phenotypes (or different distributions of migration phenotypes) during metastasis 13 .
Also, for a different parameter set where the Rac1/RhoA circuit is more responsive to Grb2 and Gab1 (the threshold for Grb2 and Gab1 signals was reduced (see supplementary information section 2)). cells can attain the {HL} phase for a larger range of parameters even in the presence of high level of μ , indicating a more metastatic behavior (the 1st-2nd panels to the left, Fig. 6b). When μ is reduced, the plasticity in cell migration is enhanced significantly (the 3rd-5th panels to the left, Fig. 6b). This difference in the response to the c-MET pathway can be attributed to different cell lines with distinct metastatic potential.
Phenotype Distribution. Even in the same cell line, cells often have non-heritable phenotypic variability or heterogeneity 39 . Such variability can often account for different model parameters for each cell. To capture this heterogeneity, we calculated the population distribution of the phenotypes for the cells with the same levels of microRNAs, but some variability in other model parameters. More specifically, we extended our simulations to a population of 5,000 cells. Each cell has different circuit parameters, which can be ± 5% different from the original parameters (details in Method section). In Fig. 7, we show the percentages of cells in one of the four different possible phenotypes-(HL, LL, LH, and HH)-for different levels of the signal μ . We found that for a high level of signal μ , a significant percentage of cells adopt the E/M phenotype (Fig. 7a,b), thereby indicating that the cells with high levels of microRNA are more robust to the parameter variability against maintaining the E/M phenotype.

Discussion
Cellular plasticity mediated by the EMT/MET and the AMT/MAT regulatory circuits has been studied extensively individually 2,5,14,15 . However, understanding the regulation of these transitions together has not received enough attention. Here, we have presented, to the best of our knowledge, the first theoretical attempt to understand the interplay between EMT/MET and AMT/MAT by investigating the dynamics of the RhoA/Rac1 circuit (the core circuit for AMT/MAT) driven by miR-200 and miR-34, the gatekeepers of epithelial phenotype and inhibitors of metastasis. Our results explain how the cancer cells can attain different phenotypes and transitions among them during metastasis-collective migration as a cluster of Circulating Tumor Cells (CTCs) 40 , and different modes of individual migration-amoeboid (A), mesenchymal (M), or hybrid A/M phenotype 5,7,[11][12][13] . We further present several predictions about how microRNAs and stromal signals such as HGF affect these transitions.
Our previous work showed that the RhoA/Rac1 circuit could perform as a three-way switch allowing three stable steady states -(high RhoA-GTP, low Rac1-GTP), (low RhoA-GTP, high Rac1-GTP), and (high RhoA-GTP, high Rac1-GTP). These three states can be associated with the A, the M and the hybrid A/M phenotypes respectively 14 , as supported by several experimental evidences. For example, amoeboid cells have been shown to have high actomyosin contractility due to high levels of RhoA-GTP 41 ; mesenchymal cells have high actin polymerization due to high levels of Rac1-GTP 41 ; the hybrid A/M phenotype displays mixed amoeboid and mesenchymal characteristics, suggesting a comparable level of both actin polymerization and actomyosin contractility due to high levels of both Rac1-GTP and RhoA-GTP 7,11-13 . Notably, similar 'hybrid' manifestations of mixed amoeboid and mesenchymal traits may not be unique to cancer, but can also be observed in other instances, such as neutrophils and leukocytes [42][43][44] .
When incorporating the microRNA signaling on the RhoA/Rac1 circuit, we observed an additional stable state (LL), when cells of which have relatively low levels of both RhoA-GTP and Rac1-GTP as compared to the other states (LH, HH, and HL). We propose to associate this new state with the hybrid E/M phenotype. Although there is no direct quantitative measurement yet of the active levels of RhoA-GTP and Rac1-GTP in the hybrid E/M cells, our association of the LL state with hybrid E/M phenotype is consistent with the following experiments. First, moderate levels of both active Rac1 and RhoA have been reported to promote wound healing 27 , a typical case of the collective migration of the hybrid E/M cells 2 . Second, a moderate level of active RhoA has been shown to induce the hybrid E/M phenotype for the human colon adenocarcinoma cell 26 , while both dominant-negative and constitutively active Rac1 has been reported to damage the collective migration of boarder cells during Drosophila early development 28 . Both evidences support the relatively low activity of Rac1 and RhoA in collectively migrating hybrid E/M cells as compared with the solitarily migrating cells in A, M and A/M phenotypes.
One limitation of our current model is the exclusion of the epithelial (E) phenotype from our analysis because this model is too simple to explain the elusive dynamics of active RhoA and Rac1 in epithelial tissue establishment and maintenance. Epithelial cells are reported to have higher levels of microRNAs (miR-200 and miR-34) than in the hybrid E/M cells 15,45 . Therefore, according to our model, these cells are likely to have even lower levels of active RhoA and Rac1 than those in hybrid E/M cells. On the other hand, experiments suggest that the active levels of RhoA and Rac1 vary significantly in different stages of epithelial establishment 46 . For instance, at the initial stage of epithelialization, de novo cell-cell adhesion and adherens junctions expansion require the activity of Rac1 and RhoA respectively. But when EMT is induced, levels of active Rac1 decrease at first and then increase again 47 . Therefore, the dynamics of Rac1/RhoA activity need to be fine-tuned either spatially and/or temporally in order to maintain the E-cadherin-mediated cell-cell adhesion. But we currently have not yet modeled the spatial heterogeneity of the two GTPases, RhoA and Rac1, as needed to allow investigation of the role of spatiotemporal dynamics of these proteins in determining the changes in cell shape during epithelialization and transitions between different phenotypes. Besides, to understand the complex temporal dynamics of the activity of Rac1 and RhoA, additional models should incorporate the interplay among the key proteins in the regulatory circuit such as E-cadherin, N-cadherin, and the Rho GTPases.
From the perspective of dynamical systems, we found that the Rac1/RhoA regulatory circuit can behave as a multi-stable switch in the regulation of different phenotypes. Such multi-stability is a hallmark of self-activating toggle switches (SATS), i.e. mutually inhibitory feedback loop with self-activations on both elements of the loop 48 as seen in RhoA/Rac1 circuit (Fig. 1b). Similar multi-stability is also seen in circuits governing Cancer Stem Cells (CSCs) 49,50 and in circuits governing cell-cell communication 51,52 , and can often be used by the cancer cells to adapt to their rapidly changing microenvironments during metastasis 4 . Restricting such cellular plasticity as enabled by these multistable systems can inhibit possible co-operation between different subpopulations and hence hamper tumor metastasis [53][54][55][56][57] .
To conclude, we present the first tractable framework towards understanding the transitions among the collective migration phenotype and the solitary migration phenotypes. We found that the microR-NAs, miR-34 and miR-200, govern various phenotypic transitions and can also mediate the effect of other signaling pathways such as Grb2/Gab1 signals on such transitions. Our framework can be further extended by incorporating many other extracellular signals to explore their impact on cellular plasticity and migration, and can possibly provide important insights into cancer therapies that target cell migration during metastasis.

Method
Deterministic analysis. The dynamics of the coupled circuit is described by a set of six nonlinear chemical rate equations as shown in equation (S3). However, this model can be further simplified into two rate equations shown in equation 1 (Details in Supplementary Information section 1). Thereafter, we performed stability analysis on the effective model. To compute all the possible steady state solutions, we calculated two nullclines -the first one satisfies Fig. 1c), and the second one satisfies Fig. 1c). The nullclines are constructed on the phase plane by the concentrations of RhoA-GTP (x-axis) and Rac1-GTP (y-axis). The intersections of these two lines are the steady states for the whole system. The stability of the steady states can be further determined by the linear approximation method 58 . The nullclines are computed by the contour-based method (contour() function in Matlab), and the bifurcation diagrams are calculated by the MatCont package 59 . Besides, we simulated the temporal dynamics of the circuit (described by equation (S3)) by using the ode15s solver in Matlab.
Similar to our previous work, we modeled the Grb2 and Gab1 signals I R c and I R h as follows: process to attain equilibrium rapidly, thus equation (S3) can be reduced to four rate equations with the rate equations for Rac1-GDI and Rho-GDI set to 0. In this study, we assumed the dynamics of the circuit to be influenced by a Gaussian white noise. So the dynamics can be described by a Langenvin equation h , ϒ is a constant representing the noise level. The stochastic differential equation is integrated by the Euler-Maruyama method 60 . The effective landscape for the state X is defined by E = − ln(P(x)) 31,32 , where P(x) is the probability for the system in the state x. The Langevin simulation gives a long time trajectory (10 6 hour) of the protein expressions of the circuit, from which we computed the probability of being in any of the steady states, and the rates of the transitions among these states (details in Supplementary Information  section 7).
Phenotype distribution. To mimic the heterogeneity of a cell population, we randomly generated 5,000 sets of circuit parameters, in which each parameter was randomly perturbed away from the original value by at most 5% (random sampling follows a uniform distribution). For each of the Hill coefficients, it follows a discrete uniform distribution within [n − 1, n + 1], where n is the original value. We evaluated the dynamics of the circuit for each cell by modeling the circuit with each randomized set of parameters. We performed stability analysis, and calculated all stable steady states, each of which was mapped to the corresponding phenotypes based on the levels of the proteins. The mapped phenotypes were weighted by the number of stable steady states for each cell to calculate the fraction of cells in each phenotype.
The relevant Matlab and C codes are accessible at https://rice.app.box.com/s/ 7vfm90cwh4lstnpsa72xhlairugdj6hh.