# A hybrid discrete-continuum modelling approach for the interactions of the
immune system with oncolytic viral infections††thanks: Corresponding author:
David Morselli<EMAIL_ADDRESS><EMAIL_ADDRESS>
This research was partially supported by the Italian Ministry of Education,
University and Research (MIUR) through the “Dipartimenti di Eccellenza”
Programme (2018-2022) – Dipartimento di Scienze Matematiche “G. L. Lagrange”,
Politecnico di Torino (CUP: E11G18000350001). MED and DM are members of GNFM
(Gruppo Nazionale per la Fisica Matematica) of INdAM (Istituto Nazionale di
Alta Matematica). ALJ acknowledges the Australian Research Council (ARC)
Discovery Project (DP) DP230100025. We also acknowledge the support of the
Australian National Health and Medical Research Council, through grant NHMRC
IDEAS 2013058. Part of this work was performed on the OzSTAR national facility
at Swinburne University of Technology. The OzSTAR program receives funding in
part from the Astronomy National Collaborative Research Infrastructure
Strategy (NCRIS) allocation provided by the Australian Government, and from
the Victorian Higher Education State Investment Fund (VHESIF) provided by the
Victorian Government.

David Morselli Department of Mathematical Sciences “G. L. Lagrange”,
Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino,
ItalyDepartment of Mathematics, School of Science, Computing and Engineering
Technologies, Swinburne University of Technology, John St, 3122, Hawthorn,
VIC, AustraliaDepartment of Mathematics “G. Peano”, Università di Torino, Via
Carlo Alberto 10, 10124 Torino, Italy Marcello Edoardo
Delitala,11footnotemark: 1 Adrianne L. Jenner School of Mathematical Sciences,
Queensland University of Technology, George St, 4000, Brisbane, QLD, Australia
Federico Frascoli22footnotemark: 2

###### Abstract

Oncolytic virotherapy, utilizing genetically modified viruses to combat cancer
and trigger anti-cancer immune responses, has garnered significant attention
in recent years. In our previous work, we developed a stochastic agent-based
model elucidating the spatial dynamics of infected and uninfected cells within
solid tumours. Building upon this foundation, we present a novel stochastic
agent-based model to describe the intricate interplay between the virus and
the immune system; the agents’ dynamics are coupled with a balance equation
for the concentration of the chemoattractant that guides the movement of
immune cells. We formally derive the continuum limit of the model and carry
out a systematic quantitative comparison between this system of PDEs and the
individual-based model in two spatial dimensions. Furthermore, we describe the
traveling waves of the three populations, with the uninfected proliferative
cells trying to escape from the infected cells while immune cells infiltrate
the tumour.

Simulations show a good agreement between agent-based approaches and numerical
results for the continuum model. Some parameter ranges give rise to
oscillations of cell number in both models, in line with the behaviour of the
corresponding nonspatial model, which presents Hopf bifurcations.
Nevertheless, in some situations the behaviours of the two models may differ
significantly, suggesting that stochasticity plays a key role in the dynamics.
Our results highlight that a too rapid immune response, before the infection
is well-established, appears to decrease the efficacy of the therapy and thus
some care is needed when oncolytic virotherapy is combined with immunotherapy.
This further suggests the importance of clinically improving the modulation of
the immune response according to the tumour’s characteristics and to the
immune capabilities of the patients.

Keywords— Oncolytic virus, Immunotherapy, Individual-based models, Continuum
models, Bifurcation analysis

MSC Classification: 35Q92, 92-08, 37N25, 37G15

##  1 Introduction

Oncolytic viruses are able to infect and kill cancer cells, while mostly
sparing healthy tissues [6, 19, 36, 42, 59]. Despite their high potential as
targeted cancer therapy, it appears now clear that, in most of the cases,
virotherapy is unable to eradicate a tumour alone; hence, most of the current
efforts are devoted towards its combination with other therapies [49]. One of
the most promising of such combinations is with immunotherapies [17], which
has been tested in several clinical trials (such as Ref. [3, 13]; we refer to
Ref. [17] for a more comprehensive review). The “avoidance of immune
destruction” is one of the hallmarks of cancer [24], therefore therapies that
contribute to the activation of the immune system may play a central role in
keeping a tumour under control and, if possible, in eradicating it. The
interplay between oncolytic viruses and immune cells is two-folded: oncolytic
viruses are able to stimulate immune cells, not only against viral particles,
but also against tumour cells; on the other hand, an immune response that
targets the oncolytic virus may prevent an effective infection in the whole
tumour, making virotherapy inefficient [61]. The complexities of these
dynamics motivate the use of mathematical models to gain a deeper
understanding, with the goal of suggesting optimal treatment schedules for the
combination of virotherapy and immunotherapies.

Several mathematical models have previously been adopted for the study of the
interactions between oncolytic viruses and immune system, including ordinary
differential equations (ODEs) [2, 15, 63, 67, 69], partial differential
equations (PDEs) [20, 38, 43, 71], stochastic agent-based models [62, 64] and
hybrid discrete-continuous multi-scale models [32]. While most of them
restrict their attention to the systemic immune response, some others also
explicitly model immunotherapies, such as immune checkpoint inhibitors [20,
38, 63, 62, 64, 67] and chimeric antigens receptor T-cells (CAR-T) [48]. In
general, individual-based models track individual cells, making it possible to
represent processes happening at single cell-scale and easily include
stochasticity; in this context, continuous fields are often used to model
molecular elements, such as nutrients, leading to a _hybrid_ modelling
approach. In contrast, deterministic continuum models describe volume
fractions, hence the biological interpretation of the terms comprised in the
model equations is less straight-forward and stochasticity cannot be included
easily; on the other hand, this approach is particularly suitable to deal with
large cell numbers for long time scales, as numerical simulations are faster
than in the case of agents-based models and sometimes analytical results may
be obtained. In order to combine the benefits of the two modelling approaches
and gain a more comprehensive understanding of the biological system under
study, a standard method is to derive a continuum macroscopic model from the
underlying discrete or hybrid stochastic model (see, for example, Ref. [8, 34,
46, 47, 57]; we refer to the introduction of Ref. [9] for a more comprehensive
literature review). The observation of significantly different behaviours of
the two modelling instances would then suggest that stochasticity plays a key
role in the phenomenon under investigation.

In our previous work [51], we have adopted this approach to model the
infection of tumour cells due to oncolytic viruses in absence of an immune
response, taking into account two alternative sets of rules governing cell
movement. Our results in the case of unrestricted cell movement show partial
tumour remission for parameter values within the biologically meaningful
range. The goal of the present work is to analyse the impact of the immune
system in this situation, with the aim to determine whether eradication or
long term control of the tumour are attainable at least in the absence of
relevant physical constraints.

Immune interactions with a tumour involve several different types of immune
cells, which are stimulated and inhibited by a large number of molecules. An
accurate description of these processes goes beyond the scopes of the present
work. In order to facilitate some theoretical understanding of the model, we
restrict our attention to a single type of immune cell, namely cytotoxic
T-cells, with the ability to kill both infected and uninfected tumour cells.
We then assume that tumour cells secrete chemoattractant and immune cells
follow the chemotactic stimuli towards the tumour (see Ref. [56] and the
references therein); this leads to a hybrid and multiscale modelling approach.
Although the derivation of this kind of model from microscopic rules is well-
known (see Ref. [1] for the specific case of immune interactions with cancer
and Ref. [7, 10] for more general situations), we are not aware of any other
work comparing agent-based and continuous models for the interactions between
immune system and oncolytic viruses.

We consider a tumour with poor immune infiltration (i.e., a _cold_ tumour in
the classification of Ref. [22]) and assume that the infection by the
oncolytic virus induces an immune anti-tumour response by increasing immune
cell inflow and improving immune recognition. We also assume that the immune
killing rate can be enhanced (e.g., by inhibition of the PD-1 and PD-L1
checkpoints [29]) and we evaluate its consequences on the therapy. First, the
spatially homogeneous ODE is considered, revealing that some parameter regions
give rise to stable limit cycles: this is not surprising, as the same
behaviour is also observed in similar models describing interactions of cancer
with immune cells [14], oncolytic virus [5, 30, 31, 58] and both together
[16]. Then, the effects of the oscillations are explored in the spatial
models: in some situations we observe the extinction of infected agents even
though the continuous model show recurrence of infection. Overall, our results
suggest that the enhancement of the immune response may either increase or
decrease the effectiveness of oncolytic virotherapy, depending on the time and
location of the viral injection.

The article is organised as follows. In Section 2, we introduce the agent-
based model and present its continuum counterpart (a formal derivation is
presented in Appendix A). In Section 3, we study the equilibria of the
spatially homogeneous ODE and the emergence of a stable limit cycle; this
analysis provides some insights on the oscillations observed in the full
system. In Section 4, we compare the results of numerical simulations of the
agent-based model and the numerical solutions of the corresponding PDEs,
comparing it with the situation in which the immune response is negligible. In
Section 5, we discuss the main findings and provide some suggestions for
future research.

##  2 Description of the agent-based model and formal derivation of the
corresponding continuum model

###  2.1 Agent-based model

We extend and improve upon the modeling framework presented in Ref. [51] by
including immune cells, which are described as agents that occupy a position
on a discrete lattice in the same way as cancer cells. We also consider a
chemoattractant secreted by cancer cells that guides the movement of immune
cells; its concentration is described as a discrete non-negative function.
Observe that we are using a hybrid discrete-continuous modeling framework,
since the chemoattractant concentration is the discretisation of a continuous
function. For ease of presentation, in this section we restrict our attention
to one spatial dimension, but there would be no additional difficulty in
considering higher spatial dimensions. In the following sections we mainly
focus on two spatial dimensions, so in Remark 2.1 we explain the small
difference of this situation.

Let us consider the temporal discretisation $t_{n}=\tau n$, with
$n\in\mathbb{N}_{0}$, $0<\tau\ll 1$, and the spatial domain
$\Omega\subseteq\mathbb{R}$ with discretisation $x_{j}=\delta j$, with
$j\in\mathbb{Z}$, $0<\delta\ll 1$; we assume $\tau$ to be small enough to
guarantee that all the probabilities defined hereafter are smaller than $1$.
We denote the number of immune cells, uninfected and infected cancer cells
that occupy position $x_{j}$ at time $t_{n}$ respectively by $Z_{j}^{n}$,
$U_{j}^{n}$ and $I_{j}^{n}$; the corresponding densities are

$z_{j}^{n}\coloneqq\frac{Z_{j}^{n}}{\delta},\qquad
u_{j}^{n}\coloneqq\frac{U_{j}^{n}}{\delta},\qquad
i_{j}^{n}\coloneqq\frac{I_{j}^{n}}{\delta}$

We then denote by $\phi_{j}^{n}$ the concentration of chemoattractant at time
$t_{n}$ and position $x_{j}$. Since the spatio-temporal scales for the
chemoattractant’s dynamics are very different from cellular ones, we describe
them with a deterministic discrete balance equation, as in Ref. [1, 12]. Table
1 summarises the variables of the hybrid agent-based model and their
macroscopic counterparts; Fig. 1 summarises the rules governing the dynamics
of the agent-based model. Cancer cells proliferate, move, become infected and
die in the same way as in Ref. [51]. The dynamics of the chemoattractant and
of the immune cells represent a novelty with respect to our previous work and
resemble some other models in the literature, as explained in the following.
We assume that the infection stimulates the immune system by increasing the
number of immune cells in the area and guiding them towards infected cells;
once an immune cell comes in contact with a cancer cell, it is able to kill it
even if it is not infected.

Quantity | Microscopic variable [Units] | Macroscopic variable [Units]
---|---|---
uninfected cancer cells |  $U_{j}^{n}$ [cells] |  $u(t,x)$ [cells/mm2]
infected cancer cells |  $I_{j}^{n}$ [cells] |  $i(t,x)$ [cells/mm2]
immune cells |  $Z_{j}^{n}$ [cells] |  $z(t,x)$ [cells/mm2]
chemoattractant |  $\phi_{j}^{n}$ [$\mu$g/mm2] |  $\phi(t,x)$ [$\mu$g/mm2]
Table 1: List of the variables for both approaches, with their units of
measurement. Figure 1: Schematic representation of the rules governing cell
dynamics in the stochastic models. Uninfected cells are represented in blue,
infected cells in red and immune cells in green. Uninfected cells may
proliferate or die according to the total density, move, become infected upon
contact with infected cells and die upon contact with immune cells. Infected
cells may move, die with constant probability and die upon contact with immune
cells. Immune cells may enter in the domain, move with the probabilities given
in Eq. (2.4) and die with constant probability.

#### Basic dynamics of uninfected cancer cells

Proliferation, undirected movement and infection of uninfected cancer cells
follow the same rules described in Ref. [51], which we recall here for the
sake of completeness.

We let an uninfected cell that occupies position $x_{j}$ at time $t_{n}$
reproduce with probability $\tau G(\rho_{j}^{n})_{+}$, die with probability
$\tau G(\rho_{j}^{n})_{-}$, and remain quiescent with probability $1-\tau
G(\rho_{j}^{n})_{+}-\tau G(\rho_{j}^{n})_{-}=1-\tau\lvert
G(\rho_{j}^{n})\rvert$. Here $\rho_{j}^{n}\coloneqq u_{j}^{n}+i_{j}^{n}$ is
the total cell density at time $t_{n}$ and position $x_{j}$ and

$G(\rho)=p\Bigl{(}1-\frac{\rho}{K}\,\Bigr{)}$ (2.1)

where $p>0$ is the maximal duplication rate and $K>0$ is the carrying
capacity.

