Integrating mechanism-based T cell phenotypes into a model of tumor–immune cell interactions

The tumor microenvironment (TME) is a complex ecosystem comprised of many different cell types, including cancer cells and immune cells. It has long been shown that the immune system plays a pivotal role in tumor development.^1–4 Although immune cells can detect malignant cells as “not-self” and then eliminate them, tumors are still able to grow and escape the immune system. In fact, in several cancer types, the TME has been shown to strongly contribute to the immunosuppressive state of the tumor.^5,6 Furthermore, several studies point to the importance of the spatial location of stromal and immune cells in predicting prognosis and treatment response.^7,8 Thus, by considering the composition and spatial structure of tumors, it may be possible to better understand the factors that drive tumor progression and the efficacy of anti-cancer treatments.⁹ 

Given its complexity, efforts to extensively explore the TME using in vitro and in vivo models alone are intractable. Fortunately, computational modeling is a useful approach to study the TME and understand how individual cell–cell interactions and changes at the molecular and cellular levels influence tumor growth and response to treatment.¹⁰ For example, computational models have been used to better understand the interactions within the TME and predict how to improve immune-based cancer therapies.¹¹ A range of computational models have been developed to study various aspects of the TME and immunotherapy.^12–15 Agent-based models (ABMs) have proven to be particularly useful,^16–18 as they predict the spatial and temporal evolution of tumors and can simulate tumor behavior for a range of different conditions that are difficult and time-consuming to study using a purely experimental approach. Many ABMs exist and enable quantitative exploration of the TME. The level of detail included in the model depends on the specific question being addressed. For example, many ABMs use an entirely phenomenological framework to determine the cells' behaviors.^17,18 There are also examples of ABMs that embed networks inside of the agents to account for the role of intracellular features on cell behaviors.^19–21 Here, we aim to study how transcriptional regulation of CD8⁺ T cells influences the cells' state and properties such as killing ability and survival. Thus, we build on our prior work using biophysical agent-based modeling to simulate tumor growth^22–24 in a hybrid modeling framework by integrating our published Boolean model²⁵ of a gene regulatory network (GRN) with our existing rule-based ABM. The Boolean model simulates T cell evolution along a pseudo time trajectory providing valuable insight into transcriptional-level processes that influence emergent outcomes. This would not be possible with an ABM that only includes a rule-based approach to inform a cell's decision to proliferate, change phenotypic states, or induce cell killing. We predict the evolution of phenotypic T cell states, informed by the GRN encoded in the Boolean model. We consider three T cell states (naïve, pro-memory, and exhausted), which each have distinct properties. We apply the integrated model to study how CD8⁺ T cell state and properties influence tumor composition and spatial organization. Furthermore, we simulate immunotherapy via blockade of programmed cell death protein 1 (PD1), which is actively being pursued as a target to control tumor growth.^26,27 By accounting for the GRN governing T cell behavior, which influences cell–cell interactions, we provide a mechanistic view of tumor evolution and a framework for predicting how tumor growth can be controlled via genetic alterations of T cells. Our model enables exploration of how targeted genetic alterations affect the TME.

We developed an ABM that captures the emergent outcomes associated with the interactions between cancer cells, cytotoxic CD8⁺ T cells, CD4⁺ T cells (T helper and T regulatory cells), and macrophages classed into three phenotypes (M0, M1, and M2). Cytotoxic CD8⁺ T cells can be in a pro-memory state or an exhausted state that has a higher probability of dying and is less cytotoxic compared to pro-memory CD8⁺ T cells. The cells and their interactions are depicted in Fig. 1. We include the cells' impact on differentiation and state changes [Fig. 1(a)]. Cancer cells and CD4⁺ T cells influence macrophage differentiation from the naïve state. Cancer cells and M2 macrophages influence differentiation of CD4⁺ T cells into T regulatory cells. The model also accounts for how other cells influence the killing ability and proliferation of CD8⁺ T cells [Figs. 1(b) and 1(c)].

