Tumor-immune metaphenotypes orchestrate an evolutionary bottleneck that promotes metabolic transformation

Abstract

Introduction:

Metabolism plays a complex role in the evolution of cancerous tumors, including inducing a multifaceted effect on the immune system to aid immune escape. Immune escape is, by definition, a collective phenomenon by requiring the presence of two cell types interacting in close proximity: tumor and immune. The microenvironmental context of these interactions is influenced by the dynamic process of blood vessel growth and remodelling, creating heterogeneous patches of well-vascularized tumor or acidic niches.

Methods:

Here, we present a multiscale mathematical model that captures the phenotypic, vascular, microenvironmental, and spatial heterogeneity which shapes acid-mediated invasion and immune escape over a biologically-realistic time scale. The model explores several immune escape mechanisms such as i) acid inactivation of immune cells, ii) competition for glucose, and iii) inhibitory immune checkpoint receptor expression (PD-L1). We also explore the efficacy of anti-PD-L1 and sodium bicarbonate buffer agents for treatment. To aid in understanding immune escape as a collective cellular phenomenon, we define immune escape in the context of six collective phenotypes (termed “meta-phenotypes”): Self-Acidify, Mooch Acid, PD-L1 Attack, Mooch PD-L1, Proliferate Fast, and Starve Glucose.

Results:

Fomenting a stronger immune response leads to initial benefits (additional cytotoxicity), but this advantage is offset by increased cell turnover that leads to accelerated evolution and the emergence of aggressive phenotypes. This creates a bimodal therapy landscape: either the immune system should be maximized for complete cure, or kept in check to avoid rapid evolution of invasive cells. These constraints are dependent on heterogeneity in vascular context, microenvironmental acidification, and the strength of immune response.

Discussion:

This model helps to untangle the key constraints on evolutionary costs and benefits of three key phenotypic axes on tumor invasion and treatment: acid-resistance, glycolysis, and PD-L1 expression. The benefits of concomitant anti-PD-L1 and buffer treatments is a promising treatment strategy to limit the adverse effects of immune escape.

1 Introduction

Metabolism plays a complex but key role in the evolution of cancerous tumors. Localized hypoxia due to vascular instability and dysfunction leads to acidification of the tumor microenvironment. Decreased pH selects for acid-resistant tumor-cell phenotypes, followed by the emergence of aerobic glycolysis [i.e., the Warburg effect (1)]. Further microenvironmental acidification by these metabolically aggressive cells foments acid-mediated invasion (2–4). This nonlinear evolutionary trajectory through a range of metabolic phenotypes has been studied clinically, experimentally, and theoretically (5–8).

The effect of metabolic processes on the immune system is a multifaceted interaction between intracellular metabolism of many varied cell types with the surrounding microenvironment. Immunometabolism is a growing area of study (9) where systems biology and mathematical approaches are highly suited to studying tumorimmune dynamics (10–16), whether using non-spatial continuum approaches (17) or spatial agent-based models (18). However, very few tumor-immune models to date have incorporated the effects of cancer metabolism on immune function (19).

1.1 Metabolism and the tumor-immune response

Cytotoxic T lymphocytes (CTL, also known as α/β CD8+ effector T-cells) are key players in adaptive immune response which are activated via antigen presentation during the body’s initial inflammatory response and subsequently rapidly proliferate. Programmed cell death-1 (PD-1) is an inhibitory immune checkpoint receptor expressed on activated CTLs, and programmed cell death ligand-1 (PD-L1) is a cell surface marker that activates PD-1 signaling (20). Some cancers constitutively express PD-L1, leading to the development of anti-PD-1/PD-L1 therapy to counter this immune escape mechanism. Immune escape or evasion mechanisms may select for subclonal populations capable of withstanding immune predation (21), often well before tumor invasion into normal tissue (22).

We investigate two key connections between tumor metabolism and immune function: acid-inactivation and glucose competition. Acidic microenvironments have been shown to inactivate otherwise viable CTLs (23), as cells rescued from low pH environments had the ability to regain effector function (24). Tumor acidity also promotes regulatory T-cell (Treg) activity as well as an increase of PD-1 expression on Tregs, indicating that PD-1 blockade may increase suppressive capacity (25). Tumor-infiltrating CD8+ T-cells require glucose to support their killing function, hence competing for glucose with cancer cells dampens their anti-cancer response (26). In contrast, Tregs avoid competition for glucose through rewired metabolism away from aerobic glycolysis, which enhances their immune-suppression function within the tumor (27).