We consider undirected random movement and assume that an uninfected cell
moves to an adjacent lattice point with probability $\theta/2$, where
$\theta\in[0,1]$, and remains at its initial position with probability
$1-\theta$.

We do not model explicitly the dynamics of the oncolytic virus and we assume
that an uninfected cell that occupies position $x_{j}$ at time $t_{n}$ becomes
infected upon contact with infected cells with probability $\tau\beta
i_{j}^{n}/K$, where $K$ is the carrying capacity and $\beta>0$ is a constant
infection rate. This approach has been commonly used for nonspatial models of
oncolytic viruses [40, 53] and allows us to avoid the complexities of viral
diffusion in the tumour microenvironment [33, 70]; we refer to the
Introduction of Ref. [51] for further details.

#### Basic dynamics of infected cancer cells

Undirected movement and death of infected cells again follow the same rules
described in Ref. [51]. We assume that the infection does not affect the cell
motility and so the probabilities are the same as the uninfected cells. We
also assume that at every time step an infected cell may die because of lysis
with probability $\tau q$, where $q>0$ is a constant death rate. We assume
that the viral replication process hijacks the cells proliferation machinery
and hence infected cells are unable to proliferate

#### Dynamics of the chemoattractant

We assume that uninfected and infected cells produce chemoattractant, at rates
$\gamma_{\phi}$ and $\alpha_{\phi}$ respectively. We choose their values so
that $\alpha_{\phi}\gg\gamma_{\phi}>0$, in line with our assumption that the
tumour is initially _cold_ and the infection by the oncolytic virus is enough
to induce an immune anti-tumour response, as often observed in vivo and in
vitro [22]. Chemoattractant density cannot grow unlimited, therefore it
saturates at $\phi^{*}>0$. The chemoattractant also decays at rate
$q_{\phi}>0$ and diffuses. The resulting balance equation is

$\phi_{j}^{n+1}=\phi_{j}^{n}+\tau
D_{\phi}\frac{\phi_{j+1}^{n}+\phi_{j-1}^{n}-2\phi_{j}^{n}}{\delta^{2}}+\tau(\alpha_{\phi}i_{j}^{n}+\gamma_{\phi}u_{j}^{n})(\phi_{j}^{n}-\phi^{*})-\tau
q_{\phi}\phi_{j}^{n}$ (2.2)

where $D_{\phi}>0$ is the diffusion coefficient. This equation closely
resembles the ones used in Ref. [1, 7] to model the evolution of a
chemoattractant concentration.

#### Dynamics of immune cells

We assume that there is a constant influx of immune cells into the
microenvironment independent of the presence of cancer cells. In addition to
this, we assume that infection by the oncolytic virus stimulates an immune
response in the whole tumour. Hence, an immune cell appears at point $x_{j}$
at time step $t_{n}$ with probability $\tau\delta S_{j}^{n}$, given by

$S_{j}^{n}=\biggl{(}S_{0}+\alpha_{z}\sum_{h}I_{h}^{n}\biggr{)}\mathbbm{1}_{\omega}(x_{j})$
(2.3)

where $\mathbbm{1}_{\omega}$ is the indicator function of the set
$\omega\subset\Omega$, $S_{0}>0$ is the base inflow rate and $\alpha_{z}>0$ is
the additional inflow rate due to the infection. In principle we could vary
$\omega$ to model the fact that some areas of the tumour are harder to reach
for immune cells (e.g. due to poor vascularisation), although this goes beyond
the scope of the present work. It is important to observe that the increase of
the inflow due to infected cells is nonlocal, as in Ref. [1]; this resembles
the recruitment of immune cells from adjacent lymph nodes and the subsequent
arrival through blood vessels.

We then assume that an immune cell that occupies position $x_{j}$ at time
$t_{n}$ moves to the lattice point $x_{j\pm 1}$ with probability $F_{j\to j\pm
1}^{n}$ and remains at its initial position with probability $1-F_{j\to
j-1}^{n}-F_{j\to j+1}^{n}$. We include both undirected, random movement and
chemotactic movement up the gradient of the chemoattractant: this latter part
depends on the difference between the chemoattractant concentration at the
initial position of the cell and the concentration of chemoattractant at the
target point. We therefore set

$F_{j\to j\pm 1}^{n}\coloneqq\frac{\theta_{z}}{2}+\nu\frac{(\phi_{j\pm
1}^{n}-\phi_{j}^{n})_{+}}{2\phi^{*}}$ (2.4)

where $z_{+}\coloneqq\max\\{z,0\\}$, $\phi^{*}$ is the saturation density of
the chemoattractant and $\theta_{z},\nu\in[0,1]$ with $\theta_{z}+\nu<1$.
Observe that, if $0\leq\phi_{j}^{n}\leq\phi^{*}$ for every $j$, then all the
probabilities are between $0$ and $1$. This kind of reasoning and the
probabilities associated have already been employed in Ref. [1, 7].

Finally, we assume that at every time step an immune cell dies with
probability $\tau q_{z}$, where $q_{z}>0$ is a constant death rate.

#### Cytotoxic action of the immune cells

We assume that cancer cells may be killed by the cytotoxic action of immune
cells upon contact; this happens at a rate proportional to the density of
immune cells. To be precise, a cancer cell that occupies position $x_{j}$ at
time $t_{n}$ dies with probability $\tau\zeta z_{j}^{n}/K$, where $K$ is the
carrying capacity and $\zeta>0$ is a constant killing rate. For the sake of
simplicity we assume that the killing rate is the same for every cancer cell,
although it could make sense to consider situations in which infected cells
are more easily recognised by immune cells and thus are killed at a higher
rate. This process is analogous to the infection of cancer cells described
above.

###  2.2 Corresponding continuum model

Letting $\tau,\delta\to 0$ in such a way that
$\frac{\delta^{2}}{2\tau}\to\tilde{D}$ and assuming that there are the
functions $u\in C^{2}([0,+\infty)\times\mathbb{R})$ such that
$u_{j}^{n}=u(t_{n},x_{j})$, $i\in C^{2}([0,+\infty)\times\mathbb{R})$ such
that $i_{j}^{n}=i(t_{n},x_{j})$, $z\in C^{2}([0,+\infty)\times\mathbb{R})$
such that $z_{j}^{n}=z(t_{n},x_{j})$ and $\phi\in
C^{2}([0,+\infty)\times\mathbb{R})$ such that $\phi_{j}^{n}=\phi(t_{n},x_{j})$
we formally obtain (see Appendix A) the following system of reaction-diffusion
PDEs

$\begin{cases}\partial_{t}u(t,x)=D\partial_{xx}^{2}u(t,x)+pu(t,x)\Bigl{(}1-\dfrac{u(t,x)+i(t,x)}{K}\Bigr{)}-\dfrac{\beta}{K}u(t,x)i(t,x)-\dfrac{\zeta}{K}u(t,x)z(t,x)\\\\[8.0pt]
\partial_{t}i(t,x)=D\partial_{xx}^{2}i(t,x)+\dfrac{\beta}{K}u(t,x)i(t,x)-qi(t,x)-\dfrac{\zeta}{K}i(t,x)z(t,x)\\\\[8.0pt]
\partial_{t}z(t,x)=D_{z}\partial_{xx}^{2}z(t,x)-\dfrac{\chi}{\phi_{\text{max}}}\partial_{x}(z(t,x)\partial_{x}\phi(t,x))-q_{z}z(t,x)+S(t,x)\\\\[8.0pt]
\partial_{t}\phi(t,x)=D_{\phi}\partial_{xx}^{2}\phi(t,x)+(\alpha_{\phi}i(t,x)+\gamma_{\phi}u(t,x))\,(\phi(t,x)-\phi^{*})-q_{\phi}\phi(t,x)\end{cases}$
(2.5)

where $D\coloneqq\theta\tilde{D}$, $D_{z}\coloneqq\theta_{z}\tilde{D}$,
$\chi\coloneqq\nu\tilde{D}$ and

$S(t,x)\coloneqq\biggl{(}S_{0}+\alpha_{z}\int_{\Omega}i(t,y)\,\mathrm{d}y\biggr{)}\mathbbm{1}_{\omega}(x)$

The first two equations are the ones studied in Ref. [51] with the addition of
the death term related to the immune system; we therefore expect to recover
similar results for small $\zeta$. This system resembles some of the models
discussed in Ref. [56] for the interactions between cancer and different kinds
of immune cells, with the relevant differences being that one of our equations
is integro-differential (as in Ref. [1]) and that the infection significantly
affects the dynamics, spatially and temporally.

In the next Section we consider the two-dimensional radially equivalent
version of this problem. Hence, we assume that

$\omega\coloneqq\Set{\boldsymbol{x}\in\Omega}{\lvert\boldsymbol{x}\rvert\leq
R}$ (2.6)

with $R>0$; this corresponds to the situation of a well-vascularised tumour in
which immune cells can easily reach any point of the domain or that of a solid
tumour that is easily accessible by the immune system both from the
hystological and topological point of view. The system of PDEs then becomes

$\begin{cases}\partial_{t}u=D\dfrac{1}{r}\partial_{r}(r\,\partial_{r}u)+pu\Bigl{(}1-\dfrac{u+i}{K}\Bigr{)}-\dfrac{\beta}{K}ui-\dfrac{\zeta}{K}uz\\\\[8.0pt]
\partial_{t}i=D\dfrac{1}{r}\partial_{r}(r\,\partial_{r}i)+\dfrac{\beta}{K}ui-
qi-\dfrac{\zeta}{K}iz\\\\[8.0pt]
\partial_{t}z=D_{z}\dfrac{1}{r}\partial_{r}(r\,\partial_{r}z)-\dfrac{\chi}{\phi_{\text{max}}}\dfrac{1}{r}\partial_{r}(rz\,\partial_{r}\phi)-q_{z}z+S\\\\[8.0pt]
\partial_{t}\phi=D_{\phi}\dfrac{1}{r}\partial_{r}(r\,\partial_{r}\phi)+(\alpha_{\phi}i+\gamma_{\phi}u)\,(\phi-\phi^{*})-q_{\phi}\phi\end{cases}$
(2.7)

with

$S(t,r)\coloneqq\biggl{(}S_{0}+2\pi\alpha_{z}\int_{0}^{r}i(t,s)s\,\mathrm{d}s\biggr{)}\mathbbm{1}_{[0,R]}(r)$

######  Remark 2.1.

When the spatial domain is the two-dimensional real plane $\mathbb{R}^{2}$
instead of the one-dimensional real line $\mathbb{R}$, the scalar index
$j\in\mathbb{Z}$ should be replaced by the vector
$\boldsymbol{j}=(j_{x},j_{y})\in\mathbb{Z}^{2}$ and the probability that a
cell moves to one of the four neighbouring lattice points is $\theta_{k}/4$,
with $k=u,i$. We then need to scale $\tau$ and $\delta$ in such a way that
$\frac{\delta^{2}}{4\tau}\to\tilde{D}$.

##  3 Corresponding ODE model and bifurcation analysis

Before comparing the agent-based and the continuous model, it is useful to
consider a homogeneous spatial configuration and analyse the equilibria of the
corresponding ODE model and their stability. The chemoattractant has the sole
purpose of guiding immune cells, therefore it can be neglected in this
nonspatial model. Hence, we now consider the system

$\begin{cases}\dfrac{\,\mathrm{d}u}{\,\mathrm{d}t}=pu\Bigl{(}1-\dfrac{u+i}{K}\Bigr{)}-\dfrac{\beta}{K}ui-\dfrac{\zeta}{K}uz\\\\[8.0pt]
\dfrac{\,\mathrm{d}i}{\,\mathrm{d}t}=\dfrac{\beta}{K}ui-
qi-\dfrac{\zeta}{K}iz\\\\[8.0pt] \dfrac{\,\mathrm{d}z}{\,\mathrm{d}t}=\alpha
i-q_{z}z+S_{z}\end{cases}$ (3.1)

It is important to observe that the inflow of immune cells in Eq. (2.5)
depends on the total number of infected cells and not just on the local
infected cell density. If we consider that $u$, $i$ and $z$ are homogeneous in
the spatial domain $\Omega$ and $\omega=\Omega$, then

$\mathbbm{1}_{\omega}(x)\int_{\Omega}i(t,y)\,\mathrm{d}y=\int_{\Omega}i(t)\,\mathrm{d}y=\lvert\Omega\rvert\,i(t)$

Hence, in this situation the parameter $\alpha$ in Eq. (3.1) corresponds to
the parameter $\alpha_{z}$ of Eqs. (2.5) and (2.3) multiplied by the measure
of the set $\Omega$ (denoted by $\lvert\Omega\rvert$).

The equilibria are $(0,0,\frac{S_{z}}{q_{z}})$, $(K-\frac{\zeta
S_{z}}{pq_{z}},0,\frac{S_{z}}{q_{z}})$, $(u^{*},i^{*},z^{*})$ and
$(0,-\frac{qq_{z}K}{\alpha\zeta}-\frac{S_{z}}{\alpha},-\frac{qK}{\zeta})$. The
latter exists only for $\alpha,\zeta\neq 0$; it is always negative, so we can
neglect it. The third one is defined by the expressions

$u^{*}\coloneqq\frac{qK}{\beta}+\zeta z^{*},\qquad
i^{*}\coloneqq\frac{Kpq_{z}(\beta-q)-S_{z}\beta(\zeta+\frac{p}{\beta}\zeta)}{\beta[q_{z}(\beta+p)+\alpha(\zeta+\frac{p}{\beta}\zeta)]},\qquad
z^{*}\coloneqq\frac{\alpha}{q_{z}}i^{*}+\frac{S_{z}}{q_{z}}$ (3.2)

