## Abstract

Glioblastoma is the most common and deadly primary brain tumor in adults. All glioblastoma patients receiving standard-of-care surgery-radiotherapy-chemotherapy (i.e., temozolomide (TMZ)) recur, with an average survival time of only 15 months. New approaches to the treatment of glioblastoma, including immune checkpoint blockade and oncolytic viruses, offer the possibility of improving glioblastoma outcomes and have as such been under intense study. Unfortunately, these treatment modalities have thus far failed to achieve approval. Recently, in an attempt to bolster efficacy and improve patient outcomes, regimens combining chemotherapy and immune checkpoint inhibitors have been tested in trials. Unfortunately, these efforts have not resulted in significant increases to patient survival. To better understand the various factors impacting treatment outcomes of combined TMZ and immune checkpoint blockade, we developed a systems-level, computational model that describes the interplay between glioblastoma, immune, and stromal cells with this combination treatment. Initializing our model to spatial resection patient samples labeled using imaging mass cytometry, our model’s predictions show how the localization of glioblastoma cells, influence therapeutic success. We further validated these predictions in samples of brain metastases from patients given they generally respond better to checkpoint blockade compared with primary glioblastoma. Ultimately, our model provides novel insights into the mechanisms of therapeutic success of immune checkpoint inhibitors in brain tumors and delineates strategies to translate combination immunotherapy regimens more effectively into the clinic.

**SIGNIFICANCE STATEMENT** Extending survival times for glioblastoma patients remains a critical challenge. Although immunotherapies in combination with chemotherapy hold promise, clinical trials have not shown much success. Here, systems models calibrated to and validated against patient samples can improve preclinical and clinical studies by shedding light on the factors distinguishing responses/failures. By initializing our model with imaging mass cytometry visualization of patient samples, we elucidate how factors such as localization of glioblastoma cells and CD8+ T cell infiltration impact treatment outcomes.

## Introduction

Glioblastoma is a rare but lethal brain cancer (Alifieris and Trafalis, 2015). Standard of care for glioblastoma is maximal surgical resection, followed by radiation and chemotherapy with temozolomide (TMZ) (Stupp et al., 2005). If left untreated, most patients succumb to glioblastoma within three months after their initial diagnosis (Becker and Yu, 2012). Standard of care improves survival (Barker et al., 2012), but due to treatment resistance and despite its intensive treatment plan, the average survival time is only 15 months (Thakkar et al., 2014).

Immunotherapies such as immune checkpoint blockade (ICB) (McGranahan et al., 2019) have shown remarkable success in treating various cancers (Bausart et al., 2022). Immune checkpoints modulate the strength of the immune response and maintain self-tolerance (Ribas and Wolchok, 2018), but tumors exploit them to evade immune pressure. For example, programmed death ligand-1 (PD-L1) receptors on glioblastoma cells bind to checkpoint proteins (programmed death protein-1 (PD-1)) on CD8+ T cells to suppress the immune response (Pardoll, 2012; Han et al., 2020). ICB works by blocking PD-L1/PD-1 interaction, reactivating CD8+ T cell-mediated antitumoral immune responses. While successful against non-central nervous system cancers (Hazarika et al., 2017), recent clinical trials of ICB in glioblastoma have been disappointing (Lim et al., 2022).

The immune cycle in glioblastoma is a multi-step process that can be activated by combination therapies to overcome multifactorial immunosuppression (Bausart et al., 2022). For example, when in combination, chemotherapy improves cancer recognition by bolstering tumor antigenicity and immunogenic cell death, while ICB reduces immunosuppression by restoring T-cell activity (Leonetti et al., 2019). Hence, there has been an increased interest in combining TMZ and ICB (Yan et al., 2018). The CheckMate 548 trial (Bristol-Myers and Ono Pharmaceutical Co, 2020) explored the efficacy of nivolumab with standard of care. However, these trials did not substantially improve patient survival, highlighting the need to identify factors contributing to poor outcomes.

Many solid tumors exhibit significant inter/intratumor heterogeneity (Galon et al., 2006; Heindl et al., 2015), particularly in spatial distribution of cell types. The heterogeneous tumor microenvironment plays a significant role in determining cancer phenotypes by exerting selective pressure through a myriad of mechanisms (Yuan, 2016; Masud et al., 2022). With the advent of highly-multiplexed imaging technologies, including imaging mass cytometry (IMC), generating detailed information on the spatial configuration of the tumor became a reality (Chang et al., 2017; Elaldi et al., 2021). IMC is a labeling strategy that detects cellular biomarkers through magnetic labeling, providing simultaneous details on cellular phenotypes and spatial loci of cells.

