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

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

Figure 1 Collective phenotypes drive acid-mediated invasion. Spatial and temporal evolution of two distinct initial spatial configurations of identical numbers of cellular phenotypes leads to differential outcomes due to context-dependent selection. A low glycolysis phenotype (blue) and a high glycolysis phenotype (purple) compete for resources according to the rules outlined in section S1. Top row: a mixed configuration leads to no evolution. Acid-mediated invasion does not occur because the volumetric concentration of acid produced by aggressive cells is not enough to cause invasion when highly glycolytic cells are seeded far apart. Bottom row: In contrast, artificially placing the aggressive high glycolysis phenotypes on the rim leads to invasion from increased volumetric acid via a group-effect. Note: this figure has shorter timescales than subsequent figures, as it is seeded with pre-existing heterogeneity.

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

Figure 2 Defining metaphenotypes in the context of immune escape. (A) Six collective cellular metaphenotypes are defined as cancer cells with a given phenotype (e.g. PD-L1), microenvironmental condition (e.g. high acid or low glucose), or neighboring cell. Immune desert is the absence of recent immune interaction. (B) PD-L1 metaphenotypes depend on the likelihood of T-cell kill as a function of PD-L1 expression of self (PD-L1 Attack) or neighbor (Mooch PD-L1). (C) Acidification metaphenotypes depend on the rate of acidification contributed by self (Self-Acidify) or neighbors (Mooch acid). (D) The rate of acid-inactivation of T-cells. (E) Data from ref. 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 3 The effect of vasculature renewal and stability on tumor size and phenotype. (A) Hybrid discrete-continuum model grids. (B) Schematic of phenotypic trajectory of weak versus intermittent hypoxia vascular conditions. (C–E) N = 10 stochastic realizations are simulated, and the average tumor area (C), acid resistance phenotype (D), and glycolytic phenotype (E) across 10 values of stability (ν_mean∈ [0,100] days), and 10 values of renewal (p_ang∈ [0,1]). (F) An example of “weak vasculature” (ν_mean= 20; p_ang= 0.1). Acidic conditions in tumor core select for acid resistant and glycolytic Warburg phenotype. (G) An example realization of “intermittent hypoxia” (ν_mean= 20; p_ang= 0.9), where selection is limited because of adequate vascularization within the tumor core. See associated Supplementary Video S2.

Figure 4

Figure 4 Immune predation induces an evolutionary bottleneck. (A, B) Tumor area over time (left) and the number of T-cells for weak vasculature (A) and intermittent hypoxia vasculature (B) conditions), shown for no T-cells (green; α_T= 0), medium (blue-gray; α_T= 10⁻³) and high (purple; α_T= 10⁻²) immune response rates. (C–J) Example simulation stochastic realizations are shown across a range of immune response from low (top) to high (bottom). Immune predation tends to suppress tumor growth in weak vasculature conditions. In contrast, immune predation induces an evolutionary bottleneck for medium immune response rates (e.g. see F, H), causing aggressive tumor growth compared to the baseline of no immune response.

Figure 5

Figure 5 Evolution of metaphenotypes under treatment. Outcomes of tumor response and immune escape can be explained by observing the evolution of metaphenotypes under treatment with anti-PD-L1 (red) and buffer (blue), given in isolation or combination (purple). (A) Tumor area over time (weak vasculature) (B) growth rate over time (weak vasculature). (C) Tumor area over time (intermittent hypoxia vasculature) (D) growth rate over time (intermittent hypoxia vasculature). (E, F) Final distribution of metaphenotypes (exlcuding Immune Desert, see Supplementary Figure S7) at t = 300, repeated for weak vasculature (E) and intermittent hypoxia (F). (G, H) Muller plots showing the frequency of metaphenotypes over time in untreated and mono- or combination therapy, with snapshots of spatial configurations during and after treatment, with moderate immune predation (α_T = 10−²). See associated Supplementary Videos S3, S4 and Supplementary Figure S7. (I) Summary schematic. Each metaphenotype is ordered from most aggressive to least aggressive in facilitating acid-mediated invasion and tumor growth under immune predation. This interaction diagram describes the role of two treatments (anti-PD-L1, buffer) in promoting (green) or inhibiting (red) each metaphenotype. Metaphenotypes names are shown on the left, and defined mathematically in Box 2. Broadly, the two treatments offset one another by inhibiting the metaphenotypes that the opposite treatment promotes.

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.

