Phenotypic variation modulates the growth dynamics and response to radiotherapy of solid tumours under normoxia and hypoxia

In cancer, treatment failure and disease recurrence have been associated with small subpopulations of cancer cells with a stem-like phenotype. In this paper, we develop and investigate a phenotype-structured model of solid tumour growth in which cells are structured by a stemness level, which varies continuously between stem-like and terminally differentiated behaviours. Cell evolution is driven by proliferation and apoptosis, as well as advection and diffusion with respect to the stemness structure variable. We use the model to investigate how the environment, in particular oxygen levels, affects the tumour's population dynamics and composition, and its response to radiotherapy. We use a combination of numerical and analytical techniques to quantify how under physiological oxygen levels the cells evolve to a differentiated phenotype and under low oxygen level (i.e., hypoxia) they de-differentiate. Under normoxia, the proportion of cancer stem cells is typically negligible and the tumour may ultimately become extinct whereas under hypoxia cancer stem cells comprise a dominant proportion of the tumour volume, enhancing radio-resistance and favouring the tumour's long-term survival. We then investigate how such phenotypic heterogeneity impacts the tumour's response to treatment with radiotherapy under normoxia and hypoxia. Of particular interest is establishing how the presence of radio-resistant cancer stem cells can facilitate a tumour's regrowth following radiotherapy. We also use the model to show how radiation-induced changes in tumour oxygen levels can give rise to complex re-growth dynamics. For example, transient periods of hypoxia induced by damage to tumour blood vessels may rescue the cancer cell population from extinction and drive secondary regrowth. Further model extensions to account for spatial variation are also discussed briefly.


Introduction
Understanding of the mechanisms by which cancer is initiated and progresses continues to increase, and, yet, cancer remains one of the leading causes of premature mortality worldwide and a major barrier to increasing average life-expectancy. For example, in 2018, 9.6 million people are estimated to have died of cancer [1]. Furthermore, treatment outcomes can differ markedly between patients with the same cancer type, with the emergence of resistance being one of the major causes of treatment failure.
Over the past twenty years, there has been a major shift in our perception of solid tumours; they are now regarded as heterogeneous tissues in which malignant cells interact with normal cells and shape their environment in ways that favour malignant growth [2]. Cancer stem cells (CSCs) were introduced to explain intra-tumour heterogeneity via the CSC hypothesis [3]. This hypothesis proposes that, while CSCs may comprise only a small fraction of the total cell population, their high clonogenic potential and their ability to produce more mature, or specialised, cancer cells enables them to create an entire tumour [4]. As CSCs are found to be resistant to standard treatments, they are recognised as a major cause of disease recurrence and treatment failure [5,4,6]. These observations have stimulated the development of novel therapeutic strategies which aim to eradicate CSCs [7,8,9,10]. In practice, the plasticity of CSCs represents a major obstacle to such treatments. Additionally, CSCs can adapt to their local micro-environment, and remodel it to create and maintain a niche which supports their survival [11].
Increasingly, researchers are turning to mathematical models in order to understand how CSCs affect the growth and composition of tumours, particularly their heterogeneity and response to treatment. These models often decompose the tumour into a series of compartments, each representing a particular cell subtype. For example, in [7], Enderling distinguishes cancer stem cells (CSCs) and cancer cells, whereas Saga and coworkers distinguish radio-resistant and radio-sensitive cells [12], and Scott and colleagues distinguish tumour-initiating cells (or CSCs), transit-amplifying cells and terminally differentiated cells (TDCs) [13]. Thus, most compartmental models are based on the CSC hypothesis which assumes that it is possible to distinguish between cancer stem cells and the tumour bulk. However, this paradigm has been challenged by recent experimental studies [14,15] that highlight the phenotypic heterogeneity and plasticity of cancer cells, whose clonogenic (or stemness) potential can be altered by the surrounding micro-environment (extrinsic forces). These findings have led to a new hypothesis for intratumoural heterogeneity, based on adaptive CSC plasticity [16]. Under this hypothesis, cancer cells move between stem-like and terminally differentiated states in response to extrinsic (environmental) and/or intrinsic (random epigenetic mutation) forces. Remarkably, the development of state-of-the-art experimental tools, such as single-cell RNA-seq, means that it is now possible to track the evolution of stemness traits [17,18], rendering this an ideal time to develop mathematical models that can explore these concepts.
Compartmental models can be used to study adaptive CSC plasticity , by allowing transitions between different compartments. However, since they assume that the tumour comprises distinct cell populations, with distinct properties, they are unable to account for continuous variation in cell properties. An increasingly popular mathematical approach for describing population heterogeneity and plasticity characterises tumour cells by their position on a continuous phenotypic axis. Position on the phenotypic axis determines cell properties such as resistance to treatment [19,20,21,22,23] and/or metabolic state [24,25]. This approach is motivated by concepts from evolutionary ecology, such as risk-spreading through spontaneous (epigenetic or genetic) variations and evolutionary pressure [26]. The resulting models are typically formulated as systems of reaction-diffusion equations [24,22,25], with an advective transport term sometimes included to account for biased mutation dynamics [21] or adaptive phenotypic switches [19,27,28].
In this paper, we formulate a mathematical model that accounts for the evolution of a cancer cell population along such a stemness axis in response to extrinsic and intrinsic stimuli. Initially, we focus on the plastic response of cells to changes in nutrient levels, in particular oxygen. This is motivated by recent experimental studies [29,30,31,32] suggesting that hypoxia (i.e. low oxygen levels) is a key driver of cell de-differentiation. From this point of view, spatial heterogeneity may introduce significant ad-ditional complications: as oxygen diffuses into a tumour and is consumed by cells, spatial gradients in the oxygen levels are established. In this way, local micro-environments characterised by normoxia, hypoxia and necrosis form as the distance to the nearest nutrient supply (i.e., blood vessels) increases [21,23,25]. For simplicity, we postpone consideration of such spatial complexity to future work and focus, instead, on a well-mixed setting where oxygen levels are homogeneous and prescribed. This idealised scenario allows us to investigate how cell properties, such as proliferation, apoptosis and adaptive response to environmental signals, contribute to the emergence of heterogeneous stemness levels in the population and the long term tumour composition. In this regard, we are interested in identifying conditions under which CSCs are favoured. We then extend the model to account for treatment via a phenotypically-modulated linear-quadratic model of radiotherapy (see, e.g., [33,34,12] for recent discussions) which accounts for differential radio-sensitivity of CSCs [10]. This allows us to investigate how different radiotherapy protocols perturb the phenotypic distribution and subsequent regrowth of the tumour.
In practice, stemness is just one of multiple traits that regulate cell behaviour and heterogeneity. We, therefore, anticipate that future models will combine multiple phenotypic axes or synthetic dimensions, such as stemness and metabolic state [24,21]. Given the complexity of such multi-dimensional models, it is important first to understand these aspects separately. Noting that considerable mathematical effort has been devoted to investigating cancer metabolism [35], we choose here to focus on population heterogeneity with respect to a continuously varying stemness axis. We hope that in the long term this work will help motivate a systematic experimental characterization of cell plasticity and phenotype.
The remainder of the article is organised as follows. In Section 2, we present a well-mixed, spatially homogeneous, model of solid tumour growth in response to a prescribed oxygen concentration. We first investigate the population dynamics in the absence of treatment, considering both normoxic and hypoxic conditions. Numerical results are presented in Section 3. As a partial validation of the numerical results, we use spectral stability analysis to characterise the long time behaviour of the solutions. Section 4 focuses on tumour cell responses to different radiotherapy protocols. As in Section 3, we simulate responses under normoxia and hypoxia, but we also consider situations in which the environment alternates between periods of hypoxia and normoxia in order to explore the different ways that radiotherapy can alter tissue oxygenation. Finally in Section 5, we summarize our key findings and propose possible directions for future work. We also present preliminary results showing how accounting for spatial and phenotypic variation may affect a tumour's growth and response to radiotherapy.