When $\alpha=S_{z}=0$ we recover the equilibria in absence of the immune
response (see also Ref. [51]). As $\alpha$ and $S_{z}$ increase, $u^{*}$
increases and $i^{*}$ decreases. Similarly, when $\zeta=0$ the equilibria are
analogous to the situation without immune response (although the value of $z$
at the equilibrium may not be $0$).

The equilibrium $(0,0,\frac{S_{z}}{q_{z}})$ has eigenvalues
$(-q_{z},p-\frac{S_{z}\zeta}{Kq_{z}},-q-\frac{S_{z}\zeta}{Kq_{z}})$ and the
equilibrium $(K-\frac{\zeta S_{z}}{pq_{z}},0,\frac{S_{z}}{q_{z}})$ has
eigenvalues
$(-q_{z},-p+\frac{S_{z}\zeta}{Kq_{z}},\beta-q-\frac{S_{z}(\beta\zeta+p\zeta)}{Kpq_{z}})$.
The first equilibrium is stable when

$Kpq_{z}<S_{z}\zeta$

corresponding to the situation in which the uninfected cell density of the
second equilibrium is negative. This means that the immune system alone may be
able to eradicate the tumour without the need of any oncolytic virus. The
second equilibrium is stable in case neither the first equilibrium is stable
nor $i^{*}>0$. In this case the oncolytic virotherapy is not effective and the
outcome of the therapy depends entirely on the immune response. Let us observe
that the density of uninfected cells $K-\frac{\zeta S_{z}}{pq_{z}}$ is
increasing in the parameters $\zeta$, $S_{z}$ and decreasing in $p$, $q_{z}$:
while a complete tumour eradication is unattainable, we may still keep the
tumour at an acceptable size if the immune response is strong enough.

Parameter | Description | Value [Units] | Reference
---|---|---|---
$p$ | maximal duplication rate of uninfected cells |  $1.87\times 10^{-2}$ [h-1] | [35]
$q$ | death rate of infected cells |  $8.34\times 10^{-3}$ [h-1] | model estimate
$D$ | diffusion coefficients of cancer cells |  $1.88\times 10^{-4}$ [mm2/h] |  estimate based on [37]
$K$ | tissue carrying capacity in two dimensions |  $10^{4}$ [cells/mm2] | [45]
$\beta$ | infection rate |  $1.02\times 10^{-1}$ [h-1] |  estimate based on [21]
$D_{\phi}$ | diffusion coefficients of chemoattractant |  $3.33\times 10^{-2}$ [mm2/h] | [50]
$\phi^{*}$ | saturation density of chemoattractant |  $2.92$ [$\mu$g/mm2] | [23]
$\alpha_{\phi}$ | secretion of chemoattractant by infected cells |  $2.50\times 10^{-4}$ [mm2/(h$\cdot$cells)] | model estimate
$\gamma_{\phi}$ | secretion of chemoattractant by uninfected cells |  $5.00\times 10^{-6}$ [mm2/(h$\cdot$cells)] | model estimate
$q_{\phi}$ | decay of chemoattractant |  $8.33\times 10^{-2}$ [h-1] | [12]
$D_{z}$ | diffusion coefficients of immune cells |  $4.20\times 10^{-3}$ [mm2/h] | [1]
$\chi$ | chemotactic coefficient of immune cells |  $1.65\times 10^{-1}$ [mm2/h] | model estimate
$q_{z}$ | death rate of immune cells |  $7.50\times 10^{-3}$ [h-1] | [25]
$S_{0}$ | base inflow rate of immune cells |  $5.00\times 10^{-2}$ [cells/(mm${}^{2}\cdot$h)] | model estimate
$\alpha_{z}$ | inflow rate of immune cells due to the infection |  $3.75\times 10^{-5}$ [(mm${}^{2}\cdot$h)-1] | model estimate
$\zeta$ | immune killing rate of cancer cells |  $0.50$ or $5.00$ [h-1] | model estimate
$R_{u}$ | initial radius of uninfected cells |  $2.60$ [mm] | [37]
$R_{i}$ | initial radius of infected cells |  $1.00$ [mm] | model estimate
Table 2: Reference parameter set used in this work.

The expressions for the eigenvalues of the third equilibrium are more
complicated. Numerical simulations show that in the parameter region where
$i^{*}>0$ either this equilibrium is stable or there appears a stable limit
cycle. Fig. 2 shows numerical bifurcation diagrams for the parameters
$\alpha$, $\beta$ and $\zeta$, which present a Hopf bifurcation; in these
simulations the other parameters of Eq. (3.1) are set to the values of Table
2. The size of the oscillations of the limit cycle increases as $\alpha$ and
$\zeta$ increase, and decreases as $\beta$ increases. As a consequence, the
enhancement of the immune response may significantly decrease the
effectiveness of the therapy. However, it is fundamental to also consider that
in some cases the oscillations have a minimum very close to zero: if we take
into account a discrete number of cells, they may go extinct when approaching
the minimum due to stochastic events and the following regrowth may not take
place. Variations of the parameter $p$ never result into a bifurcation of this
equilibrium (at least when all the other parameters are within the range of
our interest), but the size of the oscillations during the convergence
decreases as $p$ increases; a few particular parameter combinations result
into a monotone convergence towards the equilibrium (e.g., for a high value of
$p$).

Figure 2: One parameter bifurcations in $\alpha$, $\zeta$ and $\beta$ of Eq.
(3.1), with other parameters as in Table 2. The immune killing rate $\zeta$
has been set to the base value $0.50\;$h-1. In order to facilitate comparison
with the forthcoming two-dimensional simulations, we set $\alpha=\pi
r^{2}\alpha_{z}$ with $r=5\;$mm (corresponding to a late stage of tumour
growth). The green dots show the maximum and minimum values of $U/K$ during
the oscillations of the stable limit cycle. The solid lines show the value of
the equilibrium of $U$ divided by $K$; the line is red if the equilibrium is
stable and black if it is unstable. Observe that for low values of $\beta$ the
infection-free equilibrium close to carrying capacity is stable.

##  4 Comparison between agent-based and continuum models

In this section we compare numerical simulations for the agent-based model and
the corresponding system of PDEs. It is useful to recall that the equation

$\partial_{t}u=D\partial_{xx}^{2}u+f(u)u$ (4.1)

admits traveling waves with speed at least $2\sqrt{f^{\prime}(0)D}$ as
solutions and an initial condition with compact support evolves into a wave
that travels with the minimal speed [39]. Hence, in absence of infection and
immune response, the uninfected cells form a traveling wave invading the
surrounding area at the speed $2\sqrt{Dp}$. Let us also recall that, in
absence of immune response, infected cells invading a region in which
uninfected cells are at carrying capacity form a traveling wave moving at
speed $2\sqrt{D(\beta-q)}$ (this can be obtained using standard linearisation
techniques [66, §2.1]; we refer to Ref. [51] for a further discussion and
numerical simulations of this situation). These expressions only hold in one
dimension; our simulations are in two dimensions with radial symmetry, hence
the same formulas describe the asymptotic speed (see for example [52, §13.2]).

These results cannot be easily generalised to Eq. (2.5), as the chemotactic
movement of immune cells does not fall within this framework. We may still
expect that uninfected cells invade a region in which the amount of immune
cells can be neglected and thus the expression $2\sqrt{Dp}$ is still a good
approximation of the invasion speed. However, we should not expect the
considerations about the equilibria of Eq. (3.1) to be directly applicable to
the expanding central region, due to both chemotaxis and the nonlocal term in
the equation.

In our simulations we consider a spatial domain $\Omega\coloneqq[-L,L]^{2}$
with $L=10\;$mm and we adopt zero-flux boundary conditions. We define $\omega$
as in Eq. (2.6) in order to maintain the radial symmetry of the problem, with
$R=L$. The initial conditions are

$u_{0}(\boldsymbol{x})=\begin{cases}0.9\ K\quad&\text{for
}\lvert\boldsymbol{x}\rvert\leq R_{u}\\\ 0\ \quad&\text{for
}\lvert\boldsymbol{x}\rvert>R_{u}\end{cases}\qquad
i_{0}(\boldsymbol{x})=\begin{cases}0.1\ K\quad&\text{for
}\lvert\boldsymbol{x}\rvert\leq R_{i}\\\ 0\quad&\text{for
}\lvert\boldsymbol{x}\rvert>R_{i}\end{cases}$ (4.2)

where $R_{u}$ and $R_{i}$ are respectively the initial radius of uninfected
and infected cells; initial conditions for $z$ and $\phi$ are 0 across the
whole domain. The reference case with $R_{u}>R_{i}$ corresponds to the
intratumoral injection of the virus [33]. On the other hand, in some
simulations we assume that $R_{u}=R_{i}$, which corresponds to an infection of
the whole domain: since we consider a tumour that can be infiltrated by the
immune system without major obstacles, it is reasonable to assume that this
could be obtained with an intravenous administration of the virus [33].

We use the parameters listed in Table 2, whose choice is motivated in the
Appendix B. We assume $\zeta=0.50\,$h-1 as the base immune killing rate of
tumour cells; an enhancement of the immune system (e.g., due to immune
checkpoint inhibitor therapy) is then modelled by increasing it to
$\zeta=5.00\,$h-1.

###  4.1 Interactions between the tumour and the immune system in absence of
virotherapy

Figure 3: Numerical simulations of the agent-based model with the parameters
given in Table 2 and $\zeta=0.50\,$h-1 in absence of oncolytic viral infection
(i.e., with $R_{i}=0$). Panel (a) represents cell densities obtained on the
horizontal section of the domain $[-L,L]\times\\{0\\}$ by averaging 10
simulations: the density of the uninfected tumour cells is the blue solid line
and the density of the immune cells is the green solid line; observe that the
scales are different. The dotted black lines show the numerical solution of
Eq. (2.7). The vertical blue dashed lines represent the expected positions of
the uninfected invasion front, traveling at speed $2\sqrt{D_{u}p}$. The
horizontal solid black lines show respectively the equilibrium of the ODE for
the immune cell density in absence of infected cells $S_{z}/q_{z}$ (the actual
concentration is larger inside the tumour due to chemotaxis) and the expected
central tumour density $K-\frac{\zeta}{p}z(t,0)$ with $t=1000\;$h. All the
other panels show the result of a single simulation. The dashed cyan circles
in panels (b) and (c) represent the expected positions of the tumour invasion
front, traveling at speed $2\sqrt{D_{u}p}$. The dotted green circle in panel
(b) represents the internal minimum of the numerical solution of Eq. (2.7). In
panel (d) solid lines refer to the agent-based model and dotted lines refer to
the continuum model. In all the cases the maximum of the axis and the
colorbars correspond to the maximum over time of the quantity plotted. Figure
4: Numerical simulations of the agent-based model with the parameters given in
Table 2 and $\zeta=5.00\,$h-1 in absence of oncolytic viral infection (i.e.,
with $R_{i}=0$). All the graphical elements have the same meaning of Fig. 3.

We first describe the behaviour of the model without oncolytic viral
infection, in order to better understand the basic interactions between the
tumour and the immune system. Fig. 3 shows an excellent qualitative agreement
between numerical solutions of the system of PDEs (2.7) with
$i_{0}(\boldsymbol{x})=0$ and a single simulation of the agent-based model.
The number of immune cells involved is so low that stochastic fluctuations
cannot be neglected, hence the quantitative difference between the two
modelling approaches is significant; however, it is enough to consider the
average over 10 to obtain an improved quantitative agreement. At the beginning
of the simulation there are no immune cells in the domain, hence the tumour
starts to grow towards the carrying capacity and to invade the surrounding
area at the speed $2\sqrt{Dp}$ (vertical blue lines in Fig. 3a). In the
meantime, immune cells enter the domain at the constant rate $S_{0}$. Although
uninfected cells secrete much less chemoattractant than infected cells, the
high number of cells guarantees a chemoattractant secretion sufficient to
guide immune cells. Therefore, the immune cell density stabilizes around the
equilibrium value $S_{0}/q_{z}$ (lower horizontal black line in Fig. 3a) only
far from the tumour, while it is higher inside the tumour and almost 0 around
the boundary due to chemotaxis. The presence of immune cells decreases the
tumour cell density to approximately $K-\frac{\zeta}{p}z(t,\boldsymbol{x})$
(the upper horizontal black line in Fig. 3a uses the value $z(t,0)$ for
simplicity).

As time passes, the area of the tumour increases and the total number of
immune cells reaches steady state; this is due to our assumption that
uninfected tumour cells are unable to stimulate the immune system. As a
consequence, immune cell density decreases and tumour cell density increases.
The maximum density of the chemoattractant stabilizes around a value slightly
larger than $1\,\mu$g/mm2, which is a bit more than a third of $\phi^{*}$.
Overall, the total number of tumour cells increases in time, although at a
lower rate than it would do in absence of immune response.

When the immune killing action is enhanced, clearly the tumour cell density
decreases, as Fig. 4 shows. It is important to observe that in this situation,
even though the total number of immune cells does not change, the immune cell
density inside the tumour is smaller than before: this is due to the fact that
less tumour cells secrete less chemoattractant (the maximum value is now
around $0.85\,\mu$g/mm2); as a consequence, the chemotactic component of the
immune cell movement is weaker than before and the immune density is more
homogeneous in the whole domain due to diffusion. The tumour is still very far
from eradication and even a more effective immune system may not be able to
completely eradicate the mass, as a few tumour cells do not secrete enough
chemoattractant to guide immune cells. This is in line with the empirical
observation that immunotherapy alone usually cannot eradicated cold tumours
[22]. We remark that an analogous situation is observed when the immune inflow
is multiplied by a factor of 10 (not shown), i.e. no eradication is achieved
even when the immune inflow is highly enhanced.