Acid-inactivation and glucose competition may diminish immunotherapy efficacy, suggesting a potential synergy between targeting intratumoral pH and immune checkpoint blockade. For example, combining oral bicarbonate buffering with immunotherapy (adoptive T-cell transfer, anti-CTLA4, or anti PD-1) increased responses in murine cancer models, presumably due to increased immune activity in a less acidic microenvironment (24). Another study showed that targeting bicarbonate transporters (e.g. SLC4A4) known to contribute to extracellular pH during progression of pancreatic adenocarcinomas (PDAC) (28) reduces tumor acidity, increases activation, cytotoxic activity, and perfusion of CD8+ T-cells, and sensitizes PDAC-bearing mice to immune checkpoint inhibition (28). Mechanistic modeling has been used to investigate the treatment effects of systemic pH buffers (sodium bicarbonate) to limit microenvironmental selection for acid-adapted phenotypes, resulting in significantly delayed carcinogenesis in TRAMP mice (7, 29). Buffers reduce intratumoral and peritumoral acidosis, inhibiting tumor growth (5) and reducing spontaneous and experimental metastases (30, 31).

1.2 The tumor-immune gambit

The back and forth of cancer treatment and a tumor’s evolutionary response has been compared to a chess match (32). Similarly, we show that immune predation of tumors can be likened to an “immune gambit”, where a temporary sacrifice of (normal glycolytic) cells on the periphery leads to long-term acceleration of the invasion of aggressive (highly glycolytic) phenotypes into surrounding tissue. Vascular dynamics are often abnormal in tumors whereby areas of poor vascularization are prone to develop acidic niches. We show that poor vascularization selects for aggressive phenotypes while high vascularization undergoes low levels of evolution. This phenomena has a Goldilocks effect, which occurs only under moderate levels of immune response. The immune gambit is described as a collective phenotype, which critically depends on the interplay between local vascularization, immune infiltration, and immune evasive phenotypes (PD-L1).

1.3 Collective cellular phenotypes: the “Metaphenotype”

In order to describe the collective nature of phenotypes operating within the context of surrounding cells and environmental conditions, we propose the concept of a “metaphenotype”. Each of these metaphenotypes account for phenotypic traits (e.g. PD-L1 expression or glycolytic rate) as well as surrounding environmental context (e.g. local glucose or pH concentration), and competition with neighboring cell types (immune, cancer, normal). A mathematical model is the ideal testing ground for defining collective phenotypes because it enables precise characterization of local context. A simple, contrived example in Figure 1 illustrates the need to quantify context-dependent selection in this model. This figure shows the time-evolution of identical phenotypic compositions that have varied initial spatial configurations (mixed or shell). The mixed configuration of low glycolysis (blue) and high glycolysis (purple) phenotypes leads to no evolution. The volumetric concentration of acid produced by aggressive cells is not enough to cause invasion when highly glycolytic cells are seeded far apart but artificially placing the aggressive high glycolysis phenotypes on the rim leads to invasion from increased volumetric acid via a group-effect. Clearly, both tumor phenotypic composition and neighboring context are important.

Figure 1

1.4 Mathematical modeling of immune metaphenotypes

Below, we propose and define six metaphenotypes in the context of immune escape and immunotherapy (see Figure 2). Then, we present a hybrid multiscale agent-based mathematical model that incorporates phenotypic, vascular, microenvironmental, and spatial heterogeneity to investigate the evolution of aerobic glycolysis in response to immune predation, over a biologically-realistic temporal scale (Figure 3). Next, we model immune predation by T-cells in the metabolically altered tumor microenvironment, including immune escape mechanisms such as acid-mediated inactivation of T-cells, T-cell inhibition by checkpoint ligand expression on tumor cells, and T-cell glucose deprivation (Figure 4). Finally, we quantify the evolution of metaphenotypes over time, illustrating the explanatory power of collective phenotypes in describing the response to buffer therapy and anti-PD-L1 in mono- and combination therapy (Figure 5).

Figure 2

33 (blue dots) was used to parameterize T-cell death rate in low glucose, shown in Equation 20. The Starve Glucose metaphenotype expression corresponds to low glucose concentrations.

Figure 3

Figure 4

Figure 5

2 Methods

2.1 Defining collective cellular phenotypes: immune metaphenotypes

First, we define six collective phenotypes (metaphenotypes) through the lens of immune escape (see Venn diagram in Figure 2A). Each metaphenotype is contingent on a recent tumor-immune interaction and defined in the context of local microenvironment, with the exception of a “null” metaphenotype: Immune Desert. The “null” metaphenotype is the lack of collective behavior: Immune Desert are cells that do not interact with T-cells. Next, we quantify two PD-L1 metaphenotypes: a counter-attack (tumor cell with high PD-L1 expression that has recently interacted with a T-cell; PD-L1 Attack, yellow), and a mooching PD-L1 (Mooch PD-L1, blue). As seen in Figure 2B, PD-L1 Attack is high in cells with high PD-L1 expression while Mooch PD-L1 is high in cells with low PD-L1 expression, but with neighbors that are high in PD-L1 Attack. See Box 1, Equations 3, 4. Two metaphenotypes rely on acid-inactivation: self-acidifying (highly glycolytic cells which secrete acid; Self-Acidify, pink) and non-producers (reside in acidic niche but do not produce acid; Mooch Acid, green). As seen in Figure 2C, Self-Acidify is high in cells with a high glycolytic phenotype, hence high acid production (see 5). In contrast, Mooch Acid cells have low glycolytic phenotype (not producing acid) but reside in highly a acidic niche that inactivates T-cells (Figure 2D). See Box 1, Equations 5–7. We also consider a proliferative phenotype that has recently divided into empty space (Proliferate Fast; red). See Box 1, Equation 8. Tumor cells also compete with immune cells for glucose (Starve Glucose; light blue). Figure 2E illustrates that Starve Glucose reside in areas with a high probability that T-cells die due to low glucose concentration. See Box 1, Equation 9. Importantly, each of these metaphenotypes (excluding Immune Dessert; see Equation 2) is contingent on a recent tumor-immune interaction, allowing us to track effective collective phenotypes: only metaphenotypes which survive an immune interaction.