Model Formulation
We consider the temporal evolution of a heterogeneous population of tumour cells, N (z, t), where t ≥ 0 denotes time and z (0 ≤ z ≤ 1) represents their stemness or clonogenic capacity. As shown in Figure 1, z = 0 corresponds to cancer stem cells (CSCs) which have the maximum level of stemness, and z = 1 corresponds to terminally differentiated cells (TDCs), which have lost their proliferative capacity and which can either enter replicative senescence or undergo cell death [36]. We assume that the population dynamics may by described by a reaction-advection-diffusion equation (see Eq. 1a below) which accounts for two essential physical/ecological processes. First, cells move along the stemness axis (i.e., in the z-direction) in response to extrinsic (micro-environment) and intrinsic (random epimutation) forces [13], which give rise to advective and diffusive fluxes respectively. Finally, the effect of natural selection on the population is represented by the fitness function F , which models the net growth rate of the cells. Figure 1: Schematic representation of the well-mixed, phenotypic model. We associate with each cell a stemness level z, which varies continuously between the cancer stem cell state (CSCs, with z ∼ 0), the differentiated cell state (with z ∼ 0.5) and the terminally differentiated cell state (TDCs, with z ∼ 1).
While multiple nutrients and growth factors regulate the growth rate (or fitness function F ) and phenotypic adaptation (i.e., the advective velocity v z ) of the tumour cells, here, for simplicity, we focus on a single nutrient, specifically oxygen. The critical role of low oxygen levels, or hypoxia, in cancer has long been recognised due to its association with cell quiescence and poor therapeutic outcomes [21,34,12]. Recent experimental results [15] have shown that hypoxia also plays a role in de-differentiation by regulating pathways associated with a stem-like phenotype. We account for these phenomena in our model by assuming that all cells are exposed to the same level of oxygen, c = c(t), which mediates the values of the fitness function, F , and the advection velocity, v z ; the latter feature distinguishes our work from existing theoretical models in which intrinsic forces are assumed to dominate phenotypic variation (i.e., v z = 0) [24,25]. By combining the processes mentioned above, we deduce that the evolution over time t and along the phenotypic axis z of the cell concentration, N (z, t), is governed by the following non-local partial differential equation (PDE) and associated boundary and initial conditions: In Equation (1), the non-negative constant θ represents the rate at which cells diffuse along the phenotypic axis, due to random epigenetic mutations, Φ(t) denotes the density of cells in the domain at time t, and N 0 (z) is the initial distribution of cells along the phenotypic axis. In ecology, the function F is referred to as fitness landscape which is a mathematical representation of natural, or Darwinian, selection [37]. We suppose it has the following form: In Equation (1e), p = p(z, c) denotes the phenotype-dependent growth rate of the cells (see Section 2.1 for details). It is multiplied by a non-local (in the phenotypic sense) logistic term, with constant carrying capacity Φ max , to capture intra-population competition for space and other resources. We assume that oxygen levels remain sufficiently high so that necrosis can be neglected. Hence, the death rate, f , accounts only for natural cell death, or apoptosis, which is assumed to occur at a rate which is independent of the oxygen concentration, c(t). Radiotherapy (RT) also contributes to cell death and, in so doing, reduces cell fitness. We suppose that M rounds of RT are administered at discrete times t i (i = 1, 2, . . . , M ). After each treatment dose, the proportion of cells of phenotype z that survive is denoted by the survival fraction SF (z, c). By allowing SF to depend on z, we can account for phenotypic-dependent radio-sensitivity, and, for example, view the CSCs (i.e. z = 0) as the most radio-resistant tumour subpopulation [4]. Additionally, the dependence of SF (z, c) on c(t) enables us to account for differential radio-sensitivity under normoxia and hypoxia [38,39]. In contrast to [33], where the term (1 − SF ) is used to capture cell death due to radiotherapy, here we use the term log(1/SF ), to ensure that the jump in tumour cells following each dose of radiotherapy is consistent with the Linear-Quadratic (LQ) model. We now partially rescale our model by recasting the dependent variables N and Φ in the following way: where the units of time, t [hr] are preserved in a dimensional form to facilitate the interpretation of the results. Under this rescaling, equations (1) become In order to complete the model, it remains to specify several functional forms; this will be done in Sections 2.1 and 2.2. Extending the model to account for spatial variation is presented in Appendix A, and preliminary results are included in Section 5 (a full investigation of the spatially-extended model is postponed to future work).
In what follows, we assume that oxygen concentration c has been rescaled so that c = 1 corresponds to physiological oxygen levels, namely physoxia, which is about 8% oxygen [40]. When considering hypoxia, we focus on mild hypoxia, fixing c = 0.2 which corresponds to 1.6% oxygen in standard units (see Appendix A.1 for details). At this oxygen concentration, necrosis can be neglected; it typically occurs at lower oxygen tensions (approximately 0.1% oxygen in standard units).
Unless otherwise stated, we assume that the tumour initially comprises a small population of differentiated cells so that where the positive constants φ 0 and σ specify the initial size and phenotypic variance of the population. The proportion of CSCs is often used to characterise heterogeneous populations of cancer cells. CSCs are typically identified by their expression of specific markers (such as CD44/CD24 and ALDH1, depending on the tumour type [10]); thresholds in these markers are used to distinguish stem from differentiated cancer cells. Since our model treats stemness as a continuously varying cell property, we introduce a threshold z * ∈ (0, 1) in our simulations, and classify cells with 0 < z < z * as CSCs. We therefore define the proportion of stem cells at time t to be: As a further statistical feature of the cell population, we introduce the phenotypic mean, µ(t), which is defined as follows: In the absence of suitable experimental data, it is difficult to specify many of the parameters and functional forms in Equations (3). For this reason, we focus on identifying the qualitative behaviours that the model exhibits across a range of 'biologically-reasonable' situations.