Box 1. Defining metaphenotypes.
Let $T(x,y)$ be a two-dimensional grid representing the time since the last T-cell interaction has occurred within the local neighborhood of grid location $(x,y)$. We define the tumor-immune interaction grid, $I(x,y)$, where $I=1$ if an immune cell has traversed within a cancer cell’s neighborhood within the previous $T_w$ days and $I=0$ otherwise at the current timestep, $t$.

$\begin{array}{ll}I(x,y)=\{\begin{array}{ll}1,\hfill & \text{if} T(x,y)\ge t-T_w\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}& (1)\end{array}$

Metaphenotypes (MP) are defined in such a way that MP expression is scaled from zero to one and each cell can take on multiple MP: $\overrightarrow{M}=\left\{m_1,m_2,\dots ,m_7\right\}$where $m_i\in [0,1]$

2.1.1 MP1: Immune desert

We first consider the abscence of immune interaction: the immune desert metaphenotype, MP1, given by one minus I(x,y) given by Equation 1.

$\begin{array}{ll}\text{MP}1(x,y)=1-I(x,y)& (2)\end{array}$

2.1.2 MP2: PD-L1 attack

Next, we classify cells which employ the PD-L1 counter-attack, defined as high PD-L1 expression (low probability of T-cell kill; see Equation 17) with a recent T-cell interaction:

$\begin{array}{ll}\text{MP}2(x,y)=\underset{\text{Prob}.\text{ avoiding T}-\text{cell kill}}{\underbrace{(1-P_k)}}\times \underset{\text{recent T}-\text{cell interaction}}{\underbrace{I(x,y)}}& (3)\end{array}$

2.1.3 MP3: Mooch PD-L1

In contrast to MP2, cells which interact with T-cells but have low PD-L1 expression can rely on (“mooch”) neighboring cell protection. Here, the metaphenotype is proportional to neighborhood PD-L1 expression.

$\begin{array}{ll}\mathrm{MP3}(x,y)=\underset{\text{Prob}.\text{ T}-\text{cell kill}}{\underbrace{P_k}}\times \underset{\text{no} \text{T}-\text{cell interactions}}{\underbrace{1-I(x,y)}}\times \underset{j\in N_m}{\max }\text{PD}-\text{L}{1}_j& (4)\end{array}$

where $N_m$ is a Moore neighborhood of $N_m=8$ neighbors.

2.1.4 MP4: Self-acidify

As cell increase glycolytic capacity (phenotype value $p_G$), more protons are added. The per cell proton production rate is given by:

$\begin{array}{ll}{p}_R=f_H(1-B(t))& (5)\end{array}$

where proton production (see Box 2, Equation 14) rate is scaled by buffer treatment concentration, $B(t)$ (34, 35).

$\begin{array}{ll}\begin{array}{c}\mathrm{MP4}(x,y)=\hfill \\ \hspace{1em}\hspace{1em}\underset{\text{Probability of Acid}-\text{Inactivation}}{\underbrace{P_{AI}(x,y)}}\times \underset{\text{scaled proton production rate}}{\underbrace{{\overline{p}}_R}}\times \underset{\text{recent T-cell interaction}}{\underbrace{I(x,y)}}\end{array}& (6)\end{array}$

where the production rate, ${\stackrel{-}{p}}_R$, is normalized such that any value for phenotype above the buffering capability of a vessel is assumed to be mostly self-acidify metaphenotype (MP4), while below is assumed to be mostly mooch acid (MP5).

2.1.5 MP5: Mooch acid

Similarly, the mooch acidify metaphenotype occurs when the probability of T-cell acid-inactivation is high, but where the highly acidic microenvironment is not due to self-acidification.