2.2 Hybrid discrete-continuum multiscale model

We utilize this metaphenotype framework to better understand and predict tumor-immune interactions using a hybrid discrete- continuum multiscale model built using the Hybrid Automata Library framework (36). The mathematical model here is an extension of an experimentally validated multiscale model of cancer metabolism that incorporates the production of acid and acquired resistance to extracellular pH (6–8, 37) Figure 3A visualizes the model, which simulates a two-dimensional slice (panel A) through a tumor via a coupled cellular automata and partial differential equation model. A snapshot of multi-scale hybrid cellular automata model is shown (left-to-right) of the tumor spatial map, phenotypes, T-cells, diffusible molecules (oxygen, glucose, acid), PD-L1 and immune susceptibility. Values for parameterization are shown in Table 1. Values for parameters are typically identical to previous publications using the non-immune metabolism model (6, 7), except where parameter values are shown in brackets. In these cases, a parameter sweep is performed across the full range in order to determine the effect of the parameter on outcomes and test hypotheses. For convenience, we re-write the full model description, rules, and cell behaviors in Section S1.

3 Results

3.1 The effect of vasculature renewal and stability on tumor size and phenotype

In Figure 3, simulations are shown with the absence of immune predation to establish the model’s baseline dynamics, before quantifying immune predation in the next figure. The model tracks two tumor phenotypes: acid resistance and glycolysis (Figure 3B), which vary according to vascularization settings. The model contains two vascularization parameters: the rate of new vessel formation (vascular renewal) in hypoxic conditions and the average number of days before vessel collapse (vascular stability).

We compare two classifications of vasculature: weak vasculature (associated with low vessel stability and low rates of vessel renewal) and intermittent hypoxia (associated with low stability, but high renewal). Intermittent hypoxia has been shown to be an evolutionary driver of selection in tumors via environmental changes in glucose, oxygen, and acidity (see ref. 4 for a recent review).

Figures 3C–E show the average tumor area (C), and tumor phenotypes (D,E) for simulations across a range of vascular settings (no immune). Weak vasculature (low stability and renewal) typically results in larger tumors, more acid resistant phenotypes, and highly glycolytic phenotypes. Weak vasculature induces an acidic niche in the tumor core, selecting for acid-resistant phenotypes (blue). Increased turnover enables increased evolution and selection for aggressive Warburg phenotypes (pink), leading to acid-mediated invasion into surrounding normal tissue. Intermittent hypoxia (low vascular stability with high rates of renewal) generally leads to lower rates of selection and subsequently less invasion (Figure 3B).

Spatial maps of phenotypes are shown over time in Figures 3F, G along with a visualization called “phenotypic barcoding”, which visualizes the clone size, phenotype and lineage information over time (8) using the EvoFreq R package (38) (for more information on interpreting phenotypic barcoding plots, see Supplementary Figure S6). Figure 3F depicts the process by which weak vasculature selects for aggressive tumor growth. Acidic conditions in the tumor core (low glucose, low oxygen, and high pH) cause rapid death of glycolytically normal tumor cells with low levels of acid resistance. Selection for acid resistance occurs first (blue phenotypes), followed by selection for highly glycolytic tumor cells (pink phenotypes) which eventually invade into surrounding tissue. Conversely, in Figure 3G, intermittent hypoxia conditions result in little selection. The well-vascularized tumor core limits selection for aggressive phenotypes. This result underscores the importance of understanding the baseline vascular conditions before modeling the complex dynamics with the additional immune predation. A snapshot of the intratumoral oxygen, immune susceptibility (see Equation 21), phenotypes, and pH is shown at the end of each simulation.

3.2 Immune predation induces an evolutionary bottleneck

Figure 4 shows the response of two vascular conditions (weak and intermittent hypoxia) under no immune response (green; α_T= 0), medium (blue-gray; α_T= 10⁻³) and high (purple; α_T= 10⁻²) immune response rates. Immune cells are recruited in proportion to tumor size and response rate, α_T.

