An Agent-Based Model of Combination Oncolytic Viral Therapy and Anti-PD-1 Immunotherapy Reveals the Importance of Spatial Location When Treating Glioblastoma

Simple Summary A combination of oncolytic viral therapy and immunotherapy provides an alternative option to the standard of care for treating the lethal brain tumor glioblastoma (GBM). Although this combination therapy shows promise, there are many unknown questions regarding how the tumor landscape and spatial dosing strategies impact the effectiveness of the treatment. Our study aims to shed light on these questions using a novel spatially explicit computational model of GBM response to treatment. Our results suggest that oncolytic viral dosing in the location of highest tumor cell density leads to substantial tumor size reduction over viral dosing in the center of the tumor. These results can help to inform future clinical trials and more effective treatment strategies for oncolytic viral therapy in GBM patients. Abstract Oncolytic viral therapies and immunotherapies are of growing clinical interest due to their selectivity for tumor cells over healthy cells and their immunostimulatory properties. These treatment modalities provide promising alternatives to the standard of care, particularly for cancers with poor prognoses, such as the lethal brain tumor glioblastoma (GBM). However, uncertainty remains regarding optimal dosing strategies, including how the spatial location of viral doses impacts therapeutic efficacy and tumor landscape characteristics that are most conducive to producing an effective immune response. We develop a three-dimensional agent-based model (ABM) of GBM undergoing treatment with a combination of an oncolytic Herpes Simplex Virus and an anti-PD-1 immunotherapy. We use a mechanistic approach to model the interactions between distinct populations of immune cells, incorporating both innate and adaptive immune responses to oncolytic viral therapy and including a mechanism of adaptive immune suppression via the PD-1/PD-L1 checkpoint pathway. We utilize the spatially explicit nature of the ABM to determine optimal viral dosing in both the temporal and spatial contexts. After proposing an adaptive viral dosing strategy that chooses to dose sites at the location of highest tumor cell density, we find that, in most cases, this adaptive strategy produces a more effective treatment outcome than repeatedly dosing in the center of the tumor.


Introduction
Glioblastoma (GBM) is the most aggressive form of brain tumor, and its poor average survival rate of 1-2 years has remained largely constant for many years, despite significant advancements in other cancer treatment modalities [1]. The standard treatment regimen for GBM consists of surgery, followed by radiotherapy and chemotherapy. Oncolytic viral therapy (OVT) is a treatment involving the administration of a virus that infects and replicates within cancer cells, causing them to lyse. Due to oncolytic viruses' ability to selectively target tumor cells over healthy cells, OVT has the potential to perform

Materials and Methods
In this work, we develop a spatially explicit hybrid cellular automaton (CA) and partial differential equations (PDE) model. We simulate the model in a three-dimensional domain, Ω = [0, L] 3 . Each site in the lattice can be occupied by at most one tumor cell, one innate immune cell, and one T cell. Note that the average diameter of a GBM cell (line U87) is 12-14 µm [16], the average diameter of a T cell is 7-10 µm, and the average diameter of a macrophages is 21 µm [17]. Due to the similar scale between all three cell types, it is reasonable to model them on the same lattice. We define the three-dimensional neighborhood of each site to be the Moore neighborhood, i.e., the nearest 26 sites.
We model the virus using a reaction-diffusion PDE, which we describe in more detail in Section 2.3. The virus is initially injected into the center of the initial population of susceptible tumor cells, which we initialize as the center of the domain. Thus, we start with nonzero viral concentration in a spherical region at the center of the tumor, with zero concentration everywhere else. We start with 1000 tumor cells and administer a dose of 10 6 pfu, which are scaled analogously from our ODE model in [18]. Each viral dose spans a region a radius of 2 sites, corresponding to 33 total sites in the model.
To give the reader an idea of the cell actions within the model, Figure 1 displays a schematic example of a two-dimensional cross-sectional portion of the ABM before and after a time step τ. The blue shades indicate viral concentration levels, and the shapes represent tumor and immune cells within the model. The events taking place within this time step include tumor and T cell division, T cell killing of tumor cells, immune cell migration, natural cell death, and immune cell consumption of viral particles.  Figure 1. An example of a two-dimensional cross-section of the ABM before (left) and after (right) a time step τ. The blue shades in the background indicate the viral concentration levels at each site, with darker shades representing higher viral concentrations. In this example, the T cell in A4 replicates and places a daughter T cell in B4, the antitumor T cell in B2 kills the susceptible tumor cell, the T cell in B3 moves to C4, the tumor cell in B4 replicates and places a daughter cell in B3, the T cell in C2 dies naturally, the antiviral T cell in C3 kills the infected cell, the innate immune cell in C4 moves to C3, the viral concentration in D2 decreases due to T cell consumption, and the viral concentration in D3 decreases due to innate immune cell consumption.