$\begin{array}{ll}\begin{array}{c}\mathrm{MP5}(x,y)=\hfill \\ \hspace{1em}\hspace{1em}\underset{\text{Probability of Acid}-\text{Inactivation}}{\underbrace{P_{AI}(x,y)}}\times \underset{\text{scaled proton production rate}}{\underbrace{(1-{\overline{p}}_R)}}\times \underset{\text{recent T-cell interaction}}{\underbrace{I(x,y)}}\end{array}& (7)\end{array}$

This metaphenotype typically occurs early in simulations in empty regions without tumor or vasculature.

2.1.6 MP6: Proliferate fast

$\begin{array}{ll}\text{MP}6(x,y)=\left(1-\underset{\text{fraction of cell cycle completed}}{\underbrace{\frac{D_{x,y}}{T_m}}}\right)\times \underset{\text{recent T}-\text{cell interaction}}{\underbrace{I(x,y)}}& (8)\end{array}$

where $D_i$ is the time until next division for the cell at location $(x,y)$ and $T_m$ is the inter-mitotic cell division time for a metabolically normal cell.

2.1.7 MP7: Starve glucose

Tumor cells may also compete with T-cells to starve immune cells of glucose, giving rise to the following metaphenotype:

$\begin{array}{ll}\mathrm{MP7}(x,y)=\underset{\text{Prob. T}-\text{cell dies in low glucose}}{\underbrace{P_g}}\times \underset{\text{recent T}-\text{cell interaction}}{\underbrace{I(x,y)}}& (9)\end{array}$

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.

Table 1

Table 1 Model parameterization.

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.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Author contributions

JW: Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. FR: Formal analysis, Investigation, Methodology, Visualization, Writing – review & editing. CA: Formal analysis, Investigation, Methodology, Visualization, Writing – review & editing. RB: Formal analysis, Investigation, Methodology, Visualization, Writing – review & editing. KL: Conceptualization, Formal analysis, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing. MR-T: Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. AA: Conceptualization, Formal analysis, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. The authors gratefully acknowledge funding by the National Cancer Institute via the Cancer Systems Biology Consortium (CSBC) U54CA274507 and U01CA232382, the Physical Sciences Oncology Network (PSON) U54CA193489, and support from the Moffitt Center of Excellence for Evolutionary Therapy.

Conflict of interest

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

References

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

Google Scholar

2. Gatenby R, Gawlinski E, Gmitro A, Kaylor B, Gillies R. Acid-mediated tumor invasion: a multidisciplinary study. Cancer Res (2006) 66:5216.

PubMed Abstract | Google Scholar

3. Estrella V, Chen T, Lloyd M, Wojtkowiak J, Cornnell HH, Ibrahim-Hashim A, et al. Acidity generated by the tumor microenvironment drives local invasion. Cancer Res (2013) 73:1524–35.

PubMed Abstract | Google Scholar

4. Gillies RJ, Brown JS, Anderson AR, Gatenby RA. Eco-evolutionary causes and consequences of temporal changes in intratumoural blood flow. Nat Rev Cancer (2018) 18:576–85.

PubMed Abstract | Google Scholar

5. Silva AS, Yunes JA, Gillies RJ, Gatenby RA. The potential role of systemic buffers in reducing intratumoral extracellular ph and acid-mediated invasion. Cancer Res (2009) 69:2677–84.

PubMed Abstract | Google Scholar

6. Robertson-Tessi M, Gillies RJ, Gatenby RA, Anderson AR. 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

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Ibrahim-Hashim A, Robertson-Tessi M, Enriquez-Navas PM, Damaghi M, Balagurunathan Y, Wojtkowiak JW, et 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

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Damaghi M, West J, Robertson-Tessi M, Xu L, Ferrall-Fairbanks MC, Stewart PA, et 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

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Makowski L, Chaib M, Rathmell JC. Immunometabolism: From basic mechanisms to translation. Immunol Rev (2020) 295:5–14. doi: 10.1111/imr.12858

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Purohit V, Wagner A, Yosef N, Kuchroo VK. Systems-based approaches to study immunometabolism. Cell Mol Immunol (2022) 19(3):409–20. doi: 10.1038/s41423-021-00783-9

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

