Agent-based model of glioblastomaWe recently developed an ABM for glioblastoma growth and treatment with either an OV22, or an ICI and TMZ32 in PhysiCell, an open-source cell-based simulator50. In this framework, cells are considered agents and move off-lattice in a 2D domain with movement governed by a series of force-balance equations. Substrate diffusion, secretion and uptake are modelled as a system of partial differential equations (PDEs) using BioFVM51. BioFVM is a biotransport solver that simulates diffusing substrates and cell-secreted signals in the microenvironment. With BioFVM, cells also update their phenotype depending on their local microenvironment conditions. PhysiCell has extensive documentation and numerous extensions and examples of applications52,53,54,55. Here, we briefly describe our main modelling assumptions from our previous glioblastoma growth, OV, TMZ and ICI models.Our model assumes we have five main cell types: glioblastoma cells, CD8 + T cells, CD4 + T cells, macrophages, and stromal cells (Fig. 5a). There are four substrates considered in the model: OV, chemokine, ICI, and TMZ, all of which are modelled by the BioFVM PDE:$$\begin{array}{l}\frac{\partial \rho }{\partial t}=D{\nabla }^{2}\rho -\lambda \rho +S\left({p}^{* }-p\right)-{Up}\\\quad\,\,+\mathop{\sum}\limits_{{cells}}\delta \left(x-{x}_{k}\right){W}_{k}\left[{S}_{k}\left({\rho }_{k}^{* }-\rho \right)-{U}_{k}\rho \right],\end{array}$$
(1)
where \(D\) is the diffusion coefficient, \(\lambda\) the decay rate, \(S\) the secretion rate with maximum density \({\rho }^{* }\), \(U\) is the uptake rate, \({W}_{k}\) the cell volume, \({x}_{k}\) the cell position. All estimates for these parameters can be found in our previous work22,32.Fig. 5: Model and treatment schematics.a ABM describing combination treatment with TMZ, an anti-PD-1 immune checkpoint inhibitor, and an HSV-1 oncolytic virus. b We simulated the first seven days of treatment wherein the patient first receives a single administration of the oncolytic virus with the anti-PD-1 agent administered over the first hour; TMZ is given orally on days one to five. To account for the oral administration of TMZ, we use a three-compartment pharmacokinetic model to describe drug concentrations in the gastrointestinal tract, blood, and cerebrospinal fluid. Details can be found in “Treatment with temozolomide”. The oncolytic virus dosage depends on the initial cell count and is equivalent to a dose of 5 plaque-forming units (PFU) per cell for 9740 cells. PFU is a measure of virus quantity used in virology. More details on how PFU relates to viral dose size are given in “Treatment with an oncolytic virus”.As mentioned above, cells are considered individual agents that each have their own set of rules for their interactions with other cells as well as their local microenvironment. At every cell update timestep, the simulator randomly picks a cell and goes through a decision tree to determine its outcome at the next timestep (Fig. 6). Once the cell is updated, another cell is randomly selected, and this process continues until all the cells have been updated. The decision tree used in this paper is the combination of the two decision trees from our previous models22,32.Fig. 6: Flow chart of the decision tree implemented for the update of cells in the agent-based model of glioblastoma.In each cell time step, a cell agent in the simulation is chosen at random and, based on its cell type, follows a series of decision tree rules. The time step is updated only once all cell agents have been evaluated.We briefly describe the core rules of the model. Glioblastoma cells proliferate at a baseline rate \(\beta\). The probability for a glioblastoma cell to proliferate depends on the local cell density: if an individual cell is too crowded, it will not proliferate. The local cell density is correlated to the mechanical pressure from neighbouring cells. As in ref. 22, we use PhysiCell’s inbuilt cell–cell pressure calculation to compute the mechanical pressure felt by an individual cell. We define a mechanical pressure threshold \({p}_{\max }\) and cells proliferate at a rate \(\beta\) if their mechanical pressure from neighbouring cells is less or equal to this threshold. For a given cell \(i\), PhysiCell defines its set of neighbours \(\Pi \left(i\right)\) as the list of cells located in neighbouring voxels, i.e., those that share an edge with \(i\)’s current voxel, and the cells in \(i\)’s current voxel22,50.Glioblastoma cells die in three main ways: (1) a cell becomes infected by an oncolytic viral particle and dies through lysis, (2) TMZ induces cell death (see “Treatment with temozolomide”), and (3) a CD8 + T cell induces apoptosis26,32. Macrophages are activated once they phagocytose a dead cell, after which they can activate a CD4+ or a CD8 + T cell in their neighbourhood and ingest a certain number of cells; they then become exhausted, and their death rate is increased. CD4 + T cells become activated through interactions with a tumour cell in their neighbourhood. Once activated, they secrete chemokines that recruit CD8 + T cells to the site of the tumour. In the presence of an anti-PD-1 treatment, if the CD8 + T cell has more than \(N\) ICI-bound PD-1 receptors on its surface, it can attach itself to tumour cells in its neighbourhood. Once a CD8 + T cell has been attached for more than \(\tau\) minutes, it kills the tumour cell. We refer readers to ref. 32 and ref. 22 for complete details.Initialisation of the glioblastoma model with IMC dataTo initialise the tumour in our ABM, we incorporated IMC data of tumour resections from untreated glioblastoma patients from ref. 45. All data are published and publicly available; a complete description with details of ethics oversight is available in ref. 45. As described in ref. 45, written informed patient consent was obtained prior to collecting surgical specimens and clinical information. IMC on glioblastoma patient samples was performed as previously described by Karimi et al.45 and Surendran et al.32. We used spatial data from ref. 45 to localise cancer cells, macrophages, CD4 + T cells, CD8+ cells, and endothelial cells that we considered to be stromal cells56, and implemented it into our agent-based model to individualise the simulations (Fig. 7). Since microglia cells and macrophages have similar functions in the glioblastoma TME and are often recognised as one cell cluster57,58, we included microglia within the macrophage category.Fig. 7: From Imaging mass cytometry data to an individualised agent-based model.a IMC images from Karimi et al. (top) and corresponding lineage assignment (bottom), with magnified regions to the right of each image. Image reproduced under Creative Commons. b Initial spatial configurations of six chosen IMC data samples. The six IMC samples were chosen according to their adaptive immune cell proportions. Thus, we first ranked patients in our cohort by their initial CD4+ and CD8 + T cells proportions. In the case where two patients had the same initial adaptive immune cell proportions, we ranked patients according to the respective initial macrophage proportions. We then selected the bottom three patients (B1–B2–B3) and the three top patients (T1–T2–T3) as representative individuals. c Cell distributions of 145 patients in the percentage of all cells in the TME, measured in IMC. Patients are ordered according to their adaptive immune cell proportions (CD4+ and CD8 + T cells), in ascending order from left to right. Distributions corresponding to the six chosen patient samples are indicated with arrows.IMC images were all 1 mm2. Supplementary Table 1 summarises patient-related data. Samples contained between 910 and 6033 total cells and displayed a high degree of heterogeneity in the distributions of cell populations as a percentage of all cells in the TME (Fig. 7). The proportion of adaptive immune cells (CD4 + T cells and CD8 + T cells) varied between 0% and 7.38%, with an average of 0.66%. On the other hand, samples had, on average, a cancer cell proportion of 50.25%, ([0.84%; 89.54%]). Stromal cell proportions were between 0.2% and 10.54%, with an average of 2.87%. Finally, macrophages made up between 8.79% and 94.38% of samples, with an average of 46.23%.Treatment with temozolomideWe used the same TMZ pharmacokinetic/pharmacodynamic model as in ref. 32. To predict TMZ concentrations in the cerebospinal fluid (CSF), we used a three-compartment model with first-order absorption in the gastrointestinal tract (GIT) from ref. 59. TMZ pharmacokinetics were modelled according to$$\frac{d{A}_{1}}{{dt}}=-{k}_{a}{A}_{1},$$
(2)
$$\frac{d{A}_{2}}{{dt}}={k}_{a}{A}_{1}-\frac{{CL}}{{V}_{D}}{A}_{2}-{k}_{23}{A}_{2}+{k}_{32}{A}_{3},$$
(3)
$$\frac{d{A}_{3}}{{dt}}={k}_{23}{A}_{2}-{k}_{32}{A}_{3},$$
(4)
$${C}_{{CSF}}=\frac{{A}_{3}}{{V}_{P}}.$$
(5)
where \({A}_{1}\), \({A}_{2}\) and \({A}_{2}\) are the amount of TMZ in the GIT, plasma, and CSF compartments, respectively (with \({C}_{{plasma}}\) and \({C}_{{CSF}}\) the TMZ concentrations in the plasma and the CSF compartments), \({k}_{a}\) is the absorption rate, \({k}_{23}\) and \({k}_{32}\) are the rates of exchange from the plasma to CSF and from the CSF to plasma, respectively, \({CL}\) is the TMZ clearance, \({V}_{D}\) is its volume of distribution in the central compartment, and \({V}_{P}\) is its volume of distribution in the CSF compartment. The integration of this three-compartment model into the ABM and the variations in local TMZ concentrations tracked by the BioFVM biotransport system are detailed in ref. 32.TMZ is considered to be primarily cytostatic60 and to induce cell death in proliferating cells by arresting their cell cycle61. As such, we assumed that TMZ results in reduced proliferation and in increased death of glioblastoma cells. As in ref. 32, the cytostatic effect of TMZ on glioblastoma cells can be described by the following three equations:$$\frac{{dGBM}}{{dt}}=\beta \left(1-\frac{{TMZ}}{{TMZ}+I{C}_{50}}\right){GBM},$$
(6)
$$\frac{{dD}}{{dt}}=\beta \frac{{TMZ}}{{TMZ}+I{C}_{50}}{GBM},$$
(7)
$$\frac{{dTMZ}}{{dt}}=-{{\rm{\lambda }}}_{{TMZ}}{TMZ}.$$
(8)
Here, \({GBM}\) is the number of glioblastoma cells, \({TMZ}\) is the temozolomide concentration in the cerebrospinal fluid, and \(D\) is the number of dead cells after cell cycle arrest induced by TMZ. Finally, \(I{C}_{50}\) is the half effect of TMZ. Surendran et al.32. adapted this system of equations to model stochastic cell division and death with PhysiCell’s inbuilt “Live” and “Apoptosis” models50. According to the “Live” cell division model, a cell can divide with a division rate \(\varphi\) such that in a given time interval, \(\left[t,t+\Delta t\right]\), the probability of cell division is given by Ghaffarizadeh et al.50$${\mathbb{P}}\left({division}\right)=1-{e}^{-\varphi \Delta t}\,\approx \,\varphi \Delta t$$
(9)
As in ref. 32, we defined \(\varphi =\beta \left(1-\frac{{TMZ}}{{TMZ}+I{C}_{50}}\right)\). Similarly, the ”Apoptosis” model defines the probability of a cell dying in the time interval \(\left[t,t+\Delta t\right]\) to be50$${\mathbb{P}}\left({cell\; death}\right)=\psi \Delta t,$$
(10)
with \(\psi\) representing the death rate. Based on ref. 32 in our model, the death rate of glioblastoma cells is given by \(\beta \frac{{TMZ}}{{TMZ}+I{C}_{50}}\). Surendran et al.32 fitted the standard Emax effects curve to dose–response data from ref. 62 and found an estimate of \(I{C}_{50}=0.0162\frac{\mu g}{\mu L}\).Since TMZ reduces the probability of cell division of proliferating glioblastoma cells, it has no effect on cells whose proliferation is already blocked for another reason, such as infection by an oncolytic virus. Given that we lack specific data of combination HSV-1 oncolytic virus with TMZ, we consider this modelling choice to effectively capture the OV treatment in combination with TMZ while preventing our model from overpredicting the combined effect of therapy.As we did in our previous model, we considered a dosing schedule where 175 mg/m2/day of TMZ was administered for 5 consecutive days, which is consistent with standard-of-care for glioblastoma2,32,59. Supposing an average adult male has a body surface area of 1.9 m2, this dosing schedule corresponds to 332.5 mg of TMZ per day for 5 days.Treatment with immune checkpoint inhibitorsWe kept the same glioblastoma growth model with a treatment with nivolumab, an anti-PD-1 drug, as in our initial TMZ + ICI model32, which was based upon a previous model published by Storey et al.63. For each time \({t}_{d}\) (in minutes) at which nivolumab is administered,$$A\left(t\right)=\left\{\begin{array}{c}0.481,\\ 0\end{array}\right.\begin{array}{c}{t}_{d}\, < \,t \,< \,{t}_{d}+60\\ {\rm{otherwise}}\end{array},$$
(11)
where \(A(t)\) is the source of anti-PD-1 drug concentration as a function of time. Variations in the local concentration of the anti-PD-1 drug (\({ICI}\)) are tracked by the BioFVM partial differential equation$$\frac{\partial {ICI}}{{dt}}={D}_{{ICI}}{\nabla }^{2}{ICI}-{\lambda }_{{ICI}}{ICI},$$
(12)
as detailed in ref. 32. \({D}_{{ICI}}\) and \({\lambda }_{{ICI}}\) are the ICI diffusion coefficient and decay rate, respectively. The standard administration schedule of nivolumab is a single intravenous dose of \(240\) mg every two weeks administered for one hour32,63. As we performed 7-day long simulations, we administered nivolumab for the first hour of the simulations.As in ref. 32, we assumed that CD8 + T cell induced apoptosis in glioblastoma patients can only occur when the number of ICI-bound PD-1 receptors on a T cell is above a threshold \({N}_{T}\). For each CD8 + T cell, its number of ICI-bound (\({R}_{B}\)) and unbound (\({R}_{U}\)) receptors is tracked by$$\frac{d{R}_{B}}{{dt}}={k}_{{on}}{N}_{{ICI}}{R}_{U}-{k}_{{off}}{R}_{B},$$
(13)
$$\frac{d{R}_{U}}{{dt}}=-{k}_{{on}}{N}_{{ICI}}{R}_{U}+{k}_{{off}}{R}_{B},$$
(14)
where \({k}_{{on}}\) and \({k}_{{off}}\) are the binding rates, and \({N}_{{ICI}}\) is the amount of nivolumab in the local microenvironment of the CD8 + T cell tracked by Eq. (13), as detailed in ref. 32.Treatment with an oncolytic virusThe OV was administered as a single dose in the centre of the domain, in a circular region with a radius of \(20\mu m\), at the beginning of the simulations. Given that some of the patients’ IMC samples had blank spaces containing none of the cells considered in our study, to make sure that we administered the OV into a region that contains cells, we found the position of the closest cell to the centre of the domain and set its \(x\)- and \(y\)- positions as \({cente}{r}_{x}\) and \({cente}{r}_{y}\), respectively. The OV was then administered in the region$${{\mathscr{V}}}_{{admin}}=\{(x,y)|{(x-{cente}{r}_{x})}^{2}+{(y-{cente}{r}_{y})}^{2}\, < \,{20}^{2}\},$$meaning that every voxel in \({{\mathscr{V}}}_{{admin}}\) has an initial viral concentration equal to \({V}_{0}\).Details of the OV integration in the ABM can be found in ref. 22. As in ref. 22, we assumed a one-to-one relationship between the number of virions and plaque-forming units (PFU). PFU quantify the number of viral particles needed to produce one plaque by lysis after infection. In ref. 22, simulations were performed with 9740 initial cells and viral dosages of 5 plaque-forming units (PFU) per cell. Since the patient datasets used in our work had an unequal number of initial cells, we calculated the viral dosages \({PF}{U}_{x}\) (i.e., the PFU specific to patient \(x\)) according to$${PF}{U}_{x}=\frac{5\frac{{PFU}}{{cell}}}{9740{cells}}\cdot {C}_{{0}_{x}}{cells},$$
(15)
where \({C}_{{0}_{x}}\) is the initial number of cells in patient \(x\)’s sample. The initial viral concentration administered for patient \(x\) is thus given by$${V}_{{0}_{x}}=\frac{{PF}{U}_{x}\cdot {C}_{{0}_{x}}}{\pi \cdot {R}^{2}}$$
(16)
where \(R\) is the radius of the regions of administration (\(R=20\mu m\)). At the beginning of the simulations, every voxel in \({{\mathscr{V}}}_{{admin}}\) as an initial viral concentration equal to \({V}_{{0}_{x}}\frac{{virions}}{\mu {m}^{2}}\), while every voxel outside of \({{\mathscr{V}}}_{{admin}}\) have an initial viral concentration equal to \(0\frac{{virions}}{{\rm{\mu }}{{\rm{m}}}^{2}}\).Patient sample stratification and simulationsOur previous TMZ + ICI model showed that patients with higher numbers of cytotoxic T cells within their tumour respond better to treatment with ICIs32. With this in mind, we ranked patients according to their initial adaptive immune cell proportion (i.e., CD4+ and CD8 + T cells) (Fig. 7c). Adaptive immune cell proportion ties were broken according to the macrophage proportions. We then picked the three patients with the lowest adaptive immune cell proportions (patients B1, B2 and B3), and the three patients with the highest adaptive immune cell proportions (patients T1, T2 and T3), and initialised an ABM for each of these patients according to their respective IMC data (Fig. 7b). Patient data for these samples are provided in Supplementary Table 2.For each of the six patients, we simulated the model without treatment, as well as five different combination treatments: TMZ, ICI, OV, OV + ICI and TMZ + ICI + OV (Fig. 5). Since an ABM is a stochastic model, two simulations with the same input will not produce the same output. Thus, for each treatment regimen of interest, we performed ten runs of our model and used the average when analysing outcomes. Based on our previous results32, we would expect patients T1, T2 and T3 to respond better to treatment combination schemes that include immunotherapies since these patients have a higher proportion of adaptive immune system compared to patients B1, B2 and B3.Simulations of representative patientsTo study the effects of increasing the domain size and the initial number of cells, we generated two virtual patients. To generate both representative patients, we found \({\mathscr{X}}\), the average proportion of the domain that contains cells in all the 145 patients’ samples and defined \(\epsilon\) as the average cell radius between the five considered cell types (macrophages, CD4 + T cells, CD8 + T cells, tumour cells, and stroma cells). We then calculated \(K\), the circular domain radius needed such that a proportion \({\mathscr{X}}\) of \(K\) contained cells when populated with 10,000 cells of radius equal to \(\epsilon\).We ranked the patients according to their initial adaptive immune cell proportion and generated a representative patient sample with 10,000 cells (BottomRepresentative) that respects the proportions \({P}_{{x}_{{bottom}}}\) representing the initial adaptive immune cell proportion of the median patient of the bottom 10% of the cohort with respect to the order described earlier (Fig. 7c). Similarly, we also generated a representative sample consisting of 10,000 cells that maintained the proportion \({P}_{{x}_{{top}}}\) (TopRepresentative—initial adaptive immune cell proportion of the median patient of the top 10% of the cohort). For both BottomRepresentative and TopRepresentative, cells were assigned an initial position in the domain according to a uniform distribution.As with the six patients’ samples, we simulated these representative samples ten times under the same treatment combinations (i.e., No treatment, TMZ, ICI, OV, OV + ICI, TMZ + ICI + OV). For both representative samples, \({V}_{0}\) and the region of OV administration were calculated as described in “Treatments with an oncolytic virus”. Schedules of administration for our 7-day long simulations are described in Fig. 5b.Confirmation of the model’s predictive capacityBoth the TMZ + ICI and the OV models were validated with data in our previous work by Surendran et al.32 and Jenner et al.22, respectively. Given that they generally respond more to checkpoint blockade than patients with primary glioblastoma19,32,64, Surendran et al.32 used samples from metastases from patients with primary lung cancers, breast cancers, and melanomas32 and found that brain metastases exhibited significantly higher CD8 + T-cell counts, which was one of their predicted features of better responders. Their model predictions also showed that brain metastases responded better to ICIs, which is consistent with clinical observations65.The OV model was calibrated to a number of in vitro, in vivo and ex vivo datasets, and predictions from the model were further validated against clinical measurements22. For example, the OV model suggests that increasing the diffusivity of an oncolytic virus and/or decreasing its clearance would significantly improve treatment efficacy. This was consistent with evidence that arming OVs with a transgene degrading the extracellular matrix, such as hyaluronidase66,67, decorin68,69, or relaxin, improves viral spread and enhances viral potency. Also, in glioblastoma, OV clearance was linked to natural killer cells that cause diminished OV-efficacy70. Inhibiting the recruitment of natural killer cells by combining TGB-\(\beta\) and the oHSV OV showed increased viral efficacy71 .To confirm our combination model’s predictive capacity, given the lack of data on the TMZ + OV + ICI combination, in this study we selected six patients’ samples based on their adaptive immune cells proportions to make predictions for two immune system extreme scenarios: patients with a weak adaptive immune system and patients with the strongest immune response. Similarly to the six chosen patients’ samples, the two virtual domains (Top and Bottom Representatives) illustrate these two immune system extreme scenarios (i.e., weak versus strong adaptive immune system). Doing so allowed us to confirm our model’s predictions aligned with biological expectations since we know that the increased presence of adaptive immune cells (i.e., CD4+ and CD8 + T cells) enhances ICI efficacy32.Statistical analysesFor each patient sample (IMC and representative), the fold increase in cancer cells for an administered treatment was defined as$${\rm{fold\; increase}}\left({tx}\right)=\frac{C\left(7\right)}{C\left(0\right)}$$
(17)
for \({tx}\in \left\{{No\; treatment},{TMZ},{ICI},{OV},{OV}+{ICI},{TMZ}+{ICI}+{OV}\right\}\), with \(C\left(t\right)\) the number of cancer cells after \(t\) days of simulation with treatment \({tx}\). Since we were interested in the association between model inputs and predictions of the different initialisations of the ABM, we measured correlations between the eight different model initialisations (real or representative patient) and various characteristics (input or output) of the simulations. The correlation coefficient between characteristics \(i\) and \(j\), \(\rho \left(i,j\right)\) was defined as the Pearson correlation coefficient, i.e.,$$\rho \left(i,j\right)=\frac{{Cov}\left(i,j\right)}{{\sigma }_{i}{\sigma }_{j}},$$where \({Cov}\left(i,j\right)\) is the covariance between the two characteristics, and \({\sigma }_{s}\) is the standard deviation of characteristic \(s\). Statistical significance was defined at \({\rm{\alpha }}=0.05\), and correlations were measured using corrcoef in Matlab2023b72. corrcoef computes P values with a t test by transforming the correlation to create a t statistic having \({\rm{n}}-2\) degrees of freedom, where \({\rm{n}}\) is the number of rows of the correlation table. Here, \({\rm{n}}=13\). It tests the null hypothesis that there is no relationship between the observed phenomena. The P value is the probability of getting a correlation as large as the observed value by random chance when the true correlation is zero.For the input characteristics, we looked at the total initial cell count, the initial proportion of each of the different cell types, the initial percentage of blank space in the model, and the initial tumour density. For the output characteristics, we defined treatment response as the average growth speed of the tumour cells$${gs}\left(x,{tx}\right)=\frac{C\left(i\right)-C(0)}{i},$$
(18)
where \({gs}\left(x,{tx}\right)\) is the tumour cells growth speed of the simulation \(x\) under the simulated treatment \({tx}\) and \(i\) is defined as$$i=\min \left(7,j\right)$$with$$j=\left\{\begin{array}{c}7\\ \min \left\{t:C\left(t\right)=0\right\}\end{array}\begin{array}{c}{\rm{if}}C\left(7\right)\, > \,0\\ {\rm{otherwise}}\end{array}.\right.$$If \(C\left(i\right) \,< \,C\left(0\right)\), then \({gs}\left(x,{tx}\right)\, < \,0\) and \({|gs}\left(x,{tx}\right)|\) represents the average killing speed of the tumour cells. In this case, the smaller the value of \({gs}\left(x,{tx}\right)\) is, the more patient \(x\) responds to \({tx}\). Conversely, if \(g\left(x,{tx}\right)\, > \,0\), then the cancer has grown over the simulated period.We defined the percentage of white space as the percentage of the initial tumour area that does not contain any of the cell types included in this study, i.e.,$${\%}_{{white\; space}}=\frac{{\sum}_{i=1}^{5}{{\#}}({cells\; of\; type\,i})\cdot {radius}{\left({cell\; of\; type\,i}\right)}^{2}\cdot \pi }{initial\; tumour\; area}.$$
(19)
To adequately assess differences across different conditions, we min–max normalised across the columns. Let \({TC}\) and \(T{C}^{* }\) be the \(8\times n\) not normalised and normalised tables of characteristics, respectively. Then,$$T{C}^{* }\left[i,j\right]=\frac{{TC}\left[i,j\right]-\min \left\{{TC}\left[s,j\right]:s=1,\cdots ,8\right\}}{\max \left\{{TC}\left[s,j\right]:s=1,\cdots ,8\right\}-\min \left\{{TC}\left[s,j\right]:s=1,\cdots ,8\right\}},$$
(20)
where \({TC}\left[i,j\right]\left(T{C}^{* }\left[i,j\right]\right)\) refers to the (min–max normalised) value of the \({j}^{{th}}\) characteristic for the \({i}^{{th}}\) initialisation. In this way, all values in \(T{C}^{* }\) are between \(0\) and \(1\).Effects of increasing cell typesTo identify the effect of each cell type on the response to a TMZ + ICI + OV combination treatment, we generated a virtual patient with 5000 cells according to patient T1’s cell distributions. Cells were assigned an initial location uniformly within a circular domain of the same radius as BottomRepresentative and TopRepresentative. Then, we increased the number of cells up to 10,000 in increments of 1000, varying only one cell type at a time. Additional cells were also assigned an initial location uniformly within the same domain. To analyse the variation in treatment response, we looked at the fold change in cancer cells as defined by Eq. (18). Results were stored in a table \(T\) where the element in position \((i,j)\) (i.e., \(T\left(i,j\right)\)) represented the fold change in cancer cells when \(\left(j-1\right)\cdot {1000}\) cells of type i are added to the 5000 initial cells. The table was then min–max normalised, i.e.,$$T\left(i,j\right)=\frac{T\left(i,j\right)-\mathop{\min }\limits_{i,j}T(i,j)}{\mathop{\max }\limits_{i,j}T(i,j)-\mathop{\min }\limits_{i,j}T(i,j)}.$$
(21)
Spatial computational modelling illuminates the role of the tumour microenvironment for treating glioblastoma with immunotherapies