Tumor Cells
Tumor cells can become infected by the virus in a neighborhood of its location, so we first define N (x, t) to denote the Moore neighborhood of x at time t, i.e., the 26 nearest sites in the three-dimensional domain. The probability that susceptible tumor cells at site x can become infected during the time step [t, t + τ) is defined as a function of the viral concentration at that site, V(x, t), as follows: where β is the infection rate (pfu −1 h −1 ), used in the ODE model in [18], and T s (t) denotes the total number of susceptible tumor cells, T s (t) = ∑ x∈Ω T s (x, t).
For simplicity, we assume no tumor cell motility, so the tumor cells spread outward exclusively through cell division. Each new tumor cell is assigned an intrinsic cell cycle time, normally distributed with a mean of ln(2)/r t , where r t is the growth rate of tumor cells (h −1 ), and standard deviation σ cycle . At each time step, the cell cycle clock decreases by the time step length with probability P div (t), defined below. A given cell divides once its cycle clock surpasses 0. Note the mean cell cycle time, µ cycle , is converted from the fitted growth rate r t = 0.0192 h −1 in [18], and σ cycle = 3 h is an ad hoc estimate. When a susceptible tumor cell divides, it randomly chooses one of its neighboring sites, unoccupied by another tumor cell, on which to place its daughter cell. If there are no unoccupied sites, then it randomly chooses one of its 26 nearest neighbors to push outward to make room for its daughter cell, and then a chain of cells is pushed outward in the same direction. Such cell division assumptions are consistent with those used in the ABM in [19].
In order to enforce the tumor carrying capacity, we define the following probability of reducing each cell division counter at time t, as a function of the susceptible and infected tumor populations: where as above, T s (t) refers to the total number of susceptible cells, and similarly, T I (t) = ∑ x∈Ω I(x, t) denotes the total number of infected cells at time t. This assumption is consistent with those made in [13,14,20], and results from these models have shown agreement with experimental data cited in [14,21].

