Epithelial-mesenchymal transition (EMT) is an instance of cellular plasticity that plays critical roles in development, regeneration and cancer progression. While many regulatory elements have been identified to induce EMT, the complex process underlying such cellular plasticity remains poorly understood. Utilizing a systems biology approach integrating modeling and experiments, we found multiple intermediate states contributing to EMT and that the robustness of the transitions is modulated by transcriptional factor Ovol2. In particular, we observed that adding the mutual inhibition relationship between Ovol2 and EMT inducer Zeb1 generates a novel four-state system consisting of two distinct intermediate phenotypes that differ in differentiation propensities and are favored in different environmental conditions. 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 used mathematical models to show that multiple intermediate phenotypes in the EMT system help to attenuate the overall fluctuations 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 intermediate states, indicating the stem-cell population is more stable when it has more sub-states. We then attempted to bridge the gap between discrete and continuum modeling of the EMT system by incorporating the EMT core regulatory network into our heterogeneous cell population dynamics model to create a multiscale EMT model. Our model can capture the larger-scale population growth dynamics while acknowledging the intracellular interactions between proteins within each individual cell. From the two types of noise we introduced into our model, we observed that the differences in the noise design prompted distinctive behaviors in the proliferative capability of our heterogeneous population. Our findings also revealed the challenges encountered when integrating noise into a dynamic EMT model such as the multiscale model and the complex role noise plays in modulating the different phenotypic fractions of the population. Lastly, we present a class of semi-implicit integration factor methods and demonstrate its good accuracy, efficiency, and stability properties compared to existing methods. This new class of methods, which is easy to implement, will have broader applications in solving stochastic reaction-diffusion equations arising from models in biology and physical sciences.