Drug development involves lengthy clinical trials and substantial costs (Simoens and Huys, 2021). High failure rates, especially in oncology (Mohs and Greig, 2017), necessitate innovative approaches to guide drug development. Mathematical modeling is used in the development pipeline to offer insights into preclinical failures (Knight-Schrijver et al., 2016; Cassidy and Craig, 2019; Surendran et al., 2022). Previously, mathematical models have been employed to understand the effectiveness of mono/combination chemo/immunotherapies (Barazzuol et al., 2010; Lai and Friedman, 2017; Ayala-Hernández et al., 2021), as well as to infer the growth and treatment response of glioblastoma (Harpold et al., 2007; Böttcher et al., 2018; Massey et al., 2019). However, most of these are continuum models based on ordinary and partial differential equations and do not incorporate the role of spatial heterogeneity and localized cell-to-cell interactions within the tumor microenvironment.

Agent-based models (ABMs) are a computational modeling framework that consider spatial configurations of distinct cell types in the tumor (Ghaffarizadeh et al., 2018; Surendran et al., 2018; Cess and Finley, 2020). They describe realistic interactions between cells in tumors while accounting for the inherent stochasticity in cellular processes. Hence, ABMs are suitable for mapping the spatial organization of cancer, stromal, and immune cell components of glioblastomas to reflect the IMC output of patient resection data.

Building upon the previous success of ABMs in modeling glioblastoma treatments (Jenner et al., 2022), here we investigate the role of tumor spatial heterogeneity on combination chemo/immunotherapies. We construct an ABM of glioblastoma incorporating chemotherapy and immunotherapy and leverage IMC data of patient tumor resections to distinguish the spatial heterogeneity in glioblastoma tumors. By combining our data and model, we predict the influence of immune cells, particularly CD8+ T cells, as well as the impact of spatial heterogeneity on combination treatment success. Importantly, we validate our model’s predictions against IMC data from brain metastases from patients with various types of primary tumors, which generally respond better to ICB compared with glioblastoma. Together, our work highlights the role of quantitative systems pharmacology frameworks for understanding the physiologic basis behind drug efficacy and suggests new avenues to pursue experimentally and clinically to improve the use of immune checkpoint inhibitors to treat glioblastoma.

## Materials and Methods

#### Imaging Mass Cytometry Acquisition and Analyses

IMC on glioblastoma patient samples was performed as previously described (Karimi et al., 2023). Briefly, tissue microarrays containing formalin-fixed paraffin-embedded tissues were subject to deparaffinization and heat-mediated antigen retrieval according to the manufacturer’s instructions (Ventana Discovery Ultra auto-stainer, Roche Diagnostics). Slides were rinsed with 1X phosphate buffered saline and blocked for 45 minutes at room temperature (Dako Serum-free Protein Block Solution). Metal-conjugated antibodies were combined into a cocktail at the appropriate dilution (Dako Antibody Diluent) and applied to slides at 4°C overnight. Slides were washed with 0.2% Triton X and 1X phosphate buffered saline, and a secondary metal-conjugated anti-biotin was applied for 1 hour at room temperature. Slides were washed with 0.2% Triton X and 1X phosphate buffered saline and counterstained with Cell ID Intercalator-Ir as per manufacturer instructions (Fluidigm/Standard Biotools). Slides were rinsed with distilled water for 5 minutes and air dried prior to image acquisition (Hyperion Imaging System). IMC data used in this study are from Karimi et al., 2023, which details cell segmentation and lineage assignment. Briefly, cell segmentation was performed using a fully automated machine learning-based computer vision algorithm (Karimi et al., 2022). Lineage assignment was performed using a supervised approach based on canonical lineage markers.

#### Agent-Based Model of Glioblastoma

We developed an ABM of glioblastoma that incorporates imaging mass cytometry visualizations of tumor resections from untreated glioblastoma patients to initialize the tumor. To the best of our knowledge, this is the first time that patient glioblastoma tumor imaging mass cytometry data have been combined with an ABM. Since glioblastomas show significant intra- and inter-tumor heterogeneity, with diverse cell populations contained in different spatial niches (Comba et al., 2021), the use of our modeling framework allows for realistically simulating inter-patient variability in factors such as constituent cell types, their distributions and proportions, and spatial heterogeneity. The model developed here extends our recent work modeling oncolytic viral therapy in glioblastoma (Jenner et al., 2022) to include the administration of chemotherapy and immunotherapy. These extensions are detailed here, with a brief overview of the existing model.

Our agent-based model was created using the PhysiCell software framework (Ghaffarizadeh et al., 2018), an open-source platform that enables biologically accurate modeling of cellular processes, including cell cycling, cell death and cell-to-cell interactions, among others. A partial differential equation-based biotransport system, BioFVM (Ghaffarizadeh et al., 2016) is intertwined with PhysiCell to replicate diffusing substrates and cell signals within the tumor microenvironment. This biotransport system coupled with the agent-based model allows agents (cells) to dynamically update their phenotypes based on microenvironmental conditions. A detailed description of the PhysiCell platform is available in (Ghaffarizadeh et al., 2018) and (Jenner et al., 2022). In this work, we considered glioblastoma, macrophage, stromal, CD4+ and CD8+ T cells, their cell-to-cell interactions, and the effects of chemotherapy (temozolomide—TMZ) and immunotherapy (anti-PD1/PDL1 nivolumab). A schematic overview of the model is provided in Fig. 1.