Immune Cells
Innate immune cells arise, i.e., become "activated" to target viral antigens, at a site x which is unoccupied by an innate immune cell, with a probability P inn . This probability depends primarily on the viral concentration at each site, as defined below. Additionally, there is a positive feedback loop of recruitment between macrophages and natural killer cells, so the probability also depends on the current population of innate immune cells in neighboring sites of x. We define where, as above, N (x, t) denotes the neighborhood of x at time t. We estimate γ i 1 and γ i 2 so that γ i 1 V(x, t) is generally larger than γ i 2 Z(N (x, t)) in the presence of the virus. We use a similar assumption to the models used in [13,14], adapted to account for the concentration of virus at each site. Adaptive immune cells, i.e., T cells, are recruited to the tumor microenvironment (TME) by the innate immune cells. At each site occupied by an innate immune cell and currently unoccupied by a T cell, there is a probability of τa AZ that a T cell will be recruited to this site, during the time interval [t, t + τ). The type of T cell, i.e., antiviral vs. antitumor, recruited to this site is determined globally. This status depends on the proportion of total susceptible cells, T s (t), to infected tumor cells, T I (t), in the domain. If T s (t) + T I (t) > 0, then the probability that a new T cell will be antiviral is and 1 − P A V (t) is the probability that the T cell will be antitumor. Note we are assuming that the tumor vasculature exists throughout the TME, so tumor infiltration by immune cells occurs via the vasculature, as described in [22]. This allows immune cells to emerge at sites within the tumor, rather than starting on the boundary of the tumor.
We also incorporate innate and adaptive immune motility within the model. Innate and adaptive immune cells move at distinct rates and randomly choose a neighboring site that is unoccupied by an immune cell of its type. The immune cells move preferentially toward the cells that they target, modeling chemotaxis. Hence, innate immune cells and antiviral T cells are m Z and m A V times more likely, respectively, to move toward sites occupied by infected tumor cells than to other sites, and antitumor T cells are m A T times more likely to move toward sites occupied by tumor cells than to other sites. Note that m Z , m A V , m A T are all larger than 1. Innate immune cells move to an empty neighboring site, i.e., one that is not currently occupied by an innate immune cell, at rate r mov Z . Thus, the probability that such an innate immune cell moves to a neighboring site in the time interval [t, t + τ) is τr mov Z . Subsequently the innate immune cell chooses an empty site to migrate to, and it is m Z times more likely to move to a site occupied by an infected cell than any other site. Similarly, the antitumor and antiviral T cells move at rates r mov AT and r mov AV , respectively.
To simulate T cell proliferation, if an antitumor or antiviral T cell and tumor cell are located at the same site, the T cell proliferates and adds a daughter cell to a randomly chosen neighboring site during time interval [t, t + τ) with probability τa AT or τa AV , respectively. Here, a AT and a AV are the rates of the tumor cell-or infected cell-mediated proliferation of T cells, respectively, from the ODE model in [18], which we have converted to probabilities in this stochastic ABM. Similarly, innate immune cells at the same site as an infected tumor cell proliferate at rate a Z , so the probability of such a proliferation event at the site of an innate immune cell and infected tumor cell There are two mechanisms of immune cell death: one is driven by a natural death rate, causing cells to die after a certain amount of time, and the second is driven by immune activity, causing cells to die after a certain number of immune "events", modeling the exhaustion of the cells. For simplicity, we assume that the immune events are restricted to the killing of tumor cells, since the killing of viral particles is difficult to measure as individual events.
If an adaptive immune cell and tumor cell simultaneously occupy the same site, i.e., an antiviral T cell with an infected tumor cell or an antitumor T cell with any tumor cell, then the tumor cell dies at rate k I A or k TA cell −1 h −1 , respectively. Hence, the probability of such an event in [t, t + τ) is τk TA or τk I A . Similarly, the probability that an innate immune cell kills an infected tumor cell at the same site in [t, t + τ) is τk I , where k TA , k I A , k I are all rates from the ODE model in [18]. Figure 2 illustrates the three-dimensional locations of infected and susceptible tumor cells, on the left, and antitumor and antiviral T cells, on the right, on days 11, immediately after Dose 2 is administered, and on Day 80, at the end of a sample simulation.
We also incorporate PD-1/PD-L1 checkpoints in the ABM, similarly to our ODE model in [18]. The PD-1/PD-L1 checkpoint implementation is described in more detail in Appendix A.1. The immune checkpoints can be blocked by anti-PD-1 immunotherapy when it has been administered. We incorporate treatment with anti-PD-1 using a global approach, i.e., we calculate the blocking rate of PD-1 using the total concentration of PD-1 in the tumor microenvironment at time t. Further detail describing this implementation can be found in Appendix A.2.

Viral Diffusion
The viral concentration will be modeled using a reaction-diffusion equation, as follows: where D V is the diffusion coefficient for the virus, and k VZ Z(x, t) and k VA Y V (x, t) represent viral death due to consumption by innate immune cells and killing by adaptive immune cells, respectively, at rates k VZ and k VA . Additionally, α I βT s V represents viral infection of susceptible tumor cells at rate β T . The term b T δ T T I (x, t) represents the viral particles that are released when an infected cell bursts at rate δ T , with burst size b T . We assume that the virus cannot leave or enter from the boundary, so the boundary conditions are: We discretize Equation (5) using central differences for spatial derivatives and forward differences for time derivatives (FTCS).

Antitumor T Cell Killing
We utilize the model to gain a better understanding of the role of T cells in response to the combination therapy, by comparing the relative importance of hte T cell-mediated tumor killing rate with tumor antigenicity. First, we investigate the impact of the tumor cell killing rate, k TA , by antitumor T cells on the tumor response to treatment. We vary the parameter k TA in the range [0.01, 0.06], around the baseline value of 1/24 ≈ 0.42 cell −1 h −1 , estimated from [23], leaving all other parameters set at their baseline rates. Figure 3 shows the susceptible tumor population at Day 80, the end of each simulation, as a function of k TA . Interestingly, we observe a nonmonotonic relationship between tumor size and T cell killing rate. As one would expect, for the lowest killing rate in the range under consideration, k TA = 0.01 cell −1 hour −1 , the T cells are not sufficiently cytotoxic to respond effectively to the treatment, leading to a significantly larger tumor after 80 days than we see with larger k TA values. The tumor decreases to similar levels for k TA = 0.02 and 0.03, with a slightly smaller population for k TA = 0.03. However, as k TA increases further, the tumor population increases as well, suggesting that increasing the T cell cytotoxicity beyond a certain threshold will not improve the efficacy of the treatment. We hypothesize that in this high k TA range, the T cells kill cancer cells too quickly early in the treatment process, which ultimately leads to a smaller number of T cells to mount an attack on the tumor. This investigation suggests an optimal value of k TA = 0.03, at which the tumor population initially increases enough to activate a sufficient level of antitumor T cells, which in turn keep the tumor relatively controlled.