Here, we use a novel model of tumor growth that integrates a GRN of CD8⁺ T cell states (predicted by a Boolean model) into an ABM that predicts cell-level behaviors and cell–cell interactions. This model leverages our prior work in simulating T cell state transitions.²⁵ Additionally, the ABM represents a significant expansion of our published modeling of cell–cell interactions within the tumor microenvironment.^22–24 We apply the novel modeling framework to simulate tumor growth, with varying CD8⁺ T cell properties. We collect the two-dimensional position of cancer cells, T cells, and macrophages in an off-lattice setting, as well as cell-specific properties (radius, differentiation state, PDL-1 concentration in the case of cancer cells). We model TME evolution over a time period of 24 days starting from a small tumor mass located at the origin, and simulation data are logged once at the end of each simulated day. The model and its implementation are described in detail in the Methods section.

We utilized this model and considered the baseline values for the probabilities of CD8⁺ T cell death and cancer cell killing and investigated the effects of the rate of CD8⁺ T cell recruitment. Overall, the model predicts tumor escape and the presence of pro-tumor immune cells for these CD8⁺ T cell properties. Panel (a) of Figs. 2 and 3 shows the number of cancer cells over time averaged over ten simulation replicates with a baseline and high CD8⁺ T cell recruitment, respectively, as well as the standard deviation. The model predicts that higher CD8⁺ T cell recruitment leads to fewer cancer cells (26 472 and 19 591 for base and high recruitment rates, respectively, at the end of the model simulation).

As shown in Figs. 2 and 3, the final numbers of macrophages are predicted to be slightly lower for high CD8⁺ T cell recruitment (8169 and 6664 macrophages for base and high recruitment rate, respectively). Furthermore, we can analyze the proportion of different types of a given cell type at the end of the simulation. Specifically, we define the macrophage compartment as consisting of macrophages in their naïve state (M0) as well as differentiated states (M1 and M2). The T cell compartment consists of CD8⁺ T cells in their naïve, pro-memory, and exhausted states as well as CD4⁺ T cells as T-helper or T-regulatory cells. We find the final relative proportion of macrophage types is only moderately impacted by the CD8⁺ T cell recruitment rate (Fig. 4). Specifically, pro-tumor macrophages comprise 79% of the macrophage population for baseline T cell recruitment, compared to 84% for fast T cell recruitment.

Furthermore, the number of exhausted CD8⁺ T cells increases by nearly fourfold (from 668 to 2641 for base and high T cell recruitment, respectively). This corresponds to differences in the relative proportions of T cell subsets: 31% of CD8⁺ T cells are in an exhausted state for a base recruitment rate, compared to 63% for a high recruitment rate (Fig. 5). In addition, CD4⁺ T helper cells comprise 19% of the T cell compartment for baseline T cell recruitment, but only 7% for high T cell recruitment. These changes in the subsets of immune cells make sense, as faster CD8⁺ T cell recruitment enables more T cells to enter the tumor; however, those T rapidly evolve toward the exhausted state. A larger proportion of exhausted CD8⁺ T cells corresponds to fewer CD4⁺ helper cells, leading to fewer anti-tumor (M1) macrophages, due to T cell–macrophage interactions. We note that across both CD8⁺ T cell recruitment rates, naïve and pro-memory T cells represent less than 10% of the T cell compartment.

A useful feature of ABMs is their ability to predict the spatial organization of the cell populations being modeled. Thus, in addition to considering the time courses of the cell populations, we also examined the spatial layout at the end of the time course (Fig. 6) and over time (supplementary material, Movies 1 and 2). The predictions show that for both CD8⁺ T cell recruitment rates considered, there is minimal infiltration of CD8⁺ T cells, as they remain at the periphery of the tumor. This aligns with the concept of most tumors being classified as immunologically “cold,” defined as having low T cell infiltration.²⁸ The tumor periphery is also rich in macrophages. Additionally, there is moderate infiltration of macrophages into the tumor, and there is slightly less macrophage infiltration for the case of baseline CD8⁺ T cell recruitment rate. The model predicts that the average distance of macrophages from the tumor center is 1839 and 1544  $μ$m for base and high recruitment rate, respectively (Table S1). Note that shorter distance from tumor center indicates more infiltration.

The model predicts that the size of the cancer cell population is largely insensitive to the CD8⁺ T cell properties investigated. Here, we return to the analysis of the time courses shown in Fig. 2. The number of cancer cells increases exponentially for both low and high probability of CD8⁺ T cell death, as well as the two probabilities of T cell-mediated cancer cell killing [Figs. 2(b)–2(e), top row]. Furthermore, the number of cancer cells present at the end of the simulation is approximately the same across all conditions. Similarly, the spatial organization of the tumor is consistent across CD8⁺ T cell properties (supplementary material, Fig. S1).