###  4.2 Central infection by oncolytic virus and weak immune response

Figure 5: Numerical simulations of the agent-based model with the parameters
given in Table 2, $\zeta=0.50\,$h-1 and central oncolytic viral infection.
Panel (a) represents cell densities obtained on the horizontal section of the
domain $[-L,L]\times\\{0\\}$ obtained by averaging 10 simulations: the density
of the uninfected tumour cells is the blue solid line and the density of the
infected tumour cells cells is the red solid line; immune cell density is not
shown, as it would superimpose the infected tumour cell density. The dotted
black lines show the numerical solution of Eq. (2.7). The vertical blue dashed
lines represent the expected positions of the uninfected invasion front,
traveling at speed $2\sqrt{D_{u}p}$. The horizontal solid black lines show the
equilibrium of the ODE given by Eq. (3.2) with the value $\alpha=\alpha_{z}\pi
r(t)^{2}$, where $r(t)$ is the radius of infected cells at time $t=1000\;$h.
All the other panels show the result of a single simulation. The dashed cyan
circles in panels (b) and (c) represent the expected positions of the tumour
invasion front, traveling at speed $2\sqrt{D_{u}p}$. The dotted green circle
in panels (b) and (c) represents the internal minimum of the numerical
solution of Eq. (2.7). In panel (d) solid lines refer to the agent-based model
(uninfected, infected and immune cells are represented respectively in blue,
red and green) and dotted lines refer to the continuum model. The maximum of
the axis and the colorbars correspond to the maximum over time of the quantity
plotted, except panel (a), in which the maximum was scaled to enhance
readability.

Expanding on the results of the previous section, we now investigate the
effects of an oncolytic viral infection on the immune response and tumour
growth. Fig. 5, along with the video accompanying it (see Online Resource S2
and S3), shows an excellent quantitative agreement between numerical solutions
of the system of PDEs (2.7) and the average over ten numerical simulations of
the agent-based model (for the sake of clarity, the figure represents the
central section of the domain, i.e. the set $[-L,L]\times\\{0\\}$). The
dynamics in this reference situation are very similar to the ones described in
Ref. [51] in absence of an immune response, therefore we briefly present the
main takeaways. At the beginning of the simulation, the central region of the
tumour is quickly infected; the outer region starts to grow, but does not
reach carrying capacity due to the immune response. The invasion speed of
uninfected cells is $2\sqrt{Dp}$ (vertical blue lines in Fig. 5), as predicted
by theoretical results. On the other hand, the invasion speed of the infected
cells in the uninfected region is slightly less than $2\sqrt{D(\beta-q)}$
(vertical red lines in Fig. 5), which is the speed value observed in absence
of immune response; the reason is that the density of uninfected cells is
below $K$. The highest peak of immune cell density corresponds to the invasion
front of infected cells, as this is the region of the steepest gradient of
chemokines. In the central region the chemokines’ density is constantly at
$\phi^{*}$, hence the chemotactic component of the movement is weaker and the
density of immune cells is almost constant due to the diffusion, although at a
value much higher that $S_{0}/q_{z}$.

The parameters we chose are such that the infection eventually reaches the
front of the wave of uninfected cells around time $t=200\;$h; as a
consequence, the peak at the front starts to decrease for both populations and
infected cells slow down. The final peak of the uninfected cells is
approximately

$\bar{u}\coloneqq\Bigl{(}\frac{q}{\beta}+\frac{Dp}{D\beta}\Bigr{)}K$ (4.3)

which allows the infected front to move at speed $2\sqrt{Dp}$ (we refer to
Ref. [51] for more details). The fact that all the formulas related to the
invasion speeds are not affected by the immune system is not surprising, since
the linearised equations do not change. On the other hand, the central
densities are affected by immune cells: the exact values are hard to predict
due to the presence of the chemotactic term, but we can verify that uninfected
cell density is higher than in absence of immune response and infected cell
density is lower, as predicted by the analysis of Eq. (3.1). As time passes,
both front waves keep moving at the speed $2\sqrt{Dp}$ and central densities
of tumour cells and immune cells have some oscillations while converging
towards an equilibrium. In the central region the chemokines’ density
decreases to approximately $1.4\;\mu$g/mm2 due to the reduced amount of cells:
this is still enough to guide a high number of immune cells towards the
tumour, but the peak of immune cell density at the boundary of the tumour is
lower than before.

###  4.3 Emergence of oscillations

Figure 6: A numerical simulation of the agent-based model with the parameters
given in Table 2, $\zeta=0.50\,$h-1 and wide oncolytic viral infection (i.e.,
$R_{i}=R_{u}$). The dashed cyan circles in panels (a) and (c) represent the
expected positions of the tumour invasion front, traveling at speed
$2\sqrt{D_{u}p}$. The dotted green circles in panels (a), (b) and (c)
represent the internal minimum of the numerical solution of Eq. (2.7). The
dashed red circle in panel (b) represents the front given by the numerical
solution of Eq. (2.7). In panel (d) solid lines refer to the agent-based model
(uninfected, infected and immune cells are represented respectively in blue,
red and green) and dotted lines refer to the continuum model. In all the cases
the maximum of the axis and the colorbars correspond to the maximum over time
of the quantity plotted.

The discussion of Section 3 suggests that some parameter ranges may lead to
persistent oscillations in the centre of the domain. We should also take into
account that these oscillations may be biologically relevant even if they
converge towards a stable equilibrium, as the convergence may be very slow.
Fig. 6 along with the video accompanying it (see Online Resource S4), indeed
shows an example of this situation: the only difference with respect to the
reference case is the initial condition $R_{i}=R_{u}$ (i.e., the initial viral
injection covers the whole domain); we therefore expect the same asymptotic
behaviour of the reference situation, but up to time $t=1500\;$h the
difference is very significant. This can be explained by noting that the
initial number of infected cells is higher, hence more immune cells are
involved and this causes wider oscillations. It is interesting to observe that
the oscillations of the agent-based model are delayed with respect to the ones
of the continuum model (this is easy to see from the total number of cells):
when the infected cell density is very low, stochasticity becomes relevant and
in some regions infected cells go extinct; hence, the following infected cell
regrowth is at first inhomogeneous and it takes some time to diffuse in the
whole domain. The spatial inhomogeneities quickly disappear and the delay of
the oscillations is the main difference between the two modelling approaches.

This example suggests that oscillations may bring the cell density at such low
levels that the agents go extinct, even though the continuum model predicts
recurrence. Infected cells are much more likely to be eradicated then
uninfected cells, since even a small population of the latter tends to regrow;
this means that in practice most of the time after the first oscillation the
tumour keeps growing as it would do in absence of virotherapy. It makes sense
to take the bifurcation diagrams of Section 3 as a starting point and study
whether system (2.7) has the same behaviour as parameters vary. Fig. 7 indeed
shows that, as $\zeta$ increases, the oscillations that appear show an
increasing maximum value and a decreasing minimum value, suggesting that in
the agent-based model both uninfected and infected cells should become close
to eradication. The influence of the parameter $p$ is more interesting: the
stability of the equilibrium $(u^{*},i^{*},z^{*})$ appears independent of its
value, but the size of the oscillations during the convergence decreases as
$p$ increases (see Fig. 7); this leads to the counter-intuitive result that a
fast-growing tumour has a more predictable behaviour under therapy than one
with a slower growth.

Figure 7: Oscillations at the origin from numerical simulations of the PDE for
different values of $\zeta$ and $p$, with the other parameters as in Table 2,
$\zeta=0.50\;$h-1 in panel (a) and the chemotactic coefficient $\chi$ reduced
to $1.65\times 10^{-2}\;$mm2/h (i.e., one tenth of the reference value) in
order to guarantee numerical stability for all the parameter values under
investigation. The simulations are run respectively until $4000\;$h and
$3500\;$h on circular domains of radius $20\;$mm and $25\;$mm in order to
avoid boundary effects. The blue dots show the maximum and minimum values of
$u(\cdot,0)$ during the oscillations after $t=3000\;$h (when present). The red
stars show the value of $u$ in the centre at the last time in the cases in
which oscillation dampen significantly.

###  4.4 Enhanced immune response in presence of viral infection

Figure 8: A numerical simulation of the agent-based model with the parameters
given in Table 2, $\zeta=5.00\,$h-1 and central oncolytic viral infection. The
dashed cyan circles in panels (a) and (c) represent the expected positions of
the tumour invasion front, traveling at speed $2\sqrt{D_{u}p}$. The dotted
green circles in panels (a), (b) and (c) represent the internal minimum of the
numerical solution of Eq. (2.7). The dashed red circle in panel (b) represents
the front given by the numerical solution of Eq. (2.7). In panel (d) solid
lines refer to the agent-based model (uninfected, infected and immune cells
are represented respectively in blue, red and green) and dotted lines refer to
the continuum model. In all the cases the maximum of the axis and the
colorbars correspond to the maximum over time of the quantity plotted. Figure
9: A numerical simulation of the agent-based model with the parameters given
in Table 2, $\zeta=5.00\,$h-1 and wide oncolytic viral infection (i.e.,
$R_{i}=R_{u}$). The dashed cyan circles in panels (a) and (c) represent the
expected positions of the tumour invasion front, traveling at speed
$2\sqrt{D_{u}p}$. The dotted green circles in panels (a), (b) and (c)
represent the internal minimum of the numerical solution of Eq. (2.7). The
dashed red circle in panel (b) represents the front given by the numerical
solution of Eq. (2.7). In panel (d) solid lines refer to the agent-based model
(uninfected, infected and immune cells are represented respectively in blue,
red and green) and dotted lines refer to the continuum model. In all the cases
the maximum of the axis and the colorbars correspond to the maximum over time
of the quantity plotted.

The previous discussion motivates our interest in exploring different
parameter ranges. While the increase of the infection rate $\beta$ may be
challenging to be implemented biologically, the enhancement of the immune
response appears more feasible, for example through the use of immune
checkpoint inhibitors or other immune boosting techniques (T-cell transfer,
immune system modulators, etc.). Fig. 8 along with the video accompanying it
(see Online Resource S5), shows this situation, with the immune killing rate
of cancer cells increased to $5.00\;$h-1. The central region of the tumour is
quickly infected and the total cancer cell density decreases significantly
shortly after due to the action of immune cells. When the central cell
concentration is very low, the immune cells move quickly towards the outer
region of the tumour, as it secretes more chemoattractant. As the number of
infected cells decrease, the inflow of immune cells significantly reduces: as
a consequence, around time $t=200\;$h, while the cell density in the outer
region is still decreasing, the central uninfected cell density starts to
increase again due to the absence of infected and immune cells. The agent-
based model and the continuum model show an excellent quantitative agreement
up to approximately $t=490\;$h, when we observe the recurrence of the
infection only in the continuum model, which also stimulates again the immune
system. In the agent-based model, infected cells are extinct in most of the
domain, hence the second infection is less efficient because it takes more
time for the infected cells to diffuse again in the whole tumour. It is
interesting to observe that the infection remains confined in the centre of
the tumour even at later times due to the strong immune response. Overall, the
dynamics of the two models differ significantly.

The main reason for the failure of the therapy in the previous case is the
fact that the immune response is strong enough to eradicate infected cells,
but too weak to do the same with uninfected cells. A possible solution is to
try and increase the immune cell density: a natural way to do so is having
more infected cells to increase immune inflow; in the previous case the
infection was stopped before reaching the tumour boundary, therefore we now
assume that the initial infection is spread in the whole tumour so that we
avoid such issue. The result is shown in Fig. 9. In this situation, the tumour
is infected everywhere and the immune cells kill all infected cells and most
of the uninfected cells. As time passes, the number of immune cells in the
domain decreases; the few remaining uninfected cells left secrete too little
chemoattractant to guide the immune action, hence they start to regrow.
Despite the low cell numbers involved, stochasticity does not play a key role
in the process: we performed thirty identical simulations with different
random seeds and they all show the same qualitative behaviour and negligible
quantitative differences. Overall, the regrowth is only slightly delayed with
respect to the continuum model. At long times this situation is even worse
than the previous one, since the infection has completely extinguished.

It is interesting to observe that the increase of immune cell density at later
times (due for example to CAR-T therapy [26]) is not enough alone to eradicate
the few tumour cells left, as T-cells are unable to effectively detect tumour
cells (not shown). On the other hand, for higher chemoattractant secretion
rates the immune system may be able to keep the few surviving cancer cells
under control and in some cases even to completely eradicate the tumour. In
this settings we observe a significant effect of stochasticity due to the very
low uninfected cell number involved, hence we performed one hundred
simulations and compared the final cell number for three different values of
the chemoattractant secretion rates, as shown in Fig. 10. When $\alpha_{\phi}$
has the reference value and $\gamma_{\phi}=\alpha_{\phi}$ (Fig. 10a), tumour
eradication is observed only in a few simulations and in some others the
tumour appears under control for a long time (see supplementary material S6
for an example of this); nevertheless, the majority of the simulations show
recurrence. As the values of $\alpha_{\phi}$ and $\gamma_{\phi}$ increase
(Fig. 10b), eradication becomes more likely to happen; in the few cases in
which recurrence is observed, it happens later. When the secretion rates are
even higher (Fig. 10c), the immune system appears always able to either
eradicate the tumour in the wide majority of cases. As the secretion rates
increase, the number of simulations in which the tumour is not completely
eradicated become negligible. These scenarios could be interpreted as the
situation of a tumour with a very high mutational burden that is well
recognised by the immune cells in the area, despite not stimulating the immune
system by itself (see, for example, the review [55] and type 1 tumours
described there); it is therefore highly likely that such a tumour never
reaches a significant size, as any attempt to grow would be immediately
stopped by the immune cells already present in the area. Our results show that
immunovirotherapy could be efficient in the few cases in which this kind of
tumours evade the immune control. However, it appears unlikely that a direct
increase of the chemokine secretion of uninfected tumour cells could be
implemented clinically. We remark that in all these situations the agent-based
model differs significantly from the continuous model, which always shows
tumour relapse.