To mimic biologically realistic logistic growth in the absence of TMZ, glioblastoma proliferation was modeled to depend on the local cell density, as in (Jenner et al., 2022). For this, we assumed that cells which experience a net mechanical pressure from their neighboring cells above a certain threshold () do not proliferate, whereas cells that experience overall mechanical pressure below this threshold have a proliferation rate β. In this way, our growth model emulated a carrying capacity density for tumor growth in which glioblastoma cell proliferation slows to a limiting value as more and more cells attain the threshold pressure and are unable to proliferate. As temozolomide is an alkylating agent that works by arresting the cell cycle to reduce cell division, eventually leading to apoptosis (Hotchkiss and Sampson, 2021), we assumed that its presence leads to reduced proliferation and increased death of glioblastoma cells.

To model the interactions between immune cells within the tumor, we included both innate (macrophages) and adaptive (CD4+ and CD8+ T cells) immune cell subtypes. Macrophages are known to eliminate dead cells through phagocytosis and can present antigens to T cells, thereby activating the latter cells (Hirayama et al., 2017). We therefore modeled macrophages as becoming activated through the phagocytosis of apoptosed glioblastoma cells (either through the effects of chemotherapy or through cytotoxic T cells, see below). Once activated, macrophages were modeled to activate any CD4+ or CD8+ T cells in their neighborhood, defined to be twice the sum of diameters of macrophage and the corresponding T cell. CD4+ T cells can also become activated through direct interactions with a tumor cell in their neighborhood. Once primed, CD4+ T cells secrete a chemokine which recruits CD8+ T cells through chemotaxis. These cytotoxic CD8+ T cells then target glioblastoma cells, causing their death through apoptosis (Jenner et al., 2022). We consider stroma as a composite of non-cancer and non-immune content of the tumor that holds tumor tissue together. Hence, for simplicity, we model them simply as taking up space and acting as a structural component of the tumor, with no movement or proliferation.

Cell movement in PhysiCell is characterized by a net locomotive force () acting on a cell which contributes a velocity to that cell (Ghaffarizadeh et al., 2018). The cell velocity is then defined as where is the cell’s movement speed, is the movement bias direction, and b is the magnitude of the movement bias. Here, is thought to be a combination of random and biased motion, where the bias magnitude varies from (corresponding to purely Brownian motion) to , which corresponds to deterministic motion along the direction of . We set for glioblastoma, macrophage, and CD4+ T cells because these cell types are modeled to undergo solely Brownian motion (i.e., no bias). Glioblastoma cells are known to move in a “stop and go” pattern, where the cell either moves or rests for a period of time (Gallaher et al., 2020). We follow the procedure described by (Gallaher et al., 2020) and (Jenner et al., 2022), in which cells are randomly assigned a migration status (stop or go) and the persistence times (time duration over which a cell either stays stationary or continues to move) are sampled from the distribution of persistence times (Gallaher et al., 2020). Cells that are assigned a “go’’ status move with Brownian motion. Once the cell exceeds its persistence time, a new migration status is randomly assigned. As mentioned above, the stromal cells are assumed to be stationary, implying that . The CD8+ T cells are biased to move up the chemokine gradient. Therefore, we fix and to be the gradient vector of the chemokine at that cell location. i.e., , where r is the location of the cell. With this choice, CD8+ T cells are designed to move up the chemokine gradient with some random motion (since the choice of ).

Glioblastomas are known to express negative immunomodulatory surface ligands (e.g., PD-L1) that lead to reduced anti-tumor immunity and to T cell anergy (Himes et al., 2021). These PD-L1 receptors expressed on glioblastoma cells bind to checkpoint protein (e.g., PD-1) present on CD8+ T cells to evade cytotoxic T cell targeting. Thus, we assumed that CD8+ T cell-induced apoptosis of glioblastoma cells can only occur in the presence of anti-PD-1 drug nivolumab which binds to PD-1 receptors and reactivates the immune response (see Treating glioblastoma using immune checkpoint blockade, below).

#### Pharmacokinetic/Pharmacodynamic Models of Temozolomide

To model TMZ concentrations at the tumor site over time, we integrated the population pharmacokinetic model developed by Ostermann et al.*.,* 2004 to predict the concentration of TMZ in the cerebrospinal fluid (CSF). Ostermann et al. developed a three-compartment model with first-order absorption from the gastrointestinal tract to model plasma and CSF concentrations (Fig. 2A) according to:
where , and are the amounts of TMZ in the absorption, plasma, and CSF compartments, respectively, is the absorption rate constant; is the rate constant from plasma to CSF; the rate constant from CSF to plasma; CL, oral clearance; , the volume of distribution in the central compartment; the volume of the distribution in the CSF, and and are the TMZ concentrations in the plasma and CSF, respectively.

