Agent-based computational modeling of glioblastoma predicts that stromal density is central to oncolytic virus efficacy

Abstract

Oncolytic viruses (OVs) are emerging cancer immunotherapy. Despite notable successes in the treatment of some tumors, OV therapy for central nervous system cancers has failed to show efficacy. We used an ex vivo tumor model developed from human glioblastoma tissue to evaluate the infiltration of herpes simplex OV rQNestin (oHSV-1) into glioblastoma tumors. We next leveraged our data to develop a computational, model of glioblastoma dynamics that accounts for cellular interactions within the tumor. Using our computational model, we found that low stromal density was highly predictive of oHSV-1 therapeutic success, suggesting that the efficacy of oHSV-1 in glioblastoma may be determined by stromal-to-tumor cell regional density. We validated these findings in heterogenous patient samples from brain metastatic adenocarcinoma. Our integrated modeling strategy can be applied to suggest mechanisms of therapeutic responses for central nervous system cancers and to facilitate the successful translation of OVs into the clinic.

Keywords: Cancer; Computational bioinformatics; Immunology.

Graphical abstract

Figure 1

Ex vivo tumor model Patient biopsies or resections were manually cut and cultured with peripheral blood mononuclear cells (PBMCs) and grown in ex vivo tumor ecosystems (Majumder et al., 2015) in the presence of oHSV-1. Tumor tissue slices were then formalin-fixed and paraffin-embedded (FFPE), and serial sections were collected and stained for HSV-1, cleaved caspase-3, or with hemotoxylin and eosin (H&E). Multiplex immunohistochemistry (mIHC) was performed to investigate the impact of oHSV-1 on the spatial context and activity of the immune compartment in the glioblastoma TME ex vivo.

Figure 2

Agent-based model of glioblastoma patient tumor tissue slices infected with oHSV-1 (A) Schematic describing the model agents: glioblastoma (GBM) cell, stromal cell, CD4⁺ T cell, and CD8⁺ T cell. Corresponding agent rules and cell–cell interactions are in Table 1. Glioblastoma (purple) and stromal (pink) cells were randomly distributed in a 2D domain inside a circle of radius $R$ with initial numbers depending on the stromal cell to glioblastoma cell ratio. oHSV-1 particles were modeled to diffuse through the domain (Equation 1, STAR Methods), bind, and become internalized by glioblastoma or stromal cells. In infected glioblastoma cells (brown cells), the virus then undergoes replication and eventually lyses the infected cell releasing new viral particles and killing the cell (black cell) (Figure S4). CD4⁺ T cells and CD8⁺ T cells are both present in the tissue and contribute to localized clearance of infected cells. Stromal cells act as virus sinks. Through contact with an infected tumor cell, a naïve CD4⁺ T cells (Th) will become activated and secrete chemokines which attract CD8⁺ T cells (CTL) through chemotaxis. CD8+T cells kill infected cells they encounter (Figure S5). (B) A patient-derived resected glioblastoma was sub-sectioned into two smaller slices, which were then formalin-fixed, paraffin-embedded, sectioned onto slides, and stained with hematoxylin and eosin (H&E). A clinical pathologist scored the fraction of tumor cell, necrosis, immune cell, and stroma content (Table 2). Scale bar $=750\text{μm};$ further zoomed in inserts can be found in Figure S6. The two subsections were recreated in our agent-based model and designated as sparse (top) or dense (bottom) based on the pathology scores. Larger versions of the simulations can be found in Figure S7.

Figure 3

Validation of the computational model to ex vivo, in vivo, and in vitro measurements (A) The number of CD4⁺ and CD8+T cells in tumor slices (4 patients). Mean and standard deviation are given for the number of CD4+T cells (Th) cells and CD8⁺ T cells (CTLs) in each slice for patients 1098 (20 slices), 1167 (20 slices), 18IIOC-A (23 slices), and 18IIOC-B (28 slices). (B) Nearest neighbor distances between CD4⁺ (TH) and CD8⁺ (CT) T cells averaged across patients 1098, 1167, 18IIOC-A, and 18IIOC-B. Ki67+ proliferative (∗) and Ki67-non-proliferative subsets of cells were noted and for each cell–cell pair relationship, the distance from a randomly chosen cell of the first type to its nearest neighbor of the second type was calculated. Error bars: standard deviation (individual patient data in Figure S13). The average minimum distance between each cell type TH, TH∗, CT, and CT∗ after a single simulation of the Hooke’s law simulation (Figure S14) is overlaid in green. (C) To model glioblastoma cell proliferation, a logistic model was fit to in vitro cell counts for U87 cell (Mercurio et al., 2017) growth, the blue error bars are the mean and standard deviation of the cell count measurements and the black solid line is the logistic growth approximation (Figure S16). (D) Average viral infiltration in 50 μm wide bands from the periphery of ex vivo tumor slices (n = 5) measured at 24 h with an initial 5 PFU per cell by the detection of GFP label (other PFU concentrations in Figure S15). Overlaid is the fit for viral diffusion. (E) $S_{chemokine}$ and ${\rho }_{chemokine}^{\ast }$ were estimated from data from Gao et al. (Gao et al., 2014), who measured the density of IFN-$\gamma$ secreted from cells. Blue solid line and circles: data. Purple solid line: model fit. (F) Example of high-grade glioma stained for HSV-1. Scale bars represent $200\text{μm}$ and $100\text{μm}$ in the left and right panels, respectively.

