Modeling of mesenchymal hybrid epithelial state and phenotypic transitions in EMT and MET processes of cancer cells

Based on the transcriptional regulatory mechanisms between microRNA-200 and transcription factor ZEB in an individual cancer cell, a minimal dynamic model is proposed to study the epithelial-mesenchymal transition (EMT) and mesenchymal-epithelial transition (MET) processes of cancer cells. It is shown that each cancer cell can exit in any of three phenotypic states: the epithelial (E) state, the mesenchymal (M) state, and the epithelial/mesenchymal (E/M) hybrid state, and the state of cancer cell can interconvert between different states. The phase diagram shows that there are monostable, bistable, and tristable phenotypic states regions in a parameters plane. It is found that different pathway in the phase diagram can correspond to the EMT or the MET process of cancer cells, and there are two possible EMT processes. It is important that the experimental phenomenon of E/M hybrid state appearing in the EMT process but rather in the MET process can be understood through different pathways in the phase diagram. Our numerical simulations show that the effects of noise are opposite to these of time delay on the expression of transcription factor ZEB, and there is competition between noise and time delay in phenotypic transitions process of cancer cells.

understand the existence of E/M hybrid phenotypic state of cancer cell, a core regulatory network governing epithelial-mesenchymal plasticity had been proposed by Lu et al. 19 , the core regulatory network is composed of four components (two families of E-box binding transcription factors SNAIL and ZEB, and two families of microRNAs (miRs) miR-34 and miR-200). There are two double negative feedback loop between TFs and miRs in the regulatory network, that is, the miR-34/SNAIL and the miR-200/ZEB mutually inhibiting loops as shown in Fig. 1(a), and it was found that the miR-200/ZEB module functions as a ternary switch, allowing not only for the epithelial and mesenchymal phenotypes but also for a hybrid phenotype 10,19 .
On the one hand, a large number of experiments confirmed that non-coding RNAs can regulate EMT, in particular, the expression of miR-200 family is strongly associated with epithelial differentiation, and a reciprocal feedback loop between the miR-200 family and the ZEB family of transcription factors tightly controls both EMT and MET [20][21][22] . Moreover, the EMT and MET mediated by miR-200/ZEB module was found in human various tumor cells, for example, in the colon cancer cells 23,24 , in the breast cancer cells 25,26 , etc. Although there are two double negative feedback loops between TFs and miRs in the core regulatory network 19 , that is, the miR-34/ SNAIL and the miR-200/ZEB mutually inhibiting loops as shown in Fig. 1(a) with four variables, it was also showed that the miR-34/SNAIL module is a monostable circuit functions as a noise-buffering integrator, and the miR-200/ZEB module functions as a ternary switch (i.e., miR-200/ZEB acts as the decision making module for cancer cells to undergo partial or complete EMT), allowing not only for the epithelial and mesenchymal phenotypes but also for a hybrid phenotype, and these theory-based findings are consistent with several experimental findings 10,19 . Therefore, a minimal dynamic model of miR-200/ZEB module could be sufficient and appropriate to simulate EMT/MET.
On the other hand, both fluctuation and time delay are ubiquitous in various intracellular signal networks [27][28][29][30] , where the noise could originate from the random births and deaths of individual molecules or the fluctuations in biochemical reactions, and the time delay could arise from finite propagation speeds of signal molecules. The effects of both noise and time delay have been widely studied in various biological systems, such as the genetic regulatory networks 31,32 , the neuronal systems 33-36 , etc. Interesting questions now arise: is there a minimal dynamic model to understand the E/M hybrid state in the EMT process? What are the effects of noise and time delay on the minimal regulatory mechanism of multiple phenotypic states for the cancer cell?
In this paper, based on the core regulatory network proposed by Lu et al. 19 , a general minimal dynamic model with two variables is proposed to understand the EMT and MET processes and the effects of both noise and time delay. By utilizing the general minimal mathematical model, it can be found that each cancer cell may exit in any of three phenotypic states (E, M, and E/M hybrid states), there are two possible EMT pathways, and the phenomenon of E/M hybrid state in the EMT processes but rather in the MET process observed by experiments 11,12 can be revealed through different pathways in a parameters plane (phase diagram). In addition, the effects of noise and time delay on the interconversions between phenotypic states are discussed by using numerical simulations.