Immune response tends to suppress tumor growth in weak vasculature conditions (Figure 4A, left). Compared to baseline tumor growth, all levels of immune response result in greater tumor suppression. In contrast, immune predation in intermittent hypoxia conditions often leads to an initial response but fast regrowth (Figure 4B, left). This is confirmed by visual inspection of the phenotypic barcoding visualizations in Figures 4C–J. Weak vascular conditions select for aggressive phenotypes with little-to-no immune presence (Figure 4C). The addition of immune cells only serves to slow an already aggressive tumor (Figures 4E, G, I). In stark contrast, intermittent hypoxia conditions rarely select for strong growth in the absence of immune predation (Figure 4D). Immune predation serves as a selection pressure, in conditions where there would otherwise be very little selection.

Immune predation under intermittent hypoxia conditions induces an evolutionary bottleneck for medium immune response rates (e.g. see F, H), causing fast selection for aggressive growth compared to the baseline of no immune response. Interestingly, this effect occurs on a “Goldilocks” scale. The long neck of the bottleneck is associated with higher rates of tumor turnover (due to immune attack), selecting for phenotypes which are 1) inside an immune-evasive niche or 2) rapidly divide to outgrow immune kill.

Note: Figure 4 does not include immune escape mechanisms, which will be included in subsequent figures. The temporary bottleneck may be relevant to treatment with immune checkpoint inhibitors, enabling immune infiltration and predation of established tumors but leading to only a partial response (24, 39, 40).

3.3 Metaphenotypes explain immune escape under treatment

After establishing the baseline dynamics without (Figure 3) and with (Figure 4) immune predation, we next consider two treatments to mitigate immune escape and to reduce tumor growth: anti-PD-L1 and a pH buffer given in isolation or combination. A short window of treatment is simulated and results are compared to the untreated baseline. As seen in Figures 5A–D, combination therapy outperforms monotherapy in both vascular settings, but vascular dynamics drive differences in monotherapy outcomes. For example, anti-PD-L1 (red) therapy does not appreciably slow tumor evolution or growth in weak vasculature (Figures 5A, B). In contrast, anti-PD-L1 does induce large tumor remission in intermittent hypoxia (Figures 5C, D), albeit only temporarily before a strong relapse. These results are seen across a range of immune recruitment rates (Figures 5B, D).

The metaphenotypes leading to immune escape are shown in Figures 5E, F for each treatment scenario. As T-cell recruitment rate increases left-to-right, tumors evolve metaphenotypes in response to immune infiltration. Vascularization drives differential selection of metaphenotypes in baseline untreated dynamics. Weak vasculature (panel E; untreated) is associated with acidification metaphenotypes (Self-Acidify, pink; Mooch Acid, green). These are aggressive, highly glycolytic metaphenotypes that facilitate acid-mediated invasion. In contrast, intermittent hypoxia (panel F; untreated) selects for PD-L1-based immune-escape mechanisms (PD-L1 Attack, yellow; Mooch PD-L1, dark blue).

Treatment alters the type and magnitude of metaphenotype expression. Anti-PD-L1 selects for acidification metaphenotypes (Self-Acidify or Mooch Acid) in both vascularization cases. Buffer treatment eliminates the emergence of both Self-Acidify and Mooch Acid phenotypes by slowing evolution (e.g. refer to Figure 3B). But in response, PD-L1 Attack is selected (yellow). Only combination therapy targets both acidification metaphenotypes and PD-L1 phenotypes. Tracking the response of metaphenotypes to treatment explains why combination therapy is ideal for minimizing tumor growth, compared to monotherapy options. Importantly, only combination decreases the sum total of metaphenotypes expressed, and specifically targets aggressive phenotypes (Self-Acidify and Mooch Acid) across both vascularization scenarios.

3.4 Spatial configuration of metaphenotypes under treatment

The explanatory power of these defined metaphenotypes is seen most clearly by observing their spatial arrangement under high immune predation (see Figures 5G, H and associated Supplementary Videos S2, S3). For example, weak vasculature (Figure 5G) is associated with the Self-Acidify and PD-L1 Attack metaphenotypes on the invasive front of the tumor. Much of the tumor interior is unaffected by immune cells (Immune Desert), regardless of tumor phenotype. Treatment with Anti-PD-L1 selects for the aggressive Self-Acidify metaphenotype, while Buffer selects for PD-L1 Attack on the tumor rim. Combination therapy is required to achieve maximum tumor response, resulting in small tumors consisting mostly of non-aggressive metaphenotypes (Starve Glucose or Proliferate Fast).