Tumor-Mediated T Cell Proliferation
We compare the relative impact of the T cell killing rate from the previous section with the level of tumor antigenicity. In our model, this antigenicity level is dictated by the tumor-mediated proliferation rate of T cells, a AT , which we found to be a significant adaptive immune-related parameter in our ODE model parameter sensitivity analysis in [18]. Figure 4 displays the susceptible tumor population at the end of the simulation as a function of a AT . We observe a monotonic relationship between the final tumor size and a AT , with the tumor decreasing as the level of tumor antigenicity increases. Note that high levels a AT can lead to a small, well-controlled tumor 80 days after the start of treatment, around 450 cells for a AT = 0.08 and 260 cells for a AT = 0.09 after starting the simulations with about 1000 cells. We did not observe this significant reduction in tumor size at the optimal value for k TA . This suggests that the treatment is more effective when there is a critical mass of T cells to attack the tumor than when there is a smaller number of highly cytotoxic T cells. Building this critical mass of T cells seems to be particularly important during the early stages of the simulation, when the tumor increases rapidly until the T cells are able to mount a sufficient attack on the tumor cells.
In Figure 5, the graph on the right shows the tumor control achieved from the combination of six oncolytic viral doses administered in the center of the domain and anti-PD-1 immunotherapy, using a representative simulation with a AT = 0.08. The graph on the left shows the susceptible tumor population when oncolytic viral therapy is administered without anti-PD-1. In this case, the tumor size trends upward with short sporadic periods of decline, until it starts to level off toward the end of the simulation. Overall, we see no tumor control without the addition of anti-PD-1, and similar behavior is observed for other parameter sets. Thus, it is imperative to administer an anti-PD-1 immunotherapy in combination with OVT, so that the immune response stimulated by OVT can be sufficiently effective.

Viral Dose Timing
Thus far, there has been little focus in the literature on the timing of oncolytic viral doses, but timing can be crucial to treatment efficacy, especially when a multifaceted immune response plays a significant role, as it does in the case of OVT. Thus, we used our model to sequentially investigate the optimal timing of six viral doses. We fixed the initial viral dose at the start of the simulation, and then we varied the timing of viral dose 2, at Day 4, the approximate time of the population peak when a single viral dose is administered, at Day 10, Day 15, and Day 20. After comparing the means from 10 simulations for each option, we found that dosing at Day 10, about six days after the tumor population begins declining, produces the smallest tumor at the end of the simulation. Thus, we fixed Dose 2 at Day 10, and subsequently tested Dose 3 at Days 13, 20, 30, and 40, following the same procedure that we used to determine the timing of Dose 2. Our results suggested that it was optimal to administer Dose 3 shortly after Dose 2, on Day 13. By using a similar process for Doses 4-6, we determined that it was optimal to administer Dose 4 at Day 30, Dose 5 at Day 33, and Dose 6 at Day 60. Thus, our results suggest that the combination therapy may be more effective when pairs of viral doses are given close together, followed by a longer break before the next pair of doses. The timing of the doses is indicated in Figure 6, which displays a simulation of the tumor, with tumor antigenicity level a AT = 0.07 and all other parameters set at their baseline levels.

