ORIGINAL RESEARCH article
Computational Modeling of the Crosstalk Between Macrophage Polarization and Tumor Cell Plasticity in the Tumor Microenvironment
- 1Center for Theoretical Biological Physics, Rice University, Houston, TX, United States
- 2Department of Bioengineering, Rice University, Houston, TX, United States
- 3Medical Scientist Training Program, Baylor College of Medicine, Houston, TX, United States
- 4The James Buchanan Brady Urological Institute, Johns Hopkins University School of Medicine, Baltimore, MD, United States
- 5Department of Physics and Astronomy, Rice University, Houston, TX, United States
- 6Department of Physics, Northeastern University, Boston, MA, United States
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.
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 (M2) (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 (M1) (14). M1 and M2 macrophages have different roles during wound healing: while M1 macrophages initiate inflammatory responses, M2 macrophages contribute to tissue restoration (13). In the context of cancer, M1 macrophages have been generally considered anti-tumor (15–17), whereas M2 macrophages have been considered as pro-tumor (10).
However, macrophage polarization is not as rigid as the differentiation of T lymphocytes (18); instead, M1, M2, 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 M1-like macrophages (33). Forming a negative feedback loop, these M1-like macrophages can decrease the confluency of the cancer cells that polarized them (33). Moreover, pre-polarized M1 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 M1 macrophages into M2-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 M2-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 M1-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 (M1-like and M2-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.
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 M0, or macrophages M1 and M2), (b) EMT promoted by M2-like macrophages and MET promoted by M1-like macrophages, (c) polarization of monocytes (M0) to M1-like cells aided by epithelial cells, and that to M2-like cells aided by mesenchymal cells, (d) induction of senescence in epithelial cells by M1-like macrophages (Figure 1A). No inter-conversion among M1-like and M2-like macrophages or cell-death of macrophages is considered here in this model (hereafter referred to as “Model I”; see section Methods).
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 M2-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 M2-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 M2-promoted EMT or M1-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 (Mc) 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 M1-promoted MET. However, this resistance might not fully suppress MET, as MET still happens in the presence of M1 macrophages (38). Therefore, in Equation (1), we assume a Hill-like function to represent the resistance to M1-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 . 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 M2-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 M2-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 . 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 (Nmax, 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 (M1+M2+M0) to be a constant (Mc). Therefore, the maximum number of cancer cells will be Nmax -Mc.
In Model I, the final populations of E, M, M1, and M2 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 M2 macrophages in that state; since there are only three equations for four unknowns, M2 can be used as a parameter specifying (possibly discrete set of) solutions. As soon as M2/Nmax 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 M2-macrophage population. We identified that reducing the cooperativity of epithelial cancer cells in their resistance of M2-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 M2-promoted EMT or M1-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 M2/Nmax as the “control variable” is specific for Model I, because there is no interconversion between M1 and M2 macrophages. For a given set of parameters, the steady state level of M2 can vary continuously from 0 to Mc (Figures 1B–D), and in practice is determined by the initial conditions (initial number of E, M, M1, M2, and M0). For models in the following sections, M2/Nmax 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 M1 and M2, 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 M2) will first converge to one steady state in Model I (with the same M2/Nmax); on the longer time scale, as determined by the value of inter-conversion rate between M1 and M2, 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 M1- and M2-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 M1-like and M2-like macrophages, as reported in the literature (13, 35, 44, 45) (Figure 2A, see section Methods). We first assume constant interconversion between M1 and M2 (with rates denoted as and ) 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 M1/M2 ratio is simply determined by / according to Equation (2) when η12 = η21 = 0). In this model, we focused on increased cooperativity of M2-assisted EMT and M1-assisted MET (i.e., m = n = 2), because the interconversion between M1 and M2, 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 (Nmax, 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. Effects of interconversion between M1 and M2 cells, mediated by cancer cells. (A) Network of Model II. (B,C) Steady states of the epithelial (B) or M1 (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: η12 = η21 = 0 (for black lines) or 1/72 h−1 (for blue lines), m = 2, n = 2, k = 4, = 1/72 h−1. (D) In this plot, the key parameters are: η21 = 1/72 h−1 (for blue lines) or 1/36 h−1 (for green lines), m = 2, n = 2, k = 4, = = 1/72 h−1, and η12 = 1/72 h−1. In (B,D), the gray line is the maximum fraction of cancer cells (=0.7) as Mc 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 M1-assisted MET (ηem), lower M1 to M2 conversion rate (η12), or higher M2 to M1 conversion rate (η21) (Figure S4).
We next consider the case where interconversion between M1-like and M2-like macrophages is enhanced by cancer cells, i.e., mesenchymal cells enhance the M1-to-M2 transition, while epithelial cells enhance M2-to-M1 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 M1-M2 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 M2 cells makes the mesenchymal- and M2-dominated state more stable. Therefore, a higher ηme is required to compensate for the effects of low M1 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 M1 and M2, i.e., higher conversion rate from M2 to M1 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 M2 to M1. Note that the interactions in Model 2 are rather symmetric: the interconversion rates between M1 and M2 are either constants or enhanced by the corresponding cancer-cells. In the following section, we will consider the asymmetric case where only the M1 to M2 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., M1 macrophages can drive the apoptosis of epithelial cells, and factors released during cancer cell apoptosis can drive M1-to-M2 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., M1 cells driving MET and M2 cells driving EMT, and epithelial cells driving M1 maturation while mesenchymal cells driving M2 maturation. With this new interaction, the system is now biased against epithelial cells, because (a) epithelial but not mesenchymal cells, are killed by M1-macrophages, and (b) their dead counterparts may convert M1 to M2 cells that can, in turn, convert some epithelial cells to mesenchymal ones (EMT). In addition, the conversion from M2 to M1 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. Steady state solutions for a model including the M1 to M2 conversion assisted by apoptotic epithelial cancer cells. (A) Network of Model III. In this network, M1 induces the apoptosis of epithelial cancer cells. And the factors released by apoptotic cancer cells enhance the M1 to M2 conversion. These interactions are highlighted by thicker lines. (B–E) Steady states of the epithelial (B), mesenchymal (C,D) and M1 (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−1, βc = 1/1200 h−1, η12 = η21 = 1/72 h−1. When ηme is higher than a critical value , the steady state E = M = 0 becomes stable (D). When ηme is higher than a critical value , the steady state E = M = 0 is the only stable steady state of the system. (F) For a given λM, the critical value 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 M1 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 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 M1 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 0.0747 h−1, 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 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 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 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 to a lower value (Figure S6) because of a potentially higher level of cancer-killing M1-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 M1 or M2 macrophages relative to M0 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 M1 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 M2 macrophage marker CD206 (MRC1). On the other hand, out of the 52 genes investigated, many genes upregulated in M1 macrophages but downregulated in M2 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 M1 macrophages, while the other state dominated by mesenchymal cells typically has higher number of M2 macrophages.
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 M2-like macrophages will be enriched in EMT-score-high tumors and whether genes that are expressed higher in M1-like macrophages will be enriched in EMT-score-low tumors. Again, the lists of genes that are expressed higher in M1- or M2-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 M2-like macrophages (but lower in M1-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: M2-like macrophages and mesenchymal cancer cells tend to stably co-exist. However, somewhat to our surprise, the genes that expressed higher in M1-like macrophages (but lower in M2-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 M1-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. M2-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).
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 (M1 and M2) 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 M1 and M2 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 M1-M2 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 M2-dominated state to an epithelial- and M1-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 M2- into M1-like (25), such as depletion of TAMs, reprogramming of TAMs toward M1-Like macrophages, inhibition of circulating monocyte recruitment into tumor, etc. The polarization of M2 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 M1-M2 interconversion might not be effective to eliminate cancer cells; whereas in Model II increasing the M2-to-M1 conversion rate can help the system to switch to the epithelial cancer cell and M1 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 M1-like macrophages and epithelial cancer cells (and likewise M2-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 M2-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 M1-like macrophages in more mesenchymal tumors. There can be several reasons for the enrichment of M1 genes in mesenchymal tumors: (i) the list of genes used here to identify the ones preferentially expressed in M1-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 M1-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 M1-like cells, the other dominated by mesenchymal and M2-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).
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 M1-like (M2-like) macrophages. M1-like macrophages will induce senescence of epithelial cancer cells and convert mesenchymal cancer cells to epithelial ones. On the other hand, M2-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:
The description and the value of each parameter are given in Table 3.
In Model II, interconversions between M1 and M2 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:
In the third model, additional interactions were introduced as illustrated in Figure 3A. M1-like macrophages can induce apoptosis of epithelial cancer cells and factors released by apoptotic cancer cells can convert M1-like macrophages into M2-like macrophages. In order to restore the symmetry of the system, we further consider the therapeutic interaction: M2-like macrophages can be re-polarized back to M1-like macrophage by Type 1 T helper cells or IL-12. The equation to describe this system is as follows:
The description and the corresponding values of additional parameters are given in Table 3.
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).
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.
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2019.00010/full#supplementary-material
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Keywords: MET–mesenchymal-to-epithelial transition, EMT–epithelial-to-mesenchymal transition, M1-/M2-polarized macrophages, interaction network, multi-stability
Citation: Li X, Jolly MK, George JT, Pienta KJ and Levine H (2019) Computational Modeling of the Crosstalk Between Macrophage Polarization and Tumor Cell Plasticity in the Tumor Microenvironment. Front. Oncol. 9:10. doi: 10.3389/fonc.2019.00010
Received: 03 October 2018; Accepted: 03 January 2019;
Published: 23 January 2019.
Edited by:Ubaldo Emilio Martinez-Outschoorn, Thomas Jefferson University, United States
Reviewed by:Rafael Coelho Lopes De Sa, University of Massachusetts Amherst, United States
Virendra K. Chaudhri, Cincinnati Children's Hospital Medical Center, United States
Copyright © 2019 Li, Jolly, George, Pienta and Levine. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Herbert Levine, firstname.lastname@example.org
†These authors have contributed equally to this work
‡Present Address: Mohit Kumar Jolly, Centre for BioSystems Science and Engineering, Indian Institute of Science, Bangalore, India