In contrast, under intermittent hypoxia vasculature the Immune Desert comprises a much lower fraction of tumor metaphenotypes, as this improved vascularization delivers T-cells into the tumor core. PD-L1 Attack is used near blood vessels and on the tumor rim, and Self-Acidify does not occur due to low turnover in untreated conditions. Treatment with Anti-PD-L1 negates immune escape from PD-L1 Attack, inducing cellular turnover and subsequently selecting for Self-Acidify and Mooch Acid metaphenotypes. Combination therapy results in small, slow-growing tumors with less aggressive metaphenotypes (Mooch PD-L1 and Starve Glucose).

In both vasculature settings, cells slightly inset from the rim use metaphenotypes that Mooch Acid and Mooch PD-L1 from cells on the rim (see Supplementary Videos S2, S3) while cells in regions of high turnover employ the Proliferate Fast metaphenotype. Starve Glucose remains at low levels throughout all treatment modalities and vasculature settings. As seen in the (Supplementary Videos S2, S3), it is difficult to determine the major driver of immune escape from the maps of phenotypes alone, as areas of high glycolysis and high PD-L1 are each spatially heterogeneous and overlapping. There likely exists heterogeneity in vascular stability and renewal rates within a single patient’s tumor, which may drive heterogeneous metaphenotype expression (see Supplementary Video S5).

4 Discussion

Several factors contribute to a lack of responsiveness to immune checkpoint blockade, including abnormal tumor microenvironment where poor tumor perfusion hinders drug delivery and increases immunosuppression (41). Poor vascularization also leads to a hypoxic and therefore acidic microenvironment, increasing acid-mediated immunosuppression. The modeling we present here recapitulates this trend, as immune predation is less effective in weak vascularized tumors than in intermittently vascularized tumors. The importance of acidity in modulating immune response in cancer is only just beginning to be understood. Our results highlight the potential utility in buffering agents combined with immunotherapy. Whilst such buffering agents are not yet used in cancer treatment due to GI irritability and subsequent patient non-compliance, efforts continue to develop a buffer therapy that patients can tolerate and that is convenient to administer (42). Tumor acidosis can also be addressed by more indirect means. Some preclinical work has shown the potential influence of diet on acid buffering, but this remains poorly studied and may have limited effect on tumor pH (43).

Drugs that alter the vasculature are another possible indirect method for altering tumor pH. Development of agents that promote stable tumor vasculature would reduce acidosis and also increase both immune cell access and systemic-delivered drug penetration. However, there are potential risks to increasing tumor vascularity, regarding increased nutrient delivery and a higher potential for metastatic spread. Vascular renormalization can be enhanced through administration of anti-angiogenic agents (e.g., anti-vascular endothelial growth factor agents) to fortify immature blood vessels and improve tumor perfusion (44). However, our results indicate that administration of immune checkpoint blockade in tumors with increased vascularization may lead to a short-term good response but poor long-term outcomes as selection for increased glycolysis occurs. Mathematical modeling allows for direct comparison of initially identical simulations in the absence (Figure 3) and presence (Figure 4) of immune predation. We observe an immune gambit under high vascular renewal (intermittent hypoxia), due to an evolutionary bottleneck. The impact of this evolutionary bottleneck is reduced when anti-PD-L1 is combined with buffer therapy.

Characterization of collective phenotypes into metaphenotypes enables a straightforward explanation of the effect of treatment in a complex, multi-scale model. This characterization is necessary, in part, due to the fact that acid-mediated invasion is a collective phenotype phenomenon (Figure 1). Immune escape is also, by definition, a collective phenomenon by requiring the presence of two cell types in close proximity: tumor and immune. A summary schematic of the results is shown in Figure 5I. The interaction diagram describes the role of anti-PD-L1 and buffer in either promoting (green) or inhibiting (red) each metaphenotype. Broadly, the two treatments offset one another by inhibiting the metaphenotypes that the opposite treatment promotes. The two exceptions, starve glucose and immune desert, are both non-aggressive phenotypes. This summary schematic illustrates the utility of defining metaphenotypes in the context of treatment to provide insight into immune-escape dynamics. The most dominant mechanism of immune escape seen in the model is the lack of immune interactions (immune desert), especially when the tumor bed is poorly vascularized. Tumor-expressed PD-L1 is a viable immune-escape mechanism in the absence of treatment, across a range of vascularization, but treatment with anti-PD-L1 selects for the two acid-inactivation metaphenotypes (Self-Acidify and Mooch Acid). Environmental conditions must also consider neighboring (and self) cellular phenotypes. A cell in acidic conditions may rely on acid-inactivation either by self-production of acid or mooching from neighboring producer cells, a form of “public good” (45). Buffer therapy limits selection for self-acidification, driving selection toward less aggressive metaphenotypes (Glucose Starvation or Immune Desert). It’s also important to note that mooching metaphenotypes only occur in the presence of non-mooching phenotypes. Because of this, and the fact that phenotypes of individual cells change only slowly (upon division), mooching phenotypes are not expected to be a viable long-term immune escape strategy, but limited to transient, local patches co-localized with non-moochers. However, in a model where the ratio of two phenotypes is determined stochastically, for example, a population of both phenotypes could coexist for a longer period of time.

