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

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.


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 nongenetic 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 environmentalinduced 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 spatiallyresolved 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 Figure 1: Schematic representation of the model structure: we consider a fixed slice of tissue delimited by two vessels at location x = ±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).
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.

Model Development
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 [mm], where x ∈ [−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) [cells mm −1 ] indicates the number of tumour cells with phenotype s, at spatial location x and time t.We further introduce Φ = Φ(x, t) [cells mm −1 ], which represents the total number of tumour cells at location x and time t: 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): In Eq. (2a), the positive constant D xc [mm 2 hr −1 ] is the diffusion coefficient of oxygen and the function Γ = Γ(s, C) [mmHg(cell 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: where the positive constant D xn [mm 2 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: Here, the positive constant D sn [hr −1 ] represents the rate at which random epigenetic mutations take place, while the velocity 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 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, . . ., M .Under these assumptions, we define F as: In Eq. (2e), Φ max [cells/mm] represents the local carrying capacity of the tissue, so that if the initial condition is such that Φ(x, 0) ∈ [0, Φ max ] then Φ(x, t) ∈ [0, Φ max ] for all t > 0. The associated logisticgrowth 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, Φ; d i ) ∈ [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: 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 ∞ , on x = ±L.Finally, we close the governing equations by imposing the following initial conditions: 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 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 where and the evolution of the normalised oxygen concentration c is dictated by In Eqs. ( 4)-( 5), we have introduced the following (partially dimensionless) parameter groupings and functional forms: 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 B.2.

Specification of functional forms.
We now define the functional forms for v s , P, K, γ and 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 We define the proliferation rate P and natural death rate K in Eq. (4f) as follows: natural cell death where In Eq. (7a), a bimodal function describes the proliferation rate 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 ≈ 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 K 1 ≡ 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 ≈ 1), necrosis impact all cells regardless of their stemness level s only when oxygen levels are extremely low (c c N 1).To capture the effect of necrosis we introduce the continuous step function, H : where the parameter determines how rapidly H switches between its extremal values (H (x) ≈ 1 for x > and H (x) ≈ 0 for x < − ).In Eq. (7b), we assume that the transition between the hypoxic and necrotic niches is rapid by setting K = 0.01.In Fig. 3(a), we plot the net proliferation rate, in the absence of competition (i.e., P − K), and show how it depends on the local oxygen levels (c) and phenotype (s).Following [14], we propose the following functional form for the phenotypic advection velocity v s : In Eqs.(9b)-(9c), the normalizing factors V * ± are defined so that max , where V ± correspond to the maximum absolute values of the velocities with which cells differentiate (V + ) and de-differentiate (V − ).The positive constants ξ ± determine how rapidly the velocity profiles switch near s = 0 and s = 1.The exponents ω ± are positive integers: ω + = 1 and ω − = 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 v , which is defined as in Eq. ( 8), with 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 γ = γ(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 γ = γ(c).Furthermore, following [34,44], we suppose that γ drops rapidly when oxygen levels fall below the necrotic threshold c = c N : In Eq. ( 10), γ 0 is the maximum oxygen consumption rate for cells and γ = 0.05.We rewrite Eq. (10) in terms of the diffusion coefficient D −1 xc by introducing the non-dimensional parameter γ * = γ 0 D −1 xc .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]: In what follows, we suppose that α[Gy −1 ] and β[Gy −2 ] are dependent on the local environment, and, in particular, the oxygen concentration c, the cell density φ, 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 α (and, analogously, β) into three factors, which account for different mechanism of radio-resistance: 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 ≤ c R ), and propose the following functional form for α 1 (c) and β 1 (c): 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 α 2 (s) and β 2 (s): In Eq. (11d), the positive constants α − and β − represent the radio-sensitivity of CSCs (i.e., s ≈ 0) while ∆α and ∆β 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 α 3 (c, s, φ) and β 3 (c, s, φ) depend on the effective cell proliferation rate in the following way: In Eq. (11e), P max = max c,s P = p max 0 .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 α only since β 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., α ∼ α 2 (s)).Spatial competition is relevant at higher oxygen concentration (e.g, c > c H ).Where spatial competition is low (φ ≈ 0), differentiated cells (s ≈ 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.

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 B.1-B.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 ≈ 200 µm), we find that the time-scale of oxygen diffusion in the tumour region ( ) is of the order of seconds or minutes, which is much faster than the time-scale of cell proliferation (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 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 γ * = γ 0 L 2 D −1 xc , 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 (γ −1 0 ).Since we are interested in understanding how spatial variation in oxygen levels affects the tumour's phenotypic composition, we view γ * 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 B.1), oxygen consumption rates (γ 0 ) are likely to be tumour-and patient-specific.For this reason, when comparing different values of γ * , 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 γ 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 ∈ (0, 1) and variance σ s > 0. Thus, we prescribe where M 0 > 0 is the initial tumour burden, while s 0 and σ 2 s are chosen so that the mass of the Gaussian profile lies inside the stemness domain s ∈ (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 φ(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 Appendix 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): We note that M Φ max L represents the total number of tumour cells in the slice of tissue.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 (P) and their sensitivity to natural cell death (D) and treatment (α and β).Cells in P 1 (0 ≤ s ≤ s 1 ) are highly-resistant stem-like cells, while cells in P 2 (s 1 < s ≤ s 2 ) are highly proliferative and radio-sensitive.Finally, P 3 (s 2 < s ≤ 1) consists of terminally-differentiated cells that have exited the cell-cycle and have a high death rate 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 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: While, the fractions Π 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: 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.

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.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 ≈ 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.
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 ≈ 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.
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   4) with initial conditions defined by Eq. ( 13), and s 0 = 0.5.We observe long-time extinction when M 0 = 0.0001 (Fig. 7a) and evolution to a persistent tumour cell population when M 0 = 0.001 (Fig. 7b).Parameter values: as per Fig. 6.
development of CSCs, which drives further tumour growth.As shown in Fig. 8(a), CSCs appear in the tissue only at later times (t ≈ 15 weeks) and are localised in hypoxic regions (see π 1 (X, t) in Fig. 8(b)).At long times, the fraction of CSCs, Π 1 (t), increases and asymptotes to a stationary value of ≈ 25%.The TDCs remain the dominant cell phenotype, comprising ≈ 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.

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 ∞ (X, s), c ∞ (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 γ * 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 ∞ ≡ 0 and c ∞ ≡ 1, and corresponds to the "tumour extinction" (disease-free) scenario.In Appendix 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 ω, the particular wavenumber ω = 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 D xn , and considering an homogeneous oxygen distribution c ≡ c ∞ .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 B.1-B.2).Since c ≡ 1 when spatial heterogeneity is neglected, we have that v s (s) ≈ v + s (s).As shown in Fig. 9, when 0 < V + 1, the trivial steady state is unstable and a non-trivial stable solution exists (see the upper green branch).As V + increases , M ∞ decreases until V + crosses the critical point V (cr) + ≈ 6.06 × 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, We distinguish between stable (solid) and unstable (dashed) solution branches and highlight bifurcation points with circles.
this will not affect the stability of the trivial steady state.In other words, the location of the critical point is independent of D xn or the parameter associated to the oxygen distribution, i.e., γ * and D xc .

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 Appendix 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 γ * , which affects the equilibrium oxygen distribution.
The results presented in Fig. 10 show how, as γ * increases through the critical value γ (cr) ≈ 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, SN 1 and 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 (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 (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 ∞ = 0) is the only stable solution.For larger values of γ * (see Figs. 10(b) and 10(c)), the tumour burden on the upper-most branch decreases and SN 1 moves to the right, passing beyond V (cr) + .Note that in Fig. 10(c), SN 1 does not appear as it is located beyond the interval under consideration, V + ∈ [10 −4 , 10 −3 ].The second saddle-node bifurcation point, SN 2 , also moves to the right (i.e. higher V + ) but it is less sensitive to changes in γ * .In particular, comparing Fig. 10 4) with γ * being our bifurcation parameter.We compute the diagrams for three fixed values of V + : (a) V + = 2.5 × 10 −4 , (b) V + = 5 × 10 −4 and (c) V + = 8 × 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 γ * , we fix V + and compute the system steady states, viewing γ * as our continuation parameter.The corresponding diagrams, presented in Fig. 11, show that for slow advection velocities (i.e., V + = 2.5 × 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 × 10 −4 (see Fig. 11(b)), SN 2 shifts towards the bottom-right of the diagram and the size of the bistable region increases.The position of 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 ∞ initially increasing past the supercritical saddle-node point (SN 1 ) and subsequently decreasing as γ * increases further.This contrasts with the lower stable branch, which decreases monotonically until it disappears at the saddle-node bifurcation point 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 (SN 2 ).As V + → V (cr) + from below, SN 2 shifts to the right (SN 2 (γ * ) → ∞).Since SN 2 becomes infinite at a finite value of V + , our analysis suggests that the transcritical bifurcation and SN 2 annihilate each other at infinity, i.e., at γ * → +∞.This is validated by using fold continuation methods which follow the locations of SN 1 and SN 2 in the (V + ,γ * ) plane (the results are presented in Appendix C.3).Consequently, in Fig. 11(c), where V + > V (cr) + , the bifurcation diagram simplifies.For 0 < γ * < 1.93, extinction is the only possible long term solution; however as γ * increases beyond SN 1 (i.e., γ * > 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.
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 (γ * = 2.0) vs high (γ * = 4.0) oxygen consumption rates, while fixing V + = 8.0 × 10 −4 .When γ * = 2.0 (i.e., just to the right of the critical saddle node bifurcation SN 1 , see red dot in Fig. 12(a)), equilibria A and B have a similar tumour burden, M ∞ , 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 ∞ (X) < c H ). While steady state B also possesses an hypoxic region (here mild hypoxia as min X c ∞ (X) ≈ 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 φ ∞ (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 ≈ 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 ∞ , decreases monotonically for sufficiently large values of γ * .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 φ ∞ (X).For steady state C the maximum value of φ ∞ occurs at the oxygenated boundary (X ≈ 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 × 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.

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 × 10 −4 and consider two values of the oxygen consumption: γ * = 2.0 and γ * = 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 σ 2 s of the Gaussian so that σ 2 s = 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 ∈ [0.25, 0.75] so that boundary effects do not distort the Gaussian profile.
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. 13a (i.e., γ * = 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 ≈ 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 γ * .In Fig. 13b, tumour survival occurs for a wider range of values of s 0 and M 0 than when γ * = 2.0 (Fig. 13a).
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 1) initially present increases, the smaller the value of M 0 for which invasion occurs.On the other hand, as s 0 → 1, even large tumours may be eliminated.In this case,  )); 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).
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 γ * , 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 γ * , the threshold tumour burden required for the formation of hypoxia is larger, so that even large tumours may become extinct.

Dynamics in the presence of treatment
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×10 −4 ), so we know that tumour extinction is possible, and compare the response to treatment of tumours with different oxygen consumption rates: γ * = 2.0 and γ * = 4.0.For each value of γ * , 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).

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 (γ * = 2.0), while Fig. 15(b) shows the response of a tumour with a high oxygen consumption rate (γ * = 4.0).When γ * = 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 γ * = 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 γ * .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 B.1-B.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 c(t) = min X (X, t).The results presented in Fig. 16 suggest that during UHFRT the tumour composition is independent of γ * .In both cases, RT re-oxygenates the tumour (note that the minimum oxygen level recorded in the tissue, c, remains above c H throughout treatment).Therefore, although Π 1 (t) increases after each dose of radiotherapy (since RT selects for radioresistant CSCs), Π 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), Π 1 (t = 4) < Π 1 (t = 0) ≈ 0.25.By contrast, the fraction of differentiated (highly-proliferative) cells Π 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 γ * = 2.0, c remains above the threshold c H , an hypoxic region forms at t ≈ 6 weeks for γ * = 4.0.Only at this time point, does the internal composition of the two tumours start to diverge.While for γ * = 2.0, at long times all DCs become TDCs and the tumour is eliminated, when γ * = 4.0, DCs located in the hypoxic region de-differentiate, producing CSCs which enable the tumour to survive and continue to grow.Focusing now on FRT (see Fig. 17), we observe that the evolution of c differs during the early stages of treatment for the two values of γ * .For γ * = 2.0, one cycle of FRT is sufficient to raise oxygen levels above the hypoxic threshold c H , while for γ * = 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 Π 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 γ * = 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 c changes markedly with γ * , 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 c (see Fig. 16) appears to be independent of γ * during the treatment period so that it can not be used as an early indicator of outcome in response to UHFRT.
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 γ * , 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 γ * 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].

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 α and β 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 α and β 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 (γ * = 4).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 (≈ 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.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 Π 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.
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 (≈ 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.
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 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.
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.

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 γ * 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 γ * can successfully invade the tissue, while larger tumours with lower oxygen consumption rates γ * may eventually be eliminated.This suggests that the size of the basin of attraction of the "tumour-extinction" steady states decreases with γ * .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 γ * 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 (≥ 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.

Function f (n, t):
Solve for c using fix point iterations Assemble the vector f using the computed oxygen distribution c return

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 B.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 B.2. where δφ(X) = The question now is whether λ = κ is an eigenvalue for the full system Eqs.(C.2).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 ∞ .
Let us consider the case ω m ≡ 0 for which wlog Λ 0 ≡ 1.In this case, the eigenfunction associated to the eigenvalue κ is δn = W κ(s).The eigenfunction is defined up to a multiplicative constant which we set so that ) is a subset of λ m,k ≤ κ, κ will be the eigenvalue with largest real part and therefore dictates the stability of the trivial solution.From [14], we know that κ decreases from positive to negative values as the magnitude of the advection velocity V + increases.The critical value V (cr) + at which κ = 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 γ * , D xn and 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.

Appendix C.2. Numerical continuation and stability estimates
In order to compute the equilibrium states (n ∞ (x, s), c ∞ (x)) of the full system (4)-( 5) in the absence of treatment, we simply set to zero all time derivatives to obtain: 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. C.22, we see that for sufficiently small γ * , 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 γ * → ∞.Here V (cr) + corresponds to the value of V + for which the system has a transcritical bifurcation (see Fig. 10).As discussed in Appendix C.1, the location of V (cr) + is independent of the parameter γ * and divides the (γ * , 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.

Figure 3 :
Figure 3: Surface plot of the functional (a) P(s, c)−K(s, c) and (b) vs(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 B.2.

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

Figure 6 :
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 B.1-B.2, with γ * = 4 and V + = 8 × 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.7are taken.

Figure 7 :
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.7a) and evolution to a persistent tumour cell population when M 0 = 0.001 (Fig.7b).Parameter values: as per Fig.6.

Figure 12 :
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 (Π i see Eq. (15)) at points A, B, and C; (c) we illustrate the equilibrium solution, i.e., oxygen profile c∞, re-scaled cell number density φ∞, and the local distribution of cell phenotype n∞(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 × 10 −4 and γ * = 2.0 (A,B) or γ * = 4.0 (C).

Figure 13 :
Figure13: 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 γ * .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).

Figure 14 :
Figure 14: Schematic diagram showing how treatment is incorporated into the model.Treatment is repeated over a period of nrep weeks so that the final dosage delivered is 10 × nrep 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.

Figure 15 :
Figure 15: Evolution of the tumour burden, M (t), during treatment for different values of the parameter γ * : (a) γ * = 2.0 and (b) γ * = 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 B.1-B.2.

Figure 16 :
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) γ * = 2.0 and (b) γ * = 4.0.The top panels show the evolution of the minimum oxygen level within the tumour, c(t) = min 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.

Figure 17 :
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) γ * = 2.0 and (b) γ * = 4.0..In the top panels we plot the minimum oxygen level c recorded in the tissue at any time and compare it with the relevant oxygen threshold c H and c R .

Figure 18 :
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.

Figure 19 :
Figure19: 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.

Figure 20 :
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).