12. de Pillis L, Radunskaya A, Wiseman C. 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

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Mallet DG, De Pillis LG. A cellular automata model of tumor–immune system interactions. J Theor Biol (2006) 239:334–50. doi: 10.1016/j.jtbi.2005.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Robertson-Tessi M, El-Kareh A, Goriely A. A mathematical model of tumor–immune interactions. J Theor Biol (2012) 294:56–73. doi: 10.1016/j.jtbi.2011.10.027

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Robertson-Tessi M, El-Kareh A, Goriely A. 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

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Bottino D, Liu R, Bazzazi H, Venkatakrishnan K. Quantitative translation in immuno-oncology research and development. Clin Pharmacol Ther 108(3):430.

PubMed Abstract | Google Scholar

17. Mahlbacher GE, Reihmer KC, Frieboes HB. Mathematical modeling of tumor-immune cell interactions. J Theor Biol (2019) 469:47–60. doi: 10.1016/j.jtbi.2019.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Norton K-A, Gong C, Jamalian S, Popel AS. Multiscale agent-based and hybrid modeling of the tumor immune microenvironment. Processes (2019) 7:37. doi: 10.3390/pr7010037

PubMed Abstract | CrossRef Full Text | Google Scholar

19. El-Kenawi A, Gatenbee C, Robertson-Tessi M, Bravo R, Dhillon J, Balagurunathan Y, et 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

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Brown JA, Dorfman DM, Ma FR, Sullivan EL, Munoz O, Wood CR, et 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

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Rosenthal R, Swanton C, McGranahan N. 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

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Mascaux C, Angelova M, Vasaturo A, Beane J, Hijazi K, Anthoine G, et al. Immune evasion before tumour invasion in early lung squamous carcinogenesis. Nature (2019) 571:570–5. doi: 10.1038/s41586-019-1330-0

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Wu H, Estrella V, Beatty M, Abrahams D, El-Kenawi A, Russell S, et 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

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Pilon-Thomas S, Kodumudi KN, El-Kenawi AE, Russell S, Weber AM, Luddy K, et al. Neutralization of tumor acidity improves antitumor responses to immunotherapy. Cancer Res (2016) 76:1381–90. doi: 10.1158/0008-5472.CAN-15-1743

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Kumagai S, Koyama S, Itahashi K, Tanegashima T, Lin YT, Togashi Y, et al. Lactic acid promotes pd-1 expression in regulatory t cells in highly glycolytic tumor microenvironments. Cancer Cell (2022) 40(2):201–18.

PubMed Abstract | Google Scholar

26. Chang C-H, Qiu J, O’Sullivan D, Buck MD, Noguchi T, Curtis JD, et al. Metabolic competition in the tumor microenvironment is a driver of cancer progression. Cell (2015) 162:1229–41.

PubMed Abstract | Google Scholar

27. Watson MJ, Vignali PD, Mullett SJ, Overacre-Delgoffe AE, Peralta RM, Grebinoski S, et al. Metabolic support of tumour-infiltrating regulatory t cells by lactic acid. Nature (2021) 591:645–51.

PubMed Abstract | Google Scholar

28. Cappellesso F, Orban MP, Shirgaonkar N, Berardi E, Serneels J, Neveu MA, et 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

29. Ibrahim-Hashim A, Cornnell HH, Abrahams D, Lloyd M, Bui M, Gillies RJ, et al. Systemic buffers inhibit carcinogenesis in tramp mice. J Urol (2012) 188(2):624–31.

PubMed Abstract | Google Scholar

30. Robey IF, Baggett BK, Kirkpatrick ND, Roe DJ, Dosescu J, Sloane BF, et al. Bicarbonate increases tumor ph and inhibits spontaneous metastases. Cancer Res (2009) 69:2260–8.

PubMed Abstract | Google Scholar

31. Hashim AI, Cornnell HH, Coelho Ribeiro MDL, Abrahams D, Cunningham J, Lloyd M, et al. Reduction of metastasis using a non-volatile buffer. Clin Exp Metastasis (2011) 28:841–9. doi: 10.1007/s10585-011-9415-7

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Beckman RA, Yeang C-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 Full Text | Google Scholar

