Computational Modeling of the Crosstalk Between Macrophage Polarization and Tumor Cell Plasticity in the Tumor Microenvironment

Tumor microenvironments contain multiple cell types interacting among one another via different signaling pathways. Furthermore, both cancer cells and different immune cells can display phenotypic plasticity in response to these communicating signals, thereby leading to complex spatiotemporal patterns that can impact therapeutic response. Here, we investigate the crosstalk between cancer cells and macrophages in a tumor microenvironment through in silico (computational) co-culture models. In particular, we investigate how macrophages of different polarization (M1 vs. M2) can interact with epithelial-mesenchymal plasticity of cancer cells, and conversely, how cancer cells exhibiting different phenotypes (epithelial vs. mesenchymal) can influence the polarization of macrophages. Based on interactions documented in the literature, an interaction network of cancer cells and macrophages is constructed. The steady states of the network are then analyzed. Various interactions were removed or added into the constructed-network to test the functions of those interactions. Also, parameters in the mathematical models were varied to explore their effects on the steady states of the network. In general, the interactions between cancer cells and macrophages can give rise to multiple stable steady-states for a given set of parameters and each steady state is stable against perturbations. Importantly, we show that the system can often reach one type of stable steady states where cancer cells go extinct. Our results may help inform efficient therapeutic strategies.


INTRODUCTION
Cancer has been largely considered as a cell-autonomous disease, but recent investigations have highlighted the crucial role of the tumor microenvironment in determining cancer progression (1). Cancer cells can communicate bi-directionally through various mechanical and/or chemical ways with their neighboring cancer cells (2,3), and/or with other components of the tumor microenvironment such as macrophages and fibroblasts, driving aggressive malignancy (4)(5)(6). The interconnected feedback loops formed by these interactions can often generate many emergent outcomes for the tumor. Interestingly, many of the latest therapeutic innovations such as immunotherapy are aimed at targeting aspects of the tumor microenvironment instead of the cancer cells (7).
Tumor-associated macrophages (TAMs) are one of the most abundant immune cell populations in the microenvironment (8,9). They have been shown to promote cancer progression in many ways, such as promoting angiogenesis, suppressing function of cytotoxic T lymphocytes, and assisting extravasation of cancer cells (8)(9)(10)(11)(12). Generally, the secretome and functions of TAMs have been shown to be close to that of the socalled alternatively activated macrophages (M 2 ) (13). In the case of pathogen infections, macrophages that can engulf the pathogen and present processed antigens to adaptive immune cells, are generally characterized as the classically activated ones (M 1 ) (14). M 1 and M 2 macrophages have different roles during wound healing: while M 1 macrophages initiate inflammatory responses, M 2 macrophages contribute to tissue restoration (13). In the context of cancer, M 1 macrophages have been generally considered anti-tumor (15)(16)(17), whereas M 2 macrophages have been considered as pro-tumor (10).
However, macrophage polarization is not as rigid as the differentiation of T lymphocytes (18); instead, M 1 , M 2 , and any intermediate state(s) of macrophage polarization are quite plastic (13,19,20). Thus, the idea that reverting TAMs in the cancer microenvironment to its cancer-suppressing counterpart is tempting, the proof of principle of which has been demonstrated in mice models (21)(22)(23)(24)(25).
Importantly, macrophages and cancer cells can interact with and influence the behavior of one another, as shown by many in vitro experiments. Specifically, some epithelial cancer cells are capable of polarizing monocytes into M 1 -like macrophages (33). Forming a negative feedback loop, these M 1 -like macrophages can decrease the confluency of the cancer cells that polarized them (33). Moreover, pre-polarized M 1 macrophages can induce senescence and apoptosis in human cancer cell lines A549 (34) and MCF-7 (35). Intriguingly, factors released by apoptosis of cancer cells can convert M 1 macrophages into M 2 -like macrophages (35), thus switching macrophage population from being tumor-suppressive to being tumor-promoting. On the other hand, mesenchymal cancer cells can polarize monocytes into M 2 -like macrophages (33,36,37), that can in turn assist EMT (37,38). Thus, the interaction network among macrophages and cancer cells is formidably complex, and the emergent dynamics of these interactions can be non-intuitive (39) yet are often crucial in deciding the success of therapeutic strategies targeting cancer and/or immune cells. For example, even if TAMs at some time can be converted exogenously to M 1like macrophages, if most cancer cells still tend to polarize  monocytes to TAMs, other coordinated perturbations may be  needed at different time-points to restrict the aggressiveness of  the disease. Here, we develop mathematical models to capture the abovementioned set of interactions among cancer cells in varying phenotypes (epithelial and mesenchymal) and macrophages of different polarizations (M 1 -like and M 2 -like). We characterize the multiple steady states of the network that can be obtained as a function of different initial conditions and key parameters, and thus analyze various potential compositions of cellular populations in the tumor microenvironment. This in silico coculture system can not only help explain in vitro multiple experimental observations and clinical data, but also help acquire novel insights into designing effective therapeutic strategies aimed at cancer cells and/or macrophages.

