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 (M₁ vs. M₂) 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–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–12). Generally, the secretome and functions of TAMs have been shown to be close to that of the so-called alternatively activated macrophages (M₂) (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₁) (14). M₁ and M₂ macrophages have different roles during wound healing: while M₁ macrophages initiate inflammatory responses, M₂ macrophages contribute to tissue restoration (13). In the context of cancer, M₁ macrophages have been generally considered anti-tumor (15–17), whereas M₂ macrophages have been considered as pro-tumor (10).

However, macrophage polarization is not as rigid as the differentiation of T lymphocytes (18); instead, M₁, M₂, 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–25).

Not only TAMs, but also cancer cells themselves can be extremely plastic, a canonical example of which is epithelial-mesenchymal plasticity, i.e., cancer cells can undergo varying degrees of Epithelial-Mesenchymal Transition (EMT) and its reverse Mesenchymal-Epithelial Transition (MET) (26, 27). EMT/MET has been associated with metastasis (28), chemoresistance (29), tumor-initiation potential (30), resistance against cell death (31), and evading the immune system (32).

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₁-like macrophages (33). Forming a negative feedback loop, these M₁-like macrophages can decrease the confluency of the cancer cells that polarized them (33). Moreover, pre-polarized M₁ 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₁ macrophages into M₂-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₂-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₁-like 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₁-like and M₂-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 co-culture 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.

Results

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₀, or macrophages M₁ and M₂), (b) EMT promoted by M₂-like macrophages and MET promoted by M₁-like macrophages, (c) polarization of monocytes (M₀) to M₁-like cells aided by epithelial cells, and that to M₂-like cells aided by mesenchymal cells, (d) induction of senescence in epithelial cells by M₁-like macrophages (Figure 1A). No inter-conversion among M₁-like and M₂-like macrophages or cell-death of macrophages is considered here in this model (hereafter referred to as “Model I”; see section Methods).

FIGURE 1

Figure 1. Cancer-immune interplay can give rise to the co-existence of two types of steady states. (A) Interaction network for Model I. Conversions between cells are indicated by solid lines. Cell proliferation is indicated by dashed lines. Inhibition (in black) and activation (in red) is indicated by dotted lines. (B) Steady states of the epithelial population are plotted as a function of M₂-like macrophage population. Stable steady states are plotted in solid blue lines and unstable steady states are plotted in dotted red line. The key parameters are as indicated. (C) As the cooperativity of epithelial cancer cells in their resistance of M₂-promoted EMT is reduced, i.e., k = 4, 3, 2 (blue, cyan, and light-green line), the overlapped region between state I and II shrinks and then disappears. Increasing the cooperativity of M₂-promoted EMT or M₁-promoted MET, i.e., m and n increase from 1 to 2, can expand the overlapped region (between state I and II) of the system (dark-green lines). Stable steady states are plotted in solid lines and unstable steady states are plotted in dotted lines. (D) As the total number of macrophages (M_c) increases, the overlapped region (between state I and II) of the system shrinks and then disappears.

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₁-promoted MET. However, this resistance might not fully suppress MET, as MET still happens in the presence of M₁ macrophages (38). Therefore, in Equation (1), we assume a Hill-like function to represent the resistance to M₁-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 ${\text{K}}_{\text{M}}^0$. 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₂-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₂-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 ${\text{K}}_{\text{E}}^0$. Furthermore, the epithelial-cell-dependent term has a cooperativity parameter k to represent the effects of E-cadherin-β-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₁+M₂+M₀) 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₁, and M₂ cells are simply determined by the following equations:

$\begin{array}{l}{\eta }_{me}M\frac{M_1^n}{M_1^n+\frac{M}{M+K_M^0}}={\eta }_{em}E\frac{M_2^m}{M_2^m+\frac{E^k}{E^k+{\left(K_E^0\right)}^k}},\\ M_1+M_2=M_c,E+M=N_{\max }-M_c\end{array}$

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₂ macrophages in that state; since there are only three equations for four unknowns, M₂ can be used as a parameter specifying (possibly discrete set of) solutions. As soon as M₂/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₂-macrophage population. We identified that reducing the cooperativity of epithelial cancer cells in their resistance of M₂-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₂-promoted EMT or M₁-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₂/N_max as the “control variable” is specific for Model I, because there is no interconversion between M₁ and M₂ macrophages. For a given set of parameters, the steady state level of M₂ 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₁, M₂, and M₀). For models in the following sections, M₂/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₁ and M₂, 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₂) will first converge to one steady state in Model I (with the same M₂/N_max); on the longer time scale, as determined by the value of inter-conversion rate between M₁ and M₂, 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₁- and M₂-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₁-like and M₂-like macrophages, as reported in the literature (13, 35, 44, 45) (Figure 2A, see section Methods). We first assume constant interconversion between M₁ and M₂ (with rates denoted as ${\eta }_{21}^0$ and ${\eta }_{12}^0$) 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₁/M₂ ratio is simply determined by ${\eta }_{21}^0$/${\eta }_{12}^0$ according to Equation (2) when η₁₂ = η₂₁ = 0). In this model, we focused on increased cooperativity of M₂-assisted EMT and M₁-assisted MET (i.e., m = n = 2), because the interconversion between M₁ and M₂, in absence of such cooperativity (as considered in model I with m = n = 1; Figure 1B), gives rise to a narrower bi-stable 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.