Model and Method
The four components model of the miR-34/SNAIL and the miR-200/ZEB mutually inhibiting loops 19,37 of cancer cells is shown in Fig. 1(a), it can be found that the interactions between transcription factors and microRNAs have symmetrical properties. The expression of SNAIL from the miR-34/SNAIL module can been considered as a external signal action on the miR-200/ZEB module, on the contrary, the expression of ZEB from the miR-200/ ZEB module can be considered as a external signal action on the miR-34/SNAIL module. Therefore, one can simplify core regulatory network by a general minimal regulatory model with two variables as shown by Fig. 1(b), in this paper, X 1 represents the transcription factor ZEB family, and X 2 represents the miR-200, and the SNAIL signal activates the expression of ZEB and inhibits the expression of miR-200 simultaneously.
In the deterministic description, the dynamic model of the minimal regulatory motif Fig. 1(b) can be described by following ordinary differential equations in dimensionless:   where x 1 and x 2 are the expression levels of transcription factors ZEB family and miR-200. a is the activation strength of ZEB induced by both SNAIL signal and itself, b is the inhibited strength of miR-200 induced by both SNAIL signal and X 1 , c is the inhibited strength of ZEB by miR-200, and k 1 and k 2 are the self-degradation rates of ZEB and miR-200, respectively. Parameters n 1 , n 2 , and n 3 are the Hill coefficients which control the steepness of the sigmoidal function, and θ is a threshold. Here, the parameter values are n 1 = 2, n 2 = 6, and n 3 = 3 which are on the basis of the different binding sites obtained from experiments 19 , and k 1 = k 2 = 1.0 and θ = 0.5 38 . For simplicity, the inhibited strength of ZEB by mir-200 is set as constant value c = 1.0. Here, it was hypothesized that the negative regulations on X 2 by both SNAIL and X 1 (and the positive regulations on X 1 by both SNAIL and itself) can simultaneously occur, thus, the regulations of first item in the right side of Eqs (1 and 2) follow an "and" rather than "or" logic. In the stochastic description, above general dynamic model Eqs (1,2) are written by following Langevin equations where ξ t ( ) 1 and ξ t ( ) 2 are the Gaussion white noises with zero mean and auto-correlation functions . We consider a homogeneous situation (i.e., D 1 = D 2 = D), which means that the noises in Eqs (3 and 4) represent the total effect of intrinsic and extrinsic noises. The probability distribution P x x t ( , , ) 1 2 of system (3) and (4) obeys the Fokker-Planck equation 38,39 : The stationary probability P x x ( , ) st 1 2 obtained from the Fokker-Planck Eq. (5) can be used to indicate the phenotypes proportion distribution of cancer cells. In the equilibrium case, the potential function U x x ( , ) 1 2 for non-equilibrium system can be defined by the stationary probability: Each minimum of U(x 1 ; x 2 ) corresponds to one state (or phenotype) of a cancer cell, and the phenotypic switching of cancer cells is understood that the state of the cancer cell moves from one minimum of potential landscape to another.
To quantify the properties of phenotypic switching between states in the case of multiple phenotypic states coexistence, one can calculate the escape time from one steady state of U(x 1 , x 2 ) to another. A rigorous definition of escape time is provided by the mean first-passage time (MFPT) of a stochastic process y(t) to reach the point y 2 with initial condition y(t = 0) = y 1 40,41 y y st y st 1 2 In order to study the time delay arise from finite propagation speeds of regulatory molecule, here we assume that the activation of ZEB induced by itself through the post-transcriptional regulation has a time delay τ, thus, Eq. (3) is rewritten by:

Results
Although the expression levels of genes are not the only factor, they are the major influencing factor to determine the phenotypes of a cancer cell. According to the expression levels of X 1 and X 2 in the general minimal dynamic model Eqs (1, 2), one can assume that 0 represents the relatively low expression level, 1 represents the relatively high expression level, and 1/2 represents the relatively intermediate expression level. There are three states (x 1 ,    Fig. 2) was chosen to analyse the stability of deterministic trajectories of dynamic system. Figure 3 shows there is one stability point which correspond to the single E state (Fig. 3(1)) and the single M state (Fig. 3(2)) in the monostable regions. There are two stability points which correspond to the coexistence of (E, M) states ( Fig. 3(4)), and the coexistence of (E, E/M) states ( Fig. 3(5)) in the bistable region. There are three stability points which correspond the coexistence of (E, E/M, M) states ( Fig. 3(3)) in the tristable region.
The coexistence of multiple phenotypic states and the phenotypes transitions of cancer cells had been confirmed by some experiments in the EMT process [43][44][45] . It can be found from the phase diagram Fig. 2, that the epithelial (E) phenotype is appeared in all multiple phenotypic states regions, which implies that epithelial state is a fundamental phenotypic state and more robust than other states 46,47 . EMT and MET processes of cancer cells. It was demonstrated by experiments 11,12 that the E/M hybrid phenotypic state of cancer cells can be observed in the EMT process, but can not be observed in the reverse MET process. On the one hand, from the core regulatory network (see Fig. 1(a)) governing epithelial-mesenchymal The EMT process. In order to illustrate the EMT process, for example, the pink dotted line (with b = −10a + 20) in Fig. 2 represents a possible pathway that a cancer cell undergoes a EMT process with the increasing of activation strength a of ZEB. Figure 4(a) shows that, with the increasing of activation strength a of ZEB which can be induced by external signal (such as the SNAIL), the cancer cell undergoes a succession of transitions: the single E phenotypic state → the coexistence of (E, M) states → the coexistence of (E, E/M, M) states → the coexistence of (E, M) states → the single M phenotypic state. Figure 4(b) shows the stationary probability P st (x 1 , x 2 ) at the three phenotypic states of cancer cell in different region of Fig. 4(a). With the variations of parameters a and b, although the cancer cell experiences twice coexistence of (E, M) states in the EMT possess, the probability of each state is very different by comparing Fig. 4b(II) with Fig. 4b(IV). The difference between above two EMT paths is the second phase: one is the coexistence of (E, M) states, and the other is the coexistence of (E, E/M) states. However, in both EMT processes, it is interesting that the E/M hybrid phenotypic state can be observed in the EMT processes. It is also found that the E state always exist until the the activation strength a of ZEB is too large, which suggests that the E state may be a fundamental phenotype for the cancer cell which is consistent with the experimental observations 46,47 .
The MET process. To illustrate the MET process, we considered the blue dotted line (with b = −13.33a + 10) in Fig. 2, which represents a possible pathway that a cancer cell undergoes a MET process with the decreasing of activation strength a of ZEB. Figure 6(a) shows that, with the decreasing of activation strength a of ZEB, the cancer cell will undergo a succession of transitions: M → (E, M) → E. Figure 6(b) shows the stationary probability P st (x 1 , x 2 ) at the two phenotypic states of cancer cell in different region of Fig. 6(a). It can be seen that the probability of M state is decreased, however, the probability of E state is increased with the variations of parameters a and b. In above MET process, it is interesting that the cancer cell can not exist the E/M hybrid phenotypic state, that is, there is no the E/M hybrid state in the MET process.
By utilizing our general minimal mathematical model with two variables, the different phenotypic states of cancer cell as shown in ref. 19 can be reproduced. Most importantly, the phenomenon of E/M hybrid state in the EMT processes but rather in the MET process observed by experiments 11,12 can be understood through different pathways in the parameters (a, b) plane.