Crosstalk Among Cancer Cells and Macrophages Can Lead to Two Distinct Categories of Steady States
We first considered the following interactions in setting up our mathematical model: (a) proliferation of epithelial and mesenchymal cells (but not that of monocytes M 0 , or macrophages M 1 and M 2 ), (b) EMT promoted by M 2 -like macrophages and MET promoted by M 1 -like macrophages, (c) polarization of monocytes (M 0 ) to M 1 -like cells aided by epithelial cells, and that to M 2 -like cells aided by mesenchymal cells, (d) induction of senescence in epithelial cells by M 1like macrophages ( Figure 1A). No inter-conversion among M 1like and M 2 -like macrophages or cell-death of macrophages is considered here in this model (hereafter referred to as "Model I"; see section Methods). Furthermore, this model also considers that mesenchymal cells can secrete soluble factors, such as TGFβ, that can induce or maintain the mesenchymal state in autocrine or paracrine manners (40,41). Therefore, we hypothesized that mesenchymal cells can resist M 1 -promoted MET. However, this resistance might not fully suppress MET, as MET still happens in the presence of M 1 macrophages (38). Therefore, in Equation (1), we assume a Hill-like function to represent the resistance to M 1 -promoted MET and the resistance term saturates to a finite value as a function of mesenchymal population M. And the corresponding half-saturation constant is K 0 M . Similarly, epithelial cells adhere to each other via E-cadherin, sequestering β-catenin on the membrane, thus interfering with the induction of EMT (42). Thus, we hypothesized that the M 2 -promoted EMT can be inhibited by epithelial cells in a cooperative manner, because this inhibition of EMT requires direct physical cell-cell contact and hence involves multiple epithelial cells. Therefore, in Equation (1), we assume a Hill-like function to represent this resistance to M 2 -dependent EMT and the resistance term saturates to a finite value as a function of epithelial population E. And the corresponding half-saturation constant is K 0 E . Furthermore, the epithelial-cell-dependent term has a cooperativity parameter k to represent the effects of Ecadherin-β-catenin interaction between multiple epithelial cells.
Note that in all of our calculations, we assumed that there is a carrying-capacity of cells (N max , including both cancer cells and macrophages) in the co-culture system in vitro. We vary the initial number of different cells while keeping the total number of macrophages (M 1 +M 2 +M 0 ) to be a constant (M c ). Therefore, the maximum number of cancer cells will be N max -M c .
In Model I, the final populations of E, M, M 1 , and M 2 cells are simply determined by the following equations: where η em and η me are the EMT and MET rate constants, respectively. The above equations give two categories of stable steady states: (a) state I, dominated by epithelial cancer cells, and (b) state II, dominated by mesenchymal cancer cells. The final steady state of the system depends on the number of M 2 macrophages in that state; since there are only three equations for four unknowns, M 2 can be used as a parameter specifying (possibly discrete set of) solutions. As soon as M 2 /N max crosses a threshold, the system switches from state I to state II ( Figure 1B). This prediction is largely robust to parameter variation ( Figure S1). Next, we explored what factors could change the qualitative behavior of the model, i.e., enable a smoother and continuous transition of epithelial and mesenchymal percentages as a function of the M 2 -macrophage population. We identified that reducing the cooperativity of epithelial cancer cells in their resistance of M 2 -promoted EMT can lead to a loss of the feature with two-types of steady states (Figure 1C, blue, cyan, and light-green lines). Conversely, increasing the cooperativity of M 2 -promoted EMT or M 1 -promoted MET can expand the region of overlapping between state I and II of the system ( Figure 1C, green lines). Another factor that can alter the behavior of the model is the initial number of monocytes in the system. At high enough number of monocytes, the absolute number of cancer cells will be small, thus the effect of the cooperativity between epithelial cancer cells will be reduced, and consequently, as discussed above, the overlap of the two types of steady states of the system will disappear ( Figure 1D). Together, this increased propensity of multi-stability in the system upon including cooperative effects in the interactions among different species (or variables) is reminiscent of observations in models of biochemical networks (43).
Note that using the M 2 /N max as the "control variable" is specific for Model I, because there is no interconversion between M 1 and M 2 macrophages. For a given set of parameters, the steady state level of M 2 can vary continuously from 0 to M c (Figures 1B-D), and in practice is determined by the initial conditions (initial number of E, M, M 1 , M 2 , and M 0 ). For models in the following sections, M 2 /N max will be shown to be fixed to discrete allowed values (instead of continuously varying) for a given set of parameters. For example, if we add a very small conversion rate between M 1 and M 2 , the steady state of the system will collapse to only one steady state ( Figure S2): on the shorter time scale, the trajectory of the system (on the phase plane of E and M 2 ) will first converge to one steady state in Model I (with the same M 2 /N max ); on the longer time scale, as determined by the value of inter-conversion rate between M 1 and M 2 , the system will slowly evolve to the steady state (blue dot in Figure S2), following along the now-approximate steady states in Model I.