Figure 4

Computational model predicts the ex vivo infiltration of oHSV-1 in glioblastoma tissue Simulations of the computational model for fully dense sample 100:0 (first column), fully sparse sample 0:100 (second column), and patchy 50:50 dense to sparse (third column) tumor tissue configurations. The resulting viral infiltration in tumor tissue configurations at 6 days is given for three replicates of each tumor type. Uninfected glioblastoma cells (purple), infected glioblastoma tumor cells (brown), stromal cells (pink). See Figure S17 for corresponding spatial densities of chemokine and virus and total local density. The resulting % of live uninfected and infected cells relative to the total number of cells over time is plotted along with the total virion count (n = 9).

Figure 5

Systemic evaluation of parameter contributions to model behavior Parameter sensitivity was evaluated by running simulations for perturbations in the parameters $D_{virus},{\lambda }_{virus},\beta ,$ and $\alpha$. Each parameter was perturbed by $\left[0.1,0.5,1,1.5,2,10\right]\times p$, where $p$ is the parameters original value (see Table S2). The average remaining glioblastoma cells on day 3 from 10 simulations are reported in the heatmaps (left). The corresponding model spatial tumor configuration for an instance of a parameter value is given in the heatmaps for (A) fully dense sample 100:0, (B) fully sparse sample 0:100, and (C) patchy 50:50 (dense:space).

Figure 6

Effect of T cell recruitment, antigen specificity, and virus binding rates of stromal and glioblastoma cells on oHSV-1 simulations (A) The number of CD8⁺ and CD4⁺ T cells in the simulation was increased from the control and base case 2-times, 10-times, and 100-times. The control simulation considered no OV administration, whereas all other simulations included the OV. The percentage of tissue remaining at 72 h (top) and CD8⁺ T cell (CTL) percentage increase over time (bottom) were quantified for the different immune cell numbers considered (infected cells and CD4⁺ T cells (THs) over time are in Figure S22). In addition, we considered CD8⁺ T cells as either virus-specific (induce apoptosis in infected cells only) or glioblastoma and virus specific (induce apoptosis in both uninfected and infected cells). (B) The percentage of the initial tissue remaining at 72 h for different glioblastoma and stromal cell relative binding rates $\text{Ω}$, where ${\text{u}}_{\text{s}}=0.01\times \text{Ω}$ and ${\text{u}}_{\text{g}}=0.002\times \text{Ω}$ in dense and sparse tissues is given as mean and standard deviation error bars (n = 3). p-value tests for significance ($p>0.05$) with not significant (ns) is noted. (C) Simulation snapshots of (bottom) dense and (top) sparse tumor tissue at 144 h for the stromal cell relative binding rates $\text{Ω}$ of $\text{Ω}=0.5$ and $\text{Ω}=10$. (D) The percentage of the initial tissue remaining at 72 h for different stromal binding rates (${\text{u}}_{\text{s}}$) in both dense and sparse tissues. Corresponding viral infiltration measurements are in Figure S19. (E) Simulation snapshots of (bottom) dense and (top) sparse tumor tissue at 144 h for stromal cell-virus binding rate $u_s=5\times {10}^{-4}$ and $u_s=1$. Stromal cells (pink), glioblastoma cells (purple) and infected cells (brown).

Figure 7

Tumor and stromal cell density determines oHSV-1 efficacy (A) A range of tumor tissue configurations was generated with dense and sparse regions making up the proportion of total area denoted dense:sparse (Figure S11). In these simulations, the stromal binding rate was set at $u_s=1$. For different proportions of sparse and dense regions, the fraction of tumor tissue remaining was calculated at 72 h. The color of the bar represents the initial number of glioblastoma cells in that tumor tissue. For the 20:80 and 70:30 tumor configuration, the model snapshot is given on day 3. The corresponding time-evolution for 0:100, 50:50, and 100:0 tumors is in Figure 4. (B) The percentage of glioblastoma tumor cells remaining over six days for three different administration protocols: a single injection into the center of the dense slice (green), a homogeneous administration of the OV on the periphery (orange), and six OV injections each into the center of dense patches for a 50:50 slice (blue). The total amount of virus administered was conserved across the protocols. (C and D) Brain metastatic adenocarcinoma (N = 1) was isolated for ex vivo experiments and multiple independent biological replicates were then treated with oHSV1 at PFU = 5 per cell (N = 3). Following the sequential imaging strategy (described in STAR Methods), regions were characterized by clinical pathology as retaining “high” or “low” tumor cell density. Optical density of HSV-1 and caspase-3 immunostain was then quantified by automated image analysis and graphed as (E) HSV-1 expression or (F) cleaved caspase-3 per region of tissue area (${\text{mm}}^2$).