Fitness Landscape
When considering the fitness landscape, we assume that, for fixed values of c, the proliferation rate, p(z, c) has a multi-peaked profile, with local maxima centred around z = 0 and z = 0.55, representing respectively cells with stem-like (z = 0) and intermediate phenotypes (z = 0.55, this value being arbitrary). As shown in Figure 2, this choice reduces the overlap of the two Gaussian profiles while maintaining the proliferation rate at z = 1 close to zero. This asymmetry also emphasises that, under normoxia, more stem-like cells (i.e. z < 0.5) proliferate at lower rates than more differentiated cells (i.e. z > 0.5). Different environmental conditions (i.e., oxygen concentrations), will create distinct ecological niches each of which will favour a particular phenotype. We account for this effect by assuming that the amplitude of the peaks in the proliferation rate are oxygen dependent. Accordingly we write: where p 0 (c) and p 1 (c) are Hill-Langmuir type equations with fourth order exponents, so that the growth rate decays rapidly when c ∼ K i . We assume that differentiated cells are fitter than CSCs under normoxia and, therefore, choose p max 1 > p max 0 . At the same time, we note that chronic hypoxia is widely considered to favour CSCs [41,42,43]. The plasticity of CSCs enables them more easily to adapt their metabolism to changing nutrient levels than differentiated cells [29,44] and, therefore, to survive and proliferate in challenging conditions. This behaviour contrasts with that of differentiated cancer cells which tend to become quiescent when exposed to hypoxia. We account for these effects by assuming K 0 K 1 . When we consider the rate of cell death due to apoptosis, f (z), we note that apoptosis occurs predominantly when cells lose their clonogenic capacity. As such, it predominantly affects only TDCs with z ∼ 1. Motivated by the mathematical models developed in [7,13], we propose the following monotonically increasing function for f (z): Even though they may not proliferate, TDCs compete for space and resources and, thus, impact the tumour dynamics. In what follows, we consider two different cases. First, guided by experimental results reported by Driessens et al. [45], we assume that apoptosis of TDCs occurs on a much longer timescale than that on which cells proliferate so that d f << max z p(z; 1). In the second case, the rates of cell proliferation and apoptosis are assumed to be comparable. This situation represents a tumour with high cell turnover and, as we will see, gives rise to a tumour population with higher clonogenic capacity.
In Figure 2, we sketch the fitness landscape F (z, 0, t; c) for different environmental conditions in the absence of treatment and competition. In doing so, we have neglected competition and radiotherapy in Equations (3e), where p and f are defined by Equations (6)- (7).  Table 1.
Regions of positive and negative fitness are highlighted in green and red, respectively. We now consider the impact of radiotherapy on cell fitness. As mentioned above, CSCs possess protective mechanisms that enable them to withstand damage caused by radiation and oxidative stresses [46,47,48,4,6,10,49]. They are, therefore, more resistant to treatment than their differentiated counterparts. It is well known that local oxygen concentration levels also affect treatment outcomes [50,51]. While we account for this effect in the full spatial model (see Appendix A), here we focus on the role of phenotypedependent radio-sensitivity. In particular, we adapt the standard Linear-Quadratic (LQ) model so that the tissue specific coefficients, α(Gy −1 ) and β(Gy −2 ), are phenotype dependent: where d is the radiation dose in Grays (Gy). Equation (8a) is the natural, continuum extension of previous works [52,12], in which two-compartment models are used to describe the time-evolution of cancer cells and cancer stem cells exposed to radiotherapy, and CSCs are assumed to be radio-resistant. Accordingly, here, we assume α and β are increasing functions of the phenotype z [12, 6, 10] of the following form: In Equations (8b)-(8c), ξ R , α min,max and β min,max are non-negative constants with α min < α max and β min < β max . Where possible, parameter estimates are taken from the literature (see [12] for estimates of α min,max and β min,max ); the value of ξ R = 0.2 is instead chosen so that differentiated cells (i.e. z > 0.5) have maximum sensitivity to treatment (i.e., α(z) ∼ α max for z > 0.5).
We consider three different parameter sets (see Table 2); they may represent three cell populations which differ in their sensitivity to radiotherapy (RT). For cases R1 and R3, CSCs (with z ∼ 0) respond in the same way to RT, whereas differentiated cancer cells (with z > 0.5) respond differently. For case R1, the small value of α max /β max for the sensitive cells (z = 1) corresponds to a late responding tissue, whereas for case R3, the large value of α max /β max corresponds to an early responding tissue, with a low repair capacity, for which fractionation is known to be beneficial [53]. Finally, case R2 is intermediate between cases R1 and R3. By assuming heterogeneity in  the cell response to RT, we allow consideration of the selective pressure of RT. For a given dosage and LQ model, differences in the radio-sensitivity of CSCs and differentiated cells are determined by the ratios α min /α max ∈ (0, 1) and β min /β max ∈ (0, 1). When both fractions are small, CSCs are more likely to survive RT than their differentiated counterparts and, therefore, the selective pressure of RT on the population is high. By contrast, as α min /α max and β min /β max approach the value of unity, RT offers no selective advantage to CSCs as, at leading order, the response is independent of phenotype. The latter also depends on the specific dose applied. For example, for high doses the quadratic term in Equation (8a) is dominant and the selective pressure is only associated with the value of β min /β max . By contrast, for lower doses, the linear and non-linear effects contribute to cell killing and, so, the selective pressure of RT is associated with both α min /α max and β min /β max . For these reasons, we will consider two different RT protocols: either a single dose of 10 Gy is delivered or a fractionated schedule is used (here five doses of 2 Gy are delivered over five consecutive days [54,55]). While R2 is expected to have the least RT selective pressure in both scenarios, this might be higher in R1 or R3 depending on the treatment protocol considered.

Structural Flux
Plasticity is an essential feature of phenotypic adaptation to changing environmental conditions [14,37]. It assumes that cells with the same genome can acquire distinct phenotypes depending on their epigenetic status, which is also inheritable. Phenotypic variation may be mediated by random (spontaneous) epigenetic mutations [22], which we assume to be rare. We account for this effect by including in the structural flux a diffusion term with a constant diffusion coefficient θ = 5 × 10 −6 hr −1 (see Equation 1a). Such random muta-tions should not favour any specific phenotype, and Darwinian selection (i.e. the fitness function F ) drives phenotypic evolution of the population. This aspect has been widely studied in previous work in order to investigate how cells adapt to different environments [24,22,25]. At the same time, there is evidence that phenotypic switching may be mediated by environmental factors via Lamarckian selection (or induction) [37]. In this framework, cells adapt to their environment [14,56] by following a preferential (biased ) trajectory in phenotypic space. We can, therefore, envisage situations in which a subpopulation may be prevalent in a population without being the fittest subpopulation (i.e. the population with the highest proliferation rate). For example, recent studies have identified cell de-differentiation and CSC maintenance as stress responses to harsh environmental conditions [37], including hypoxia. More specifically, cells respond to hypoxic stress by up-regulating Hypoxia Inducible Factors (HIFs) which, in turn, promote the expression of stem-related genes [29,30,31,32]. HIF suppression has also been linked to cell differentiation and reduced levels of stemness [57]. We account for such micro-environment mediated adaptation by incorporating an advective term in the structural flux. Cells are assumed to evolve along the stemness axis with a velocity v z = v z (z, c), that depends on the oxygen concentration c and cell phenotype z. Under normoxia, cells tend to differentiate, and v z > 0. From this point of view, the model is similar to classical age-structured models [58,59], with v z being analogous to a maturation velocity. In our model, however, ageing (i.e. differentiation or loss of clonogenic potential [13]) may be reversible. For example, under hypoxia (i.e. c ≤ c H ), we assume v z < 0 (see Figure 3) and a more stem-like character is promoted. Figure 3: Series of sketches showing how v + z and v − z , as defined by Equations (9b) and (9c) respectively, change as the parameters ξ ± and ω ± vary.
Combining the above observations, and motivated in part by recent, similar considerations [21], we propose the following functional forms for the phenotypic drift term, v z : where H is a smooth variant of the Heaviside function approaching the latter in the limit of → 0 (i.e., H (x) = (1 + tanh( −1 x))/2). In Equations (9), the normalising factors V * ± ensure that (max z v ± z ) /V ± = 1 and V ± (hr −1 ) corresponds to the magnitude of the velocity. Further, by controlling the advection speed along the stemness axis, V −1 ± determines the timescales for maturation and de-differentiation. The parameters ξ ± regulate the slopes of v z at the boundaries z = 0, 1. As shown in Figure 3a, when ξ ± 1, the advection velocity is steep when z ∼ 0, 1 and flatter elsewhere. This functional form is similar to that proposed in [21]. For larger values of ξ ± , the variation is more gradual, with a single maximum (or minimum) near z ∼ 0.5 (see Figure 3b). The exponents ω ± allow us to tune the symmetry/asymmetry in v z and also to modulate the flux at the boundaries (see Figure 3). For example, if ω + = 2, then v(0) = ∂ z v(0) = 0 which means that CSCs will be less likely to differentiate compared to the case ω + = 1. In the absence of experimental data with which to specify the parameters in the phenotypic drift velocity, we consider combinations of the following parameter sets: In summary, our phenotype-structured model for the growth and response to radiotherapy of a solid tumour is defined by Equations (3)- (9). A list of the model parameters and estimates of their values can be found in Table A.4 in Appendix A.

Population Dynamics in the Absence of Treatment
In this section, we present numerical solutions of Equations (3)- (7) and (9) showing how, in the absence of treatment, the tumour cell distribution along the stemness axis evolves under normoxia and hypoxia. Our numerical solutions are generated using the method of lines, with discretisation performed in the z-direction. In more detail and following [60], we use a finite volume scheme, opting for a Koren limiter to control the advection component of the structural flux. In this way, we reduce (3) to a system of time-dependent, ordinary differential equations which can be solved in MATLAB using ode15s, an adaptive solver for stiff equations. The numerical simulations are validated in Section 3.3 where we perform a linear stability analysis. The associated eigenvalue problem is solved numerically using MATLAB's chebfun package [61].