Cancer-Cell Enhanced Interconversion Between M 1 -and M 2 -Like Macrophages Lead to Bi-Stability
In the next iteration of our model (hereafter referred to as "Model II"), we included the possibility of interconversion between M 1like and M 2 -like macrophages, as reported in the literature (13,35,44,45) (Figure 2A, see section Methods). We first assume constant interconversion between M 1 and M 2 (with rates denoted as η 0 21 and η 0 12 ) and investigate the effects of varying the rate of conversion of mesenchymal cells to epithelial cells (η me ). At small values of η me , the system has small number of epithelial cells, which is equivalent to state II in Model I; with increasing η me , a threshold is crossed, and the system can switch to states with larger number of epithelial cells, which is equivalent to state I in Model I (solid black lines, Figure 2B). Thus, we can observe bi-stability of cancer cell population in this system. However, the populations of macrophages stay constant as a function of η me , since the M 1 /M 2 ratio is simply determined by η 0 21 /η 0 12 according to Equation (2) when η 12 = η 21 = 0). In this model, we focused on increased cooperativity of M 2 -assisted EMT and M 1 -assisted MET (i.e., m = n = 2), because the interconversion between M 1 and M 2 , in absence of such cooperativity (as considered in model I with m = n = 1; Figure 1B), gives rise to a narrower bistable region ( Figure S3). Note again that in our calculations, we assumed that there is a carrying-capacity of cells (N max , including both cancer cells and macrophages) in the co-culture system in vitro. Again, the total number of macrophages is constant, since no cell death or cell division of macrophages is considered.
We chose η me as a control parameter, because it can act as a bottleneck for the transition between the two types of stable steady states of the system. Lowering the transition threshold of η me can be helpful in a sense that the system can potentially stay at state I (high epithelial state, supposed to be less aggressive) only. Due to the inherent symmetry in the network, the effect of lowering the threshold of η me can be recapitulated via other perturbations, such as a smaller rate of M 1 -assisted MET (η em ), lower M 1 to M 2 conversion rate (η 12 ), or higher M 2 to M 1 conversion rate (η 21 ) (Figure S4).
We next consider the case where interconversion between M 1 -like and M 2 -like macrophages is enhanced by cancer cells, i.e., mesenchymal cells enhance the M 1 -to-M 2 transition, while epithelial cells enhance M 2 -to-M 1 transition. In this case, we observe the existence of a bi-stable region, and the range of the epithelial-cell-low solution (equivalent to state II) is wider than that for the previous case without any effects of epithelial (E) and mesenchymal (M) cancer cells on M 1 -M 2 interconversion (Figures 2B,C, solid blue lines). The reason for the expanded epithelial-cell-low region is that the positive feedback loop between mesenchymal and M 2 cells makes the mesenchymaland M 2 -dominated state more stable. Therefore, a higher η me is required to compensate for the effects of low M 1 population. For therapeutic purposes, state I is favored, for it is believed that epithelial cells are typically less aggressive than mesenchymal ones. Therefore, symmetrical mutual enhancement might not be helpful for the therapy because of the expanded epitheliallow (mesenchymal-high) region. For the same reason, the lower threshold of η me to switch from state II to state I also shift to a higher value, which means that for a small η me , the system will only have one stable steady state, i.e., state II with high number of mesenchymal cells.
In order to drive the system to be bi-stable at a smaller η me , we tested the case with asymmetrical interconversion rates between M 1 and M 2 , i.e., higher conversion rate from M 2 to M 1 enhanced by epithelial cells (η 21 ). In this case, the lower threshold of η me can indeed shift to a smaller value (Figure 2D), which makes it theoretically possible to switch the system from state II (low epithelial cells) to state I (high epithelial cells) at a fairly small η me value.
In summary, our results for Model II suggest that one can switch the system from state II (low epithelial cells) to state I (high epithelial cells) by increasing both η me and the effective conversion from M 2 to M 1 . Note that the interactions in Model 2 are rather symmetric: the interconversion rates between M 1 and M 2 are either constants or enhanced by the corresponding cancer-cells. In the following section, we will consider the asymmetric case where only the M 1 to M 2 conversion is enhanced by factors released by apoptotic epithelial cancer cells.
Cell Apoptosis-Induced M1-to-M2 Conversion Leads to Symmetry Breaking in the Cancer-Immune Interaction Network Finally, we incorporated one other set of experimentally documented interactions, i.e., M 1 macrophages can drive the apoptosis of epithelial cells, and factors released during cancer cell apoptosis can drive M 1 -to-M 2 conversion (referred as Model III, see section Methods) (35). The newly introduced interactions are highlighted in thick lines in Figure 3A. This interaction induces "symmetry breaking" (46) in the system, as previously most of the interactions considered were of a "symmetric" nature, i.e., M 1 cells driving MET and M 2 cells driving EMT, and epithelial cells driving M 1 maturation while mesenchymal cells driving M 2 maturation. With this new interaction, the system is now biased against epithelial cells, because (a) epithelial but not mesenchymal cells, are killed by M 1 -macrophages, and (b) their dead counterparts may convert M 1 to M 2 cells that can, in turn, convert some epithelial cells to mesenchymal ones (EMT). In addition, the conversion from M 2 to M 1 is essential to bring the system back from the mesenchymal-cancer-cell biased state and it can be enhanced by introducing IL-12 (45) or Type 1 T helper cells (44).
Thus, in the parameter region investigated, there are 3 types of stable steady states, which are represented by solid blue, red and black lines, respectively. At small values of the control parameter η me (rate of M 1 macrophage assisted MET), there is only one type of stable steady state: the system is biased toward the mesenchymal-dominated state (solid blue curves in Figures 3B-E), whereas the cancer-extinction state (E = M = 0, dashed black curves in Figures 3B-E) is unstable. As η me increases and goes across a critical value η a me =0.0626 h −1 ), the extinction state (E = M = 0) becomes stable as well as a new set of steady states emerges (red lines in Figures 3B-E).
Between η me =0.0663 h −1 and η me =0.0747 h −1 , the steady states depicted as a solid red line (Figures 3B-E) are stable; here both populations of epithelial and mesenchymal cancer cells are at a low level. This phenomenon can be understood in the following sense: at higher η me , proliferating mesenchymal cells are continuously being converted to epithelial cells, which will be killed by M 1 macrophages. This effect brings down the overall fraction of cancer cells. Note that in this region, three types of stable steady states co-exist in the system for a given set of parameters. As η me increases and goes across η b me =0.0747 h −1 , the stable steady states in solid red lines disappears and the other two types of stable steady states coexisted (solid blue and black lines in Figures 3B-E). As η me further increases and goes across η c me =0.1532 h −1 , there is only one stable steady state: cancer cells are necessarily eliminated from the system. We note in passing that the instability that drives the red solution unstable as η me is lowered past 0.0663 h −1 is a Hopf bifurcation to an unstable limit cycle.
Furthermore, for the breast cancer cell line MDA-MB-231 used in Yang et al. (38), η me is estimated to be around 1/120 h −1 (∼0.0083 h −1 ). In order to explore the conditions for the absolute extinction of cancer cells around the estimated η me , we varied different parameters (such as λ M , η em, and η 21 ) to get the corresponding η c me where the extinction state is the only stable steady state of the system. We found that lowering the growth rate of mesenchymal cells (λ M ) can reduce η c me to the nominal experimental value 1/120 h −1 (Figure 3F) whereas lowering down η em or increasing η 21 might not ( Figure S5). Therefore, for therapeutic purposes, increasing η me and decreasing λ M can be a promising combination to help to eliminate cancer cell populations. In addition, increasing the number of macrophages in the co-culture system can also shift η c me to a lower value ( Figure S6) because of a potentially higher level of cancer-killing M 1 -like macrophages.