It would be of interest to test the predictions of this model, and as described above, previous experimental exploration of acid-mediated invasion and immune suppression aligns with the findings of our paper (3, 7, 19, 24). New technology has enabled spatially resolved transcriptomics which can quantify cellular heterogeneity in context of spatial information (46). However, the metaphenotype as a dynamic spatiotemporal metric is a challenge to measure (47), since most spatially resolved methods of interrogating a tumor in vivo are either destructive or of low resolution, and therefore lack the needed temporal component. Novel in vivo live-imaging technologies are one possible route to investigate how the different cellular phenotypes in the tumor environment change over time in response to emergent physiological changes and to therapeutic interventions (48). Another option would be the use of organoid cultures(e.g. 49), especially in conjunction with a 3D printer that could initialize different spatial configurations of the cells and environment for testing of hypotheses (50, 51).

The modeling framework presented here is not without its limitations. For example, 1) it is a two-dimensional representation of a three-dimensional tumors, 2) tumors may be heterogeneous in vasculature conditions, 3) we ignore directed motion of immune cells, and 4) parameters are an estimation based on literature values but may in reality be patient-specific. We also note limitations regarding the fundamental biological assumptions made in the model. For one, we are modeling the dominant form of cellular metabolism, namely the glycolytic and aerobic respiration pathway, but this is not the only source of cellular energy. Other forms involving glutamate, lactate, and more have been observed in tumor cells and would potentially alter how cells evolve in different environments. Another simplification involves the dynamics of immune response. T-cell activation, recruitment, engagement, and tolerance are all highly complex processes involving numerous cell types and cytokines, dynamic expression of different surface markers, and processes that work on many different timescales. Here we have limited ourselves to the influx of active T cells and not modeled upstream processes, nor immunosuppression generated by factors other than checkpoints and pH (e.g., regulatory T cells and other suppressive cells, TGF-beta and other suppressive secreted factors, and variable tumor antigenicity). These additional elements of immune response and tumor escape would certainly be worthy of further investigation in future work but will of course add significant complexity.

The intimate feedback between a growing tumor and the homeostatic tissue it is invading drives both ecological and evolutionary dynamics that should not be ignored in modern cancer therapy. The results we presented here indicate that treatments that modulate context may turn out to be just as important as those that target the tumor.

References

• 1

  WarburgO. Uber den stoffwechsel der karzinomezellen. Biochem Z (1924) 152:309–44.

  • Google Scholar
  • View reference in article

• 2

  GatenbyRGawlinskiEGmitroAKaylorBGilliesR. Acid-mediated tumor invasion: a multidisciplinary study. Cancer Res (2006) 66:5216.

  • Google Scholar
  • View reference in article

• 3

  EstrellaVChenTLloydMWojtkowiakJCornnellHHIbrahim-HashimAet al. Acidity generated by the tumor microenvironment drives local invasion. Cancer Res (2013) 73:1524–35.

  • Google Scholar
  • View reference in article

• 4

  GilliesRJBrownJSAndersonARGatenbyRA. Eco-evolutionary causes and consequences of temporal changes in intratumoural blood flow. Nat Rev Cancer (2018) 18:576–85.

  • Google Scholar
  • View reference in article

• 5

  SilvaASYunesJAGilliesRJGatenbyRA. The potential role of systemic buffers in reducing intratumoral extracellular ph and acid-mediated invasion. Cancer Res (2009) 69:2677–84.

  • Google Scholar
  • View reference in article

• 6

  Robertson-TessiMGilliesRJGatenbyRAAndersonAR. Impact of metabolic heterogeneity on tumor growth, invasion, and treatment outcomes. Cancer Res (2015) 75:1567–79. doi: 10.1158/0008-5472.CAN-14-1428

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  Ibrahim-HashimARobertson-TessiMEnriquez-NavasPMDamaghiMBalagurunathanYWojtkowiakJWet al. Defining cancer subpopulations by adaptive strategies rather than molecular properties provides novel insights into intratumoral evolution. Cancer Res (2017) 77:2242–54. doi: 10.1158/0008-5472.CAN-16-2844

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  DamaghiMWestJRobertson-TessiMXuLFerrall-FairbanksMCStewartPAet al. The harsh microenvironment in early breast cancer selects for a warburg phenotype. Proc Natl Acad Sci (2021) 118(3):e2011342118. doi: 10.1073/pnas.2011342118

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  MakowskiLChaibMRathmellJC. Immunometabolism: From basic mechanisms to translation. Immunol Rev (2020) 295:5–14. doi: 10.1111/imr.12858

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  PurohitVWagnerAYosefNKuchrooVK. Systems-based approaches to study immunometabolism. Cell Mol Immunol (2022) 19(3):409–20. doi: 10.1038/s41423-021-00783-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  KirschnerDPanettaJ. Modeling immunotherapy of the tumor-immune interaction. J Math Biol (1998) 37:235–52. doi: 10.1007/s002850050127

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  de PillisLRadunskayaAWisemanC. A validated mathematical model of cell-mediated immune response to tumor growth. Cancer Res (2005) 65:7950. doi: 10.1158/0008-5472.CAN-05-0564

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  MalletDGDe PillisLG. A cellular automata model of tumor–immune system interactions. J Theor Biol (2006) 239:334–50. doi: 10.1016/j.jtbi.2005.08.002

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  Robertson-TessiMEl-KarehAGorielyA. A mathematical model of tumor–immune interactions. J Theor Biol (2012) 294:56–73. doi: 10.1016/j.jtbi.2011.10.027

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  Robertson-TessiMEl-KarehAGorielyA. A model for effects of adaptive immunity on tumor response to chemotherapy and chemoimmunotherapy. J Theor Biol (2015) 380:569–84. doi: 10.1016/j.jtbi.2015.06.009

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  BottinoDLiuRBazzaziHVenkatakrishnanK. Quantitative translation in immuno-oncology research and development. Clin Pharmacol Ther108(3):430.

  • Google Scholar
  • View reference in article