Density-Based Adaptive Viral Dosing
The location of oncolytic viral doses has also not been well-explored, and can be consequential for developing an effective treatment regimen. Our spatially explicit ABM proves particularly useful for investigating the effect of the viral dosing locations on the treatment response. We model the direct intratumoral delivery of the virus, which generally controls the viral concentration more effectively in a specific location than intravenous delivery [24]. Typically, the intratumoral dose is administered roughly in the center of the tumor region, which we simulated in the model by administering each viral dose in the center of the three-dimensional domain. We then compared the center dosing results with an adaptive dosing strategy, in which viral doses were administered at the locations of maximum tumor cell density at the time of each dose. We implemented this in the model by calculating the cell density in a three-dimensional neighborhood with a two-site radius around each site in the domain. We chose the site with maximum tumor cell density to administer the viral dose; if there is more than one site with maximal cell density, then we split the viral dose evenly between these sites.  (41, 56), and Dose 6 is split between the three neighboring sites at (34,73), (34,74), and (35,74). We note that even when multiple sites have neighborhoods with the same maximum tumor cell density, they are always neighboring sites in this example, and they are nearly always neighboring sites in the other realizations that we have simulated. Thus, in a clinical setting, it is reasonable to choose a single location of estimated highest tumor cell density using tumor imaging technology.
We compared the adaptive viral dosing strategy to center dosing in the range of antigenicity levels leading to a responsive or stable tumor, characterized by at least 50% tumor reduction from its maximum size, by the simulation end time, 80 days after the start of treatment. These response levels are reached with a AT = 0.07, 0.08, so we simulated 40 realizations of the ABM for each of these antigenicity levels to account for the stochasticity within the model. Figure 8 displays the mean tumor size resulting from these realizations, with the center dosing mean shown in blue and the adaptive dosing mean shown in yellow. The graph on the right is a zoomed in version of the left graph to highlight the difference in tumor size at the end of the simulation. We also investigated the benefit conferred by increasing the viral infection rate to β = 1 × 10 −7 pfu −1 h −1 , with the mean tumor size shown by the green curve in this figure. The average tumor size 80 days after the start of treatment with center dosing is 2.425 × 10 3 cells, and adaptive viral dosing yields a 24.5% reduction in this average tumor size. Increasing the viral infection rate to 1 × 10 −7 pfu −1 h −1 yields an additional 8.6% reduction in average tumor size over the center dosing case. Hence, we observe a sizable improvement in treatment response due to adaptive dosing in locations of high tumor cell density, with a marginal improvement when the virus is highly infectious.
Similarly, Figure 9 shows the tumor population means for the same cases, with a higher level of tumor antigenicity, a AT = 0.08 cell −1 h −1 . In this case, we see even more improvement resulting from the adaptive dosing strategy. The mean tumor size 80 days after the start of treatment with center dosing is 855.3 cells, and the adaptive dosing strategy reduces this mean by 33.4%, which is comparable to the reduction resulting from both adaptive dosing and increased viral infectivity in the a AT = 0.07 case. Increasing the viral infection rate when a AT = 0.08 provides an additional 8.8% reduction in tumor size over the center dosing mean. These results suggest that choosing viral dosing locations in the regions of highest cell density produces significant improvement over repeatedly dosing in the same location, and we observe from Figures 8 and 9 that this improvement appears to increase over time. Note that for a AT ≥ 0.09, both center and adaptive viral dosing produce near tumor clearance, so adaptive dosing does not provide a significant improvement when center dosing yields such a strong treatment response.
In order to make model simulations computationally tractable, we were required to scale down the initial tumor population to start around 1000 cells, which is two orders of magnitude smaller the initial tumor population we used in the ODE model in [18]. Similarly we scaled the standard viral dose used in the ODE model by 10 −2 to obtain a proportional viral dose to use in the ABM. To gauge whether we can expect comparable results in a larger tumor, we tested the a AT = 0.08 case in a tumor starting with 5000 cells, so about 5 times our typical simulation size. In a single simulation, we found the ratio of tumor sizes between center dosing and adaptive dosing to be 0.665, or a nearly identical tumor size reduction of 33.5% with adaptive dosing to the average reduction we observed with the smaller tumors. The results of this simulation are shown in Figure 10, with the graph on the right showing a zoomed in version of the full trajectory on the left. Although this tumor is still significantly smaller than a human tumor in vivo, we are encouraged by the fact that the percent size reduction scales up nearly identically from our typical simulation size in this larger example.   We also investigate the chosen locations in the adaptive dosing procedure. In order to compare the dosing locations for each dose number, we plot the mean adaptive viral dosing distances from the center of the domain for Doses 2-6, generated from 40 realizations of each parameter set, in Figure 11. We plotted these means for a AT = 0.07, 0.08 and β = 2.5 × 10 −8 , 1 × 10 −7 , and we observed that, in almost all cases, the dosing locations move farther away from the domain center as the dose number increases. Due to the growing difference between the adaptive and center dosing strategies as the dose number increases, this suggests that adaptive viral dosing may confer a more significant benefit in later doses, when the high-density regions have spread farther from the initial tumor center.