EMT Scores Correlate With Multiple Genes Upregulated in M2 Macrophages
As a proof of principle for the predictions of our model, we investigated multiple TCGA datasets (see section Methods for the source of datasets), using our previously developed EMT scoring metric (47). This metric quantifies the extent of EMT in a particular sample, and correspondingly assigns a score between 0 (fully epithelial) and 2 (fully mesenchymal). We calculated the correlation coefficients for EMT scores with various genes reported to be differently regulated in M 1 or M 2 macrophages relative to M 0 macrophages. The list of those genes is taken from Jablonski et al. (48). Out of 52 genes investigated, we observed that many genes upregulated in M2 and downregulated in M 1 macrophages-ACTN1, FLRT2, MRC1, PTGS1, RHOJ, TMEM158-correlated positively with the EMT scores ( Table 1, first 6 rows) across multiple cancer types. The higher the EMT scores, the higher the levels of those genes, including the canonical M 2 macrophage marker CD206 (MRC1). On the other hand, out of the 52 genes investigated, many genes upregulated in M 1 macrophages but downregulated in M 2 macrophages-FIIR, STAT1, RSAD2, TUBA4A, and XAF1-showed either a negative or an overall weak positive correlation with EMT scores ( Table 1, last 5 rows). The only exception observed in this trend was that for ARHGP24 which correlated positively with EMT scores across cancer types. Together, these correlation results in multiple TCGA datasets offer a promising initial validation of our model predictions that a state dominated by epithelial cells typically has higher number of M 1 macrophages, while the other state dominated by mesenchymal cells typically has higher number of M 2 macrophages. Furthermore, to go beyond the correlation analysis between a single gene and the EMT-score, we investigated whether genes that are expressed higher in M 2 -like macrophages will be enriched in EMT-score-high tumors and whether genes that are expressed higher in M 1 -like macrophages will be enriched in EMT-score-low tumors. Again, the lists of genes that are expressed higher in M 1 -or M 2 -like macrophages are from Jablonski et al. (48) and the EMT-scores were calculated for each sample of the four TCGA dataset investigated in this work. We define a tumor to be EMT-score-high if its EMT-score is higher than the median score of the specific dataset and vice versa. Then, the expressions of a gene in EMT-score-high and low tumors in each dataset are used to determine whether this gene has a significantly high/low expression in EMT-score-high/low tumors. As expected, the genes that are expressed higher in M 2 -like macrophages (but lower in M 1 -like macrophages) are enriched in EMT-score-high tumors across the 4 types of cancers investigated ( Table 2). This result is consistent with our model predictions: M 2 -like macrophages and mesenchymal cancer cells tend to stably co-exist. However, somewhat to our surprise, the genes that expressed higher in M 1 -like macrophages (but lower in M 2 -like macrophages) are also enriched in EMT-score high tumors across the 4 types of cancers investigated (Table S1). This result is not in line with our model prediction where M 1 -like macrophages are expected to stably co-exist with epithelial cancer cells. Possible reasons for this inconsistency will be discussed later. Nevertheless, the gene enrichment analysis in multiple TCGA datasets offer a promising initial validation of (most of) our model predictions.