FIGURE 2

Figure 2. Effects of interconversion between M₁ and M₂ cells, mediated by cancer cells. (A) Network of Model II. (B,C) Steady states of the epithelial (B) or M₁ (C) populations are plotted as a function of η_me. Stable steady states are plotted in solid lines and unstable steady states are plotted in dotted lines. The key parameters in (B,C) are: η₁₂ = η₂₁ = 0 (for black lines) or 1/72 h⁻¹ (for blue lines), m = 2, n = 2, k = 4, ${\eta }_{21}^0$ = ${\eta }_{12}^0=$1/72 h⁻¹. (D) In this plot, the key parameters are: η₂₁ = 1/72 h⁻¹ (for blue lines) or 1/36 h⁻¹ (for green lines), m = 2, n = 2, k = 4, ${\eta }_{21}^0$ = ${\eta }_{12}^0$ = 1/72 h⁻¹, and η₁₂ = 1/72 h⁻¹. In (B,D), the gray line is the maximum fraction of cancer cells (=0.7) as M_c is set as 0.3.

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₁-assisted MET (η_em), lower M₁ to M₂ conversion rate (η₁₂), or higher M₂ to M₁ conversion rate (η₂₁) (Figure S4).

We next consider the case where interconversion between M₁-like and M₂-like macrophages is enhanced by cancer cells, i.e., mesenchymal cells enhance the M₁-to-M₂ transition, while epithelial cells enhance M₂-to-M₁ 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₁-M₂ 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₂ cells makes the mesenchymal- and M₂-dominated state more stable. Therefore, a higher η_me is required to compensate for the effects of low M₁ 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 epithelial-low (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₁ and M₂, i.e., higher conversion rate from M₂ to M₁ enhanced by epithelial cells (η₂₁). 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₂ to M₁. Note that the interactions in Model 2 are rather symmetric: the interconversion rates between M₁ and M₂ are either constants or enhanced by the corresponding cancer-cells. In the following section, we will consider the asymmetric case where only the M₁ to M₂ 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₁ macrophages can drive the apoptosis of epithelial cells, and factors released during cancer cell apoptosis can drive M₁-to-M₂ 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₁ cells driving MET and M₂ cells driving EMT, and epithelial cells driving M₁ maturation while mesenchymal cells driving M₂ maturation. With this new interaction, the system is now biased against epithelial cells, because (a) epithelial but not mesenchymal cells, are killed by M₁-macrophages, and (b) their dead counterparts may convert M₁ to M₂ cells that can, in turn, convert some epithelial cells to mesenchymal ones (EMT). In addition, the conversion from M₂ to M₁ 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).

FIGURE 3