Phenotypic transitions in tristable region.
The regulation of cancer cell state decisions is of particular significance in tumor pathobiology. The cancer cells frequently exist in any of possible phenotypic states, and cancer cells can interconvert between different states. There are extensive studies 2-4 on the phenotypes E and M of cancer cells. However, the E/M hybrid state of cancer cells is more complex and also important in physiology and pathology 10 since the hybrid phenotype with the mixed traits can lead to cancer cells collective migration. The hybrid phenotype can induce the cancer cells maximum cellular plasticity 6 .
In  With the increasing of noise intensity, Fig. 7(a) shows that the stationary probability of M state is gradually increased, and the potential landscapes corresponding to Fig. 7(a) are shown in Fig. 7(b). It is found that the cancer cell has a large probability of transition into the M state through the E/M hybrid phenotype. When the cancer cell is in the M state, it becomes more migratory and invasive 4,6 . To quantify the properties of phenotypic switching between different states, one can calculate the mean first passage time (MFPT) from one steady state to another according to Eq (7). In addition, the barrier height of minima of potential function Eq (6) can also be used to imply the properties of phenotypic switching, the height of barriers of two attractors (e.g., the A and B states) is defined by: , where S is the saddle point between A and B states. With the increase of the noise intensity, Fig. 7(c) shows that the barrier heights of E, E/M hybrid, and M states are decreased, which means the larger noise intensity can lead to the steady states more instability, and the noise lets cancer cell easily escapes from the steady states. Figure 7(d) shows the MFPT from E state to M state or from M state to E state, respectively. Although all MFPTs are decreased with the increasing of noise intensity, the MFPT from E to M state is smaller than that from M to E state, which means that the M phenotype is much more stable than the E phenotype for cancer cell.
Interestingly, with the increase of noise intensity, it is found that the barrier height ϕ SE/M of E/M hybrid phenotype (see the blue line in Fig. 7(c)) is decreased, and then reaches a platform value when the noise intensity is larger than a critical value D c = 0.03, that is, the E/M hybrid state always exists when D > D c . The cancer cell in the E/M hybrid state can acquire the capabilities of resistance and robustness to against the high pressure environment. It was reported that cancer cells in the E/M hybrid phenotype exhibit stemness characters during EMT in developmental, regenerative, as well as pathological contexts [48][49][50][51] , and related to drug resistance 52-54 . Effects of time delay and noise on phenotypic transitions. The multi-phenotypic states of cancer cells arise in a cancer cell with different genes expression levels, for example, the low ZEB expression corresponds to E phenotype, the high ZEB expression corresponds to M phenotype, and the intermediate ZEB expression corresponds to E/M hybrid phenotype 19 . In order to discuss the effects of time delay and noise on the phenotypic states transitions, two points in the phase diagram Fig. 2   In the bistable region. In the absence of noise, when the cancer cell starts in the M state or in the E state (i.e., the steady state of expression level of x 1 over the time, not the initial value of x 1 ). With the increasing of time delays, Fig. 8(a,b) (the top) show that the expression level x 1 of ZEB is independent of the time delay τ, only depends on the initial state of cancer cell.
In the presence of noise, however, Fig. 8(a,b) (the bottom) show that the expression level x 1 of ZEB is not only depend on the time delay but also on the initial state of cancer cell. When the cancer cell starts in the M state ( Fig. 8(a)), the cancer cell transforms from M to E state (i.e., M → E) with the increasing of time delay. When the cancer cell starts in the E state ( Fig. 8(b)), with the increasing of time delay, the cancer cell firstly transforms from E to M state, and then to E state (E → M → E).
In the tristable region. In the absence of noise, when the cancer cell starts in the M state, Fig. 9(a) (the top) shows that, with the increasing of time delay, the expression level x 1 is firstly oscillation, and then reaches to a middle constant level, and finally to low constant level. It implies that the cancer cell undergoes a succession of transitions: from M to E/M hybrid, and then to E state (M → E/M → E). When the cancer cell starts in the E/M hybrid state, Fig. 9(b) (the top) shows that, with the increasing of time delay, the expression level x 1 is firstly oscillation, and then reaches to a low constant level. It implies that the cancer cell undergoes a transition from E/M hybrid to E state (E/M → E). When the cancer cell starts in the E state, Fig. 9(c) (the top) shows that, with the increasing of time delay, the expression level x 1 is firstly oscillation, and then reaches to a low constant level, it implies that the cancer cell stays in the E state at all time.
In the presence of noise, with the increasing of time delay, when the cancer cell starts in the M state, Fig. 9