In comparison, the number and relative proportions of macrophages depend on CD8⁺ T cell properties. The model predicts increases in the number of pro-tumor macrophages over time for all four conditions [Figs. 2(b)–2(e), middle row]. When CD8⁺ T cells have a low probability of cell-mediated cancer killing, M2 cells comprise 68%–80% of the macrophage compartment (Fig. 4). However, there are more M2 cells when the CD8⁺ T cells have a high probability of T cell-mediated cancer cell killing, comprising 83%–85% of macrophages, for high or low probability of CD8⁺ T cell death, respectively.

As anticipated, the probability of CD8⁺ T cell death and the probability of inducing cancer cell killing influence the number of T cells [Figs. 2(b)–2(e), bottom row] and relative proportion of T cells in the tumor (Fig. 5). When CD8⁺ T cells have a low probability of cell death, exhausted CD8⁺ T cells, Tregs, and CD4⁺ T helper cells make up approximately 62%, 21%–25%, and 9%–12% of the T cell population, respectively. In contrast, with high probability of CD8⁺ T cell death, the proportion of exhausted CD8⁺ T cell death is limited to 6% of all T cells. Having a smaller relative proportion of exhausted T cells corresponds to greater fractions of Tregs and T helper cells, ranging 61%–69% and 21%–28%, respectively, depending on the probability of T cell-mediated cancer cell killing. Again, these results can be traced to the cell–cell interactions included in the model. Having fewer pro-tumor macrophages (in the condition of high probability of T cell death) corresponds to more CD4⁺ T helper cells, which promotes differentiation of naïve macrophages into the anti-tumor state (M1).

Overall, the model predicts that with a baseline recruitment of CD8⁺ T cells, varying the properties of the CD8⁺ T cells influences the number and proportions of immune cells, with only a minimal effect on the number of cancer cells. This is not an intuitive result, demonstrating the utility of a systems-level model of the tumor microenvironment.

Across a range of values for the probability of CD8⁺ T cell death and cancer cell killing probability, the final number of cancer cells is lower for fast CD8⁺ T cell recruitment, compared to the baseline recruitment rate (Figs. 2 and 3, top row). One exception is the condition of a high probability of CD8⁺ T cell death combined with a low probability of T cell-mediated killing [panel (d) in Figs. 2 and 3], as compared to the baseline rate of CD8⁺ T cell recruitment.

The faster rate of CD8⁺ T cell recruitment is predicted to influence the absolute number and relative proportions of macrophage and T cell subsets in the tumor. For all four conditions, the number of pro-tumor macrophages is lower when T cell recruitment is high (Figs. 2 and 3, middle row). Additionally, anti-tumor macrophages comprise a larger fraction of the macrophage population (4%–10% and 3%–8% for high and base recruitment rates, respectively), as shown in Fig. 4. Furthermore, increasing CD8⁺ T cell recruitment leads to a larger number (Figs. 2 and 3, middle row) and a greater proportion (Fig. 5) of exhausted CD8⁺ T cells across all conditions. For example, exhausted CD8⁺ T cells comprise 31% for the base T cell recruitment rate, compared to 63% for the case of high recruitment. These changes correspond to a smaller fraction of Tregs and greater percentage of pro-memory CD8⁺ T cells. Altogether, the model predicts that higher CD8⁺ T cell recruitment increases the proportion of anti-tumor immune cells (M1 and pro-memory CD8⁺ T cells), concomitant with a decrease in the relative size of the anti-tumor T cell population (Tregs) and a nominal average decrease in anti-tumor macrophages (M2).

The spatial layout of the tumor further demonstrates the effect of increased CD8⁺ T cell recruitment (Fig. 7). A comparison of the spatial organization of the tumor for base and high recruitment rates for a representative model output shows fewer cancer cells, with the exception of the case for the high probability of CD8⁺ T cell death combined with low probability of T cell-mediated killing [Figs. 7(c) and 7(g)]. The tumors appear smaller, and we can visually appreciate the presence of substantially more exhausted CD8⁺ T cells (colored blue). Additionally, there is greater infiltration of immune cells into the tumor when CD8⁺ T cells are recruited at a faster rate and have a low probability of death [Figs. 7(e) and 7(g) and Table S1].