33. Jacobs SR, Herman CE, MacIver NJ, Wofford JA, Wieman HL, Hammen JJ, et al. Glucose uptake is limiting in t cell activation and requires cd28-mediated aktdependent and independent pathways. J Immunol (2008) 180:4476–86.

PubMed Abstract | Google Scholar

34. Mendler AN, Hu B, Prinz PU, Kreutz M, Gottfried E, Noessner E. Tumor lactic acidosis suppresses ctl function by inhibition of p38 and jnk/c-jun activation. Int J Cancer (2012) 131:633–40.

PubMed Abstract | Google Scholar

35. Rothstein TL, Mage M, Jones G, McHugh LL. Cytotoxic t lymphocyte sequential killing of immobilized allogeneic tumor target cells measured by time-lapse microcinematography. J Immunol (1978) 121:1652–6.

PubMed Abstract | Google Scholar

36. Bravo RR, Baratchart E, West J, Schenck RO, Miller AK, Gallaher J, et 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 Full Text | Google Scholar

37. Anderson AR. 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 Full Text | Google Scholar

38. Gatenbee CD, Schenck RO, Bravo RR, Anderson AR. 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 Full Text | Google Scholar

39. Nishino M, Ramaiya NH, Hatabu H, Hodi FS. Monitoring immune-checkpoint blockade: response evaluation and biomarker development. Nat Rev Clin Oncol (2017) 14:655–68. doi: 10.1038/nrclinonc.2017.88

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Miao Y, Yang H, Levorse J, Yuan S, Polak L, Sribour M, et al. Adaptive immune resistance emerges from tumor-initiating stem cells. Cell (2019) 177:1172–1186. doi: 10.1016/j.cell.2019.03.025

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Munn LL, Jain RK. Vascular regulation of antitumor immunity. Science (2019) 365:544–5. doi: 10.1126/science.aaw7875

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Gillies RJ, Ibrahim-Hashim A, Ordway B, Gatenby RA. Back to basic: Trials and tribulations of alkalizing agents in cancer. Front Oncol (2022) 12:981718.

PubMed Abstract | Google Scholar

43. Pilot C, Mahipal A, Gillies R. Buffer therapy—buffer diet. J Nutr Food Sci (2018) 8(02):684. doi: 10.4172/2155-9600.1000685

CrossRef Full Text | Google Scholar

44. Mpekris F, Voutouri C, Baish J.W, Duda D.G, Munn L.L, Stylianopoulos T, et al. Combining microenvironment normalization strategies to improve cancer immunotherapy. Proc Natl Acad Sci (2020) 117:3728–37. doi: 10.1073/pnas.1919764117

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Archetti M, Pienta KJ. Cooperation among cancer cells: applying game theory to cancer. Nat Rev Cancer (2019) 19:110–7. doi: 10.1038/s41568-018-0083-7

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

47. West J, Robertson-Tessi M, Anderson AR. Agent-based methods facilitate integrative science in cancer. Trends Cell Biol (2023) 33:300–11. doi: 10.1016/j.tcb.2022.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Alieva M, Wezenaar AK, Wehrens EJ, Rios AC. Bridging live-cell imaging and next-generation cancer treatment. Nat Rev Cancer (2023) 23:731–45. doi: 10.1038/s41568-023-00610-5

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Strelez C, Jiang HY, Mumenthaler SM. Organs-on-chips: A decade of innovation. Trends Biotechnol (2023) 41:278–80. doi: 10.1016/j.tibtech.2023.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Moghimi N, Hosseini SA, Dalan AB, Mohammadrezaei D, Goldman A, Kohandel M. Controlled tumor heterogeneity in a co-culture system by 3d bio-printed tumor-onchip model. Sci Rep (2023) 13:13648.

PubMed Abstract | Google Scholar

51. Nguyen DT, Famiglietti JE, Smolchek RA, Dupee Z, Diodati N, Pedro DI, et al. 3d in vitro platform for cell and explant culture in liquid-like solids. Cells (2022) 11:967.

PubMed Abstract | Google Scholar