• 17

  MahlbacherGEReihmerKCFrieboesHB. Mathematical modeling of tumor-immune cell interactions. J Theor Biol (2019) 469:47–60. doi: 10.1016/j.jtbi.2019.03.002

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  NortonK-AGongCJamalianSPopelAS. Multiscale agent-based and hybrid modeling of the tumor immune microenvironment. Processes (2019) 7:37. doi: 10.3390/pr7010037

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  El-KenawiAGatenbeeCRobertson-TessiMBravoRDhillonJBalagurunathanYet al. Acidity promotes tumour progression by altering macrophage phenotype in prostate cancer. Br J Cancer (2019) 121:556–66. doi: 10.1038/s41416-019-0542-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  BrownJADorfmanDMMaFRSullivanELMunozOWoodCRet al. Blockade of programmed death-1 ligands on dendritic cells enhances t cell activation and cytokine production. J Immunol (2003) 170:1257–66. doi: 10.4049/jimmunol.170.3.1257

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  RosenthalRSwantonCMcGranahanN. Understanding the impact of immune-mediated selection on lung cancer evolution. Br J Cancer (2021) 124(10):1615–7. doi: 10.1038/s41416-020-01232-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  MascauxCAngelovaMVasaturoABeaneJHijaziKAnthoineGet al. Immune evasion before tumour invasion in early lung squamous carcinogenesis. Nature (2019) 571:570–5. doi: 10.1038/s41586-019-1330-0

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  WuHEstrellaVBeattyMAbrahamsDEl-KenawiARussellSet al. T-cells produce acidic niches in lymph nodes to suppress their own effector functions. Nat Commun (2020) 11:1–13. doi: 10.1038/s41467-020-17756-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  Pilon-ThomasSKodumudiKNEl-KenawiAERussellSWeberAMLuddyKet al. Neutralization of tumor acidity improves antitumor responses to immunotherapy. Cancer Res (2016) 76:1381–90. doi: 10.1158/0008-5472.CAN-15-1743

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  KumagaiSKoyamaSItahashiKTanegashimaTLinYTTogashiYet al. Lactic acid promotes pd-1 expression in regulatory t cells in highly glycolytic tumor microenvironments. Cancer Cell (2022) 40(2):201–18.

  • Google Scholar
  • View reference in article

• 26

  ChangC-HQiuJO’SullivanDBuckMDNoguchiTCurtisJDet al. Metabolic competition in the tumor microenvironment is a driver of cancer progression. Cell (2015) 162:1229–41.

  • Google Scholar
  • View reference in article

• 27

  WatsonMJVignaliPDMullettSJOveracre-DelgoffeAEPeraltaRMGrebinoskiSet al. Metabolic support of tumour-infiltrating regulatory t cells by lactic acid. Nature (2021) 591:645–51.

  • Google Scholar

• 28

  CappellessoFOrbanMPShirgaonkarNBerardiESerneelsJNeveuMAet al. Targeting the bicarbonate transporter slc4a4 overcomes immunosuppression and immunotherapy resistance in pancreatic cancer. Nat Cancer 1–20 (2022) 3(12):1464–83.

  • Google Scholar
  • View reference in article

• 29

  Ibrahim-HashimACornnellHHAbrahamsDLloydMBuiMGilliesRJet al. Systemic buffers inhibit carcinogenesis in tramp mice. J Urol (2012) 188(2):624–31.

  • Google Scholar
  • View reference in article

• 30

  RobeyIFBaggettBKKirkpatrickNDRoeDJDosescuJSloaneBFet al. Bicarbonate increases tumor ph and inhibits spontaneous metastases. Cancer Res (2009) 69:2260–8.

  • Google Scholar
  • View reference in article