Finally, we considered the effects of altering the intracellular GRN of CD8⁺ T cells compared to the baseline model (termed wildtype, or “WT”). In particular, we simulated inhibition of PD1. Experimental studies demonstrate that blocking the PD1 pathway in naïve CD8⁺ T cells during differentiation can alter the cell state transitions. As in our previous work, we implement the PD1 blockade by removing the ability of NFATC1 to activate PD1 in the Boolean model of the CD8⁺ T cell GRN.

When CD8⁺ T cells are recruited to the tumor at the base rate, the model predicts that inhibiting PD1 strongly decreases the absolute number of exhausted CD8⁺ T cells, versus the WT case (compare the bottom rows of Figs. 2 and S1). Additionally, PD1 blockade significantly increases the proportion of pro-memory CD8⁺ T cells [compare Fig. 5(a) to Fig. S1(b)]. Despite altering the composition of the T cell population, when CD8⁺ T cell recruitment is at a baseline, a PD1 blockade does not affect the number of tumor cells or macrophages (compare the top two rows in Figs. 2 and S2). PD1 blockade also fails to alter the relative proportions of macrophages at the end of the model simulation [compare Fig. 4(a) to Fig. S1(a)].

In contrast, with higher CD8⁺ T cell recruitment, the model predicts that PD1 blockade can lead to smaller tumors compared to WT with high T cell recruitment. In particular, for the base and low values of the probability of CD8⁺ T cells considered, the numbers of tumor cells and pro-tumor macrophages are predicted to be lower than WT (compare the top two rows in Figs. 3 and 8). For example, with the baseline values for the probability of CD8⁺ T cell death and T cell-mediated cancer cell killing, the number of cancer cells is 19 591 and 17 036 for the WT and PD1 blockade, respectively. However, the number of cancer cells also depends on CD8⁺ T cell properties, even leading to more cancer cells with PD1 blockade. For example, when the CD8⁺ T cell death probability and killing probability are low, there are 18 245 and 20 867 cancer cells for WT and PD1 blockade, respectively. Interestingly, while the absolute numbers of M2 cells can decrease with PD1 blockade (compare the bottom row of Figs. 3 and 8), this does not always correspond to a change in the final relative proportion of macrophages [compare Figs. 4(b) and 9(a)].

In contrast, PD1 blockade consistently increases the fraction of pro-memory CD8⁺ T cells [compare Figs. 5(b) and 9(b)]. The condition of having a low probability of CD8⁺ T cell death and high probability of T cell-mediated killing is an exemplary case. Here, the percentages of T cells in the pro-memory state at the end of the model simulation are 6% and 40%, respectively, for the WT and PD1 blockade cases. Visual inspection of the spatial layouts also demonstrates the effect of PD1 blockade (Fig. S4). We note that for a high probability of CD8⁺ T cell death, PD1 blockade is not able to reduce the number of cancer cells [the top row of Figs. 8(d) and 8(e)] despite there being more pro-memory CD8⁺ T cells. In summary, a combination of high CD8⁺ T cell recruitment and low CD8⁺ T cell death is needed for PD1 blockade to effectively control the size of the cancer cell population.

We have presented an agent-based modeling framework that predicts the spatiotemporal behavior of the tumor-immune ecosystem comprised of cancer cells and distinct subsets of macrophages and T cells. The model highlights that CD8⁺ T cell recruitment, their death, and cancer cell-mediated killing can influence tumor growth. Specifically, the combination of low T cell death probability and high T cell-mediated cancer killing most strongly controls tumor size with high T cell recruitment. In addition to predicting how the number of cancer cells evolves with time, our results consider the absolute number and relative proportions of macrophage and T cell subsets, representing the immune component of the tumor. Experimental and clinical studies demonstrate the efficacy of targeting the immune checkpoint inhibitor, PD1, to inhibit tumor growth.²⁷ Therefore, we applied the model to predict the effects of altering gene regulation within CD8⁺ T cells to inhibit activation of the PD1 gene. This strategy is predicted to enhance the composition of the immune compartment, in agreement with published experimental studies. Given the mechanistic basis of our model, we can explain the efficacy of PD1 blockade, finding that it drastically shifts CD8⁺ T cells from an exhausted state to a pro-memory state. Importantly, the model shows that the efficacy of this PD1 blockade in controlling the number of cancer cells depends on CD8⁺ T cell properties. This result also aligns with experimental studies, where engineering CD8⁺ T cell characteristics is used to enhance the effects of immune checkpoint blockade.²⁹ 