Within the agent-based model we developed, variations in local TMZ concentrations are tracked by the BioFVM biotransport system (Ghaffarizadeh et al., 2016). Hence, the TMZ concentration at the tumor site ) is given by the equation:
where and are the diffusion coefficient and the decay rate of TMZ, respectively. Here is the Dirac delta function, **^{2}** is the k th cell’s position, is its volume and is its uptake rates. Eq. 7 is solved using first order implicit operating splitting on a discretized grid of voxels with length . We took the periphery of the tumor (denoted ) to be the set of points such that +

**20**)

^{2}, where R is the tumor radius. We imposed the boundary condition such that at , and 0 elsewhere.

The TMZ diffusion coefficient within the tumor can be estimated by using Stokes-Einstein equation (Valencia and González, 2011): where k is the Boltzmann constant, T is the absolute temperature, η is the solvent viscosity, and r is the solute hard sphere molecular ratio. The parameter n is an integer whose value (6, 4, or 2) is dependent of the volume relationship between the solvent and solute and we set . The molecular mass of TMZ is 194.151 g/mol (Ananta et al., 2016), which is equivalent to 0.194151 kDa. Converting this to a hydrodynamic radius gives 0.521 nm. The value of the Boltzmann constant is . The absolute temperature of the cerebellum in the human brain is 37.3°C (Kiyatkin, 2010) or . The dynamic viscosity of the tumor was estimated using that of the extracellular matrix of glioblastomas (), estimated by Santaguiliana et al. Thus, the diffusion coefficient of TMZ was taken to be (eq. 8). As TMZ has a half-life of 1.8–2 hours (Svec et al., 2018), we set the decay rate to be . The uptake rate of TMZ by glioblastoma tumor cells (U) was estimated using the TMZ brain tumor influx and efflux rates determined by Ballesta et al. (Ballesta et al., 2014). Assuming that the volume of the brain is 360 ml and using their estimate for the volume fraction of the interstitial compartment gives an uptake rate of /min.

TMZ induces cell death in proliferating cells by arresting the cell cycle. As such, to model the effect of TMZ and the induction of apoptosis in cells, we created a simple deterministic system describing the apoptotic effects of this interaction: where GBM is the number of glioblastoma cells, TMZ is the temozolomide concentration, D is the number of dead cells, and is the half effect of TMZ on the death of glioblastoma cells. To adapt this system of equations for modeling stochastic cell division and death in the ABM, we leveraged PhysiCell’s inbuilt “Live” and “Apoptosis” models (Ghaffarizadeh et al., 2018). The “Live” cell division model allows each cell to stochastically divide according to its division rate φ such that the probability of cell division in a given time interval, is given by (Ghaffarizadeh et al., 2018)

We set . Similarly, in the PhysiCell “Apoptosis” cell death model, the probability for a cell to die in the time interval is , where ψ represents the death rate. Hence, we set the death rate of glioblastoma cells as As mentioned previously, tumor cells in the absence of TMZ proliferate at the exponential rate β. (Saha et al., 2020) computed dose-responses of the glioblastoma stem cells by measuring the viability of these cells 4 days after an initial dose of TMZ. Fitting the standard Emax effects curve, i.e., to their measurements gave an estimate of using the molecular mass of TMZ (Supplemental Fig. 1).

Finally, we considered the dosing schedule where of TMZ was administered for 5 days, consistent with standard-of-care for glioblastoma (Ostermann et al., 2004; Stupp et al., 2005). Considering an average adult male has a body surface area of , this schedule corresponds to of TMZ per day for 5 days (Fig. 2B).

#### Treating Glioblastoma Using Immune Checkpoint Blockade

The approved flat dosage regimen for the anti-PD-1 drug nivolumab is 240 mg every two weeks. A study by Lee et al. (Lee et al., 2018) established that the flat dosage results in similar exposure to regimens dosing 3 mg/kg once every two weeks. To model the concentration of anti-PD-1 drug at the tumor site, we leveraged the method developed by Storey et al., 2020 who used pharmacokinetic data from the phase 1 study (Brahmer et al., 2010) to derive a linear relationship between the dose of anti-PD-1 drug (, in mg/kg) to plasma concentrations () (Supplemental Fig. 2). This relationship is given by,

Converting the units of to using the molecular mass of nivolumab as 143599.39 g/mol (Chaudhari, 2017), we get

The change in concentration of the anti-PD-1 immunotherapy drug at the periphery of the tumor , denoted by APD, is represented by the equation: where represents the source of the anti-PD-1 drug as a function of time. As mentioned above, the standard administration schedule of nivolumab is a single intravenous dose of 3 mg/kg, administered for one hour every two weeks. Computing for this dose gives . For simplicity, we assumed this value as the baseline estimate for during the one hour of anti-PD-1 drug administration (Storey et al., 2020). This implies that for each time (in the unit of minutes) at which anti-PD-1 is administered,

Variations in the local concentration of the anti-PD-1 drug through diffusion and decay are tracked by the BioFVM partial differential equation,