Conclusions and Discussions
In this paper, based on the core regulatory network proposed by Lu et al. 19 , a general minimal dynamic model with two variables is proposed to understand the EMT and MET processes and the effects of both noise and time delay on phenotypic transitions of cancer cell. By utilizing the stability analysis of dynamic model, it is shown that single cancer cell exits three phenotypic states: the epithelial state, the mesenchymal state, and the epithelial hybrid mesenchymal state. The phase diagram of dynamic model shows that there are monostable, bistable, and tristable states regions in a parameters plane, that is, the cancer cell exists the coexistence of multiple phenotypic states, and interconverts between different states.
By using the bifurcation analysis and the stationary probability of genes expression level, it is found that different pathways in the phase diagram correspond to the EMT or the MET process of cancer cells, and there are two EMT paths. The phenomenon of E/M hybrid state in the EMT processes but rather in the MET process observed by experiments 11,12 can be understood through different pathways in the (a, b) parameters plane. In the multiple states (E, E/M, M) coexistence region, it is found that the transitions in phenotypes of cancer cells can be induced by noise. With the increasing of noise intensity, the mean first passage time and the height of barriers of two attractors show that the cancer cell has a large probability of transition into the M state through the E/M hybrid phenotype. It is found that the barrier height of E/M hybrid phenotype is firstly decreased, and then reaches a platform value when the noise intensity is larger than D c = 0.03, that is, when D > D c , the E/M hybrid state always exists no matter how large the noise intensity is. The cancer cell in the hybrid state can acquire the capabilities of resistance and robustness to against the high pressure environment. It was reported that cancer cells in the E/M hybrid phenotype exhibit stemness characters during EMT in developmental, regenerative, as well as pathological contexts [48][49][50][51] , and related to drug resistance [52][53][54] .
The numerical simulations show that the effects of noise on multiple phenotypic states transitions are opposite to these of time delay, and there is competitive mechanism between noise and time delay in the multiple phenotypic transitions process. In conclusion, a general minimal dynamic model with two variables is proposed to understand the EMT and MET processes and the effects of both noise and time delay, and the phenomenon of E/M hybrid state in the EMT process but rather in the reverse MET process observed by experiments is revealed through different pathways in a parameters plane (phase diagram).
Epithelial to mesenchymal transition (EMT) plays an important role in embryonic development, tissue regeneration, and cancer metastasis. By using our minimal model, it was found that there are three states of each cancer cell, the states could be reproduced by using the model of miR-200/ZEB module. The three states are (high miR-200/low ZEB) or (1, 0), (low miR-200/high ZEB) or (0, 1), and (medium miR-200/medium ZEB) or (1/2, 1/2), which correspond to the epithelial phenotype (E), mesenchymal phenotype (M) and epithelial/mesenchymal hybrid (E/M), respectively.
Compared our results with previous work 19 , the phases and bifurcation diagrams looked differences, the reasons for the differences are (i) the parameters in our general minimal model are dimensionless, and the parameters in the model of core regulatory network 19 are dimension. (ii) The assumptions of genes switching functions (or Hill functions) in both models are different. However, the biological significance of phases and bifurcation diagrams are the same. For example, by using our minimal dynamic model, the phase diagram of phenotypic states of single cancer cell can be drawn in parameters (a, b) plane as shown by Fig. 2, the phase diagram can reproduce the coexistence of multiple phenotypic states and the phenotypes transitions of cancer cells: two monostable regions which correspond to the E sate or the M state, two bistable regions which correspond to the coexistence of E and M states or the coexistence of E and E/M hybrid states, and a tristable region which corresponds to the coexistence of E, E/M hybrid, and M states. Importantly, by using our general minimal model, it is found that the E state always exist until the the activation strength a of ZEB is too large in the bifurcation diagrams (see Figs 4(a)-6(a)). Our results suggest that the E state may be a fundamental phenotype for the cancer cell which is consistent with the experimental observations 46,47 . Our general minimal model is involved in negative regulation between microRNA and transcription factor such as miR-200/ZEB. It is well known that RNA interference and microRNA are new technologies in drug development. For example, Olmeda et al. 55 disclosed that specific silencers of endogenous miRNAs, antagomirs, are powerful tools to silence specific miRNAs in vivo. Iorio and Croce 56 reviewed the involvement of microRNAs in cancer and their potential as diagnostic, prognostic, and therapeutic tools. Wang et al. 57 discussed the CSCs, EMT and the role of regulation by miRNAs in the context of drug resistance, tumor recurrence and metastasis. Therefore, microRNAs associated with EMT such as miR-200 family could be exploited as therapeutic strategies in the future. Therefore, our general minimal model for the miR-200/ZEB module could provide new insights into the roles of microscopic regulatory mechanisms, noise, and time delay in single cancer cell, and the dynamic model might give some insights for various tumors clinical therapy strategies.