A defining aspect of our work is the use of a Boolean model to inform the state transitions of CD8⁺ T cells. Specifically, we incorporate a Boolean representation of the GRN that mediates progression of CD8⁺ T cells from a naïve state to a terminally differentiated exhausted or pro-memory state.²⁵ The T cell's state influences its phenotype (probability of cell death) and cell–cell interactions (probability of promoting cancer cell killing). Thus, integration of the Boolean model allows the behaviors of CD8⁺ T cells in the ABM to be influenced by the intracellular gene regulation rather than taking a purely rule-based approach. The model also includes the effects of other cells on CD8⁺ T cell-mediated killing and death probability. We find that cell–cell interactions contribute to the final tumor composition (proportions of cancer and immune cells). With this modeling framework, we can simulate the effects of targeting interactions between cells. Furthermore, this integrated approach enables the investigation of targeted genetic strategies to alter a T cell's state. We have previously incorporated an intracellular network within an ABM to study the population-level effects of therapies that target macrophage intracellular signaling.²² Similarly, some ABMs consider intracellular models to inform the agent's behaviors.^19,30,31 However, most ABMs use phenomenological rules to determine cell state changes and phenotypes. Thus, having a mechanistic basis for the CD8⁺ T cell behavior helps advance the field of agent-based modeling.

We acknowledge some areas for future improvements of this work. First, our model includes three classes of macrophages and five types of T cells. This is a vast simplification of the tumor-immune ecosystem. Thus, we can expand the model to include other immune subsets. For example, recent work reveals that interactions between T cells and dendritic cells influence T cell state changes and tumor growth.^32,33 Related to this, we include differentiation of macrophages into just two broad classes (M1 and M2, anti- and pro-tumor, respectively); however, future work can consider a range of macrophage states and behaviors. Aligning the time steps simulated by the Boolean model is a second limitation of our work. We model that CD8⁺ T cells become exhausted within hours of stimulation by a cancer cell based on the findings of Rudloff and coworkers.³⁴ An issue with Boolean modeling is the use of pseudo-time, where the discrete time steps do not directly match real time. If the time required to reach an attractor state is known, one must assume what the interval between Boolean pseudo-time step corresponds to. For simplicity, we assume linear pseudo-time steps and that the Boolean model time steps correspond to a time interval of one hour. While the current work focuses on phenotypic features of CD8⁺ T cells, the effects of the time intervals within the Boolean model and aligning the timing between Boolean and ABM can be explored in future work. We can leverage the work of others in aligning timescales.^30,35 Third, we note that altering the GRN that mediates CD8⁺ T cell state transitions may not happen with 100% efficiency in a heterogeneous population of cells. We can account for this by sampling from 10⁴ simulated trajectories that are a mixture of WT and PD1 blockade. For example, assuming an alteration to the GRN within a population of cells is 80% effective, we would sample from a set of trajectories comprised of 2 × 10³ WT and 8 × 10³ PD1 blockade Boolean simulations. Finally, this work focuses on development of the integrated Boolean–ABM framework and explores how CD8⁺ T cell properties affect tumor growth. We have not performed a full analysis of the effect of varying all model parameters, a rigorous statistical analysis of model outputs, or a rigorous calibration of the model to tumor images from experimental models such as organoids or in vivo mouse tumors. Future work will leverage our recently developed machine learning-based approach to calibrate ABMs to such imaging data.²⁴ 

Taking an ecological view of a biological system by considering its spatial organization provides a framework for understanding the evolution of the system and the factors that control system behavior.³⁶ Similarly, in the context of cancer, by considering the composition and spatial structure of tumors, it may be possible to better understand the factors that drive tumor progression and response to treatment.⁹ We have demonstrated that biophysical agent-based modeling can produce biological insights into the tumor-immune ecosystem. Furthermore, accounting for the evolution of T cell states as determined by an intracellular model of CD8⁺ T cell gene activation simulates experimentally observed behaviors. Overall, our work provides a solid foundation for future studies to explore the role of cell–cell interactions in the TME and how they influence the response to treatment.