Similar to eq. 7, we used the boundary condition and initial condition and to solve eq. 18.

Using the Stokes-Einstein equation given in eq. 8 and taking the molecular mass of nivolumab as 143599.39 g/mol (Chaudhari, 2017), we computed the diffusion coefficient of the immune check point inhibitor to be As nivolumab has a half-life of 15 days (Brahmer et al., 2010), we estimated the decay rate to be

We assumed that CD8+ T cells induce apoptosis in glioblastoma cells only when the number of ICB-bound PD-1 receptors on CD8+ T cells is above the threshold . For each CD8+ T cell, we tracked the number of ICB-bound PD-1 receptors ( and the number of the unbound PD-1 receptor . The immune checkpoint inhibitor and PD-1 receptor binding kinetics were modeled by the system of equations: where and are the binding and unbinding rates and is the amount of anti-PD-1 drug present at the location of CD8+ T cell. We converted the concentration of ICB in the unit of to by multiplying it by the Avogadro number . We set the initial number of unbound PD-1 receptors on CD8+ T cells as 3096, which is found to be the average number of PD-1 receptors expressed on CD8+ T cell (Cheng et al., 2013).

#### Quantification of Spatial Structure Using Pair Correlation Function

Since we were interested in exploring how the spatial heterogeneity of glioblastoma tumors affect combination chemotherapy and ICB therapies, we sought to define a metric to compare and quantify the spatial configuration of the patient samples measured by IMC. For this, we leveraged the pair correlation function (PCF) () expressed as a function of separation distance between cells (r) as in previous work in ecology (Agnew et al., 2014; Surendran et al., 2019, 2020a,b). A PCF is defined as the ratio of density of pairs of cells with the square of the density of cells such that, in a population in which cells are randomly distributed (complete absence of spatial structure), . Here, the subscripts i and j represent the cell types of the two respective cells that form the pairs. To calculate the PCF for any population of cells, a reference cell at a location is randomly chosen and the distances from all other cells to the reference cell are computed. This procedure is repeated for all the cells in the sample (i.e., all cells act once as the reference cell). Once all possible pair distances are calculated, the PCF is constructed by enumerating the distances between pairs of agents that fall into the interval ensuring a bin width of. Normalizing the bin count by a factor of ensures that in the absence of spatial structure, where and are the number of cells of type i and j, respectively, and L is the length of the computational domain.

In a clustered spatial configuration, where pairs of cells are more likely to be found in close proximity, the PCF will have values for short distances of *r*. Conversely, when cells are less likely to be found close together, resulting in a segregated spatial pattern, for short distances of *r.* Thus, we measured the PCFs in each of the IMC images to quantify the extent of spatial structure in live cell populations which adopt uniform, clustered, or segregated spatial structure to varying degrees. Note that we only considered clustering of cancer cells, given that CD8+ T cells were sparse in our patient samples.

However, plotting and comparing PCFs as a function of distance for a large cohort of patients can be computationally challenging. To compare various possible combinations of spatial structure among multiple tumors, we also defined a simplified measure of spatial structure that expresses the type and extent of spatial structure as a scalar quantity. Hence, we defined a summary statistic based on the area under the PCF curve, i.e.,
as a convenient way to simply express the nature and extent of the spatial structure as a single number (Surendran et al., 2019). Thus, if cells are more (less) likely to be in proximity to others, then . If there is no correlation between cell locations, then . Supplemental Fig. 3 shows the distribution of (i.e., the metric measuring the correlation between glioblastoma cells) for the patient tumor samples in the glioblastoma cohort. We found that is insensitive to maximum distance between the cells, *H,* when we choose *H* to be sufficiently large since the PCFs approach unity for large separation distances. In our simulations, we set .

## Results

#### Increasing CD8+ T Cell Counts Can Improve Immune Checkpoint Blockade Efficacy

We first investigated how mono- and combination therapy impact tumor growth. For this, we compared the increase in the number of glioblastoma cells under various treatment conditions: 1) no treatment (baseline), 2) TMZ monotherapy, 3) ICB monotherapy, and 4) combination treatment with TMZ and ICB. To simulate each of the four cases, the respective therapies not being administered were set to zero. In Jenner et al., two model initializations were established based on random cell distributions within the tumor (according to cell counts established in ex vivo tumor slices) or “patchy” distributions (for more details, see the methods in Jenner et al., 2022). Before integrating patient IMC data to understand individual heterogeneity, we first studied the effects of treatment using these model initializations as a reference case to compare effectiveness of each of the different treatments. A snapshot of the initial configuration of all constituent glioblastoma, stromal, macrophage, CD4+, and CD8+ T cells for the random distribution of cells within the tumor is shown in Fig. 3. To create this, we placed cells randomly on the computational domain by sampling an angle () and a radius (), where R is the radius of the tumor.