Discussion
Over the past several decades, efforts to develop new therapies for GBM tumors have intensified, yet none have significantly impacted patient mortality [25]. OVs, equipped with specific oncolytic properties, have the potential to confer therapeutic benefits to GBM patients due to their specificity, potency, tolerability, and potential to be combined to with multiple immunotherapies [26,27]. Questions remain surrounding OVT dosing strategies and concerning the precise role of the immune system in response to OVT, in a quest to increase the efficacy of OVT. In order to address these questions, we developed an agent-based model that allows us to consider immune interactions and to investigate the impact of the spatial location of viral dosing.
The agent-based approach allowed us to extend our mechanistic ODE model in [18] to consider more complex spatial interactions. We previously developed the ODE model to study the interactions between GBM cancer cells and innate and adaptive immune cells in response to OVT and anti-PD-1 therapy. A limitation of that mechanistic model was that it could not capture spatial heterogeneity inherent in GBM tumors, including the spatial distribution of various cell types and the diffusion of the virus and anti-PD-1 drug. This motivated us to develop a spatially explicit ABM of tumor response to a combination of oncolytic viral therapy and anti-PD-1 immunotherapy. ABMs are important tools in translational systems biology, due to their ability to incorporate spatial structure and stochasticity on multiple spatial and time scales, and they have become much more prevalent in recent years [9].
Both the ABM and our previous ODE model suggest that anti-PD-1 immunotherapy is necessary to allow the OVT to work effectively, with the body's T cells serving as the primary antitumor weapon. When OVT is not combined with anti-PD-1, the PD-1/PD-L1 complex prevents many T cells from targeting the tumor cells, resulting in an insufficient immune response. The ABM also confirms our finding using the ODE model, that the tumor antigenicity level is more consequential for OVT treatment response than the T cell killing rate of tumor cells. Thus, if the treatment does not yield a reduction in tumor size shortly after the start of treatment, it may be advantageous to administer an IL-2 immunotherapy in combination with the OVT and anti-PD-1, in order to increase the T cell proliferation rate. Modeling the combination of these three treatment options will be the subject of future work.
An important benefit of the ABM approach is that it allows us to explore spatial strategies for oncolytic viral dosing, which is not possible with mechanistic ODE models. We compare an adaptive dosing strategy that chooses the sites with the highest tumor cell density with a fixed dosing strategy that administers each dose in the same location throughout the treatment period. Given a sufficiently antigenic tumor, the density-based adaptive dosing strategy is more effective than the fixed dosing strategy, yielding a smaller post-treatment tumor. When using the adaptive strategy, each successive viral dosing location is chosen to be farther from the center, as the tumor spreads outward. Hence, the utility of the adaptive dosing strategy is more pronounced for doses later in the treatment schedule. We find that there is a more significant benefit conferred by this adaptive dosing strategy than from using a highly infectious virus. We note that we explored optimal dose timing before our investigation of optimal dose location, so it is possible that the adaptive dosing strategy could be improved even further with the reconsideration of optimal dose timing. In a clinical setting, the dosing locations can be determined using tumor imaging, e.g., computed tomography (CT) scans, which display relative tumor density. Such scans can be used to estimate the location of maximum tumor cell density before each viral dose. Thus, our results suggest that if clinicians have the ability to collect CT scans before each viral dose, then this practice will prove valuable, particularly for doses later in the treatment schedule.
Our ABM allowed us to more accurately approximate tumor growth and local interactions between viral particles, immune cells, and tumor cells. We observed similar overarching trends between the ODE model in [18] and the ABM, particularly with respect to the tumor-immune interactions. Additionally, the spatial features of the ABM allowed us to consider more complex viral dosing strategies. A limitation of the model is that although we did incorporate T cell exhaustion due to tumor cell killing, we did not include a limit on the number of times T cells can proliferate. While we included the PD-1/PD-L1 checkpoint within the model, we did not incorporate other checkpoints like CTLA-4, CD28, and TIM3, which can also impact T cell proliferation. This may result in the model overestimation of the T cell population within the TME. Our model can be easily adapted to explore the impact of additional checkpoints, so in future work, we would like to incorporate these checkpoints and to investigate the effect of a proliferation limit for T cells on the response to this combination therapy. We hope to use this work as motivation to design experiments that test adaptive immune cell killing and proliferation rates and exhaustion-related limits on these processes. This will help to provide further clarity about the feasible ranges for these parameters in the model. Additionally, the model developed here is parameterized using data from mouse models, so it cannot be directly translated to human patients. We would like to validate our computational results by designing future murine experiments that compare the adaptive and fixed viral dosing strategies, for example, by administering the same combination therapy to each mouse, with the oncolytic virus administered in the tumor center in half of the mice, and the virus administered adaptively in the site of estimated maximum tumor cell density in the other half. The proportional difference between tumor sizes after 80 days can then be compared to our ABM results. We hope that the results will motivate a subsequent clinical trial in GBM patients. In addition, we plan to model the administration of a new type of oncolytic virus that is engineered to express a PD-L1 inhibitor, which provides an opportunity to combine OVT and immunotherapy in a single treatment [28]. The ABM framework is flexible enough to investigate novel therapies such as this one, allowing us to adapt our current ABM for this future work.