We develop an ABM to simulate the actions and interactions between cellular species in the TME. Specifically, we are able to capture the emergent outcomes associated with the interactions between cancer cells, cytotoxic CD8⁺ T cells (also referred to as CTLs), helper CD4⁺ T cells, and macrophages classed into three phenotypes (M0, M1, and M2). The cell–cell interactions are depicted in Fig. 1. This model provides a framework to develop sophisticated simulations of the TME and capture complex biophysical phenomena in simplified forms. We describe the model features below.

When simulating the model, it is computationally expensive to determine $Ni(t)$ for each cell, as the distances between each cell need to be calculated. To improve computational efficiency, we use a larger simulation time step of 1 h. This time step is both biologically realistic and computationally efficient. The reported time it takes for macrophages to respond to cytokines is on the order of hours.³⁸ However, to retain the accuracy of solving the force equation for every cell, we allow macrophage re-differentiation to occur every hour. At the beginning of each simulation step, we define a cell neighborhood for each cell. This is a list of cells with a $10xmax$ distance of cell i. We then determine $Ni(t)$ only from cells in that neighborhood and proceed to solve the above-mentioned equation at the small time step. This hastens computational time while preserving the accuracy of the force calculations.

To simplify proliferation and cell death, these processes are modeled as probabilities of occurring at each time point. Cell proliferation is also influenced by the physical presence of other cells, assuming that cell proliferation is inhibited by excess physical forces.²³ We model the excess physical stress as being due to the “overlap” between the radii of adjacent cells. If the sum of the total overlap with other cells is above a threshold, a cell's proliferation is inhibited until this overlap decreases. When a cell proliferates, the daughter cell is placed the distance of the cell radius away from the mother cell at a randomly selected angle. When a cell dies, it is removed from the simulation. Only cancer cells and CD8⁺ T cells can undergo proliferation.

Immune cells are recruited to the environment at a rate proportional to the number of cancer cells. They are recruited at a random angle around the tumor center, under the assumption that the area surrounding the tumor is well vascularized. When recruited, these cells are placed a random distance away from the tumor edge, sampled from a uniform distribution (Table S2).

In the model, immune cells migrate toward tumor cells, mimicking chemotaxis due to tumor-secreted chemokines.^39–41 Cells migrate with a specified amount of bias, which is modeled as the probability at each migration step of moving in the direction of the tumor center versus a random direction. This represents the impact of the extracellular matrix surrounding the tumor.^42–46 However, there are many barriers to immune cell infiltration into the tumor, leading to three general tumor-immune states: ignored (immune cells cannot detect the tumor), inflamed (immune cells can infiltrate the tumor), and excluded (immune cells are restricted to the outside of the tumor).⁴⁷ Some of these are due to interactions between immune cells, while others are due to physical barriers of the tumor. To capture this phenomenon in a simplistic manner with as few additional parameters as possible, we assume that both migration speed and bias decrease upon entry into the tumor. By adjusting these decreases, this model can capture different phenotypic tumor behaviors. Migration for all immune cell types is modeled in the same way; however, each cell type has its own bias and speed parameters. CD4⁺ and CD8⁺ T cells have a comparable migration speed and bias, whereas macrophages have a reduced speed and increased bias compared to the T cells (Table S2).

Within the TME, cells can communicate via diffusible factors, such as cytokines.⁴¹ These are often modeled using PDEs, with cells acting as point sources/sinks.^12,15,48 PDEs, however, require a great amount of computational time to solve, especially as environment size and the number of factors increases. Additionally, cells secrete several different cytokines, many with overlapping effects. This complexity makes it difficult to model the explicit biological effects of all the cytokines. Often, the effects of these factors are modeled either as a probability or a threshold value, instead of explicit modeling of downstream intracellular signaling. Because of this, researchers lump different cytokines into generic factors,⁴⁸ or replace diffusible factors with a distance threshold, where an effect occurs if one cell is within a specified distance from another cell.⁴⁹ 

In this model, cancer cells have three main actions: proliferation, death, and expression of the PD1 ligand (PD-L1). Proliferation and death occur as described above. Cancer cells gain PD-L1 at a probability proportional to the influence from CD4⁺ helper cells (via T cell-secreted IFN-γ, which promotes PD-L1 expression) and M2 macrophages,^12,52 dictated by the influence equation mentioned earlier. After proliferation, the daughter cells inherit the PD-L1 expression of the mother.