Beginning from this random initialization, we sought to understand how each of the four treatment scenarios described above impact on overall tumor dynamics (Fig. 4, A–D). While our model predicted an overall reduction in cancer cell growth under TMZ monotherapy (Fig. 4B) compared with the case with no treatment (Fig. 4A), this did not result in total removal of glioblastoma cells, and there is still an overall positive growth in the presence of TMZ. This finding is consistent with clinical observations that patients receiving TMZ do not experience significant tumor reductions, particularly over short treatment periods. The ICB therapy was predicted to have a negligible effect on overall tumor growth (Fig. 4C), most likely due to the low number of adaptive immune cells initialized in these simulations. While the notion that the central nervous system and brain are completely immune-privileged is increasingly being challenged, the overall presence of CD8+ T cells in the brain is low compared with other tissues (Urban et al., 2020). Our results suggest that the failure of ICB monotherapy in glioblastoma can be in part attributed to the low number of CD8+ T cells present in the tumor (Fig. 3). Lastly, in the case of combination therapy, although we found reduced glioblastoma growth, it was largely driven by the effects of TMZ and not the presence of ICB (Fig. 4D), again due to the lack of cytotoxic T cells in the tumor field.

To better understand how CD8+ T cell recruitment to the tumor site affects treatment outcomes, we next varied the initial number of CD8+ T cells in the tumor domain by 5-, 10-, and 20-fold the initial amount (22 CD8+ T cells, see Jenner et al., 2022). We again simulated the administration of TMZ and ICB for each of these cases and looked at the change in glioblastoma cell counts in each of these scenarios. Note that initial parameters in ABM simulations, such as the number of glioblastoma cells and the tumor radius, can influence the individual tumor growth trajectories due to the pressure-dependent proliferation. However, our focus here is on looking at how varying initial CD8+ T cell count for a specific tumor sample impacts its growth dynamics. Our results showed a reduction in glioblastoma cell growth as the number of CD8+ T cells in the tumor increased (Fig. 5). This result suggests that patients with higher numbers of cytotoxic T cells within their tumor would benefit more greatly from treatment with ICB, and further, that this biomarker may act as a differentiator between response and non-response.

#### Spatial Configuration of the Tumor Impacts on Therapeutic Success

Glioblastomas are significantly spatially heterogeneous, as shown in the histopathological analyses in Jenner et al. (Jenner et al., 2022) which revealed critical differences in the size and distribution of cancer cell-enriched regions within their ex vivo tumor spheroids. Jenner et al. showed that this spatial heterogeneity determines the efficacy of oncolytic virotherapy beyond the impact of CD8+ T cell numbers. Therefore, to similarly elucidate how glioblastoma cell clustering and spatial heterogeneity impact on the therapeutic success of combination TMZ and ICB therapy, we next used our agent-based glioblastoma model to simulate two tumors with different (uniform and clustered) spatial configurations under combination TMZ and nivolumab. To generate the clustered configuration, we used the method described in Jenner et al. (Jenner et al., 2022) where they considered 50% of the tumor was comprised of cancer cell-enriched regions. We then randomly assigned five points as the center of the clusters and sampled the size (radius) of each of the clusters from a distribution informed by the approximate radius of clusters in the patient samples described in Jenner et al., 2022. Snapshots of the initial tumor samples with the uniform and clustered configuration are shown in Fig. 6, A and B, respectively.

For the uniform and clustered tumor samples, we simulated the agent-based glioblastoma model under treatment by the combination regimen and quantified the change in glioblastoma cell counts. We found that the clustered tumor results in worse treatment outcomes (Fig. 6D) compared with a more homogeneous cell distribution (Fig. 6C). As seen in Fig. 6E, the uniformly randomized tumor showed drastically decreased tumor growth over the short treatment period, suggesting that glioblastoma cells within clusters may get shielded from the drug action and thereby evade therapy within heterogenous tumors. This result may also explain the failure of mono-ICB in glioblastomas. Additional simulation results of tumor growths in absence of treatment in both clustered and uniform tumors (Supplemental Fig. 4) further verify this observation since both types of tumors have similar growth dynamics when no treatment is applied.

#### Interpatient Spatial Heterogeneity and Treatment Outcome

Given our previous results, we returned to our IMC data to predict how interpatient variability in tumor spatial configuration influences treatment outcomes. We used the pair correlation function expressed as a scalar metric of spatial structure (, see eq. 21) to rank patient samples according to the extent of spatial clustering. We then compared and quantified the extent of spatial structure and observed significant variability in the spatial organization these samples (Fig. 7).

We then initialized our agent-based model according to four representative patient samples (Fig. 7, A–D) and simulated the administration of combination chemo-immunotherapy, focusing on the change in glioblastoma cell counts. Note that the number of cells of each type is not identical among these samples since they correspond to distinct patients. Hence for comparison, we computed fold increases to glioblastoma cell numbers to accommodate differences in initial glioblastoma cell counts across samples. Consistent with the findings from our simulations based on the synthetic patient tumor sample (see Section: *Spatial configuration of the tumor impacts on therapeutic success*), we found that treatment outcomes worsened with increased clustering (Fig. 8). These findings suggest that the variability in the spatial configuration of the tumors is also a significant factor in determining treatment success, and as in oncolytic virotherapy, efforts should be made to account for spatial heterogeneity for clinical treatment decision making.

