# Mathematical modeling of glioma invasion and therapy approaches

Martina Conte Corresponding author<EMAIL_ADDRESS>Department of
Mathematical Sciences ”G. L. Lagrange”, Politecnico di Torino
Corso Duca degli Abruzzi 24 - 10129 Torino, Italy Yvonne Dzierma  Saarland
University Medical Centre, Department of Radiotherapy and Radiation Oncology
Kirrberger Str. Geb. 6.5 - 66421 Homburg/Saar, Germany Sven Knobe  Saarland
University Medical Centre, Department of Radiotherapy and Radiation Oncology
Kirrberger Str. Geb. 6.5 - 66421 Homburg/Saar, Germany Christina Surulescu
Technische Universität Kaiserslautern, Felix-Klein-Zentrum für Mathematik
Paul-Ehrlich-Str. 31 - 67663 Kaiserslautern, Germany

###### Abstract

A multiscale model for glioma spread in brain tissue under the influence of
vascularization and vascular endothelial growth factors is proposed. It
accounts for the interplay between the different components of the neoplasm
and the healthy tissue and it investigates and compares various therapy
approaches. Precisely, these involve radio- and chemotherapy in a concurrent
or adjuvant manner together with anti-angiogenic therapy affecting the
vascular component of the system. We assess tumor growth and spread on the
basis of DTI data, which allows us to reconstruct a realistic brain geometry
and tissue structure, and we apply our model to real glioma patient data. In
this latter case, a space-dependent radiotherapy description is considered
using data about the corresponding isodose curves.

##  1 Introduction

Glioma is the most prevalent, aggressive, and invasive subtype of primary
brain cancer, characterized by rapid cell proliferation and great infiltration
capacity. Its infiltrative spread is strongly related to very poor patient
survival prognosis, considering the difficulty to properly assess tumor
margins for effective treatment [5, 45]. The tumor evolution takes advantage
of the patient-specific anisotropic brain structure, in particular white
matter tracts and blood vessels [22, 23, 38]. In particular, the anisotropic
fiber structure characterizing the brain tissue influences the direction of
cell migration determining preferential paths along which the cells move using
specific receptor-binding mechanisms. Thus, the relevance of these specific
structures in the overall tumor evolution motivates the need of developing
personalized treatment planning, which should include such characteristics.

The first approach to treat newly diagnosed gliomas is usually a maximal
surgical resection of the neoplasm, whose benefits, for both high- and low-
grade gliomas, have been extensively debated since it could cause neurologic
deficits and, often, a complete resection is impossible [47]. In fact, the
tumor outer rims, composed of cells infiltrating the surrounding normal brain
tissue and protected by normal vasculature, are often not visible with the
common medical imaging techniques. Therefore, several supporting treatments
are typically used in combination with surgical resection. The Stupp protocol
has become standard of care for the treatment of gliomas. It consists of
radiotherapy and concomitant chemotherapy with temozolomide. Precisely,
radiotherapy (RT) is provided with the standard dose of 60 Gy delivered in 30
daily fractions of 2 Gy (Monday to Friday) over 6 weeks. RT planning is based
on the definition of specific tumor volumes allowing to assess the margins of
the gross tumor volume and the area to treat. These could be identified using
advanced imaging techniques. The administration of temozolomide changes from 7
days per week at 75 mg per square meter of body-surface area per day, during
radiotherapy, to six cycles of 150-200 mg per square meter for 5 days during
each 28-day cycle, after the radiotherapy. This second part of the
chemotherapy treatment is also named adjuvant chemotherapy. Together with
these standard approaches, new treatments are emerging to target molecules
involved in various tumor mechanisms. Among them, angiogenesis inhibitors have
been tested as anti-migratory agents for glioblastoma (GB). For instance,
bevacizumab, a neutralizing monoclonal antibody against vascular endothelial
growth factors (VEGFs), demonstrated a promising response in patients with
recurrent malignant gliomas in combination with other drugs. However, the role
and efficacy of bevacizumab are still debated, especially when combined with
the classical treatment approaches [42].

From the mathematical viewpoint, several classes of models have been proposed
during the last two decades with the aim of using modern biomedical
visualization techniques to help in predicting the tumor volumes (CTV=clinical
target volume, PTV=planning target volume) for therapy planning. Some of those
models [27, 51] give purely macroscopic descriptions of glioma evolution in
interaction with the tumor environment, others, more recent, have a multiscale
character and obtain effective equations on the macroscale of the whole tumor
from settings describing microscopic and/or mesoscale dynamics [6, 7, 8, 9,
12, 13, 14, 15, 26, 31, 39, 50, 55]. Among these models, [7, 8, 51] suggest
ways to assess tumor growth and grading. Settings taking into account the
interplay of glioma with vascularization and other components (e.g., acidity)
of the peritumoral space have been proposed within both model classes [7, 12,
34, 51]. Models that pay attention to various therapy approaches, in a more or
less detailed way, are, instead, relatively scarce in this context, see e.g.
[25, 26, 43]. Therefore, there is still a need for mathematical formulations
aiming at describing and comparing personalized treatment schedules, which
involve state-of-the-art methods, with the purpose of identifying the most
effective ones. This paper aligns with this aim.

The basis for the mesoscopic description of glioma and endothelial cells (ECs)
dynamics can be found in [6, 7]. With respect to these works, here we extend
those settings considering a more detailed description at the microscopic
level, where subcellular dynamics of tumor receptor binding to anisotropic
brain tissue and blood vessels are taken into account together with the
dynamics of VEGF receptors located on the EC membrane. With respect to [7], we
consider here a different coupling between endothelial and tumor cell dynamics
at the mesoscopic level and we mainly focus on the description and inclusion
of different terms characterizing the effects of radio-, chemo-, and anti-
angiogenic therapy. The upscaling of the corresponding kinetic transport PDEs
for glioma and endothelial cells leads in the parabolic limit to a system of
reaction-advection-diffusion equations, which we couple with the macroscopic
descriptions of VEGFs, healthy, and necrotic tissue evolution.

The subsequent content is organized as follows. In Section 2, we describe the
model setup on the subcellular and mesoscopic scales and we deduce the
corresponding macroscopic PDEs. In Section 3, we concretize the coefficient
functions and parameters, while the numerical simulations are presented in
Section 4. In particular, the last part of this section is dedicated to the
application of the proposed approach to real glioma patient data. A brief
discussion is proposed in Section 5.

##  2 Modeling

In this note, we propose a multiscale model for glioma development and spread
which interconnects the dynamics of glioma cells, vasculature, and vascular
endothelial growth factors. Glioma cells mainly spread according to the
anisotropy of brain tissue or moving along the blood vessels therein [17, 38].
In order to provide nutrient and oxygen supply to its growing cell population,
the tumor produces and releases in the extracellular environment growth
factors, such as VEGFs, that attract endothelial cells. The resulting
angiogenesis boosts the growth of the neoplasm, leading to an abundantly
vascularized and progressed tumor mass. Standard treatment usually comprises
chemo- and radiotherapy, in a concurrent or an adjuvant manner. Thereby, both
tumor and normal tissue are affected, although to a different degree. A
necrotic region is formed within the tumor mass, due to impaired tumor-
associated vascularization. We also consider here the effect of an anti-
angiogenic drug (e.g., bevacizumab [28]), aimed at reducing the affinity of
VEGF receptors on ECs toward their ligand released in the extracellular space
by the tumor cells. This therapeutic approach affects the EC tactic motility
by reducing their bias towards (increasing gradients of) VEGFs and impairing
their proliferation, which, in turn, leads to degradation of glioma through
diminished nutrient supply. Starting from the subcellular level of interaction
between cells, VEGFs, and tissue, we set up the corresponding kinetic
transport equations (KTEs) for glioma cells and endothelial cells, and perform
a (non-rigorous) parabolic limit to deduce the macroscopic system of reaction-
advection-diffusion PDEs for the involved cell populations: glioma and
endothelial cells. These equations are then coupled to the evolution of VEGF
concentration, healthy tissue, and necrotic matter, which are considered
directly on the macroscopic scale.

###  2.1 Subcellular level

On the microscopic scale, we describe the interaction of glioma cells with the
surrounding tissue including extracellular matrix (ECM), neuron bundles, and
other components, thereby also considering the blood vessels (described via EC
density), and the interaction of ECs with VEGFs diffusing in the extracellular
space. The exchange of information between cells and their environment happens
through receptors located on the cell membranes.

Considering such interactions follows the idea employed in [30, 37] to build a
micro-meso model for tumor invasion with chemo- and haptotaxis and revisited
in [6, 8, 9, 12, 13, 14, 15, 26, 33, 34], where upscalings from the KTEs on
mesoscale have been performed, leading to PDEs with cross diffusion and
various kinds of taxis.

We denote by $y_{1}(t)$ the amount of tumor receptors bound to the blood
vessels, by $y_{2}(t)$ the amount of tumor receptors bound to the tissue, and
by $\zeta(t)$ the amount of EC receptors bound to VEGF molecules. The
corresponding binding dynamics is characterized by simple mass action
kinetics:

$\begin{split}(\bar{R}_{y}-y_{1}-y_{2})\,\,+\,\,W\stackrel{{\scriptstyle[}}{{k}}_{W}^{-}]{k_{W}^{+}}{\rightleftharpoons}y_{1}\\\\[5.69046pt]
(\bar{R}_{y}-y_{1}-y_{2})\,\,+\,\,Q\stackrel{{\scriptstyle[}}{{k}}_{Q}^{-}]{k_{Q}^{+}}{\rightleftharpoons}y_{2}\\\\[5.69046pt]
(\bar{R}_{\zeta}-\zeta)\,\,+\,\,H\stackrel{{\scriptstyle[}}{{k}}_{H}^{-}]{k_{H}^{+}}{\rightleftharpoons}\zeta\,.\end{split}$

Here, $\bar{R}_{y}$ is the total amount of receptors on a tumor cell, and
$\bar{R}_{\zeta}$ is that on an endothelial cell; both are assumed to be
constant. We denote by $Q(t,x)$ the macroscopic brain tissue density and by
$W(t,x)$ the macroscopic density of ECs. They depend on the position
$x\in\mathbb{R}^{n}$ and on the time $t>0$. Further, $H(t,x)$ denotes the
concentration of VEGFs. Accordingly, we get the ODE system

$\begin{split}&\dot{y}_{1}=k_{W}^{+}\dfrac{W}{W_{c,0}}S_{W}(d_{r},\alpha_{W})\,(\bar{R}_{y}-y_{1}-y_{2})-k_{W}^{-}y_{1}\\\\[2.84544pt]
&\dot{y}_{2}=k_{Q}^{+}\dfrac{Q}{Q^{*}}S_{Q}(d_{r},\alpha_{Q})(\bar{R}_{y}-y_{1}-y_{2})-k_{Q}^{-}y_{2}\\\\[2.84544pt]
&\dot{\zeta}=k_{H}^{+}(d_{b})\dfrac{H}{H_{c,0}}(\bar{R}_{\zeta}-\zeta)-k_{H}^{-}(d_{b})\zeta,\end{split}$

where $k_{W}^{+}$ and $k_{W}^{-}$, $k_{Q}^{+}$ and $k_{Q}^{-}$ represent
attachment and detachment rates of glioma cells to vessels and tissue,
respectively, while $k_{H}^{+}$ and $k_{H}^{-}$ are the corresponding rates in
the process of EC binding to VEGFs. The parameters $W_{c,0}$ and $Q^{*}$ are
reference density values for endothelial cell and healthy tissue,
respectively. At this subcellular level we also include the effect of therapy.
In particular, we model the effect of radiotherapy on the normal tissue and on
the ECs, with which glioma cells interact by way of their receptors. The
function $S(d_{r},\alpha):=exp(-\alpha d_{r}-\beta d_{r}^{2})$ is used to
model the fraction of cells that survive to a dose $d_{r}$ of radiotherapy:
this is the so-called linear quadratic (LQ) formulation. Thereby, each cell
population affected by radiotherapy has different values for the parameters
$\alpha$ and $\beta$, referring to the lethal lesions produced by a single
radiation track or by two radiation tracks [18]. The parameters $\alpha_{W}$
and $\alpha_{Q}$ refer to the effect of radiotherapy on ECs and on brain
tissue, respectively. In the following we will indicate with $S_{W}$ and
$S_{Q}$ the radiotherapy survival fractions related to $W$ and $Q$, without
explicitly mentioning their dependency on the parameters $\alpha$, $\beta$,
and $d_{r}$. On the other hand, the interaction between ECs and VEGF is
affected by the anti-angiogenic drug at dose $d_{b}$ that reduces the binding
affinity of VEGF molecules to the corresponding receptors on ECs. In
particular, $k_{H}^{+}:=\bar{k_{H}^{+}}l^{+}(d_{b})$ and
$k_{H}^{-}:=\bar{k_{H}^{-}}l^{-}(d_{b})$, with $l^{+}(d_{b})$ and
$l^{-}(d_{b})$ being a decreasing and an increasing function of $d_{b}$,
respectively.

As in [16], we define $y:=y_{1}+y_{2}$ to be the total amount of transmembrane
entities occupied by tissue or ECs, which allows us to lump together the two
corresponding ODEs, into

$\dot{y}=\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}W+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}Q\right)\bar{R}_{y}-y\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}W+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}Q\right)-k_{W}^{-}y_{1}-k_{Q}^{-}y_{2}\,.$

Assuming that $k_{W}^{-}=k_{Q}^{-}=k^{-}_{1}$, we get the microscopic equation
for the subcellular dynamics of glioma cells

$\dot{y}=\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}W+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}Q\right)(\bar{R}_{y}-y)-k^{-}_{1}y\,.$
(1)

Rescaling $y/{\bar{R}}_{y}\leadsto y$ and
$\zeta/{\bar{R}}_{\zeta}\leadsto\zeta$ will further simplify the notation.
Then, the unique steady states of the above equations are given by:

$\begin{split}&y^{*}=\dfrac{\left(k_{W}^{+}\frac{S_{W}}{W_{c,0}}W+k_{Q}^{+}\frac{S_{Q}}{Q^{*}}Q\right)}{k_{W}^{+}\frac{S_{W}}{W_{c,0}}W+k_{Q}^{+}Q\frac{S_{Q}}{Q^{*}}+k^{-}_{1}}=:\bar{f}(W,Q)\\\\[2.84544pt]
&\zeta^{*}=\dfrac{k_{H}^{+}\frac{H}{H_{c,0}}}{k_{H}^{+}(d_{b})\frac{H}{H_{c,0}}+k_{H}^{-}(d_{b})}=:\bar{e}(H)\,.\end{split}$

The variable $y$ characterizes the ’internal’ state of tumor cells (called
’activity’ in the kinetic theory of active particles, see e.g. [4]). Likewise,
$\zeta$ is the activity variable for ECs. In the sequel we will consider the
mesoscopic densities $p(t,x,v,y)$ and $w(t,x,\vartheta,\zeta)$ of tumor cells
and ECs, respectively, hence both depending on such activity variables.
Thereby, $v\in V\subset\mathbb{R}^{n}$ and
$\vartheta\in\Theta\subset\mathbb{R}^{n}$ are the velocity vectors of glioma
and ECs, respectively, with the spaces $V$ and $\Theta$ to be closer explained
in Subsection 2.2 below.

For the glioma cells, we assume that they follow the gradients of tissue and
vasculature, therefore we look at the path of a single cell starting at
position $x_{0}$ and moving to position $x$ with velocity $v$ in the density
fields $Q$ and $W$. Denoting by $z:=y^{*}-y$ the deviation of $y$ from its
steady state, we have:

$\begin{split}\dot{z}&=\dfrac{k_{1}^{-}}{\left(k_{W}^{+}\frac{S_{W}}{W_{c,0}}W+k_{Q}^{+}\frac{S_{Q}}{Q^{*}}Q+k_{1}^{-}\right)^{2}}\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}\,v\cdot\nabla_{x}W+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}\,v\cdot\nabla_{x}Q+\bar{F}_{W}(t)+\bar{F}_{Q}(t)\right)\\\\[8.5359pt]
&-z\left(k_{W}^{+}\dfrac{W}{W_{c,0}}S_{W}+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}Q+k^{-}_{1}\right):=G(z,W,Q)\end{split}$
(2)

with

$\begin{split}&\bar{F}_{W}(t):=k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}\partial_{t}W,\\\\[5.69046pt]
&\bar{F}_{Q}(t):=k_{Q}^{+}\dfrac{S_{Q}}{Q_{c,0}}\partial_{t}Q\,.\end{split}$

For simplicity we denote
$B_{p}(W,Q):=\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}W+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}Q+k_{1}^{-}\right)$.
With an analogous argument, denoting by $u=\zeta^{*}-\zeta$ the deviation of
$\zeta$ from its steady state, we assume that the ECs follow the gradient of
VEGF. Therefore, we look at the path of a single cell starting at position
$x_{0}$ and moving to position $x$ with velocity $\vartheta$ in the
concentration field $H$ so that $H(t,x)=H(t,x_{0}+\vartheta t)$.

Thus, the equation for $u$ is given by:

$\dot{u}=\dfrac{\frac{k_{H}^{+}}{H_{c,0}}k_{H}^{-}}{\left(\frac{k_{H}^{+}}{H_{c,0}}H+k_{H}^{-}\right)^{2}}\left(\partial_{t}\,H+\vartheta\cdot\nabla_{x}H\right)-u\left(\dfrac{k_{H}^{+}}{H_{c,0}}H+k_{H}^{-}\right)=:\Gamma(u,H)$
(3)

For simplicity we denote
$B_{w}(H):=\left(\dfrac{k_{H}^{+}}{H_{c,0}}H+k_{H}^{-}\right)$.

###  2.2 Mesoscopic level

We model the mesoscopic behavior of glioma and endothelial cells with the aid
of KTEs describing velocity-jump processes and taking into account the
subcellular dynamics. Concretely, we consider the following cell density
functions:

  * •

$p(t,x,v,y)$ for glioma cells;

  * •

$w(t,x,\vartheta,\zeta)$ for endothelial cells (ECs) forming capillaries,

with the time and space variables $t>0$ and $x\in\mathbb{R}^{n}$, velocities
$v\in V=s\mathbb{S}^{n-1}$ and $\vartheta\in\Theta:=\sigma\mathbb{S}^{n-1}$,
and activity variables $y\in\mathcal{Y}=(0,1)$ and
$\zeta\in\mathcal{Z}=(0,1)$. These choices mean that we assume for glioma and
ECs constant speeds $s>0$ and $\sigma>0$, respectively, and we denote by
$\mathbb{S}^{n-1}$ the unit sphere in $\mathbb{R}^{n}$. As in [13, 14, 15,
26], in the sequel we will work with the deviations
$z=y^{*}-y\in\mathit{Y}\subseteq(y^{*}-1,y^{*})$ rather then with $y$ for
glioma cells, and with
$u=\zeta^{*}-\zeta\in\mathit{Z}\subseteq(\zeta^{*}-1,\zeta^{*})$ rather then
with $\zeta$ for ECs. The corresponding macroscopic cell densities are denoted
by $M(t,x)$ and $W(t,x)$, respectively.

The kinetic transport equation for the glioma cell density function is given
by:

$\partial_{t}p+\nabla_{x}\cdot\,(vp)+\partial_{z}(G(z,W,Q,S_{W},S_{Q})p)=\mathcal{L}_{p}[\lambda(z)]p+\mathcal{P}\,p-(R_{M}(d_{r})+C_{M}(k_{c}))p$
(4)

where $\mathcal{L}_{p}[\lambda(z)]p$ denotes the turning operator, describing
tumor velocity changes. In particular, such changes are due to contact
guidance: glioma cells tend to align according to the brain tissue anisotropy,
mainly associated with white matter tracts. As in previous models of this
type, $\mathcal{L}_{p}[\lambda(z)]p$ is a Boltzmann-like integral operator of
the form

$\mathcal{L}_{p}[\lambda(z)]p=-\lambda(z)p+\lambda(z)\int_{V}K(x,v)\,p(v^{\prime})\,dv^{\prime}\,,$
(5)

where $\lambda(z):=\lambda_{0}-\lambda_{1}z\geq 0$ is the cell turning rate
depending on the microscopic variable $z$, while $\lambda_{0}$ and
$\lambda_{1}$ are positive constants. The integral term describes the
reorientation of cells from any previous velocity $v^{\prime}$ to a new
velocity $v$ upon interacting with the tissue. The turning kernel $K(x,v)$ is
assumed here to be independent on the incoming velocity $v^{\prime}$ and it
describes the dominating directional cue as given by the orientation of tissue
fibers. Following [13, 39] we take $K(x,v):=\frac{q(x,\theta)}{\omega}$, where
$\theta=\frac{v}{|v|}$ and $q(x,\theta)$, with $\theta\in\mathbb{S}^{n-1}$, is
the orientational distribution of the fibers, normalized by $\omega=s^{n-1}$.
This function encodes the information about the individual brain structure,
obtained by diffusion tensor imaging (DTI). A concrete choice of $q(x,\theta)$
will be provided in Section 3. We also assume the tissue to be
undirected111for a motivation refer to [33], where the model developed and its
numerical simulations hint on the brain tissue being undirected, hence
$q(x,\theta)=q(x,-\theta)$ for all $x\in\mathbb{R}^{n}$. For later reference
we introduce the notations

$\displaystyle\mathbb{E}_{q}(x):=\int_{\mathbb{S}^{n-1}}\theta
q(x,\theta)d\theta$
$\displaystyle\mathbb{V}_{q}(x):=\int_{\mathbb{S}^{n-1}}(\theta-\mathbb{E}_{q})\otimes(\theta-\mathbb{E}_{q})\
q(x,\theta)d\theta$

for the mean fiber orientation and the variance-covariance matrix for the
orientation distribution of tissue fibers, respectively. Notice that the above
symmetry of $q$ implies $\mathbb{E}_{q}=0$.

As in [14, 26], the term $\mathcal{P}\,p$ is used to describe the
proliferation process triggered by cell receptor binding with the brain tissue
and enhanced in the present formulation by the capillaries that provide the
necessary nutrients and oxygen. We assume it to be given by:

$\mathcal{P}\,p:=\mu_{M}(M,W,N_{e})\int_{Z}\chi(x,z,z^{{}^{\prime}})\dfrac{Q(x)}{Q^{*}}p(t,x,v,z^{{}^{\prime}})dz^{{}^{\prime}}\,.$

The proliferation rate $\mu(M,W,N_{e})$, depends on the macroscopic densities
of tumor cells and ECs ($M$ and $W$, respectively), but is also negatively
influenced by necrotic matter, of density $N_{e}$. In the integral operator,
the kernel $\chi(x,z,z^{\prime})$ characterizes the transition from state
$z^{\prime}$ to state $z$ during such proliferation-initiating interaction at
position $x$. No further conditions are required on $\chi$; we only assume
that the nonlinear proliferative operator is uniformly bounded in the
$L^{2}$-norm, which is reasonable because of space-limited cell division.

The last terms in (4) describe the decay of glioma cells due to radio- and
chemotherapy. Thereby, $R_{M}(d_{r})$ describes the radiation effect on tumor
cells; its expression relates to the survival fraction $S(d_{r},\alpha)$ as
follows:

$R_{M}(d_{r})=1-S(d_{r},\alpha_{M})$

with $d_{r}$ denoting the radiation dose and $\alpha_{M}$ being a parameter
which refers (within the LQ model mentioned in Section 2.1) to the lesions
caused by irradiation. Chemotherapeutic effects inducing glioma depletion are
quantified by $C_{M}(k_{c})$, where the parameter $k_{c}$ represents the
killing rate, related to the effectiveness of chemotherapy on the tumor cells.
We consider that the same quantity of drug is administered once a day and that
it has each time the same effectiveness, thus

$C_{M}(k_{c}):=k_{c}(t)=\sum_{i=1}^{\nu}k_{c}\eta_{\delta}(t-t_{i}),$ (6)

with $t_{i}$ denoting the time moments at which the drug is given and $\nu$
representing the number of times (days) allocated to this chemotherapy.
Thereby, $\eta_{\delta}$ is a smooth function with unit mass and support in
$(-\delta,\delta)$, with $\delta\ll 1$. In particular, we indicate with
$L_{M}(d_{r},k_{c}):=(R_{M}(d_{r})+C_{M}(k_{c}))$ the overall loss term for
the tumor population.

The KTE for ECs is given by:

$\partial_{t}w+\nabla_{x}\cdot\,(\vartheta
w)+\partial_{u}(\Gamma(u,H,d_{b})w)=\mathcal{L}_{w}[\gamma(u)]w+\mathcal{P}_{w}\,w-R_{W}(d_{r})w\,,$
(7)

where the turning operator $\mathcal{L}_{w}[\gamma(u)]w$ describes changes in
the orientation of ECs due to a chemotactic response to the concentration of
growth factors produced by tumor cells. It is well-established, in fact, (see
e.g. [24]) that tumor cells produce angiogenic signals acting as
chemoattractants for ECs. These growth factors (mainly VEGF) stimulate
endothelial cells to migrate toward the tumor and proliferate, a process
called angiogenesis. Typically, the blood vessels produced within tumors and
in their immediate proximity are distorted and enlarged, with erratic blood
flux or microhemorrhages [24]. A switch from low to high vascularization is
considered to be an essential indicator of tumor progression [44].

Similarly to the description of the turning operator for glioma cells,
$\mathcal{L}_{w}[\gamma]w$ is an integral operator of the form

$\mathcal{L}_{w}[\gamma(u)]w=-\gamma(u)w+\gamma(u)\int_{\Theta}\dfrac{1}{|\Theta|}w(\vartheta^{{}^{\prime}})d\vartheta^{{}^{\prime}}\,.$
(8)

In this case, for the turning kernel characterizing ECs reorientations we
simply took a uniform distribution, while
$\gamma(u):=\gamma_{0}-\gamma_{1}u\geq 0$, which represents the EC turning
rate, depends on the microscopic variable $u$ and, consequently, on the
concentration of VEGF, with $\gamma_{0}$ and $\gamma_{1}$ being some positive
constants.

$\mathcal{P}_{w}\,w$ is the proliferation term related to ECs. It describes
proliferation as a result of interactions between ECs and their
chemoattractant VEGF

$\mathcal{P}_{w}\,w:=\mu_{W}(W,Q,d_{b})\int_{U}\chi_{W}(x,u,u^{{}^{\prime}})\dfrac{H(t,x)}{H_{c,0}}w(t,x,\vartheta,u^{{}^{\prime}})du^{{}^{\prime}}\,.$

The proliferation rate $\mu_{W}(W,Q,d_{b})$ depends on the macroscopic density
of ECs and on the density of brain tissue $Q$. We also consider that the anti-
angiogenic therapy might also affect the proliferation of ECs more directly,
in addition to the effect on the receptor binding kinetics. Thus, we include
the dependency of the proliferation rate on the anti-angiogenic dose $d_{b}$.

Finally, $R_{W}(d_{r})$ describes the depletion of ECs due to radiotherapy. We
take $R_{W}(d_{r}):=1-S(d_{r},\alpha_{W})$, with $\alpha_{W}>0$ a parameter
related to the lesions caused on ECs by ionizing radiation.

More details about the concrete choices of the coefficient functions involved
in (4) and (7) will be provided in Section 3.

###  2.3 Parabolic scaling of the mesoscopic model

We deduce effective equations for the macroscopic dynamics of ECs and glioma
cells, represented by the quantities $W$ and $M$, respectively. The PDEs for
the VEGF concentration $H$, the density of healthy tissue $Q$ and of necrotic
matter $N_{e}$ are stated directly at the macroscopic level. We rescale the
time and space variables as $t\to\varepsilon^{2}t$ and $x\to\varepsilon x$ and
we assume that the proliferation and death terms in (4) and (7) will be scaled
by $\varepsilon^{2}$ in order to account for the mitotic and apoptotic events
taking place on a much larger time scale than migration. Hence, (4) and (7)
become:

$\begin{split}\varepsilon^{2}\partial_{t}p+&\varepsilon\nabla_{x}\cdot(vp)-\partial_{z}\left(\left(zB_{p}(W,Q)-\frac{k_{1}^{-}}{B_{p}(W,Q)^{2}}\left(\varepsilon
k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}v\cdot\nabla W+\varepsilon
k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}v\cdot\nabla
Q+\varepsilon^{2}\left(\bar{F}_{W}+\bar{F}_{Q}\right)\right)\right)p\right)\\\\[5.69046pt]
&=\mathcal{L}_{p}[\lambda(z)]p+\varepsilon^{2}(\mathcal{P}-L_{M}(d_{r},k_{c}))p\end{split}$
(9) $\varepsilon^{2}\partial_{t}w+\varepsilon\nabla_{x}\cdot\,(\vartheta
w)-\partial_{u}\left(\left(uB_{w}(H)-\frac{\frac{k_{H}^{+}}{H_{c,0}}k_{H}^{-}}{B_{w}(H)^{2}}\left(\varepsilon\vartheta\cdot\nabla
H+\varepsilon^{2}\partial_{t}H\right)\right)w\right)=\mathcal{L}_{w}[\gamma(u)]w+\varepsilon^{2}(\mathcal{P}_{w}-R_{W}(d_{r}))w\,.$
(10)

As in [13, 14, 15], we define the following moments:

$\begin{split}&m(t,x,v)=\int_{Z}p(t,x,v,z)dz,\quad\quad
m^{z}(t,x,v)=\int_{Z}z\,p(t,x,v,z)dz,\\\\[5.69046pt]
&M(t,x)=\int_{V}m(t,x,v)dv,\quad\quad\quad
M^{z}(t,x)=\int_{V}m^{z}(t,x,v)dv\\\\[5.69046pt]
&\bar{w}(t,x,\vartheta)=\int_{U}w(t,x,\vartheta,u)du,\quad\,\,\,\bar{w}^{u}(t,x,\vartheta)=\int_{U}u\,w(t,x,\vartheta,u)du\\\\[5.69046pt]
&W(t,x)=\int_{\Theta}\bar{w}(t,x,\vartheta)d\vartheta,\quad\,\quad\quad
W^{u}(t,x)=\int_{\Theta}\bar{w}^{u}(t,x,\vartheta)d\vartheta\end{split}$

and neglect the higher order moments w.r.t. the variables $z$ and $u$, in
virtue of the subcellular dynamics being much faster than the events on the
higher scales, hence $z\ll 1$ and $u\ll 1$. We assume the function $p$ to be
compactly supported in the phase space $\mathbb{R}^{n}\times
V\times\mathit{Y}$ and $w$ to be compactly supported in
$\mathbb{R}^{n}\times\Theta\times\mathit{Z}$.

We first integrate equation (9) with respect to $z$, getting the following
equation for $m(t,x,v)$:

$\begin{split}\varepsilon^{2}\partial_{t}m+\varepsilon\,\nabla_{x}\cdot(vm)&=-\lambda_{0}\left(m-\dfrac{q}{\omega}M\right)+\lambda_{1}\left(m^{z}-\dfrac{q}{\omega}M^{z}\right)+\varepsilon^{2}\int_{Z}\mu_{M}(M,W,N_{e})\int_{Z}\chi(z,z^{\prime},x)p(z^{\prime})\dfrac{Q}{Q^{*}}dz^{\prime}dz\\\\[5.69046pt]
&-\varepsilon^{2}L_{M}(d_{r},k_{c})m\,.\end{split}$

Using the fact that $\chi(z,z^{\prime},x)$ is a probability kernel with
respect to $z$ for all $(x,z^{{}^{\prime}})$, the previous equation for
$m(t,x,v)$ reduces to:

$\varepsilon^{2}\partial_{t}m+\varepsilon\nabla_{x}\cdot(vm)=-\lambda_{0}\left(m-\dfrac{q}{\omega}M\right)+\lambda_{1}\left(m^{z}-\dfrac{q}{\omega}M^{z}\right)+\varepsilon^{2}\mu_{M}(M,W,N_{e})\dfrac{Q}{Q^{*}}m-\varepsilon^{2}L_{M}(d_{r},k_{c})m.$
(11)

Then, we multiply equation (9) by $z$ and integrate it w.r.t. $z$, obtaining

$\begin{split}\varepsilon^{2}\partial_{t}m^{z}&+\varepsilon\nabla_{x}\cdot(vm^{z})\\\\[5.69046pt]
&-\int_{Z}z\partial_{z}\left[\left(zB_{p}(W,Q)-\dfrac{k^{-}_{1}}{B_{p}(W,Q)^{2}}\left(\varepsilon\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}v\cdot\nabla
W+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}v\cdot\nabla
Q\right)+\varepsilon^{2}(\bar{F}_{W}+\bar{F}_{Q})\right)\right)p(z)\right]dz\\\\[5.69046pt]
&=\int_{Z}z\mathcal{L}_{p}[\lambda(z)]p(z)dz+\varepsilon^{2}\int_{Z}z\,\mu_{M}(M,W,N_{e})\int_{Z}\chi(z,z^{\prime},x)p(z^{\prime})\dfrac{Q}{Q^{*}}dz^{\prime}dz\,-\varepsilon^{2}L_{M}(d_{r},k_{c})m^{z}\,.\end{split}$

The calculation of the integral term on the left hand side leads to the
following equation for $m^{z}(t,x,v)$:

$\displaystyle\varepsilon^{2}\partial_{t}m^{z}+\varepsilon\nabla_{x}\cdot(vm^{z})=$
$\displaystyle-
B_{p}(W,Q)m^{z}+\dfrac{k^{-}_{1}}{B_{p}(W,Q)^{2}}\left(\varepsilon\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}v\cdot\nabla
W+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}v\cdot\nabla
Q\right)+\varepsilon^{2}\left(\bar{F}_{W}+\bar{F}_{Q}\right)\right)m$
$\displaystyle-\lambda_{0}\left(m^{z}-\dfrac{q}{\omega}M^{z}\right)+\varepsilon^{2}\mu_{M}(M,W,N_{e})\int_{Z}\int_{Z}\,z\chi(z,z^{\prime},x)p(z^{\prime})\dfrac{Q}{Q^{*}}dz^{\prime}dz-\varepsilon^{2}L_{M}(d_{r},k_{c})m^{z}.$
(12)

Applying the same procedure to equation (10), but integrating with respect to
$u$, we obtain the following equation for the moment $\bar{w}$:

$\varepsilon^{2}\partial_{t}\bar{w}+\varepsilon\nabla_{x}\cdot(\vartheta\bar{w})=-\gamma_{0}\left(\bar{w}-S_{n}^{\sigma}W\right)+\gamma_{1}\left(\bar{w}^{u}-S_{n}^{\sigma}W^{u}\right)+\varepsilon^{2}\mu_{W}(W,Q,d_{b})\dfrac{H}{H_{c,0}}\,\bar{w}-\varepsilon^{2}R_{W}(d_{r})\bar{w}\,,$
(13)

assuming $\chi_{W}(x,u,u^{{}^{\prime}})$ to be a probability kernel with
respect to $u$ for all $(x,u^{{}^{\prime}})$. Multiplying by $u$ and
integrating (10) again w.r.t. $u$, we obtain the equation for $\bar{w}^{u}$:

$\begin{split}\varepsilon^{2}\partial_{t}\bar{w}^{u}&+\varepsilon\nabla_{x}\cdot(\vartheta\bar{w}^{u})-\int_{U}u\partial_{u}\left[\left(uB_{w}(H)-\dfrac{\frac{k_{H}^{+}}{H_{c,0}}k_{H}^{-}}{B_{w}(H)^{2}}(\varepsilon\vartheta\cdot\nabla
H+\varepsilon^{2}\bar{E})\right)w(u)\right]du=\\\\[5.69046pt]
&=\int_{U}u\mathcal{L}_{w}[\gamma(u)]w(u)du+\varepsilon^{2}\int_{U}\,u\,\mu_{W}(W,Q,d_{b})\int_{U}\chi_{W}(u,u^{\prime},x)\dfrac{H}{H_{c,0}}w(u^{\prime})du^{\prime}du-\varepsilon^{2}R_{W}(d_{r})\bar{w}^{u}\,,\end{split}$

After the calculation of the integral on the left hand side, it reduces to

$\begin{split}\varepsilon^{2}\partial_{t}\bar{w}^{u}+&\varepsilon\nabla_{x}\cdot(\vartheta\bar{w}^{u})+B_{w}(H)\bar{w}^{u}-\dfrac{\frac{k_{H}^{+}}{H_{c,0}}k_{H}^{-}}{B_{w}(H)^{2}}(\varepsilon\vartheta\cdot\nabla
H\,+\varepsilon^{2}\bar{E})\bar{w}=-\gamma_{0}\left(\bar{w}^{u}-S_{n}^{\sigma}W^{z}\right)+\\\\[5.69046pt]
&+\varepsilon^{2}\mu_{W}(W,Q,d_{b})\int_{U}\int_{U}\chi_{W}(u,u^{\prime},x)\dfrac{H}{H_{c,0}}w(u^{\prime})du^{\prime}du-\varepsilon^{2}R_{W}(d_{r})\bar{w}^{u}\,.\end{split}$
(14)

We consider Hilbert expansions for the previously introduced moments:

$\begin{split}&m(t,x,v)=\sum_{k=0}^{\infty}\varepsilon^{k}m_{k},\quad\quad
m^{z}(t,x,v)=\sum_{k=0}^{\infty}\varepsilon^{k}m^{z}_{k}\quad\quad
M(t,x)=\sum_{k=0}^{\infty}\varepsilon^{k}M_{k},\quad\quad
M^{z}(t,x)=\sum_{k=0}^{\infty}\varepsilon^{k}M^{z}_{k}\\\\[5.69046pt]
&\bar{w}(x,t,\vartheta)=\sum_{k=0}^{\infty}\varepsilon^{k}\bar{w}_{k},\quad\quad\bar{w}^{u}(x,t,\vartheta)=\sum_{k=0}^{\infty}\varepsilon^{k}\bar{w}^{u}_{k},\quad\quad
W(x,t)=\sum_{k=0}^{\infty}\varepsilon^{k}W_{k},\quad\quad
W^{u}(t,x)=\sum_{k=0}^{\infty}\varepsilon^{k}W^{u}_{k}\,.\end{split}$

For the subsequent calculations it will be useful to Taylor-expand the
proliferation coefficient functions involving $W$ or $M$ in the scaled
equations (11)- (14):

$\displaystyle\mu(M,W,N_{e})=\mu(M_{0},W_{0},N_{e})+\partial_{M}\mu(M_{0},W_{0},N_{e})(M-M_{0})+\partial_{W}\mu(M_{0},W_{0},N_{e})(W-W_{0})+O(\varepsilon^{2}),$
$\displaystyle\mu_{W}(W,Q,d_{b})=\mu_{W}(W_{0},Q,d_{b})+\partial_{W}\mu_{W}(W_{0},Q,d_{b})(W-W_{0})+O(\varepsilon^{2}).$

Moreover, we observe that

$B_{p}(W,Q)=B_{p}(W_{0},Q)+k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}(W-W_{0})+O(|W-W_{0}|^{2})$

and

$\displaystyle\dfrac{1}{(B_{p}(W,Q))^{2}}$
$\displaystyle=\dfrac{1}{(B_{p}(W_{0},Q))^{2}}-2k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}\dfrac{1}{(B_{p}(W_{0},Q))^{3}}(W-W_{0})+O(|W-W_{0}|^{2})\,.$

Then, equating the powers of $\varepsilon$ in the scaled equations (11)-(14),
we obtain:

$\varepsilon^{0}$ terms:

$\displaystyle
0=-\lambda_{0}\left(m_{0}-\dfrac{q}{\omega}M_{0}\right)+\lambda_{1}\left(m_{0}^{z}-\dfrac{q}{\omega}M_{0}^{z}\right)$
(15) $\displaystyle
0=-B_{p}(W_{0},Q)m_{0}^{z}-\lambda_{0}\left(m_{0}^{z}-\dfrac{q}{\omega}M_{0}^{z}\right)$
(16) $\displaystyle
0=-\gamma_{0}(\bar{w}_{0}-S_{n}^{\sigma}W_{0})+\gamma_{1}(\bar{w}_{0}^{u}-S_{n}^{\sigma}W_{0}^{u})$
(17) $\displaystyle
0=-B_{w}(H)\bar{w}_{0}^{u}-\gamma_{0}\left(\bar{w}_{0}^{u}-S_{n}^{\sigma}W_{0}^{u}\right)$
(18)

where
$S_{n}^{\sigma}:=\frac{1}{|\Theta|}=\frac{\sigma^{1-n}}{|\mathbb{S}^{n-1}|}$.

$\varepsilon^{1}$ terms:

$\displaystyle\nabla_{x}\cdot(vm_{0})=-\lambda_{0}\left(m_{1}-\dfrac{q}{\omega}M_{1}\right)+\lambda_{1}\left(m_{1}^{z}-\dfrac{q}{\omega}M_{1}^{z}\right)$
(19)
$\displaystyle\nabla_{x}\cdot(vm_{0}^{z})+B_{p}(W_{0},Q)m_{1}^{z}+k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}W_{1}m_{0}^{z}-\dfrac{k_{1}^{-}}{(B_{p}(W_{0},Q)^{2}}\,\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}v\cdot\nabla
W_{0}+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}v\cdot\nabla Q\right)m_{0}=$ (20)
$\displaystyle=-\lambda_{0}\left(m_{1}^{z}-\dfrac{q}{\omega}M_{1}^{z}\right)$
(21)
$\displaystyle\nabla_{x}\cdot(\vartheta\bar{w}_{0})=-\gamma_{0}(\bar{w}_{1}-S_{n}^{\sigma}W_{1})+\gamma_{1}(\bar{w}_{1}^{u}-S_{n}^{\sigma}W_{1}^{u})$
(22)
$\displaystyle\nabla_{x}\cdot(\vartheta\bar{w}_{0}^{u})+B_{w}(H)\bar{w}_{1}^{u}-\dfrac{\frac{k_{H}^{+}}{H_{c,0}}k_{H}^{-}}{B_{w}(H)^{2}}\,\vartheta\cdot\nabla
H\bar{w}_{0}=-\gamma_{0}\left(\bar{w}_{1}^{u}-S_{n}^{\sigma}W_{1}^{u}\right)$
(23)

$\varepsilon^{2}$ terms:

$\displaystyle\partial_{t}m_{0}+\nabla_{x}\cdot(vm_{1})=-\lambda_{0}\left(m_{2}-\dfrac{q}{\omega}M_{2}\right)+\lambda_{1}\left(m_{2}^{z}-\dfrac{q}{\omega}M_{2}^{z}\right)+\mu_{M}(M_{0},W_{0},N_{e})\dfrac{Q}{Q^{*}}\,m_{0}-L_{M}(d_{r},k_{c})m_{0}$
(24) $\displaystyle\partial_{t}\bar{w}_{0}+\nabla_{x}\cdot(\vartheta
w_{1})=-\gamma_{0}(\bar{w}_{2}+S_{n}^{\sigma}W_{2})+\gamma_{1}(\bar{w}^{u}_{2}-S_{N}^{\sigma}W_{2}^{u})+\mu_{W}(W_{0},Q,d_{b})\dfrac{H}{H_{c,0}}\bar{w}_{0}-R_{W}(d_{r})\bar{w}_{0}$
(25)

Integrating (16) with respect to $v$, we have

$0=-M_{0}^{z}B_{p}(W_{0},Q)$

i.e.,

$M_{0}^{z}=0\,.$ (26)

This, together with (16), leads to

$m_{0}^{z}=0\,.$ (27)

From (15) and (16) we obtain:

$0=-\lambda_{0}m_{0}+\lambda_{0}\dfrac{q}{\omega}M_{0}\quad\Longrightarrow
m_{0}=\dfrac{q}{\omega}M_{0}\,.$ (28)

From equation (18) we have, by integrating w.r.t. $\vartheta$,

$0=-W_{0}^{u}B_{w}(H)$

i.e.,

$W_{0}^{u}=0\,.$ (29)

This, together with (18), leads to

$\bar{w}_{0}^{u}=0\,.$ (30)

From (17) and (18) we obtain:

$0=-\gamma_{0}\bar{w}_{0}+\gamma_{0}S_{n}^{\sigma}W_{0}\quad\Longrightarrow\bar{w}_{0}=S_{n}^{\sigma}W_{0}\,,$
(31)

thus $\bar{w}_{0}$ depends on the (constant) speed $\sigma$, but not on the
direction $\vartheta\in\mathbb{S}^{n-1}$.

Now we turn to the equations stemming from the $\varepsilon^{1}$-terms.
Integrating (23) w.r.t $\vartheta$ and using (31) and (30), we have

$0=-B_{w}(H)W_{1}^{u}+\dfrac{k_{H}^{-}\frac{k_{H}^{+}}{H_{c,0}}}{B_{w}(H)^{2}}\int_{\Theta}\vartheta\,\bar{w}_{0}\,\nabla
H\,d\vartheta\,,$

i.e.,

$W_{1}^{u}=0$ (32)

due to (31). Plugging this result into (23) we get

$\begin{split}&0=-(B_{w}(H)+\gamma_{0})\bar{w}_{1}^{u}+\dfrac{k_{H}^{-}\frac{k_{H}^{+}}{H_{c,0}}}{B_{w}(H)^{2}}\vartheta\cdot\nabla
H\,\bar{w}_{0}\\\\[8.5359pt]
&\Rightarrow\quad\bar{w}_{1}^{u}=\dfrac{k_{H}^{-}\frac{k_{H}^{+}}{H_{c,0}}}{(B_{w}(H)+\gamma_{0})\,B_{w}(H)^{2}}\,\vartheta\cdot\nabla
H\,\bar{w}_{0}\,.\end{split}$ (33)

Considering (22) and using (32), we can derive

$\displaystyle\nabla_{x}\cdot(\vartheta\bar{w}_{0})=-\gamma_{0}\left(\bar{w}_{1}-S_{n}^{\sigma}W_{1}\right)+\gamma_{1}\bar{w}_{1}^{u}$
$\displaystyle\Rightarrow-\gamma_{0}\bar{w}_{1}=\nabla_{x}\cdot(\vartheta\bar{w}_{0})-\gamma_{0}S_{n}^{\sigma}W_{1}-\gamma_{1}\bar{w}_{1}^{u}$
$\displaystyle\Rightarrow\bar{\mathcal{L}}_{W}[\gamma_{0}]\bar{w}_{1}:=-\gamma_{0}\bar{w}_{1}+\gamma_{0}S_{n}^{\sigma}W_{1}=\nabla_{x}\cdot(\vartheta\bar{w}_{0})-\gamma_{1}\bar{w}_{1}^{u}\,.$
(34)

In order to get an explicit expression for $w_{1}$, we would like to invert
the operator $\bar{\mathcal{L}}_{W}[\gamma_{0}]$. Similarly to [39] we define
it on the weighted $L^{2}$-space $L^{2}_{S_{n}^{\sigma}}(\Theta)$, in which
the measure $d\vartheta$ is weighted by $S_{n}^{\sigma}(\vartheta)$. In
particular, $L^{2}_{S_{n}^{\sigma}}(\Theta)$ can be decomposed as
$L^{2}_{S_{n}^{\sigma}}(\Theta)=<S_{n}^{\sigma}>\oplus<S_{n}^{\sigma}>^{\perp}$.
Due to the properties of the chosen turning kernel,
$\bar{\mathcal{L}}_{W}[\gamma_{0}]$ is a compact Hilbert-Schmidt operator with
kernel $<S_{n}^{\sigma}>$. We can therefore calculate its pseudo-inverse on
$<S_{n}^{\sigma}>^{\perp}$. Thus, to determine $\bar{w}_{1}$ from (34) we
check the solvability condition, which holds for the above results on
$\bar{w}_{0}$. Hence we obtain from (34) and (33)

$\bar{w}_{1}=-\dfrac{1}{\gamma_{0}}\left[\nabla_{x}\cdot(\vartheta\bar{w}_{0})-\gamma_{1}\dfrac{\frac{k_{H}^{+}}{H_{c,0}}k_{H}^{-}}{(B_{w}(H)+\gamma_{0})B_{w}(H)^{2}}\vartheta\cdot\nabla
H\,\bar{w}_{0}\right]$ (35)

and

$W_{1}=0\,.$ (36)

On the other hand, integrating (21) with respect to $v$, we have upon using
(28) and (27):

$0=-B_{p}(W_{0},Q)M_{1}^{z}+\dfrac{k^{-}_{1}}{B_{p}(W_{0},Q)^{2}}M_{0}s^{n}\mathbb{E}_{q}\,\cdot\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}\nabla
W_{0}+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}\nabla Q\right)\,.$

Using the fact that the tissue is undirected, i.e. $\mathbb{E}_{q}=0$, the
previous equation leads to:

$M_{1}^{z}=0\,.$ (37)

Now plugging these results into (21) we get

$\begin{split}&0=-(B_{p}(W_{0},Q)+\lambda_{0})m_{1}^{z}+\dfrac{k^{-}_{1}}{B_{p}(W_{0},Q)^{2}}\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}v\cdot\nabla
W_{0}+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}v\cdot\nabla Q\right)\,m_{0}\\\\[8.5359pt]
&\Rightarrow\quad
m_{1}^{z}=\dfrac{k^{-}_{1}}{(B_{p}(W_{0},Q)+\lambda_{0})\,B_{p}(W_{0},Q)^{2}}\,\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}v\cdot\nabla
W_{0}+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}v\cdot\nabla Q\right)\,m_{0}\,.\end{split}$
(38)

Using (19) we derive

$\displaystyle\nabla_{x}\cdot(vm_{0})=-\lambda_{0}\left(m_{1}-\dfrac{q}{\omega}M_{1}\right)+\lambda_{1}m_{1}^{z}$
$\displaystyle\Rightarrow-\lambda_{0}m_{1}=\nabla_{x}\cdot(vm_{0})-\lambda_{0}\dfrac{q}{\omega}M_{1}-\lambda_{1}m_{1}^{z}$
$\displaystyle\Rightarrow\bar{\mathcal{L}}_{M}[\lambda_{0}]m_{1}:=-\lambda_{0}m_{1}+\lambda_{0}\dfrac{q}{\omega}M_{1}=\nabla_{x}\cdot(vm_{0})-\lambda_{1}m_{1}^{z}\,.$
(39)

As done for the turning operator of the endothelial cells, in order to get an
explicit expression for $m_{1}$, we would like to invert the operator
$\bar{\mathcal{L}}_{M}[\lambda_{0}]$. Again, similarly to [39], we define it
on the weighted $L^{2}$-space $L^{2}_{q}(V)$, in which the measure $dv$ is
weighted by $q(x,\theta)/\omega$. In particular, $L^{2}_{q}(V)$ can be
decomposed as $L^{2}_{q}(V)=<q/\omega>\oplus<q/\omega>^{\perp}$. Due to the
properties of the chosen turning kernel, $\bar{\mathcal{L}}_{M}[\lambda_{0}]$
is a compact Hilbert-Schmidt operator with kernel $<q/\omega>$. We can
therefore calculate its pseudo-inverse on $<q/\omega>^{\perp}$.

Thus, to determine $m_{1}$ from (39) we need to check the solvability
condition

$\int_{V}\left[\nabla_{x}\cdot(vm_{0})-\lambda_{1}m_{1}^{z}\right]dv=0\,.$

This holds thanks to the above results and to the symmetry of $q(x,\theta)$.
Therefore, we obtain from (39) and (38)

$m_{1}=-\dfrac{1}{\lambda_{0}}\left[\nabla_{x}\cdot(vm_{0})-\lambda_{1}\dfrac{k^{-}_{1}}{(B_{p}(W_{0},Q)+\lambda_{0})B_{p}(W_{0},Q)^{2}}\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}v\cdot\nabla
W_{0}+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}v\cdot\nabla Q\right)\,m_{0}\right]$ (40)

and

$M_{1}=0\,.$ (41)

Integrating equation (24) with respect to $v$ we get

$\partial_{t}M_{0}+\int_{V}\nabla_{x}\cdot(vm_{1})dv=\mu_{M}(M_{0},W_{0},N_{e})\dfrac{Q(x)}{Q^{*}}M_{0}-L_{M}(d_{r},k_{c})M_{0}\,,$
(42)

where

$\begin{split}\int_{V}\nabla_{x}\cdot(vm_{1})dv&=\int_{V}\nabla_{x}\cdot\Bigg{[}\,v\Bigg{(}-\dfrac{1}{\lambda_{0}}\bigg{(}\nabla_{x}\cdot(vm_{0})\\\\[8.5359pt]
&-\lambda_{1}\dfrac{k^{-}_{1}}{(B_{p}(W_{0},Q)+\lambda_{0})B_{p}(W_{0},Q)^{2}}\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}v\cdot\nabla
W_{0}+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}v\cdot\nabla
Q\right)\,m_{0}\bigg{)}\Bigg{)}\Bigg{]}dv\\\\[8.5359pt]
&=\nabla_{x}\cdot\left[\int_{V}-\dfrac{1}{\lambda_{0}}v\otimes
v\,\nabla_{x}\cdot\left(\dfrac{q}{\omega}M_{0}\right)\right]dv\\\\[8.5359pt]
&+\nabla_{x}\cdot\Bigg{[}\dfrac{\lambda_{1}k_{1}^{-}}{\omega\lambda_{0}(B_{p}(W_{0},Q)+\lambda_{0})B_{p}(W_{0},Q)^{2}}\int_{V}v\otimes
v\,q(\theta,x)dv\,\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}\nabla
W_{0}+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}\nabla
Q\right)\,M_{0}\Bigg{]}\,.\end{split}$

With the notation

$\mathbb{D}_{T}(x):=\dfrac{1}{\omega}\int_{V}v\otimes
v\,q(x,\theta)dv=s^{2}\int_{\mathbb{S}^{n-1}}\theta\otimes\theta\,q(x,\theta)d\theta=s^{2}\mathbb{V}_{q}(x)$
(43)

we obtain from (42) the following macroscopic equation for $M_{0}(t,x)$:

$\begin{split}\partial_{t}M_{0}&-\nabla_{x}\cdot\left[\dfrac{1}{\lambda_{0}}\nabla_{x}\cdot\left(\mathbb{D}_{T}(x)M_{0}\right)\right]+\nabla_{x}\cdot\left[\dfrac{\lambda_{1}k_{1}^{-}}{\lambda_{0}(B_{p}(W_{0},Q)+\lambda_{0})B_{p}(W_{0},Q)^{2}}\mathbb{D}_{T}(x)\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}\nabla
W_{0}\,+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}\nabla Q\right)M_{0}\right]\\\\[8.5359pt]
&=\left(\mu_{M}(M_{0},W_{0},N_{e})\dfrac{Q}{Q^{*}}-L_{M}(d_{r},k_{c})\right)M_{0}\,.\end{split}$
(44)

Integrating (25) with respect to $\vartheta\in\Theta$ gives

$\partial_{t}W_{0}+\int_{\Theta}\nabla_{x}\cdot(\vartheta\bar{w}_{1})d\vartheta=\mu_{W}(W_{0},Q,d_{b})\dfrac{H}{H_{c,0}}W_{0}-R_{W}(d_{r})W_{0}$

where

$\begin{split}\int_{\Theta}\nabla_{x}\cdot(\vartheta\bar{w}_{1})d\vartheta&=\int_{\Theta}\nabla_{x}\cdot\left[\vartheta\left(-\dfrac{1}{\gamma_{0}}\left(\nabla_{x}\cdot(\vartheta\bar{w}_{0})-\gamma_{1}\dfrac{\frac{k_{H}^{+}}{H_{c,0}}k_{H}^{-}}{(B_{w}(H)+\gamma_{0})B_{w}(H)^{2}}\vartheta\cdot\nabla
H\,\bar{w}_{0}\right)\right)\right]\\\\[8.5359pt]
&=\nabla_{x}\cdot\left[\int_{\Theta}-\dfrac{1}{\gamma_{0}}\vartheta\otimes\vartheta\,\nabla_{x}\bar{w}_{0}\right]d\vartheta+\nabla_{x}\cdot\Bigg{[}\dfrac{\gamma_{1}\frac{k_{H}^{+}}{H_{c,0}}k_{H}^{-}}{\gamma_{0}(B_{w}(H)+\gamma_{0})B_{w}(H)^{2}}\int_{\Theta}\vartheta\otimes\vartheta\,d\vartheta\,\nabla
H\,W_{0}\Bigg{]}\,.\end{split}$

Recalling (31), this leads to the following macroscopic equation for the
density $W_{0}$ of endothelial cells:

$\partial_{t}W_{0}-\,\nabla_{x}\cdot\left(\mathbb{D}_{EC}\nabla_{x}W_{0}\right)+\nabla_{x}\cdot\left(\mathbb{D}_{EC}\dfrac{\gamma_{1}k_{H}^{+}k_{H}^{-}}{H_{c,0}(B_{w}(H)+\gamma_{0})B_{w}(H)^{2}}\nabla
H\,W_{0}\right)=\left(\mu_{W}(W_{0},Q,d_{b})\dfrac{H}{H_{c,0}}-R_{W}(d_{r})\right)W_{0}\,,$
(45)

where

$\mathbb{D}_{EC}:=\dfrac{\sigma^{2}}{n\gamma_{0}}\mathbb{I}_{n}\,.$ (46)

The two macroscopic equations obtained in (44) and (45) for the evolution of
glioma cells and ECs, respectively, are coupled with a PDE for the dynamics of
VEGF concentration:

$\partial_{t}H=D_{H}\Delta H+f(M_{0},W_{0},H)\,,$ (47)

where $D_{H}>0$ is the constant VEGF diffusion coefficient, while
$f(M_{0},W_{0},H)$ is a reaction term describing the processes of VEGF
production by tumor cells in hypoxic situations and its uptake by ECs. The
concrete form of the reaction term is defined in the next section. The
degradation of the healthy tissue is due to the acidity produced by tumor
cells and to the effects of radiotherapy. Thus, we describe the evolution of
healthy tissue by means of the following ODE:

$\partial_{t}Q=-d_{Q}\dfrac{M_{0}}{M_{c,0}}Q-R_{Q}(d_{r})Q\,,$ (48)

where $M_{c,0}$ denotes the reference density for tumor cells and $d_{Q}$ is
the tissue degradation rate. The term $R_{Q}(d_{r})Q$ collects the effects of
the radiation on the healthy tissue and is described with the linear quadratic
(LQ) model, which we previously introduced and used to formulate the survival
fractions $S_{W}$ and $S_{Q}$. Finally, the ODE describing the evolution of
necrotic matter takes into account both tissue degradation and therapy effects
on tumor, endothelial cells, and normal tissue, i.e,

$\partial_{t}N_{e}=L_{M}(d_{r},k_{c})M+R_{W}(d_{r})W_{0}+d_{Q}\dfrac{M_{0}}{M_{c,0}}Q+R_{Q}(d_{r})Q\,.$
(49)

In view of (41) and (36), the $\varepsilon$-correction terms for $M$ and $W$
can be left out and, neglecting the higher order moments, we get the following
PDE system characterizing the macroscopic evolution of the tumor under the
influence of tissue, vasculature, and growth factors:

$\left\\{\begin{array}[]{@{}l@{}}\partial_{t}M-\nabla_{x}\cdot\left[\dfrac{1}{\lambda_{0}}\nabla_{x}\cdot\left(\mathbb{D}_{T}(x)M\right)\right]+\nabla_{x}\cdot\left[\dfrac{\lambda_{1}k_{1}^{-}}{\lambda_{0}(B_{p}(W,Q)+\lambda_{0})B_{p}(W,Q)^{2}}\mathbb{D}_{T}(x)\left(k_{W}^{+}\dfrac{S_{W}}{W_{c,0}}\nabla
W\,+k_{Q}^{+}\dfrac{S_{Q}}{Q^{*}}\nabla Q\right)M\right]\vspace{0.7cm}\\\
\hskip
28.45274pt=\left(\mu_{M}(M,W,N_{e})\dfrac{Q(x)}{Q^{*}}-L_{M}(d_{r},k_{c})\right)M,\vspace{0.7cm}\\\
\partial_{t}W-\,\nabla_{x}\cdot\left(\mathbb{D}_{EC}\nabla_{x}W\right)+\nabla_{x}\cdot\Big{(}\mathbb{D}_{EC}\dfrac{\gamma_{1}k_{H}^{+}k_{H}^{-}}{H_{c,0}(B_{w}(H)+\gamma_{0})B_{w}(H)^{2}}\nabla
H\,W\Big{)}=\left(\mu_{W}(W,Q,d_{b})\dfrac{H}{H_{c,0}}-R_{W}(d_{r})\right)W,\vspace{0.7cm}\\\
\partial_{t}H=D_{H}\Delta H+f(M,W,H)\,,\vspace{0.7cm}\\\
\partial_{t}Q=-d_{Q}\dfrac{M}{M_{c,0}}Q-R_{Q}(d_{r})Q\,,\vspace{0.7cm}\\\
\partial_{t}N_{e}=L_{M}(d_{r},k_{c})M+R_{W}(d_{r})W+d_{Q}\dfrac{M}{M_{c,0}}Q+R_{Q}(d_{r})Q\,.\end{array}\right.$
(50)

with the tumor diffusion tensor $\mathbb{D}_{T}$ given in (43) and the EC
diffusion tensor $\mathbb{D}_{EC}$ given by (46). This system has been deduced
considering $x\in\mathbb{R}^{n}$ and it has to be supplemented with adequate
initial conditions, For the numerical simulations performed in the next
Section 4 we will consider the system to be set in a bounded, sufficiently
regular domain $\Omega\subset\mathbb{R}^{2}$ and endow it with no-flux
boundary conditions.

The next section is dedicated to specifying the remaining coefficients of the
above macroscopic system.

##  3 Assessment of coefficients

To determine the tumor diffusion tensor $\mathbb{D}_{T}(x)$ in (43) we need to
provide a concrete form for the (mesoscopic) orientational distribution of
tissue fibers $q(x,\theta)$. Several different choices are available in the
literature, see e.g. [6, 39]. As in [26] we use the orientation distribution
function (ODF)

$q(x,\theta)=\dfrac{1}{4\pi|\mathbb{D}_{W}(x)|^{\frac{1}{2}}(\theta^{T}(\mathbb{D}_{W}(x))^{-1}\theta)^{\frac{3}{2}}}\,,$
(51)

where $\mathbb{D}_{W}(x)$ is the water diffusion tensor obtained from
processing the (patient-specific) DTI data.

The growth rate $\mu(M,Q)$ can be also defined in different ways and it should
be motivated by biological evidence. As in e.g., [14], we choose a logistic-
like function to describe the growth self-limitation, enhancing this
limitation with a factor depending on the necrotic tissue $Ne$

$\mu(M,W,N_{e})=\mu_{M,0}\,\left(1-\dfrac{M}{K_{M}}-\dfrac{N_{e}}{K_{Ne}}\right)g(W)\,.$
(52)

Here, $K_{N}$ is the tumor carrying capacity, $K_{Ne}$ is the carrying
capacity of the necrotic matter, and $\mu_{M,0}$ is the proliferation rate.
The function $g(W)$ describes the influence of vasculature on tumor growth. It
can be simply a linear function of $W$, such as $g(W)=\frac{W}{W_{c,0}}$
describing tumor proliferation as an effect of interactions between tumor and
ECs, or can be assumed to have an expression like
$g(W)=(1+\frac{W}{W_{c,0}})$, describing ECs enhancing tumor proliferation in
better-vascularized regions. For the simulations presented in Section 4, we
take into account this second choice for the function $g(W)$.

Similarly, for the term $\mu_{W}(W,Q,d_{b})$ describing the proliferation rate
of ECs, apart from a limited growth term, we have to include the influence of
the healthy tissue on EC growth. In particular, the healthy tissue can be seen
as necessary to drive ECs proliferation and, therefore, binary interactions
between $W$ and $Q$ could be favorable

$\mu_{W}(W,Q,d_{b})=\mu_{W,0}(d_{b})\left(1-\dfrac{W}{K_{W}}\right)\dfrac{Q}{Q^{*}}\,,$
(53)

Here, $Q^{*}$ is the reference density for the tissue, $K_{W}$ is the carrying
capacity of the ECs, while $\mu_{W}(d_{b})$ is the proliferation rate,
depending on the dose of anti-angiogenic therapy $d_{b}$. This treatment, in
fact, also impairs the ability of ECs to grow and, therefore, $\mu_{W}(d_{b})$
can be described as a decreasing function of $d_{b}$.

Concerning radiotherapy, we already introduced the general expression related
to the LQ model $S(d_{r},\alpha):=\exp{(-\alpha d_{r}-\beta^{2}d_{r})}$;
precisely, the radiation-associated death for the two cell populations is
described as

$\begin{split}&R_{M}(d_{r}):=\sum_{i=1}^{\xi}(1-S(d_{r},\alpha_{M}))\eta_{\delta}(t-t_{i})\\\\[8.5359pt]
&R_{W}(d_{r}):=\sum_{i=1}^{\xi}(1-S(d_{r},\alpha_{W}))\eta_{\delta}(t-t_{i}),\end{split}$
(54)

where the total dose $d_{r}$ of the drug is given in smaller fractions, to
avoid toxic effects on healthy tissue. Thereby, $\xi$ is the number of
fractions, $\eta_{\delta}$ is a $C_{0}^{\infty}$ function with unit mass and
support in $(-\delta,\delta)$, $\delta\ll 1$, and $t_{i}$ denotes the time
instants at which ionizing radiation is applied to the patient. Moreover, the
effects of anti-angiogenic therapy affecting ECs migration and proliferation
are included in the binding/unbinding functions and in the proliferation term
as

$\displaystyle
k_{H}^{+}:=\bar{k}_{H}^{+}l^{+}(d_{b}):=\bar{k}_{H}^{+}\left(\dfrac{1}{2}+\dfrac{1}{2(1+d_{b}^{2})}\right)\,,$
(55) $\displaystyle
k_{H}^{-}:=\bar{k}_{H}^{-}l^{-}(d_{b}):=\bar{k}_{H}^{-}\left(1+d_{b}\right)\,,$
(56)
$\displaystyle\mu_{W,0}(d_{b}):=\bar{\mu}_{W,0}\left(\dfrac{1}{2}+\dfrac{1}{2(1+d_{b}^{2})}\right)\,.$
(57)

Finally, the reaction term in the PDE for the VEGF dynamics is chosen as

$f(M_{0},W_{0},H):=c_{1}\frac{M_{0}}{M_{c,0}}\left(\dfrac{W_{T}}{W_{c,0}}-\dfrac{W}{W_{c,0}}\right)_{+}\,\,-c_{2}\,\frac{W_{0}}{W_{c,0}}H\,,$
(58)

where the parameters $c_{1}>0$ and $c_{2}>0$ represent the production and
uptake rates of VEGFs, respectively, while $(\cdot)_{+}$ represents the
positive part. The parameter $W_{T}$ represents a threshold of EC density
referring to hypoxic condition. This means that when EC density is below such
threshold, the tumor cells keep producing VEGFs. The latter is produced in
order to boost angiogenesis and therewith associated nutrient supply. The ECs
are attracted by the VEGF gradient and uptake VEGF to sustain the growth of
the capillary network.

In Table 1 we report the (range of) values for the constant parameters
involved in system (50), as well as the references they were drawn from. In
particular, concerning the parameters $\lambda_{0}$ and $k^{+}_{Q,S}$ we
observe that the estimations are in good agreement with the ranges provided in
[31], where the clinical applicability of a similar multiscale model in
radiotherapy is studied for a collective of 3 GBM patients. Instead, following
the results presented in [7], we consider for the speed of tumor cells a
smaller range of variability than that provided in [31]. More details about
the therapy are provided in the explanations of the corresponding simulations.

Parameter | Description | Value (unit) | Source
---|---|---|---
$\lambda_{0}$ | turning frequency in $\mathcal{L}_{p}[\lambda(z)]$ |  $0.001$ (s${}^{-1})$ | [7, 46]
$\lambda_{1}$ | turning frequency in $\mathcal{L}_{p}[\lambda(z)]$ |  $0.001$ (s${}^{-1})$ | [7, 46]
$s$ | tumor cells speed |  $[0.0042,0.0084]\cdot 10^{-3}$ (mm $\cdot$ s${}^{-1})$ | [7, 11]
$k^{+}_{Q}$ | tumor-ECM attachment rate |  $0.034$ (s${}^{-1})$ | [35, 7]
$k^{+}_{W}$ | tumor-EC attachment rate |  $0.034$ (s${}^{-1})$ | [35]
$k^{-}$ | detachment rate |  $0.01$ (s${}^{-1})$ | [35]
$\mu_{M,0}$ | tumor proliferation rate |  $[3.5,9]\cdot 10^{-6}$ (s${}^{-1})$ | [29]
$K_{M}$ | tumor carrying capacity |  $\sim 10^{5}$ (cells $\cdot$ mm${}^{-3})$ | [1]
$K_{Ne}$ | necrotic matter carrying capacity |  $\sim 10^{5}$ (cells $\cdot$ mm${}^{-3})$ | [1]
$K_{Q}$ | healthy tissue carrying capacity |  $\sim 10^{5}$ (cells $\cdot$ mm${}^{-3})$ | [51]
$\alpha_{M}$ | single radiation track lesion on tumor |  $[0.04,0.11]$ (Gy${}^{-1})$ | [41, 3]
$\alpha_{Q}$ | single radiation track lesion on tissue |  $0.0025$ (Gy${}^{-1})$ | [32]
$\alpha_{W}$ | single radiation track lesion on ECs |  $0.0025$ (Gy${}^{-1})$ | this work
$\beta_{M}$ | two radiation tracks lesions on tumor |  $[0.006,0.019]$ (Gy${}^{-2})$ | [41, 3]
$\beta_{Q}$ | two radiation tracks lesions on tissue |  $0.00005$ (Gy${}^{-2})$ | [32]
$\beta_{W}$ | two radiation tracks lesions on ECs |  $0.0008$ (Gy${}^{-2})$ | [28]
$k_{c}$ | chemotherapy killing rate |  $0.23\cdot 10^{-6}$ (s${}^{-1})$ | [40]
$d_{r}$ | radiotherapy dose |  $60$ (Gy) | this work
$\gamma_{0}$ | turning frequency in $\mathcal{L}_{w}[\gamma(u)]$ |  $2\cdot 10^{-4}$ (s${}^{-1})$ | [52, 7]
$\gamma_{1}$ | turning frequency in $\mathcal{L}_{w}[\gamma(u)]$ |  $0.001$ (s${}^{-1})$ | this work
$\sigma$ | ECs speed |  $[0.0028,0.0069]\cdot 10^{-3}$ (mm $\cdot$ s${}^{-1})$ | [10, 7]
$\bar{k}_{H}^{+}$ | EC-VEGF attachment rate |  $0.03$ (s${}^{-1})$ | [20]
$\bar{k}_{H}^{-}$ | EC-VEGF detachment rate |  $0.001$ (s${}^{-1})$ | [20]
$K_{W}$ | EC carrying capacity |  $\sim 10^{5}$ (cells $\cdot$ mm${}^{-3})$ | [2]
$\bar{\mu}_{W,0}$ | EC proliferation rate |  $9.6\cdot 10^{-6}$ (s${}^{-1})$ | [7]
$W_{T}$ | EC threshold for hypoxia | $60\,\%\,K_{W}$ | this work
$d_{b}$ | anti-angiogenic therapy dose |  $10$ (mg $\cdot$ kg${}^{-1})$ | this work
$D_{H}$ | VEGF diffusion coefficient |  $[10^{-8},10^{-7}]$ (mm${}^{2}\cdot$ s${}^{-1})$ | [36]
$c_{1}$ | VEGF production rate |  $1.4\cdot 10^{-5}$ (s${}^{-1})$ | [21]
$c_{2}$ | VEGF consumption rate |  $2.3\cdot 10^{-4}$ (s${}^{-1})$ | [54]
$d_{Q}$ | tissue degradation rate (by acidity) |  $0.58\cdot 10^{-7}$ (s${}^{-1})$ | [7]
Table 1: Model parameters (dimensional)

##  4 Numerical simulations

We perform numerical simulations of the non-dimensionalized version of the
macroscopic setting (50) (see Appendix A for further details about the non-
dimensionalization) in 2D with the parameters listed in Table 1. For the
initial conditions we choose a Gaussian-like aggregate of tumor cells centered
in $(x_{0,M},y_{0,M})=(-17,5)$, in the right hemisphere of the brain slice
representing our illustrative computational domain $\Omega$ (in the following,
$(x,y)\in\Omega\subseteq\mathbb{R}^{2}$):

$M_{0}(x,y)=e^{-\frac{(x-x_{0,M})^{2}+(y-y_{0,M})^{2}}{8}}\,.$

For the initial distribution of normal tissue we consider the expression
proposed in [14]

$Q_{0}(x,y)=1-\dfrac{l_{c}^{3}(x,y)}{h^{3}}\,,$ (59)

where $h$ is the side length of one voxel of the DTI dataset and $l_{c}$ a
characteristic length estimated as

$l_{c}(x)=\sqrt{\frac{tr(\mathbb{D}_{W}(x))h^{2}}{4l_{1}}}$

with $l_{1}$ the leading eigenvalue of the diffusion tensor $\mathbb{D}_{W}$.
For the ECs, we consider the situation of several blood vessels located in
regions close to the tumor and where there is a high tissue density. From
usual MRI and DTI scans, in fact, it is not possible to distinguish between
vessels and ECM. Therefore, we assume the vasculature to be denser in the
areas where more tissue is available. Defining

$\displaystyle
W_{1}(x,y)=2\,e^{-\frac{(x+2)^{2}}{0.5}}\,\sin^{6}{\left(\dfrac{\pi}{8}\,y\right)}\qquad\forall\,y\in[-5,15]\,,$
(60) $\displaystyle
W_{2}(x,y)=2\,e^{-\frac{(x+12)^{2}}{0.5}}\,\sin^{4}{\left(\dfrac{\pi}{8}\,y\right)}\qquad\forall\,y\in[-35,-15]\,,$
(61) $\displaystyle
W_{3}(x,y)=2\,e^{-\frac{(x+5)^{2}+(y-15)^{2}}{0.75}}+2\,e^{-\frac{(x+5)^{2}+(y-20)^{2}}{0.75}}+2\,e^{-\frac{(x+10)^{2}+(y-20)^{2}}{0.75}}\,,$
(62)

we set $W_{0}(x,y)=W_{1}(x,y)+W_{2}(x,y)+W_{3}(x,y)$. For the VEGF initial
profile, we consider a Gaussian distribution, centered in the same point as
tumor cells ${(x_{0,H},y_{0,H})=(-17,5)}$ and given by:

$H_{0}(x,y)=0.5\,e^{-\frac{(x-x_{0,H})^{2}+(y-y_{0,H})^{2}}{4}}\,.$

Moreover, we assume that initially there is no necrotic matter, setting its
initial value to zero everywhere in the domain. Figure 1 shows the initial
conditions on the entire 2D brain slice, zooming then in the region
$\bar{\Omega}=[-35,5]\times[-15,25]\subseteq\Omega$.

Figure 1: Initial conditions of system (50). The five columns refer to
densities of glioma, ECs, VEGFs, necrotic matter, and healthy tissue,
respectively. The initial densities in the bottom row are visualized on the
zoomed domain $\bar{\Omega}=[-35,5]\times[-15,25]$.

The numerical simulations are performed with a self-developed code in Matlab
(MathWorks Inc., Natick, MA). The computational domain is a horizontal brain
slice reconstructed from the processing of an MRI scan. The macroscopic tensor
$\mathbb{D}_{T}(x)$ is pre-calculated using DTI data and the ODF for the fiber
distribution function (51) (see Section 3). The dataset was acquired at the
Hospital Galdakao-Usansolo (Galdakao, Spain), and approved by its Ethics
Committee: all the methods employed were in accordance with approved
guidelines. We consider a Galerkin finite element scheme for the spatial
discretization of the equations and an implicit Euler scheme for the time
discretization. We run several simulations to test different scenarios, which
can be summarized as follows.
System evolution without any therapy

  * (A)

We study the dynamics of the five species involved in the macroscopic setting
(50) in the absence of any treatment. This means that the therapy doses
($d_{r}$ and $d_{b}$) and the chemotherapy killing rate ($k_{c}$) are set to
zero and the tumor is let to grow and spread without any external impediment.

  * (B)

We analyze the effect of the subcellular dynamics included in the EC equation,
comparing the migration of these cells in three different cases. Precisely, we
deduce for ECs two possible different equations: a simple diffusion model
deduced from a kinetic equation where no subcellular dynamics are considered;
a diffusion-advection model, where the transport term is obtained as in [7].
We, then, compare the evolution of the ECs described in these two alternative
ways with the approach proposed in these notes.

Effects of treatment

  * (C)

We study the effect of concurrent radio- and chemotherapy, which is the
standard in clinical practice after tumor resection. The protocol we consider
accounts for irradiation with 60 Gy in fractions of 2 Gy, simultaneously
chemotherapy with temozolomide. The dosage of temozolomide (which is given
orally) is according to the body surface of the patient: 75 mg/m2 surface. The
treatment schedule is sketched in Figure 2. With respect to the complete Stupp
protocol [49, 48], we neglect the adjuvant chemotherapy, i.e. the so-called
”5/23 cycles”, which consists of 5 days of temozolomide (150-200 mg/m2) out of
every 28 days and repeated for 6 times. Since our main aim is to compare the
different treatment schemes and to point out the effects of the anti-
angiogenic therapy, the adjuvant chemotherapy would only be an additional (and
perhaps confusing) factor. Moreover, in the real case reported in Section 4.3,
the example patient actually did not receive the adjuvant chemotherapy.

Figure 2: Schematic representation of the therapeutic plan in experiment (C).
The blue arrows refer to the radiotherapy administration, which is given in
fractionated doses of $d_{r}=2$ Gy per day; the red arrow identifies the
period of chemotherapy administration, with a daily dose of $75$ mg $\cdot$
m-2.

  * (D)

Following the protocol reported in the trials NCT01209442 and
NCT01390948222https://www.clinicaltrials.gov, we combined the treatment based
on radio- and chemotherapy with an anti-angiogenic therapy with bevacizumab
given every two weeks. The dosage of bevacizumab is accordingly to the
bodyweight of the patient: 10 mg/kg. The treatment schedule is sketched in
Figure 3.

Figure 3: Schematic representation of the therapeutic plan in experiment (D).
The blue arrows refer to the radiotherapy administration, which is given in
fractionated doses of $d_{r}=2$ Gy per day; the yellow arrows refer to the
anti-angiogenic therapy administration, given on a 2-weeks base at dose
$d_{b}=10$ mg$\cdot$ kg-1; the red arrow identifies the period of chemotherapy
administration, with a daily dose of $75$ mg $\cdot$ m-2.

  * (E)

Considering the trial AVF3708, whose results have also been reported in [19],
we test an adjuvant, second-line therapy with bevacizumab (10 mg/kg) every 2
weeks on patients received prior radiotherapy and temozolomide. The treatment
schedule for this experiment is sketched in Figure 4.

Figure 4: Schematic representation of the therapeutic plan in experiment (E).
The blue arrows refer to the radiotherapy administration, which is given in
fractionated doses of $d_{r}=2$ Gy per day; the yellow arrows refer to the
adjuvant anti-angiogenic therapy administration, given on a 2-weeks base at
dose $d_{b}=10$ mg$\cdot$ kg-1; the red arrow identifies the period of
chemotherapy administration, with a daily dose of $75$ mg $\cdot$ m-2.

###  4.1 System evolution with no therapy

Starting with experiment (A), we simulate system (50) for the parameters
listed in Table 1. To investigate the behavior of the different cell
populations in the absence of therapy, we preliminary set
$k_{c}=d_{r}=d_{b}=0$. This choice allows us to observe the tumor evolution
over time without external intervention and to study the possible progression
of the disease. The corresponding simulation results are shown in Figure 5,
where the five columns report the evolution of glioma density $M$, endothelial
cell density $W$, VEGF concentration $H$, necrotic matter density $N_{e}$, and
normal tissue density $Q$.

Figure 5: Experiment (A). Numerical simulation of system (50) with the
parameters listed in Table 1 and without any treatment.

Tumor growth and spread, which mainly depend on the choice of the parameters
$\lambda_{0}$ and $s$ (as also observed in [7]), seems to be reasonably fast.
Unlike in [7, 13, 26], here we are not relying on the go-or-grow hypothesis
for the description of tumor evolution, whence we assume a phenotypically
homogeneous tumor. In Figure 5, we observe the growth and spread of the tumor
mass in accordance with the distribution of healthy tissue (shown in the last
column), which is crucially influencing the tumor tactic term, and of the
vasculature. The EC dynamics also influence tumor proliferation through the
function $g(W)$, which here is chosen to describe ECs as enhancing tumor
proliferation and not directly controlling it (as considered, for instance, in
[7]). This choice is mainly related to the fact that tumor cells are still
able to undergo mitosis (at a different rate) when there is a lack of
nutrients; however, the presence of ECs substantially enhances cell
proliferation. In Figure 5, we also observe how ECs grow and diffuse, showing
an evident migration toward the tumor population, due to the produced VEGFs.
The latter is expressed by tumor cells and represents the chemotactic agent
guiding ECs migration and, thus, mediating angiogenesis. In this case, since
no therapy is included, the healthy tissue (last column of Figure 5) is
degraded only by glioma cells activity 333more precisely by the acidity and
matrix-degrading agents they produce, but we wanted to avoid overloading the
model by explicitly accounting for such variables; refer, however, to [7, 34]
for closely related problems with explicit acidity dynamics., and the necrotic
compartment collects this portion of degraded tissue. The tissue degradation
shows an overall behavior comparable to the one observed in [7].

To study the impact of including microscopic dynamics also in the EC equation,
which led to the tactic term in (45), we compare the behavior of system (50)
with that of two similar models where we change the EC description: one of
them simply excludes the microscopic dynamics for ECs, deducing a macroscopic
diffusion equation for these cells, while a second one considers EC migration
as biased by the gradient of tumor cells (instead of the gradient of tumor-
produced VEGF), similar to what is proposed in [7]. In particular, in order to
observe the impact on cell migration, we exclude EC proliferation from the
three compared settings. In the first case, system (50) (without EC
proliferation) is compared with an analogous system where the tumor equation
(44) is coupled with the following reduced equation for ECs:

$\partial_{t}W_{1}-\nabla\cdot\left(\mathbb{D}_{EC}\nabla W_{1}\right)=0$ (63)

We refer to this system as $MC_{1}$ with $\mathbb{D}_{EC}$ given in (46). In
the second case, instead, system (50) (without EC proliferation) is compared
with an analogous one, where the glioma equation (44) is coupled with the EC
equation

$\partial_{t}W_{2}-\nabla\cdot\left(\mathbb{D}_{EC}\nabla
W_{2}\right)+\nabla\cdot\left(\chi_{a}(M_{2})W_{2}\nabla M_{2}\right)=0\,.$
(64)

We refer to this system as $MC_{2}$ and, here,
$\chi_{a}(M_{2}):=\gamma_{0}a(M_{2})\mathbb{D}_{EC}$, where $a(M_{2})$ is the
function accounting for the chemotactic effect of pro-angiogenic signals on
ECs. Following [7], this function can be given as

$a(M_{2}):=\chi_{a,0}\dfrac{K_{M_{2}}}{(K_{M_{2}}+M_{2})^{2}}\,,$

with the parameter $\chi_{a,0}$ describing the duration between EC turnings
damped by large glioma densities. This experiment (B) is depicted in Figure 6,
which collects the results of these comparisons.

Figure 6: Experiment (B). Differences between the respective solution
components of Model (50) and of the analogous one with equation (63) (first
and second column) or equation (64) (third and fourth column) for ECs at 18
weeks (upper row), 24 weeks (middle row), and 30 weeks (lower row).

Precisely, Figure 6 shows the two different comparisons. The first two columns
refer to the comparison between the model given by (50) and the one we named
$MC_{1}$, while the third and fourth columns refer to the comparison between
system (50) and model $MC_{2}$. We plot the differences between the solutions
for the evolution of tumor ($M-M_{i}$) and endothelial cells ($W-W_{i}$)
obtained with the different models. In the first comparison, when the equation
for ECs is the simple diffusion equation (63), we observe - as expected - a
stronger diffusive spread of ECs, which determines the negative values of the
plotted difference (blue regions in the column $W-W_{1}$), while including the
subcellular level for ECs leads to a more directed migration of these cells
towards the tumor location, as shown by the positive values of the plotted
difference (yellow/red regions in the column $W-W_{1}$). Correspondingly, as
ECs enhance tumor proliferation, when their subcellular dynamics are included,
we observe a higher density of tumor cells at the interface between the
regions occupied by ECs and tumor cells (yellow/red regions in the column
$M-M_{1}$). In the second comparison, we observe a much stronger directed
migration of the ECs described by (64). In fact, higher values of EC density
(leading to negative values of the plotted difference in the column $W-W_{2}$)
are shown mainly in the area affected by the tumor mass, whose gradient is
driving the migration, when compared to the case in which EC motion is driven
by the concentration of (tumor-produced) VEGFs. Once these growth factors are
produced, they diffuse in the microenvironment, determining a lower tactic
force on the endothelial cells. Tumor cell growth is, then, correspondingly
affected, as shown in the column $M-M_{2}$. Precisely, there is enhanced tumor
proliferation in the core region in the case described by (64), while there is
more proliferation around the external border in the case described by (50).

###  4.2 Effects of treatment

We study here the effects of the various therapeutic strategies described in
the experiments (C), (D) and (E) on the different system components. We denote
with the index $S_{1}$ the solution components of experiment (C), with $S_{2}$
the solution components of experiment (D), and with $S_{3}$ the ones of
experiment (E). Starting from experiment (C), we first assess the impact of
standard radio- and chemotherapy on the densities of glioma, ECs, necrotic
matter, healthy tissue, and on the concentration of VEGFs. This numerical
experiment is motivated by the fact that in the real case reported in Section
4.3 the patient was treated with this specific combined treatment. Figure 2
shows a scheme of this treatment schedule.

With the same choice of parameters as that employed to generate Figure 5 for
experiment (A), we test the effect of the treatment plan depicted in Figure 2.
The results of this experiment (C) are shown in Figure 7. The first row of
Figure 7 (that also corresponds to the third row of Figure 5) represents the
state of the five species at the beginning of the treatment, while the second
and third rows correspond to the situation after 3 and 6 weeks, respectively.
Precisely, the third row represents the situation at the end of the combined
treatment. In the last row, we show the system evolution after the resting
period of 10 weeks, during which no therapy is applied.

Figure 7: Experiment (C). Numerical simulation of system (50) with the
parameters listed in Table 1 and for the therapeutic plan schematized in
Figure 2.

Comparing Figures 5 and 7, we immediately grasp the effects of the radio- and
chemotherapy on the tumor population, whose density strongly decreases during
treatment, while the density of the necrotic matter increases, as this
component collects the effects of the therapy on tumor, ECs, and healthy
tissue. The reduction in glioma density consequently affects VEGF production,
with a decrease in its expression. The impact of the treatment is also evident
in the evolution of the healthy tissue. In turn, the reduction of normal
tissue affects the proliferation of cancer and endothelial cells. In fact,
this depends on the availability of healthy tissue and, when the latter is
excessively degraded, then proliferation is impaired. This effect can be
observed by comparing the last row of Figures 5 and the third row of Figure 7,
which relates to the system evolution at 30 weeks (at the end of the
treatment).

In experiment (D) we investigate the overall effects of the treatment as given
in experiment (C), but supplemented with anti-angiogenic therapy. Precisely,
the tumor is let to grow and spread for 24 weeks, allowing it to form a
considerable neoplastic mass. Afterward, we apply a therapeutic plan
consisting of a combination of radio- and chemotherapy, the latter involving
temozolomide for inducing apoptosis of glioma cells, and anti-angiogenic
therapy with bevacizumab for 6 weeks. Temozolomide is administered orally
every day at a standard constant dose of $75$ mg/m2. Radiotherapy is applied 5
days per week (from Monday to Friday) for 6 weeks, and the total dose of $60$
Gy is fractionated into smaller doses of $2$ Gy per day. Anti-angiogenic
therapy is administered at a dose of $10$ mg/kg intravenously every 2 weeks,
thus, providing a total of 3 doses in the whole treatment period of 6 weeks.
At the end of the treatment, we let the patient rest for 10 weeks and analyze
how the tumor eventually reorganizes and evolves. Figure 3 shows the specific
scheme of this treatment schedule. We also prolonged the simulations by 10
more therapy-free weeks, in order to better observe a possible tumor relapse.
We consider the same choice of parameters as in Figure 7 for experiment (C) to
test the effect of the described combined therapy plan and we show the
differences in the evolution of the solution components in Figure 8.
Precisely, we consider the differences between the species indicated with
index $S_{1}$ (referring to the schedule in Figure 2) and those indicated with
the index $S_{2}$ (referring to the schedule in Figure 3). Results are shown
at 27 weeks (after three weeks of combined treatment), at 30 weeks (end of the
treatment), at 40 weeks (after 10 weeks of no therapy), and finally at 50
weeks (after 10 more weeks without treatment).

Figure 8: Differences between experiment (C) and experiment (D). Differences
between the respective solution components of system (50) with the two
different treatment schedules sketched in Figures 2 and 3 at 24 weeks (first
row), 27 weeks (second row), 40 weeks (third row), and 50 weeks (last row).

In Figure 8, the effects of the additional anti-angiogenic therapy
administered in the case $S_{2}$ can be observed. Firstly, the EC density is
reduced when bevacizumab is administered, as shown by the positive values of
the plotted difference (yellow/red regions in the column
$W_{S_{1}}-W_{S_{2}}$). The consequent decrease in nutrient availability
provided by the vasculature affects the tumor cell population and impacts the
VEGF dynamics. Its expression is actually enhanced, thus leading to the
negative values of the plotted difference (blue regions in the column
$H_{S_{1}}-H_{S_{2}}$. The situation after 10 weeks of rest (third row of
Figure 8) suggests that in the long term the differences in the EC and VEGF
densities reduce and, for the latter, they become almost negligible. Moreover,
as a consequence of the decreasing EC density, when bevacizumab is
administered the tumor mass is more confined, although featuring a slightly
higher cell density at the tumor core (observed by the emergence of blue areas
in those regions along with a yellow outer rim in the column
$M_{S_{1}}-M_{S_{2}}$). The emergence of the latter induced us to simulate 10
more weeks of tumor evolution. The results are shown in the last row of Figure
8 and they confirm our guess that the difference between the $S_{1}$ and
$S_{2}$ therapy plans becomes more visible after a longer intervention-free
period. This is put in evidence by the clear enhancement of the yellow tumor
rim: more glioma cells are spreading into the tumor-adjacent area if no anti-
angiogenic treatment has been applied. Especially for glioma, this might be
relevant, since very often the (infiltrative) margins of the treated tumor are
the source for relapse. As with every therapy, one should also take into
account the side effects, which should be balanced against reducing the tumor
spread. The last column of the figure suggests that the effect of the anti-
angiogenic drug on normal tissue is negligible. Concerning the necrotic
matter, it correspondingly infers a substantial decrease when the treatment
plan of Figure 3 is administered, as observed in the column
$N_{S_{1}}-N_{S_{2}}$. These simulations suggest that scheme $S_{2}$ might
have some advantage against $S_{1}$, in line with some clinical results (refer
e.g., to trial NCT01390948) showing a slight increase in the 1-year
progression-free survival (from 67% to 74%) after the end of the treatment
combining radio-, chemo-, and anti-angiogenic therapy. However, the clinical
debate about the effects of bevacizumab is still open, as it did not seem to
improve considerably the overall therapy outcome, despite its clear impact
during the period of administration.

Finally, in experiment (E), we analyze the effects of the treatment schedule
sketched in Figure 4. Precisely, after letting the tumor grow for 24 weeks, we
apply a combined treatment of radio- and chemotherapy. The former is applied 5
days per week (from Monday to Friday) for 6 weeks, at a fractionated dose of
$2$ Gy per day (total dose of $60$ Gy), while the latter is based on
temozolomide administration orally every day at a standard constant dose of
$75$ mg/m2. Then, after a resting period of 4 weeks, adjuvant anti-angiogenic
therapy is applied at a standard dose of $10$ mg/kg intravenously every 2
weeks for other 6 weeks, thus providing a total of 3 doses during the whole
treatment period. With the same choice of parameters as in Figure 7 for
experiment (C), we test the effects of the described combined therapy against
the alternative therapy plan proposed in experiment (D). We show the
differences in the evolution of the solution components for the two schedules
in Figure 9. Precisely, we consider the differences between the populations
indicated with index $S_{2}$ (referring to the schedule in Figure 3) and the
one indicated with the index $S_{3}$ (referring to the schedule in Figure 4).
Results are shown at 27 weeks (after three weeks of radio- and chemotherapy),
at 30 weeks (end of this treatment), at 40 weeks (after the resting period for
the case $S_{2}$ or after the adjuvant anti-angiogenic treatment for the case
$S_{3}$), and at 50 weeks, allowing for 10 more weeks without therapy for the
follow-up.

Figure 9: Differences between experiment (D) and experiment (E). Differences
between the respective solution components of system (50) with the two
different treatment schedules sketched in Figures 3 and 4 at 24 weeks (first
row), 27 weeks (second row), 40 weeks (third row), and 50 weeks (last row).

From Figure 9 we can immediately notice the effects of the different schedules
involving anti-angiogenic therapy on the EC population. Precisely, a
concurrent administration of radio-, chemo-, and anti-angiogenic therapy seems
to have a stronger effect on the EC density, which appears more reduced in
case $S_{2}$ than in case $S_{3}$, as shown by the negative values of the
plotted difference (blue regions in the column $W_{S_{2}}-W_{S_{3}}$). This
has a clear impact on the tumor cell population, whose spread is reduced.
Instead, the adjuvant anti-angiogenic therapy given after the resting period
of 4 weeks seems to allow a stronger tumor infiltration into the surrounding
areas, as shown by the blue outer ring of tumor cells in the column
$M_{S_{2}}-M_{S_{3}}$. The reduction in EC density for the case $S_{2}$
correspondingly affects the VEGF component (yellow/red regions in the column
$H_{S_{2}}-H_{S_{3}}$), which infers less uptake by ECs. The necrotic
component sees a substantial increase in a much extended region when the
treatment in Figure 4 is administered, as tumor and EC activities are less
controlled. The overall effect on the healthy tissue is, instead, minimal, in
line with the results shown in Figure 8.

###  4.3 Application to real glioma patient data

We perform further numerical simulations starting from real data of a
glioblastoma patient in order to assess the predicted numerical outcome of the
proposed therapeutic schemes.

The patient aged 75 years first presented at the Department of Neurosurgery
(Saarland University Medical Centre) with word-finding disorders, unsteady
gait, and disorientation regarding time and situation. Computed tomography
scans and subsequent magnetic resonance imaging revealed a left temporo-
occipital brain lesion suspect for glioblastoma, which was confirmed
pathologically after gross total resection. The histological work-up revealed
glioblastoma cells, necrotic streaks, and prominent, partially hyperplastic
endothelial cells. A combined radio- and chemotherapy was recommended.
Precisely, the treatment followed the standard protocol, applying 60 Gy in 30
fractions of 2 Gy over 6 weeks with simultaneous daily temozolomide (75 mg/m2
surface area), analogous to the scheme shown in Figure 2. Due to severe
thrombopenia, chemotherapy was stopped two days before the end of radiotherapy
and was not resumed as adjuvant therapy.

We simulate a non-dimensionalized version of the macroscopic setting (50) on a
horizontal brain slice reconstructed from the processing of the patient-
specific MRI scans. The macroscopic tensor $\mathbb{D}_{T}(x)$ is pre-
calculated using the patient DTI data and the bimodal von Mises-Fisher
distribution [39] for the fiber distribution function, whose expression is
given by:

$q(x,\theta)=\dfrac{\delta}{2\pi}+(1-\delta)\dfrac{1}{4\pi
I_{0}(k(x))}(e^{k(x)v_{1}\cdot\theta}+e^{-k(x)v_{1}\cdot\theta}).$

Here, the parameter $\delta\in[0,1]$ describes the inherent degree of
randomized turning of the cell, while the function $k(x):=\kappa
FA(\mathbb{D}_{W}(x))$, refers to the level of concentration around the
dominant direction, with the parameter $\kappa$ quantifying the sensitivity of
cells to the directional information given by the environment. $I_{0}(k(x))$
is the modified Bessel function of first kind and of order $0$ and $v_{1}$ is
the leading eigenvector of $\mathbb{D}_{W}$. In particular, the parameters
$\delta$ and $\kappa$ are optimized using the patient data as in [6], while
the values for the other parameters are listed in Table 1.

We reconstruct the initial condition for tumor cells and necrotic matter from
the provided patient segmentation data, where the different components of the
tissue are classified into cerebrospinal fluid, gray matter, white matter,
necrotic matter, edema, nonenhancing and enhancing tumor. In particular, for
the tumor population we consider the areas occupied by both enhancing and
nonenhancing tumor. For the healthy tissue, we used the formulation given in
(59), with $l_{c}(x,y)$ and $h$ estimated from the patient DTI. For the
initial density of ECs, which is neither explicitly available from the DTI
brain data set nor the segmentation, we follow the choice proposed in [12].
Precisely, we choose the random variable

$W_{0}(x,y)=(10^{-3}+u_{r})Q_{0}(x,y),$

where $u_{r}\sim U(0,10^{-3})$ denotes a random array derived from a
continuous uniform distribution on $[0,10^{-3}]$. Finally, for the VEGF
initial profile, we assume a low concentration located in the region occupied
by necrotic matter, where the lack of nutrients would favor VEGF production,
i.e.,

$H_{0}(x,y)=0.1\,N_{e,0}\,.$

Figure 10 shows the initial conditions on the entire 2D brain slice.

Figure 10: Initial conditions of system (50) reconstructed from real patient
data. The five columns refer to densities of glioma, ECs, VEGFs, necrotic
matter, and healthy tissue, respectively.

We can notice that there are some errors or missed data in the border regions
of the initial distributions derived from the DTI data (i.e., $Q_{0}$ and
$W_{0}$). In fact, the main problem is that DTI measurements are usually
recorded to analyze diffusion in white matter and, thus, data acquisition is
terminated as soon as it reaches margins of such regions. As a consequence,
some values in these border areas can be missed and set to zero. However, the
region of greater interest for our study, i.e., the one affected by the tumor
mass, does not show particular problems to deal with. Concerning the treatment
schedule, we consider the one described in Figure 2, i.e., a combination of
chemo- and radiotherapy, as this was the plan actually administered to the
patient. Precisely, we reconstructed from the data the dose distribution with
which the patient was actually irradiated and we modified the terms (54)
modeling radiotherapy assuming a spatial-dependent radiotherapy dose
$d_{r}(x,y)$. Figure 11 shows the isodose curves corresponding to radiation
fractions of 10, 20, 30, 40, 50, and 60 Gy.

Figure 11: Isodose curves of the radiotherapy treatment. Level set curves
referring to the specific dose distribution $d_{r}(x,y)$ with which the
patient was irradiated.

Starting from the initial conditions given in Figure 10, we simulate the
system evolution for 10 weeks (6 weeks of combined therapy + 4 weeks of rest).
Results are given in Figure 12.

Figure 12: Experiment (C) on real patient data. Numerical simulation of system
(50) with the parameters listed in Table 1, for the therapeutic plan
schematized in Figure 2, and on the basis of the provided patient data.

Firstly, we observe the tumor dynamics (first column of Figure 12), which is
characterized by a strong reduction of the tumor density in the area
originally occupied by the tumor mass. This shows the efficacy of the combined
treatment in reducing the overall tumor volume. Moreover, we notice that, even
after the resting period of 4 weeks, no evident regrowth of the tumor mass
emerges. Concerning VEGFs and ECs, we observe that the production of VEGFs by
tumor cells and their accumulation in the region originally affected by the
tumor determines a correspondingly increasing concentration of ECs in the same
area. In turn, the increasing density of ECs, together with the reduction of
tumor density, leads to a noticeable decrease of VEGF concentration, as shown
in the third column of Figure 12. In particular, the EC density (especially
after three weeks of therapy) shows a certain level of heterogeneity, which
may be due to the random distribution we chose as initial condition for this
population and which determines a corresponding heterogeneity in the VEGFs
concentration (first plot in the third column of Figure 12). An outer rim of
VEGFs seems to remain after the 10 weeks and it may be due to a corresponding
low-dense outer rim of tumor cells surviving the treatment. Looking at the
evolution of necrotic matter and healthy tissue, we notice the effects of
having assumed a space-dependent radiotherapy dose distribution. In fact, the
degradation of the healthy tissue is not homogenous and the further the region
from the tumor area, the smaller the effect of the radiations on healthy
tissue. In fact, the most distant regions (e.g. those in the right hemisphere)
do not seem to be affected by the therapy. The effects of both heterogeneous
tissue degradation and tumor cell death are, then, collected in the necrotic
compartment.

##  5 Discussion

The model proposed in this note is a therapy-oriented development of that
considered in [7] and aligns to the approaches in [6, 8, 9, 13, 14, 15, 26].
Unlike [7, 15] we do not differentiate here between moving and proliferating
cancer cells, but account for a single population of cells forming the tumor.
In fact, the analysis in [7] showed (at least in our framework) that the main
effect of splitting the tumor population into moving and resting (and
proliferating) cells was the slower evolution of tumor mass, but qualitatively
it behaved similarly to the situation with one population of moving and
proliferating cancer cells. Therefore, in this note we do not differentiate
between the two cell phenotypes and focus instead on therapy effects. We
consider tumor cells interacting with their environment, whose main physical
components are the anisotropic brain tissue as well as brain-own and tumor-
generated vasculature. We also modeled explicitly the evolution of VEGF
controlling angiogenesis and the necrotic component of the neoplasm, which is
relevant for assessing the tumor stage and for the segmentation needed for
treatment planning. The model investigates the impact of including microscopic
dynamics also in the equation for endothelial cells, which has not been done
before, comparing our approach with other two possible descriptions of EC
dynamics. Then, we analyze via simulations various therapy approaches and
their effects on tumor development and on normal tissue, thereby following
treatment schedules addressed in clinical studies available in the literature.
Although this has been done in [26] in a rather rudimentary way, we do not
include tumor resection as part of the therapy; the problems have already been
stated in that reference: it is hardly possible to provide a proper
mathematical characterization of the tissue and neoplasm dynamics in the
resected region in a continuous manner and the solution infers sharp
discontinuities, which are difficult to handle. Considering the real data
employed in Subsection 4.3, the simulations showed qualitatively reasonable
results, however, a quantitative assessment could not be performed, due to
missing data, as that patient’s therapy was stopped.

System (50) with its initial and boundary conditions is also interesting from
a mathematical viewpoint: it describes a novel sort of taxis cascade, in which
tumor cells are chemotactically following ECs, which in turn bias their motion
towards gradients of a chemical signal (VEGF) produced by tumor cells and
depleted by themselves. The global well-posedness of a system featuring this
complexity has not been investigated and raises manifold challenges, due to
the intricate couplings and nonlinearities involved in the source and motility
terms. The passage from KTEs to the macroscopic dynamics has been done here in
a merely formal way, the rigorous convergence is still open and might probably
use some of the ideas in [55], where a much simpler system has been
considered.

In [31] we commented about the feasibility of using such multiscale models to
predict tumor spread and establish CTV and PTV margins for treatment planning.
Those observations still apply here - the main issue remains the relatively
large number of parameters and therewith related uncertainties. However, the
increasingly fast development of biomedical imaging, computing power, and
technology for necessary cell biology experiments will provide a means to
assess at least some of the missing quantitative information. On the other
hand, such multiscale models seem to offer an adequate frame for studying the
effects of various therapy ansatzes. In fact, chemo- and radiotherapy
primarily act on the level of single cells and ultimately lead to the observed
effects on the whole tumor and this mathematical approach allows us to account
for dynamics on both levels in a reasonably detailed manner. Our numerical
experiments also suggest that letting the clinical studies have a longer
follow-up might provide useful information about the tumor behavior after
ceasing the actual therapy. That would presumably lead to higher costs, but
these might be justified by the achieved understanding. Mathematical models
could help in identifying the adequate duration of such studies. Moreover,
they have the potential to investigate a great variety of therapeutic
scenarios (of which we showed here just a few examples) in an unprecedented
complexity and accuracy - provided the necessary quantitative information
becomes available. Intra- and interdisciplinary studies including such models
are called upon to shed light on the intricate biological processes associated
with tumor growth, expansion, and treatment response.

## Acknowledgments

MC acknowledges funding by the by Ministry of Education, Universities and
Research, through the MIUR grant Dipartimento di Eccellenza 2018-2022, Project
no. E11G18000350001. This work was also supported by National Group of
Mathematical Physics” (GNFM-INdAM). MC acknowledges also funding by the Basque
Government through the BERC 2018- 2021 program, by the Spanish State Research
Agency through BCAM Severo Ochoa excellence accreditation SEV-2017-0718 and by
the MINECO- Feder (Spain) through the research grant number
RTI2018-098850-B-I00. MC has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No. 713673. The project that gave rise to these results
received the support of a fellowship from “a Caixa” Foundation (ID 100010434).
The fellowship code is LCF/BQ/IN17/11620056. CS was partially supported by the
Federal Ministry of Education and Research BMBF, project GlioMaTh 05M2016.

##  Appendix A Non-dimensionalization

To proceed with the non-dimensionalization, we firstly observe that the
variables $M$, $W$, $N_{e}$ involved in systems (50) are expressed in
cells/mm3, $Q$ in mg/mm3, while the concentration of VEGF $H$ is given in
mol/liter (=:M). The reference values we use for the adimensionalization are
listed in Table 2. In particular, we rescale the tumor, EC, and necrotic
matter densities with respect to their carrying capacities, i.e.,
$M_{c,0}=K_{N}$, $W_{c,0}=K_{W}$, and $N_{e,0}=K_{Ne}$.

Parameter | Description | Value (units) | Source
---|---|---|---
T | time |  $1$ (d) |
L | length |  $0.875$ (mm) |
$M_{c,0}$ | tumor cell density |  $10^{5}$ (cell$\cdot$mm-3) | [7]
$W_{c,0}$ | ECs density |  $10^{5}$ (cell$\cdot$mm-3) | [7]
$N_{e,0}$ | Necrotic matter density |  $10^{5}$ (cell$\cdot$mm-3) | [7]
$H_{c,0}$ | VEGF concentration |  $10^{-9}$ (M) | [53]
$Q^{*}$ | healthy tissue density |  $10^{-3}$ (mg$\cdot$mm-3) | [7]
Table 2: Reference variables

We nondimensionalize the partial differential equations introduced above as
follows:

$\tilde{t}=\dfrac{t}{T}\quad\quad\tilde{x}=\dfrac{x}{L}\quad\quad\tilde{M}=\dfrac{M}{M_{c,0}}\quad\quad\tilde{W}=\dfrac{W}{W_{c,0}}\quad\quad\tilde{N_{e}}=\dfrac{N_{e}}{N_{e,0}}\quad\quad\tilde{H}=\dfrac{H}{H_{c,0}}\quad\quad\tilde{Q}=\dfrac{Q}{Q^{*}}\,.$

A proper scaling of the parameters involved in the macroscopic setting leads
to:

$\begin{split}&\tilde{k}_{W}^{+}=\dfrac{k_{W}^{+}}{\lambda_{0}}\quad\quad\tilde{k}_{Q}^{+}=\dfrac{k_{Q}^{+}}{\lambda_{0}}\quad\quad\tilde{k}_{1}^{-}=\dfrac{k_{1}^{-}}{\lambda_{0}}\quad\quad\tilde{k}_{H}^{+}=\dfrac{\bar{k}_{H}^{+}}{\gamma_{0}}\quad\quad\tilde{\mu}_{M,0}=\mu_{M,0}\,T\\\\[8.5359pt]
&\tilde{\lambda_{1}}=\dfrac{\lambda_{1}}{k_{1}^{-}}\quad\quad\tilde{k}_{H}^{-}=\dfrac{\bar{k}_{H}^{-}}{\gamma_{0}}\quad\quad\tilde{\gamma}_{1}=\dfrac{\gamma_{1}}{\gamma_{0}}\quad\quad\tilde{\mu}_{W,0}=\bar{\mu}_{W,0}\,T\\\\[8.5359pt]
&\tilde{\mathbb{D}}_{T}=\dfrac{1}{\lambda_{0}}\dfrac{T}{L^{2}}\mathbb{D}_{T}\quad\quad\tilde{\mathbb{D}}_{EC}(x)=\dfrac{T}{L^{2}}\,\mathbb{D}_{EC}(x)\quad\quad\tilde{D}_{H}=\dfrac{T}{L^{2}}D_{H}\quad\quad\tilde{c}_{1}=\dfrac{c_{1}}{H_{c,0}}T\quad\quad\tilde{c}_{2}=c_{2}\,T.\end{split}$

If we drop the tilde (” $\tilde{}$ ”) in the new variables and parameters,
then the differential equations in sistem (50) keep the same form, with the
properly rescaled expressions for following functions:

$\begin{split}&\tilde{\mu}(\tilde{M},\tilde{W},\tilde{N}_{e}):=\tilde{\mu}_{M,0}\,\left(1-\tilde{M}-\tilde{N}_{e}\right)g(\tilde{W})\,,\qquad\tilde{\mu}_{W,0}(d_{b})=\tilde{\mu}_{W,0}\left(1+\dfrac{1}{1+d_{b}^{2}}\right)\,,\\\\[8.5359pt]
&\tilde{\mu}_{W}(\tilde{W},\tilde{Q},\tilde{d_{b}})=\tilde{\mu}_{W,0}(d_{b})(1-\tilde{W})\tilde{Q}\,,\qquad\tilde{B}_{p}(\tilde{W},\tilde{Q})=\left(\tilde{k}_{W}^{+}S_{W}\tilde{W}+\tilde{k}_{Q}^{+}S_{Q}\tilde{Q}+\tilde{k}_{1}^{-}\right)\,,\qquad\tilde{B}_{w}(\tilde{H})=(\tilde{k}_{H}^{+}\tilde{H}+\tilde{k}_{H}^{-}).\\\
\end{split}$

## References

  * [1] Estimation taken from.  https://bionumbers.hms.harvard.edu/bionumber.aspx?s=n&v=0&id=108941.
  * [2] Estimation taken from.  http://www.lab.anhb.uwa.edu.au/mb140/MoreAbout/Endothel.htm.
  * [3] L. Barazzuol, N.G. Burnet, R. Jena, B. Jones, S.J. Jefferies, and N.F. Kirkby.  A mathematical model of brain tumour response to radiotherapy and chemotherapy considering radiobiological aspects.  Journal of theoretical biology, 262(3):553–565, 2010.
  * [4] N. Bellomo.  Modeling Complex Living Systems.  Birkhäuser Boston, 2008.
  * [5] M.E. Berens and A. Giese.  “…those left behind.” Biology and oncology of invasive glioma cells.  Neoplasia, 1(3):208–219, 1999.
  * [6] M. Conte, L. Gerardo-Giorda, and M. Groppi.  Glioma invasion and its interplay with nervous tissue and therapy: A multiscale model.  Journal of Theoretical Biology, 486:110088, 2020.
  * [7] M. Conte and C. Surulescu.  Mathematical modeling of glioma invasion: acid- and vasculature mediated go-or-grow dichotomy and the influence of tissue anisotropy.  Applied Mathematics and Computation, 407:126305, 2021.
  * [8] G. Corbin, C. Engwer, A. Klar, J. Nieto, J. Soler, C. Surulescu, and M. Wenske.  Modeling glioma invasion with anisotropy- and hypoxia-triggered motility enhancement: From subcellular dynamics to macroscopic pdes with multiple taxis.  Mathematical Models and Methods in Applied Sciences, 31(01):177–222, 2021.
  * [9] G. Corbin, A. Hunt, A. Klar, F. Schneider, and C. Surulescu.  Higher-order models for glioma invasion: From a two-scale description to effective equations for mass density and momentum.  Mathematical Models and Methods in Applied Sciences, 28(09):1771–1800, 2018.
  * [10] A. Czirok.  Endothelial cell motility, coordination and pattern formation during vasculogenesis.  Wiley Interdisciplinary Reviews: Systems Biology and Medicine, 5(5):587–602, 2013.
  * [11] W. Diao, X. Tong, C. Yang, F. Zhang, C. Bao, H. Chen, L. Liu, M. Li, F. Ye, Q. Fan, and J. Wang.  Behaviors of glioblastoma cells in in vitro microenvironments.  Scientific reports, 9(1):1–9, 2019.
  * [12] A. Dietrich, N. Kolbe, N. Sfakianakis, and C. Surulescu.  Multiscale modeling of glioma invasion: from receptor binding to flux-limited macroscopic PDEs.  Multiscale Modeling and Simulation, accepted, arXiv:2010.03277.
  * [13] C. Engwer, T. Hillen, M. Knappitsch, and C. Surulescu.  Glioma follow white matter tracts: a multiscale DTI-based model.  Journal of Mathematical Biology, 71(3):551–582, 2014.
  * [14] C. Engwer, A. Hunt, and C. Surulescu.  Effective equations for anisotropic glioma spread with proliferation: a multiscale approach.  IMA Journal of Mathematical Medicine and Biology, 33:435–459, 2016\.
  * [15] C. Engwer, M. Knappitsch, and C. Surulescu.  A multiscale model for glioma spread including cell-tissue interactions and proliferation.  Mathematical Biosciences and Engineering, 13(2):443–460, 2016.
  * [16] C. Engwer, C. Stinner, and C. Surulescu.  On a structured multiscale model for acid-mediated tumor invasion: The effects of adhesion and proliferation.  Mathematical Models and Methods in Applied Sciences, 27:1355–1390, 2017.
  * [17] V. Pereira Ferrer, N. Vivaldo Moura, and R. Mentlein.  Glioma infiltration and extracellular matrix: key players and modulators.  Glia, 66(8):1542–1565, 2018.
  * [18] J.F. Fowler.  The linear-quadratic formula and progress in fractionated radiotherapy.  The British journal of radiology, 62(740):679–694, 1989.
  * [19] H.S. Friedman, M.D. Prados, P.Y. Wen, T. Mikkelsen, D. Schiff, L.E. Abrey, W.K.A. Yung, N. Paleologos, M.K. Nicholas, R. Jensen, et al.  Bevacizumab alone and in combination with irinotecan in recurrent glioblastoma.  Journal of clinical oncology, 27(28):4733–4740, 2009.
  * [20] F. Mac Gabhann and A.S. Popel.  Targeting neuropilin-1 to inhibit vegf signaling in cancer: comparison of therapeutic approaches.  PLoS computational biology, 2(12):e180, 2006.
  * [21] J.L. Gevertz and S. Torquato.  Modeling the effects of vasculature evolution on early brain tumor growth.  Journal of Theoretical Biology, 243(4):517–531, 2006.
  * [22] A. Giese, R. Bjerkvig, M.E. Berens, and M. Westphal.  Cost of migration: Invasion of malignant gliomas and implications for treatment.  Journal of Clinical Oncology, 21(8):1624–1636, 2003.
  * [23] A. Giese, L. Kluwe, H. Meissner, E. Michael, and M. Westphal.  Migration of human glioma cells on myelin.  Neurosurgery, 38:755–764, 1996.
  * [24] D. Hanahan and R.A. Weinberg.  Hallmarks of cancer: the next generation.  Cell, 144(5):646–674, 2011.
  * [25] C.H. Holdsworth, D. Corwin, R.D. Stewart, R. Rockne, A.D. Trister, K.R. Swanson, and M. Phillips.  Adaptive IMRT using a multiobjective evolutionary algorithm integrated with a diffusion–invasion model of glioblastoma.  Physics in Medicine and Biology, 57(24):8271–8283, 2012.
  * [26] A. Hunt and C. Surulescu.  A multiscale modeling approach to glioma invasion with therapy.  Vietnam Journal of Mathematics, 45(1-2):221–240, 2016.
  * [27] S. Jbabdi, E. Mandonnet, H. Duffau, L. Capelle, K.R. Swanson, M. Pélégrini-Issac, R. Guillevin, and H. Benali.  Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging.  Magnetic Resonance in Medicine, 54(3):616–624, 2005.
  * [28] M.C. Joiner and A. Van der Kogel.  Basic clinical radiobiology fourth edition.  CRC press, 2009.
  * [29] L.D. Ke, Y.X. Shi, S.A. Im, X. Chen, and W.K.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(6):2562–2572, 2000.
  * [30] J. Kelkel and C. Surulescu.  A multiscale approach to cell migration in tissue networks.  Mathematical Models and Methods in Applied Sciences, 22(03):1150017, 2012.
  * [31] S. Knobe, Y. Dzierma, M. Wenske, C. Berdel, J. Fleckenstein, P. Melchior, J. Palm, F. Nuesken, A. Hunt, C. Engwer, et al.  Feasibility and clinical usefulness of modelling glioblastoma migration in adjuvant radiotherapy.  Zeitschrift für Medizinische Physik, 2021.
  * [32] J.M. Kroos, C. Stinner, C. Surulescu, and N. Surulescu.  SDE-driven modeling of phenotypically heterogeneous tumors: The influence of cancer cell stemness.  Discrete and Continuous Dynamical Systems - Series B, 24(8):4629, 2019.
  * [33] P. Kumar, J. Li, and C. Surulescu.  Multiscale modeling of glioma pseudopalisades: contributions from the tumor microenvironment.  Journal of Mathematical Biology, 82(6), 2021.
  * [34] P. Kumar and C. Surulescu.  A flux-limited model for glioma patterning with hypoxia-induced angiogenesis.  Symmetry, 12(11):1870, 2020.
  * [35] D.A. Lauffenburger and J.L. Lindermann.  Receptors. Models for binding, trafficing and signaling. Oxford University Press, 1993.
  * [36] H.A. Levine, S. Pamuk, B.D. Sleeman, and M. Nilsen-Hamilton.  Mathematical modeling of capillary formation and development in tumor angiogenesis: Penetration into the stroma.  Bulletin of Mathematical Biology, 63(5):801–863, 2001.
  * [37] T. Lorenz and C. Surulescu.  On a class of multiscale cancer cell migration models: Well-posedness in less regular function spaces.  Mathematical Models and Methods in Applied Sciences, 24(12):2383–2436, 2014.
  * [38] M. Onishi, T. Ichikawa, K. Kurozumi, and I. Date.  Angiogenesis and invasion in glioma.  Brain tumor pathology, 28(1):13–24, 2011.
  * [39] K. Painter and T. Hillen.  Mathematical modelling of glioma growth: the use of diffusion tensor imaging (DTI) data to predict the anisotropic pathways of cancer invasion.  Journal of Theoretical Biology, 323:25–39, 2013.
  * [40] G. Powathil, M. Kohandel, S. Sivaloganathan, A. Oza, and M. Milosevic.  Mathematical modeling of brain tumors: effects of radiotherapy and chemotherapy.  Physics in Medicine & Biology, 52(11):3291, 2007.
  * [41] X. S. Qi, C.J. Schultz, and X.A. Li.  An estimation of radiobiologic parameters from clinical outcomes for radiation treatment planning of brain tumor.  International Journal of Radiation Oncology, Biology, Physics, 64(5):1570–1580, 2006.
  * [42] D.A. Reardon, S. Turner, K.B. Peters, A. Desjardins, S. Gururangan, J.H. Sampson, R.E. McLendon, J.E. Herndon, L.W. Jones, J.P. Kirkpatrick, et al.  A review of vegf/vegfr-targeted therapeutics for recurrent glioblastoma.  Journal of the National Comprehensive Cancer Network, 9(4):414–427, 2011.
  * [43] R. Rockne, J.K. Rockhill, M. Mrugala, A.M. Spence, I. Kalet, K. Hendrickson, A. Lai, T. Cloughesy, E.C. Alvord, and K.R. Swanson.  Predicting the efficacy of radiotherapy in individual glioblastoma patientsin vivo:a mathematical modeling approach.  Physics in Medicine and Biology, 55(12):3271–3285, 2010.
  * [44] T. Seyfried.  Cancer as a metabolic disease: on the origin, management, and prevention of cancer.  John Wiley & Sons, 2012.
  * [45] W.R. Shapiro.  Current Therapy for Brain Tumors: Back to the Future.  Archives of Neurology, 56(4):429–432, 1999.
  * [46] M. Sidani, D. Wessels, G. Mouneimne, M. Ghosh, S. Goswami, C. Sarmiento, W. Wang, S. Kuhl, M. El-Sibai, J-M. Backer, R. Eddy, D. Soll, and J. Condeelis.  Cofilin determines the migration behavior and turning frequency of metastatic cancer cells.  Journal of Cell Biology, 179(4):777–791, 2007.
  * [47] W. Stummer.  Extent of resection and survival in glioblastoma multiforme.  Neurosurgery, 64(6):E1206, 2009.
  * [48] R. Stupp, M.E. Hegi, W.P. Mason, M.J. Van Den Bent, M.J.B. Taphoorn, R.C. Janzer, S.K. Ludwin, A. Allgeier, B. Fisher, K. Belanger, et al.  Effects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase iii study: 5-year analysis of the eortc-ncic trial.  The lancet oncology, 10(5):459–466, 2009.
  * [49] R. Stupp, W.P. Mason, M.J. Van Den Bent, M. Weller, B. Fisher, M.J.B. Taphoorn, K. Belanger, A.A. Brandes, C. Marosi, U. Bogdahn, et al.  Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma.  New England journal of medicine, 352(10):987–996, 2005.
  * [50] A. Swan, T. Hillen, J.C. Bowman, and A.D. Murtha.  A patient-specific anisotropic diffusion model for brain tumour spread.  Bulletin of Mathematical Biology, 80(5):1259–1291, 2017.
  * [51] K.R. Swanson, R.C. Rockne, J. Claridge, M.A. Chaplain, E.C. Alvord, and A.R.A. Anderson.  Quantifying the role of angiogenesis in malignant progression of gliomas: In silico modeling integrates imaging and histology.  Cancer Research, 71(24):7366–7375, 2011.
  * [52] A. Szabó, R. Unnep, E. Méhes, W.O. Twal, W.S. Argraves, Y. Cao, and A. Czirók.  Collective cell motion in endothelial monolayers.  Physical biology, 7(4):046007, 2010.
  * [53] S. Takano, Y. Yoshii, S. Kondo, H. Suzuki, T. Maruno, S. Shirai, and T. Nose.  Concentration of vascular endothelial growth factor in the serum and tumor tissue of brain tumor patients.  Cancer research, 56(9):2185–2190, 1996.
  * [54] P. Vempati, A.S. Popel, and F. Mac Gabhann.  Extracellular regulation of vegf: isoforms, proteolysis, and vascular patterning.  Cytokine & growth factor reviews, 25(1):1–19, 2014.
  * [55] A. Zhigun and C. Surulescu.  A novel derivation of rigorous macroscopic limits from a micro-meso description of signal-triggered cell migration in fibrous environments.  SIAM Journal on Applied Mathematics, 82(1):142–167, 2022.