Classical in vitro characterization of macrophages leads to their classification into what have been called M2 pro-tumor and M1 anti-tumor macrophages. While it is now known that macrophages exist on a continuum, whereby their functions vary, and they can exhibit a range of phenotypes and functions,^53,54 we retain these broad M1 and M2 groupings for the purpose of exploring the effects of anti- and pro-tumor macrophages.

M2 cells promote the differentiation of CD4⁺ T cells into the regulatory state (via IL-10 and TGF-β) and promote M0 differentiation into the M2 state (via TGF-β).^56,58,59 We model all of these effects via the cell influence function described earlier. M2 cells also express PD-L1.

CD8⁺ T cells serve the main function of killing cancer cells. In the TME, we consider a CD8⁺ T cell to be in a naïve (N), exhausted (E), or pro-memory (M) state. CD8⁺ T cells can undergo cell death, kill a neighboring cancer cell, and migrate.⁶³ The probability of a CD8⁺ T cell dying or killing a cancer cell is influenced by (1) the CD8⁺ T cell state and (2) the cells in its neighborhood. Specifically, we set the model parameter for each of these cell behaviors relative to a baseline set of behaviors. CD8⁺ T cells in the N state are assumed to show no deviance from the base probabilities. In comparison, cells in the M state are half as likely to die and twice as likely to trigger cancer cell death, relative to baseline. Finally, CD8⁺ T cells in the E state maintain the base probability of death; however, their ability to promote cancer cell killing is reduced by an order of magnitude, relative to baseline. By varying cell death and cytotoxicity, we capture CD8⁺ T cell state-specific behaviors. The way in which the CD8⁺ T cell's intracellular GRN influences the probability of a CD8⁺ T cell dying or killing a cancer cell is detailed in the section “Integrating CD8⁺ T cell state decision into ABM.” Here, we describe how other cells in the TME influence CD8⁺ T cell behaviors.

The simulation loop proceeds as follows (Fig. 10):

1. Immune cells are recruited to the environment.

2. Cell neighborhoods are determined.

3. Influence from cell–cell interactions due to proximity (mimicking the effects of diffusible factors) is calculated, direct contact effects (CD8⁺ T cell-mediated cancer cell killing) are accounted for, CD8⁺ T cell states are updated based on intracellular Boolean model.

4. Forces between cells are solved.

5. Species proliferation and spontaneous cell death are modeled.

CD8⁺ T cell state decisions were determined directly from the output of our previously published Boolean model of the GRN that defines CD8⁺ T cell states.²⁵ This model predicts the expression profile of CD8⁺ T cell genes over pseudo-time. Although CD8⁺ T cells occupy a continuum of phenotypes, we broadly class the cell's state into one of three main phenotypes: naïve (N), pro-memory (M), and exhausted (E), based on the output of the GRN predicted by the Boolean model. The model logs the on/off expression of 24 genes, where each gene is classified as promoting a proliferative pro-memory (PP) state or an effector exhausted (EE) state. The Boolean model predicts the state of each CD8⁺ T cell at every time step based on the relative expression of PP and EE genes. Specifically, if more PP genes are on compared to EE, the phenotype is assumed to be M; the phenotype is assumed to be E if there are more EE genes on than PP. The first three pseudo-time steps are assumed to be in the N state; following stimulation of the T cell receptor, the predicted state in subsequent time steps is based on the expression of PP and EE genes.

For a single simulation, the output from a Boolean model simulation is given by an n × k matrix, where k represents the CD8⁺ T cell state (N, P, or E) given by the rules described above, and n represents the number of time steps. Thus, moving down the rows of this matrix represents the “trajectory” that a cell takes in its evolution from a naïve state to the terminally differentiated state. Once a cell reaches its terminal state, it stays in that state for the remaining n time steps. Our prior work demonstrated that the Boolean model captures both linear differentiation of a T cell from M to E, and circular differentiation where the cell oscillates between M and E before reaching the final state.^64,65

Ten thousand simulations of the Boolean model were performed, representing a population of 10⁴ distinct CD8⁺ T cells. We first consider the baseline or wildtype (WT). We also considered the case where the CD8⁺ T cell GRN encoded by the Boolean model reflects a PD1 blockade. Literature evidence shows that nuclear factor of activated T cells 1 (NFATC1) can promote antitumoral effector functions and memory CD8⁺ T cell differentiation through regulation of PD1 expression upon T cell activation.^66,67 Therefore, as described in our prior work,²⁵ we implemented indirect inhibition of PD1 via NFATC1. In the WT model, NFATC1 can activate PD1 (NFATC1 → PD1), and here, we removed NFATC1 from the model, inhibiting PD1 activation. With this alteration of the GRN, 10⁴ additional simulations of the Boolean model were performed to produce a population of 10⁴ distinct CD8⁺ T cells in which PD1 is inhibited.