#### Predicted Attributes of Better Responses Remain Consistent in Brain Metastases

In the results described thus far, we identified that the extent of T cell recruitment and infiltration, and the spatial architecture of the tumor have a significant impact on responses to combined TMZ and ICB. To validate our findings, we used data from IMC analyses (see Methods) of brain metastases from primary lung cancers, breast cancers, and melanomas (Karimi et al., 2023) to establish whether their architectures and degrees of T cell infiltration differed from our glioblastoma samples, and whether these potential differences lead to better predicted responses. Brain metastases typically respond better than primary brain tumors to ICBs (Aquilanti and Brastianos, 2020; Lim et al., 2022). Hence, we investigated whether brain metastases had attributes similar to the predicted features of better responders, such as the presence of more T cells (Fig. 5).

Leveraging the IMC spatial snapshots from our glioblastoma and brain metastases samples, we first computed CD8+ T cell counts and the cross-correlation between cancer cells and CD8+ T cells expressed as a scalar metric, i.e., (Fig. 9, A and B). Note that can be computed by following the method described in section quantification of spatial structure using pair correlation function, and by specifically considering the distances between cancer cells and CD8+ T cells. Computing the correlation between CD8+ T cells () may provide insight into the possible occurrence of CD8+ T cell clusters and their potential impact on treatment efficacy. However, the number of CD8+ T cells in all the samples we considered are relatively low, and under these conditions, accurately computing , is not possible. A two sample T-test performed using *ttest2* function in MATLAB with α-significance level of 0.05 revealed a significant difference in the mean number of CD8+ T cells between these two cohorts, with brain metastases exhibiting overall higher CD8+ T cell counts (*P* value of 0.000475). No statistically significant difference in the mean value of was found between these two cohorts (*P* value of 0.8403).

To explore potential differences in treatment responses between glioblastoma and brain metastases, we again chose a representative subset of patients and simulated the administration of combination TMZ and ICB therapy. We selected three brain metastases (Fig. 9, C–E) and two glioblastoma patients (Fig. 9, F–G) based on the density of CD8+ T cells and the relationship between CD8+ T cell and cancer cell localization calculated as the product of a sample’s CD8+ T cell count and (with calculated from the cross-correlation function, , as described above—see Fig. 9, H–L). Put otherwise, this metric distinguishes patients with higher CD8+ T cell counts positioned closer to cancer cells. In all the cases considered here, we choose samples that maximized this metric.

Our results showed the fold increase in cancer cell counts after combination treatment to be consistent with our previous predictions, namely that tumors with higher recruitment/infiltration of CD8+ T cells respond better to combination chemotherapy and ICB, and clinical observations that brain metastases have better responses to immune checkpoint inhibitors as compared with primary glioblastomas (Fig. 9M), thus validating our model’s predictions. These results further underline the role of tumor architecture and CD8+ T cell infiltration as a biomarker of response and non-response.

## Discussion

Glioblastoma continues to be one of the most lethal brain tumors, despite intensive standard-of-care treatment. New treatment modalities, such as ICB, have been extensively tested both experimentally and in clinical trials with the goal of improving treatment outcomes for this difficult-to-treat cancer. Regrettably, ICB has thus far failed to show efficacy with respect to survival in clinical trials. Consequently, a better understanding of the tumor microenvironment and spatial effects determining responses to ICB, particularly in combination with approved treatment modalities like temozolomide, is paramount. In response, we developed a computational agent-based model that incorporates chemotherapy and PD1/PD-L1 immune checkpoint inhibition to treat glioblastoma to predict the key factors affecting the effectiveness of combination therapeutic approaches. Notably, our model describes key immunologic (e.g., macrophages, CD4+ and CD8+ T cells), phenotypic (e.g., stromal cells), and spatial interactions to provide a rational modeling framework for preclinical study.

Consistent with our previous study using an agent-based model of glioblastoma (Jenner et al., 2022), we first explored the role of CD8+ T cell numbers and localization on short-term therapeutic outcomes after monotherapies (TMZ or ICB alone) and in combination. Our model predicted that the efficacy of the immune response is a function of the amount of recruited CD8+ T cells. Since glioblastoma tumors are highly immunosuppressive and the ICB works not by directly targeting the cancer but rather enabling the patient’s CD8+ T cells to attack and remove cancerous cells, a sufficient number of CD8+ T cells must be present in the tumor for ICB to be effective. Indeed, these results reiterate the fact that both CD4+ and CD8+ T cell recruitment and specificity play a major role in ICB therapy. A potential avenue for further research here is to test combinations of ICB with immunostimulatory agents or oncolytic viruses that kill tumor cells and cause the release of danger signals to elicit an effective and sufficiently strong antitumoral immune response.

