# Spatio-temporal modelling of phenotypic heterogeneity in tumour tissues and
its impact on radiotherapy treatment.

Giulia L. Celora  Helen M. Byrne  P.G. Kevrekidis  Mathematical Institute,
University of Oxford, Oxford, UK  Department of Mathematics & Statistics,
University of Massachusetts, Amherst 01003 USA 11footnotemark: 1

###### Abstract

We present a mathematical model that describes how tumour heterogeneity
evolves in a tissue slice that is oxygenated by a single blood vessel.
Phenotype is identified with the stemness level of a cell, $s$, that
determines its proliferative capacity, apoptosis propensity and response to
treatment. Our study is based on numerical bifurcation analysis and dynamical
simulations of a system of coupled non-local (in phenotypic “space”) partial
differential equations that links the phenotypic evolution of the tumour cells
to local oxygen levels in the tissue. In our formulation, we consider a 1D
geometry where oxygen is supplied by a blood vessel located on the domain
boundary and consumed by the tumour cells as it diffuses through the tissue.
For biologically relevant parameter values, the system exhibits multiple
steady states; in particular, depending on the initial conditions, the tumour
is either eliminated (“tumour-extinction”) or it persists (“tumour-invasion”).
We conclude by using the model to investigate tumour responses to radiotherapy
(RT), and focus on establishing which RT strategies can eliminate the tumour.
Numerical simulations reveal how phenotypic heterogeneity evolves during
treatment and highlight the critical role of tissue oxygen levels on the
efficacy of radiation protocols that are commonly used clinically.

###### keywords:

cancer stem cells , radio-resistance , heterogeneity

††journal:  label4label4footnotetext<EMAIL_ADDRESS>
##  1 Introduction

In recent years, there has been a paradigm shift in how we understand cancer.
New technologies, such as single-cell sequencing, have showcased the complex
composition of tumours where cells with different genotypes or phenotypes
coexist in the same tissue region [11, 48, 56]. While cancer was initially
defined as a genetic disease, recent evidence suggests that phenotypic
heterogeneity arises from both genetic and non-genetic sources of variability,
such as epigenetic and environmental factors [11, 13, 49]. Viewed through this
lens, tumours are complex ecosystems in which competition and cooperation
between cells, combined with environmental pressures, shape the system
dynamics [33, 49] and its response to treatment. Poor treatment outcomes have
indeed been related to tumour heterogeneity [18], with failure of therapies
being associated with the presence and emergence of resistant phenotypes [8,
49]. Resistant cells are indeed selected for during treatment and can drive
subsequent relapse. Understanding the forces that drive such heterogeneity is
therefore important for the design of more effective treatment strategies [8,
13, 49].

The concept of cancer stem cells (CSCs) was originally introduced to describe
the hierarchical organization of tumours, motivated by analogies with normal
tissue development [1, 33]. While differentiated cells (DCs) have a finite
clonogenic capacity, (_i.e._ , after dividing a finite number of times they
exit the cell-cycle and become terminally differentiated (TDCs)), CSCs can
give rise to cells with the same clonogenic capacity when they divide. This
capacity for self-renewal makes CSCs key drivers of tumour initiation,
progression and relapse. The coexistence of CSCs, DCs and TDCs in tumours is a
source of intra-tumour heterogeneity. Furthermore, CSCs can contribute to
heterogeneity because they can adapt to stressful environmental conditions
(such as nutrient starvation or treatment) by changing their metabolism and/or
entering a quiescent state [1, 24, 66]. From this point of view, CSCs have the
ability to survive standard treatment by transiently changing their behaviour
[53] and can transmit these changes to their differentiated progeny resulting
in a differentiated tumour bulk that is also resistant to treatment.

Hypoxia, i.e., abnormally low oxygen levels, is commonly found in solid
tumours due to aberrant vascularisation [38, 57]. As oxygen diffuses into a
tumour and is consumed by cells, spatial gradients in oxygen levels are
established and regions of chronic or transient hypoxia may emerge. The
clinical implications of hypoxia and its link to several hallmarks of cancer
[36, 38, 60], such evasion of apoptosis, genetic instabilities, induction of
angiogenesis and metastasis, have long been known. The response of cells to
hypoxia is in part regulated by _hypoxia inducible factors_ (HIFs) [60], which
can influence cancer cell metabolism, proliferation and “stemness” (a term
that we use hereafter when referring to phenotypic variation of cancer cells).
It has been shown that hypoxia facilitates the formation of stem-like cells by
up-regulating stem-related genes in cancer cells in a HIF-dependent manner [6,
42, 37]. _De novo_ formation of CSCs can therefore occur in tumours in
response to hypoxia. Indeed, hypoxia-mediated stress can drive cells to
acquire a more aggressive, stem-like phenotype [59].

Mathematical models are a useful tool for studying the complex interplay
between spatio-temporal variation in the tumour micro-environment and
phenotypic heterogeneity and how this interplay impacts a tumours’ growth
dynamics and response to treatment. The critical role of space in shaping
tumour composition and/or treatment outcome has been studied in previous
works, where a wide range of mathematical formalisms has been used: from
discrete agent based models [7, 12, 19, 54, 62, 68] to continuum models of
different degree of complexity [5, 23, 29, 34, 44, 75]. In particular,
phenotypic-structured models have been shown to be a useful tool for
understanding intra-tumour phenotypic heterogeneity [2, 16, 17, 25, 39, 45,
47, 63, 67, 75]. These models comprise partial integro-differential equations
in which changes in the tumour heterogeneity are mediated by proliferation,
competition and epigenetic mutations. Such frameworks have been used to
understand how oxygen gradients shape metabolic heterogeneity [23, 75] as well
as the consequences of resistance to chemotherapy [46, 76].

In our previous work [14], we developed a stemness-structured model to
investigate how heterogeneity in stemness levels affects the growth dynamics
of tumours in a well-mixed environment. In our model, a cell’s phenotypic
state is described by a continuous variable, $s$, such that $s=0$ correspond
to CSCs and $s=1$ to TDCs. Furthermore, cells may change their state in
response to random and environmental-induced epigenetic changes, which are
represented respectively via diffusion and advection along the $s$-axis. When
local oxygen levels are sufficiently high, the advection flux favours
differentiation; on the contrary, when oxygen levels are low the advection
term favours de-differentiation towards a stem-like phenotype. The
differential effects effects of radiotherapy (RT) are also included by
assuming that cell radio-sensitivity depends on the stemness level of cells;
in line with experimental evidence [3, 21, 30], CSCs are the most radio-
resistant subpopulation. In this paper, we extend our earlier work to account
for oxygen diffusion and consumption in the tumour and investigate how spatial
variation in local oxygen levels affects phenotypic (_i.e._ ,“stemness”)
heterogeneity within a tumour. Using a combination of numerical techniques, we
find that the interplay between physical and phenotypic space gives rise to
more complex dynamics than those observed in our earlier work [14]. In
particular, features of the spatial model, such as multistability and cusp
bifurcations, disappear when spatial effects are neglected. We also
investigate the combined effect of spatial and phenotypic heterogeneity on
radiation outcome. We use numerical simulations to compare the efficacy of
different treatment strategies and, in particular, their ability to prevent
CSC-driven relapse. In doing so, we provide insight into the role that dose
fractionation and scheduling may play in preventing relapse.

The remainder of this paper is organised as follows. In §2 we present our
stemness-structured spatially-resolved model and recast it in a partially non-
dimensional form. In §3, we present typical numerical simulations of the
model, which highlight the influence of the initial conditions in determining
the long-time behaviour of the system in the absence of treatment. The
observed multistability of the model is further investigated in §4, where we
numerically compute steady state solutions and use continuation techniques to
identify how equilibria and their linear stability change as key model
parameters vary. In §6, we investigate tumour characterised by different
levels of heterogeneity responds to standard radiotherapy protocols. how
applying treatment can drive the eradication of the tumour by perturbing an
initially stable state where the tumour persists. We find that the rate at
which tumour cells consume oxygen plays a critical role in dictating the
response to treatment, with high oxygen consumption rates, _i.e._ , higher
spatial-heterogenity of the oxygen levels in the tumour, correlating with poor
outcome. In §7, we summarise our findings and outline future research
directions. The appendices offer a number of details about the numerical
methods, the parameter values and the bifurcation analysis.

##  2 Model Development

Figure 1: Schematic representation of the model structure: we consider a fixed
slice of tissue delimited by two vessels at location $x=\pm L$. We reduce the
model to a 1D Cartesian geometry by assuming that there is no dependency on
the $y$ direction. Tumour cells living in this region are classified according
to their stemness level $s$ where $s=0$ corresponds to fully stem-like cells
while $s=1$ to cells that are terminally differentiated. At each location
$x^{*}$ in space the density of cell with phenotype $s$ is denoted by
$N(x^{*},s,t)$ whose time evolution is described by the system (2).

As shown in Fig. 1, we consider a fixed tissue region in which oxygen is
supplied by blood vessels located on the domain boundaries. For simplicity, we
consider a 1D Cartesian domain of length $2L$ and assume that the model
variables depend on the spatial location $x\ [\mbox{mm}]$, where $x\in[-L,L]$,
and time $t\,[hr]$. The state of a tumour cell in phenotypic space is
described by its stemness level $s$ (we view $s$ as a dimensionless
independent variable); stemness levels influence a cell proliferation
capacity, resistance to treatment (here radiotherapy) and response to
physiological stresses (here _hypoxia_). The local cell phenotypic
distribution $N=N(x,s,t)\ [\mbox{cells}\mbox{ mm}^{-1}]$ indicates the number
of tumour cells with phenotype $s$, at spatial location $x$ and time $t$. We
further introduce $\Phi=\Phi(x,t)\ [\mbox{cells}\mbox{ mm}^{-1}]$, which
represents the total number of tumour cells at location $x$ and time $t$:

$\Phi(x,t)=\int_{0}^{1}N(x,s,t)\,ds.$ (1)

The evolution of $N$ is dictated by cell proliferation and death as well as
spatial ($J_{x}=J_{x}(x,s,t)$) and phenotypic ($J_{s}=J_{s}(x,s,t)$) fluxes
along the $x$ and $s$ axis, respectively. Cell behaviour is mediated by local
levels of oxygen $C=C(x,t)\ [mmHg]$. As oxygen diffuses into the tissue, it is
consumed by cells at a rate that may depend on their phenotype, causing
different oxygen environments, or _niches_ , to form. Under these assumptions,
we have that the time evolution of $N$ and $C$ is determined by the following
system of coupled non-linear and non-local parabolic partial differential
equations (PDEs):

$\displaystyle\frac{\partial C}{\partial
t}=D_{xc}\frac{\partial^{2}C}{\partial
x^{2}}-\underbrace{\int_{0}^{1}\Gamma(s,C)N(x,s,t)\,ds}_{\mbox{net oxygen
consumption}},$ (2a) $\displaystyle\frac{\partial N}{\partial
t}+\frac{\partial J_{x}}{\partial x}+\frac{\partial J_{s}}{\partial
s}=F(s,C,\Phi,t)N.$ (2b) In Eq. (2a), the positive constant $D_{xc}\
[\mbox{mm}^{2}\mbox{hr}^{-1}]$ is the diffusion coefficient of oxygen and the
function $\Gamma=\Gamma(s,C)\ [\mbox{mmHg}(\mbox{cell}\mbox{ hr})^{-1}]$
represents the rate at which cells of phenotype $s$ exposed to local oxygen
levels $C$ consume oxygen. Focusing now on Eq. (2b), and following [25, 75],
we assume that cells simply move randomly in space, in which case the spatial
flux $J_{x}$ is given by: $\displaystyle J_{x}=-D_{xn}\frac{\partial
N}{\partial x},$ (2c) where the positive constant $D_{xn}\
[\mbox{mm}^{2}\mbox{hr}^{-1}]$ represents the diffusion coefficient of cells.
What distinguishes our model from previous works is the inclusion of the
phenotypic flux $J_{s}$, which comprises two effects: random epigenetic
mutations and environment-mediated changes in cell stemness. These are
represented, respectively, via diffusion and advection terms so that $J_{s}$
is given by: $\displaystyle J_{s}=-D_{sn}\frac{\partial N}{\partial
s}+V_{s}(s,C)N.$ (2d) Here, the positive constant $D_{sn}\ [\mbox{hr}^{-1}]$
represents the rate at which random epigenetic mutations take place, while the
velocity $V_{s}=V_{s}(s,C)\ [\mbox{hr}^{-1}]$ captures the rate at which cells
of phenotype $s$ alter their stemness level in response to local oxygen
levels. Focusing on the right-hand side of Eq. (2b), the function
$F=F(s,C,\Phi,t)\ [\mbox{hr}^{-1}]$ describes the net rate of cell
proliferation. As such, it also captures the effect of radiotherapy, which we
assume to be administered in $M$ doses $d_{i}\ [Gy]$ at discrete times $t_{i}$
with $i=1,\ldots,M$. Under these assumptions, we define $F$ as:
$\displaystyle\begin{aligned}
F(s,C,\Phi,t)=\underbrace{P(s,C)\left(1-\frac{\Phi}{\Phi_{max}}\right)}_{\mbox{proliferation}}&-\underbrace{K(s,C)}_{\mbox{cell
death}}\\\
&-\underbrace{\sum^{M}_{i=1}\log\left(\frac{1}{S_{RT}(s,C,\Phi;d_{i})}\right)\delta(t-t_{i})}_{\mbox{radiotherapy}}.\end{aligned}$
(2e) In Eq. (2e), $\Phi_{max}\ [\mbox{cells}/\mbox{mm}]$ represents the local
carrying capacity of the tissue, so that if the initial condition is such that
$\Phi(x,0)\in[0,\Phi_{max}]$ then $\Phi(x,t)\in[0,\Phi_{max}]$ for all $t>0$.
The associated logistic-growth term captures growth inhibition due to
competition for space. Locally, cells proliferate at a rate $P=P(s,C)\
[hr^{-1}]$ and die at a rate $K=K(s,C)\ [hr^{-1}]$ both rates being modulated
by the stemness levels $s$ and the local oxygen levels, $C$. The final sink
term in Eq. (2e) is associated with cell death due to radiotherapy (RT); here
$S_{RT}=S_{RT}(s,C,\Phi;d_{i})\,\in[0,1]$ denotes the fraction of cells that
survive exposure to a dose $d_{i}$ of RT. As discussed in [14], we use the
term $\log(1/S_{RT})$ although other forms have been proposed in the
literature. For example, in [43], the term $(1-S_{RT})$ is used. The improved
accuracy of the logarithmic functional form was argued in [14]. We postpone
the specification of functional forms in the model (_i.e._ , $V$,$P$, $K$, and
$SF$) until we have rescaled the model equations.

In line with the geometry presented in Fig. 1, we impose the following
boundary conditions:

$\displaystyle\left.D_{sn}\frac{\partial N}{\partial
s}-NV_{s}\right|_{s\in\left\\{0,1\right\\}}=0,$ $\displaystyle\quad
x\in[-L,L],\,t>0,$ (2f) $\displaystyle\left.\frac{\partial N}{\partial
x}\right|_{x=\pm L}=0,$ $\displaystyle\quad s\in(0,1),\,t>0,$ (2g)
$\displaystyle C(\pm L,t)=C_{\infty},$ $\displaystyle\quad t>0.$ (2h)

Eqs. (2f)-(2g) state that the phenotypic and spatial fluxes of the cells
vanish on the relevant domain boundaries. In Eq. (2h), the supply of oxygen
(from vessels located at tissue boundaries) is modelled by prescribing a
constant oxygen concentration, $C_{\infty}$, on $x=\pm L$. Finally, we close
the governing equations by imposing the following initial conditions:

$\displaystyle N(x,s,0)=N_{0}(x,s)$ $\displaystyle\quad
x\in[-L,L],\,s\in(0,1),$ (2i) $\displaystyle C(x,0)=C_{0}(x)$
$\displaystyle\quad x\in[-L,L].$ (2j)

#### Non-dimensional model formulation

We now recast the model in a partially dimensionless form; all variables are
rescaled except for time $t$ which is measured in hours. We continue to view
time as a dimensional variable to facilitate interpretation of model
simulations, particularly when RT is included. We choose the following
scalings for the other model variables

$X=\frac{x}{L},\quad
n=\frac{N}{\Phi_{max}},\quad\phi=\frac{\Phi}{\Phi_{max}},\quad
c=\frac{C}{C_{\infty}}.$ (3)

Substituting (3) into (2) and exploiting the assumed symmetry about the axis
$X=0$, we deduce that the evolution of the re-scaled cell phenotypic
distribution $n$ is governed by

$\displaystyle\frac{\partial n}{\partial
t}=\mathcal{D}_{xn}\frac{\partial^{2}n}{\partial
X^{2}}+\frac{\partial}{\partial s}\left(\mathcal{D}_{sn}\frac{\partial
n}{\partial s}-v_{s}(s,c)n\right)+\mathcal{F}(s,c,\phi,t)n,$
$\displaystyle\quad X\in(0,1),\ s\in(0,1),\ t>0,$ (4a)
$\displaystyle\left.\mathcal{D}_{sn}\frac{\partial n}{\partial
s}-nv_{s}\right|_{s\in\left\\{0,1\right\\}}=0,$ $\displaystyle\quad
X\in[0,1],\,t>0,$ (4b) $\displaystyle\left.\frac{\partial n}{\partial
X}\right|_{X\in\left\\{0,1\right\\}}=0,$ $\displaystyle\quad s\in(0,1),\,t>0,$
(4c) $\displaystyle n(X,s,0)=n_{0}(X,s),$ $\displaystyle\quad
X\in[0,1],\,s\in(0,1),$ (4d) where $\displaystyle\phi(X,t)$
$\displaystyle=\int_{0}^{1}n(X,s,t)\,ds,$ (4e)
$\displaystyle\mathcal{F}(s,c,\phi,t)$
$\displaystyle=\mathcal{P}(s,c)\left(1-\phi\right)-\mathcal{K}(s,c)-\sum^{M}_{i=1}\log\left(\frac{1}{\mathcal{S}_{RT}(s,c,\phi;d_{i})}\right)\delta(t-t_{i}),$
(4f)

and the evolution of the normalised oxygen concentration $c$ is dictated by

$\displaystyle\frac{\partial c}{\partial
t}=\mathcal{D}_{xc}\frac{\partial^{2}c}{\partial
X^{2}}-\int_{0}^{1}\gamma(s,c)n(X,s,t)\,ds,$ $\displaystyle\quad X\in(0,1),\
t>0,$ (5a) $\displaystyle c(1,t)=1,\quad\frac{\partial c}{\partial X}(0,t)=0,$
$\displaystyle\quad t>0,$ (5b) $\displaystyle c(x,0)=c_{0}(x),$
$\displaystyle\quad X\in[0,1].$ (5c)

In Eqs. (4)-(5), we have introduced the following (partially dimensionless)
parameter groupings and functional forms:

$\displaystyle\mathcal{D}_{xn}$
$\displaystyle=\frac{D_{xn}}{L^{2}},\quad\mathcal{D}_{xc}=\frac{D_{xc}}{L^{2}},\quad\gamma(s,c)=\frac{\Gamma(s,c\,C_{\infty})\Phi_{max}}{C_{\infty}},$
(6) $\displaystyle v_{s}(s,c)$
$\displaystyle=V_{s}(s,c\,C_{\infty}),\quad\mathcal{P}(s,c)=P(s,c\,C_{\infty}),\quad\mathcal{K}(s,c)=K(s,c\,C_{\infty}),$
$\displaystyle\mathcal{F}(s,c,\phi,t)$
$\displaystyle=F(s,c\,C_{\infty},\phi\,\Phi_{max},t),\qquad\mathcal{S}_{RT}(s,c,\phi;d_{i})=S_{RT}(s,c\,C_{\infty},\phi\,\Phi_{max};d_{i}).$

###  2.1 Specification of functional forms.

We now define the functional forms for $v_{s}$, $\mathcal{P}$, $\mathcal{K}$,
$\gamma$ and $\mathcal{SF}$. Some are similar to those used in [14] where we
have studied the spatially well-mixed counterpart of Eqs. (4)-(5). Therefore
we refer to [14] for a more detailed discussion and here focus on the key
differences.

To facilitate understanding of our model, we summarise in Fig. 2 the ways in
which we assume local oxygen levels influence tumour cell behaviours. This, in
turn, leads us to subdivide the tumour into four _niches_ according to local
oxygen levels: normoxia ($c_{H}<c<c_{\infty}$), hypoxia ($c_{N}<c<c_{H}$),
radio-biological hypoxia ($c_{N}<c<c_{R}$) and necrosis ($0\leq c<c_{N}$).

Figure 2: Schematic representing the role of oxygen in the model. Local oxygen
levels mediate proliferation as well as death, the differentiation of cells
and their sensitivity to radiotherapy. Note that the scaling of oxygen levels
is here for illustrative purposes. Values associated with the thresholds can
be found in Table 2.

We define the proliferation rate $\mathcal{P}$ and natural death rate
$\mathcal{K}$ in Eq. (4f) as follows:

$\displaystyle\mathcal{P}(s,c)=p_{0}(c)\exp\left[-\frac{s^{2}}{g_{0}}\right]+p_{1}(c)\exp\left[-\frac{(s-s_{max})^{2}}{g_{1}}\right],$
(7a)
$\displaystyle\mathcal{K}(s,c)=\underbrace{d_{f}\,e^{-k_{f}(1-s)}}_{\footnotesize\begin{subarray}{c}\mbox{natural
cell}\\\\[2.0pt]
\mbox{death}\end{subarray}}+\underbrace{d_{n}H_{\epsilon_{K}}(c_{N}-c)}_{\footnotesize\mbox{necrosis}},$
(7b) $\displaystyle\mbox{where
}p_{i}(c)=p_{i}^{max}\genfrac{}{}{}{}{\raisebox{2.0pt}{$c^{4}$}}{\raisebox{-2.0pt}{$K_{i}^{4}+c^{4}$}},\quad
i\in\left\\{0,1\right\\}.$ (7c) In Eq. (7a), a bimodal function describes the
proliferation rate $\mathcal{P}$; the two peaks correspond to CSCs ($s=0$) and
differentiated cells (DCs) ($s_{max}=0.55$) which proliferate at rates
$p_{0}(c)$ and $p_{1}(c)$ respectively. We assume further that both $p_{0}$
and $p_{1}$ are increasing, saturating functions of $c$ (see Eq. (7c)). In
physiological oxygen conditions ($c\approx 1$), CSCs and DCs proliferate at
their maximum rates, $p^{max}_{0}$ and $p^{max}_{1}$, respectively. The
parameters $K_{0}$ and $K_{1}$ represent the oxygen concentrations at which
CSCs and DCs proliferate at half their maximum rates (see Fig. 2). We fix
$p^{max}_{0}<p^{max}_{1}$ and $K_{0}<K_{1}$ so that DCs are more proliferative
than CSCs in normoxia, in line with [20, 40], but their growth is more
significantly inhibited in hypoxia ($c\lessapprox K_{1}\equiv c_{H}$). In
contrast, CSCs proliferate more slowly than DCs in normoxia, but they are less
sensitive hypoxia (since $K_{0}<K_{1}$) and continue to proliferate albeit at
a low rate. In this way, we capture the higher plasticity of CSCs and, in
particular, their ability to adapt to hypoxia [15, 31, 64]. In Eq. (7b), we
account for two mechanisms of cell death: natural cell death and necrosis.
While natural cell death is assumed to be independent of oxygen levels $c$ and
only to affect TDCs ($s\approx 1$), necrosis impact all cells regardless of
their stemness level $s$ only when oxygen levels are extremely low
($c\lessapprox c_{N}\ll 1$). To capture the effect of necrosis we introduce
the continuous step function, $H_{\epsilon}$:
$H_{\epsilon}(x)=\frac{1}{2}+\frac{1}{2}\tanh\left(x\epsilon^{-1}\right),$ (8)

where the parameter $\epsilon$ determines how rapidly $H_{\epsilon}$ switches
between its extremal values ($H_{\epsilon}(x)\approx 1$ for $x>\epsilon$ and
$H_{\epsilon}(x)\approx 0$ for $x<-\epsilon$). In Eq. (7b), we assume that the
transition between the hypoxic and necrotic niches is rapid by setting
$\epsilon_{K}=0.01$. In Fig. 3(a), we plot the net proliferation rate, in the
absence of competition (_i.e._ , $\mathcal{P}-\mathcal{K}$), and show how it
depends on the local oxygen levels ($c$) and phenotype ($s$).

(a)

(b)

Figure 3: Surface plot of the functional (a)
$\mathcal{P}(s,c)-\mathcal{K}(s,c)$ and (b) $v_{s}(s,c)$ with indication of
the oxygen threshold that dictates their behaviour. Note that we use a
logarithmic scale for the oxygen $c$. For the parameters used in the figure
see Table 2.

Following [14], we propose the following functional form for the phenotypic
advection velocity $v_{s}$:

$\displaystyle
v_{s}(s,c)=v^{+}_{s}(s)H_{\epsilon_{v}}(c-c_{H})-v^{-}_{s}(s)H_{\epsilon_{v}}(c_{H}-c),$
(9a) $\displaystyle
v^{+}_{s}(s)=\genfrac{}{}{}{}{\raisebox{2.0pt}{$V_{+}$}}{\raisebox{-2.0pt}{$V^{*}_{+}$}}\tanh\left(\genfrac{}{}{}{}{\raisebox{3.0pt}{$s^{\omega_{+}}$}}{\raisebox{-3.0pt}{$\xi_{+}$}}\right)\tanh\left(\genfrac{}{}{}{}{\raisebox{2.0pt}{$1-s$}}{\raisebox{-2.0pt}{$\xi_{+}$}}\right),$
(9b) $\displaystyle
v^{-}_{s}(s)=\genfrac{}{}{}{}{\raisebox{2.0pt}{$V_{-}$}}{\raisebox{-2.0pt}{$V^{*}_{-}$}}\tanh\left(\genfrac{}{}{}{}{\raisebox{2.0pt}{$s$}}{\raisebox{-2.0pt}{$\xi_{-}$}}\right)\tanh\left(\genfrac{}{}{}{}{\raisebox{3.0pt}{$(1-s)^{\omega_{-}}$}}{\raisebox{-3.0pt}{$\xi_{-}$}}\right),$
(9c)

In Eqs. (9b)-(9c), the normalizing factors $V_{\pm}^{*}$ are defined so that
$\max\limits_{s\in[0,1]}(v_{s}^{\pm})=V_{\pm}$, where $V_{\pm}$ correspond to
the maximum absolute values of the velocities with which cells differentiate
($V_{+}$) and de-differentiate ($V_{-}$). The positive constants $\xi_{\pm}$
determine how rapidly the velocity profiles switch near $s=0$ and $s=1$. The
exponents $\omega_{\pm}$ are positive integers: $\omega_{+}=1$ and
$\omega_{-}=2$ (this choice is based on the analysis performed in [14]). As in
[58], we assume that there is a threshold oxygen concentration $c_{H}$ that
separates the tumour into two _niches_. As illustrated in Fig. 3(b), when
$c>c_{H}$ cells differentiate and progress towards $s=1$ (i.e., $v_{s}>0$);
when $c<c_{H}$, cells instead de-differentiate and assume stem-like traits
(i.e., $v_{s}<0$). There is a smooth transition between the two niches,
capture by the function $H_{\epsilon_{v}}$, which is defined as in Eq. (8),
with $\epsilon_{v}=0.05$.

To complete the model formulation for tumour growth in the absence of
treatment, we must specify the functional form of the oxygen consumption rate
$\gamma=\gamma(s,c)$. As already mentioned, in practice oxygen consumption
might depend on the stemness variable $s$. However, in the absence of suitable
experimental data to define this dependence, we follow [62] and assume
$\gamma=\gamma(c)$. Furthermore, following [34, 44], we suppose that $\gamma$
drops rapidly when oxygen levels fall below the necrotic threshold $c=c_{N}$:

$\gamma=\gamma_{0}H_{\epsilon_{\gamma}}(c-c_{N}).$ (10)

In Eq. (10), $\gamma_{0}$ is the maximum oxygen consumption rate for cells and
$\epsilon_{\gamma}=0.05$. We rewrite Eq. (10) in terms of the diffusion
coefficient $\mathcal{D}_{xc}^{-1}$ by introducing the non-dimensional
parameter $\gamma^{*}=\gamma_{0}\mathcal{D}_{xc}^{-1}$. The latter gives the
ratio between the timescales for oxygen diffusion and oxygen consumption. As
we will see in §4, this parameter plays a key role in determining the long
time behaviour of the system.

#### The survival fraction

Sensitivity of cells to radiotherapy is modulated by environmental conditions
[65, 69] as well as cell-intrinsic factors [3, 30, 72, 73]. We start by
considering the standard linear-quadratic model which states that the fraction
of cells that survive a single dose $d$ of radiation is given by [51]:

$\displaystyle\ln\left(\mathcal{S}_{RT}(c,s,\phi;d)\right)=-\alpha(c,s,\phi)d-\beta(c,s,\phi)d^{2}.$
(11a) In what follows, we suppose that $\alpha[\mbox{Gy}^{-1}]$ and
$\beta[\mbox{Gy}^{-2}]$ are dependent on the local environment, and, in
particular, the oxygen concentration $c$, the cell density $\phi$, and the
stemness of the cell, $s$. Consequently, in our model the overall tumour
radio-sensitivity depends on the tissue composition and the environment, both
of which are themselves affected by treatment.

We decompose $\alpha$ (and, analogously, $\beta$) into three factors, which
account for different mechanism of radio-resistance:

$\displaystyle\alpha(s,c,\phi)=\alpha_{1}(c)\alpha_{2}(s)\alpha_{3}(c,s,\phi),\quad\beta(s,c,\phi)=\beta_{1}(c)\beta_{2}(s)\beta_{3}(c,s,\phi).$
(11b)

First, it is well known that low oxygen levels ($c<c_{R}$) can adversely
impact radiotherapy efficacy; oxygen is an essential reactant for the fixation
of radio-induced DNA damage and subsequent cell death [78, 51, 65, 79]. This
effect is commonly incorporated into the survival fraction by introducing the
_oxygen enhancement ratio_ (OER)[78, 43]. The OER is a constant parameter
whose value usually ranges between $1$ and $4$; it is defined to be the ratio
of the radiation dose that needs to be given in anoxia compared to normoxia to
achieve the same biological effects (i.e., the same overall survival fraction)
[78]. As in [43], we assume that the OER plays a role only in regions of
_radio-biological hypoxia_ (i.e., regions where $c\leq c_{R}$), and propose
the following functional form for $\alpha_{1}(c)$ and $\beta_{1}(c)$:

$\displaystyle\alpha_{1}(c)=\begin{cases}1,&\quad c>c_{R}\\\
\mbox{OER}^{-1},&\quad c\leq
c_{R}\end{cases},\quad\beta_{1}(c)=\begin{cases}1,&\quad c>c_{R}\\\
\mbox{OER}^{-2},&\quad c\leq c_{R}\end{cases}$ (11c)

where the threshold $c_{R}$ is chosen such that $c_{N}<c_{R}<c_{H}$ (see Fig.
2). As in [14], we account for the intrinsic radio-resistance of CSCs [3, 21,
30, 55] by proposing the following functional forms for $\alpha_{2}(s)$ and
$\beta_{2}(s)$:

$\displaystyle\alpha_{2}(s)=\alpha_{-}+\Delta\alpha\tanh(\kappa
s),\quad\beta_{2}(s)=\beta_{-}+\Delta\beta\tanh(\kappa s),$ (11d)

In Eq. (11d), the positive constants $\alpha_{-}$ and $\beta_{-}$ represent
the radio-sensitivity of CSCs (_i.e._ , $s\approx 0$) while $\Delta\alpha$ and
$\Delta\beta$ denote the difference in radio-sensitivity between CSCs and more
differentiated cell phenotypes ($s>0.5$). In contrast to [14], we include an
additional term to account for the increased sensitivity of rapidly
proliferating cells compared to growth-arrested cells [10, 22, 28] by assuming
that $\alpha_{3}(c,s,\phi)$ and $\beta_{3}(c,s,\phi)$ depend on the effective
cell proliferation rate in the following way:

$\displaystyle\alpha_{3}(c,s,\phi)=\beta_{3}(c,s,\phi)=\left(1+\frac{\mathcal{P}(c,s)(1-\phi)}{\mathcal{P}_{max}}\right).$
(11e)

In Eq. (11e), $\mathcal{P}_{max}=\max_{c,s}\mathcal{P}=p_{0}^{max}$. Equation
(11e) implies that rapidly proliferating cells are more sensitive to radiation
than cells which proliferate less rapidly.

The combined effect of the three radio-resistance mechanisms is illustrated in
Fig. 4. Note that we show plots of $\alpha$ only since $\beta$ has a
qualitatively similar effect. When oxygen levels drop below $c_{R}$ (see
lighter line in Fig 4), the radio-sensitivity of the cells also falls. For
milder levels of hypoxia (e.g., $c_{R}<c<c_{H}$), cell proliferation is
inhibited for all phenotypes $s$ so that the dominant mechanism modulating
radio-sensitivity is the intrinsic radio-sensitivity of the cells (_i.e._ ,
$\alpha\sim\alpha_{2}(s)$). Spatial competition is relevant at higher oxygen
concentration (e.g, $c>c_{H}$). Where spatial competition is low ($\phi\approx
0$), differentiated cells ($s\approx 0.5$) are highly radio-sensitive; by
contrast, TDCs are less sensitive as they are not proliferating. However, as
the population grows, competition for space inhibits cell proliferation so
that TDCs and differentiated cells have similar responses to radiotherapy.

Figure 4: Series of plot showing how, for fixed value of $\phi$,
$\alpha=\alpha(s,c,\phi)$ (defined by Eq. (11a)) changes as the phenotypic
variable $s$ and the oxygen concentration $c$ vary. We plot $\alpha$ for
different values of the re-scaled cell density $\phi$ to show also how
competition for space modulates the radio-sensitivity. Parameter values set as
in Tables 1-2.

###  2.2 Numerical methods and interpretation of the numerical simulations

When performing dynamic simulations of Eqs (4)-(5), most parameters are fixed
to the values listed in Tables 1-2; where possible these values have been
taken from the literature.

Based on standard values reported in the literature and considering a tissue
region that is several micrometers wide ($L\approx 200\,\mu m$), we find that
the time-scale of oxygen diffusion in the tumour region
($\mathcal{D}_{xc}^{-1}=L^{2}D_{xc}^{-1}$) is of the order of seconds or
minutes, which is much faster than the time-scale of cell proliferation
($\mathcal{P}^{-1}$ of the order of hours and days) which is the timescale of
interest here. We exploit this separation of timescales to simplify the
equations governing the oxygen dynamics. As is standard in mathematical models
of tumour growth [12, 44], we assume that the oxygen distribution is in quasi-
equilibrium so that $c$ satisfies

$\frac{\partial^{2}c}{\partial X^{2}}-\gamma^{*}\phi(X,t)H(c-c_{N})=0,\quad
c(1,t)=0,\quad\frac{\partial c}{\partial X}=0.$ (12)

Under the quasi-steady-state approximation, the oxygen quickly adapts is
spatial distribution to changes in the tumour composition. As a result, we do
not need to impose initial conditions for $c$. Under this assumption, the
oxygen profile is primarily characterised by the non-dimensional parameter
$\gamma^{*}=\gamma_{0}L^{2}D_{xc}^{-1}$, which we recall measures the ratio of
the time scales for oxygen diffusion in a tissue of size $L$ to the time scale
for oxygen consumption ($\gamma_{0}^{-1}$). Since we are interested in
understanding how spatial variation in oxygen levels affects the tumour’s
phenotypic composition, we view $\gamma^{*}$ as a free parameter, which varies
between different tumour cell populations. While the oxygen diffusion
coefficient is typically considered as a constant $D_{xc}$ (see Table 1),
oxygen consumption rates ($\gamma_{0}$) are likely to be tumour- and patient-
specific. For this reason, when comparing different values of $\gamma^{*}$, we
will implicitly assume that the tumours have the same ‘typical’ length scale,
$L$, and oxygen diffusion coefficient but differ in their oxygen consumption
rates $\gamma_{0}$.

Unless otherwise stated, we initialise the simulations by assuming that, at
$t=0$, the cells are distributed uniformly in space and normally along the
stemness axis, with mean stemness $s_{0}\in(0,1)$ and variance $\sigma_{s}>0$.
Thus, we prescribe

$n(s,X,t=0)=n_{0}(s,X)=\frac{M_{0}}{\sqrt{2\pi}\sigma_{s}}\exp\left[-\frac{(s-s_{0})^{2}}{2\sigma^{2}_{s}}\right].$
(13)

where $M_{0}>0$ is the initial tumour burden, while $s_{0}$ and
$\sigma_{s}^{2}$ are chosen so that the mass of the Gaussian profile lies
inside the stemness domain $s\in(0,1)$ and the Neumann boundary conditions
(4b) are satisfied, up to a negligible (exponentially small) correction.

We use the method of lines to solve Eqs. (4)-(12). In more detail, we first
discretise Eqs. (4) and (12) using finite volume and central finite
differences to discretise the phenotypic ($s$) and spatial ($X$) axis,
respectively. We then advance Eq. (4a) in time by using an implicit multi-step
time-stepping method based on backward differentiation formulae (BDF)
implemented in the _scipy_ python package, which is suitable for stiff
equations. At each time point $t$, given the discrete form of the solution
$n(s,X,t)$, we can compute $\phi(X,t)$. We use the latter to solve for the
stationary distribution of the oxygen concentration $c$ via a fixed point
iteration method. A detailed description of the numerical method used can be
found in A.

#### Interpretation of the numerical results

The model has several independent variables which can make it difficult to
simultaneously visualise the evolution of the tumour’s spatial and phenotypic
distributions.

As a measure of tumour load, we introduce the total tumour burden $M=M(t)$:

$M(t)=\int_{0}^{1}\int_{0}^{1}n(X,s,t)dXds.$ (14)

We note that $M\Phi_{max}L$ represents the total number of tumour cells in the
slice of tissue.

Figure 5: Illustration of the relationship between the definition of the three
discrete sub-population, (namely $P_{1}$ (CSCs), $P_{2}$ (DCs) and $P_{3}$
(TDCs)) and the cells maximum proliferation potential (_i.e._ ,
$\mathcal{P}$), natural death (_i.e._ , $\mathcal{K}$) radio-sensitivity
(_i.e._ , $\alpha$). The yellow and blue dotted line indicate, respectively
$s=s_{1}=0.3$ and $s=s_{2}=0.8$.

When studying the tumour phenotypic composition, it will be convenient to
group cells into three phenotypic subpopulations: $P_{1}$, $P_{2}$ and
$P_{3}$. As shown in Fig. 5, the subdivision is based on the cells’
proliferative capacity ($\mathcal{P}$) and their sensitivity to natural cell
death ($\mathcal{D}$) and treatment ($\alpha$ and $\beta$). Cells in $P_{1}$
($0\leq s\leq s_{1}$) are highly-resistant stem-like cells, while cells in
$P_{2}$ ($s_{1}<s\leq s_{2}$) are highly proliferative and radio-sensitive.
Finally, $P_{3}$ ($s_{2}<s\leq 1$) consists of terminally-differentiated cells
that have exited the cell-cycle and have a high death rate $\mathcal{D}$. For
brevity, henceforth, we refer to $P_{1}$ as CSCs, $P_{2}$ as DCs and $P_{3}$
as TDCs. We emphasise that this definition of the sub-populations is not a
“firm distinction”; the thresholds $s_{1}$ and $s_{2}$ are chosen arbitrarily
to showcase the dynamics of the respective populations and simplify
interpretation of the simulation results. In particular, we note that small
changes in the values of $s_{1}$ and $s_{2}$ do not significantly affect our
results and/or conclusions.

Given the solution $n(X,s,t)$ of Eqs. (4), the proportion of cells in each of
the three classes $P_{1}$, $P_{2}$ and $P_{3}$ at time $t$ is given by

$\displaystyle\Pi_{i}(t)=\frac{1}{M(t)}\int_{s_{i-1}}^{s_{i}}\int_{0}^{1}n(X,s,t)dXds,\quad
i=1,2,3,$ (15)

where, without loss of generality, we fix $s_{0}=0$, $s_{1}=0.3$, $s_{2}=0.8$
and $s_{3}=1$. We note that, by definition, the three fractions are not
independent since they satisfy the condition: $\sum_{m=1}^{3}\Pi_{m}(t)\equiv
1$.

While, the fractions $\Pi_{i}$ provide information about the overall tumour
phenotypic composition, they lack information about the spatial location of
the different sub-populations. Therefore, we also define the local cell
fractions:

$\displaystyle\pi_{i}(t,X)=\frac{1}{M(t)}\int_{s_{i-1}}^{s_{i}}n(X,s,t)ds,\quad
i=1,2,3,$ (16)

to study the composition of cells in the different niches that forms within
the tumour as a result of spatial heterogeneity in the oxygen levels.

##  3 Dynamics in the absence of treatment

We start our analysis by presenting results from dynamic simulations of our
tumour growth model, in the absence of treatment. These results highlight one
of the key-features of our model: multi-stability.

Figure 6: Numerical results generated from solving the full model, defined by
Eqs. (4) and (13) with $s_{0}=0.5$ and two different values of $M_{0}$
($M_{0}=0.001$ on the blue curve, and $M_{0}=0.0001$ on the red curve). All
other model parameters are fixed at the default values specified in Tables
1-2, with $\gamma^{*}=4$ and $V_{+}=8\times 10^{-4}$. When $M_{0}=0.001$ (blue
curve), the tumour persists and eventually evolves to a non-trivial steady
state; when the initial tumour burden is decreased ($M_{0}=0.0001$, red
curve), the tumour grows initially but at long times becomes extinct. The
yellow dots indicate the time points at which the snapshots in Fig. 7 are
taken.

In Fig. 6, we show how the tumour burden $M(t)$ evolves for two simulations in
which all parameters are the same except for the initial tumour burden
$M_{0}$. When $M_{0}$ is sufficiently large (e.g., $M_{0}=0.001$), the tumour
invades the tissue; for smaller values of $M_{0}$ (e.g., $M_{0}=0.0001$) the
tumour cells become extinct (see also Fig. 7).

Despite the different long-time behaviour of $M$ for these two examples their
initial dynamics (until about $t=9$ weeks) are qualitatively similar. In both
cases, the tumour grows until $t\approx 6$ weeks. Thereafter, the tumour with
the smaller initial burden, decreases monotonically towards extinction, while
the tumour with the higher burden resumes growth and eventually attains a non-
trivial steady state.

(a) $M_{0}=0.0001$

(b) $M_{0}=0.001$

Figure 7: Time series plots showing the local cell distributions $n(x,s,t)$
for the two scenarios in Fig. 6. Solutions are computed by numerically solving
Eqs. (4) with initial conditions defined by Eq. (13), and $s_{0}=0.5$. We
observe long-time extinction when $M_{0}=0.0001$ (Fig. 7(a)) and evolution to
a persistent tumour cell population when $M_{0}=0.001$ (Fig. 7(b)). Parameter
values: as per Fig. 6.

To better understand why the dynamics diverge, we consider the time evolution
of the cell distribution $n(X,s,t)$ (see Fig. 7). In both examples, the tumour
is initially composed of highly-proliferative differentiated cells thus the
increase in the tumour burden $M(t)$. However, as these cells proliferate,
they also become terminally differentiated. At about 6 weeks, the tumour is
primarily composed of TDCs (note the large value of $n(X,1,t=6)$ in Figs.
7(a.2) and 7(b.2)). Since the death rate of TDCs exceeds their proliferation
rate (see Fig. 3(a)), $M(t)$ start to decrease. Further, while TDCs do not
proliferate, they still consume oxygen. Consequently, if sufficiently many
TDCs are present, an hypoxic region forms at the centre of the tumour (where
$x\approx 0$) as observed in Fig. 7(b.2). In the hypoxic _niche_ , cell
differentiation slows down and, if oxygen levels are sufficiently low, cells
de-differentiate. The different signals experienced by tumour cells in the
normoxic and hypoxic niches select for different phenotypes and drive spatial
variation in the phenotypic composition of the tumour which is amplified over
time (compare the cell distributions at $t=6$ and $t=60$ weeks, for example).
By contrast, in Fig. 7(a.2), the number of TDCs are low so that oxygen levels
remain above the hypoxic threshold $c_{H}$. In this case, $n(X,s,t)$ is
spatially homogeneous and the tumour comprises only TDCs which are exposed to
normoxia. As a result, the tumour is eventually driven to extinction, and the
rate at which $M(t)$ decays matches $d_{f}$, the rate of at which TDCs die
under normoxia.

Figure 8: Series of plots showing how, for the simulation illustrated in Fig.
7(b), the tumour composition changes over time. All parameters are fixed at
the default values stated in Tables 1-2 with $\gamma^{*}=4$ and $V_{+}=8\times
10^{-4}$ and $M_{0}=0.001$. (a) Cumulative plot indicating the fractions
$\Pi_{i}(t)$, with $i=1,2,3$, as defined in Eq. (15); (b)-(d) local
subpopulation distribution $\pi_{i}(x,t)$ as defined by Eqs. (16).

Having understood the origin of the divergent behaviour observed in the two
examples, we focus now on the invasion scenario (_i.e_ , $M_{0}=0.001$). As
mentioned above, it is the presence of de-differentiated cells in the hypoxic
region that drives secondary regrowth of $M(t)$ for $t>9$ weeks. This is
apparent, when looking at Fig. 8, where we show how the tumour’s phenotypic
composition evolves over time. When the hypoxic region first forms, the DCs
localise in the hypoxic niche (see Fig. 8(c)) while the normoxic region
contains only TDCs (see Fig. 8(d)). While DCs can proliferate, under hypoxia
they stop proliferating and become quiescent. Consequently, between 6 and 9
weeks, TDCs in the outer normoxic region tend to die while DCs in the hypoxic
niche survive but do not proliferate. Hence the net growth rate is negative.
At later times, the DCs diffuse from the hypoxic to the normoxic region, and
in the presence of increased oxygen levels they resume proliferation and drive
the secondary increase in $M(t)$. This initiates a feedback loop: as $M(t)$
increases, the size of the hypoxic region also increases, creating a niche
which is favourable for the development of CSCs, which drives further tumour
growth. As shown in Fig. 8(a), CSCs appear in the tissue only at later times
($t\approx 15$ weeks) and are localised in hypoxic regions (see $\pi_{1}(X,t)$
in Fig. 8(b)). At long times, the fraction of CSCs, $\Pi_{1}(t)$, increases
and asymptotes to a stationary value of $\approx 25\%$. The TDCs remain the
dominant cell phenotype, comprising $\approx 50\%$ of the tumour cells. While
CSCs and TDCs are located in distinct tumour regions (the hypoxic and normoxic
niches, respectively), DCs eventually spread uniformly across the tissue (Fig.
8(c)). Therefore, our model predicts coexistence of different cell populations
not only within the same tumour, due to the formation of different niches, but
also at the same spatial location, due to the cells’ ability to move. The
interplay between these two effects sustains tumour heterogeneity and is what
eventually drives not only the formation but also the growth of the tumour.

##  4 Equilibrium states in the absence of treatment and their stability.

The numerical results presented in Fig. 6 suggest that the model admits
multiple stable steady state solutions, at least for some parameter sets.
Consequently, the long time behaviour of the system depends on the initial
conditions, and their positioning with respect to the basins of attraction of
the stable steady states. To investigate the possible model outcomes, we
compute the steady states $(n_{\infty}(X,s),c_{\infty}(X))$ of the full non-
dimensional model, Eqs. (4)-(5), and study their stability. We note that, in
so doing, we consider the full, time-dependent model, with no quasi-steady
assumption for the oxygen.

Figure 9: Bifurcation diagram for the system (4) as a function of the
parameter $\gamma^{*}$ for fixed value of $V_{+}$. We differentiate between
stable (solid) and unstable (dashed) branches of the diagram and highlight
with a circle the relevant bifurcation point.

It is straightforward to show that Eqs. (4)-(5) admit a unique homogeneous
steady state (independent of $s$ and $X$). This is the trivial solution, with
$n_{\infty}\equiv 0$ and $c_{\infty}\equiv 1$, and corresponds to the “tumour
extinction” (disease-free) scenario. In C.1, we show that its linear stability
can be determined from that of the simpler, well-mixed model formulation of
the model which was investigated in [14]. This is because, when decomposing
the spatial perturbation into Fourier modes with wavenumber $\omega$, the
particular wavenumber $\omega=0$ (corresponding, effectively, to the well-
mixed case) is present. Further, since it is the most unstable wavenumber, it
determines the stability of the trivial steady state. The well-mixed model is
obtained by neglecting spatial variation in Eqs. (4), i.e., setting
$\mathcal{D}_{xn}$, and considering an homogeneous oxygen distribution
$c\equiv c_{\infty}$. While all other model parameters influence the stability
of the trivial steady state, we focus on parameters associated with the
advection velocity $v_{s}$ and fix all other parameters at values based on
estimates from the literature (see Tables 1-2). Since $c\equiv 1$ when spatial
heterogeneity is neglected, we have that $v_{s}(s)\approx v_{s}^{+}(s)$. As
shown in Fig. 9, when $0<V_{+}\ll 1$, the trivial steady state is unstable and
a non-trivial stable solution exists (see the upper green branch). As $V_{+}$
increases , $M_{\infty}$ decreases until $V_{+}$ crosses the critical point
$V_{+}^{(cr)}\approx 6.06\times 10^{-4}$ where the two solution branches cross
and exchange stability via a transcritical bifurcation. For $V>V_{+}^{(cr)}$
the trivial steady state is the unique, non-negative (i.e., physically-
realistic) solution and it is stable.

The bifurcation diagram becomes more complex when we account for spatial
heterogeneity. However, this will not affect the stability of the trivial
steady state. In other words, the location of the critical point
$V_{+}^{(cr)}$ is independent of $\mathcal{D}_{xn}$ or the parameter
associated to the oxygen distribution, i.e., $\gamma_{*}$ and
$\mathcal{D}_{xc}$.

###  4.1 Bifurcation diagrams: bistability

Computing heterogeneous stationary steady state solutions is challenging and,
therefore, we rely on numerical continuation techniques. In particular, we use
Julia’s package BifurcationKit [74] to compute solutions, study their
stability and identify bifurcation points (BifurcationKit is an open-source
software package for computing bifurcation diagrams). More details about the
numerical techniques used can be found in A. Before continuing, we note that
we present below all of the positive steady state that were identified. These
solutions were also observed when we performed dynamic numerical simulations;
however we cannot, a priori, exclude the possibility that other stationary
steady state solutions, not detected by our analysis, or other type of
attractors (such as limit cycles) might exist.

As mentioned above, the parameter $V_{+}$ plays a critical role in the long
time system dynamics. We now show how the diagram in Fig. 9 changes when we
account for spatial structure. In doing so, we compute bifurcation diagrams
for different values of the dimensionless oxygen consumption rate
$\gamma^{*}$, which affects the equilibrium oxygen distribution.

Figure 10: Series of diagrams showing how the bifurcation structure of system
(4) changes as the parameter $V_{+}$ varies, for three fixed values of
$\gamma^{*}$: (a) $\gamma^{*}$ = 1.65, (b) $\gamma^{*}$ =2.00 and (c)
$\gamma^{*}$ =4.00. We distinguish between stable (solid) and unstable
(dashed) solution branches and highlight bifurcation points with circles.

The results presented in Fig. 10 show how, as $\gamma^{*}$ increases through
the critical value $\gamma^{(cr)}\approx 1.62$, the qualitative properties of
the bifurcation diagram change. Comparison of Fig. 10 with the bifurcation
diagram for the well-mixed scenario (i.e., Fig. 9), we observe new steady
state solutions and new bifurcation points. In particular, the two saddle-node
bifurcation points, $\mbox{SN}_{1}$ and $\mbox{SN}_{2}$, define the boundary
of a bistable region. Focusing on Fig. 10(a), we see that, for small $V_{+}$,
there exists a single stable steady state, which is characterised by a large
tumour burden. As $V_{+}$ increases, the system undergoes a saddle node
bifurcation ($\mbox{SN}_{2}$). For larger values of $V_{+}$, the system admits
four non-negative steady state solutions; two of these solutions are stable
and characterised by non-zero tumour burdens. As $V_{+}$ increases further,
there is a second saddle-node bifurcation ($\mbox{SN}_{1}$), beyond which a
unique, non-trivial, stable steady state exists whose mass decreases with
$V_{+}$ until it crosses the trivial branch when $V_{+}=V_{+}^{(cr)}$ (the
yellow circle in Fig. 10(a)). Beyond the transcritical bifurcation, the
trivial steady state (with $M_{\infty}=0$) is the only stable solution. For
larger values of $\gamma^{*}$ (see Figs. 10(b) and 10(c)), the tumour burden
on the upper-most branch decreases and $\mbox{SN}_{1}$ moves to the right,
passing beyond $V_{+}^{(cr)}$. Note that in Fig. 10(c), $\mbox{SN}_{1}$ does
not appear as it is located beyond the interval under consideration,
$V_{+}\in[10^{-4},10^{-3}]$. The second saddle-node bifurcation point,
$\mbox{SN}_{2}$, also moves to the right (i.e. higher $V_{+}$) but it is less
sensitive to changes in $\gamma^{*}$. In particular, comparing Fig. 10(b) and
Fig. 10(c), we see that $\mbox{SN}_{2}$ approaches the transcritical
bifurcation point. Further numerical analysis of the fold points (see Fig. 22
in C.3) suggests that the transcritical and fold bifurcations intersect at
infinity (i.e., $\mbox{SN}_{2}(\gamma^{*})\rightarrow V_{+}^{(cr)}$ as
$\gamma^{*}\rightarrow\infty$).

Figure 11: Series of bifurcation diagrams of the system (4) with $\gamma^{*}$
being our bifurcation parameter. We compute the diagrams for three fixed
values of $V_{+}$: (a) $V_{+}=2.5\times 10^{-4}$, (b) $V_{+}=5\times 10^{-4}$
and (c) $V_{+}=8\times 10^{-4}$. We distinguish between stable (solid) and
unstable (dashed) branches of solutions and highlight with circles the
bifurcation points.

To better understand the role played by $\gamma^{*}$, we fix $V_{+}$ and
compute the system steady states, viewing $\gamma^{*}$ as our continuation
parameter. The corresponding diagrams, presented in Fig. 11, show that for
slow advection velocities (i.e., $V_{+}=2.5\times 10^{-4}$), the bi-stability
region is narrow, and the unstable branch connecting the two fold points is
almost vertical. As $V_{+}$ increases to $5\times 10^{-4}$ (see Fig. 11(b)),
$\mbox{SN}_{2}$ shifts towards the bottom-right of the diagram and the size of
the bistable region increases. The position of $\mbox{SN}_{1}$ moves in a
similar manner, but the effect is less pronounced. We note that the upper
stable branch is non-monotonic, with $M_{\infty}$ initially increasing past
the supercritical saddle-node point ($\mbox{SN}_{1}$) and subsequently
decreasing as $\gamma^{*}$ increases further. This contrasts with the lower
stable branch, which decreases monotonically until it disappears at the
saddle-node bifurcation point $\mbox{SN}_{2}$. In Fig. 11(c), we take
$V_{+}>V_{+}^{(cr)}$ (i.e., $V_{+}$ is beyond the transcritical bifurcation
point shown in Fig. 10). Comparing Figs. 11(b) and. 11(c), we see that the
middle stable branch disappears and the trivial steady state becomes stable.
This raises the question of what happens to the second saddle-node point
($\mbox{SN}_{2}$). As $V_{+}\rightarrow V_{+}^{(cr)}$ from below,
$\mbox{SN}_{2}$ shifts to the right
($\mbox{SN}_{2}(\gamma^{*})\rightarrow\infty$). Since $\mbox{SN}_{2}$ becomes
infinite at a finite value of $V_{+}$, our analysis suggests that the
transcritical bifurcation and $\mbox{SN}_{2}$ annihilate each other at
infinity, i.e., at $\gamma^{*}\rightarrow+\infty$. This is validated by using
fold continuation methods which follow the locations of $\mbox{SN}_{1}$ and
$\mbox{SN}_{2}$ in the ($V_{+}$,$\gamma^{*}$) plane (the results are presented
in C.3). Consequently, in Fig. 11(c), where $V_{+}>V_{+}^{(cr)}$, the
bifurcation diagram simplifies. For $0<\gamma^{*}<1.93$, extinction is the
only possible long term solution; however as $\gamma^{*}$ increases beyond
$\mbox{SN}_{1}$ (i.e., $\gamma^{*}>1.93$), two outcomes are possible: either
the tumour goes extinct, or it persists, giving rise to what we term the
“tumour-invasion” steady state.

Figure 12: Characterization of the steady state solution of system (4). (a) We
reproduce the bifurcation diagram from Fig. 11(c) and indicate the equilibrium
points of interest (A,B,C); (b) we compare the tumour phenotypic composition
($\Pi_{i}$ see Eq. (15)) at points A, B, and C; (c) we illustrate the
equilibrium solution, i.e., oxygen profile $c_{\infty}$, re-scaled cell number
density $\phi_{\infty}$, and the local distribution of cell phenotype
$n_{\infty}(X,s)$. In the surface plot for $n$, the vertical orange and red
lines indicate the boundaries of the hypoxic ($X_{H}$) and necrotic regions
($X_{N}$), respectively. The parameter used to generate these results are as
follows: $V_{+}=8\times 10^{-4}$ and $\gamma^{*}=2.0$ (A,B) or
$\gamma^{*}=4.0$ (C).

#### Characterisation of the computed equilibria

Thus far, we have identified the number of steady state solutions of Eqs. (4)
and determined their stability. We now characterise these solutions by
focusing on a couple of characteristic examples. In particular in Fig. 12, we
compare the cases of low ($\gamma^{*}=2.0$) vs high ($\gamma^{*}=4.0$) oxygen
consumption rates, while fixing $V_{+}=8.0\times 10^{-4}$. When
$\gamma^{*}=2.0$ (i.e., just to the right of the critical saddle node
bifurcation $\mbox{SN}_{1}$, see red dot in Fig. 12(a)), equilibria A and B
have a similar tumour burden, $M_{\infty}$, but their internal phenotypic
compositions differ significantly (see Fig. 12(b)). While approximately $20\%$
of cells for the stable steady state A are stem-like, CSCs are completely
absent from the unstable steady state B. Comparison of the spatial
distributions (see Fig. 12(c)), reveals that steady state A is characterised
by a large hypoxic niche in which CSCs accumulate (here
$\min_{X}c_{\infty}(X)<c_{H}$). While steady state B also possesses an hypoxic
region (here mild hypoxia as $\min_{X}c_{\infty}(X)\approx c_{H}$), this is
smaller and does not contain CSCs; we find instead coexistence of TDCs and DCs
in the hypoxic niche, where DCs are uniquely localised in the hypoxic region
and absent from the normoxic tissue regions. Plots of $\phi_{\infty}(X)$,
indicate that the rescaled cell density is almost constant for steady state A,
whereas for B it increases towards the tumour interior (where $X\approx 0$).
This is because cell growth is concentrated in the hypoxic niche and dominated
by slowly proliferating DCs.

Focusing on the upper stable branch in Fig. 12(a), we see that the tumour
burden, $M_{\infty}$, decreases monotonically for sufficiently large values of
$\gamma^{*}$. This is due to the formation of a necrotic region in which cells
die (see solution C). Overall, the tumour phenotypic composition for steady
state C is similar to that of A (see Fig. 12(b)), with steady state C
containing a slightly higher fraction of CSCs. Nonetheless, there is a marked
difference in the equilibrium cell densities $\phi_{\infty}(X)$. For steady
state C the maximum value of $\phi_{\infty}$ occurs at the oxygenated boundary
($X\approx 1$) and it decreases monotonically towards the interior, where the
oxygen levels drops below $c_{N}$ (i.e., the necrotic region). This suggests
that, in our model, while necrosis decreases the overall tumour burden, it
does not significantly impact the phenotypic heterogeneity of the tumour.
However, as we will show in the next sections, the presence of a necrotic
region in the tumour still correlates with poor treatment outcomes in contrast
to a scenario like A, where only hypoxic niches are present.

Our analysis pertains to the local stability of the different steady states of
the system. From a mathematical point of view, it is possible to drive the
tumour to extinction, if the trivial solution is stable (as for $V_{+}=8\times
10^{-4}$) and if the system enters the basin of attraction of the trivial
steady state. The question of interest is, therefore, whether it is possible
to apply treatment, specifically radiotherapy, in a way that modifies the long
time system dynamics and drives the tumour to extinction. Given the properties
of the unstable solution branch (see B in Fig. 12), we might expect that the
transition between the “tumour-invasion” and “tumour-extinction” regimes is
mediated by the formation of a sufficiently large hypoxic region. Even if
resistant CSCs are present, they would differentiate under normoxia. We
conclude that a successful treatment should avoid the formation of CSC-
favourable niches (i.e, hypoxia). In other words, cell numbers should be
maintained at low enough levels to prevent the formation of steep oxygen
gradients.

##  5 Dynamic simulations: switching between basins of attraction

Before introducing treatment with radiotherapy into the model, we pause to
investigate how the initial composition of the tumour affects its long time
behaviour. As in Fig. 12, we fix $V_{+}=8\times 10^{-4}$ and consider two
values of the oxygen consumption: $\gamma^{*}=2.0$ and $\gamma^{*}=4.0$.
Recall that we impose initial conditions in which the cell distribution in
uniform in space and normally distributed along the stemness axis (see Eq.
(13)). In what follows, we fix the variance $\sigma_{s}^{2}$ of the Gaussian
so that $\sigma_{s}^{2}=0.02$ and investigate how the long time behaviour
changes as the initial tumour burden $M_{0}$ and the mean phenotype $s_{0}$
vary. We restrict attention to $s_{0}\in[0.25,0.75]$ so that boundary effects
do not distort the Gaussian profile.

(a) $\gamma^{*}=2.0$

(b) $\gamma^{*}=4.0$

Figure 13: Series of plots showing how the long time behaviour of the tumour
depends on its initial composition for the cases of (a) low and (b) high
$\gamma^{*}$. We plot the evolution of the tumour burden as predicted by Eqs.
(4) for the initial condition (13) (see (a.1) and (b.1)). Each curve in (a.1)
(or (b.1)) corresponds to a point in the $(M_{0},s_{0})$ diagram (a.2) (or
(b.2)); we label with red crosses initial conditions that result in tumour
extinction and with blue dots the ones corresponding to tumour invasion. In
the subplots (a.3) and (b.3), we indicate how the $(M_{0},s_{0})$ plane
partitions into region that map to tumour extinction (red areas) and those
that map to tumour invasion (blue areas).

In Fig. 13, we present results from dynamic simulations. Each red curve in
Fig. 13(a.1) or 13(b.1) depicts tumour extinction and corresponds to a
specific cross in $(s_{0},M_{0})$ plane in Figs. 13(a.2) or 13(b.2); the blue
curves correspond to “tumour-invasion” and are denoted by blue circles in
Figs. 13(a.2) or 13(b.2). From Fig. 13(a) (i.e., $\gamma^{*}=2.0$), we note
that, for most initial conditions considered, the tumour is driven to
extinction. Focusing on Fig. 13(a.1), we see that the tumour burden $M$ does
not decay monotonically; rather, it grows initially, to relatively large
values ($M\approx 0.75$) in some cases, and only starts to decay at longer
times. Focusing on Fig. 13(b.1), we observe again non-monotonic evolution of
$M$, and the decay to extinction (see red curves) is delayed, as in Fig.
13(a.1). However, the maximum tumour burden prior to the start of the decay
phase for the red curves is much smaller than those observed in Fig. 13(a.1).
When considering the $(s_{0},M_{0})$ diagrams (Fig. 13(a.2) and 13(b.2)), we
see that “tumour-invasion” scenarios are concentrated in the upper left
corner. This suggests that only sufficiently large tumours, comprised of more
stem-like cells ($s_{0}<0.5$), can persist. This effect is more pronounced for
larger values of $\gamma^{*}$. In Fig. 13(b), tumour survival occurs for a
wider range of values of $s_{0}$ and $M_{0}$ than when $\gamma^{*}=2.0$ (Fig.
13(a)).

Based on the numerical simulations, we decompose the $(s_{0},M_{0})$ plane
into two distinct regions: one that maps to tumour-invasion (blue) and one
that maps to tumour-extinction (red). As shown in Fig. 13(a.3) and (b.3), the
boundary between these regions can be approximated by a straight line
$M_{0}=M_{cr}(s_{0})$ (note we are using a semi-logarithmic scale).
Accordingly, for any mean phenotype $s_{0}$, all tumours with initial mass
$M_{0}>M_{cr}(s_{0})$ will persist. Interestingly, the model predicts that
$M_{cr}$ increases with $s_{0}$, suggesting that as the proportion of stem-
like cells ($s_{0}\ll 1$) initially present increases, the smaller the value
of $M_{0}$ for which invasion occurs. On the other hand, as $s_{0}\rightarrow
1$, even large tumours may be eliminated. In this case, the tumour initially
comprises TDCs which are unable to drive tumour growth. This is consistent
with the idea that CSCs have the unique capacity to drive the formation and
maintenance of tumour colonies.

When applying treatment, the aim is to drive the tumour into the extinction
region. Therefore, we expect that it will be easier to eradicate a tumour with
a larger extinction region. From this point of view, comparing panels (a.3)
and (b.3), we see that the size of the extinction region increases with
$\gamma^{*}$, suggesting that tumours with higher consumption rates are more
resilient, even though they are smaller (see Fig. 12). This is because, as the
rate at which cells consume oxygen increases, hypoxic regions are initiated at
smaller tumour volumes (i.e., smaller numbers of tumour cells), creating an
environment in which CSCs can emerge and drive tumour growth (as in Fig. 12,
case C). On the other hand, for smaller values of $\gamma^{*}$, the threshold
tumour burden required for the formation of hypoxia is larger, so that even
large tumours may become extinct.

##  6 Dynamics in the presence of treatment

Figure 14: Schematic diagram showing how treatment is incorporated into the
model. Treatment is repeated over a period of $n_{rep}$ weeks so that the
final dosage delivered is $10\times n_{rep}$ Gy. We consider two different
protocols: ultra-hypo-fractionated radiotherapy (UHFRT), and fractionated
radiotherapy (FRT). For UHFRT, a single dose of 10 Gy is delivered on Monday
of each week during treatment; for FRT, 5 doses of 2 Gy are delivered on
Monday to Friday, with a break from treatment on Saturday and Sunday.

Having studied the model predictions in the absence of treatment, we now use
it to investigate tumour responses to radiotherapy. Due to the large number of
unknown model parameters, it is difficult to apply our model to patient
specific data and personalised radiation scheduling. However, it can provide
insight into the ways in which interactions between the different radio-
resistance mechanisms determine treatment responses and, in doing so, increase
understanding of the possible causes of relapse.

As in §5, we fix ($V_{+}=8\times 10^{-4}$), so we know that tumour extinction
is possible, and compare the response to treatment of tumours with different
oxygen consumption rates: $\gamma^{*}=2.0$ and $\gamma^{*}=4.0$. For each
value of $\gamma^{*}$, we consider as initial conditions the corresponding
non-trivial equilibrium solutions (i.e, steady state A or C in Fig. 12). As
shown in Fig. 14, we assume that the tumours receive 10 Gy of radiotherapy
every week for $n_{rep}$ weeks. Treatment is delivered either as a single,
high dose (Mon) or as 5 doses (Mon-Fri) with a two-day break (Sat-Sun). We
refer to the first treatment strategy as ultra-hypo-fractionated radiotherapy
(UHFRT) and the second one, which corresponds to standard-of-care, as
fractionated radiotherapy (FRT).

###  6.1 Numerical simulations

In Fig. 15, we compare the evolution of the tumour burden $M(t)$ for the two
treatment protocols. Fig. 15(a) shows the response of a tumour with a low rate
of oxygen consumption ($\gamma^{*}=2.0$), while Fig. 15(b) shows the response
of a tumour with a high oxygen consumption rate ($\gamma^{*}=4.0$). When
$\gamma^{*}=2.0$, a 5-week course of treatment is sufficient to eradicate the
tumour, regardless of the protocol used (see Fig. 15(a)). During the early
stages of treatment, UHFRT reduces the tumour burden more rapidly than FRT,
but after about 10 weeks, the dynamics of $M(t)$ are indistinguishable, with
the tumour burden decreasing monotonically to zero in both cases. By contrast,
Fig. 15(b) shows that when $\gamma^{*}=4.0$ treatment fails for both
protocols: while the tumour initially responds to treatment, it eventually
regrows to its original size. Interestingly, the tumour burdens at the end of
treatment, with either FRT or UHFRT, are similar for both values of
$\gamma^{*}$.

Figure 15: Evolution of the tumour burden, $M(t)$, during treatment for
different values of the parameter $\gamma^{*}$: (a) $\gamma^{*}=2.0$ and (b)
$\gamma^{*}=4.0$ . We compute $M(t)$ numerically by solving Eqs. (4)-(12) and
applying 50 Gy radiotherapy over a period of 5 weeks according to the two
protocols described in Fig. 14. The duration of treatment is highlighted by
the yellow area; for each panel, we zoom on this time frame to highlight
differences in the evolution of $M$. The remaining parameters are fixed at the
default values listed in Tables 1-2.

We now consider how the internal phenotypic composition of the tumours changes
during treatment, in order to determine whether there is any early indication
of different treatment outcomes. We also plot the minimum oxygen concentration
recorded in the tissue $\underline{c}(t)=\min_{X}(X,t)$. The results presented
in Fig. 16 suggest that during UHFRT the tumour composition is independent of
$\gamma^{*}$. In both cases, RT re-oxygenates the tumour (note that the
minimum oxygen level recorded in the tissue, $\underline{c}$, remains above
$c_{H}$ throughout treatment). Therefore, although $\Pi_{1}(t)$ increases
after each dose of radiotherapy (since RT selects for radio-resistant CSCs),
$\Pi_{1}$ decreases between doses as CSCs tend to differentiate when exposed
to sufficiently high oxygen levels. As a result, at $4$ weeks (i.e., just
after the last dose), $\Pi_{1}(t=4)<\Pi_{1}(t=0)\approx 0.25$. By contrast,
the fraction of differentiated (highly-proliferative) cells $\Pi_{2}(t)$
increases during treatment which explains the rapid regrowth of the tumour
when treatment ends. As the number of tumour cells increases, tissue oxygen
levels decrease. While for $\gamma^{*}=2.0$, $\underline{c}$ remains above the
threshold $c_{H}$, an hypoxic region forms at $t\approx 6$ weeks for
$\gamma^{*}=4.0$. Only at this time point, does the internal composition of
the two tumours start to diverge. While for $\gamma^{*}=2.0$, at long times
all DCs become TDCs and the tumour is eliminated, when $\gamma^{*}=4.0$, DCs
located in the hypoxic region de-differentiate, producing CSCs which enable
the tumour to survive and continue to grow.

Figure 16: Numerical results showing how, for the simulations presented in
Fig. 15, the tumour’s internal composition changes following 5 weeks of UHFRT
when (a) $\gamma^{*}=2.0$ and (b) $\gamma^{*}=4.0$. The top panels show the
evolution of the minimum oxygen level within the tumour,
$\underline{c}(t)=\min\limits_{X}c(X,t)$. The latter is compared to the oxygen
thresholds $c_{H}$ (orange horizontal line) and $c_{R}$ (red horizontal line).
In the lower panels, the white vertical line indicates the time at which the
last dose is administered.

Focusing now on FRT (see Fig. 17), we observe that the evolution of
$\underline{c}$ differs during the early stages of treatment for the two
values of $\gamma^{*}$. For $\gamma^{*}=2.0$, one cycle of FRT is sufficient
to raise oxygen levels above the hypoxic threshold $c_{H}$, while for
$\gamma^{*}=4.0$ two cycles are required. As a result of differences in the
oxygen levels, a larger fraction of CSCs are present in Fig. 17(b) during the
initial treatment phase. However, for both scenarios, the value of $\Pi_{1}$
at the end of treatment (here indicated by the vertical white line) is
negligible, and significantly smaller than the value recorded after the last
dose for UHFRT (see Fig. 16). However, FRT is less effective than UHFRT in re-
oxygenating the tumour because fewer cells are killed. For $\gamma^{*}=4.0$,
this results in the rapid formation of an hypoxic region after the end of
treatment so that a small fraction of CSCs persist in the tumour, in contrast
to the behaviour observed in Fig. 16(b). As mentioned above, in the case of
FRT, the dynamics of $\underline{c}$ changes markedly with $\gamma^{*}$, even
during the early phase of treatment, making tumour oxygen levels a possible
biomarker for predicting outcomes in response to FRT. This is in contrast to
UHFRT, where the dynamics of $\underline{c}$ (see Fig. 16) appears to be
independent of $\gamma^{*}$ during the treatment period so that it can not be
used as an early indicator of outcome in response to UHFRT.

Figure 17: Numerical results showing how, for the simulations presented in
Fig. 15, the tumour’s internal composition changes following 5 weeks of FRT
when (a) $\gamma^{*}=2.0$ and (b) $\gamma^{*}=4.0$.. In the top panels we plot
the minimum oxygen level $\underline{c}$ recorded in the tissue at any time
and compare it with the relevant oxygen threshold $c_{H}$ and $c_{R}$.

The simulation results presented in this section, indicate that oxygen plays a
key role in defining the response of a tumour to RT due in large part to its
(eventual) impact on the dynamics of CSCs. From this viewpoint, the
dimensionless oxygen consumption rate $\gamma^{*}$, which influences the
distribution of oxygen in the tumour, plays a key role in determining the
treatment outcome. Overall, we see that low values of $\gamma^{*}$ correlate
with better responses to both UHFRT and FRT. Our findings support the possible
benefit of combining RT and metabolic inhibitors, which can reduce oxygen
consumption by tumour cells by inhibiting mitochondrial activity, as recently
proposed by several experimental studies [4, 26, 52].

###  6.2 Alternative protocols: the possible benefit of waiting

The results presented in the previous section suggest that tumours with higher
oxygen consumption rates are more aggressive and might not respond well to
standard treatment. Here we investigate how RT scheduling might be adapted to
improve treatment efficacy. Classical approaches for RT scheduling look at the
total number of doses and the level of fractionation based on the concept of a
Biologically Effective Dose (BED) [9, 27]. Under the assumption that $\alpha$
and $\beta$ in the LQ model (11a) are constant parameters, the BED is regarded
as a measure of the total biological dose delivered to the tissue for a
certain dose fractionation schedule [41]. However, our model illustrates how
$\alpha$ and $\beta$ may change over time, as treatment affects the tumour’s
internal composition and its micro-environment (i.e, oxygen levels).
Consequently, the time at which a dose is administered can impact the overall
outcome and must therefore be taken into consideration when planning treatment
schedules. We therefore use our model to test whether waiting longer periods
between treatment cycles (for example, we might repeat treatment every two
weeks instead of every week) might be beneficial for more aggressive tumours
($\gamma^{*}=4$).

Figure 18: Comparison of different schedules for UHFRT where the interval
between successive doses and the total dose given are changed. (a) Summary of
the effect of different treatments (red cross: tumour goes extinct; blue
circle: treatment fails and the tumour persists); (b)-(e) series of plots
showing how the tumour burden $M(t)$ evolves for the different scenarios in
(a); while colour indicates the total dose, each subplot corresponds to a
different dose interval: (b) $1$ week; (c) $1.5$ weeks; (d) $2$ weeks; (e)
$2.5$ weeks.

The results for UHFRT are summarised in Fig. 18. Here the interval between
cycles of treatment corresponds to the interval between consecutive doses.
When the interval between doses is 1 week, as in Fig. 15, treatment fails,
even when we increase the total dose delivered to 80 Gy. By allowing more time
between doses, we observe better responses with a total dose of 70 Gy being
sufficient to prevent recurrence when there is an interval of 1.5 weeks
between doses. However, if the interval between doses is too long as in Fig.
18(e), cells have time to grow between doses so that, even after a total dose
of 80 Gy, the tumour burden $M$ is relatively large ($\approx 0.1$). Overall
we observe that the shorter the interval between doses, the smaller the tumour
burden at the end of treatment. Yet, recurrence is observed for inter-dose
interval of 1 week and not for intervals of 1.5 weeks (the latter is true for
70 and 80 Gy total doses). These results suggest that the different outcomes
cannot be (solely) understood on the basis of the overall tumour burden.

Figure 19: Plots showing how, for UHFRT, the evolution of the internal
phenotypic composition of the tumour changes as the period between successive
treatment cycles is increased from (a) 1 week to (b) 1.5 weeks. We apply a
total dose of 70 Gy and the vertical white line indicates the end of
treatment. On the top panel we plot the evolution of the minimum oxygen level
recorded in the tissue.

In Fig. 19, we show how the internal tumour composition changes over time when
a total of $70$ Gy is delivered waiting either 1 (a) or 1.5 (b) weeks between
doses. While tumour re-oxygenation is less effective when waiting longer
between doses, a dose interval of 1.5 weeks is more effective in reducing the
overall fraction of CSCs. This is evident when comparing the values of
$\Pi_{1}$ for the two protocols at the end of treatment (indicated by the
vertical white line). Given that a single dose of 10 Gy is sufficient to
reoxygenate the tissue, waiting longer between doses gives more time for CSCs
to differentiate into radio-sensitive DCs. However, if the time between doses
is too long (such as 2.5 weeks), the regrowth of the tumour between doses
would result in the formation of an hypoxic _niche_ in which CSCs would start
to accumulate resulting in increased resistance of the tissue. This explains
the non-monotonic dependence of treatment outcomes on inter-dose period
reported in Fig. 18.

When we increase the interval between treatment cycles for FRT, rather than
UHFRT, our model predicts a different trend. In this case, the interval
between cycles corresponds to the time between the administration of the first
dose of the five $2$ Gy doses, which are still delivered over 5 consecutive
days. As shown in Fig. 20(a), the efficacy of treatment decreases
monotonically with the inter-cycle interval. While for a standard interval of
1 week a total of 100 Gy is effective in preventing relapse, other schedules
fail. Considering the growth curves in Fig. 20(d) (corresponding to a 2 weeks
interval between cycles), we see that, by waiting longer, the tumour acquires
resistance to treatment. Indeed, the value of $M$ after the last dose of each
cycle increases which indicates strong selection for resistant phenotypes.

Figure 20: Comparison of different schedules for FRT where we change the
interval between successive treatment cycles and the total dose given over the
treatment period. (a) Summary of the effect of different treatments (red
cross: tumour goes extinct; blue circle: treatment fails and tumour persists);
(b)-(e) series of plots showing how the tumour burden $M(t)$ evolves for the
different scenarios in (a); while colour indicates the total dose, each
subplot corresponds to a different interval: (b) $1$ week; (c) $1.5$ weeks;
(d) $2$ weeks. The insets in panel (d) highlight that cycles of treatment
become less effective at diminishing the tumour burden due to the accumulation
of radio-resistance cells (see black arrow interpolating the tumour burden at
the end of each treatment cycle).

This selection is evident when we look at the evolution of the internal tumour
composition n Fig. 21. For an interval of 1 week, oxygen levels in the tumour
increase during treatment and the lack of an hypoxic niche allows CSCs to de-
differentiation. As a consequence, CSCs are negligible at the end of
treatment. On the contrary, for intervals of two weeks, the treatment fails to
re-oxygenate the tissue, the hypoxic _niche_ persists and oxygen levels remain
below $c_{R}$. As a result, large numbers of CSCs accumulate in the tissue,
well above their equilibrium value ($\approx 25\%$). Consequently, subsequent
rounds of RT become less effective in killing cells. Once treatment ends, the
large number of CSCs in the tumour enables it to quickly evolve to a non-zero
steady state.

Figure 21: Plots showing how, for FRT, the evolution of the internal
phenotypic composition of the tumour changes as the period between successive
treatment cycles is increased from (a) 1 week to (b) 2 weeks. We apply a total
dose of 100 Gy and the vertical white line indicates the end of treatment. On
the top panel we plot the evolution of the minimum oxygen level recorded in
the tissue.

In summary, our model suggests that the optimal schedules for UHFRT and FRT
may differ markedly due to the different ways in which the composition of the
tumours evolve during the two treatment protocols. For UHFRT, if doses are too
close to each other, then treatment selects for CSCs which drive subsequent
disease relapse. If, instead, the doses are too far apart, then the tumour has
time to regrow, and to generate a severely hypoxic region, which reduces the
overall treatment efficacy and drives treatment failure. By choosing an
intermediate interval between doses (i.e., 1.5 weeks), the model predicts
successful extinction of the tumour for physically realistic protocols (with
total RT doses of 70 Gy). The absence of an hypoxic _niche_ in this case, in
combination with the conversion of potential CSCs into radio-sensitive DCs,
enables the radiation to eradicate the tumour. By contrast, FRT becomes less
effective as the time between treatment cycles increases. This is because,
tumour re-oxygenation in FRT is less pronounced than in UHFRT, so that hypoxic
_niches_ form more rapidly after the end of treatment. Consequently shorter
intervals between treatment cycles are predicted to be more effective.
Nevertheless, this approach leads to the early appearance of a larger fraction
of CSCs and, thereafter, (under the conditions considered here) to the
survival of the tumour.

##  7 Conclusion

We have proposed and analysed a new mathematical model, structured by
phenotype and space, to describe the evolution of intra-tumour heterogeneity
in the presence and absence of radiotherapy. We have focused on the role
played by tissue oxygen levels in shaping the tumour’s internal phenotypic
composition, with hypoxia (i.e., abnormally low oxygen levels) being a key
driver of cell de-differentiation towards an aggressive, stem-like phenotype.
The limited susceptibility of the CSCs to treatment with radiotherapy is found
to be a key factor in driving the tumour’s survival.

After introducing our model in §2, we analysed its long time behaviour in the
absence of treatment, combining analytical and numerical techniques to
investigate its bifurcation structure as two model parameters vary: the
magnitude of the differentiation velocity $V_{+}$ and the dimensionless oxygen
(and, by extension, nutrient) consumption rate $\gamma^{*}$ which measures the
ratio of the timescales for oxygen diffusion and consumption i the tissue. We
find that for biologically reasonable values of these parameters the system
possesses two non-negative steady states which we term “tumour-extinction” and
“tumour-invasion”. Based on extensive numerical simulations, we conclude that
relatively small tumours characterised by large values of $\gamma^{*}$ can
successfully invade the tissue, while larger tumours with lower oxygen
consumption rates $\gamma^{*}$ may eventually be eliminated. This suggests
that the size of the basin of attraction of the “tumour-extinction” steady
states decreases with $\gamma^{*}$. This result can be better understood by
focusing on the internal composition of a tumour and establishing whether it
contains sufficiently large hypoxic regions. Hypoxia niches are essential for
the accumulation and persistence of cancer stem cells, which, due to their
high clonogenic potential, are responsible for the growth and successful
invasion of the tumour.

In §6, we have used the model to study tumour responses to radiotherapy (RT).
We find that large values of $\gamma^{*}$ correlates with more aggressive
tumours that respond poorly to treatment. This is mainly due to the difficulty
in reoxygenating the tumour rather than the presence, prior to treatment, of
radio-resistant CSCs. In particular, our model predicts that treatment will be
successful only if the tumour burden at the end of treatment is low enough to
prevent the formation of a durable hypoxic region after the end of treatment.
However, as shown in §5, the maximum tumour burden that prevents relapse
depends on the internal phenotypic composition of the tumour at the end of
treatment and is, therefore, case (or patient) specific. When comparing
different treatment protocols, we find that larger doses administered less
often can be more effective than standard of care fractionation (2 Gy for 5
days per week). In §6.2, we further investigate how different scheduling of
the same dose fractionation protocol can impact the outcome of treatment. We
find that, in certain cases, “sooner” is not always “better” and waiting
longer between treatment cycles can be beneficial. This is because RT selects
for radio-resistant phenotypes and therefore, if doses are too close to each
other, resistant cells accumulate until eventually the tumour stops responding
to treatment. At the same time, waiting too long can also be harmful, as the
tumour has the time to grow sufficiently large and to develop extremely
hypoxic regions where radio-resistant CSCs accumulates and all tumour cells
are less radio-sensitive due to the lack of oxygen.

A natural question is whether our conclusions generalise to more realistic,
higher-dimensional spatial geometries, e.g., the framework of tumour spheroids
[12]. When symmetry arguments are used to reduce 3D-problems to 1D geometries,
we find that the qualitative behaviour of the solution is similar to that
observed for the 1D-Cartesian problem investigated here (results not shown).
However, fully 2D and 3D simulations may exhibit additional behaviours not
captured by the simplified 1D-geometry.

An interesting prediction of our model is that the rate of oxygen consumption
by the cells plays a key role in determining treatment (here, radiotherapy)
outcomes. This emphasises the possible benefit of combining RT and metabolic
inhibitors to improve response to treatment, as recently proposed by several
studies [4, 26, 52]. By decreasing the rate at which the cells consume oxygen,
oxygen levels in the tumour can be increased to prevent the formation of
regions of radio-biological hypoxia and increase the tumour’s overall
sensitivity to radiotherapy. Furthermore, our model suggests that metabolic
inhibitors can also prevent the emergence of CSCs, by re-oxygenation of large
hypoxic _niches_ in which stem-like cells initially might accumulate. A
natural extension of the model would be to introduce the effect of metabolic
inhibitors and investigate how they should be combined with RT to optimise
treatment outcomes. Furthermore, when considering the use of high RT doses
($\geq 10$ Gy) in a clinical setting, RT side-effects should also be included
into the model.

While our model captures the higher radio-resistance of CSCs, recent findings
suggest that radio-therapy can directly impact the differentiation status of
cells [77]. In particular, RT-induced de-differentiation has been observed _in
vitro_ and linked to the accumulation of senescent-type cells that can support
the _de novo_ formation of CSCs by influencing the tumour micro-environment.
Once more information about the mechanism linking senescence and stemness
become available, they could be easily incorporated into our model by for
example allowing the phenotypic advection velocity itself (i.e., $v_{s}$) to
depend on other environmental factors besides oxygen. Such additions would
affect the evolution of the tumour heterogeneity during treatment, shifting
the system into the basin of attraction of the “tumour-invasion” steady state
solution. This raises the interesting question of what treatment strategies
may be more effective when RT-induced radio-resistance is included and whether
these may differ from the discussion herein.

Ultimately, one can envision enlarging the cells’ phenotypic state space by
introducing additional “synthetic” axes that can capture more fully the cancer
cells heterogeneity and evolution. However, the complexity of the model would
make it difficult to gain insight into the role of the different “phenotypic”
variables. Hence the need to start analysing simpler model, where a subset of
these “synthetic” axes are included. From this point of view, given the amount
of experimental evidence linking metabolism, stemness and treatment outcomes,
a natural next step would be to include in the model two structure variables
accounting for the stemness and metabolic status of cells. While previous
theoretical studies have focused on the ways in which nutrient levels and/or
treatment shape metabolic heterogeneity in tumours (e.g, [2], [39],[75]), the
role of stemness is yet to be investigated. A model that captures all of these
aspects could provide insight into the evolution of complex tumour ecosystem
and help in the design of effective treatment strategies. While these goals
are intriguing from a theoretical viewpoint, a different, but equally
important, direction involves attempting to apply the model to experimental
data. The development of RNA-seq techniques has resulted in a wealth of data
being collected on the evolution of cells phenotypes (see, e.g., [56]). The
question then is how these rich datasets can be used to inform phenotypic-
structured models. From this point of view, important practical aspects have
to be addressed, such as parameter identifiability and estimation, that can be
challenging when considering highly-parameterised and computationally-
expensive models. Addressing these limitations will be crucial for translating
our theoretical study into clinically-useful insight.

## Funding

The authors wish to thank Prof Christos E. Zois for the helpful and
illuminating discussions on the biology underpinning the model. G.L.C. is
supported by from EPSRC and MRC Centre for Doctoral Training in Systems
Approaches to Biomedical Science and Cancer Research UK. G.L.C. acknowledges
also the support from Worcester College via the Martin Senior Scholarships.
P.G.K. acknowledges support from the Leverhulme Trust via a Visiting
Fellowship and thanks the Mathematical Institute of the University of Oxford
for its hospitality during the early stages of this work.

##  Appendix A Numerical methods for dynamical simulations

We discretise the spatial variable $X$ using standard central finite
difference on a stensil of $m_{x}+1$ points, where the Dirichlet condition in
Eq. (12) is strongly imposed, while ghost points are used to enforce the
Neumann boundary conditions. For the phenotypic axis $s$, we use instead
finite volume where we divide the axis into $m_{s}$ cells. Following [32], we
control the advection component using a Koren limiter. As discussed in [14],
the use of a finite volume scheme allows us to properly account for the
advective structural flux. After having discretised both the $X$ and the $s$
axis, we obtain from Eq. (12) a set of $m_{x}+1$ non-linear algebraic
constraints ($\boldsymbol{g}(\boldsymbol{n}(t),\boldsymbol{c}(t))=0$) coupled
to a system of $(m_{x}+1)\times m_{s}$ non-linear odes of the form
$\dot{\boldsymbol{n}}(t)=\boldsymbol{f}(\boldsymbol{c}(t),\boldsymbol{n}(t),t)$
derived from Eq. (4), where $\boldsymbol{n}(t)\in\mathcal{R}^{(m_{x}+1)m_{s}}$
and $\boldsymbol{c}(t)\in\mathcal{R}^{(m_{x}+1)}$. Note that, given
$\boldsymbol{n}$, the implicit function $\boldsymbol{g}$ uniquely defines
$\boldsymbol{c}$, i.e., we can write
$\boldsymbol{c}=\boldsymbol{c}(\boldsymbol{n})$. We define
$\tilde{\boldsymbol{f}}(\boldsymbol{n},t)=\boldsymbol{f}(\boldsymbol{c}(\boldsymbol{n}(t)),\boldsymbol{n}(t),t)$,
so that the evolution of the vector function $\boldsymbol{n}(t)$ can be
written as
$\dot{\boldsymbol{n}}(t)=\tilde{\boldsymbol{f}}(\boldsymbol{n}(t),t)$. To
advance in time, we use an implicit multistep _bdf15_ method for stiff
equation from the _scipy_ library in python. Practically the input function
$\tilde{\boldsymbol{f}}$ to the _bdf15_ method reads

Function  _$\tilde{\boldsymbol{f}}$(_$\boldsymbol{n}$ , $t$_)_:

Solve for $\boldsymbol{c}$ using fix point iterations

Assemble the vector $\tilde{\boldsymbol{f}}$ using the computed oxygen
distribution $\boldsymbol{c}$

return

When including treatment, we want to properly capture the fast time scale of
RT delivery. To do so, at each time $t_{i}$ at which a dose of RT is
delivered, we integrate the solution ${\boldsymbol{n}}$ up to time $t_{i}$
using _bdf15_ method (as described above). In the interval
$[t_{i},t_{i}+\delta T]$, with $\delta T=0.5$, we instead use a simple
explicit Euler scheme with small fixed time step $dt=5\times 10^{-3}$; we then
switch back to _bdf15_ method up to time $t_{i+1}$, when the next radiation
dose is delivered.

##  Appendix B Parameter values

We summarise here the parameter values used in the simulation of Eqs. (4)-(12)
presented in the main text. In Table 1, we list the value of dimensional
parameter values that could be estimated from the literature. The values of
non-dimensional parameters are instead given in Table 2.

Table 1: Summary of typical dimensional values for some of the parameters in Eq. (2), togehter with supporting references. Parameter | Reference value(s) | Refs
---|---|---
$D_{xn}$ |  $10^{-10}-10^{-8}$ [mm2/s] | [70, 75]
$p^{max}_{0}$, $p^{max}_{1}$ | 0.005, 0.02 [1/hr] | [14, 71]
OER | 3 | [43]
$\alpha_{min}$, $\alpha_{max}$ | 0.005, 0.1 [1/Gy] | [61]
$\beta_{min}$, $\beta_{max}$ | 0.002, 0.01 [1/Gy2] | [61]
$\gamma_{0}C_{\infty}$ |  $5–50$ [mmHg/s] | [35]
$D_{xc}$ |  $2\times 10^{-3}$ [mm2/s] | [35]
$C_{\infty}$ |  $50$ [mmHg] | [50, 79]
$C_{H}$ |  $15$ [mmHg] | [50, 79]
$C_{R}$ |  $5$ [mmHg] | [50, 79]
$C_{N}$ |  $\approx 0$ [mmHg] | [79]
$L$ |  $\approx 200$ [$\mu$m] |
Table 2: Summary of non-dimensional parameters groups that appears in Eq. (4). Parameter values are consistent with those reported in Table 1 and/or in [14]. Parameter | Value(s)
---|---
$\mathcal{D}_{sn}$ |  $5\times 10^{-6}$ [1/hr]
$V_{+}$, $V_{-}$ |  $8$, $2$ [$\times 10^{-4}$/hr]
$\xi_{+}$, $\xi_{-}$ |  $0.1$, $0.1$
$\omega_{+}$, $\omega_{-}$ | 1, 2
$\mathcal{D}_{xn}$ |  $10^{-4}$ [1/hr]
$K_{0}$, $K_{1}$ |  $0.05$, $0.3$
$g_{0}$, $g_{1}$ | 0.01, 0.04
$d_{f}$ |  $0.001$ [1/hr]
$k_{f}$ | $10$
$\xi_{R}$ | 0.2
$\gamma^{*}$ |  $2-20$ [1/hr]
$\mathcal{D}_{xc}$ |  $180$ [1/hr]
$c_{\infty}$ | $1.0$
$c_{H}$ | 0.3
$c_{R}$ | 0.1
$c_{N}$ | 0.0125

##  Appendix C Bifurcation Analysis

###  C.1 Linear stability analysis of the trivial steady state

In this section, we determine the linear stability of the trivial steady state
by seeking solutions of Eq. (4)-(5) of the form:

$n=\delta ne^{\lambda t}+O(\delta^{2}),\quad c=c_{\infty}+\delta ce^{\lambda
t}+O(\delta^{2}),$ (17)

where $0<\delta c,\delta n\ll 1$. Substituting (17) into Eqs. (4)-(5) and
retaining terms which are linear in $\delta$ we obtain the following
linearised problem:

$\displaystyle\lambda\delta n=\mathcal{D}_{xn}\frac{\partial^{2}\delta
n}{\partial X^{2}}+\frac{\partial}{\partial
s}\left(\mathcal{D}_{sn}\frac{\partial\delta n}{\partial s}-\delta
nv_{s}(s,c_{\infty})\right)+\mathcal{F}(s,c_{\infty},0)\delta n,$ (18a)
$\displaystyle\lambda\delta c=\mathcal{D}_{xc}\frac{\partial^{2}\delta
c}{\partial X^{2}}-\gamma\mathcal{D}_{xc}\delta\phi,$ (18b)
$\displaystyle\mathcal{D}_{xs}\frac{\partial\delta n}{\partial s}-\delta
nv_{s}=0,\qquad s\in\left\\{0,1\right\\},\,X\in[0,1],\,t>0,$ (18c)
$\displaystyle\left.\frac{\partial\delta c}{\partial
X}\right|_{X=0}=0,\quad\delta c(1,t)=0,\quad t>0,$ (18d)
$\displaystyle\left.\frac{\partial\delta n}{\partial
x}\right|_{X=0}=\left.\frac{\partial\delta n}{\partial x}\right|_{X=1}=0,\quad
s\in(0,1),\,t>0.$ (18e)

where $\delta\phi(X)=\int_{0}^{1}\delta n(X,s)ds$. We note that Eq. (18b) for
$\delta c$ decouples from Eq. (18a) for $\delta n$. We can therefore solve Eq.
(18a) independently of $\delta c$, and seek a separable solution of the form:

$\delta n=\Lambda(x)W(s)$ (19)

where $\Lambda$ and $W$ needs to satisfy:

$\displaystyle\begin{aligned}
\lambda=\frac{\mathcal{D}_{xn}}{\Lambda}\frac{d^{2}\Lambda}{dX^{2}}+\frac{1}{W}\frac{d}{ds}\left(\mathcal{D}_{sn}\frac{dW}{ds}-v_{s}W\right)+(p(s,c_{\infty})-f(s)),\end{aligned}$
(20a) $\displaystyle\mathcal{D}_{sn}\frac{dW}{ds}=0,\qquad
s\in\left\\{0,1\right\\},$ (20b)
$\displaystyle\left.\frac{d\Lambda}{dX}\right|_{X=0}=\left.\frac{d\Lambda}{dX}\right|_{X=1}=0,$
(20c)

where in writing Eq. (20b) we have used the fact that the velocity $v_{s}$
vanishes at the boundaries.

Using standard arguments, we find that the eigenfunctions $\Lambda_{m}$
associate to the eigenvalue $\omega_{m}$ are given by

$\displaystyle\Lambda_{m}=\cos\left(\omega_{m}X\right),\quad\omega_{m}=\pi
m,\quad m=0,1,2,\ldots.$ (21a) so that the eigenfunctions $W$ associated to
the phenotypic axis, $s$, must satisfy the following ODE:
$\displaystyle\left(\omega_{m}^{2}\mathcal{D}_{xn}+\lambda\right)W=\frac{d}{ds}\left(\mathcal{D}_{sn}\frac{dW}{ds}-v_{s}W\right)+\left[p(s,c_{\infty})-f(s)\right]W\equiv\mathcal{M}(W),$
(21b) $\displaystyle\mathcal{D}_{sn}\frac{dW}{ds}=0,\qquad
s\in\left\\{0,1\right\\}.$ (21c)

Here, the operator $\mathcal{M}$ is equivalent to the differential operator
obtained by linearising the well-mixed formulation of Eqs. (4), which is
obtained by setting all spatial derivatives to zero and considering an
homogeneous oxygen concentration. The spectrum of the operator $\mathcal{M}$
has been investigated in more detail in [14]. While we have no explicit
solution for the eigenvalues $\kappa_{k}$ and eigenfunctions $W_{k}$ of
$\mathcal{M}$, we can show that $\kappa_{k}$ are real and bounded from above.
If $\mathcal{D}_{xn}$ is constant and $W_{k}$ denotes the eigenfunctions of
$\mathcal{M}$ associated to the eigenvalue $\kappa_{k}$, then $\delta
n_{m,k}=\Lambda_{m}(X)W_{k}(s)$ is an eigenfunction of the operator, Eq.
(20a), with real eigenvalue
$\lambda_{m,k}=\kappa_{k}-\omega_{m}^{2}\mathcal{D}_{xn}$. Since
$\omega^{2}\mathcal{D}_{xn}\geq 0$, we find that $\lambda_{m,k}\leq\kappa_{k}$
and $\lambda_{m,k}<\lambda_{m^{\prime},k}$ if $m^{\prime}<m$. Consequently,
$\lambda_{m,k}<0$ for all $k$ and $m$ if and only if
$\bar{\kappa}=\max_{k}\kappa_{k}<0$.

The question now is whether $\lambda=\bar{\kappa}$ is an eigenvalue for the
full system Eqs. (18). If this is the case, then the stability of the trivial
steady-state for the model Eqs. (4)-(5) boils down to studying the linear
stability of the steady state for the well-mixed model for homogeneous oxygen
concentration $c_{\infty}$.

Let us consider the case $\omega_{m}\equiv 0$ for which wlog
$\Lambda_{0}\equiv 1$. In this case, the eigenfunction associated to the
eigenvalue $\bar{\kappa}$ is $\delta n=W_{\bar{\kappa}}(s)$. The eigenfunction
is defined up to a multiplicative constant which we set so that
$\int_{0}^{1}W_{k}(s)=\delta\phi\equiv\Delta>0$. In order for $\bar{\kappa}$
to be an eigenvalue of Eqs. (18), Eq. (18b) must have a solution for
$\lambda=\bar{\kappa}$ and $\delta\phi\equiv\Delta$. Substituting for
$\lambda$ and $\delta\phi$, we obtain:

$\displaystyle\bar{\kappa}\delta
c+\gamma^{*}\mathcal{D}_{xc}\Delta=\mathcal{D}_{xc}\frac{\partial^{2}\delta
c}{\partial X^{2}},$ (22a) $\displaystyle\delta c(1)=0,\quad\partial_{x}\delta
c(0)=0,$ (22b) which has the solution: $\displaystyle\delta
c=\bar{C}(X)=\begin{cases}\frac{\gamma^{*}\mathcal{D}_{xc}\Delta}{\bar{\kappa}}\left[\frac{\cosh(\bar{\tau}X)}{\cosh(\bar{\tau})}-1\right],&\quad\bar{\kappa}>0,\\\
\frac{\gamma^{*}\Delta}{2}(X^{2}-2),&\quad\bar{\kappa}=0,\\\
\frac{\gamma^{*}\mathcal{D}_{xc}\Delta}{\bar{\kappa}}\left[\frac{\cos(\bar{\tau}X)}{\cos(\bar{\tau})}-1\right],&\quad\bar{\kappa}<0,\end{cases}$
(22c) where $\bar{\tau}=\sqrt{|\bar{\kappa}/\mathcal{D}_{xn}|}$.

Given that a solution exists, we have that $\bar{\kappa}$ belongs to the
spectrum of the full operator, Eqs. (18), with eigenfunction
$(\bar{W}(s),\bar{C}(X))$. Further, since the spectrum of Eqs. (18) is a
subset of $\lambda_{m,k}\leq\bar{\kappa}$, $\bar{\kappa}$ will be the
eigenvalue with largest real part and therefore dictates the stability of the
trivial solution. From [14], we know that $\bar{\kappa}$ decreases from
positive to negative values as the magnitude of the advection velocity $V_{+}$
increases. The critical value $V_{+}^{(cr)}$ at which $\bar{\kappa}=0$,
corresponds to a transcritical bifurcation. This bifurcation must therefore
occur also in the spatial model, and, based on the above result, the location
of the transcritical bifurcation $V_{+}=V_{+}^{cr}$ is independent of
$\gamma^{*}$, $\mathcal{D}_{xn}$ and $\mathcal{D}_{xc}$. We conclude therefore
that the critical value $V_{+}^{(cr)}$ at which the trans-critical bifurcation
occurs will not been affected by the inclusion of spatial effects, in line
with our numerical observations in the text.

###  C.2 Numerical continuation and stability estimates

In order to compute the equilibrium states $(n_{\infty}(x,s),c_{\infty}(x))$
of the full system (4)-(5) in the absence of treatment, we simply set to zero
all time derivatives to obtain:

$\displaystyle E_{n}(n_{\infty},c_{\infty})=0,\quad
E_{c}(n_{\infty},c_{\infty})=0,\quad\boldsymbol{E}_{bc}(n_{\infty},c_{\infty})=\boldsymbol{e}_{6},$
(23a) where $\boldsymbol{e}_{6}=[0,0,0,0,0,1]^{T}$ and $\displaystyle
E_{n}(n,c)$ $\displaystyle=\mathcal{D}_{xn}\frac{\partial^{2}n}{\partial
X^{2}}+\frac{\partial}{\partial s}\left(\mathcal{D}_{sn}\frac{\partial
n}{\partial s}-v_{s}(s,c)n\right)+\mathcal{F}(s,c,\phi)n,$ (23b)
$\displaystyle E_{c}(n,c)$
$\displaystyle=\mathcal{D}_{xc}\frac{\partial^{2}c}{\partial
X^{2}}-\gamma^{*}\mathcal{D}_{xc}\phi H_{\epsilon_{\gamma}}(c-c_{N}),$ (23c)
$\displaystyle\boldsymbol{E}_{bc}(n,c)$
$\displaystyle=\left[\left.\frac{\partial n}{\partial
s}\right|_{s=0},\left.\frac{\partial n}{\partial
s}\right|_{s=1},\left.\frac{\partial n}{\partial
X}\right|_{X=0},\left.\frac{\partial n}{\partial
X}\right|_{X=1},\left.\frac{\partial c}{\partial
X}\right|_{X=0},\left.c\right|_{X=1}\right].$ (23d)

Due to the non linearities in the system, we cannot solve for Eqs. (23)
explicitly. We instead rely on numerical approaches in particular continuation
techniques to investigate the system equilibria and their stability. Let us
consider a perturbation of the equilibrium points:

$n=n_{\infty}+e^{\lambda t}\delta n,\quad c=c_{\infty}+e^{\lambda t}\delta c,$
(24)

where the perturbations are small, i.e. $|\delta c|,|\delta n|\ll 1$, and the
growth rate $\lambda\in\mathbb{C}$ is a complex number.

Substituting Eqs. (24) into the system (4)-(5), we obtain:
$\displaystyle\lambda\begin{bmatrix}\delta n\\\ \delta c\end{bmatrix}$
$\displaystyle=J(n_{\infty},c_{\infty})\begin{bmatrix}\delta n\\\ \delta
c\end{bmatrix},$ (25a) $\displaystyle J(n_{\infty},c_{\infty})$
$\displaystyle=\begin{bmatrix}\frac{\delta E_{n}}{\delta
n}(n_{\infty},c_{\infty})&\frac{\delta E_{n}}{\delta
c}(n_{\infty},c_{\infty})\\\\[5.0pt] \frac{\delta E_{c}}{\delta
n}(n_{\infty},c_{\infty})&\frac{\delta E_{c}}{\delta
c}(n_{\infty},c_{\infty})\end{bmatrix}$ (25b) with the constraint:
$\displaystyle\boldsymbol{E}_{bc}(\delta n,\delta c)=\boldsymbol{0}.$ (25c)
The stability of the steady state can therefore be related to the spectrum of
the Jacobian $J$ in particular $(n_{\infty},c_{\infty})$ will be stable if
$\Re(\lambda)<0$ for all $\lambda$ such that Eq. (25) admits a nontrivial
solution.

We solve the nonlinear problem (23) numerically via parameter continuation
techniques using Julia’s package BifurcationKit [74], which allows us to draw
the bifurcation diagram for the system as a function of a given parameter. To
do so we discretise Eqs. (23) as described in A and we compute the discretised
approximation of $E_{n},E_{c}$ and $J$. By computing the eigenvalues of the
matrix $J$ we can approximate the eigenvalues of the differential operator
(25) and study the stability of the computed steady states. To characterise
the solution we define the total number of cells in the domain:

$M_{\infty}=\int_{0}^{1}\phi_{\infty}(X)dX.$ (26)

###  C.3 Fold continuation

Figure 22: Continuation of the fold points $\mbox{SN}_{1}$ and $\mbox{SN}_{2}$
in the ($\gamma^{*}$,$V_{+}$) parameter space. Here $\mbox{SN}_{1}$ and
$\mbox{SN}_{2}$ are as identified in the bifurcation diagrams in Figs. 10 and
11. The points $\mbox{SN}_{1}$ and $\mbox{SN}_{2}$ divide the parameter spaces
in regions where Eqs. (23) admits a different number of solutions (as
indicated by the different colours). The two saddle nodes originates from a
cusp point.

Using Julia’s package BifurcationKit [74], we are also able to compute the
location of saddle-node bifurcation points and follow their location as we
vary the values of model parameters. This, however, is highly computationally
expensive as it augments the system (23), with an additional $(m_{x}+1)m_{s}$
constraints to define the fold points. We were therefore constrained in the
number of points $m_{x}$ and $m_{s}$ that we could use to discretise the
system. This implies that the results are not highly accurate but can still be
informative in understanding the bifurcation structure trends of the system
under study.

Looking at Fig. 22, we see that for sufficiently small $\gamma^{*}$, there are
no fold points (i.e., the bifurcation diagram is analogous to Fig. 9). The two
fold points $SN_{1}$ and $SN_{2}$ originates from a singular cusp (i.e., a
codimension two) point, and then depart from each other. While $SN_{2}$
steeply increase reaching large values of $V_{+}$, the trajectory of $SN_{1}$
asymptotes to a finite value of the advection velocity, $V_{+}^{(cr)}$, as
$\gamma^{*}\rightarrow\infty$. Here $V_{+}^{(cr)}$ corresponds to the value of
$V_{+}$ for which the system has a transcritical bifurcation (see Fig. 10). As
discussed in C.1, the location of $V_{+}^{(cr)}$ is independent of the
parameter $\gamma^{*}$ and divides the $(\gamma^{*},V_{+})$ plane into two
regions: the one where the trivial steady state is stable (see white dotted
region) and the one where it is unstable.

## References

  * Aponte and Caicedo [2017] Aponte, P.M., Caicedo, A., 2017\.  Stemness in cancer: Stem cells, cancer stem cells, and their microenvironment.  Stem Cells International 2017.  doi:10.1155/2017/5619472.
  * Ardaševa et al. [2020] Ardaševa, A., Gatenby, R.A., Anderson, A.R.A., Byrne, H.M., Maini, P.K., Lorenzi, T., 2020\.  A Mathematical Dissection of the Adaptation of Cell Populations to Fluctuating Oxygen Levels.  Bulletin of Mathematical Biology 2020 82:6 82, 1–24.  doi:10.1007/S11538-020-00754-7.
  * Arnold et al. [2020] Arnold, C.R., Mangesius, J., Skvortsova, I.I., Ganswindt, U., 2020\.  The Role of Cancer Stem Cells in Radiation Resistance.  Frontiers in Oncology 10, 164\.  doi:10.3389/FONC.2020.00164.
  * Ashton et al. [2016] Ashton, T.M., Fokas, E., Kunz-Schughart, L.A., Folkes, L.K., Anbalagan, S., Huether, M., Kelly, C.J., Pirovano, G., Buffa, F.M., Hammond, E.M., Stratford, M., Muschel, R.J., Higgins, G.S., McKenna, W.G., 2016\.  The anti-malarial atovaquone increases radiosensitivity by alleviating tumour hypoxia.  Nature Communications 7, 1–13.  doi:10.1038/ncomms12308.
  * Astanin and Preziosi [2009] Astanin, S., Preziosi, L., 2009\.  Mathematical modelling of the Warburg effect in tumour cords.  Journal of Theoretical Biology 258, 578–590.  doi:10.1016/J.JTBI.2009.01.034.
  * Axelson et al. [2005] Axelson, H., Fredlund, E., Ovenberger, M., Landberg, G., Påhlman, S., 2005.  Hypoxia-induced dedifferentiation of tumor cells – A mechanism behind heterogeneity and aggressiveness of solid tumors.  Seminars in Cell & Developmental Biology 16, 554–563.  doi:10.1016/J.SEMCDB.2005.03.007.
  * Bacevic et al. [2017] Bacevic, K., Noble, R., Soffar, A., Ammar, O.W., Boszonyik, B., Prieto, S., Vincent, C., Hochberg, M.E., Krasinska, L., Fisher, D., 2017\.  Spatial competition constrains resistance to targeted cancer therapy.  Nature Communications 2017 8:1 8, 1–15.  doi:10.1038/s41467-017-01516-1.
  * Baliu-Piqué et al. [2020] Baliu-Piqué, M., Pandiella, A., Ocana, A., 2020.  Breast Cancer Heterogeneity and Response to Novel Therapeutics.  Cancers 12, 1–16.  doi:10.3390/CANCERS12113271.
  * Barendsen [1982] Barendsen, G., 1982.  Dose fractionation, dose rate and iso-effect relationships for normal tissue responses.  International Journal of Radiation Oncology*Biology*Physics 8, 1981–1997.  doi:https://doi.org/10.1016/0360-3016(82)90459-X.
  * Barendsen et al. [2001] Barendsen, G.W., Van Bree, C., Franken, N.A., 2001.  Importance of cell proliferative state and potentially lethal damage repair on radiation effectiveness: Implications for combined tumor treatments (Review).  International Journal of Oncology 19, 247–256.  URL: http://www.spandidos-publications.com/10.3892/ijo.19.2.247/abstracthttps://www.spandidos-publications.com/10.3892/ijo.19.2.247, doi:10.3892/IJO.19.2.247.
  * Black and McGranahan [2021] Black, J.R.M., McGranahan, N., 2021\.  Genetic and non-genetic clonal diversity in cancer evolution.  Nature Reviews Cancer , 1–14doi:10.1038/s41568-021-00336-2.
  * Bull et al. [2020] Bull, J.A., Mech, F., Quaiser, T., Waters, S.L., Byrne, H.M., 2020.  Mathematical modelling reveals cellular dynamics within tumour spheroids.  PLOS Computational Biology 16, e1007961.  doi:10.1371/JOURNAL.PCBI.1007961.
  * Cajal et al. [2020] Cajal, S.R.y., Sesé, M., Capdevila, C., Aasen, T., Mattos-Arruda, L.D., Diaz-Cano, S.J., Hernández-Losa, J., Castellví, J., 2020.  Clinical implications of intratumor heterogeneity: challenges and opportunities.  Journal of Molecular Medicine (Berlin, Germany) 98, 161.  doi:10.1007/S00109-020-01874-2.
  * Celora et al. [2021] Celora, G., Byrne, H., Zois, C., Kevrekidis, P., 2021\.  Phenotypic variation modulates the growth dynamics and response to radiotherapy of solid tumours under normoxia and hypoxia.  Journal of Theoretical Biology , 110792doi:https://doi.org/10.1016/j.jtbi.2021.110792.
  * Chae and Kim [2018] Chae, Y.C., Kim, J.H., 2018\.  Cancer stem cell metabolism: target for cancer therapy.  BMB Reports 51, 319\.  doi:10.5483/BMBREP.2018.51.7.112.
  * Chisholm et al. [2016] Chisholm, R., Lorenzi, T., Clairambault, J., 2016.  Cell population heterogeneity and evolution towards drug resistance in cancer: Biological and mathematical assessment, theoretical treatment optimisation.  Biochimica et Biophysica Acta 1860, 2627–2645.  doi:10.1016/j.bbagen.2016.06.009.
  * Chisholm et al. [2015] Chisholm, R.H., Lorenzi, T., Lorz, A., et al., 2015\.  Emergence of Drug Tolerance in Cancer Cell Populations: An Evolutionary Outcome of Selection, Nongenetic Instability, and Stress-Induced Adaptation.  Cancer Research 75, 930–939.  doi:10.1158/0008-5472.CAN-14-2103.
  * Dagogo-Jack and Shaw [2018] Dagogo-Jack, I., Shaw, A.T., 2018\.  Tumour heterogeneity and resistance to cancer therapies.  Nature Reviews Clinical Oncology 15, 81–94.  doi:10.1038/nrclinonc.2017.166.
  * Damaghi et al. [2021] Damaghi, M., West, J., Robertson-Tessi, M., Xu, L., Ferrall-Fairbanks, M.C., Stewart, P.A., Persi, E., Fridley, B.L., Altrock, P.M., Gatenby, R.A., Sims, P.A., Anderson, A.R.A., Gillies, R.J., 2021.  The harsh microenvironment in early breast cancer selects for a Warburg phenotype.  Proceedings of the National Academy of Sciences 118\.  doi:10.1073/PNAS.2011342118.
  * Dhawan et al. [2016] Dhawan, A., Tonekaboni, S.A.M., Taube, J.H., Hu, S., Sphyris, N., Mani, S.A., Kohandel, M., 2016.  Mathematical modelling of phenotypic plasticity and conversion to a stem-cell state under hypoxia.  Scientific Reports 2016 6:1 6, 1–10.  doi:10.1038/srep18074.
  * Diehn et al. [2009] Diehn, M., Cho, R.W., Lobo, N.A., Kalisky, T., Dorie, M.J., Kulp, A.N., Qian, D., Lam, J.S., Ailles, L.E., Wong, M., Joshua, B., Kaplan, M.J., Wapnir, I., Dirbas, F.M., Somlo, G., Garberoglio, C., Paz, B., Shen, J., Lau, S.K., Quake, S.R., Brown, J.M., Weissman, I.L., Clarke, M.F., 2009\.  Association of reactive oxygen species levels and radioresistance in cancer stem cells.  Nature 2009 458:7239 458, 780–783.  URL: https://www.nature.com/articles/nature07733, doi:10.1038/nature07733.
  * Enderling et al. [2009] Enderling, H., Park, D., Hlatky, L., Hahnfeldt, P., 2009\.  The Importance of Spatial Distribution of Stemness and Proliferation State in Determining Tumor Radioresponse.  Mathematical Modelling of Natural Phenomena 4, 117–133.  URL: https://www.mmnp-journal.org/articles/mmnp/abs/2009/03/mmnp20093p117/mmnp20093p117.html, doi:10.1051/MMNP/20094305.
  * Fadai and Simpson [2020] Fadai, N.T., Simpson, M.J., 2020\.  Population Dynamics with Threshold Effects Give Rise to a Diverse Family of Allee Effects.  Bulletin of Mathematical Biology 82, 74.  doi:10.1007/s11538-020-00756-5.
  * Fanelli et al. [2020] Fanelli, G.N., Naccarato, A.G., Scatena, C., 2020.  Recent Advances in Cancer Plasticity: Cellular Mechanisms, Surveillance Strategies, and Therapeutic Optimization.  Frontiers in Oncology 10, 569\.  URL: https://www.frontiersin.org/article/10.3389/fonc.2020.00569/full, doi:10.3389/fonc.2020.00569.
  * Fiandaca et al. [2021] Fiandaca, G., Delitala, M., Lorenzi, T., 2021.  A Mathematical Study of the Influence of Hypoxia and Acidity on the Evolutionary Dynamics of Cancer.  Bulletin of Mathematical Biology 2021 83:7 83, 1–29.  doi:10.1007/S11538-021-00914-3.
  * Fiorillo et al. [2016] Fiorillo, M., Lamb, R., Tanowitz, H.B., Mutti, L., Krstic-Demonacos, M., Cappello, A.R., Martinez-Outschoorn, U.E., Sotgia, F., Lisanti, M.P., 2016.  Repurposing atovaquone: Targeting mitochondrial complex III and OXPHOS to eradicate cancer stem cells.  Oncotarget 7, 34084–34099.  doi:10.18632/oncotarget.9122.
  * Fowler [1989] Fowler, J.F., 1989.  The linear-quadratic formula and progress in fractionated radiotherapy.  The British Journal of Radiology 62, 679–694.  doi:10.1259/0007-1285-62-740-679.
  * Franken et al. [2013] Franken, N.A., Oei, A.L., Kok, H.P., Rodermond, H.M., Sminia, P., Crezee, J., Stalpers, L.J., Barendsen, G.W., 2013\.  Cell survival and radiosensitisation: modulation of the linear and quadratic parameters of the LQ model (Review).  International Journal of Oncology 42, 1501–1515.  URL: http://www.spandidos-publications.com/10.3892/ijo.2013.1857/abstracthttps://www.spandidos-publications.com/10.3892/ijo.2013.1857, doi:10.3892/IJO.2013.1857.
  * Frieboes et al. [2010] Frieboes, H.B., Jin, F., Chuang, Y.L., Wise, S.M., Lowengrub, J.S., Cristini V, V., 2010\.  Three-dimensional multispecies nonlinear tumor growth—II: Tumor invasion and angiogenesis.  Journal of Theoretical Biology 264, 1254–1278.  doi:10.1016/J.JTBI.2010.02.036.
  * Galeaz et al. [2021] Galeaz, C., Totis, C., Bisio, A., 2021.  Radiation Resistance: A Matter of Transcription Factors.  Frontiers in Oncology 0, 2055\.  doi:10.3389/FONC.2021.662840.
  * Garnier et al. [2019] Garnier, D., Renoult, O., Alves-Guerra, M.C., Paris, F., Pecqueur, C., 2019.  Glioblastoma stem-like cells, Metabolic strategy to kill a challenging target.  Frontiers in Oncology 9, 118\.  doi:10.3389/fonc.2019.00118.
  * Gerisch and Chaplain [2006] Gerisch, A., Chaplain, M.A., 2006\.  Robust numerical methods for taxis-diffusion-reaction systems: Applications to biomedical problems.  Mathematical and Computer Modelling 43, 49–75.  doi:10.1016/j.mcm.2004.05.016.
  * Greaves and Maley [2012] Greaves, M., Maley, C.C., 2012\.  Clonal evolution in cancer.  Nature 481, 306–313.  doi:10.1038/nature10762.
  * Greenspan [1972] Greenspan, H.P., 1972.  Models for the Growth of a Solid Tumor by Diffusion.  Studies in Applied Mathematics 51, 317–340.  doi:10.1002/sapm1972514317.
  * Grimes et al. [2014] Grimes, D.R., Fletcher, A.G., Partridge, M., 2014.  Oxygen consumption dynamics in steady-state tumour models.  Royal Society Open Science 1, 140080\.  URL: https://royalsocietypublishing.org/doi/abs/10.1098/rsos.140080, doi:10.1098/rsos.140080, arXiv:https://royalsocietypublishing.org/doi/pdf/10.1098/rsos.140080.
  * Hanahan and Weinberg [2000] Hanahan, D., Weinberg, R.A., 2000\.  The hallmarks of cancer.  doi:10.1016/S0092-8674(00)81683-9.
  * Heddleston et al. [2010] Heddleston, J.M., Li, Z., Lathia, J.D., Bao, S., Hjelmeland, A.B., Rich, J.N., 2010\.  Hypoxia inducible factors in cancer stem cells.  British Journal of Cancer 2010 102:5 102, 789–795.  doi:10.1038/sj.bjc.6605551.
  * Höc̈kel, Michael and Vaupel, Peter [2001] Höc̈kel, Michael and Vaupel, Peter, 2001\.  Tumor Hypoxia: Definitions and Current Clinical, Biologic, and Molecular Aspects.  JNCI: Journal of the National Cancer Institute 93, 266–276.  doi:10.1093/jnci/93.4.266, arXiv:https://academic.oup.com/jnci/article-pdf/93/4/266/10889684/266.pdf.
  * Hodgkinson et al. [2019] Hodgkinson, A., Le Cam, L., Trucu, D., et al., 2019\.  Spatio-genetic and phenotypic modelling elucidates resistance and re-sensitisation to treatment in heterogeneous melanoma.  Journal of Theoretical Biology 466, 84–105.  doi:10.1016/j.jtbi.2018.11.037.
  * Hoffmann et al. [2008] Hoffmann, M., Chang, H.H., Huang, S., Ingber, D.E., Loeffler, M., Galle, J., 2008\.  Noise-Driven Stem Cell and Progenitor Population Dynamics.  PLOS ONE 3, e2922.  doi:10.1371/JOURNAL.PONE.0002922.
  * Jones et al. [2001] Jones, B., Dale, R., Deehan, C., Hopkins, K., Morgan, D., 2001.  The role of biologically effective dose (bed) in clinical oncology.  Clinical Oncology 13, 71–81.  doi:https://doi.org/10.1053/clon.2001.9221.
  * Koh et al. [2011] Koh, M.Y., Lemos, R., Liu, X., Powis, G., 2011.  The Hypoxia-Associated Factor Switches Cells from HIF-1$\alpha$\- to HIF-2$\alpha$-Dependent Signaling Promoting Stem Cell Characteristics, Aggressive Tumor Growth and Invasion.  Cancer Research 71, 4015–4027.  doi:10.1158/0008-5472.CAN-10-4142.
  * Lewin et al. [2020] Lewin, T.D., Byrne, H.M., Maini, P.K., Caudell, J.J., Moros, E.G., Enderling, H., 2020\.  The importance of dead material within a tumour on the dynamics in response to radiotherapy.  Physics in Medicine & Biology 65, 015007.  doi:10.1088/1361-6560/ab4c27.
  * Lewin, Thomas D. et al. [2020] Lewin, Thomas D., Maini, Philip K., Moros, Eduardo G., Enderling, Heiko, Byrne, Helen M., 2020.  A three phase model to investigate the effects of dead material on the growth of avascular tumours.  Math. Model. Nat. Phenom. 15, 22\.  doi:10.1051/mmnp/2019039.
  * Lorenzi et al. [2016] Lorenzi, T., Chisholm, R., Clairambault, J., 2016.  Tracking the evolution of cancer cell populations through the mathematical lens of phenotype-structured equations.  Biology Direct 11, 43\.  doi:10.1186/s13062-016-0143-4.
  * Lorenzi et al. [2018] Lorenzi, T., Venkataraman, C., Lorz, A., Chaplain, M.A., 2018\.  The role of spatial variations of abiotic factors in mediating intratumour phenotypic heterogeneity.  Journal of Theoretical Biology 451, 101–110.  doi:10.1016/J.JTBI.2018.05.002.
  * Lorz et al. [2015] Lorz, A., Lorenzi, T., Clairambault, J., et al., 2015\.  Modeling the Effects of Space Structure and Combination Therapies on Phenotypic Heterogeneity and Drug Resistance in Solid Tumors.  Bulletin of Mathematical Biology 77, 1–22.  doi:10.1007/s11538-014-0046-4.
  * Marusyk et al. [2012] Marusyk, A., Almendro, V., Polyak, K., 2012.  Intra-tumour heterogeneity: A looking glass for cancer?  Nature Reviews Cancer 12, 323–334.  doi:10.1038/nrc3261.
  * Marusyk et al. [2020] Marusyk, A., Janiszewska, M., Polyak, K., 2020.  Intratumor Heterogeneity: The Rosetta Stone of Therapy Resistance.  Cancer Cell 37, 471–484.  doi:10.1016/J.CCELL.2020.03.007.
  * McKeown [2014] McKeown, S.R., 2014.  Defining normoxia, physoxia and hypoxia in tumours—implications for treatment response.  The British Journal of Radiology 87, 20130676.  doi:10.1259/bjr.20130676. pMID: 24588669.
  * McMahon [2018] McMahon, S.J., 2018.  The linear quadratic model: usage, interpretation and challenges.  Physics in Medicine & Biology 64, 01TR01.  doi:10.1088/1361-6560/AAF26A.
  * Mudassar et al. [2020] Mudassar, F., Shen, H., O’Neill, G., Hau, E., 2020\.  Targeting tumor hypoxia and mitochondrial metabolism with anti-parasitic drugs to improve radiation response in high-grade gliomas.  Journal of Experimental and Clinical Cancer Research 39, 208.  doi:10.1186/s13046-020-01724-6.
  * Najafi et al. [2019] Najafi, M., Mortezaee, K., Majidpoor, J., 2019.  Cancer stem cell (CSC) resistance drivers.  Life Sciences 234, 116781\.  doi:10.1016/J.LFS.2019.116781.
  * Noble et al. [2019] Noble, R., Burri, D., Kather, J.N., Beerenwinkel, N., 2019\.  Spatial structure governs the mode of tumour evolution.  bioRxiv , 586735doi:10.1101/586735.
  * Olivares-Urbano et al. [2020] Olivares-Urbano, M.A., Griñán-Lisón, C., Marchal, J.A., Núñez, M.I., 2020\.  CSC Radioresistance: A Therapeutic Challenge to Improve Radiotherapy Effectiveness in Cancer.  Cells 9, 1651.  URL: https://www.mdpi.com/2073-4409/9/7/1651, doi:10.3390/cells9071651.
  * Patel et al. [2014] Patel, A.P., Tirosh, I., Trombetta, J.J., Shalek, A.K., Gillespie, S.M., Wakimoto, H., Cahill, D.P., Nahed, B.V., Curry, W.T., Martuza, R.L., Louis, D.N., Rozenblatt-Rosen, O., Suvà, M.L., Regev, A., Bernstein, B.E., 2014.  Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma.  Science 344, 1396–1401.  doi:10.1126/science.1254257.
  * Petrova et al. [2018] Petrova, V., Annicchiarico-Petruzzelli, M., Melino, G., Amelio, I., 2018\.  The hypoxic tumour microenvironment.  Oncogenesis 2018 7:1 7, 1–13.  doi:10.1038/s41389-017-0011-9.
  * Picco et al. [2017] Picco, N., Gatenby, R.A., Anderson, A.R., 2017.  Stem Cell Plasticity and Niche Dynamics in Cancer Progression.  IEEE Transactions on Biomedical Engineering 64, 528–537.  doi:10.1109/TBME.2016.2607183.
  * Pisco and Huang [2015] Pisco, A.O., Huang, S., 2015\.  Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: ’What does not kill me strengthens me’.  British Journal of Cancer 112, 1725–1732.  doi:10.1038/bjc.2015.146.
  * Ruan et al. [2009] Ruan, K., Song, G., Ouyang, G., 2009.  Role of hypoxia in the hallmarks of human cancer.  Journal of Cellular Biochemistry 107, 1053–1062.  doi:10.1002/JCB.22214.
  * Saga et al. [2019] Saga, R., Matsuya, Y., Takahashi, R.et al.., 2019.  Analysis of the high-dose-range radioresistance of prostate cancer cells, including cancer stem cells, based on a stochastic model.  Journal of Radiation Research 60, 298–307.  doi:10.1093/jrr/rrz011.
  * Scott et al. [2014] Scott, J.G., Hjelmeland, A.B., Chinnaiyan, P., Anderson, A.R.A., Basanta, D., 2014.  Microenvironmental Variables Must Influence Intrinsic Phenotypic Parameters of Cancer Stem Cells to Affect Tumourigenicity.  PLOS Computational Biology 10, e1003433.  doi:10.1371/JOURNAL.PCBI.1003433.
  * Shen and Clairambault [2020] Shen, S., Clairambault, J., 2020\.  Cell plasticity in cancer cell populations.  F1000Research 9.  doi:10.12688/f1000research.24803.1.
  * Snyder et al. [2018] Snyder, V., Reed-Newman, T.C., Arnold, L., Thomas, S.M., Anant, S., 2018.  Cancer stem cell metabolism and potential therapeutic targets.  Frontiers in Oncology 8.  doi:10.3389/fonc.2018.00203.
  * Sørensen and Horsman [2020] Sørensen, B.S., Horsman, M.R., 2020\.  Tumor Hypoxia: Impact on Radiation Therapy and Molecular Pathways.  Frontiers in Oncology 10, 562\.  doi:10.3389/FONC.2020.00562.
  * Sousa et al. [2019] Sousa, B., Ribeiro, A.S., Paredes, J., 2019.  Heterogeneity and Plasticity of Breast Cancer Stem Cells.  Advances in Experimental Medicine and Biology 1139, 83–103.  doi:{10.1007/978-3-030-14366-4\\_5}.
  * Stace et al. [2020] Stace, R.E., Stiehl, T., Chaplain, M.A., Marciniak-Czochra, A., Lorenzi, T., 2020.  Discrete and continuum phenotype-structured models for the evolution of cancer cell populations under chemotherapy.  Mathematical Modelling of Natural Phenomena 15, 14.  doi:10.1051/MMNP/2019027.
  * Strobl et al. [2021] Strobl, M.A.R., Gallaher, J., West, J., Robertson-Tessi, M., Maini, P.K., Anderson, A.R.A., 2021\.  Spatial structure impacts adaptive therapy by shaping intra-tumoral competition.  bioRxiv 2020.11.03, 2020.11.03.365163.  doi:10.1101/2020.11.03.365163.
  * Suwa et al. [2021] Suwa, T., Kobayashi, M., Nam, J.M., Harada, H., 2021\.  Tumor microenvironment and radioresistance.  Experimental & Molecular Medicine 53, 1029–1035.  doi:10.1038/s12276-021-00640-9.
  * Swanson et al. [2000] Swanson, K.R., Alvord Jr, E.C., Murray, J.D., 2000.  A quantitative model for differential motility of gliomas in grey and white matter.  Cell Proliferation 33, 317–329.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-2184.2000.00177.x, doi:https://doi.org/10.1046/j.1365-2184.2000.00177.x, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1046/j.1365-2184.2000.00177.x.
  * Sweeney et al. [1998] Sweeney, K.J., Swarbrick, A., et al., 1998.  Lack of relationship between CDK activity and G1 cyclin expression in breast cancer cells.  Oncogene 16, 2865–2878.  doi:10.1038/sj.onc.1201814.
  * Syljuåsen, Randi G. [2019] Syljuåsen, Randi G., 2019.  Cell Cycle Effects in Radiation Oncology. Springer International Publishing, Cham.  pp. 1–8.  doi:{10.1007/978-3-319-52619-5\\_101-1}.
  * Tang et al. [2018] Tang, L., Wei, F., Wu, Y., He, Y., Shi, L., Xiong, F., Gong, Z., Guo, C., Li, X., Deng, H., Cao, K., Zhou, M., Xiang, B., Li, X., Li, Y., Li, G., Xiong, W., Zeng, Z., 2018.  Role of metabolism in cancer cell radioresistance and radiosensitization methods.  Journal of Experimental & Clinical Cancer Research 37.  doi:10.1186/S13046-018-0758-7.
  * Veltz [2020] Veltz, R., 2020.  BifurcationKit.jl.  URL: https://hal.archives-ouvertes.fr/hal-02902346.
  * Villa et al. [2021] Villa, C., Chaplain, M.A., Lorenzi, T., 2021.  Modeling the emergence of phenotypic heterogeneity in vascularized tumors.  SIAM Journal on Applied Mathematics 81, 434–453.  doi:10.1137/19M1293971, arXiv:https://doi.org/10.1137/19M1293971.
  * Villa et al. [2020] Villa, C., Chaplain, M.A.J., Lorenzi, T., 2020.  Evolutionary Dynamics in Vascularised Tumours under Chemotherapy: Mathematical Modelling, Asymptotic Analysis and Numerical Simulations.  Vietnam Journal of Mathematics 2020 49:1 49, 143–167.  doi:10.1007/S10013-020-00445-9.
  * Walcher et al. [2020] Walcher, L., Kistenmacher, A.K., Suo, H., Kitte, R., Dluczek, S., Strauß, A., Blaudszun, A.R., Yevsa, T., Fricke, S., Kossatz-Boehlert, U., 2020\.  Cancer stem cells—origins and biomarkers: Perspectives for targeted personalized therapies.  Frontiers in Immunology 11, 1280\.  doi:10.3389/fimmu.2020.01280.
  * Wenzl and Wilkens [2011] Wenzl, T., Wilkens, J.J., 2011\.  Theoretical analysis of the dose dependence of the oxygen enhancement ratio and its relevance for clinical applications.  Radiation Oncology 2011 6:1 6, 1–9.  doi:10.1186/1748-717X-6-171.
  * West and Slevin [2019] West, C.M., Slevin, F., 2019\.  Tumour Hypoxia.  Clinical Oncology 31, 595–599.  doi:10.1016/j.clon.2019.06.008.