The predicted CD8⁺ T cell state given by the Boolean model transitions through pseudo-time. It is straightforward to initialize every CD8⁺ T cell that is recruited to the tumor by sampling from the 10⁴ trajectories produced by the Boolean model and assigning the selected trajectory to the CD8⁺ T cell. However, we must assign a time interval for each pseudo-time step to align the transition of CD8⁺ T cells states with the timescale of the ABM. Experimental evidence indicates that CD8⁺ T cells become exhausted within hours following stimulation by a tumor cell.³⁴ We reflected this finding by updating the CD8⁺ T cell state every hour. This would involve an assumption that each pseudo-time step of the Boolean model corresponds to one hour. Thus, the predicted CD8⁺ T cell state given by the Boolean model at every time step is passed into the ABM. In the cases where the ABM time point is longer than the trajectory predicted by the Boolean (i.e., the CD8⁺ T cell has reached its terminally differentiated state before the current ABM simulation timepoint), we assume the T cell remains in its final state (the last value in its trajectory predicted by the Boolean model) for the rest of the simulation until cell death. For PD1 inhibition, we instead initialize and update the states of CD8⁺ T cells using the Boolean simulations corresponding to the altered GRN in which PD1 activation by NFATC1 is removed.

The “Model description” section details that CD8⁺ T cells can die and kill cancer cells. The probability of each of these behaviors is determined by the state predicted by the Boolean model, relative to a base probability. CD8⁺ T cells in the N state are assumed to show no deviance from the base probabilities. In comparison, cells in the M state are half as likely to die and twice as likely to trigger cancer cell death, relative to baseline. Finally, CD8⁺ T cells in the E state maintain the base probability of death; however, their ability to promote cancer cell killing is reduced by an order of magnitude, relative to baseline. By varying cell death and cytotoxicity, we capture CD8⁺ T cell state-specific behaviors.

This model is built in C++ using CMake for build processes. Scripts that develop the parameter grid as well as plotting scripts were implemented in Python. Two builds are presented, one developed for Apple Silicon and the other for a Unix-based computing cluster. The code is available at: https://github.com/FinleyLabUSC/Boolean-TME-ABM.

The parameter values used in the model are given in Table S2. We perform simulations to explore the effects of CD8⁺ T cell properties: recruitment rate (k_recr), probability of cell death (DP), and cancer cell killing probability (KP). We further investigate increased T cell recruitment (fivefold higher than the baseline value); varying the probability of cell death: fivefold increase and reduction in DP (high DP and low DP, respectively); and varying the probability of cancer cell killing: fivefold increase and reduction in KP (high KP and low KP, respectively). To account for variability in the predictions due to the probabilistic nature of the agent-based modeling, we performed ten replicates for each parameter combination considered.

See the supplementary material for supplementary figures and files. Movie 1 is representative of the time course of spatial layout over time for the baseline probabilities of CD8+ T cell death and cancer cell killing for base CD8+ T cell recruitment. Movie 2 is representative of the time course of spatial layout over time for the baseline probabilities of CD8+ T cell death and cancer cell killing for high CD8+ T cell recruitment. Table S1 provides the number and proportion of cell types, as well as mean infiltration distance for immune cells at the end of model simulations. Table S2 details the model parameters.

We acknowledge the Finley research group for critical evaluation of this manuscript. This work was supported by the USC Center for Computational Modeling of Cancer.

The authors have no conflicts to disclose.

Ethics approval is not required.

Neel Tangella: Formal analysis (lead); Investigation (lead); Methodology (lead); Software (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Colin G. Cess: Methodology (supporting); Software (equal); Writing – review & editing (supporting). Geena V. Ildefonso: Methodology (supporting); Writing – review & editing (supporting). Stacey D. Finley: Conceptualization (lead); Funding acquisition (lead); Resources (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available within the article and openly available on GitHub at https://github.com/FinleyLabUSC/Boolean-TME-ABM, Ref. [68].