Normoxic Conditions
In well-oxygenated environments, the advection velocity is positive and cells are driven towards a terminally differentiated phenotype, with z = 1. Depending on the balance between the advective flux and cell renewal (i.e., Darwinian selection and Lamarckian induction), the model predicts a variety of long-time behaviours: the system relaxes to its steady state via damped fluctuations or monotonically. We start by considering symmetric velocity profiles (see Figure 3a). As summarised in Figures 4 and 5, as the magnitude of the advection velocity, V + , and its steepness, ξ + , are varied, the system exhibits different long time behaviours, even though the dynamics at early times are similar for all parameter sets considered (see Figure 4). If simulations are initialised with a small population of cells with z ∼ 0.5, then the dynamics are initially dominated by proliferation. Over time, as φ increases, competition slows the cell proliferation rate and phenotypic advection becomes more important. As the cells mature, they accumulate near z = 1, and the rate of natural cell death exceeds the rate of cell proliferation. From this time onwards, the growth curves corresponding to different parameter sets start to deviate.
For example, in case A.2, the system rapidly relaxes to a non-zero steady state distribution characterised by cells with medium clonogenic capacity (i.e., a mix of highly proliferating and terminally differentiated cells or TDCs). Similarly, for cases C.1 and C.2, the cell density, φ(t), decays exponentially  7) and (9), showing how the cell distribution, n(z, t), the phenotypic mean, µ(t), and the cell density, φ(t), change over time when we use a symmetric velocity profile (i.e., ω + = 1 in Equation (9)). As V + increases and ξ + decreases, the system can be driven to extinction. See Figure 5 for the values of the other model parameters.
to extinction at a rate dictated by d f . In other parameter regimes, the relaxation phase is characterised by damped fluctuations. In case A.1, for example, fluctuations are driven by the interplay between apoptosis, competition and advection. As TDCs are eliminated, the reduction in competition allows re-growth of highly proliferative cancer cells (i.e., z ∼ 0.55). As these cells proliferate, competition slows growth and advection becomes dominant, resulting in the alternating pattern of red and white stripes observed in the surface plot for n(z, t) shown in Figure 4 for case A.1. Over time, the fluctuations decay and the system relaxes to its steady state distribution. In Section 3.3, we present a complementary investigation of this behavior, relating the damped oscillations to a complex eigenvalue in the linearisation about the equilibrium solution.
Focusing on the long time behaviour, the symmetric advective profile gives rise to a population with a unimodal equilibrium distribution where the location of the peak is dictated by the values of the other parameters. For example, for small values of the maximum death rate, d f (see case A.1), the distribution is skewed towards z = 1, while for higher values of d f the peak is shifted towards the centre of the domain. These observations are summarized in Figure 5, where we have further analysed how the properties of the equilibrium distribution depend on other parameters in the model. We note that as the advective velocity increases (i.e., larger V + ) the value of ξ + determines whether total extinction occurs. This suggests that there is a bifurcation as V + and ξ + vary, with the system transitioning from a trivial to a non-zero steady state (this behaviour will be investigated in Section 3.3). Figure 5: Series of phase diagrams characterising the steady state distribution predicted by the model as properties of the advection velocity, v z , vary (i.e., for different values of the parameters V + , ξ + and ω + ), and the rate of apoptosis, d f . At each point in (V + , ξ + ) parameter space, we characterise the equilibrium distribution based on the number of peaks and the dominant phenotype (i.e., the z-locations of the local maxima) for different values of the parameters ω + and d f . For parameter sets that give rise to a significant fraction of CSCs (i.e., % CSCs ≥ 1%), the value of φ CSC (0.3, t ∞ ), as defined by Equation (4), is also indicated.
By contrast, the equilibrium distribution for an asymmetric velocity profile (i.e., ω + = 2, as in Figure 3b), has a multimodal distribution, typically characterised by two peaks. In this case, since the CSCs have a lower propensity to mature, they accumulate and persist in the population, even under normoxia. The second column of Figure 5 shows that the proportion of CSCs at long time increases as the death rate, d f , the steepness parameter, ξ + , and the maturation velocity, V + , increase, until the CSCs become the dominant subpopulation (see, for example, Case B.3 in Figure 6). Varying the death rate, d f , does not significantly affect whether extinction occurs; rather, it determines the location of the maximum peak in the equilibrium distribution (see, for example, case B.2 in Figure 6). For low death rates, cells are predominantly in a terminally differentiated state. As the death rate increases, the peak moves to the left, producing an equilibrium distribution in which a higher proportion of rapidly proliferating cells balances the high death rate. Figure 6 shows how the system relaxes to its steady state when ω + = 2. Comparison with Figure 4 reveals that in this case the dynamics are characterised by secondary regrowth, driven by the accumulation of CSCs. For example, in case B.1, phenotypic diffusion enables the cancer cells to de-differentiate, acquire a stem-like phenotype and, therefore, contribute to population growth.  7) and (9), showing how the cell distribution, n(z, t), the phenotypic mean, µ(t), and the cell density, φ(t), change over time. For these results, we use an asymmetric velocity profile (i.e., ω + = 2 in Equation (9)). See Figure 5 for the values of the other model parameters.
To summarise, the properties of the advection velocity v z , determine whether the model predicts extinction or persistence of CSCs, regardless of whether they are present initially. When ω + = 2, random mutations (i.e., diffusion), may dominate the advective force near z = 0, allowing CSCs first to form, then to proliferate and ultimately to comprise a significant proportion of the equilibrium population. CSCs have been observed in normoxic regions; for example, they have been found in perivascular tumour regions, where endothelial cells secrete factors that inhibit CSC maturation [62]. By contrast, when ω + = 1 (i.e., for symmetric velocity profiles), all cells mature over time, leading to the eventual extinction of CSCs. This behaviour could describe that of tumours which lack CSCs, or the effect of drugs which induce stem cell differentiation and, thereby, reduce the incidence of resistance to other treatments, such as radiotherapy. We conclude that targetting V + and ξ + may be effective for eliminating CSCs, increasing tumour sensitivity to treatment and, in certain scenarios, driving tumour extinction.

Hypoxic Conditions
Under hypoxia, the advection velocity in our model is negative and cells will be driven to de-differentiate. In this case, the equilibrium distribution is unimodal, with the dominant phenotype at z = 0. Although varying the death rate d f does not effect the equilibrium distribution (compare cases H3 and H4 in Figure 7), the values of ω − and ξ influence the width of the peak (compare cases H1 and H2 in Figure 7) and, therefore, the variability in the population. Differences in the system dynamics also arise as the initial conditions n 0 (z) vary. The results in Figure 7a indicate little variation in the system dynamics when the initial conditions from Section 3.1 are used. By contrast, in Figure 7b we observe marked differences when the initial conditions are centred around the TDCs. In this case, population regrowth is delayed, the delay depending on the choice of parameter values. For example, when ω − = 2, the velocity in a neighbourhood of z = 1 is smaller than when ω − = 1. Consequently, cells de-differentiate more slowly, delaying tumour regrowth. Similarly, increasing the death rate, d f , reduces the number of cells that can de-differentiate and, subsequently, delays regrowth. Therefore, while d f does not affect the equilibrium distribution, it influences the system dynamics. These results show how the formation of hypoxic regions can shape the development of a tumour. In particular, the emergence of hypoxia maintains and enhances the pool of CSCs, preventing population extinction (see, for example, scenario D in Section 3.1).

Linear Stability Analysis
We now validate some of the above numerical results by performing a linear stability analysis which enables us to characterise the equilibrium states. We denote byn =n(z) a steady state for the (untreated) system (3)-(9), with a total cell densityφ = 1 0n (z)dz and let δn represent a small perturbation to this solution. Then we can approximate the solution n in a neighbourhood ofn as: n(z, t) =n + δn(z, t), δn 1 ∀t > 0.
Substituting this ansatz into (3) and retaining linear terms, we obtain the following equation for δn: where M is the following integro-differential operator The solutionn is spectrally stable if the spectrum of the operator, σ(M), does not contain eigenvalues with positive real part, i.e., Moreover, the dynamics of the system will be dominated by the fastest growing mode (i.e., the eigenfunction corresponding to the eigenvalue with the largest real part, λ 0 ). In Appendix B we transform the above eigen-problem so that it does not include any first order derivatives. For a non-zero steady state, we retain a non-local term in the eigenvalue problem and this can give rise to a spectrum with a pair of complex eigenvalues. Recalling case A.1 from Section 3.1 (see Figure 4), the numerically estimated value of λ 0 is indeed complex (λ 0 = −1.535 × 10 −4 ± i 2.24 × 10 −3 , where i 2 = −1). This result, in turn, explains why damped fluctuations are observed in the numerical simulations.
By contrast, when considering the trivial steady state,n ≡ 0, which is always a fixed point for the system, the non-local term vanishes and we obtain the standard form analysed by Sturm-Liouville theory. Using well known results, we can identify sufficient conditions for the stability/instability of the trivial steady state (see Lemmas 1-3 in Appendix B). Under hypoxia, where v z < 0, we find that the trivial steady state is unstable (for the parameter sets in Table A.4) and the system evolves to a non-zero distribution, which is consistent with the numerical results from Section 3.2. We note that the results relate only to the behaviour of the fitness function and advection velocity near the boundary z = 0, suggesting that the most relevant parameters are p max 0 , V − , θ and ξ − . By contrast, under normoxia, and for the range of parameter considered here, the system undergoes a bifurcation. For sufficiently small V + , the trivial steady state is unstable; for sufficiently large V + and for large values of the death rate, d f , the trivial steady state is stable, (see, for example, case C2 in Section 3.1). To investigate other parameter regimes that we can not tackle analytically, we rely on numerical estimation of the largest eigenvalue, λ 0 . As shown in Figure 8, it is possible to identify the boundary of the region of stability in (ξ + , V + ) space. This diagram does not change significantly as the death rate varies in the range from d f = 0.001 to d f = 0.015 (results not shown). However, the results are highly sensitive to the value of ω + . Comparing Figures 8a and 8b, we see that setting ω + = 2 favours the formation of a non-trivial equilibrium distribution, with the curve shifting to the far right of the parameter space (i.e., The diagrams are obtained for d f = 0.001. We note that changing ω + has a significant impact on the size of the region of (V + , ξ + ) parameter space in which the non-trivial steady state is stable (compare (a) and (b)).
small values of ξ + and large values of V + ). In the latter case, this implies that even higher velocities V + are needed to stabilise the tumour elimination solution. This is consistent with the numerical results in Section 3.1, where setting ω + = 2 (see scenario B in Section 3.1) favoured the accumulation of CSCs which acted a reservoir for tumour cells.

