CONTROLLING STOCHASTICITY IN EPITHELIAL-MESENCHYMAL TRANSITION THROUGH MULTIPLE INTERMEDIATE CELLULAR STATES

. Epithelial-mesenchymal transition (EMT) is an instance of cellular plasticity that plays critical roles in development, regeneration and cancer progression. Recent studies indicate that the transition between epithelial and mesenchymal states is a multi-step and reversible process in which several in- termediate phenotypes might coexist. These intermediate states correspond to various forms of stem-like cells in the EMT system, but the function of the multi-step transition or the multiple stem cell phenotypes is unclear. Here, we use mathematical models to show that multiple intermediate phenotypes in the EMT system help to attenuate the overall ﬂuctuations of the cell population in terms of phenotypic compositions, thereby stabilizing a heterogeneous cell population in the EMT spectrum. We found that the ability of the system to attenuate noise on the intermediate states depends on the number of inter- mediate states, indicating the stem-cell population is more stable when it has more sub-states. Our study reveals a novel advantage of multiple intermediate EMT phenotypes in terms of systems design, and it sheds light on the general design principle of heterogeneous stem cell population.


1.
Introduction. Epithelial-mesenchymal transition (EMT) is an extreme form of cellular plasticity that is involved in morphogenesis, tissue regeneration and cancer progression. During EMT, epithelial cells undergo dramatic changes in cell morphology and behavior to form mesenchymal cells. These changes include loss of cell-to-cell junction, loss of cell polarity and acquisition of migratory and invasive properties [5] [16]. The migratory behavior of newly formed mesenchymal cells is critical for the formation of internal organs during embryonic development [5], and it is also involved in cancer metastasis, which often requires the dissemination of cancerous epithelial cells from the primary tumors [25]. After arriving to their destinations, mesenchymal cells sometimes revert to epithelial cells via a process called mesenchymal-epithelial transition (MET), suggesting that EMT is a reversible process [2].
Previous theories and experiments showed that EMT does not give rise to terminal mesenchymal phenotype in some biological scenarios [24] [30] [36]. In other words, intermediate epithelial-mesenchymal phenotype may be generated by partial EMT. Interestingly, it has been shown that partial EMT, or intermediate EMT phenotype, is associated with the stemness of epithelial and mesenchymal cells [9].
In addition, it was hypothesized that partial EMT is responsible for collective migration and invasiveness of cancer cells [24]. These observations and hypothesis strongly suggest that EMT is a multi-step transition, and various degrees of the transition may have distinct physiological outcomes.
Using modeling and experimental approaches, we recently demonstrated the existence of multiple intermediate phenotypes in the EMT spectrum [12]. This is consistent with the observations which show four epigenetic states in multi-step EMT [29], multiple epithelial/mesenchymal stem cell phenotypes [1] and four distinct types of ovarian cancer cells in the EMT spectrum [14]. We also showed that the multiple intermediate states arise from complex gene regulatory networks consisting of interconnected feedback loops [12]. Previous theoretical studies suggested the role of heterogeneity in cell population and its feedback control in stem cell regeneration [20]. However, it is unclear why the epithelial/mesenchymal cell population needs these multiple intermediate states, or why the stem cells in the EMT spectrum need to have several subtypes. In particular, what performance objective does the system evolve to achieve by having multiple intermediate states at the expense of simplicity of cellular phenotypes?
To address these questions, we built and compared a series of mathematical models containing various number of intermediate states in EMT spectrum. By analyzing the behaviors of the system in facing noise in cell populations, we found that multiple intermediate phenotypes in the system help to attenuate the overall fluctuations of the cell population in terms of phenotypic compositions, suggesting their ability to stabilize a heterogeneous cell population. When we compared the models with at least one intermediate state, we found that the number of intermediate states positively correlates with the ability to attenuate noise on the intermediate states, indicating that the stem-cell population is more stable when it has more sub-states. These results suggest a performance objective that the EMT system might have evolved to achieve by having multiple intermediate phenotypes, and it sheds light on the general principles of heterogeneous stem cell population in term of systems design.
2. Mathematical models and stochastic simulation of multiple state EMT models.
2.1. Model construction. We explored the functional significance of the intermediate states in the EMT process by adapting the population dynamics to modelling the transition from epithelial to mesenchymal phenotypes. Each steady state in the transition was portrayed by a population of cells that assumed the phenotype reflective of that state. In other words, an n-state EMT system has n different populations of cell, with one of them being epithelial cells and another being mesenchymal cells. Hence there is a total of n − 2 groups that represent the cell populations of different intermediate stages. In this paper, we discussed the population dynamics for four different systems: two-state, three-state, four-state, and five-state EMT processes ( Figure 1). In each system, each population of cells was characterized by death rate and transition rates from/to another population. The populations that corresponded to the intermediate states were given stem cell features, described by self-renewal rates and lower death-rates than other types of cells. We modeled the population of each cell phenotype with a linear ordinary differential equation (ODE), where we introduced a constant influx of cells to the mesenchymal population under the assumption that mesenchymal cells are mobile and invasive. With this influx, all of the matrix systems in our models are nonsingular. The ODE that depicts the change in the population size of a particular cell phenotype P i , i = 1, 2, ..., n, is as followed: where α ki is the cellular conversion rate from P k to P i , d i is the death rate, and s i is the self-renewal rate where s i = 0 if the phenotype of the cell population P i is epithelial or mesenchymal. The cellular conversion rates, death rates, and selfrenewal rates all assume constant values. In addition, if P i describes the population of mesenchymal cell, then with I M being the constant influx of cells into the mesenchymal population. For our simulations, we chose the range for our parameter values to be between 1e−5 and 10, with the death rates d i of the "stem cell" populations to be between 1e−5 and 1e−2 to reflect the lower death rates of stem cells. Meanwhile, since the epithelial and mesenchymal populations were assumed to have higher death rates, we chose the death rates to be between 0.05 and 10. The quantity of cell population is in arbitrary unit. One time unit in our model corresponds to one day. In order to explore the noise attenuation properties of each population dynamics model, we introduced multiplicative noise σ i P i dW i to our ODE system via multiple means, where σ i is the noise coefficient and W i describes the Standard Brownian motion. This noise is referred to as "input noise" in this paper. First, we added noise to the epithelial and mesenchymal populations, then observed the effects that noise introduction had on each individual population as well as overall. Next, we added noise to one intermediate population besides the two existing noises on epithelial and mesenchymal cells, and finally we introduced multiplicative noise to each cell population. To ensure consistency between the population fluctuations, we chose the noise coefficient σ i = 1 for all i = 1, 2, ...n.

2.2.
Model comparison through differential evolution algorithm. To compare various models with distinct number of parameters, we used Differential Evolution (DE) to optimize these models with respect to their ability to attenuate noise, and we subsequently compared the optimized parameter sets in terms of the noise attenuation property. DE is a method for improving a set of parameter values (i.e., a parameter vector x, as above) with respect to an objective function, Φ(x), by generating a sequence of trial parameter vectors by processes of reproduction and selection [26]. Reproduction generates an offspring parameter vector u from a parent parameter vector x by diversification. If the offspring performs better than the parent, then the cost value produced by the objective function for the offspring is less than that of the parent, i.e. Φ(u) < Φ(x). In this case, the parent vector x is replaced by the offspring vector u in the next generation. The DE algorithm was previously shown to be efficient in optimizing models for biological systems [13]. With this algorithm, we iteratively searched for a newer set of parameter values that provided a better cost to the objective function Φ(x) than the previous set until the cost converged to the most optimal value. We chose the cost to be the average coefficient of variation of the population size for all phenotypes, i.e. the sum of the coefficients of variation of all the cellular populations divided by the number of states in a particular EMT system that we study. As a result, a lower cost value implies an overall reduction in fluctuations in each cell population and is thus desirable. We applied the optimization algorithm ten times to obtain ten sets of parameters with their corresponding cost values.
During DE, parameter vectors evolve from generation to generation. Each generation (indexed by t = 0, 1, ...) consists of N parameter vectors x j , j = 1, ..., N . Hence, the real number x ij (t) is the value of the i th parameter in the j th parent in the t th generation. Let u j (t) be the parameter vector for the single offspring of the j th parent in the t th generation. The components of this vector, u ij (t) for i = 1, ..., D, are constructed in two steps (mutation and crossover). Then, given the two parameter vectors x j (t) and u j (t), a decision is made as to which one is propagated to generation t + 1. The following three operations propagate parameter vectors from one generation to the next: 1. Mutation: First, we create a "mutant" vector v j (t) by perturbing a parental vector x j (t). In our DE approach, we let the perturbation vector be the difference between the parameter vectors of two additional parents, j and j , chosen at random from the t th generation of parents. All three parents must be different. The "mutant" vector is defined by: where F is a scalar, 0 < F < 1, that determines how "aggressive" the mutation is. 2. Crossover: Next we allow for crossover between the parental parameter set x j (t) and the mutant parameter set v j (t). Component-wise, the offspring vector u j (t) receives a parameter value from the mutant vector with probability C (the crossover probability) or from the parent vector with probability 1 − C: i = 1, 2, ..., D and j = 1, 2, ..., N.
We choose C = 0.5 so that parental values and mutant values have equal chances to be in the offspring vector. 3. Selection: Depending on their relative performances, x j (t) or u j (t) passes on to the next generation. There are two possibilities here. For "greedy" selection, the offspring replaces its parent if it is superior: For "non-greedy" selection, the condition for replacement is: We implement DE using the Matlab code available in Price and Storn's algorithm for DE [38].

Application of fluctuation dissipation theorem.
To measure the fluctuations on each cell population, we calculated the coefficient of variation of the population size for each cell phenotype from the covariance matrix that we found using the Fluctuation Dissipation Theorem [17] [18]. The Fluctuation Dissipation Theorem helps to predict the behaviors of a system in thermal equilibrium by establishing a relation between the frictional force and random force created by the Brownian particle's impacts with surrounding molecules. We applied the version of the Fluctuation Dissipation Theorem mentioned in [32] to our system of linear stochastic ODEs: where A is a constant matrix and Γ(P(t))dW t represents the multiplicative noise term. Γ(P(t)) is a diagonal matrix whose i th entry is σ i P i (t). With this application of the Fluctuation Dissipation Theorem, we found the covariance matrix Σ for the equilibrium system through the formula with P ss (t) being the steady states of the deterministic system dP/dt = AP(t).
To quantify the noise attenuation performance of each EMT system, we aggregate the results of ten different optimization simulations using the DE algorithm [26] and calculate the mean and standard deviation of the coefficients of variation that describe the fluctuations in a specific population. Figures 2-8. 2.4. Sensitivity analysis. In our sensitivity analysis, we defined our parameters as followed: k ij describes the cellular transition rate from the population in the i-state to the population in the j-state where i, j = E, M, I, 1, 2, and 3. Here, E denotes epithelial state, M mesenchymal state, I the intermediate state in the three-state EMT system, and 1, 2, and 3 denote the three intermediate states I 1 , I 2 , and I 3 in the four-state and five-state systems. The parameter k sI is the selfrenewal rate of the intermediate population in the three-state EMT system while the parameters k s1 , k s2 , and k s3 are the self-renewal rates of the I 1 , I 2 , and I 3 populations in the five-state process.
The sensitivity of each parameter is measured by the average change in the average coefficient of variation upon perturbation by 50%. For some transition rates, a positive change in the parameter value due to perturbation results in a significantly smaller average coefficient of variation but this parameter value can be discounted because it lies outside the designated value range that we used in our optimization algorithm  Figure 2B). We performed similar analysis on the fraction of E cells, and we found that when the input noise was only in E and M states, the three-state model performed significantly better than the two-state in terms of minimizing the fluctuations in E cells ( Figure 3A We quantified the average fluctuations of all cellular states with the two-state and the three-state models. When we added noise on E and M states, the average cellular fluctuations were markedly reduced in the presence of the intermediate state ( Figure 4A), and there was still a moderate reduction of average fluctuations when the noise was on all cellular states ( Figure 4B). This suggests that the overall performance of the three-state model is better than that of the two-state model in terms of stability of phenotypic compositions.
We next asked which cellular state transition rates are critical for maintaining the stability. By perturbing the parameters representing the rates transitions among the three cellular states (input noise was added to all three states), we found that the performance of noise attenuation in the three-state system is most sensitive to three parameters: the transition rate from I to E (k IE ), the transition rate from I to To obtain the normalization, we perform stochastic simulations on the steady state population of mesenchymal cells, then divide the mesenchymal population size at each time point by its mean value obtained over the 10-day period. In the middle two panels, the normalized mesenchymal population size is plotted against the number of times that particular size occurs. Yellow: two-state EMT, green: three-state EMT. Bottom two panels display the quantification of the noise attenuation performance of the two-state and three-state systems using the mean and standard deviation of the coefficients of variation (CV) of the mesenchymal population. The mean is plotted here in the form of a bar graph (blue), while the standard deviation is described by the red error bar.  M (k IM ) and the transition rate from M to I (k M I ) ( Figure 4C). These results suggest that the transitions involving the intermediate states are important for noise attenuation, and they further corroborate the advantage of the intermediate state.
They correlation between these transition rates and the ability of the system to attenuate noise may be tested in future experiments. Here, the normalized mesenchymal population size is plotted against the number of times that particular size occurs. Green: three-state EMT, blue: four-state EMT, red: five-state EMT. Bottom three panels display the quantification of the noise attenuation performance of the three-state, four-state, and five-state systems. Here, the normalized average intermediate population size is plotted against the number of times that particular size occurs. The color coding scheme for each system is similar to that of Figure 5. Bottom three panels display the quantification of the noise attenuation performance of the three-state, four-state, and five-state systems.

3.2.
Multiple intermediate states enhance noise attenuation. Next, we explored the advantages of modeling the EMT system with multiple steady states. Using identical input noise on E and M populations, we examined the noise buffering property of three different EMT models: three-state, four-state, and five-state transitions. Analysis of the fluctuations in M cells demonstrates a reduction in noise effects in systems with more than one intermediate states, as evident by the normalized trajectories of the fraction of M cells and the coefficient of variation that reflects the variability of the mesenchymal population size ( Figure 5A). In addition, fivestate system filters noise more efficiently than four-state system. Similarly, when we introduced noise to one intermediate population in addition to the E and M Figure 7. Comparison of the noise attenuation property between the three-state, four-state, and five-state EMT systems in terms of the effects on the epithelial population size. A) Multiplicative noise is introduced to epithelial and mesenchymal populations. B) Multiplicative noise is introduced to the epithelial, mesenchymal, and one intermediate populations. C) Multiplicative noise is introduced to all the populations. Top three panels are the time-course trajectories that represent the normalized number of epithelial cells N E /µ(N E ) over a period of 10 days. Middle three panels illustrate the distribution of the different population sizes of the normalized epithelial population. Here, the normalized epithelial population size is plotted against the number of times that particular size occurs. The color coding scheme for each system is similar to that of Figure 5. Bottom three panels display the quantification of the noise attenuation performance of the three-state, four-state, and five-state systems.
populations, we observed the same correlation between the number of intermediate populations in an EMT system and the noise attenuation property of that system ( Figure 5B). We also confirmed this result by adding noise to all the populations in every EMT system, where the five-state system achieved the best noise filtering results, followed closely by the four-state system ( Figure 5C). From analyzing the responses of the average intermediate population to noise with similar approach, we found that the populations of the five-state and four-state systems outperformed Figure 8. Comparison of the noise attenuation property between the three-state, four-state, and five-state EMT systems in terms of the effects on the average population size. A) Multiplicative noise is introduced to epithelial and mesenchymal populations. B) Multiplicative noise is introduced to the epithelial, mesenchymal, and one intermediate populations. C) Multiplicative noise is introduced to all the populations. Top three panels are the time-course trajectories that represent the normalized average number of cells N avg /µ(N avg ) over a period of 10 days. Middle three panels illustrate the distribution of the different population sizes of the normalized average population. The color coding scheme for each system is similar to that of Figure 5. Bottom three panels display the quantification of the noise attenuation performance of all the EMT systems. D) Sensitivity analysis of the parameters representing unique cell transition rates in the five-state EMT system. Here, we plot the mean of the average change in the average CV as bar graphs (blue) accompanied by red error bars that describe the standard deviation of the average change.
that of the three-state system in attenuating noise at the intermediate levels, regardless of how noise was introduced to the different populations ( Figure 6A-C). This suggests that the number of stem cell sub-states correlates with their ability to reduce the fluctuations of stem cell population. When we studied the effects of noise on the fraction of E cells, we did not observe a similar trend to that of the M and I cells. When input noise was added to E and M populations only, the epithelial population of the five-state system was more susceptible to the fluctuations ( Figure  7A) while the epithelial population of the same system showed the least variability when noise was also added to one intermediate population ( Figure 7B). With noise in every state, the difference in the performance of each system in filtering noise at the epithelial level was not discernible ( Figure 7C). These results suggest that when the number of intermediate state varies, there might be a tradeoff between attenuating noise in E cells and other cells.
Lastly, we quantified the overall performance of each EMT system by considering the fluctuations on the average population of all cellular states. In the scenario where noise is only introduced to the E and M states, the five-state system performs the best at suppressing noise, while the four-state system displays better noise attenuation than the three-state system ( Figure 8A). With the addition of noise in one intermediate state, the four-state and the five-state EMT processes demonstrate improvement in noise attenuation from their three-state counterpart by 6.9% and 9.6% respectively ( Figure 8B). Likewise, we confirmed our results with the addition of noise to all cellular populations and concluded that the EMT system achieves better noise attenuation property when more intermediate states are taken into consideration ( Figure 8C).
Using the five-state EMT process with noise in every state, we analyzed the sensitivity of the parameters that described the cellular transition/self-renewal rates. For our analysis, we examined the average change in the average coefficient of variation of all five populations upon perturbing each parameter. We performed the optimization procedure ten times, whereupon, we carried out the perturbations each time and presented our results in Figure 8D. We found that the parameters that described the transition rate from any intermediate population I 1 , I 2 , or I 3 to the epithelial population (k 1E , k 2E , k 3E ) are most sensitive to their own perturbations, resulting in a marked increase in the average coefficient of variation, therefore a decrease in noise buffering ability of this EMT system. Besides those three transition rates, the cellular transitions from the mesenchymal population to the intermediate , also show considerable sensitivity to perturbations. We thus recognize the crucial roles that the intermediate states such as the I 1 , I 2 , and I 3 states play in ensuring the better noise attenuation performance of the fivestate system. The significance of the transitions between intermediate states and the others reinforce the notion that having multiple intermediate states is beneficial to reducing the effects of the fluctuations on the overall population.
4. Discussion. Studying the multi-step and reversible dynamics of EMT is essential for understanding the roles of EMT in various biological processes such as development and cancer progression. Previous studies revealed multiple intermediate states lying between terminal epithelial and mesenchymal states [12] [14] [29], and this is consistent with the observations that epithelial cell populations show remarkable heterogeneity in normal and tumor tissues [19] [35]. These intermediate states are associated with stemness and invasiveness during cancer progression [9] [15], but the functions of these multiple intermediate states and the complex transitions among all the phenotypic states are unclear. In this study, we used mathematical models to show that systems with multiple intermediate (stem cell) states have advantages in attenuating fluctuations of the heterogeneous cell population, and this effect is particularly beneficial for maintaining stable fractions of stem cells in cell populations, thus providing a strategy for maintaining homoeostasis in facing stochasticity of cell fate changes. The possible performance objective at systems level suggests a new design principle for multi-step EMT and other multi-state systems involving transitions among diverse types of cells.
Complex networks in biological systems have been extensively studied in terms of design principles and performance objectives [21] [31]. However, most of these studies focused on intracellular networks of molecular interactions and influence, and much less is known for networks formed by phenotypic transitions among cell types [10]. We showed that the later type of networks can carry rich information in terms of performance objectives such as facilitating noise attenuation. This finding improves our understanding of systems design at tissue level, and suggests that the formation of cellular states and the transitions among them may have evolved to reach performance objectives that are beneficial to the organisms. It is likely that multiple performance objectives are influencing the evolution of the cellular networks, and further studies are needed to identify these objectives.
Due to the stochastic nature of biochemical reactions and fluctuating environmental conditions, many biological systems are designed to attenuate noise such that stable response to signals can be attained [6] [33]. Intermediate states and regulators have been previously found to improve the robustness of developmental patterning [3]. Our finding that adding intermediate cell states can reduce the fluctuations in cell population provides another strategy for noise attenuation at tissue level. The fluctuations that we included in our simulations represent the stochasticity of cell division [7], cell death [8] and phenotypic transitions [4], all of which were previously observed in experiments. Nonetheless, the sources of these types of noise are mostly molecular fluctuations that affect various types of rates of the cellular activities [37]. In order to describe the dynamics of the noise in more details, future models will need to incorporate the stochastic dynamics of biomolecules into the framework.
Multiple intermediate cell states may arise from complex gene regulatory networks with interconnected positive feedback loops [12] [13]. It is conceivable that the formation of more intermediate cells would require more complex gene regulatory networks, which in turn need some other strategies and/or more energy to control. Therefore, the correlation between the number of states and the ability to attenuate noise suggests that a tradeoff between the simplicity of the gene regulatory network and noise attenuation may exist when the system is subject to design via evolution.
Dynamic equilibrium of sub-states of stem cells is observed not only in the EMT system, but also in hematopoietic and embryonic stem cell populations [11] [27] [34]. Moreover, it was suggested that heterogeneity of cell population might be related to disseminated cancer cell dormancy [28]. Although these sub-states of stem cells might represent distinct functional entities in specific contexts, the recurring phenomena indicate that there might be advantages of this dynamic behavior at systems level. Our modeling work provides a general framework to study the cell differentiation systems with multiple sub-states of stem cells, and our conclusions with respect to the noise attenuation property of multiple stem cell states can apply to other systems, such as embryonic stem cell and disseminated cancer cell [11] [28]. In addition, feedback control among multiple cellular states have been shown to be critical to the growth and stability of cell populations [20] [22] [23]. Future works are warranted to investigate the role of feedback controls involving the multiple intermediate states.
In conclusion, our computational analysis shows that the existence of multiple intermediate states between epithelial and mesenchymal states and the transitions among these states are able to attenuate the fluctuations of the fractions of cell population. We found a general correlation between the number of states and the ability to attenuate fluctuations. In particular, the fluctuations of the stem cell populations are reduced by increasing the number of intermediate states. These results improve the understanding of the intermediate states in the EMT system in terms of performance objectives, and provide insights into the stem cell systems with multiple sub-states in general.