Figure 10: Violin plots of the final cell number at time $t=1500\;$ obtained
from one hundred simulations of the agent-based model with the parameters
given in Table 2, $\zeta=5.00\,$h-1, wide oncolytic viral infection (i.e.,
$R_{i}=R_{u}$) and different chemoattractant secretion rates. The blue dots
show the single results, not shown when they are at $0$; this happens in 6% of
simulations of panel (a), 72% of simulations in panel (b) and 92% of
simulations in panel (c). The white dots show the median and the blue
horizontal lines show the average. The dark blue areas show the region between
the first and the third quartile; in panel (c) this is not shown, as both
quartiles coincide with 0. Randomness plays a very important role in
therapeutic outcomes.

###  4.5 Different treatment protocols

Figure 11: A numerical simulation of the agent-based model with the parameters
given in Table 2, central oncolytic viral infection, $\zeta=0.50\,$h-1 up to
time $t=200\;$h and $\zeta=5.00\,$h-1 afterwards. The dashed cyan circles in
panels (a) and (c) represent the expected positions of the tumour invasion
front, traveling at speed $2\sqrt{D_{u}p}$. The dotted green circles in panels
(a), (b) and (c) represent the internal minimum of the numerical solution of
Eq. (2.7).The dashed red circle in panel (b) represents the front given by the
numerical solution of Eq. (2.7). In panel (d) solid lines refer to the agent-
based model (uninfected, infected and immune cells are represented
respectively in blue, red and green) and dotted lines refer to the continuum
model. In all the cases the maximum of the axis and the colorbars correspond
to the maximum over time of the quantity plotted. Figure 12: A numerical
simulation of the agent-based model with the parameters given in Table 2,
$\zeta=5.00\,$h-1 and wide oncolytic viral infection (i.e., $R_{i}=R_{u}$); a
second viral injection is performed at time $t=900\;$h, infecting $30$% of
cells everywhere. The dashed cyan circles in panels (a) and (c) represent the
expected positions of the tumour invasion front, traveling at speed
$2\sqrt{D_{u}p}$. The dotted green circles in panels (a), (b) and (c)
represent the internal minimum of the numerical solution of Eq. (2.7). The
dashed red circle in panel (b) represents the front given by the numerical
solution of Eq. (2.7). In panel (d) solid lines refer to the agent-based model
(uninfected, infected and immune cells are represented respectively in blue,
red and green) and dotted lines refer to the continuum model. In all the cases
the maximum of the axis and the colorbars correspond to the maximum over time
of the quantity plotted.

The best outcome that we have achieved so far without changing the
chemoattractant secretion rate is a temporary tumour remission in the
situation of an initial infection that affects the whole area of the tumour.
It is natural to wonder whether the same result can be obtained when the
initial infection is only in the centre. Fig. 11, along with the video
accompanying it (see Online Resource S7), shows that this is possible if the
immune killing rate is increased at time $t=200\;$h, corresponding to the
moment in which the infected front reaches the tumour boundary. This is
because the amount of infected cells at that time is enough to stimulate an
appropriate immune response and its spread configurations allows the immune
system to attack every area of the tumour. The following dynamics are similar
to the ones of Fig. 9, except that the relapse is slightly faster.

Our goal is now to exploit this temporary remission and try to achieve a
better therapeutic outcome. A possible solution would be to have a second
viral injection when the tumour starts to relapse, but this is challenging for
several reasons. First of all, it is not clear when this should be done: on
the one hand, we do not want to wait until the tumour is too big as this is
inconvenient for the patient; on the other hand, if we do it early when the
number of cells is too low, the infection quickly dies away. A second
challenge is the location of the injection: the spatial configuration during
the remission is very sparse and it is impossible to predict where the
infection should start in order to be effective. We ignore the second issue,
assuming that the virus may easily reach any area of the tumour, and focus on
the first one. We assume that the immune system is enhanced all the time and
the first wide injection is followed by a second one at time $t=900\;$h; we
assume that this second viral injection causes $30\%$ of the tumour cells to
become infected, irrespective of their location. Fig. 12, along with the video
accompanying it (see Online Resource S8), shows that we can indeed keep the
tumour under control for a longer period of time, although in the end we
always observe recurrence.

This suggests that a good therapeutical approach could be to perform periodic
repeated viral injections: let us now focus on optimising the schedule. Fig.
12 shows a very good quantitative agreement between numerical solutions of the
system of PDEs (2.7) and single numerical simulations of the agent-based
model. We exploit this fact to reduce our attention to the continuum model and
simulate an automatic viral injection when the cell count reaches a fixed
threshold decided a priori. Fig. 13 shows that, as this threshold increases,
the minimum tumour size achieved decreases, but the total area under the curve
increases. An ideal treatment would require very frequent viral injections,
but its implementation in real life may be inconvenient. Nonetheless, for some
chronic conditions where the patient requires a lifelong, periodic monitoring,
this approach should not be completely discarded.

Figure 13: Numerical solutions of Eq. (2.7) with the parameters given in Table
2, $\zeta=5.00\,$h-1 and wide oncolytic viral infection (i.e., $R_{i}=R_{u}$);
after the first $500\;$h, a new wide viral injection is performed in the whole
tumour as soon as the cell count reaches a given threshold (decided a priori)
and as a result $30$% of cells in every location become infected. Panel (a)
shows the total tumour cell numbers for some thresholds. Panels (b) and (c)
show respectively the total tumour cell number from $t=0\;$h to $t=3000\;$h
and the average time between two consecutive injections for different values
of the threshold at which the injection is performed.

##  5 Conclusions

Figure 14: Comparison between the sum of tumour cells obtained as the average
of five agent-based simulations in different regimes in case of central viral
injection (a) and wide viral injection (b). It is here evident that the
enhancement of the immune system during the virotherapy may worsen the
outcome. The delay of the enhancement in case of central injection is very
similar to the case of a wide injection with the immune system enhanced from
the beginning.

A minimal, hybrid discrete-continuum model for the interactions between tumour
cells, oncolytic viruses and immune system has been developed. The
deterministic continuum counterpart is formally derived and the numerical
results of the two approaches are compared. The main assumption is that the
tumour under investigation is immunologically cold (i.e., its immmunogenicity
is very low) and the viral infection stimulates an immune response.

The continuum model is an excellent approximation of the underlying
microscopic model in several cases. This allows to improve our understanding
of the therapeutical outcome in different settings, relying on some analytical
insights coming from the analysis of the nonspatial model and performing
extensive numerical simulations in a reasonable amount of time. On the other
hand, in some situations we observe significantly different behaviours between
the two models. The main explanations for this is the appearance of
oscillations that bring the cell density at very low levels: the continuum
model may then exhibit a quick regrowth, which does not take place in the same
way in the agent-based model due to the extinction of a cell population in
some locations. This extinction is more likely to happen for infected cells,
since uninfected cells have the ability to regrow; therefore, it appears as a
major obstacle to the effectiveness of immunovirotherapy. Even though
stochasticity plays a key role in the process, it is still possible to predict
when the phenomenon may be observed based on the bifurcations of the
corresponding ODE (which mirrors the oscillations of the spatial continuum
model). We remark the importance of taking into account also the transient
behaviour of the system, since oscillations may dampen on a time scale much
longer than the biologically meaningful one.

Our results show that, according to the continuum model, any immune response
has the tendency to decrease the effectiveness of the virotherapy. This holds
true for the agent-based model whenever oscillations are absent or too weak to
drive the uninfected cell number very low; in the latter situation, the
partial extinction of the sole infected populations may result into a complete
failure of virotherapy. On the other hand, stronger oscillations are sometimes
able to lead all the cancer cells close to extinction in the individual-based
model. This happens when the infection has the possibility to propagate in the
whole tumour before the enhancement of the immune system, hence it can only be
achieved if the time and location of the therapies are correctly calibrated.
Fig. 14 summarises the total number of cells of the agent-based model in
different scenarios, clearly showing that the combination of treatments may
worsen the outcome; this result is in line with several experimental evidences
[18]. Even though our model shows recurrence of the tumour at later times in
most cases, such a low number of cells suggests that it would be possible to
completely eradicate the tumour (e.g., by implementing an additional therapy)
or at least to keep it under control (e.g., by repeated viral injections, as
shown in Fig. 13).

It is fundamental to keep in mind that there are several factors that may
hinder this partial success. First, we assume that immune cells kill
uninfected and infected cells at exactly the same rate, but it would be
reasonable to assume that infected cells are more easily recognised and
killed: this highlights the risk that the immune response could stop the
infection and the need to enhance it only when a sufficient number of infected
cells are present. Furthermore, our model neglects spatial constraints or
immunorefractory aspects of some tumour microenvironments that could affect
the viral diffusion and the immune infiltration inside the tumour, which are
well-known obstacles to the success of both virotherapy [51] and immunotherapy
[1] administered alone. On the other hand, in the present work we only take
into account an increase of the immune killing rate that resembles a generic
immune boosting therapy, such as immune checkpoint inhibition. Several
immunotherapies have shown their success when combined to oncolytic
virotherapy [61], such as, but not limited to, adoptive T-cell transfer [41],
CAR-T [26], CAR-T and bispecific T-cell engager (BiTE) [68]. It could
therefore be interesting to analyse whether the combination of different
immunotherapies could partially overcome the above mentioned obstacles, by
refining and augmenting our model.

Another potentially interesting addition to the model is the role of hypoxia,
whose interaction with oncolytic virotherapy is not straight-forward: while
some viruses are able to specifically target hypoxic cells, some others are
unable to effectively infect hypoxic regions due to the reduced protein
translation of those cells [60]. Hypoxia-driven inflammation constitutes an
additional challenge when considering the interactions with the immune system.

From the mathematical point of view, it could be interesting to perform a
rigorous analysis of the PDE that we have obtained. It is well-known that
chemotactic models may lead to blow-up in finite time [27], hence the well-
posedness for long times may not be completely trivial.

##  Appendix A Formal derivation of continuum models

In this Appendix we describe how to derive the model discussed in the main
text.

###  A.1 Uninfected cells

Uninfected cells can first move, then reproduce or die based on the pressure
value, then become infected and finally be killed by immune cells, as
explained in Section 2. The principle of mass balance gives the equation

$\begin{split}u_{j}^{n+1}&=\Bigl{[}\frac{\theta}{2}u_{j-1}^{n}+\frac{\theta}{2}u_{j+1}^{n}+(1-\theta)u_{j}^{n}\Bigr{]}\Bigl{[}1+\tau
G(\rho_{j}^{n})_{+}-\tau G(\rho_{j}^{n})_{-}\Bigl{]}\\\
&\phantom{=}\times\Bigl{(}1-\tau\frac{\beta}{K}i_{j}^{n}\Bigr{)}\Bigl{(}1-\tau\frac{\zeta}{K}z_{j}^{n}\Bigr{)}\end{split}$

and using the algebraic relation $x_{+}-x_{-}=x$, this simplifies to

$\begin{split}u_{j}^{n+1}&=\Bigl{[}\frac{\theta}{2}u_{j-1}^{n}+\frac{\theta}{2}u_{j+1}^{n}+(1-\theta)u_{j}^{n}\Bigr{]}\Bigl{[}1+\tau
G(\rho_{j}^{n})\Bigl{]}\\\
&\phantom{=}\times\Bigl{(}1-\tau\frac{\beta}{K}i_{j}^{n}\Bigr{)}\Bigl{(}1-\tau\frac{\zeta}{K}z_{j}^{n}\Bigr{)}\end{split}$

Let us define

$\Phi\coloneqq\frac{\theta}{2}u_{j-1}^{n}+\frac{\theta}{2}u_{j+1}^{n}-\theta
u_{j}^{n}=\frac{\theta}{2}\delta^{2}\frac{u_{j-1}^{n}+u_{j+1}^{n}-2u_{j}^{n}}{\delta^{2}}$

so that the previous equation becomes

$\displaystyle u_{j}^{n+1}$ $\displaystyle=(u_{j}^{n}+\Phi)\Bigl{[}1+\tau
G(\rho_{j}^{n})\Bigl{]}\Bigl{(}1-\tau\frac{\beta}{K}i_{j}^{n}\Bigr{)}\Bigl{(}1-\tau\frac{\zeta}{K}z_{j}^{n}\Bigr{)}$
$\displaystyle=u_{j}^{n}+\tau
G(\rho_{j}^{n})u_{j}^{n}-\tau\frac{\beta}{K}u_{j}^{n}i_{j}^{n}-\tau\frac{\zeta}{K}u_{j}^{n}z_{j}^{n}+\Phi$
$\displaystyle-\tau^{2}G(\rho_{j}^{n})\frac{\beta}{K}u_{j}^{n}i_{j}^{n}-\tau^{2}G(\rho_{j}^{n})\frac{\zeta}{K}u_{j}^{n}z_{j}^{n}+\tau^{2}\frac{\beta\zeta}{K^{2}}u_{j}^{n}i_{j}^{n}z_{j}^{n}+\tau^{3}G(\rho_{j}^{n})\frac{\beta\zeta}{K^{2}}u_{j}^{n}i_{j}^{n}z_{j}^{n}$
$\displaystyle\phantom{=}+\tau\Phi\Bigl{[}G(\rho_{j}^{n})-\frac{\beta}{K}i_{j}^{n}-\frac{\zeta}{K}z_{j}^{n}-\tau
G(\rho_{j}^{n})\frac{\beta}{K}i_{j}^{n}-\tau
G(\rho_{j}^{n})\frac{\zeta}{K}z_{j}^{n}+\tau\frac{\beta\zeta}{K^{2}}i_{j}^{n}z_{j}^{n}$
$\displaystyle\phantom{=}+\tau^{2}G(\rho_{j}^{n})\frac{\beta\zeta}{K^{2}}i_{j}^{n}z_{j}^{n}\Bigr{]}$