Population Dynamics in the Presence of Treatment
In the previous section, we found that the system possesses a stable steady state to which the dynamics converge for the range of parameter values considered. Therefore, we anticipate that, while treatment can perturb the system from its equilibrium, it will eventually relax to its stable steady state once treatment ends. Thus we expect extinction to occur for parameter values lying in the stability region of the trivial steady state (see Figure 8). From this point of view, we are interested in understanding how different environmental conditions (i.e. normoxia and hypoxia), different treatment protocols and different tumour compositions affect the relaxation phase and, in particular, the time to recurrence.
To account for variability in tumour responses, we consider the different advection velocities used in our earlier analysis (see Table 3). Starting from the initial condition (3f), cells follow different pre-treatment protocols as specified in Table 3. Without loss of generality, we shift time so that t = 0 corresponds to 24 hours before treatment begins. While attention will focus on tumour responses in constant environmental conditions, we also consider briefly treatment responses in changing environments. For each scenario, we simulate the response to treatment for the range of values of the radiation parameters listed in Table 2. We denote by n (S1,R1) (z, t) the solutions corresponding to scenario S1 from Tables 3 and radio-sensitivity parameter set R1 from Table 2.

Treatment Response in Normoxic Conditions
The simulation results presented in Figure 9 illustrate the different regrowth dynamics that can arise when well-oxygenated tumour cells are exposed to a single dose of RT. We identify three distinct behaviours: instantaneous regrowth (S1), decay and extinction (S2) and initial remission with subsequent regrowth (S3). While the cell survival fraction immediately post-treatment depends on the parameter values used in the LQ-model (see Equation (8)), the qualitative population regrowth dynamics post-treatment do not depend on these values.
In more detail, for scenario S1, the cell density increases rapidly after treatment, driving the system towards its (asymptotic) equilibrium. By contrast, for scenarios S2 and S3, the growth curves initially decrease at similar rates until about 40 days after treatment. Thereafter, for scenario S3 the tumour exhibits rapid regrowth to the equilibrium distribution, whereas for scenario S2, the tumour continues to shrink, until it is eventually eliminated.
The origin of such differences can be understood from the time evolution of n(z, t) post-radiotherapy. Figure 9 shows that for case R1 of Table 2, the balance between cell proliferation and advection drives the system dynamics. The reduction in the cell density φ(t) post-radiotherapy reduces intra-population competition and allows the cells to resume proliferation. Depending on the magnitude of the advection velocity (which is positive), the cells either regrow (S3) or they are driven to a terminally differentiated state and, thereafter, become extinct (S2). For scenario S3, the presence of radioresistant CSCs post treatment and a small positive velocity at z = 0 together drive regrowth. As the CSCs start to mature, there is a continuous source of highly proliferative cells which, in turn, drive rapid regrowth of the tumour. As the total cell number increases, intra-population competition slows cell proliferation until eventually advection becomes dominant, driving Scenario Protocol Parameters Subsection     For each scenario S1, S2 and S3 (see Table 3) we consider the dynamics of the total cell number, φ(t), and compare the responses for the radio-sensitivity parameter sets R1, R2 and R3 (see Table 2) to the control, untreated case. For each scenario we also present plots of the phenotypic cell distribution, n(z, t), at different times for radiotherapy protocol R1. The vertical line indicates the time of irradiation, while a line is also shown that follows the evolution of the control (i.e., in the absence of treatment).
the cells to de-differentiate. By contrast, for scenario S2, advection dominates proliferation along the entire phenotypic axis. Additionally, CSCs are absent so that all cells are rapidly terminally differentiated and, thereafter, undergo cell death. Comparison of scenarios S2 and S3 reveals how different phenotypic compositions can generate treatment responses which are initially qualitatively similar, but differ markedly at long times. This finding is reinforced in Figure  10 where we plot the mean phenotypes, µ = µ(t), as defined by Equation (5). For scenarios S2 and S3, the dynamics of the mean phenotype are indistinguishable at short times and do not start to diverge until approximately 20 days after treatment.
More generally, the results presented in Figure 10 reveal three charac- Figure 10: Series of plots showing the evolution of the phenotypic mean, µ(t), for scenarios S1, S2 and S3 (see Figure 9). We note that the scales used on the vertical axes are different.
teristic behaviours for the evolution of the phenotypic mean following radiotherapy. The dynamics of µ may be the same as those prior to treatment, with negligible deviation from the control (see scenario S2). A discontinuity in µ may be induced by radiotherapy (see scenario S1). In this case, CSCs comprise a significant proportion of the population prior to RT and the effect of radioresistance is pronounced (see Figure 9). As CSCs are more likely to survive radiotherapy than more mature cells, we observe an "instantaneous" shift in µ towards less mature phenotypes. The size of the discontinuity depends on the relative sensitivity of CSCs and TDCs to RT, or, using the terminology introduced in Section 2.1, the selective power of RT. Since we are considering high radiation dosages, the discontinuity is determined by the ratio β min /β max . In order for the selective pressure of treatment to be apparent, CSCs must comprise a significant fraction of the population prior to treatment. This explains why, for scenario S3, there is an initial transient period during which, as for scenario S2, there is no discernible deviation from the control. Only at later times does the difference in the evolution of µ(t) for the different parameter sets become apparent. We note that other factors, in addition to stemness, influence cell radiosensitivity. It is natural to expect cells that have permanently exited the cell-cycle will be less radio-sensitive than cycling cells, as the DNA damage response may already be active in such cells [36]. The functional forms for α and β defined by Equations (8b)-(8c) assume that radio-sensitivity increases monotonically with cell phenotype, z. In order to investigate situations in which TDCs have lower radio-sensitivity than proliferating cancer cells, we Figure 11: Series of numerical results showing how the growth dynamics and the phenotypic mean evolves following exposure to a single dose of radiotherapy when cell radiosensitivity is a non-monotonic function of cell phenotype. The simulations are analogous to those presented in Figure 9 and 10, except that Equations (14) are used in place of Equations (8b)-(8c). now the following, non-monotonic functional forms: where H is defined in §2.2, and we arbitrarily fix = 0.075 (all other parameters are as defined in §2.1). When the single dose experiment is repeated with the new radio-sensitivity profile, we observe an overall increase in the population survival fraction (compare Figures 11 and 9) and changes in the dynamics of the population mean µ(t) (compare Figures 11 and 10). The differences are most pronounced for scenarios S2 and S3 where TDCs, localised near z = 1, are dominant in the population prior to treatment. The qualitative growth dynamics (i.e., φ(t)) is similar for both cases. Further investigation of these differences is beyond the scope of the current study and is postponed for future work.
In practice, delivery of a single (high) dose of 10 Gy may not be practical Figure 12: Simulation results for fractionated radiotherapy protocols, showing how the total cell number φ(t) and the phenotypic mean µ(t) evolve for scenarios S1 and S3 (see Figure 9 for details). In all plots, the light purple shaded area indicates the variability in responses when a single dose of 10 Gy is administered and is included for comparison with the fractionated treatments (see Figure 10). The yellow shaded area indicates the duration of the treatment for the fractionated case.
for treating patients, due to adverse side effects [63]. Therefore, we now consider tumour responses to fractionated RT protocols. The trends for fractionated RT are similar to those for single doses for all scenarios in Table  3. Typically, the proportion of cells that survive fractionated therapy is larger than for the single-dose case, by a factor of about 100. Consequently, for scenarios S1 and S3, the time to return to the equilibrium population distributions is reduced. For S2, while treatment causes a monotonic decrease in the cell density φ, since more cells survive fractionated RT, it takes longer for the cell population to become extinct. For scenarios S1 and S3, we recall that for high doses of RT, the phenotypic mean was markedly affected by the specific LQ model parameters considered; this is not the case when lower doses are applied (see Figure 12). Figure 13: Phenotypic distribution n (S1,R1) (z, t) for the control (light blue), the colony exposed to a single dose (dark blue) and the one treated with fractionated dose 2 Gy ×5 (green). The orange and yellow lines indicate the phenotypic mean for the single dose (orange) and fractionated (yellow) therapy respectively. Note that the first panel corresponds to the end of the treatment so that t pt is 24 hr and 120 hr for the 10 Gy and fractionated protocol respectively. On the other hand, the remaining panels are measured relative to the beginning of the treatment, which is at t = 24 hr for both protocols.
The variability in responses for scenarios S1 and S3 following a single dose of radiotherapy can be attributed to the temporary advantage CSCs have post treatment. When using a fractionated protocol, intra-population competition is maintained at the cost of fewer cells being killed. This is apparent when we compare the phenotypic distribution at different times for the two treatment protocols (see Figure 13). When 10 Gy is administered in one dose (first panel, dark blue region), the peak of the distribution is at z = 0. On the other hand, after 5 doses of 2Gy per day (first panel, green region), the proportions of differentiated and cancer stem cells are approximately equal. Given that the former proliferate faster than the latter, the differentiated cells quickly become the dominant phenotype. Consequently, one month after treatment ends (third panel in Figure 13), the proportion of CSCs in the population is the same for both protocols. We conclude further that the single dose protocol outperforms the fractionated protocol when we compare the total number of cells (the blue curve is below the green one for all values of z).