Conclusions
In this study, we developed an ABM to model the response of GBM to a combination of oncolytic viral therapy and anti-PD-1 immunotherapy. Our model simulations suggest that the tumor antigenicity level, which impacts T cell proliferation, is more consequential for treatment efficacy than the T cell killing rate. In addition, we determine an optimal viral dosing schedule and consider an alternative spatial dosing strategy. We show that an adaptive viral dosing strategy that chooses to dose in the locations of highest tumor cell density is substantially more effective than administering each dose in the same location in the center of the tumor. These results suggest that collecting CT scans during treatment for GBM patients can help to inform the location of oncolytic viral dosing, thereby improving the effectiveness of this treatment.

Data Availability Statement:
The authors confirm that the data supporting the findings of this study are available within the article. MATLAB code for model generation and simulation can be found at https://github.com/kstorey90/ABM_OVT_antiPD1 (accessed on 10 October 2021).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
For computational feasibility, we initialized the tumor population with about 1000 cells, corresponding to a radius of six sites at the center of the lattice domain. HSV particles have diameters ranging from 155 to 240 nm [29], so we approximate the diameter as 200 nm, which converts to 2 × 10 −5 cm. Using our baseline estimate of l = 0.0014 cm for site length, the initial viral dose spans a radius of two sites, corresponding to 33 total sites.

Appendix A.1. PD-1/PD-L1 Checkpoints
To incorporate PD-1/PD-L1 checkpoints within the model, we multiply the probabilities of the innate immune-mediated activation of adaptive immune cells and of tumor cell-mediated T cell proliferation by a factor F(P, L). This factor represents the suppression of T cell activation and proliferation via the PD-1/PD-L1 checkpoint, where P, L denote the molar concentrations of PD-1 and PD-L1, respectively, expressed by cells within the model. The molar concentrations (in µmol/cm 2 ) are obtained by first calculating the PD-1 expression on all T cells and the PD-L1 expression on all T cells, tumor cells, and innate immune cells, as outlined below. As P and L increase, so does the number of PD-1/PD-L1 complexes within the tumor region. This increase corresponds to a smaller F(P, L) value, modeling the inhibition of T cell activity.
In [18], we calculate estimates for the parameters used in the factor F(P, L). We estimate the concentration of PD-1 and PD-L1 per T cell to be ρ P = 1.259 × 10 −11 µM and ρ L = 2.510 × 10 −11 µM, respectively.

Appendix A.2. Anti-PD-1
In the model, we track the concentration of anti-PD-1 in the tumor microenvironment, using µmol/cm 3 .
For the decay rate of anti-PD-1, we use δ A in the range 0.0019-0.01 h −1 , from [30,31]. Just as in [18], we estimate the anti-PD-1 blocking rate of PD-1, µ PA , by assuming that 10% of the anti-PD-1 drug is used in blocking PD-1 and that 90% degrades naturally. Thus, so we assume that for a steady stateP and δ A = 0.0019, In [18], we derived the following conversion from the dosage, D, of nivolumab, in mg/kg, to plasma concentration C max , in µM = µmol/cm 3 : A typical nivolumab regimen consists of a single intravenous dose of 3 mg/kg nivolumab, administered for one hour, once every two weeks. We incorporate this treatment regimen within our model, by adding nivolumab every two weeks, at concentration C max (3 mg/kg) = 0.481 µmol/cm 3 .