Figure 3. Steady state solutions for a model including the M₁ to M₂ conversion assisted by apoptotic epithelial cancer cells. (A) Network of Model III. In this network, M₁ induces the apoptosis of epithelial cancer cells. And the factors released by apoptotic cancer cells enhance the M₁ to M₂ conversion. These interactions are highlighted by thicker lines. (B–E) Steady states of the epithelial (B), mesenchymal (C,D) and M₁ (E) populations are plotted as a function of η_me. (D) is a zoomed-in version of (C). Stable steady states are plotted in solid lines and unstable steady states are plotted in dotted lines. The model-specific parameters are: β = 1/36 h⁻¹, β_c = 1/1200 h⁻¹, η₁₂ = η₂₁ = 1/72 h⁻¹. When η_me is higher than a critical value ${\eta }_{\text{me}}^{\text{a}}$, the steady state E = M = 0 becomes stable (D). When η_me is higher than a critical value ${\eta }_{\text{me}}^{\text{c}}$, the steady state E = M = 0 is the only stable steady state of the system. (F) For a given λ_M, the critical value ${\eta }_{\text{me}}^{\text{c}}$ is plotted.

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₁ 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 ${\eta }_{\text{me}}^{\text{a}}=$0.0626 h⁻¹), 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⁻¹ and η_me = 0.0747 h⁻¹, 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₁ 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 ${\eta }_{\text{me}}^{\text{b}}=$0.0747 h⁻¹, the stable steady states in solid red lines disappears and the other two types of stable steady states co-existed (solid blue and black lines in Figures 3B–E). As η_me further increases and goes across ${\eta }_{\text{me}}^{\text{c}}=$0.1532 h⁻¹, 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⁻¹ 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⁻¹ (~0.0083 h⁻¹). 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 η₂₁) to get the corresponding ${\eta }_{\text{me}}^{\text{c}}$ 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 ${\eta }_{\text{me}}^{\text{c}}$ to the nominal experimental value 1/120 h⁻¹ (Figure 3F) whereas lowering down η_em or increasing η₂₁ 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 ${\eta }_{\text{me}}^{\text{c}}$ to a lower value (Figure S6) because of a potentially higher level of cancer-killing M₁-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₁ or M₂ macrophages relative to M₀ 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₁ 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₂ macrophage marker CD206 (MRC1). On the other hand, out of the 52 genes investigated, many genes upregulated in M₁ macrophages but downregulated in M₂ 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₁ macrophages, while the other state dominated by mesenchymal cells typically has higher number of M₂ macrophages.

TABLE 1

Table 1. Correlation coefficients of expression levels of specific genes with EMT scores, across many TCGA datasets.

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₂-like macrophages will be enriched in EMT-score-high tumors and whether genes that are expressed higher in M₁-like macrophages will be enriched in EMT-score-low tumors. Again, the lists of genes that are expressed higher in M₁- or M₂-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₂-like macrophages (but lower in M₁-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₂-like macrophages and mesenchymal cancer cells tend to stably co-exist. However, somewhat to our surprise, the genes that expressed higher in M₁-like macrophages (but lower in M₂-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₁-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.

TABLE 2

Table 2. M₂-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).

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–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–58).

Here, we focused on the interactions among macrophages of different polarizations (M₁ and M₂) 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₁ and M₂ 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₁-M₂ interconversion), with a given set of parameters, the system 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₂-dominated state to an epithelial- and M₁-dominated state.

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₂- into M₁-like (25), such as depletion of TAMs, reprogramming of TAMs toward M₁-Like macrophages, inhibition of circulating monocyte recruitment into tumor, etc. The polarization of M₂ 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₁-M₂ interconversion might not be effective to eliminate cancer cells; whereas in Model II increasing the M₂-to-M₁ conversion rate can help the system to switch to the epithelial cancer cell and M₁ 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 model-predicted stable co-existence of M₁-like macrophages and epithelial cancer cells (and likewise M₂-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₂-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₁-like macrophages in more mesenchymal tumors. There can be several reasons for the enrichment of M₁ genes in mesenchymal tumors: (i) the list of genes used here to identify the ones preferentially expressed in M₁-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₁-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, MET-inducing 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₁-like cells, the other dominated by mesenchymal and M₂-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).

Methods

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₁-like (M₂-like) macrophages. M₁-like macrophages will induce senescence of epithelial cancer cells and convert mesenchymal cancer cells to epithelial ones. On the other hand, M₂-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:

$\begin{array}{l}\frac{dE}{dt}={\lambda }_{E}E\left(1-\frac{E+M+M_1+M_2+M_0}{N_{max}}\right)\frac{1}{1+\alpha \frac{M_1}{M_1+K_1}}\\ +{\eta }_{me}M\frac{M_1^n}{M_1^n+\frac{M}{M+K_M^0}}-{\eta }_{em}E\frac{M_2^m}{M_2^m+\frac{E^k}{E^k+{\left(K_E^0\right)}^k}}\end{array}$

$\begin{array}{l}\frac{dM}{dt}={\lambda }_{M}M\left(1-\frac{E+M+M_1+M_2+M_0}{N_{max}}\right)\\ -{\eta }_{me}M\frac{M_1^n}{M_1^n+\frac{M}{M+K_M^0}}+{\eta }_{em}E\frac{M_2^m}{M_2^m+\frac{E^k}{E^k+{\left(K_E^0\right)}^k}}& (1)\end{array}$

$\begin{array}{l}\frac{d{M}_1}{dt}={\eta }_1\frac{E}{E+K_E}{M}_0\\ \frac{d{M}_2}{dt}={\eta }_2\frac{M}{M+K_M}{M}_0\\ \frac{d{M}_0}{dt}=-{\eta }_1\frac{E}{E+K_E}{M}_0-{\eta }_2\frac{M}{M+K_M}{M}_0\end{array}$