We now divide both sides of the previous equation by $\tau$ and rearrange the
terms to get

$\frac{u_{j}^{n+1}-u_{j}^{n}}{\tau}=G(\rho_{j}^{n})u_{j}^{n}-\frac{\beta}{K}u_{j}^{n}i_{j}^{n}-\frac{\zeta}{K}u_{j}^{n}z_{j}^{n}+\frac{1}{\tau}\Phi+H_{1}$
(A.1)

and observe that every term of $H_{1}$ is multiplied either by $\tau$ or by
$\Phi$.

Let us now assume that there is a function $u\in
C^{2}([0,+\infty)\times\mathbb{R})$ such that $u_{j}^{n}=u(t_{n},x_{j})=u$
(from now on we omit the arguments of functions computed at $(t_{n},x_{j})$);
thus, we can use Taylor expansions for $u$ in time and space as follows

$\displaystyle
u_{j}^{n+1}=u(t_{n}+\tau,x_{j})=u+\tau\partial_{t}u+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau^{2})$
$\displaystyle u_{j\pm
1}^{n}=u(t_{n},x_{j}\pm\delta)=u\pm\delta\partial_{x}u+\frac{1}{2}\delta^{2}\partial_{xx}^{2}u+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})$

This implies that

$\Phi=\frac{\theta}{2}\delta^{2}\partial_{xx}^{2}u+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})$

and thus
$H_{1}=\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau)+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{2})$.
Eq. (A.1) then becomes

$\partial_{t}u+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau^{2})=\theta\frac{\delta^{2}}{2\tau}\partial_{xx}^{2}u+G(\rho)u-\frac{\beta}{K}ui-\frac{\zeta}{K}uz+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau)+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{2})$

Letting $\tau,\delta\to 0$ in such a way that
$\frac{\delta^{2}}{2\tau}\to\tilde{D}$, we obtain

$\partial_{t}u=\theta\tilde{D}\partial_{xx}^{2}u+G(\rho)u-\frac{\beta}{K}ui-\frac{\zeta}{K}uz$

###  A.2 Infected cells

Infected cells can first move, then die, as explained in Section 2. Also,
uninfected cells may be infected. The principle of mass balance gives the
equation

$\displaystyle i_{j}^{n+1}$
$\displaystyle=\Bigl{[}\frac{\theta}{2}i_{j-1}^{n}+\frac{\theta}{2}i_{j+1}^{n}+(1-\theta)i_{j}^{n}\Bigr{]}(1-\tau
q)\Bigl{(}1-\tau\frac{\zeta}{K}z_{j}^{n}\Bigr{)}$
$\displaystyle\phantom{=}+\tau\frac{\beta}{K}i_{j}^{n}(1+\tau
G(\rho_{j}^{n}))\Bigl{(}1-\tau\frac{\zeta}{K}z_{j}^{n}\Bigr{)}(u_{j}^{n}+\Phi)$

which simplifies to

$\displaystyle i_{j}^{n+1}$ $\displaystyle=(1-\tau
q)\Bigl{(}1-\tau\frac{\zeta}{K}z_{j}^{n}\Bigr{)}(i_{j}^{n}+\Psi)+\tau\frac{\beta}{K}u_{j}^{n}i_{j}^{n}+\tau
H_{2}$ $\displaystyle=i_{j}^{n}-\tau
qi_{j}^{n}-\tau\frac{\zeta}{K}i_{j}^{n}z_{j}^{n}+\Psi+\tau H_{3}$
$\displaystyle\phantom{=}+\tau\frac{\beta}{K}u_{j}^{n}i_{j}^{n}+\tau H_{2}$

where

$\Psi\coloneqq\frac{\theta}{2}i_{j-1}^{n}+\frac{\theta}{2}i_{j+1}^{n}-\theta
i_{j}^{n}=\frac{\theta}{2}\delta^{2}\frac{i_{j-1}^{n}+i_{j+1}^{n}-2i_{j}^{n}}{\delta^{2}}$

and

$\displaystyle H_{2}\coloneqq\tau
G(\rho_{j}^{n})\frac{\beta}{K}u_{j}^{n}i_{j}^{n}-\tau\frac{\beta\zeta}{K^{2}}u_{j}^{n}i_{j}^{n}z_{j}^{n}-\tau^{2}G(\rho_{j}^{n})\frac{\beta\zeta}{K^{2}}u_{j}^{n}i_{j}^{n}z_{j}^{n}$
$\displaystyle+\frac{\beta}{K}i_{j}^{n}(1+\tau
G(\rho_{j}^{n}))\Bigl{(}1-\tau\frac{\zeta}{K}z_{j}^{n}\Bigr{)}\Phi$
$\displaystyle H_{3}\coloneqq-q\Psi-\frac{\zeta}{K}z_{j}^{n}\Psi+\tau
q\frac{\zeta}{K}i_{j}^{n}z_{j}^{n}+\tau q\frac{\zeta}{K}z_{j}^{n}\Psi$

Let us observe that every term of $H_{2}$ and $H_{3}$ is multiplied either by
$\tau$, $\Phi$ or $\Psi$. Dividing both sides by $\tau$ and rearranging the
terms we get

$\frac{i_{j}^{n+1}-i_{j}^{n}}{\tau}=\frac{1}{\tau}\Psi-
qi_{j}^{n}+\frac{\beta}{K}u_{j}^{n}i_{j}^{n}-\frac{\zeta}{K}i_{j}^{n}z_{j}^{n}+H_{2}+H_{3}$
(A.2)

Let us now assume that there is a function $i\in
C^{2}([0,+\infty)\times\mathbb{R})$ such that $i_{j}^{n}=i(t_{n},x_{j})=i$, so
that

$\displaystyle
i_{j}^{n+1}=i(t_{n}+\tau,x_{j})=i+\tau\partial_{t}i+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau^{2})$
$\displaystyle i_{j\pm
1}^{n}=i(t_{n},x_{j}\pm\delta)=i\pm\delta\partial_{x}i+\frac{1}{2}\delta^{2}\partial_{xx}^{2}i+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})$

This implies that

$\Psi=\frac{\theta}{2}\delta^{2}\partial_{xx}^{2}i+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})$

and thus
$H_{2}+H_{3}=\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau)+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{2})$.
Eq. (A.2) then becomes

$\partial_{t}i+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau^{2})=\theta\frac{\delta^{2}}{2\tau}\partial_{xx}^{2}i+\frac{\beta}{K}ui-\frac{\zeta}{K}iz-
qi+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau)+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{2})$

Letting $\tau,\delta\to 0$ in such a way that
$\frac{\delta^{2}}{2\tau}\to\tilde{D}$ we obtain

$\partial_{t}i=\theta\tilde{D}\partial_{xx}^{2}i+\frac{\beta}{K}ui-\frac{\zeta}{K}iz-
qi$

###  A.3 Chemoattractant

The chemoattractant is produced by cancer cells, decays at a constant rate and
diffuses, as explained in Section 2; its dynamics are described by Eq. (2.2),
which can be written as

$\frac{\phi_{j}^{n+1}-\phi_{j}^{n}}{\tau}=D_{\phi}\frac{\phi_{j+1}^{n}+\phi_{j-1}^{n}-2\phi_{j}^{n}}{\delta^{2}}+(\alpha_{\phi}i_{j}^{n}+\gamma_{\phi}u_{j}^{n})(\phi_{j}^{n}-\phi^{*})-q_{\phi}\phi_{j}^{n}$
(A.3)

Let us now assume that there is a function $\phi\in
C^{2}([0,+\infty)\times\mathbb{R})$ such that
$\phi_{j}^{n}=\phi(t_{n},x_{j})=\phi$. Then clearly the discrete diffusion
term converges to the second derivative of $\phi$, hence for $\tau,\delta\to
0$ we obtain

$\partial_{t}\phi=D_{\phi}\partial_{xx}^{2}\phi+(\alpha_{\phi}i+\gamma_{\phi}u)(\phi-\phi^{*})-q_{\phi}\phi$

###  A.4 Immune cells

Immune cells can first move, then die, as explained in Section 2. Also, new
immune cells may enter into the domain. The principle of mass balance gives
the equation

$z_{j}^{n+1}=\Bigl{[}F_{j-1\to j}^{n}\,z_{j-1}^{n}+F_{j+1\to
j}^{n}\,z_{j+1}^{n}+(1-F_{j\to j-1}^{n}-F_{j\to
j+1}^{n})z_{j}^{n}\Bigr{]}(1-\tau q_{z})+\tau S_{j}^{n}$

with

$F_{j\to j\pm
1}^{n}\coloneqq\frac{\theta_{z}}{2}+\underbrace{\nu\frac{(\phi_{j\pm
1}^{n}-\phi_{j}^{n})_{+}}{2\phi^{*}}}_{\eqqcolon\tilde{F}_{j\to j\pm
1}^{n}},\qquad
S_{j}^{n}=\biggl{(}S_{0}+\alpha_{z}\sum_{h}I_{h}^{n}\biggr{)}\mathbbm{1}_{\omega}(x_{j})$

Observe that the source term is no longer multiplied by $\delta$, since we are
considering the cell density. The previous equation can be written as

$\begin{split}z_{j}^{n+1}&=(z_{j}^{n}+\Xi)(1-\tau q_{z})+\tau S_{j}^{n}\\\
&=z_{j}^{n}+\Xi-\tau q_{z}z_{j}^{n}-\tau q_{z}\Xi+\tau S_{j}^{n}\end{split}$

where

$\displaystyle\Xi\coloneqq\Xi_{1}+\Xi_{2}$
$\displaystyle\Xi_{1}\coloneqq\theta_{z}\frac{z_{j-1}^{n}+z_{j+1}^{n}-2z_{j}^{n}}{2}$
$\displaystyle\Xi_{2}\coloneqq-(\tilde{F}_{j\to j-1}^{n}+\tilde{F}_{j\to
j+1}^{n})z_{j}^{n}+\tilde{F}_{j-1\to j}^{n}z_{j-1}^{n}+\tilde{F}_{j+1\to
j}^{n}z_{j+1}^{n}$

Dividing both sides by $\tau$ and rearranging the terms we get

$\frac{z_{j}^{n+1}-z_{j}^{n}}{\tau}=\frac{\Xi}{\tau}-q_{z}z_{j}^{n}+S_{j}^{n}+q_{z}\Xi$
(A.4)

Let us now assume that there is a function $z\in
C^{2}([0,+\infty)\times\mathbb{R})$ such that $z_{j}^{n}=z(t_{n},x_{j})=z$, so
that

$\displaystyle
z_{j}^{n+1}=z(t_{n}+\tau,x_{j})=z+\tau\partial_{t}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau^{2})$
$\displaystyle z_{j\pm
1}^{n}=z(t_{n},x_{j}\pm\delta)=z\pm\delta\partial_{x}z+\frac{1}{2}\delta^{2}\partial_{xx}^{2}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})$

This implies that

$\Xi_{1}=\frac{\theta_{z}}{2}\delta^{2}\partial_{xx}^{2}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})$

Furthermore, the assumptions on $\phi$ imply that

$\tilde{F}_{j\to j\pm
1}^{n}=\frac{\nu}{2\phi^{*}}\Bigl{(}\pm\delta\partial_{x}\phi+\frac{1}{2}\delta^{2}\partial_{xx}^{2}\phi+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})\Bigr{)}_{+}=\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta)$

and we can easily conclude that
$\Xi=\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta)$.
We then use the Taylor expansion of $z$ in the definition of $\Xi_{2}$ to get

$\begin{split}\Xi_{2}&=-(\tilde{F}_{j\to j-1}^{n}+\tilde{F}_{j\to
j+1}^{n})z+\tilde{F}_{j-1\to
j}^{n}\Bigl{(}z-\delta\partial_{x}z+\frac{1}{2}\delta^{2}\partial_{xx}^{2}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})\Bigr{)}\\\
&\phantom{=}+\tilde{F}_{j+1\to
j}^{n}\Bigl{(}z+\delta\partial_{x}z+\frac{1}{2}\delta^{2}\partial_{xx}^{2}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})\Bigr{)}\\\
&=(\tilde{F}_{j-1\to j}^{n}-\tilde{F}_{j\to j-1}^{n}+\tilde{F}_{j+1\to
j}^{n}-\tilde{F}_{j\to j+1}^{n})z+\delta(-\tilde{F}_{j-1\to
j}^{n}+\tilde{F}_{j+1\to j}^{n})\partial_{x}z\\\
&\phantom{=}+\frac{1}{2}\delta^{2}(\tilde{F}_{j-1\to j}^{n}+\tilde{F}_{j+1\to
j}^{n})\partial_{xx}^{2}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})\end{split}$

Now, let us observe that

$\displaystyle\tilde{F}_{j\pm 1\to j}^{n}-\tilde{F}_{j\to j\pm 1}^{n}$
$\displaystyle=\frac{\nu}{2\phi^{*}}[(\phi_{j}^{n}-\phi_{j\pm
1}^{n})_{+}-(\phi_{j\pm 1}^{n}-\phi_{j}^{n})_{+}]$
$\displaystyle=\frac{\nu}{2\phi^{*}}(\phi_{j}^{n}-\phi_{j\pm
1}^{n})=\frac{\nu}{2\phi^{*}}\Bigl{(}\mp\delta\partial_{x}\phi-\frac{1}{2}\delta^{2}\partial_{xx}^{2}\phi+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})\Bigr{)}$