• 31

  HashimAICornnellHHCoelho RibeiroMDLAbrahamsDCunninghamJLloydMet al. Reduction of metastasis using a non-volatile buffer. Clin Exp Metastasis (2011) 28:841–9. doi: 10.1007/s10585-011-9415-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  BeckmanRAYeangC-H. Nonstandard personalized medicine strategies for cancer may lead to improved patient outcomes. Personalized Med (2014) 11:705–19. doi: 10.2217/pme.14.57

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  JacobsSRHermanCEMacIverNJWoffordJAWiemanHLHammenJJet al. Glucose uptake is limiting in t cell activation and requires cd28-mediated aktdependent and independent pathways. J Immunol (2008) 180:4476–86.

  • Google Scholar
  • View reference in article

• 34

  MendlerANHuBPrinzPUKreutzMGottfriedENoessnerE. Tumor lactic acidosis suppresses ctl function by inhibition of p38 and jnk/c-jun activation. Int J Cancer (2012) 131:633–40.

  • Google Scholar
  • View reference in article

• 35

  RothsteinTLMageMJonesGMcHughLL. Cytotoxic t lymphocyte sequential killing of immobilized allogeneic tumor target cells measured by time-lapse microcinematography. J Immunol (1978) 121:1652–6.

  • Google Scholar
  • View reference in article

• 36

  BravoRRBaratchartEWestJSchenckROMillerAKGallaherJet al. Hybrid automata library: A flexible platform for hybrid modeling with real-time visualization. PloS Comput Biol (2020) 16:1–28. doi: 10.1371/journal.pcbi.1007635

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  AndersonAR. A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion. Math Med Biology: A J IMA (2005) 22:163–86. doi: 10.1093/imammb/dqi005

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  GatenbeeCDSchenckROBravoRRAndersonAR. Evofreq: visualization of the evolutionary frequencies of sequence and model data. BMC Bioinf (2019) 20:1–4. doi: 10.1186/s12859-019-3173-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  NishinoMRamaiyaNHHatabuHHodiFS. Monitoring immune-checkpoint blockade: response evaluation and biomarker development. Nat Rev Clin Oncol (2017) 14:655–68. doi: 10.1038/nrclinonc.2017.88

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  MiaoYYangHLevorseJYuanSPolakLSribourMet al. Adaptive immune resistance emerges from tumor-initiating stem cells. Cell (2019) 177:1172–1186. doi: 10.1016/j.cell.2019.03.025

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  MunnLLJainRK. Vascular regulation of antitumor immunity. Science (2019) 365:544–5. doi: 10.1126/science.aaw7875

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  GilliesRJIbrahim-HashimAOrdwayBGatenbyRA. Back to basic: Trials and tribulations of alkalizing agents in cancer. Front Oncol (2022) 12:981718.

  • Google Scholar
  • View reference in article

• 43

  PilotCMahipalAGilliesR. Buffer therapy—buffer diet. J Nutr Food Sci (2018) 8(02):684. doi: 10.4172/2155-9600.1000685

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  MpekrisFVoutouriCBaishJ.WDudaD.GMunnL.LStylianopoulosTet al. Combining microenvironment normalization strategies to improve cancer immunotherapy. Proc Natl Acad Sci (2020) 117:3728–37. doi: 10.1073/pnas.1919764117

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  ArchettiMPientaKJ. Cooperation among cancer cells: applying game theory to cancer. Nat Rev Cancer (2019) 19:110–7. doi: 10.1038/s41568-018-0083-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  MarxV. Method of the year: spatially resolved transcriptomics. Nat Methods (2021) 18:9–14. doi: 10.1038/s41592-020-01033-y

  • CrossRef
  • Google Scholar

• 47

  WestJRobertson-TessiMAndersonAR. Agent-based methods facilitate integrative science in cancer. Trends Cell Biol (2023) 33:300–11. doi: 10.1016/j.tcb.2022.10.006

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  AlievaMWezenaarAKWehrensEJRiosAC. Bridging live-cell imaging and next-generation cancer treatment. Nat Rev Cancer (2023) 23:731–45. doi: 10.1038/s41568-023-00610-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  StrelezCJiangHYMumenthalerSM. Organs-on-chips: A decade of innovation. Trends Biotechnol (2023) 41:278–80. doi: 10.1016/j.tibtech.2023.01.004

  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  MoghimiNHosseiniSADalanABMohammadrezaeiDGoldmanAKohandelM. Controlled tumor heterogeneity in a co-culture system by 3d bio-printed tumor-onchip model. Sci Rep (2023) 13:13648.

  • Google Scholar
  • View reference in article

• 51

  NguyenDTFamigliettiJESmolchekRADupeeZDiodatiNPedroDIet al. 3d in vitro platform for cell and explant culture in liquid-like solids. Cells (2022) 11:967.

  • Google Scholar
  • View reference in article