Treatment Response in Hypoxic Conditions
Cell populations that are continuously exposed to hypoxia, exhibit instantaneous re-growth following RT, as shown in Figure 14. Compared with the treatment outcome under normoxia, a higher percentage of cells survive radiation, because there is a larger proportion of radio-resistant cells in the population under hypoxia. Even though a smaller fraction of cells are killed, re-growth is also usually slower under hypoxia than under normoxia. We note also that, following exposure to the single and fractionated protocols, the phenotypic mean µ(t) shifts toward z = 0 under hypoxia, favouring CSCs as the dominant phenotype (see Figure 14). The drift in µ is less pronounced for the fractionated case, suggesting the latter protocol is less favourable for the immediate accumulation of resistant subpopulation of CSCs than the single dose.  Table 3). Simulation results showing the time evolution of the cell density, φ(t), and phenotypic mean, µ(t), are presented. For comparison, the light purple shaded areas in the fractionated plots indicate the variability in the response when a single dose of 10 Gy is administered. The yellow shaded areas indicate the duration of treatment for the fractionated case.
Taken together, our simulation results suggest that, under hypoxia, RT may accelerate the accumulation of resistant cells, while significantly reducing the overall growth rate of the population.

Treatment Response in a Changing Environment
Thus far we have assumed that the oxygen concentration remains constant throughout treatment. While this may accurately describe RT responses for cells cultured in vitro, such control is likely to be absent in vivo [64,65,66]. There is currently no consensus about the impact of RT on tumour vasculature and, hence, tissue re-oxygenation. On the one hand, high doses of radiotherapy may damage the vasculature [67], and decrease nutrient availability post radiotherapy. On the contrary, moderate RT may transiently increase tissue oxygenation by normalising the tumour vasculature (vessel normalisation is a phenomenon that has been observed when tumours are exposed to vascular-targetting agents which destroy some of the blood vessels in a way that increases blood flow through the network and, thereby, tissue oxygen levels [68,69]).
Moreover, as tumour cells are killed, the pressure on immature vessels, not damaged by the radiation, decreases, and oxygen supply to the surviving cells may increase. Equally, hypoxic regions may form at later times as the tumour regrows. From this point of view, radiotherapy may impact both the phenotypic distribution of the cell population (and, thereby, its radio sensitivity), and oxygen levels post-treatment. We can use our mathematical model to investigate these scenarios, by assuming that oxygen levels change post radiotherapy.
Based on the results presented in Sections 4.1 and 4.2, we anticipate that reoxygenation of a hypoxic tumour will be beneficial in certain cases, driving CSC maturation, and even leading to tumour eradication. The results presented in Figure 15 show that the long-term tumour regression is preceded by an initial phase of regrowth during which CSCs that survive treatment de-differentiate and proliferate. Such a treatment might initially be considered unsuccessful, although the stability of the trivial steady state upon re-oxygenation leads to extinction at longer times.
As mentioned previously, when high radiation doses are applied in vivo, it is likely that the vessel network is also damaged, potentially inducing hypoxia [64]. Figure 15b shows that such environmental changes may negatively impact the outcome. The formation of an hypoxic region favours the development and maintenance of radioresistant CSCs, reducing the treatment effi-  Table 3, respectively. Different response to treatment are compared based on parameter values from Table 2.
(c) Growth curve φ(t) and phenotypic mean µ(t) evolution for model R1 from Table 2, when exposed to transient post-treatment hypoxia. We denote by T R the time at which re-oxygenation occur (indicated by the arrows in the plot). If T R is sufficiently small than re-oxygenation does not drive re-growth of the cells population. If we waited for a sufficiently long time (as in case T R = 1000) then re-oxygenation would first drive regrowth. Areas in blue and pink correspond to intervals of normoxia and hypoxia respectively.
cacy and making it more difficult to eradicate the tumour. At the same time, environmental changes may be transient: damaged blood vessels are likely to be replaced by new vessels which form via angiogenesis and re-oxygenate the damaged regions. As shown in Figure 15c, depending on the time-scale required for vessel regrowth (indicated by T R ), different behaviours may arise. If the duration of RT-induced periods of hypoxia is sufficiently short, then the size of the cell population remains low. By contrast, if there is sufficient time for cells to de-differentiate (see T R = 1000), then re-oxygenation leads to a rapid increase in cell number, although eventually the cells die out. These results highlight the complex interplay between tumour growth and treatment response in vivo and the importance of environmental factors in determining the eventual outcome of radiotherapy treatment.