Facilitated by imaging mass cytometry visualizations of patient tumor resections, we quantified the spatial configuration of the tumor microenvironment to assess the role of spatial heterogeneity on treatment outcomes. To better understand the effect of the tumor microenvironment, we considered random (i.e., largely spatially homogeneous) and “patchy” (i.e., spatially heterogeneous) tumors and used our ABM to predict how spatial heterogeneity determines therapeutic efficacy. Varying microenvironmental conditions, such as immune infiltration, nutrients, and drug diffusion, exert selective pressure and shape cancer development. Treatment resistance due to selection is a major barrier to effective cancer treatments. Considering our results from a cancer evolution perspective, spatial heterogeneity can amplify these effects. For example, we found that the glioblastoma cells within clusters may get shielded from the drug action and thus evade therapy within heterogenous tumors. Thus, together, our results suggest that the varied spatial structure of glioblastomas contributes greatly to the challenges of treating these tumors. This is due to therapies not being equally effective throughout heterogeneous tumors (Figs. 7 and 8), highlighting the need for customizing the treatments to target specific tumor microenvironments. However, taking account of tumor heterogeneity when designing treatment is practically difficult to achieve and should be an area of focus going forward, particularly in preclinical models for which we can obtain better quantification. Crucially, we validated our model’s predictions using IMC data from brain metastases, which have shown greater responsiveness to immune checkpoint inhibition (Aquilanti and Brastianos, 2020). Comparing representative samples from patients with glioblastoma to patients with brain metastases, we found that the number of recruited and infiltrating CD8+ T cells was determinant of the response to treatment, corroborating our model’s predictions of why ICB has not shown efficacy for the treatment of glioblastoma and underlining the role of CD8+ T cells for immune checkpoint inhibition success.

There are limitations to our approach. Most notably, we considered only the short-term (i.e., up to 7 days) effects of therapy on glioblastoma due to the computational complexity of long-term simulations of our model (which models several tens of thousands of cells/agents). In future work, we will integrate both homogenous (ordinary differential equations) models of combination glioblastoma treatment with TMZ and ICBs and our ABM to provide more comprehensive and long-term predictions of therapeutic scheduling while still accounting for the spatial effects of treatment. Further, our model simplifies the immunologic response within the tumor (e.g., a thresholding effects for the ICB efficacy, simplified accounting of immune cell subsets acting within the tumor). Future work should explore the impact of increasing the complexity of immunologic interactions on model predictions.

Quantitative systems pharmacology and computational modeling are increasingly recognized for their contributions to preclinical drug development and for improving our understanding of the mechanisms underlying drug responses. Given the severity of glioblastoma, better insight into the drivers of response to new treatment modalities, like ICB, have the potential to provide significant benefits to patients. The modeling framework developed here is general enough to be adapted and extended to include other treatments, such as oncolytic viral therapy, in addition to the ones considered in this study. Thus, our study provides a rationale for mechanistic modeling in the preclinical space, and motivates the continued integration of experimental, clinical, and predictive models.

## Acknowledgments

The authors wish to thank the patients who participated in this study.

## Data Availability

The computational code to implement our agent-based model is available on GitHub at https://github.com/Anudeep-Surendran/GBM_TMZ_ICB.

## Authorship Contributions

*Participated in research design:* Surendran, Jenner, Quail, Walsh, Craig.

*Conducted experiments:* Karimi, Quail, Walsh.

*Contributed new reagents or analytic tools:* Surendran, Jenner, Craig.

*Performed data analysis:* Surendran, Jenner, Fiset, Craig.

*Wrote or contributed to the writing of the manuscript:* Surendran, Jenner, Karimi, Fiset, Quail, Walsh, Craig.

## Footnotes

- Received December 15, 2022.
- Accepted June 13, 2023.
AS was funded by a Centre de recherches mathématiques and Institut de sciences mathématiques postdoctoral fellowship and NSERC Discovery Grant [RGPIN-2018-04546] (to MC). ALJ was funded by a Fonds de recherche du Québec-Santé Programme de bourse de formation postdoctorale pour les citoyens d’autres pays and a Centre for Applied Mathematics in Bioscience and Medicine (CAMBAM) postdoctoral fellowship. IMC was supported by the Brain Tumor Funders’ Collaborative (DFQ and LAW). MC is an FRQS J1 Research Scholar and was supported by the Fondation du CHU Sainte-Justine.

No author has an actual or perceived conflict of interest with the contents of this article.

↵This article has supplemental material available at jpet.aspetjournals.org.

## Abbreviations

- ABM
- agent-based model
- ICB
- immune checkpoint blockade
- IMC
- imaging mass cytometry
- PCF
- pair correlation function
- PD-1
- programmed death protein-1
- PD-L1
- programmed death ligand-1
- TMZ
- temozolomide

- Copyright © 2023 by The Author(s)

This is an open access article distributed under the CC BY-NC Attribution 4.0 International license.