1 0
δn(X, s)ds.We note that Eq. (C.2b) for δc decouples from Eq. (C.2a) for δn.We can therefore solve Eq. (C.2a) independently of δc, and seek a separable solution of the form:δn = Λ(x)W (s) (C.3)where Λ and W needs to satisfy: s W + (p(s, c ∞ ) − f (s)Eq.(C.4b) we have used the fact that the velocity v s vanishes at the boundaries.Using standard arguments, we find that the eigenfunctions Λ m associate to the eigenvalue ω m are given byΛ m = cos (ω m X) , ω m = πm, m = 0, 1, 2, . . . .(C.5a)so that the eigenfunctions W associated to the phenotypic axis, s, must satisfy the following ODE:ω 2 m D xn + λ W = d ds D sn dW ds − v s W + [p(s, c ∞ ) − f (s)] W ≡ M(W ), (C.5b) operator 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 M has been investigated in more detail in[14].While we have no explicit solution for the eigenvalues κ k and eigenfunctions W k of M, we can show that κ k are real and bounded from above.If D xn is constant and W k denotes the eigenfunctions of M associated to the eigenvalue κ k , then δn m,k = Λ m (X)W k (s) is an eigenfunction of the operator, Eq. (C.4a), with real eigenvalue λ m,k = κ k − ω 2 m D xn .Since ω 2 D xn ≥ 0, we find that λ m,k ≤ κ k and λ m,k < λ m ,k if m < m.Consequently, λ m,k < 0 for all k and m if and only if κ = max k κ k < 0.

Table B .
[14]ummary of typical dimensional values for some of the parameters in Eq. (2), togehter with supporting references.Table B.2: Summary of non-dimensional parameters groups that appears in Eq. (4).Parameter values are consistent with those reported in TableB.1 and/or in[14].