Conclusion and Future Challenges
We have developed a structured model to investigate how clonogenic heterogeneity affects the growth and treatment response of a population of tumour cells. Cell heterogeneity is incorporated via an independent and continuous structural variable which represents stemness. As proposed by [13,19], we view stemness as a plastic trait, with cells becoming more, or less, stemlike depending on their environmental conditions. Our mathematical model accounts for cell proliferation and apoptosis, inter-cell competition, and phenotypic movement along the stemness axis, via diffusion and advection. Studies of the population dynamics in the absence of treatment revealed that, under normoxia, a variety of qualitative behaviours may arise depending on the functional forms used to represent the structural flux and fitness landscape. When advection dominates movement along the stemness axis, its magnitude, relative to the rates of proliferation and cell death, determines whether the population is driven to extinction. Multimodal distributions, which allow for the formation and maintenance of CSCs pools, are observed for asymmetric velocity profiles. Under hypoxia, the population distribution is unimodal and skewed toward stem-like phenotypes, with little intra-population variability. The resulting cell distribution is highly resistant to radiotherapy, the tumour will typically regrow following treatment. By contrast, under normoxia (or re-oxygenated hypoxia), and for suitable parameter values, the tumour may become extinct following radiotherapy.
There are many ways in which the work presented in this paper could be extended. A first, natural extension would be to incorporate structural and spatial heterogeneity (i.e., both phenotypic and spatial dimensions) [21]. This would enable us to consider in vivo situations, where spatial gradients in oxygen levels emerge naturally, due to oxygen consumption by the cells as it diffuses from blood vessels. As outlined in Appendix A, in such a model oxygen consumption rates may vary with cell phenotype, and spatial fluxes may account for random movement of the cells. Preliminary results for such a model are presented in Figure 16. We consider a 1D Cartesian geometry and focus on a tumour region of width L, in which a blood vessel located at x = 0 provides a continuous supply of oxygen to the tissue. If the tumour initially comprises a spatially homogeneous distribution of terminally differ-entiated cells (see Equation (3f)), then the oxygen rapidly relaxes to a steady state and a hypoxic region forms at distance from x = 0. In contrast to the well-mixed model, cells are now able to move, by random motion, between normoxic and hypoxic regions. While terminally differentiated cancer cells are dominant in the well-oxygenated region, a small fraction persists in the hypoxic region (in particular, near the boundary of the hypoxic region, orange line in the plots in Figure 16). This is due to the influx of cells from the well-oxygenated portion of the domain. Similarly, CSCs are dominant in the hypoxic region, but a small fraction of hypoxic CSCs migrate towards x = 0, where re-oxygenation induces their maturation, creating a differentiated and highly proliferative cell phenotype, alongside terminally differentiated cancer cells. These results illustrate how the interplay between space, resources and phenotypic adaptation may give rise to complex behaviours; their investigation is the focus of ongoing work. Figure 16: Series of plots showing how, in the absence of treatment, the cancer cell population n(x, z, t) and the oxygen concentration c(x, t) change over time t when we account for spatial and phenotypic variation (see Equations (A.1)). We indicate the threshold c = c H which defines the boundary of the hypoxic region with a horizontal red line in the upper plots and with a vertical orange line in the lower plots. We fix V ± = 4 × 10 −4 , ξ ± = 0.1, ω + = 1, ω − = 2 and d f = 0.001, while the remaining model parameters are fixed at the values stated in Table A.

4.
A significant challenge of the modelling approach presented in this paper is the determination of model parameters and functional forms. In the longer term, techniques such as single-cell RNA sequencing [17,18] will make it be possible to quantify specific aspects of our model, such as the dependence of the proliferation and apoptosis rates on cell stemness and the dependence on the tumour micro-environment of the (phenotypic) advection velocity associated with cell maturation and de-differentiation. In spite of their current limitations, we believe that studies of such models can increase understanding of the ways in which specific physical processes may influence the phenotypic distribution of cell populations in different environments. At the same time, we acknowledge that it remains a matter of debate as to whether asymmetric cell distributions are driven by micro-environmental signals (as in the model presented here), asymmetric division, or a combination of the two [70]. By using a non-local proliferation kernel to account for asymmetric division, we could investigate these alternative hypotheses and identify conditions under which they lead to different outcomes.
A important feature of our model is the way in which the response to radiotherapy (RT) varies with cell stemness (i.e., z). Our analysis shows how the functional forms used to describe the advection velocity and fitness functions can affect the system dynamics post-RT. While unimodal phenotypic distributions lead to monotonic growth curves post-treatment, more complex behaviour is observed when heterogeneous populations, with a pool of CSCs, are considered. For example, under normoxia, the presence of radio-resistant CSCs can drive recurrence, despite an initial phase of tumour regression. As the CSCs mature into highly-proliferating cancer cells, rapid re-growth is accompanied by re-sensitisation of the population to RT. Under hypoxia, CSCs maintain their stemness, leading to a slowly growing, radio-resistant cell population. More complex outcomes arise when we consider the effect that treatment might have on the environment. As noted in 4.3, changes in the vasculature induced by radiotherapy can result in either post-treatment reoxygenation or hypoxia. While re-oxygenation increases the radio-sensitivity of the population, hypoxia increases their radio-resistance. In practice, such environmental changes are likely to be transient. Even in an untreated tumour, fluctuations in oxygen levels can occur. Consider, for example, cells in a neighborhood of immature blood vessels. As the cells proliferate, they exert mechanical pressure on the vessels, causing them to collapse and local oxygen levels to fall. Under hypoxia, the tumour cells stimulate the growth of new blood vessels from pre-existing ones, via angiogenesis. In this way, tumour regions may cycle between periods of hypoxia and normoxia. It would be of interest to extend the model to account explicitly for the tumour vasculature and its interaction with tumour cells. This could be achieved at a "high level" of description, via simple ODE models such as [71,72], or via more complex, multi-phase [73] or multi-scale approaches [74,75,76,77].
This would enable us to better capture the different time-scales on which the oxygen dynamics and cell adaptation velocity change. As shown in Figure 17, variations in oxygen levels emerge naturally within spatially-resolved models. Here, cell killing leads to tissue re-oxygenation which, in turn, disrupts the CSC niche. Depending on the time scale over which the cells adapt to their new environmental conditions, this may increase the overall radiosensitivity. Understanding and accounting for such phenomena is particularly relevant for predicting responses to RT and comparing alternative treatment protocols. Figure 17: Evolution of the population n(x, z, t) in the spatial and phenotypic dimensions following a cycle of fractionated radiotherapy (5 × 2 Gy). The parameter values are the same as those used in Figure 16 and the initial cell distribution is the same as the final distribution in Figure 16. For the LQ-model we used parameter set R3 in Table 2.
In the extinction scenario, or post administration of high RT doses, the number of cells in the population can become low and our continuum model may cease to be valid. In such conditions, stochastic effects which are neglected herein may become important. As in [78,79,80], stochastic and mean field approaches may be combined with hybrid discrete-continuum techniques to account for small population effects and to study their impact on the probability of tumour extinction.
In this paper, we considered only single dose and fractionated treatment protocols. In future work, we could investigate alternative strategies, such as adaptive therapeutic protocols [81] and/or multi-drug treatments, which have been proposed as an effective way to overcome radio-resistance. From this point of view, considerable efforts have been invested in designing treatments that exploit features of CSCs, such as their metabolic plasticity [10]. Motivated by recent metabolically-structured models [24,21,25], a natural extension of our model would be to include a "metabolic dimension" in order to investigate the interplay between stemness, metabolic switching and resistance. A biologically informed model that incorporates metabolic and phenotypic effects, together with the tumour micro-environment and vascular remodelling lies at the heart of a mathematical program that would enable systematic comparison with in vivo observations. The framework and results outlined in this work represent a first step towards achieving this long-term goal.

Appendix A. Spatial Model
We outline here the set up for the 1D simulations presented in Section 5. As a full description of the spatial model goes beyond the scope of the present work, we focus on the main changes to (3)- (9). We now view the oxygen concentration c as a dependent variable, rather than a prescribed function. We suppose that oxygen is supplied to the region by blood vessels on the domain boundary ∂Ω 2 (see Figure A.18). Oxygen diffuses from the boundary into the tissue where it is consumed by the tumour cells at rates which depend on their phenotype and the local oxygen concentration. The evolution of the dimensionless cell density, n = n(x, z, t), is driven by a phenotypic flux of the same form as in Equation (3) but a spatial flux is included to account for random motion in the spatial dimension.  Figure A.18, we consider a fixed tissue slice where the oxygen supply (i.e. vasculature) is confined to one of the tissue boundaries. Given the assumed symmetry of the problem, we can consider a 1D Cartesian geometry with x ∈ [0, L]. The spatial model is defined by the following system of coupled PDEs:

As shown in
In Equation (A.1), D N and D C are the assumed constant spatial diffusion coefficient for the cells and oxygen, respectively, while γ denotes the rate at which cells of phenotype z consume oxygen and Γ the net rate of oxygen consumption at position x and time t. The advection velocity v z is as defined by Eq. (9), while the fitness function F is analogous to that defined in Section 2.1, with an additional term to account for necrosis. The latter is assumed to occur at a constant rate g ≥ 0, independent of cell phenotype, when the oxygen concentration falls below a threshold value, c N ≥ 0. We also modify the definition of the survival fraction SF given in §2.1 (see Equation (8a)) to account for the oxygen-enhancement ratio (OER) [34,33]. According to the oxygen fixation hypothesis [82], part of the biological damage induced by radiation is indirect, being mediated by the presence of free radicals. Thus, when oxygen is limited, radio-sensitivity is accordingly reduced. Based on experiments, the range of oxygen concentrations at which this effect is relevant corresponds to more severe levels of hypoxia (where c ∼ 0.5% or lower). We do not consider such situations for the well-mixed model, where we consider mild hypoxia. However, accounting for the OER will be important for the spatially extended model. Recall from Section 3.2 that hypoxia is a favourable niche for CSCs. Therefore the OER will endow them with additional protection from radiation. Denoting by c R H the oxygen threshold at which the OER becomes active, we use the following functional form for the survival fraction when simulating the spatially-extended model: In Equation (A.2), α and β are defined by (8). We note that in the main text, we consider c = 1 (normoxia) and c = 0.2 (hypoxia), so that the OER does not impact cell responses to RT. For the well-mixed model, the oxygen concentration is typically maintained at a prescribed, constant value. By contrast, for the spatially extended model, we suppose that the tumour cells consume oxygen at a rate γ which depends on their phenotype, z. As mentioned previously, stem cells are known to have a glycolytic metabolism and, thus, we assume that they consume less oxygen than cancer cells. Consequently, we consider γ to be a monotonically increasing function of the phenotypic variable z which asymptotes to its maximum value for z > 0.5: In   (6b). As most data reported in the literature refer to processes, such as cell proliferation, at the population/cell colony and do not account for phenotypic variation, it was difficult to estimate parameters that characterise the phenotypic variation in these processes. We based our estimates of the proliferation rate on the doubling times reported by [84] for two breast cancer cell lines, MCF-7 and BT-549. The former belong to the class of laminal -like cells which are characterised by low stemness levels [89] and high proliferation rates (doubling time 1.8 days, i.e., growth rate 0.016 hr −1 ). On the other hand, BT-549 belong to the class of triple-negative cells whose population is dominated by highly aggressive but slowly proliferating stem-like cells [89] (doubling time 3.7 days [84], i.e., growth rate 0.008 hr − 1). Given the variability in the phenotypic distribution of these cell lines, we have rounded the values to those presented in Table  A.4. As is common in the literature, we have chosen the source of oxygen (i.e. c ∞ ) to be at a pressure of 100 mmHg [87]. Given that atmospheric pressure corresponds to 760 mmHg with 21% O 2 , the oxygen tension corresponding to c ∞ is about 8% O 2 . The hypoxic and necrotic thresholds (c H and c N ) are then equivalent to oxygen pressures of 2.5% O 2 and 0.1% O 2 in line with [90,88]. These values can be converted into oxygen concentrations by use of Henry's law [87], see Table A.4. where q(z; c,φ) = p(z; c) ( and k(s, z) = exp 1 2θ We conclude that the trivial steady state is either a stable node (if λ 0 < 0) or a saddle (if λ 0 > 0). In addition to numerical estimation of λ 0 , analytical approximations and bounds can be obtained via the so-called Rayleigh quotient R(y). If we multiply Equation (B.6a) by y and integrate by parts, then we obtain: where y also satisfies the Neumann boundary conditions B.6b. We deduce that the following therefore holds: where E is the set of twice differentiable functions that satisfy condition (B.6b).

Lemma 1.
If the functionq is such that max z∈(0,1)q < 0 then the null steady state is stable.
Proof. Consider the numerator of the quotient defining R(y): We deduce that It is therefore apparent that if the function q is negative throughout the domain, then R up is negative for any choice of y ∈ E. In such a case, we have that: We now show that under normoxia, q < 0 if the death rate is high and the magnitude of the phenotypic advection velocity is sufficiently large.
Lemma 2. If the model proliferation rate, apoptosis rate and phenotypic advection velocity and diffusion coefficient are chosen such that: then the trivial steady state is unstable.
Remark. Note that for (B.10) to hold we require 1 0 (p − f )dz > 0 so that cell proliferation dominates apoptosis. Based on the functional form defined in Section 2.1, we have that: where Z denotes the cumulative distribution function for the normal distribution. We note that I(1; d f ) > 0 while I(0.2; d f ) < 0 for all values of the parameters listed in Table 1. We conclude that under normoxia there is a threshold V + (ξ + , ω) such that the system is unstable for all choices of V + < V + (ξ + , ω): where I v (ξ + , ω + ) = We note also that higher values of θ favour instability of the trivial solution as V + increases with θ. By inspecting Figure 3, we note qualitatively that I v is expected to decrease for increasing values of ξ + and ω + .
To analyse other regions of parameter space, where neither of the sufficient conditions holds, we rely on numerical estimates of the eigenvalue λ 0 . As shown in Figure B.19, and as expected based on the above findings, when the magnitude of the velocity V + is small, λ 0 > 0 for all ξ and the trivial solution is unstable. By contrast, as the magnitude of the advection velocity increases, its steepness, ξ, determines the stability of the trivial solution. Using this estimate, we can identify the region of stability of the trivial steady state (see Figure 8 in Section 3.3). We remark that the boundary between the regions of stability is non-smooth. This is because λ 0 = λ 0 (ξ) plateaus as ξ 1 (see Figure B.19). By computing the second largest eigenvalue, λ 1 (ξ), we observe that the sharp change in the profile of λ 0 as ξ decreases occurs where |λ 0 − λ 1 | attains its minimum value. It is possible to show that the two eigenvalues do not cross, as expected by the Sturm Oscillation theorem. A similar phenomenon occurs in quantum physics [92] where it is known as avoided crossing.
Finally, we consider the stability of the trivial solution in an hypoxic environment. We confirm the numerical simulations from §3.2 by showing that, under hypoxia, the trivial solution is always unstable.  Table A.4, the trivial steady state is always unstable.
Proof. Let us consider as a trial function: where a small parabolic correction is added to the standard Gaussian, the constant A being chosen to ensure that the boundary condition at z = 1 is satisfied: A = e − 1 2κ 2 2π 1/4 κ 5/2 ; (B.14) the derivative y at z = 0 vanishes, by construction. We now want to show that the Rayleigh quotient is positive for such a choice of the test function y which implies that the trivial steady state is unstable.
Given that the denominator of R(y) is always positive, its sign will be determined by the numerator R n (y) that is: Computing the derivative of y and denoting the Gaussian by y 0 , we obtain: y 2 = y 2 0 + 2Az 2 y 0 + A 2 z 4 , (B.16a) Recalling that the constant A is exponentially small in κ while y 0 grows only as a power law of κ −1 , the terms multiplied by A will be negligible and the sign of R n (y) will be determined also by the leading term: Proving instability therefore reduces to show that I 0 is positive for the range of parameters and functional forms considered in hypoxic condition. We do so finding a lower bound on the value on I 0 , exploiting the quick decay of the function y 0 , whose mass is concentrated in a neighborhood of z = 0. Given that p(0) − f (0) > 0 and m(0) = 0, provided that m does not grow too fast near z = 0, we can intuitively see that the major contribution to the integral I 0 will be positive. We will now expand this intuitive argument with a more rigorous calculation. We first focus on I (1) 0 = 1 0 (p − f )y 2 0 dz, the contribution in (B.16e) due to cell proliferation. We can compute this integral exactly as the integrand comprises products of exponentials, that can be re-written as the integral of Gaussian distribution: Let us reiterate that we want z 0 and z 1 to be such that m 2 (z 0 κ) and m 2 (z 1 κ) are not too large while [Z] √ 2z 0 0 and [Z] ∞ √ 2z 1 are sufficiently small. In this way, the growth of m is balanced by the exponential decay of the Gaussian function y 2 0 . In particular, we choose z 0 = √ 2 and z 1 = 5/ √ 2. Combining the above with the estimate from Equation (B.17), we obtain: magnitude V − and steepness ξ − considered in the paper (without loss of generality, we only consider d f = 0.015 as I low 0 decreases with d f ). As shown in Figure B.20, for all such values, we have that I low 0 > 0. As I 0 > I low 0 , we therefore have that generically I 0 is also positive. We estimate A ≤ O(10 −13 ) which justifies us dropping the O(A) in (B.16d). Consequently, we conclude that R n (y) is positive and so is the quotient R. Hence, in hypoxia, the trivial steady state is always unstable.