using the relation $x_{+}-(-x)_{+}=x_{+}-x_{-}=x$. We therefore have

$\Xi_{2}=\frac{\nu}{2\phi^{*}}\Bigl{\\{}-\delta^{2}\partial_{xx}^{2}\phi
z+\delta[-(\delta\partial_{x}\phi+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{2}))_{+}+(-\delta\partial_{x}\phi+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{2}))_{+}]\partial_{x}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta^{3})\Bigr{\\}}$

Finally, Eq. (A.4) becomes

$\displaystyle\partial_{t}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\tau^{2})$
$\displaystyle=\theta_{z}\frac{\delta^{2}}{2\tau}\partial_{xx}^{2}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\frac{\delta^{3}}{\tau})+\frac{\nu}{\phi^{*}}\
\frac{\delta^{2}}{2\tau}\Bigl{\\{}-\partial_{xx}^{2}\phi z$
$\displaystyle\phantom{=}+[-(\partial_{x}\phi+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta))_{+}+(-\partial_{x}\phi+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta))_{+}]\partial_{x}z+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta)\Bigr{\\}}-q_{z}z+S_{j}^{n}+\mathchoice{{\mathcal{O}}}{{\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptstyle\mathcal{O}$}}(\delta)$

Let us now observe that

$\sum_{h}I_{h}^{n}=\sum_{h}i_{h}^{n}\,\delta=\sum_{h}i(t_{n},x_{h})\delta\xrightarrow{\delta\to
0}\int_{\Omega}i(t_{n},y)\,\mathrm{d}y$

which means that

$S_{j}^{n}=\biggl{(}S_{0}+\alpha_{z}\sum_{h}I_{h}^{n}\biggr{)}\mathbbm{1}_{\omega}(x_{j})\xrightarrow{\delta\to
0}\biggl{(}S_{0}+\alpha_{z}\int_{\Omega}i(t_{n},y)\,\mathrm{d}y\biggr{)}\mathbbm{1}_{\omega}(x_{j})\eqqcolon
S(t_{n},x_{j})$

Letting $\tau,\delta\to 0$ in such a way that
$\frac{\delta^{2}}{2\tau}\to\tilde{D}$ we arrive at the final result:

$\begin{split}\partial_{t}z&=\theta_{z}\tilde{D}\partial_{xx}^{2}z+\frac{\nu\tilde{D}}{\phi^{*}}\\{-\partial_{xx}^{2}\phi
z+[-(\partial_{x}\phi)_{+}+(-\partial_{x}\phi)_{+}]\partial_{x}z\\}-q_{z}z+S\\\
&=\theta_{z}\tilde{D}\partial_{xx}^{2}z+\frac{\nu\tilde{D}}{\phi^{*}}(-\partial_{xx}^{2}\phi
z-\partial_{x}\phi\partial_{x}z)-q_{z}z+S\\\
&=\theta_{z}\tilde{D}\partial_{xx}^{2}z-\frac{\nu\tilde{D}}{\phi^{*}}\partial_{x}(z\partial_{x}\phi)-q_{z}z+S\end{split}$

##  Appendix B Details of numerical simulations

#### Parameter values

The majority of the parameters of the model has been estimated from the
empirical literature, while a few others are specific of our formulation of
the model and have been set to reasonable values in order to reproduce
plausible dynamics. The parameters $p,D,K$ and $\beta$ assume the values
listed in Table 2, which are the same used in Ref. [51]; for the sake of
brevity, we omit further comments on them and report only the references in
the aforementioned table. On the other hand, the basic death rate of infected
cells $q$ has been decreased to $8.34\times 10^{-3}\;$h-1, which is one fifth
of the value used in Ref. [51]. This is due to the fact that in the current
model $q$ does not take into account the death of infected cells due to immune
killing, which is considered separately.

The diffusion coefficient of the chemoattractant $D_{\phi}$ has been taken
equal to $3.33\cdot 10^{-2}\;$mm2/h following Ref. [50]. The saturation
density of the chemoattractant $\phi^{*}$ and the secretion rates of
chemoattractant by infected cells $\alpha_{\phi}$ have been adapted from the
values reported in Ref. [32], which fit the data of IFN $\gamma$ taken from
Ref. [23]. The value of $\phi^{*}$ has been obtained by rescaling to our two-
dimensional setting the estimate of Ref. [32], yielding a value of
$2.92\;\mu$g/mm2. The secretion rate reported in Ref. [32] refers to a single
CD4+T and these cells are assumed to be stimulated by infected tumour cells;
we have adapted their value to our setting by dividing it by the carrying
capacity $K$, obtaining the value $\alpha_{\phi}=2.50\times
10^{-4}\;$mm2/(h$\cdot$cells). Since we are assuming that immune cells are
much less stimulated by uninfected tumour cells, we have set the secretion
rate of chemoattractant by uninfected cells $\gamma_{\phi}$ to $5.00\times
10^{-6}\;$mm2/(h$\cdot$cells), which still allows to obtain a considerable
reduction of the tumour load when the immune system is stimulated by the
infection. In the main text we also discuss how different values of
$\alpha_{\phi}$ and $\gamma_{\phi}$ affect the outcome. The decay of
chemoattractant $q_{\phi}$ has been taken equal to $8.3\cdot 10^{-2}\;$h-1, as
in Ref. [12].

While it is clear from experimental results that the speed of an immune cell
is around $1.08\;$mm/h [65], the estimate of the diffusion and chemotactic
coefficients $D_{z}$ and $\chi$ from this consideration constitute a
particular challenge. The diffusion coefficient of immune cells has been set
to $D_{z}=4.20\cdot 10^{-3}\;$mm2/h, as in Ref. [1], noting that similar
values are used elsewhere in the literature (such as in Ref. [4]). We remark
that, following the approach of Ref. [28, 54], one could estimate a value in
the same order of magnitude relying on reasonable biological assumptions
(although the precise quantities needed are hard to estimate). We performed
several simulations of our agent-based model to conclude that for
$\chi=1.65\;$mm2/h immune cells move toward a gradient of chemoattractant
similar to the one present in our simulations (i.e., the stationary profile of
the chemokines for the initial condition of our simulations) with an average
speed of approximately $0.6\;$mm/h. We decrease this value to
$\chi=0.165\;$mm/h2 to avoid an excessive concentration of immune cells; the
resulting average speed is around $0.06\;$mm/h, which is plausible considering
that immune cells face many physical obstacles in penetrating the tumour
microenvironment.

The death rate of immune cells $q_{z}$ has been taken equal to $7.5\cdot
10^{-3}\;$h-1, which is the value used in [25] for T4 cells. The base inflow
rate of immune cells $S_{0}$ has been set to $5.00\times 10^{-2}$
cells/(mm${}^{2}\cdot$h) in order to have a density of immune cells inside the
tumour coherent with experimental observations [11, 72]. The additional inflow
rate of immune cells due to the infection $\alpha_{z}$ is one of the main
peculiarities of our model and summarises several biological processes, hence
it is hard to find meaningful estimates in the literature. We have set it to
$3.75\times 10^{-5}\;$(mm${}^{2}\cdot$h)-1, which allows to have an immune
cell density comparable with the aforementioned experimental references. The
same problem arises with the immune killing rate of cancer cells $\zeta$:
since the model is very sensitive to this parameter, in the main text we
compare the differences between setting it to $0.50\;$h-1 and $5.00\;$h-1,
corresponding respectively to weak and enhanced immune responses.

#### Numerical simulations for the discrete models

We used a temporal step $\tau=0.02\;$h and a spatial step of $\delta=0.1\;$mm,
as already done in Ref. [51]. All simulations have been performed in Matlab
2021b.

At every iteration the cell numbers and the chemoattractant density are
updated according to the rules described in Section 2. We first consider
movement, then reproductions and deaths of all cell populations, inflow of
immune cells, infections and finally chemoattractant dynamics. Since we only
need to keep track of the collective fate of cells in the same lattice point,
we used the built-in Matlab functions binornd and mnrnd, which compute random
arrays according to binomial and multinomial distributions. Zero-flux boundary
conditions for cell populations are implemented by not allowing cells at the
boundary to leave the domain. The density of the chemoattractant is updated
through the two-dimensional analogue of Eq. (2.2); Neumann boundary conditions
are then implemented by considering additional grid points outside the domain,
with the same density value as the boundary points of the grid.

The one-dimensional plots in Figs. 3a, 4 and 5a are obtained by averaging ten
simulation. The cell sums of Fig. 14 are obtained by averaging five
simulations (although the cell sum obtained from a single simulation does not
show any significant difference). All the two-dimensional plots show a single
simulation. We remark that in all cases we performed at least five simulations
and did not observe any relevant qualitative difference with respect to the
result shown; the only exception is Fig. 10 and the Online Resource S6, as
explained in the main text.

In order to allow reproducibility, a random seed has been set at the beginning
of each new simulation. In the figures representing a single simulation only
the one with random seed equal to 1 is shown (with the exception of the Online
Resource S6, in which it has been set to 4).

#### Numerical simulations for the continuum models

The system of equations (2.7) has been solved with a finite difference scheme
explicit in time, using the discretisations $\Delta t=10^{-4}$, $\Delta
x=0.01$; such a low space step allows to appropriately describe the high peaks
of immune cells of some simulations. The only exceptions are the simulations
of Fig. 7, which are run for a very long time in a bigger domain: here the
discretisation for space is $\Delta x=0.02$, which guarantees stability at
late times without the need to decrease the time step. We used a forward
upwind scheme for the chemotactic term in the equation for immune cell
density, following Ref. [44]; this is a common strategy to deal with this kind
of equations [1, 7]. In the plots of the supplementary material we truncated
the solutions at a value $\frac{1}{\delta^{2}}$ to be consistent with the
representation of the agent-based model. We also use the same threshold
$\frac{1}{\delta^{2}}$ to identify the wave front of infected cells.

## References

  * [1] L. Almeida, C. Audebert, E. Leschiera, and T. Lorenzi, A hybrid discrete–continuum modelling approach to explore the impact of T-cell infiltration on anti-tumour immune response, Bulletin of Mathematical Biology, 84 (2022).
  * [2] N. Almuallem, D. Trucu, and R. Eftimie, Oncolytic viral therapies and the delicate balance between virus-macrophage-tumour interactions: A mathematical approach, Mathematical Biosciences and Engineering, 18 (2021), p. 764 – 799.
  * [3] R. H. I. Andtbacka, F. Collichio, K. J. Harrington, M. R. Middleton, G. Downey, K. Öhrling, and H. L. Kaufman, Final analyses of optim: a randomized phase III trial of talimogene laherparepvec versus granulocyte-macrophage colony-stimulating factor in unresectable stage III–IV melanoma, Journal for ImmunoTherapy of Cancer, 7 (2019).
  * [4] K. Atsou, F. Anjuère, V. M. Braud, and T. Goudon, A size and space structured model describing interactions of tumor cells with immune cells reveals cancer persistent equilibrium states in tumorigenesis, Journal of Theoretical Biology, 490 (2020), p. 110163.
  * [5] A. A. Baabdulla and T. Hillen, Oscillations in a spatial oncolytic virus model.  preprint, 2023.
  * [6] P. Blanchette and J. G. Teodoro, A renaissance for oncolytic adenoviruses?, Viruses, 15 (2023).
  * [7] F. Bubba, T. Lorenzi, and F. R. Macfarlane, From a discrete model of chemotaxis with volume-filling to a generalized Patlak–Keller–Segel model, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, (2020).
  * [8] N. Champagnat and S. Méléard, Invasion and adaptive evolution for individual-based spatially structured populations, Journal of Mathematical Biology, 55 (2007), pp. 147–188.
  * [9] M. A. J. Chaplain, T. Lorenzi, and F. R. Macfarlane, Bridging the gap between individual-based and continuum models of growing cell populations, Journal of Mathematical Biology, 80 (2020), pp. 343–371.
  * [10] N. Charteris and E. Khain, Modeling chemotaxis of adhesive cells: Stochastic lattice approach and continuum description, New Journal of Physics, 16 (2014).
  * [11] K. Chatzopoulos, V. Kotoula, K. Manoussou, K. Markou, K. Vlachtsis, N. Angouridakis, A. Nikolaou, M. Vassilakopoulou, A. Psyrri, and G. Fountzilas, Tumor infiltrating lymphocytes and CD8+ T cell subsets as prognostic markers in patients with surgically treated laryngeal squamous cell carcinoma, Head and Neck Pathology, 14 (2020), p. 689 – 700.
  * [12] A. K. Cooper and P. S. Kim, A cellular automata and a partial differential equation model of tumor–immune dynamics and chemotaxis, in Mathematical models of tumor-immune system dynamics, Springer, 2014, pp. 21–46.
  * [13] A. Desjardins, M. Gromeier, J. E. Herndon, N. Beaubier, D. P. Bolognesi, A. H. Friedman, H. S. Friedman, F. McSherry, A. M. Muscat, S. Nair, K. B. Peters, D. Randazzo, J. H. Sampson, G. Vlahovic, W. T. Harrison, R. E. McLendon, D. Ashley, and D. D. Bigner, Recurrent glioblastoma treated with recombinant poliovirus, New England Journal of Medicine, 379 (2018), p. 150 – 161.
  * [14] R. Eftimie, J. L. Bramson, and D. J. Earn, Interactions between the immune system and cancer: A brief review of non-spatial mathematical models, Bulletin of Mathematical Biology, 73 (2011), p. 2 – 32.
  * [15] R. Eftimie and G. Eftimie, Tumour-associated macrophages and oncolytic virotherapies: A mathematical investigation into a complex dynamics, Letters in Biomathematics, 5 (2018), p. S6 – S35.
  * [16] R. Eftimie, C. MacNamara, J. Dushoff, J. Bramson, and D. Earn, Bifurcations and chaotic dynamics in a tumour-immune-virus system, Mathematical Modelling of Natural Phenomena, 11 (2016), p. 65 – 85.
  * [17] C. E. Engeland, J. P. Heidbuechel, R. P. Araujo, and A. L. Jenner, Improving immunovirotherapies: The intersection of mathematical modelling and experiments, ImmunoInformatics, 6 (2022), p. 100011.
  * [18] A. C. Filley and M. Dey, Immune system, friend or foe of oncolytic virotherapy?, Frontiers in Oncology, 7 (2017).
  * [19] C. Fountzilas, S. Patel, and D. Mahalingam, Review: Oncolytic virotherapy, updates and future directions, Oncotarget, 8 (2017), pp. 102617–102639.
  * [20] A. Friedman and X. Lai, Combination therapy for cancer with oncolytic virus and checkpoint inhibitor: A mathematical model, PLOS ONE, 13 (2018), pp. 1–21.
  * [21] A. Friedman, J. P. Tian, G. Fulci, E. A. Chiocca, and J. Wang, Glioma virotherapy: Effects of innate immune suppression and increased viral replication capacity, Cancer Research, 66 (2006), pp. 2314–2319.
  * [22] J. Galon and D. Bruni, Approaches to treat immune hot, altered and cold tumours with combination immunotherapies, Nature Reviews Drug Discovery, 18 (2019), p. 197 – 218.
  * [23] Y. Gao, Q. Zhou, Z. Matharu, Y. Liu, T. Kwa, and A. Revzin, A mathematical method for extracting cell secretion rate from affinity biosensors continuously monitoring cell activity, Biomicrofluidics, 8 (2014).
  * [24] D. Hanahan and R. A. Weinberg, Hallmarks of cancer: The next generation, cell, 144 (2011), pp. 646–674.
  * [25] W. Hao, E. D. Crouser, and A. Friedman, Mathematical model of sarcoidosis, Proceedings of the National Academy of Sciences, 111 (2014), pp. 16065–16070.
  * [26] J. He, F. Munir, D. Ragoonanan, W. Zaky, S. J. Khazal, P. Tewari, J. Fueyo, C. Gomez-Manzano, and H. Jiang, Combining CAR T cell therapy and oncolytic virotherapy for pediatric solid tumors: A promising option, Immuno, 3 (2023), p. 37 – 56.
  * [27] T. Hillen and K. Painter, A user’s guide to PDE models for chemotaxis, Journal of Mathematical Biology, 58 (2009), p. 183 – 217.
  * [28] T. Hillen and A. Swan, The diffusion limit of transport equations in biology, Lecture Notes in Mathematics, 2167 (2016), p. 73 – 129.
  * [29] Y. Iwai, M. Ishida, Y. Tanaka, T. Okazaki, T. Honjo, and N. Minato, Involvement of PD-L1 on tumor cells in the escape from host immune system and tumor immunotherapy by PD-L1 blockade, Proceedings of the National Academy of Sciences of the United States of America, 99 (2002), p. 12293 – 12297\.
  * [30] A. Jenner, A. Coster, P. Kim, and F. Frascoli, Treating cancerous cells with viruses: Insights from a minimal model for oncolytic virotherapy, Letters in Biomathematics, 5 (2018), pp. S117–S136.
  * [31] A. L. Jenner, P. S. Kim, and F. Frascoli, Oncolytic virotherapy for tumours following a Gompertz growth law, Journal of Theoretical Biology, 480 (2019), pp. 129–140.
  * [32] A. L. Jenner, M. Smalley, D. Goldman, W. F. Goins, C. S. Cobbs, R. B. Puchalski, E. A. Chiocca, S. Lawler, P. Macklin, A. Goldman, and M. Craig, Agent-based computational modeling of glioblastoma predicts that stromal density is central to oncolytic virus efficacy, iScience, 25 (2022).
  * [33] K.-T. Jin, W.-L. Du, Y.-Y. Liu, H.-R. Lan, J.-X. Si, and X.-Z. Mou, Oncolytic virotherapy in solid tumors: The challenges and achievements, Cancers, 13 (2021), pp. 1–28.
  * [34] S. T. Johnston, M. J. Simpson, and R. E. Baker, Modelling the movement of interacting cell populations: A moment dynamics approach, Journal of Theoretical Biology, 370 (2015), pp. 81–92.
  * [35] L. D. Ke, Y.-X. Shi, S.-A. Im, X. Chen, and W. A. Yung, The relevance of cell proliferation, vascular endothelial growth factor, and basic fibroblast growth factor production to angiogenesis and tumorigenicity in human glioma cell lines, Clinical Cancer Research, 6 (2000), pp. 2562–2572.
  * [36] E. Kelly and S. Russell, History of oncolytic viruses: Genesis to genetic engineering, Molecular Therapy, 15 (2007).
  * [37] J.-H. Kim, Y.-S. Lee, H. Kim, J.-H. Huang, A.-R. Yoon, and C.-O. Yun, Relaxin expression from tumor-targeting adenoviruses and its intratumoral spread, apoptosis induction, and efficacy, Journal of the National Cancer Institute, 98 (2006), pp. 1482–1493.
  * [38] Y. Kim, J. Y. Yoo, T. J. Lee, J. Liu, J. Yu, M. A. Caligiuri, B. Kaur, and A. Friedman, Complex role of NK cells in regulation of oncolytic virus–bortezomib therapy, Proceedings of the National Academy of Sciences of the United States of America, 115 (2018), p. 4927 – 4932.
  * [39] A. N. Kolmogorov, I. G. Petrovsky, Petrovsky, and N. S. Piskunov, Étude de l’équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique, Bull. Univ. Moskow, Ser. Internat., Sec. A, 1 (1937), pp. 1–25.
  * [40] N. L. Komarova and D. Wodarz, ODE models for oncolytic virus dynamics, Journal of theoretical biology, 263 (2010), pp. 530–543.
  * [41] T. Krabbe, J. Marek, T. Groll, K. Steiger, R. M. Schmid, A. M. Krackhardt, and J. Altomonte, Adoptive T cell therapy is complemented by oncolytic virotherapy with fusogenic VSV-NDV in combination treatment of murine melanoma, Cancers, 13 (2021), p. 1 – 18.
  * [42] S. E. Lawler, M.-C. Speranza, C.-F. Cho, and E. A. Chiocca, Oncolytic viruses in cancer treatment: A review, JAMA Oncology, 3 (2017), pp. 841–849.
  * [43] T. Lee, A. L. Jenner, P. S. Kim, and J. Lee, Application of control theory in a delayed-infection and immune-evading oncolytic virotherapy, Mathematical Biosciences and Engineering, 17 (2020), p. 2361 – 2383.
  * [44] R. J. LeVeque, Finite difference methods for ordinary and partial differential equations: Steady-state and time-dependent problems, SIAM, Philadelphia, 2007.
  * [45] H. Lodish, A. Berk, C. A. Kaiser, C. Kaiser, M. Krieger, M. P. Scott, A. Bretscher, H. Ploegh, P. Matsudaira, et al., Molecular cell biology, Macmillan, New York, 2008.
  * [46] T. Lorenzi, P. J. Murray, and M. Ptashnyk, From individual-based mechanical models of multicellular systems to free-boundary problems, Interfaces and Free Boundaries, 22 (2020), pp. 205–244.
  * [47] F. R. Macfarlane, X. Ruan, and T. Lorenzi, Individual-based and continuum models of phenotypically heterogeneous growing cell populations, AIMS Bioengineering, 9 (2022), pp. 68–92.
  * [48] K. J. Mahasa, R. Ouifki, A. Eladdadi, and L. de Pillis, A combination therapy of oncolytic viruses and chimeric antigen receptor T cells: A mathematical model proof-of-concept, Mathematical Biosciences and Engineering, 19 (2022), p. 4429 – 4457.
  * [49] N. T. Martin and J. C. Bell, Oncolytic virus combination therapy: Killing one bird with two stones, Molecular Therapy, 26 (2018), p. 1414 – 1422\.
  * [50] A. Matzavinos, M. A. Chaplain, and V. A. Kuznetsov, Mathematical modelling of the spatio-temporal response of cytotoxic T-lymphocytes to a solid tumour, Mathematical Medicine and Biology, 21 (2004), pp. 1–34.
  * [51] D. Morselli, M. E. Delitala, and F. Frascoli, Agent-based and continuum models for spatial dynamics of infection by oncolytic viruses, Bulletin of Mathematical Biology, 85 (2023).
  * [52] J. D. Murray, Mathematical biology: I. An introduction, Springer, New York, 2002.
  * [53] A. S. Novozhilov, F. S. Berezovskaya, E. V. Koonin, and G. P. Karev, Mathematical modeling of tumor therapy with oncolytic viruses: Regimes with complete tumor eination within the framework of deterministic models, Biology Direct, 1 (2006), pp. 1–18.
  * [54] H. G. Othmer and T. Hillen, The diffusion limit of transport equations II: Chemotaxis equations, SIAM Journal on Applied Mathematics, 62 (2002), p. 1222 – 1250.
  * [55] J. S. O’Donnell, M. W. L. Teng, and M. J. Smyth, Cancer immunoediting and resistance to t cell-based immunotherapy, Nature Reviews Clinical Oncology, 16 (2018), p. 151–167.
  * [56] K. J. Painter, Mathematical models for chemotaxis and their applications in self-organisation phenomena, Journal of Theoretical Biology, 481 (2019), p. 162 – 182.
  * [57] C. J. Penington, B. D. Hughes, and K. A. Landman, Building macroscale models from microscale probabilistic models: A general probabilistic approach for nonlinear diffusion and multispecies phenomena, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 84 (2011).
  * [58] P. Pooladvand, C.-O. Yun, A.-R. Yoon, P. S. Kim, and F. Frascoli, The role of viral infectivity in oncolytic virotherapy outcomes: A mathematical study, Mathematical Biosciences, 334 (2021), p. 108520.
  * [59] L. Russell and K.-W. Peng, The emerging role of oncolytic virus therapy against cancer, Chinese Clinical Oncology, 7 (2018).
  * [60] Z. Sheng Guo, The impact of hypoxia on oncolytic virotherapy, Virus Adaptation and Treatment, 3 (2011), p. 71 – 82.
  * [61] T. Shi, X. Song, Y. Wang, F. Liu, and J. Wei, Combining oncolytic viruses with cancer immunotherapy: Establishing a new generation of cancer treatment, Frontiers in Immunology, 11 (2020).
  * [62] K. M. Storey and T. L. Jackson, An agent-based model of combination oncolytic viral therapy and anti-PD-1 immunotherapy reveals the importance of spatial location when treating glioblastoma, Cancers, 13 (2021).
  * [63] K. M. Storey, S. E. Lawler, and T. L. Jackson, Modeling oncolytic viral therapy, immune checkpoint inhibition, and the complex dynamics of innate and adaptive immunity in glioblastoma treatment, Frontiers in Physiology, 11 (2020).
  * [64] A. Surendran, A. L. Jenner, E. Karimi, B. Fiset, D. F. Quail, L. A. Walsh, and M. Craig, Agent-based modelling reveals the role of the tumor microenvironment on the short-term success of combination temozolomide/immune checkpoint blockade to treat glioblastoma, Journal of Pharmacology and Experimental Therapeutics, 387 (2023), p. 66 – 77.
  * [65] J. Textor, A. Peixoto, S. E. Henrickson, M. Sinn, U. H. Von Andrian, and J. Westermann, Defining the quantitative limits of intravital two-photon lymphocyte tracking, Proceedings of the National Academy of Sciences of the United States of America, 108 (2011), p. 12401 – 12406.
  * [66] W. Van Saarloos, Front propagation into unstable states, Physics Reports, 386 (2003), pp. 29–222.
  * [67] G. Vithanage, H.-C. Wei, and S. R.-J. Jang, The role of tumor activation and inhibition with saturation effects in a mathematical model of tumor and immune system interactions undergoing oncolytic viral therapy, Mathematical Methods in the Applied Sciences, 46 (2023), p. 10787 – 10813.
  * [68] A. Wing, C. A. Fajardo, A. D. Posey, C. Shaw, T. Da, R. M. Young, R. Alemany, C. H. June, and S. Guedan, Improving CART-cell therapy of solid tumors with oncolytic virus–driven production of a bispecific T-cell engager, Cancer Immunology Research, 6 (2018), p. 605 – 616.
  * [69] D. Wodarz, Viruses as antitumor weapons: Defining conditions for tumor remission, Cancer Research, 61 (2001), pp. 3501–3507.
  * [70] J. Wojton and B. Kaur, Impact of tumor microenvironment on oncolytic viral therapy, Cytokine and Growth Factor Reviews, 21 (2010), pp. 127–134.
  * [71] J. T. Wu, D. H. Kirn, and L. M. Wein, Analysis of a three-way race between tumor growth, a replication-competent virus and an immune response, Bulletin of Mathematical Biology, 66 (2004), pp. 605–625.
  * [72] K. Yasuda, T. Nirei, E. Sunami, H. Nagawa, and J. Kitayama, Density of CD4(+) and CD8(+) T lymphocytes in biopsy samples can be a predictor of pathological response to chemoradiotherapy (CRT) for rectal cancer, Radiation Oncology, 6 (2011).