The description and the value of each parameter are given in Table 3.

TABLE 3

Table 3. Parameters used in our models.

In Model II, interconversions between M₁ and M₂ macrophages are included. Furthermore, the interconversions can be enhanced by corresponding cancer cells. The figure that illustrate this model is in Figure 2A. The equation to describe this system is as follows:

$\begin{array}{l}\frac{dE}{dt}={\lambda }_{E}E\left(1-\frac{E+M+M_1+M_2+M_0}{N_{max}}\right)\frac{1}{1+\alpha \frac{M_1}{M_1+K_1}}\\ +{\eta }_{me}M\frac{M_1^n}{M_1^n+\frac{M}{M+K_M^0}}-{\eta }_{em}E\frac{M_2^m}{M_2^m+\frac{E^k}{E^k+{\left(K_E^0\right)}^k}}\end{array}$

$\begin{array}{l}\frac{dM}{dt}={\lambda }_{M}M\left(1-\frac{E+M+M_1+M_2+M_0}{N_{max}}\right)\\ -{\eta }_{me}M\frac{M_1^n}{M_1^n+\frac{M}{M+K_M^0}}+{\eta }_{em}E\frac{M_2^m}{M_2^m+\frac{E^k}{E^k+{\left(K_E^0\right)}^k}}& (2)\end{array}$

$\begin{array}{l}\frac{d{M}_1}{dt}={\eta }_1\frac{E}{E+K_E}{M}_0-\left({\eta }_{12}^0+{\eta }_{12}\frac{M}{M+K_M}\right)M_1\\ +\left({\eta }_{21}^0+{\eta }_{21}\frac{E}{E+K_E}\right)M_2\end{array}$

$\begin{array}{l}\frac{d{M}_2}{dt}={\eta }_2\frac{M}{M+K_M}{M}_0+\left({\eta }_{12}^0+{\eta }_{12}\frac{M}{M+K_M}\right)M_1\\ -\left({\eta }_{21}^0+{\eta }_{21}\frac{E}{E+K_E}\right)M_2\end{array}$

$\frac{d{M}_0}{dt}=-{\eta }_1\frac{E}{E+K_E}{M}_0-{\eta }_2\frac{M}{M+K_M}{M}_0$

In the third model, additional interactions were introduced as illustrated in Figure 3A. M₁-like macrophages can induce apoptosis of epithelial cancer cells and factors released by apoptotic cancer cells can convert M₁-like macrophages into M₂-like macrophages. In order to restore the symmetry of the system, we further consider the therapeutic interaction: M₂-like macrophages can be re-polarized back to M₁-like macrophage by Type 1 T helper cells or IL-12. The equation to describe this system is as follows:

$\begin{array}{l}\frac{dE}{dt}={\lambda }_{E}E\left(1-\frac{E+M+M_1+M_2+M_0}{N_{max}}\right)\frac{1}{1+\alpha \frac{M_1}{M_1+K_1}}\\ +{\eta }_{me}M\frac{M_1^n}{M_1^n+\frac{M}{M+K_M^0}}-{\eta }_{em}E\frac{M_2^m}{M_2^m+\frac{E^k}{E^k+{\left(K_E^0\right)}^k}}\\ -\beta E\frac{M_1}{M_1+K_2}\end{array}$

$\begin{array}{l}\frac{dM}{dt}={\lambda }_{M}M\left(1-\frac{E+M+M_1+M_2+M_0}{N_{max}}\right)\\ -{\eta }_{me}M\frac{M_1^n}{M_1^n+\frac{M}{M+K_M^0}}+{\eta }_{em}E\frac{M_2^m}{M_2^m+\frac{E^k}{E^k+{\left(K_E^0\right)}^k}}& (3)\end{array}$

$\begin{array}{l}\frac{dC}{dt}=\beta E\frac{M_1}{M_1+K_2}-{\beta }_{c}C\\ \frac{d{M}_1}{dt}={\eta }_1\frac{E}{E+K_E}{M}_0-{\eta }_{12}\frac{C}{C+K_C}{M}_1+{\eta }_{21}{M}_2\\ \frac{d{M}_2}{dt}={\eta }_2\frac{M}{M+K_M}{M}_0+{\eta }_{12}\frac{C}{C+K_C}{M}_1-{\eta }_{21}{M}_2\\ \frac{d{M}_0}{dt}=-{\eta }_1\frac{E}{E+K_E}{M}_0-{\eta }_2\frac{M}{M+K_M}{M}_0\end{array}$

The description and the corresponding values of additional parameters are given in Table 3.

Correlation Analysis