DISCUSSION
The tumor microenvironment involves multiple cell types that interact among each other in diverse ways, thus giving rise to a complex adaptive ecological system (49)(50)(51). For such a system, mathematical models can help reveal the mechanisms underlying the emergent behavior and eventually aid in designing effective therapeutic strategies to modulate those aspects of the microenvironment that fuel disease aggressiveness (52)(53)(54)(55)(56)(57)(58).
Here, we focused on the interactions among macrophages of different polarizations (M 1 and M 2 ) and cancer cells with different phenotypes (epithelial and mesenchymal). Based on the literature, we focused on two types of models: with (Model II and III) or without (Model I) interconversion between M 1 and M 2 macrophages. All three models share a common feature: with a given set of parameters, multiple types of stable steady states can co-exist. More specifically, in Model I (without M 1 -M 2 interconversion), with a given set of parameters, the system 2 | M 2 -associated gene set enrichment: Reported are all genes (out of 31 total) having statistically significant differential expression (samples segregated based on median EMT score) across TCGA colorectal adenocarcinoma (Adeno COL; n = 286), prostate adenocarcinoma (Adeno PCA; n = 497), invasive breast cancer (BCA invasive; n = 1097), and lung adenocarcinoma (Lung Adeno; n = 515). can converge to continuous range of states depending on the initial condition. However, these steady states belong to two categories: state I with a higher epithelial population and state II with a lower epithelial population. After the system reaches a steady state, perturbations applied only to cancer cell populations cannot drive the system out of the original steady state whereas perturbations on macrophage populations will drive the system out of the original steady state. However, the perturbation on macrophage populations might not drive the system out of the original category of steady states unless the perturbation is sufficiently strong. In Model II and III, with a given set of parameters, the system can again converge to two types of steady states depending on the initial condition: state I with a higher epithelial population and state II with a lower (or even zero for Model III) epithelial population. Now, however, the states are discrete. After the system reaches a steady state, perturbations of any single cell population might not drive the system out of the original steady state. Therefore, in general, it is not easy to switch the cancer-macrophage system from a mesenchymal-and M 2 -dominated state to an epithelial-and M 1 -dominated state.

Adeno COL
Mathematical approaches similar to ours may be useful in explaining, and even predicting, the efficacy of different therapeutic intervention(s) and their combinations. For example, various efforts have been made to switch populations of macrophages from M 2 -into M 1 -like (25), such as depletion of TAMs, reprogramming of TAMs toward M 1 -Like macrophages, inhibition of circulating monocyte recruitment into tumor, etc. The polarization of M 2 macrophages can also be altered via inhibiting the endothelial-mesenchymal transition (EndMT), a process similar to EMT (59). However, it is unclear that how effective these types of strategies would be. Our modeling studies suggest that when we consider the interactions in Model III, the efforts on manipulating the M 1 -M 2 interconversion might not be effective to eliminate cancer cells; whereas in Model II increasing the M 2 -to-M 1 conversion rate can help the system to switch to the epithelial cancer cell and M 1 macrophage dominated state, which is believed to be less aggressive. Therefore, the effective therapeutic strategies strongly depend on the type of interactions present between cancer cells and macrophages.
Furthermore, we attempted to investigate whether the modelpredicted stable co-existence of M 1 -like macrophages and epithelial cancer cells (and likewise M 2 -like macrophages and mesenchymal cancer cells) can be seen in the TCGA dataset. Indeed, we can observe the enrichment of many genes that are expressed higher in M 2 -like macrophages in more mesenchymal tumors. This observation is consistent with our mathematical modeling analysis. However, we also detected the enrichment of genes that are expressed higher in M 1 -like macrophages in more mesenchymal tumors. There can be several reasons for the enrichment of M 1 genes in mesenchymal tumors: (i) the list of genes used here to identify the ones preferentially expressed in M 1 -like macrophages is from a study of mouse macrophages (48). For human macrophages, not as many markers have been validated as those in mice; future experiments can help characterize these better; (ii) for the TCGA sample used, the gene expression data that are used to calculate the EMT-score could be not exclusively that of cancer cells but may incorporate other types of cells, such as fibroblasts. Future experiments that can separate and measure gene expression of cancer cells separately, through techniques such as laser microdissection, will be more precise in quantifying the EMT status of cancer cells.
In addition, it is also important to recognize that our model suffers from limitations. For instance: (a) it does not consider spatial aspects of these interactions, for instance, mesenchymal cells may migrate and invade through the matrix, thus changing the interactions among the cells considered in our framework; (b) it does not consider the effects of senescence on epithelial cell growth (60); and (c) it considers EMT and macrophages polarization as a binary process, whereas emerging reports support the notion that in both cases there is likely to exist a spectrum of states/ phenotypes (61,62). A more realistic model that can overcome the abovementioned assumptions and thus reflect the dynamics of tumor microenvironment more accurately can be used to help designing a more effective way to switch and stably maintain the system into a less aggressive state.
Despite the abovementioned limitations, our model can contribute to identifying key parameters of the system. For example, it suggests that to design an effective therapy to maintain the system in a M 1 -dominated and cancer-free steady state, not only the conversion rate from mesenchymal to epithelial cells should be significant, but also the growth rate of mesenchymal cells should be low enough. In other words, METinducing and cell-growth-suppressing mechanisms can together synergistically restrict disease aggressiveness.
In summary, our results show that the interaction network between tumor cells and macrophages may lead to multi-stability in the network: one state dominated by epithelial and M 1 -like cells, the other dominated by mesenchymal and M 2 -like cells. We also identify that inducing MET and inhibiting cancer-cell growth might be much more efficient in "normalizing" (1) the tumor microenvironment that can otherwise be engineered by cancer cells to support their growth (63).

Computational Models
According to in vitro experiments in literature, we construct three mathematical models. In Model I, factors secreted by epithelial (mesenchymal) cancer cells will polarize monocytes into M 1like (M 2 -like) macrophages. M 1 -like macrophages will induce senescence of epithelial cancer cells and convert mesenchymal cancer cells to epithelial ones. On the other hand, M 2 -like macrophages will convert epithelial cancer cells to mesenchymal ones. In addition, a mesenchymal cancer cell can secrete soluble factors to help maintain the mesenchymal state of itself and its neighbors, whereas an epithelial cancer cell can maintain the epithelial state through being in contact with its neighboring epithelial cancer cells. There is assumed to be no cell growth or death for macrophages and there is a maximum "carrying capacity" of cells in the system. The figure that illustrate this model is in Figure 1A. The equations to describe this system is as follows: