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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Glioblastoma TMZ and immune checkpoint inhibitor agent-based model. Schematic summarizing the interactions between different cell types as well as drugs. The agent-based model describes and tracks glioblastoma tumor cells, stromal cells, macrophages, and CD4+ and CD8+ T cells (see inset legend). The chemotherapy drug temozolomide (yellow) reduces proliferation and enhances the death of glioblastoma cells (purple). Macrophages (green) become activated through the phagocytosis of apoptosed glioblastoma cells (black). CD4+ T cells (orange) become activated through direct interactions with a glioblastoma cell. Once activated, CD4+ T cells secrete chemokines that attract CD8+ T cells (blue) through chemotaxis. CD8+ T cells are activated either by CD4+T cells or macrophages. CD8+ T cell-induces apoptosis of glioblastoma cells in the presence of the immune checkpoint inhibitor drug nivolumab which binds to PD-1 receptors and reactivates the immune response.

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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Summary of the three-compartment pharmacokinetic TMZ model developed by Ostermann et al. A) The amount of TMZ in the GIT, blood, and brain (CSF) is modeled with eqs. 2–6. Details on the population pharmacokinetic model can be found in Ostermann et al., 2004. B) Model simulations predicting the amount of TMZ in the GIT and plasma (top), and the concentration in the CSF (bottom).

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, ² 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)², 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 .