For a fixed dataset, linear regression was performed for each gene of interest against the predicted EMT score (47). In each case samples were discretized into EMT-high or EMT-low based on median EMT score for hypothesis testing. The linear correlation coefficient was recorded in each case. We performed statistical analysis under the null hypothesis of zero correlation between gene expression fold-change and EMT score and recorded the corresponding p-values at significance level α = 0.05 using unpaired t-test. All datasets (colon adenocarcinoma, n = 286; lung adenocarcinoma, n = 525; prostate adenocarcinoma, n = 497; breast invasive carcinoma, n = 1,097) were obtained from the R2: Genomics Analysis and Visualization Platform (http://r2.amc.nl).

Author Contributions

XL, MKJ, KJP, and HL designed research. XL, MKJ, and JTG performed research. MKJ and JTG analyzed data. XL, MKJ, KJP, and HL wrote the paper.

Funding

This work was sponsored by the National Science Foundation NSF grant PHY-1427654 (Center for Theoretical Biological Physics). XL was supported by Stand Up to Cancer and The V Foundation. MKJ was also supported by a training fellowship from Gulf Coast Consortium as computational cancer biology training grant (CPRIT RP1705593). JTG was supported by National Cancer Institute of NIH (F30CA213878). KJP is supported by National Cancer Institute NCI grant CA093900. HL is also a Cancer Prevention and Research Institute of Texas Scholar in Cancer Research of the State of Texas at Rice University.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2019.00010/full#supplementary-material

References

1. Jain RK. Normalizing tumor microenvironment to treat cancer: bench to bedside to biomarkers. J Clin Oncol. (2013) 31:2205–18. doi: 10.1200/JCO.2012.46.3653

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Bocci F, Jolly MK, Tripathi SC, Aguilar M, Hanash SM, Levine H, et al. Numb prevents a complete epithelial-mesenchymal transition by modulating Notch signaling. J R Soc Interface (2017) 14:20170512. doi: 10.1098/rsif.2017.0512

CrossRef Full Text | Google Scholar

3. Peng DH, Ungewiss C, Tong P, Byers LA, Wang J, Canales JR, et al. ZEB1 induces LOXL2-mediated collagen stabilization and deposition in the extracellular matrix to drive lung cancer invasion and metastasis. Oncogene (2016) 36:1925–38. doi: 10.1038/onc.2016.358

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Wang T, Liu G, Wang R. The intercellular metabolic interplay between tumor and immune cells. Front Immunol. (2014) 5:358. doi: 10.3389/fimmu.2014.00358

PubMed Abstract | CrossRef Full Text | Google Scholar

5. del Pozo Martin Y, Park D, Ramachandran A, Ombrato L, Calvo F, Chakravarty P, et al. Mesenchymal cancer cell-stroma crosstalk promotes niche activation, epithelial reversion, and metastatic colonization. Cell Rep. (2015) 13:2456–69. doi: 10.1016/j.celrep.2015.11.025

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Nagarsheth N, Wicha MS, Zou W. Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy. Nat Rev Immunol. (2017) 17:559–72. doi: 10.1038/nri.2017.49

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Pitt JM, Marabelle A, Eggermont A, Soria JC, Kroemer G, Zitvogel L. Targeting the tumor microenvironment: removing obstruction to anticancer immune responses and immunotherapy. Ann Oncol. (2016) 27:1482–92. doi: 10.1093/annonc/mdw168

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Mantovani A, Bottazzi B, Colotta F, Sozzani S, Ruco L. The origin and function of tumor-associated macrophages. Immunol Today (1992) 13:265–70. doi: 10.1016/0167-5699(92)90008-U

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Noy R, Pollard JW. Tumor-associated macrophages: from mechanisms to therapy. Immunity (2014) 41:49–61. doi: 10.1016/j.immuni.2014.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Mantovani A, Schioppa T, Porta C, Allavena P, Sica A. Role of tumor-associated macrophages in tumor progression and invasion. Cancer Metastasis Rev. (2006) 25:315–22. doi: 10.1007/s10555-006-9001-7

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Allavena P, Sica A, Solinas G, Porta C, Mantovani A. The inflammatory micro-environment in tumor progression: The role of tumor-associated macrophages. Crit Rev Oncol Hematol. (2008) 66:1–9. doi: 10.1016/j.critrevonc.2007.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Solinas G, Germano G, Mantovani A, Allavena P. Tumor-associated macrophages (TAM) as major players of the cancer-related inflammation. J Leukoc Biol. (2009) 86:1065–73. doi: 10.1189/jlb.0609385

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Sica A, Mantovani A. Macrophage plasticity and polarization: in vivo veritas. J Clin Invest. (2012) 122:787–95. doi: 10.1172/JCI59643

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Mantovani A, Sozzani S, Locati M, Allavena P, Sica A. Macrophage polarization: tumor-associated macrophages as a paradigm for polarized M2 mononuclear phagocytes. Trends Immunol. (2002) 23:549–55. doi: 10.1016/S1471-4906(02)02302-5

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Edin S, Wikberg ML, Dahlin AM, Rutegård J, Öberg Å, Oldenborg PA, et al. The Distribution of macrophages with a M1 or M2 phenotype in relation to prognosis and the molecular characteristics of colorectal cancer. PLoS ONE (2012) 7:e47045. doi: 10.1371/journal.pone.0047045

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ohri CM, Shikotra A, Green RH, Waller DA, Bradding P. Macrophages within NSCLC tumour islets are predominantly of a cytotoxic M1 phenotype associated with extended survival. Eur Respir J. (2009) 33:118–26. doi: 10.1183/09031936.00065708

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Ma J, Liu L, Che G, Yu N, Dai F, You Z. The M1 form of tumor-associated macrophages in non-small cell lung cancer is positively associated with survival time. BMC Cancer (2010) 10:112. doi: 10.1186/1471-2407-10-112

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Murphy KM, Stockinger B. Effector T cell plasticity: flexibility in the face of changing circumstances. Nat Immunol. (2010) 11:674–80. doi: 10.1038/ni.1899

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Das A, Sinha M, Datta S, Abas M, Chaffee S, Sen CK, et al. Monocyte and macrophage plasticity in tissue repair and regeneration. Am J Pathol. (2015) 185:2596–606. doi: 10.1016/j.ajpath.2015.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Moghaddam AS, Mohammadian S, Vazini H, Taghadosi M, Esmaeili S-A, Mardani F, et al. Macrophage plasticity, polarization and function in health and disease. J Cell Physiol. (2018) 233:6425–40. doi: 10.1002/jcp.26429

CrossRef Full Text | Google Scholar

21. Guiducci C, Vicari AP, Sangaletti S, Trinchieri G, Colombo MP. Redirecting in vivo elicited tumor infiltrating macrophages and dendritic cells towards tumor rejection. Cancer Res. (2005) 65:3437–46. doi: 10.1158/0008-5472.CAN-04-4262

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Rolny C, Mazzone M, Tugues S, Laoui D, Johansson I, Coulon C, et al. HRG inhibits tumor growth and metastasis by inducing macrophage polarization and vessel normalization through downregulation of PlGF. Cancer Cell (2011) 19:31–44. doi: 10.1016/j.ccr.2010.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Xu M, Liu M, Du X, Li S, Li H, Li X, et al. Intratumoral delivery of IL-21 overcomes anti-Her2/Neu resistance through shifting tumor-associated macrophages from M2 to M1 phenotype. J Immunol. (2015) 194:4997–5006. doi: 10.4049/jimmunol.1402603

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Geeraerts X, Bolli E, Fendt SM, Van Ginderachter JA. Macrophage metabolism as therapeutic target for cancer, atherosclerosis, and obesity. Front Immunol. (2017) 8:289. doi: 10.3389/fimmu.2017.00289

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Genard G, Lucas S, Michiels C. Reprogramming of tumor-associated macrophages with anticancer therapies: radiotherapy versus chemo- and immunotherapies. Front Immunol. (2017) 8:828. doi: 10.3389/fimmu.2017.00828

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Nieto MA, Huang RY, Jackson RA, Thiery JP. EMT: 2016. Cell (2016) 166:21–45. doi: 10.1016/j.cell.2016.06.028

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Jolly MK, Ware KE, Gilja S, Somarelli JA, Levine H. EMT and MET: necessary or permissive for metastasis? Mol Oncol. (2017) 11:755–69. doi: 10.1002/1878-0261.12083

CrossRef Full Text | Google Scholar

28. Krebs AM, Mitschke J, Losada ML, Schmalhofer O, Boerries M, Busch H, et al. The EMT-activator Zeb1 is a key factor for cell plasticity and promotes metastasis in pancreatic cancer. Nat Cell Biol. (2017) 19:518–29. doi: 10.1038/ncb3513

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Siebzehnrubl FA, Silver DJ, Tugertimur B, Deleyrolle LP, Siebzehnrubl D, Sarkisian MR, et al. The ZEB1 pathway links glioblastoma initiation, invasion and chemoresistance. EMBO Mol Med. (2013) 5:1196–212. doi: 10.1002/emmm.201302827

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Jolly MK, Jia D, Boareto M, Mani SA, Pienta KJ, Ben-Jacob E, et al. Coupling the modules of EMT and stemness: A tunable “stemness window” model. Oncotarget (2015) 6:25161–74. doi: 10.18632/oncotarget.4629

CrossRef Full Text | Google Scholar

31. Huang RYJ, Wong MK, Tan TZ, Kuay KT, Ng AHC, Chung VY, et al. An EMT spectrum defines an anoikis-resistant and spheroidogenic intermediate mesenchymal state that is sensitive to e-cadherin restoration by a src-kinase inhibitor, saracatinib (AZD0530). Cell Death Dis. (2013) 4:e915. doi: 10.1038/cddis.2013.442

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Tripathi SC, Peters HL, Taguchi A, Katayama H, Wang H, Momin A, et al. Immunoproteasome deficiency is a feature of non-small cell lung cancer with a mesenchymal phenotype and is associated with a poor outcome. Proc Natl Acad Sci USA. (2016) 113:E1555–64. doi: 10.1073/pnas.1521812113

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Hollmén M, Roudnicky F, Karaman S, Detmar M. Characterization of macrophage–cancer cell crosstalk in estrogen receptor positive and triple-negative breast cancer. Sci Rep. (2015) 5:9188. doi: 10.1038/srep09188

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Yuan A, Hsiao YJ, Chen HY, Chen HW, Ho CC, Chen YY, et al. Opposite effects of M1 and M2 macrophage subtypes on lung cancer progression. Sci Rep. (2015) 5:14273. doi: 10.1038/srep14273

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Weigert A, Tzieply N, von Knethen A, Johann AM, Schmidt H, Geisslinger G, et al. Tumor cell apoptosis polarizes macrophages role of sphingosine-1-phosphate. Mol Biol Cell (2007) 18:3810–9. doi: 10.1091/mbc.E06-12-1096

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Sousa S, Brion R, Lintunen M, Kronqvist P, Sandholm J, Mönkkönen J, et al. Human breast cancer cells educate macrophages toward the M2 activation status. Breast Cancer Res. (2015) 17:101. doi: 10.1186/s13058-015-0621-0

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Su S, Liu Q, Chen J, Chen J, Chen F, He C, et al. A Positive feedback loop between mesenchymal-like cancer cells and macrophages is essential to breast cancer metastasis. Cancer Cell (2014) 25:605–20. doi: 10.1016/j.ccr.2014.03.021

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Yang M, Ma B, Shao H, Clark AM, Wells A. Macrophage phenotypic subtypes diametrically regulate epithelial-mesenchymal plasticity in breast cancer cells. BMC Cancer (2016) 16:419. doi: 10.1186/s12885-016-2411-1

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Mobius W, Laan L. Physical and mathematical modeling in experimental papers. Cell (2015) 163:1577–83. doi: 10.1016/j.cell.2015.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Gregory PA, Bracken CP, Smith E, Bert AG, Wright JA, Roslan S, et al. An autocrine TGF-beta/ZEB/miR-200 signaling network regulates establishment and maintenance of epithelial-mesenchymal transition. Mol Biol Cell (2011) 22:1686–98. doi: 10.1091/mbc.E11-02-0103

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Scheel C, Eaton EN, Li SHJ, Chaffer CL, Reinhardt F, Kah KJ, et al. Paracrine and autocrine signals induce and maintain mesenchymal and stem cell states in the breast. Cell (2011) 145:926–40. doi: 10.1016/j.cell.2011.04.029

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Schmalhofer O, Brabletz S, Brabletz T. E-cadherin, beta-catenin, and ZEB1 in malignant progression of cancer. Cancer Metastasis Rev. (2009) 28:151–66. doi: 10.1007/s10555-008-9179-y

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Jia D, Jolly MK, Harrison W, Boareto M, Ben-Jacob E, Levine H. Operating principles of tristable circuits regulating cellular differentiation. Phys Biol. (2017) 14:035007. doi: 10.1088/1478-3975/aa6f90

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Heusinkveld M, de Vos van Steenwijk PJ, Goedemans R, Ramwadhdoebe TH, Gorter A, Welters MJ, et al. M2 macrophages induced by prostaglandin E2 and IL-6 from cervical carcinoma are switched to activated M1 macrophages by CD4+ Th1 cells. J Immunol. (2011) 187:1157–65. doi: 10.4049/jimmunol.1100889

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wang Y, Lin Y-X, Qiao S-L, An H-W, Ma Y, Qiao Z-Y, et al. Polymeric nanoparticles promote macrophage reversal from M2 to M1 phenotypes in the tumor microenvironment. Biomaterials (2017) 112:153–63. doi: 10.1016/J.BIOMATERIALS.2016.09.034

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Frost JJ, Pienta KJ, Coffey DS. Symmetry and symmetry breaking in cancer: a foundational approach to the cancer problem. Oncotarget (2017) 9:11429–40doi: 10.18632/oncotarget.22939

PubMed Abstract | CrossRef Full Text | Google Scholar

47. George JT, Jolly MK, Xu S, Somarelli JA, Levine H. Survival outcomes in cancer patients predicted by a partial EMT gene expression scoring metric. Cancer Res. (2017) 77:6415–28. doi: 10.1158/0008-5472.CAN-16-3521

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Jablonski KA, Amici SA, Webb LM, Ruiz-Rosado JDD, Popovich PG, Partida-Sanchez S, et al. Novel markers to delineate murine M1 and M2 macrophages. PLoS ONE (2015) 10:e0145342. doi: 10.1371/journal.pone.0145342

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Whiteside TL. The tumor microenvironment and its role in promoting tumor growth. Oncogene (2008) 27:5904–12. doi: 10.1038/onc.2008.271

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Gajewski TF, Schreiber H, Fu YX. Innate and adaptive immune cells in the tumor microenvironment. Nat Immunol. (2013) 14:1014–22. doi: 10.1038/ni.2703

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Augsten M. Cancer-associated fibroblasts as another polarized cell type of the tumor microenvironment. Front Oncol. (2014) 4:62. doi: 10.3389/fonc.2014.00062

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Kuznetsov VA, Knott GD. Modeling tumor regrowth and immunotherapy. Math Comput Model (2001) 33:1275–87. doi: 10.1016/S0895-7177(00)00314-9

CrossRef Full Text | Google Scholar

53. Knútsdóttir H, Pálsson E, Edelstein-Keshet L. Mathematical model of macrophage-facilitated breast cancer cells invasion. J Theor Biol. (2014) 357:184–99. doi: 10.1016/j.jtbi.2014.04.031

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Knutsdottir H, Condeelis JS, Palsson E. 3-D individual cell based computational modeling of tumor cell–macrophage paracrine signaling mediated by EGF and CSF-1 gradients. Integr Biol. (2016) 8:104–19. doi: 10.1039/C5IB00201J

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Serre R, Benzekry S, Padovani L, Meille C, Andre N, Ciccolini J, et al. Mathematical modeling of cancer immunotherapy and its synergy with radiotherapy. Cancer Res. (2016) 76:4931–40. doi: 10.1158/0008-5472.CAN-15-3567

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Guo Y, Nie Q, MacLean AL, Li Y, Lei J, Li S. Multiscale modeling of inflammation-induced tumorigenesis reveals competing oncogenic and oncoprotective roles for inflammation. Cancer Res. (2017) 77:6429–41. doi: 10.1158/0008-5472.CAN-17-1662

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Mahlbacher G, Curtis LT, Lowengrub J, Frieboes HB. Mathematical modeling of tumor-associated macrophage interactions with the cancer microenvironment. J Immunother Cancer (2018) 6:10. doi: 10.1186/s40425-017-0313-7

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Michor F, Beal K. Improving cancer treatment via mathematical modeling: surmounting the challenges is worth the effort. Cell (2015) 163:1059–63. doi: 10.1016/j.cell.2015.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Choi SH, Kim AR, Nam JK, Kim JM, Kim J-Y, Seo HR, et al. Tumour-vasculature development via endothelial-to-mesenchymal transition after radiotherapy controls CD44v6+ cancer cell and macrophage polarization. Nat Commun. (2018) 9:5108. doi: 10.1038/s41467-018-07470-w

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Storer M, Mas A, Robert-Moreno A, Pecoraro M, Ortells MC, Di Giacomo V, et al. Senescence is a developmental mechanism that contributes to embryonic growth and patterning. Cell (2013) 155:1119–30. doi: 10.1016/j.cell.2013.10.041

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Jolly MK, Ward C, Eapen MS, Myers S, Hallgren O, Levine H, et al. Epithelial–mesenchymal transition, a spectrum of states: Role in lung development, homeostasis, and disease. Dev Dyn (2018) 247:346–58. doi: 10.1002/dvdy.24541

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Aras S, Raza Zaidi M. TAMeless traitors: macrophages in cancer progression and metastasis. Br J Cancer (2017) 117:1583–91. doi: 10.1038/bjc.2017.356

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Amend SR, Roy S, Brown JS, Pienta KJ. Ecological paradigms to understand the dynamics of metastasis. Cancer Lett. (2016) 380:237–42. doi: 10.1016/j.canlet.2015.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar