Phenotypic variation modulates the growth dynamics and response to radiotherapy of solid tumours under normoxia and hypoxia
Giulia L. Celora, Helen M. Byrne, Christos Zois, Panos G. Kevrekidis
PPhenotypic variation modulates the growth dynamicsand response to radiotherapy of solid tumours undernormoxia and hypoxia
G.L. Celora , H. Byrne a , C. Zois b , P.G. Kevrekidis c a Mathematical Institute, University of Oxford, Oxford, UK b Molecular Oncology Laboratories, Department of Oncology, Oxford University,Weatherall Institute of Molecular Medicine, John Radcliffe Hospital, Oxford, UnitedKingdom. c Department of Mathematics & Statistics, University of Massachusetts, Amherst 01003USA
Abstract
In cancer, treatment failure and disease recurrence have been associated withsmall subpopulations of cancer cells with a stem -like phenotype. In this pa-per, we develop and investigate a phenotype-structured model of solid tumourgrowth in which cells are structured by a stemness level, which varies con-tinuously between stem-like and terminally differentiated behaviours. Cellevolution is driven by proliferation and apoptosis, as well as advection anddiffusion with respect to the stemness structure variable. We use the modelto investigate how the environment, in particular oxygen levels, affects the tu-mour’s population dynamics and composition, and its response to radiother-apy. We use a combination of numerical and analytical techniques to quantifyhow under physiological oxygen levels the cells evolve to a differentiated phe-notype and under low oxygen level (i.e., hypoxia) they de-differentiate. Un-der normoxia, the proportion of cancer stem cells is typically negligible andthe tumour may ultimately become extinct whereas under hypoxia cancerstem cells comprise a dominant proportion of the tumour volume, enhancingradio-resistance and favouring the tumour’s long-term survival. We then in-vestigate how such phenotypic heterogeneity impacts the tumour’s responseto treatment with radiotherapy under normoxia and hypoxia. Of particularinterest is establishing how the presence of radio-resistant cancer stem cells [email protected] Preprint submitted to Journal of Theoretical Biology January 15, 2021 a r X i v : . [ q - b i o . CB ] J a n an facilitate a tumour’s regrowth following radiotherapy. We also use themodel to show how radiation-induced changes in tumour oxygen levels cangive rise to complex re-growth dynamics. For example, transient periods ofhypoxia induced by damage to tumour blood vessels may rescue the cancercell population from extinction and drive secondary regrowth. Further modelextensions to account for spatial variation are also discussed briefly. Keywords: cancer stem cells, phenotypic variability, radio-resistance
1. Introduction
Understanding of the mechanisms by which cancer is initiated and progressescontinues to increase, and, yet, cancer remains one of the leading causes ofpremature mortality worldwide and a major barrier to increasing averagelife-expectancy. For example, in 2018, 9.6 million people are estimated tohave died of cancer [1]. Furthermore, treatment outcomes can differ markedlybetween patients with the same cancer type, with the emergence of resistancebeing one of the major causes of treatment failure.Over the past twenty years, there has been a major shift in our perceptionof solid tumours; they are now regarded as heterogeneous tissues in whichmalignant cells interact with normal cells and shape their environment inways that favour malignant growth [2]. Cancer stem cells (CSCs) were intro-duced to explain intra-tumour heterogeneity via the
CSC hypothesis [3]. Thishypothesis proposes that, while CSCs may comprise only a small fraction ofthe total cell population, their high clonogenic potential and their ability toproduce more mature, or specialised, cancer cells enables them to create anentire tumour [4]. As CSCs are found to be resistant to standard treatments,they are recognised as a major cause of disease recurrence and treatment fail-ure [5, 4, 6]. These observations have stimulated the development of noveltherapeutic strategies which aim to eradicate CSCs [7, 8, 9, 10]. In practice,the plasticity of CSCs represents a major obstacle to such treatments. Ad-ditionally, CSCs can adapt to their local micro-environment, and remodel itto create and maintain a niche which supports their survival [11].Increasingly, researchers are turning to mathematical models in order tounderstand how CSCs affect the growth and composition of tumours, par-ticularly their heterogeneity and response to treatment. These models oftendecompose the tumour into a series of compartments, each representing aparticular cell subtype. For example, in [7], Enderling distinguishes cancer2tem cells (CSCs) and cancer cells, whereas Saga and coworkers distinguishradio-resistant and radio-sensitive cells [12], and Scott and colleagues distin-guish tumour-initiating cells (or CSCs), transit-amplifying cells and termi-nally differentiated cells (TDCs) [13]. Thus, most compartmental models arebased on the CSC hypothesis which assumes that it is possible to distinguishbetween cancer stem cells and the tumour bulk. However, this paradigm hasbeen challenged by recent experimental studies [14, 15] that highlight thephenotypic 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 intra-tumoural heterogeneity, based on adaptive CSC plasticity [16]. Under thishypothesis, cancer cells move between stem-like and terminally differentiatedstates in response to extrinsic (environmental) and/or intrinsic (random epi-genetic mutation) forces. Remarkably, the development of state-of-the-artexperimental tools, such as single-cell RNA-seq, means that it is now possi-ble to track the evolution of stemness traits [17, 18], rendering this an idealtime 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, sincethey assume that the tumour comprises distinct cell populations, with dis-tinct properties, they are unable to account for continuous variation in cellproperties. An increasingly popular mathematical approach for describingpopulation heterogeneity and plasticity characterises tumour cells by theirposition on a continuous phenotypic axis. Position on the phenotypic axisdetermines cell properties such as resistance to treatment [19, 20, 21, 22, 23]and/or metabolic state [24, 25]. This approach is motivated by concepts fromevolutionary ecology, such as risk-spreading through spontaneous (epigeneticor genetic) variations and evolutionary pressure [26]. The resulting modelsare typically formulated as systems of reaction-diffusion equations [24, 22, 25],with an advective transport term sometimes included to account for biasedmutation dynamics [21] or adaptive phenotypic switches [19, 27, 28].In this paper, we formulate a mathematical model that accounts forthe evolution of a cancer cell population along such a stemness axis in re-sponse to extrinsic and intrinsic stimuli. Initially, we focus on the plasticresponse of cells to changes in nutrient levels, in particular oxygen. Thisis motivated by recent experimental studies [29, 30, 31, 32] suggesting thathypoxia (i.e. low oxygen levels) is a key driver of cell de-differentiation.From this point of view, spatial heterogeneity may introduce significant ad-3itional complications: as oxygen diffuses into a tumour and is consumedby cells, spatial gradients in the oxygen levels are established. In this way,local micro-environments characterised by normoxia, hypoxia and necrosisform as the distance to the nearest nutrient supply (i.e., blood vessels) in-creases [21, 23, 25]. For simplicity, we postpone consideration of such spatialcomplexity to future work and focus, instead, on a well-mixed setting whereoxygen levels are homogeneous and prescribed. This idealised scenario al-lows us to investigate how cell properties, such as proliferation, apoptosisand adaptive response to environmental signals, contribute to the emergenceof heterogeneous stemness levels in the population and the long term tumourcomposition. In this regard, we are interested in identifying conditions underwhich CSCs are favoured. We then extend the model to account for treat-ment via a phenotypically-modulated linear-quadratic model of radiotherapy(see, e.g., [33, 34, 12] for recent discussions) which accounts for differentialradio-sensitivity of CSCs [10]. This allows us to investigate how differentradiotherapy protocols perturb the phenotypic distribution and subsequentregrowth of the tumour.In practice, stemness is just one of multiple traits that regulate cell be-haviour and heterogeneity. We, therefore, anticipate that future models willcombine multiple phenotypic axes or synthetic dimensions , such as stemnessand metabolic state [24, 21]. Given the complexity of such multi-dimensionalmodels, it is important first to understand these aspects separately. Notingthat considerable mathematical effort has been devoted to investigating can-cer metabolism [35], we choose here to focus on population heterogeneity withrespect to a continuously varying stemness axis. We hope that in the longterm this work will help motivate a systematic experimental characterizationof cell plasticity and phenotype.The remainder of the article is organised as follows. In Section 2, wepresent a well-mixed, spatially homogeneous, model of solid tumour growthin response to a prescribed oxygen concentration. We first investigate thepopulation dynamics in the absence of treatment, considering both normoxicand hypoxic conditions. Numerical results are presented in Section 3. As apartial validation of the numerical results, we use spectral stability analysisto characterise the long time behaviour of the solutions. Section 4 focuses ontumour cell responses to different radiotherapy protocols. As in Section 3,we simulate responses under normoxia and hypoxia, but we also considersituations in which the environment alternates between periods of hypoxiaand normoxia in order to explore the different ways that radiotherapy can4lter tissue oxygenation. Finally in Section 5, we summarize our key findingsand propose possible directions for future work. We also present preliminaryresults showing how accounting for spatial and phenotypic variation mayaffect a tumour’s growth and response to radiotherapy.
2. Model Formulation
We consider the temporal evolution of a heterogeneous population of tumourcells, N ( z, t ), where t ≥ z (0 ≤ z ≤
1) represents theirstemness or clonogenic capacity . As shown in Figure 1, z = 0 correspondsto cancer stem cells (CSCs) which have the maximum level of stemness, and z = 1 corresponds to terminally differentiated cells (TDCs), which have losttheir proliferative capacity and which can either enter replicative senescenceor undergo cell death [36]. We assume that the population dynamics mayby 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], whichgive rise to advective and diffusive fluxes respectively. Finally, the effect ofnatural 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 associatewith each cell a stemness level z , which varies continuously between the cancer stem cellstate (CSCs, with z ∼ z ∼ .
5) and the terminallydifferentiated cell state (TDCs, with z ∼ While multiple nutrients and growth factors regulate the growth rate (orfitness 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 , incancer has long been recognised due to its association with cell quiescence5nd poor therapeutic outcomes [21, 34, 12]. Recent experimental results [15]have shown that hypoxia also plays a role in de-differentiation by regulat-ing pathways associated with a stem-like phenotype. We account for thesephenomena in our model by assuming that all cells are exposed to the samelevel 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 fromexisting theoretical models in which intrinsic forces are assumed to domi-nate phenotypic variation (i.e., v z = 0) [24, 25]. By combining the processesmentioned above, we deduce that the evolution over time t and along thephenotypic axis z of the cell concentration, N ( z, t ), is governed by the fol-lowing non-local partial differential equation (PDE) and associated boundaryand initial conditions: ∂N∂t = ∂∂z (cid:18) θ ∂N∂z − N v z ( z, c ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) structural flux + F ( z, Φ , t ; c ) (cid:124) (cid:123)(cid:122) (cid:125) fitness N, (1a) θ ∂N∂z − N v z = 0 , z = { , } , t > , (1b) N ( z,
0) = N ( z ) z ∈ (0 , , (1c)Φ( t ) = (cid:90) N ( z, t ) dz. (1d)In Equation (1), the non-negative constant θ represents the rate at whichcells diffuse along the phenotypic axis, due to random epigenetic mutations,Φ( t ) denotes the density of cells in the domain at time t , and N ( z ) is theinitial distribution of cells along the phenotypic axis. In ecology, the function F is referred to as fitness landscape which is a mathematical representationof natural, or Darwinian , selection [37]. We suppose it has the followingform: F ( z, Φ , t ; c ) = p ( z, c ) (cid:18) − ΦΦ max (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) proliferation − natural celldeath (cid:122)(cid:125)(cid:124)(cid:123) f ( z ) − M (cid:88) i =1 log (cid:18) SF ( z, c ) (cid:19) δ ( t − t i ) (cid:124) (cid:123)(cid:122) (cid:125) radiotherapy . (1e)In Equation (1e), p = p ( z, c ) denotes the phenotype-dependent growth rateof the cells (see Section 2.1 for details). It is multiplied by a non-local (in6he phenotypic sense) logistic term, with constant carrying capacity Φ max ,to capture intra-population competition for space and other resources. Weassume that oxygen levels remain sufficiently high so that necrosis can beneglected. Hence, the death rate, f , accounts only for natural cell death,or apoptosis, which is assumed to occur at a rate which is independent ofthe oxygen concentration, c ( t ). Radiotherapy (RT) also contributes to celldeath and, in so doing, reduces cell fitness. We suppose that M rounds of RTare administered at discrete times t i ( i = 1 , , . . . , M ). After each treatmentdose, the proportion of cells of phenotype z that survive is denoted by thesurvival fraction SF ( z, c ). By allowing SF to depend on z , we can account forphenotypic-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 differentialradio-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 cellsfollowing 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: n = N Φ max , φ = ΦΦ max , (2)where the units of time, t [hr] are preserved in a dimensional form to facilitatethe interpretation of the results. Under this rescaling, equations (1) become ∂n∂t = ∂∂z (cid:18) θ ∂n∂z − nv z ( z, c ) (cid:19) + F ( z, φ, t ; c ) n, (3a) θ ∂n∂z − nv z = 0 , z ∈ { , } , t > , (3b) n ( z,
0) = n ( z ) z ∈ (0 , , (3c) φ ( t ) = (cid:90) n ( z, t ) dz, (3d) F ( z, φ, t ; c ) = p ( z, c ) (1 − φ ) − f ( z ) − N (cid:88) i log (cid:18) SF ( z, c ) (cid:19) δ ( t − t i ) . (3e)In order to complete the model, it remains to specify several functionalforms; this will be done in Sections 2.1 and 2.2. Extending the model to7ccount for spatial variation is presented in Appendix A, and preliminaryresults are included in Section 5 (a full investigation of the spatially-extendedmodel is postponed to future work).In what follows, we assume that oxygen concentration c has been rescaledso that c = 1 corresponds to physiological oxygen levels, namely physoxia ,which is about 8% oxygen [40]. When considering hypoxia, we focus on mildhypoxia, fixing c = 0 . .
6% oxygen in standard units(see Appendix A.1 for details). At this oxygen concentration, necrosis can beneglected; it typically occurs at lower oxygen tensions (approximately 0 . n ( z ) = φ √ πσ e − ( z − . σ , (3f)where the positive constants φ and σ specify the initial size and phenotypicvariance of the population.The proportion of CSCs is often used to characterise heterogeneous pop-ulations of cancer cells. CSCs are typically identified by their expression ofspecific markers (such as CD44/CD24 and ALDH1, depending on the tu-mour type [10]); thresholds in these markers are used to distinguish stemfrom differentiated cancer cells. Since our model treats stemness as a con-tinuously varying cell property, we introduce a threshold z ∗ ∈ (0 ,
1) in oursimulations, and classify cells with 0 < z < z ∗ as CSCs. We therefore definethe proportion of stem cells at time t to be: φ CSC ( t, z ∗ ) = (cid:82) z ∗ n ( t, s ) dsφ ( t ) . (4)As a further statistical feature of the cell population, we introduce the phe-notypic mean, µ ( t ), which is defined as follows: µ ( t ) = 1 φ ( t ) (cid:90) zn ( z, t ) dz. (5)In the absence of suitable experimental data, it is difficult to specify manyof the parameters and functional forms in Equations (3). For this reason,we focus on identifying the qualitative behaviours that the model exhibitsacross a range of ‘biologically-reasonable’ situations.8 .1. 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 maximacentred around z = 0 and z = 0 .
55, representing respectively cells withstem-like ( z = 0) and intermediate phenotypes ( z = 0 .
55, this value beingarbitrary). As shown in Figure 2, this choice reduces the overlap of the twoGaussian profiles while maintaining the proliferation rate at z = 1 close tozero. This asymmetry also emphasises that, under normoxia, more stem-likecells (i.e. z < .
5) proliferate at lower rates than more differentiated cells (i.e. z > . niches each of which will favour a particularphenotype. We account for this effect by assuming that the amplitude of thepeaks in the proliferation rate are oxygen dependent. Accordingly we write: p ( z ; c ) = p ( c ) exp (cid:20) − z g (cid:21) + p ( c ) exp (cid:20) − ( z − . g (cid:21) , (6a) p i ( c ) = p maxi c K i + c , i = 0 , , (6b)where p ( c ) and p ( c ) are Hill–Langmuir type equations with fourth orderexponents, so that the growth rate decays rapidly when c ∼ K i . We assumethat differentiated cells are fitter than CSCs under normoxia and, therefore,choose p max > p max . At the same time, we note that chronic hypoxia iswidely considered to favour CSCs [41, 42, 43]. The plasticity of CSCs en-ables them more easily to adapt their metabolism to changing nutrient levelsthan differentiated cells [29, 44] and, therefore, to survive and proliferate inchallenging conditions. This behaviour contrasts with that of differentiatedcancer cells which tend to become quiescent when exposed to hypoxia. Weaccount for these effects by assuming K (cid:28) K .When we consider the rate of cell death due to apoptosis, f ( z ), we notethat apoptosis occurs predominantly when cells lose their clonogenic capac-ity. As such, it predominantly affects only TDCs with z ∼
1. Motivatedby the mathematical models developed in [7, 13], we propose the followingmonotonically increasing function for f ( z ): f ( z ) = d f e − k f (1 − z ) . (7)Even though they may not proliferate, TDCs compete for space and resourcesand, thus, impact the tumour dynamics. In what follows, we consider two9ifferent cases. First, guided by experimental results reported by Driessens etal. [45], we assume that apoptosis of TDCs occurs on a much longer timescalethan that on which cells proliferate so that d f << max z p ( z ; 1). In thesecond case, the rates of cell proliferation and apoptosis are assumed to becomparable. This situation represents a tumour with high cell turnover and,as we will see, gives rise to a tumour population with higher clonogeniccapacity.In Figure 2, we sketch the fitness landscape F ( z, , t ; c ) for different envi-ronmental conditions in the absence of treatment and competition. In doingso, we have neglected competition and radiotherapy in Equations (3e), where p and f are defined by Equations (6)-(7). (a) high d f (b) low d f (c) low d f Figure 2: Series of sketches showing how the maximum growth rate p ( z, c ) − f ( z ), as definedby Equations (6)-(7) changes in different micro-environments: (a)-(b) under normoxia( c = 1), the progenitor cells ( z = 0 .
55) are the fittest phenotype, and the death ratemay be either high (a) or low (b); (c) under hypoxia ( c = 0 . z = 0) are thefittest phenotype. The parameter values used to produce the sketches are listed in Table 1.Regions of positive and negative fitness are highlighted in green and red, respectively. p maxi ( hr − ) K i g i i =0 0 .
005 0 .
05 0 . i =1 0 .
02 0 . . (a) proliferation d f ( hr − ) k f { . , . } (b) apoptosis Table 1: Range of parameter values used in the sensitivity analysis. More information onthe specific parameter choice can be found in Appendix A.
We now consider the impact of radiotherapy on cell fitness. As mentionedabove, CSCs possess protective mechanisms that enable them to withstand10amage caused by radiation and oxidative stresses [46, 47, 48, 4, 6, 10, 49].They are, therefore, more resistant to treatment than their differentiatedcounterparts. It is well known that local oxygen concentration levels alsoaffect treatment outcomes [50, 51]. While we account for this effect in thefull spatial model (see Appendix A), here we focus on the role of phenotype-dependent radio-sensitivity. In particular, we adapt the standard Linear-Quadratic (LQ) model so that the tissue specific coefficients, α ( Gy − ) and β ( Gy − ), are phenotype dependent: − log( SF ) = α ( z ) d + β ( z ) d , (8a)where d is the radiation dose in Grays (Gy). Equation (8a) is the natural,continuum extension of previous works [52, 12], in which two-compartmentmodels are used to describe the time-evolution of cancer cells and cancer stemcells exposed to radiotherapy, and CSCs are assumed to be radio-resistant.Accordingly, here, we assume α and β are increasing functions of the pheno-type z [12, 6, 10] of the following form: α ( z ) = α min + ( α max − α min ) tanh (cid:18) zξ R (cid:19) , (8b) β ( z ) = β min + ( β max − β min ) tanh (cid:18) zξ R (cid:19) . (8c)In Equations (8b)-(8c), ξ R , α min,max and β min,max are non-negative con-stants with α min < α max and β min < β max . Where possible, parameterestimates are taken from the literature (see [12] for estimates of α min,max and β min,max ); the value of ξ R = 0 . z > .
5) have maximum sensitivity to treatment (i.e., α ( z ) ∼ α max for z > . z ∼
0) respond in the same wayto RT, whereas differentiated cancer cells (with z > .
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 valueof α max /β max corresponds to an early responding tissue, with a low repaircapacity, for which fractionation is known to be beneficial [53]. Finally, caseR2 is intermediate between cases R1 and R3. By assuming heterogeneity in11 α min , α max ]( Gy − ) [ β min , β max ]( Gy − ) α min β min ( Gy ) α max β max ( Gy )R1 [0 . , .
15] [0 . , .
10] 2.5 1.5R2 [0 . , .
20] [0 . , .
05] 2.5 4R3 [0 . , .
40] [0 . , .
05] 2.5 8
Table 2: Summary of the parameter values used in Equation (8) to describe the threedifferent RT responses used in model simulations. In all cases, we fix ξ R = 0 . the cell response to RT, we allow consideration of the selective pressure ofRT. For a given dosage and LQ model, differences in the radio-sensitivity ofCSCs and differentiated cells are determined by the ratios α min /α max ∈ (0 , β min /β max ∈ (0 , α min /α max and β min /β max approach the value of unity, RT offers no selective advantageto CSCs as, at leading order, the response is independent of phenotype. Thelatter also depends on the specific dose applied. For example, for high dosesthe quadratic term in Equation (8a) is dominant and the selective pressureis 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 selectivepressure of RT is associated with both α min /α max and β min /β max . For thesereasons, we will consider two different RT protocols: either a single dose of10 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 tohave the least RT selective pressure in both scenarios, this might be higherin R1 or R3 depending on the treatment protocol considered. Plasticity is an essential feature of phenotypic adaptation to changing en-vironmental conditions [14, 37]. It assumes that cells with the same genomecan acquire distinct phenotypes depending on their epigenetic status, which isalso inheritable. Phenotypic variation may be mediated by random (sponta-neous) epigenetic mutations [22], which we assume to be rare. We account forthis effect by including in the structural flux a diffusion term with a constantdiffusion coefficient θ = 5 × − hr − (see Equation 1a). Such random muta-12ions should not favour any specific phenotype, and Darwinian selection (i.e.the fitness function F ) drives phenotypic evolution of the population. Thisaspect has been widely studied in previous work in order to investigate howcells adapt to different environments [24, 22, 25]. At the same time, thereis evidence that phenotypic switching may be mediated by environmentalfactors via Lamarckian selection (or induction) [37]. In this framework, cellsadapt to their environment [14, 56] by following a preferential ( biased ) tra-jectory in phenotypic space. We can, therefore, envisage situations in whicha subpopulation may be prevalent in a population without being the fittestsubpopulation (i.e. the population with the highest proliferation rate). Forexample, recent studies have identified cell de-differentiation and CSC main-tenance as stress responses to harsh environmental conditions [37], includinghypoxia. More specifically, cells respond to hypoxic stress by up-regulatingHypoxia Inducible Factors (HIFs) which, in turn, promote the expression ofstem-related genes [29, 30, 31, 32]. HIF suppression has also been linked tocell differentiation and reduced levels of stemness [57]. We account for suchmicro-environment mediated adaptation by incorporating an advective termin the structural flux. Cells are assumed to evolve along the stemness axiswith a velocity v z = v z ( z, c ), that depends on the oxygen concentration c andcell phenotype z . Under normoxia, cells tend to differentiate, and v z > v z being analogous to a maturation velocity. In our model,however, ageing (i.e. differentiation or loss of clonogenic potential [13]) maybe reversible. For example, under hypoxia (i.e. c ≤ c H ), we assume v z < (a) ξ + = 0 .
05 (b) ξ + = 0 . ξ − = 0 . 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. v z : v z ( z ; c ) = v + z ( z ) H (cid:15) ( c − c H ) − v − z ( z ) H (cid:15) ( c H − c ) , (9a) v + z ( z ) = V + V ∗ + tanh (cid:32) z ω + ξ + (cid:33) tanh (cid:32) (1 − z ) ξ + (cid:33) , (9b) v − z ( z ) = V − V ∗− tanh (cid:18) zξ − (cid:19) tanh (cid:32) (1 − z ) ω − ξ − (cid:33) . (9c)where H (cid:15) is a smooth variant of the Heaviside function approaching the latterin the limit of (cid:15) → H (cid:15) ( x ) = (1 + tanh( (cid:15) − x )) / V ∗± ensure that (max z v ± z ) /V ± = 1 and V ± ( hr − )corresponds to the magnitude of the velocity. Further, by controlling theadvection speed along the stemness axis, V − ± determines the timescales formaturation and de-differentiation. The parameters ξ ± regulate the slopes of v z at the boundaries z = 0 ,
1. As shown in Figure 3a, when ξ ± (cid:28)
1, the ad-vection velocity is steep when z ∼ , ξ ± , the variationis more gradual, with a single maximum (or minimum) near z ∼ . ω ± allow us to tune the symmetry/asymmetryin v z and also to modulate the flux at the boundaries (see Figure 3). Forexample, if ω + = 2, then v (0) = ∂ z v (0) = 0 which means that CSCs will beless likely to differentiate compared to the case ω + = 1. In the absence ofexperimental data with which to specify the parameters in the phenotypicdrift velocity, we consider combinations of the following parameter sets: • V ± ∈ { , , } × − [ hr − ], • ξ ± ∈ { . , . , . } , and • ω ± ∈ { , } .In summary, our phenotype-structured model for the growth and responseto radiotherapy of a solid tumour is defined by Equations (3)-(9). A list ofthe model parameters and estimates of their values can be found in Table A.4in Appendix A. 14 . 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 alongthe stemness axis evolves under normoxia and hypoxia. Our numerical solu-tions are generated using the method of lines, with discretisation performedin the z -direction. In more detail and following [60], we use a finite volumescheme, opting for a Koren limiter to control the advection component of thestructural flux. In this way, we reduce (3) to a system of time-dependent, or-dinary differential equations which can be solved in MATLAB using ode15s ,an adaptive solver for stiff equations. The numerical simulations are validatedin Section 3.3 where we perform a linear stability analysis. The associatedeigenvalue problem is solved numerically using MATLAB’s chebfun package[61]. In well-oxygenated environments, the advection velocity is positive and cellsare driven towards a terminally differentiated phenotype, with z = 1. De-pending on the balance between the advective flux and cell renewal (i.e.,Darwinian selection and Lamarckian induction), the model predicts a vari-ety of long-time behaviours: the system relaxes to its steady state via dampedfluctuations or monotonically. We start by considering symmetric velocityprofiles (see Figure 3a). As summarised in Figures 4 and 5, as the magnitudeof the advection velocity, V + , and its steepness, ξ + , are varied, the systemexhibits different long time behaviours, even though the dynamics at earlytimes are similar for all parameter sets considered (see Figure 4). If simula-tions are initialised with a small population of cells with z ∼ .
5, then thedynamics are initially dominated by proliferation. Over time, as φ increases,competition slows the cell proliferation rate and phenotypic advection be-comes more important. As the cells mature, they accumulate near z = 1,and the rate of natural cell death exceeds the rate of cell proliferation. Fromthis time onwards, the growth curves corresponding to different parametersets start to deviate.For example, in case A.2, the system rapidly relaxes to a non-zero steadystate 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 exponentially15 igure 4: Results from a series of numerical simulations of Equations (3)-(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 Equa-tion (9)). As V + increases and ξ + decreases, the system can be driven to extinction. SeeFigure 5 for the values of the other model parameters. to extinction at a rate dictated by d f . In other parameter regimes, the re-laxation phase is characterised by damped fluctuations. In case A.1, forexample, fluctuations are driven by the interplay between apoptosis, compe-tition and advection. As TDCs are eliminated, the reduction in competitionallows re-growth of highly proliferative cancer cells (i.e., z ∼ . n ( z, t ) shown in Figure 4 for case A.1. Over time, the fluc-tuations decay and the system relaxes to its steady state distribution. InSection 3.3, we present a complementary investigation of this behavior, re-lating the damped oscillations to a complex eigenvalue in the linearisationabout the equilibrium solution.Focusing on the long time behaviour, the symmetric advective profilegives rise to a population with a unimodal equilibrium distribution wherethe 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 thepeak is shifted towards the centre of the domain. These observations aresummarized in Figure 5, where we have further analysed how the propertiesof 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 abifurcation as V + and ξ + vary, with the system transitioning from a trivial toa non-zero steady state (this behaviour will be investigated in Section 3.3).16 igure 5: Series of phase diagrams characterising the steady state distribution predictedby the model as properties of the advection velocity, v z , vary (i.e., for different values ofthe 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 peaksand the dominant phenotype (i.e., the z -locations of the local maxima) for different valuesof the parameters ω + and d f . For parameter sets that give rise to a significant fraction ofCSCs (i.e., % CSCs ≥ φ CSC (0 . , t ∞ ), as defined by Equation (4), is alsoindicated. By contrast, the equilibrium distribution for an asymmetric velocity pro-file (i.e., ω + = 2, as in Figure 3b), has a multimodal distribution, typicallycharacterised by two peaks. In this case, since the CSCs have a lower propen-sity to mature, they accumulate and persist in the population, even under17ormoxia. The second column of Figure 5 shows that the proportion ofCSCs at long time increases as the death rate, d f , the steepness parameter, ξ + , and the maturation velocity, V + , increase, until the CSCs become thedominant subpopulation (see, for example, Case B.3 in Figure 6). Varyingthe death rate, d f , does not significantly affect whether extinction occurs;rather, it determines the location of the maximum peak in the equilibriumdistribution (see, for example, case B.2 in Figure 6). For low death rates,cells are predominantly in a terminally differentiated state. As the death rateincreases, the peak moves to the left, producing an equilibrium distributionin which a higher proportion of rapidly proliferating cells balances the highdeath 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 arecharacterised by secondary regrowth, driven by the accumulation of CSCs.For example, in case B.1, phenotypic diffusion enables the cancer cells tode-differentiate, acquire a stem-like phenotype and, therefore, contribute topopulation growth. Figure 6: Results from a series of numerical simulations of Equations (3)-(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 , determinewhether the model predicts extinction or persistence of CSCs, regardlessof whether they are present initially. When ω + = 2, random mutations (i.e.,diffusion), may dominate the advective force near z = 0, allowing CSCs firstto form, then to proliferate and ultimately to comprise a significant propor-tion of the equilibrium population. CSCs have been observed in normoxicregions; for example, they have been found in perivascular tumour regions,where endothelial cells secrete factors that inhibit CSC maturation [62]. By18ontrast, when ω + = 1 (i.e., for symmetric velocity profiles), all cells matureover time, leading to the eventual extinction of CSCs. This behaviour coulddescribe that of tumours which lack CSCs, or the effect of drugs which inducestem cell differentiation and, thereby, reduce the incidence of resistance toother treatments, such as radiotherapy. We conclude that targetting V + and ξ + may be effective for eliminating CSCs, increasing tumour sensitivity totreatment and, in certain scenarios, driving tumour extinction. Under hypoxia, the advection velocity in our model is negative and cells willbe driven to de-differentiate. In this case, the equilibrium distribution isunimodal, with the dominant phenotype at z = 0. Although varying thedeath rate d f does not effect the equilibrium distribution (compare cases H3and 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 thepopulation.Differences in the system dynamics also arise as the initial conditions n ( z ) vary. The results in Figure 7a indicate little variation in the systemdynamics when the initial conditions from Section 3.1 are used. By contrast,in Figure 7b we observe marked differences when the initial conditions arecentred 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 tumourregrowth. Similarly, increasing the death rate, d f , reduces the number ofcells that can de-differentiate and, subsequently, delays regrowth. Therefore,while d f does not affect the equilibrium distribution, it influences the systemdynamics. These results show how the formation of hypoxic regions canshape the development of a tumour. In particular, the emergence of hypoxiamaintains and enhances the pool of CSCs, preventing population extinction(see, for example, scenario D in Section 3.1). We now validate some of the above numerical results by performing a linearstability analysis which enables us to characterise the equilibrium states. Wedenote by ¯ n = ¯ n ( z ) a steady state for the (untreated) system (3)-(9), with atotal cell density ¯ φ = (cid:82) ¯ n ( z ) dz and let δn represent a small perturbation to19 a)(b) Figure 7: Numerical results under hypoxic condition for four parameter sets, all with V − =4 × − . In (a) we use the standard initial condition defined by Equation (3f) while in (b)the population is centred around z = 1. The other parameter values are as follows: (H1) ξ − = 0.05, ω − = 1 and d f = 0.001; (H2) ξ − = 0 . ω − = 1 and d f = 0 . ξ − = 0 . ω − = 2 and d f = 0 . ξ − = 0 . ω − = 2 and d f = 0 . this solution. Then we can approximate the solution n in a neighbourhoodof ¯ n as: n ( z, t ) = ¯ n + δn ( z, t ) , (cid:107) δn (cid:107) (cid:28) ∀ t > . (10)Substituting this ansatz into (3) and retaining linear terms, we obtain thefollowing equation for δn : ∂δn∂t = M δn, (11a) ∂δn∂z = 0 , z = 0 , , (11b) δn ( z, ≡ , (11c)where M is the following integro-differential operator M δn ≡ ∂∂z (cid:18) θ ∂δn∂z − v z δn (cid:19) + (cid:2) p (cid:0) − ¯ φ (cid:1) − f (cid:3) δn − p ¯ n (cid:90) δndz. (12)20he solution ¯ n is spectrally stable if the spectrum of the operator, σ ( M ),does not contain eigenvalues with positive real part, i.e., σ ( M ) (cid:92) { λ ∈ C : (cid:60) ( λ ) > } = ∅ . (13)Moreover, the dynamics of the system will be dominated by the fastest grow-ing mode (i.e., the eigenfunction corresponding to the eigenvalue with thelargest real part, λ ).In Appendix B we transform the above eigen-problem so that it does notinclude any first order derivatives. For a non-zero steady state, we retain anon-local term in the eigenvalue problem and this can give rise to a spectrumwith a pair of complex eigenvalues. Recalling case A.1 from Section 3.1 (seeFigure 4), the numerically estimated value of λ is indeed complex ( λ = − . × − ± i . × − , where i = − n ≡
0, which is al-ways a fixed point for the system, the non-local term vanishes and we obtainthe standard form analysed by Sturm-Liouville theory. Using well knownresults, we can identify sufficient conditions for the stability/instability ofthe 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 pa-rameter 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 notethat the results relate only to the behaviour of the fitness function and ad-vection velocity near the boundary z = 0, suggesting that the most relevantparameters are p max , V − , θ and ξ − . By contrast, under normoxia, and forthe range of parameter considered here, the system undergoes a bifurcation.For sufficiently small V + , the trivial steady state is unstable; for sufficientlylarge V + and for large values of the death rate, d f , the trivial steady stateis stable, (see, for example, case C2 in Section 3.1). To investigate otherparameter regimes that we can not tackle analytically, we rely on numericalestimation of the largest eigenvalue, λ . As shown in Figure 8, it is possibleto identify the boundary of the region of stability in ( ξ + , V + ) space. Thisdiagram does not change significantly as the death rate varies in the rangefrom d f = 0 .
001 to d f = 0 .
015 (results not shown). However, the results arehighly sensitive to the value of ω + . Comparing Figures 8a and 8b, we seethat setting ω + = 2 favours the formation of a non-trivial equilibrium distri-bution, with the curve shifting to the far right of the parameter space (i.e.,21 a) ω + = 1 (b) ω + = 2 Figure 8: Series of phase diagrams partitioning the ( V + , ξ + ) parameter space into regionswhere the trivial steady state is linearly stable (green regions) and unstable (white regions).The diagrams are obtained for d f = 0 . ω + has a significantimpact on the size of the region of ( V + , ξ + ) parameter space in which the non-trivial steadystate is stable (compare (a) and (b)). small values of ξ + and large values of V + ). In the latter case, this impliesthat even higher velocities V + are needed to stabilise the tumour eliminationsolution. This is consistent with the numerical results in Section 3.1, wheresetting ω + = 2 (see scenario B in Section 3.1) favoured the accumulation ofCSCs which acted a reservoir for tumour cells.
4. Population Dynamics in the Presence of Treatment
In the previous section, we found that the system possesses a stable steadystate to which the dynamics converge for the range of parameter values con-sidered. Therefore, we anticipate that, while treatment can perturb the sys-tem from its equilibrium, it will eventually relax to its stable steady stateonce treatment ends. Thus we expect extinction to occur for parameter val-ues lying in the stability region of the trivial steady state (see Figure 8).From this point of view, we are interested in understanding how differentenvironmental conditions (i.e. normoxia and hypoxia), different treatmentprotocols 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 differentadvection velocities used in our earlier analysis (see Table 3). Starting fromthe initial condition (3f), cells follow different pre-treatment protocols as22pecified in Table 3. Without loss of generality, we shift time so that t =0 corresponds to 24 hours before treatment begins. While attention willfocus on tumour responses in constant environmental conditions, we alsoconsider briefly treatment responses in changing environments. For eachscenario, we simulate the response to treatment for the range of values ofthe radiation parameters listed in Table 2. We denote by n ( S ,R ( z, t ) thesolutions corresponding to scenario S R The simulation results presented in Figure 9 illustrate the different regrowthdynamics that can arise when well-oxygenated tumour cells are exposed toa single dose of RT. We identify three distinct behaviours: instantaneous re-growth (S1), decay and extinction (S2) and initial remission with subsequentregrowth (S3). While the cell survival fraction immediately post-treatmentdepends on the parameter values used in the LQ-model (see Equation (8)),the qualitative population regrowth dynamics post-treatment do not dependon these values.In more detail, for scenario S1, the cell density increases rapidly aftertreatment, driving the system towards its (asymptotic) equilibrium. By con-trast, for scenarios S2 and S3, the growth curves initially decrease at similarrates until about 40 days after treatment. Thereafter, for scenario S S
2, the tumour continues to shrink, until it is eventually eliminated.The origin of such differences can be understood from the time evolutionof 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 dy-namics. The reduction in the cell density φ ( t ) post-radiotherapy reducesintra-population competition and allows the cells to resume proliferation.Depending on the magnitude of the advection velocity (which is positive),the cells either regrow ( S
3) or they are driven to a terminally differentiatedstate and, thereafter, become extinct ( S S
3, the presence ofradioresistant CSCs post treatment and a small positive velocity at z = 0together drive regrowth. As the CSCs start to mature, there is a continuoussource of highly proliferative cells which, in turn, drive rapid regrowth of thetumour. As the total cell number increases, intra-population competitionslows cell proliferation until eventually advection becomes dominant, driving23cenario Protocol Parameters Subsection S V + [10 − ] , ξ + , ω + , d f )=(4 , . , , . S V + [10 − ] , ξ + , ω + , d f )=(8 , . , , . S V + [10 − ] , ξ + , ω + , d f )=(8 , . , , . S V − [10 − ] , ξ − , ω − , d f )=(2 , . , , . S V ± [10 − ] , ξ + , ξ − , ω + , ω − , d f )=(8 , . , . , , , . S V ± [10 − ] , ξ + , ξ − , ω + , ω − , d f )=(8 , . , . , , , . Table 3: Parameter sets used to generate the numerical simulations in Section 2.1, togetherwith the corresponding environmental conditions pre- and post-treatment (blue: normoxia,red: hypoxia). Simulations are initialised using equation 3f at different times t = − t s asindicated in the second column. Radiotherapy is administered at time t = 24 hours. Theparameter values have been chosen to illustrate the range of qualitative behaviours thatthe model exhibits. igure 9: Different treatment outcomes under normoxia. For each scenario S1, S2 andS3 (see Table 3) we consider the dynamics of the total cell number, φ ( t ), and comparethe responses for the radio-sensitivity parameter sets R1, R2 and R3 (see Table 2) to thecontrol, untreated case. For each scenario we also present plots of the phenotypic celldistribution, n ( z, t ), at different times for radiotherapy protocol R1. The vertical lineindicates the time of irradiation, while a line is also shown that follows the evolution ofthe control (i.e., in the absence of treatment). the cells to de-differentiate. By contrast, for scenario S
2, advection domi-nates proliferation along the entire phenotypic axis. Additionally, CSCs areabsent so that all cells are rapidly terminally differentiated and, thereafter,undergo cell death.Comparison of scenarios S2 and S3 reveals how different phenotypic com-positions can generate treatment responses which are initially qualitativelysimilar, but differ markedly at long times. This finding is reinforced in Figure10 where we plot the mean phenotypes, µ = µ ( t ), as defined by Equation (5).For scenarios S2 and S3, the dynamics of the mean phenotype are indistin-guishable at short times and do not start to diverge until approximately 20days after treatment.More generally, the results presented in Figure 10 reveal three charac-25 igure 10: Series of plots showing the evolution of the phenotypic mean, µ ( t ), for scenariosS1, 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 radio-therapy. The dynamics of µ may be the same as those prior to treatment,with negligible deviation from the control (see scenario S2). A discontinuityin µ may be induced by radiotherapy (see scenario S µ towards less mature phenotypes. The size of the discontinuitydepends on the relative sensitivity of CSCs and TDCs to RT, or, using theterminology introduced in Section 2.1, the selective power of RT. Since weare considering high radiation dosages, the discontinuity is determined bythe ratio β min /β max . In order for the selective pressure of treatment to beapparent, CSCs must comprise a significant fraction of the population priorto treatment. This explains why, for scenario S3, there is an initial transientperiod during which, as for scenario S2, there is no discernible deviation fromthe 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 radio-sensitivity. It is natural to expect cells that have permanently exited thecell-cycle will be less radio-sensitive than cycling cells, as the DNA damageresponse may already be active in such cells [36]. The functional forms for α and β defined by Equations (8b)-(8c) assume that radio-sensitivity increasesmonotonically with cell phenotype, z . In order to investigate situations inwhich TDCs have lower radio-sensitivity than proliferating cancer cells, we26 igure 11: Series of numerical results showing how the growth dynamics and the pheno-typic mean evolves following exposure to a single dose of radiotherapy when cell radio-sensitivity is a non-monotonic function of cell phenotype. The simulations are analogousto those presented in Figure 9 and 10, except that Equations (14) are used in place ofEquations (8b)-(8c). now the following, non-monotonic functional forms: α ( z ) = α min + ( α max − α min ) tanh (cid:18) zξ R (cid:19) H . (1 − z ) , (14a) β ( z ) = β min + ( β max − β min ) tanh (cid:18) zξ R (cid:19) H . (1 − z ) , (14b)where H (cid:15) is defined in § (cid:15) = 0 .
075 (all other pa-rameters are as defined in § µ ( t ) (compare Figures 11 and 10). The differences are most pronouncedfor scenarios S S z = 1, are dominant inthe population prior to treatment. The qualitative growth dynamics (i.e., φ ( t )) is similar for both cases. Further investigation of these differences isbeyond 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 practical27 igure 12: Simulation results for fractionated radiotherapy protocols, showing how thetotal cell number φ ( t ) and the phenotypic mean µ ( t ) evolve for scenarios S1 and S3 (seeFigure 9 for details). In all plots, the light purple shaded area indicates the variabilityin responses when a single dose of 10 Gy is administered and is included for comparisonwith the fractionated treatments (see Figure 10). The yellow shaded area indicates theduration of the treatment for the fractionated case. for treating patients, due to adverse side effects [63]. Therefore, we nowconsider tumour responses to fractionated RT protocols. The trends forfractionated RT are similar to those for single doses for all scenarios in Table3. Typically, the proportion of cells that survive fractionated therapy islarger than for the single-dose case, by a factor of about 100. Consequently,for scenarios S S
3, the time to return to the equilibrium populationdistributions is reduced. For S2, while treatment causes a monotonic decreasein the cell density φ , since more cells survive fractionated RT, it takes longerfor the cell population to become extinct. For scenarios S S
3, we recallthat for high doses of RT, the phenotypic mean was markedly affected bythe specific LQ model parameters considered; this is not the case when lower28oses are applied (see Figure 12).
Figure 13: Phenotypic distribution n ( S ,R ( z, t ) for the control (light blue), the colonyexposed to a single dose (dark blue) and the one treated with fractionated dose 2 Gy × t pt is 24 hr and 120 hr for the 10 Gy andfractionated protocol respectively. On the other hand, the remaining panels are measuredrelative 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 singledose of radiotherapy can be attributed to the temporary advantage CSCshave post treatment. When using a fractionated protocol, intra-populationcompetition is maintained at the cost of fewer cells being killed. This is ap-parent when we compare the phenotypic distribution at different times forthe two treatment protocols (see Figure 13). When 10 Gy is administeredin one dose (first panel, dark blue region), the peak of the distribution is at z = 0. On the other hand, after 5 doses of 2 Gy per day (first panel, greenregion), the proportions of differentiated and cancer stem cells are approx-imately equal. Given that the former proliferate faster than the latter, thedifferentiated cells quickly become the dominant phenotype. Consequently,one month after treatment ends (third panel in Figure 13), the proportion ofCSCs in the population is the same for both protocols. We conclude furtherthat the single dose protocol outperforms the fractionated protocol when wecompare the total number of cells (the blue curve is below the green one forall values of z ). Cell populations that are continuously exposed to hypoxia, exhibit instanta-neous re-growth following RT, as shown in Figure 14. Compared with the29reatment outcome under normoxia, a higher percentage of cells survive ra-diation, because there is a larger proportion of radio-resistant cells in thepopulation under hypoxia. Even though a smaller fraction of cells are killed,re-growth is also usually slower under hypoxia than under normoxia. Wenote also that, following exposure to the single and fractionated protocols,the phenotypic mean µ ( t ) shifts toward z = 0 under hypoxia, favouring CSCsas the dominant phenotype (see Figure 14). The drift in µ is less pronouncedfor the fractionated case, suggesting the latter protocol is less favourable forthe immediate accumulation of resistant subpopulation of CSCs than thesingle dose. Figure 14: Comparison of the tumour cell responses to single and fractionated radiotherapyprotocols under hypoxia for scenario S4 (See Table 3). Simulation results showing thetime evolution of the cell density, φ ( t ), and phenotypic mean, µ ( t ), are presented. Forcomparison, the light purple shaded areas in the fractionated plots indicate the variabilityin the response when a single dose of 10 Gy is administered. The yellow shaded areasindicate the duration of treatment for the fractionated case. Thus far we have assumed that the oxygen concentration remains constantthroughout treatment. While this may accurately describe RT responses forcells 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 vascu-lature and, hence, tissue re-oxygenation. On the one hand, high doses ofradiotherapy may damage the vasculature [67], and decrease nutrient avail-ability post radiotherapy. On the contrary, moderate RT may transientlyincrease tissue oxygenation by normalising the tumour vasculature (vesselnormalisation is a phenomenon that has been observed when tumours areexposed to vascular-targetting agents which destroy some of the blood ves-sels 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 survivingcells may increase. Equally, hypoxic regions may form at later times as thetumour regrows. From this point of view, radiotherapy may impact boththe phenotypic distribution of the cell population (and, thereby, its radiosensitivity), and oxygen levels post-treatment. We can use our mathematicalmodel to investigate these scenarios, by assuming that oxygen levels changepost radiotherapy.Based on the results presented in Sections 4.1 and 4.2, we anticipate thatreoxygenation of a hypoxic tumour will be beneficial in certain cases, drivingCSC maturation, and even leading to tumour eradication. The results pre-sented in Figure 15 show that the long-term tumour regression is precededby an initial phase of regrowth during which CSCs that survive treatmentde-differentiate and proliferate. Such a treatment might initially be con-sidered unsuccessful, although the stability of the trivial steady state uponre-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 hy-poxia [64]. Figure 15b shows that such environmental changes may negativelyimpact the outcome. The formation of an hypoxic region favours the devel-opment and maintenance of radioresistant CSCs, reducing the treatment effi-31 a) S5(b) (c) S6
Figure 15: Growth curves for changing environmental conditions: (a) re-oxygenation and(b) post-radiation hypoxia for the parameter values S5 and S6 in 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 R T R the time at whichre-oxygenation occur (indicated by the arrows in the plot). If T R is sufficiently smallthan re-oxygenation does not drive re-growth of the cells population. If we waited for asufficiently 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 tobe replaced by new vessels which form via angiogenesis and re-oxygenate thedamaged regions. As shown in Figure 15c, depending on the time-scale re-quired for vessel regrowth (indicated by T R ), different behaviours may arise.If the duration of RT-induced periods of hypoxia is sufficiently short, thenthe size of the cell population remains low. By contrast, if there is sufficienttime for cells to de-differentiate (see T R = 1000), then re-oxygenation leadsto a rapid increase in cell number, although eventually the cells die out.These results highlight the complex interplay between tumour growth and32reatment response in vivo and the importance of environmental factors indetermining the eventual outcome of radiotherapy treatment.
5. Conclusion and Future Challenges
We have developed a structured model to investigate how clonogenic hetero-geneity affects the growth and treatment response of a population of tumourcells. Cell heterogeneity is incorporated via an independent and continuousstructural variable which represents stemness . As proposed by [13, 19], weview stemness as a plastic trait, with cells becoming more, or less, stem-like depending on their environmental conditions. Our mathematical modelaccounts for cell proliferation and apoptosis, inter-cell competition, and phe-notypic movement along the stemness axis, via diffusion and advection.Studies of the population dynamics in the absence of treatment revealedthat, under normoxia, a variety of qualitative behaviours may arise depend-ing on the functional forms used to represent the structural flux and fitnesslandscape. When advection dominates movement along the stemness axis,its magnitude, relative to the rates of proliferation and cell death, deter-mines whether the population is driven to extinction. Multimodal distribu-tions, which allow for the formation and maintenance of CSCs pools, areobserved for asymmetric velocity profiles. Under hypoxia, the populationdistribution is unimodal and skewed toward stem-like phenotypes, with littleintra-population variability. The resulting cell distribution is highly resis-tant to radiotherapy, the tumour will typically regrow following treatment.By contrast, under normoxia (or re-oxygenated hypoxia), and for suitableparameter values, the tumour may become extinct following radiotherapy.There are many ways in which the work presented in this paper couldbe extended. A first, natural extension would be to incorporate structuraland spatial heterogeneity (i.e., both phenotypic and spatial dimensions) [21].This would enable us to consider in vivo situations, where spatial gradientsin oxygen levels emerge naturally, due to oxygen consumption by the cells asit diffuses from blood vessels. As outlined in Appendix A, in such a modeloxygen consumption rates may vary with cell phenotype, and spatial fluxesmay account for random movement of the cells. Preliminary results for sucha model are presented in Figure 16. We consider a 1D Cartesian geometryand 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 tumourinitially comprises a spatially homogeneous distribution of terminally differ-33ntiated cells (see Equation (3f)), then the oxygen rapidly relaxes to a steadystate and a hypoxic region forms at distance from x = 0. In contrast to thewell-mixed model, cells are now able to move, by random motion, betweennormoxic and hypoxic regions. While terminally differentiated cancer cellsare dominant in the well-oxygenated region, a small fraction persists in thehypoxic region (in particular, near the boundary of the hypoxic region, or-ange line in the plots in Figure 16). This is due to the influx of cells from thewell-oxygenated portion of the domain. Similarly, CSCs are dominant in thehypoxic region, but a small fraction of hypoxic CSCs migrate towards x = 0,where re-oxygenation induces their maturation, creating a differentiated andhighly proliferative cell phenotype, alongside terminally differentiated cancercells. These results illustrate how the interplay between space, resources andphenotypic adaptation may give rise to complex behaviours; their investiga-tion is the focus of ongoing work. Figure 16: Series of plots showing how, in the absence of treatment, the cancer cell popula-tion n ( x, z, t ) and the oxygen concentration c ( x, t ) change over time t when we account forspatial 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 upperplots and with a vertical orange line in the lower plots. We fix V ± = 4 × − , ξ ± = 0 . ω + = 1, ω − = 2 and d f = 0 . A significant challenge of the modelling approach presented in this paperis the determination of model parameters and functional forms. In the longerterm, techniques such as single-cell RNA sequencing [17, 18] will make it bepossible to quantify specific aspects of our model, such as the dependence ofthe proliferation and apoptosis rates on cell stemness and the dependence onthe tumour micro-environment of the (phenotypic) advection velocity asso-ciated with cell maturation and de-differentiation. In spite of their currentlimitations, we believe that studies of such models can increase understanding34f the ways in which specific physical processes may influence the phenotypicdistribution of cell populations in different environments. At the same time,we acknowledge that it remains a matter of debate as to whether asymmetriccell distributions are driven by micro-environmental signals (as in the modelpresented here), asymmetric division, or a combination of the two [70]. Byusing a non-local proliferation kernel to account for asymmetric division, wecould investigate these alternative hypotheses and identify conditions underwhich they lead to different outcomes.A important feature of our model is the way in which the response to ra-diotherapy (RT) varies with cell stemness (i.e., z ). Our analysis shows howthe functional forms used to describe the advection velocity and fitness func-tions can affect the system dynamics post-RT. While unimodal phenotypicdistributions lead to monotonic growth curves post-treatment, more complexbehaviour is observed when heterogeneous populations, with a pool of CSCs,are considered. For example, under normoxia, the presence of radio-resistantCSCs can drive recurrence, despite an initial phase of tumour regression. Asthe CSCs mature into highly-proliferating cancer cells, rapid re-growth is ac-companied by re-sensitisation of the population to RT. Under hypoxia, CSCsmaintain their stemness, leading to a slowly growing, radio-resistant cellpopulation. More complex outcomes arise when we consider the effect thattreatment might have on the environment. As noted in 4.3, changes in thevasculature induced by radiotherapy can result in either post-treatment re-oxygenation or hypoxia. While re-oxygenation increases the radio-sensitivityof the population, hypoxia increases their radio-resistance. In practice, suchenvironmental changes are likely to be transient. Even in an untreated tu-mour, fluctuations in oxygen levels can occur. Consider, for example, cellsin a neighborhood of immature blood vessels. As the cells proliferate, theyexert mechanical pressure on the vessels, causing them to collapse and localoxygen levels to fall. Under hypoxia, the tumour cells stimulate the growthof new blood vessels from pre-existing ones, via angiogenesis. In this way, tu-mour regions may cycle between periods of hypoxia and normoxia. It wouldbe of interest to extend the model to account explicitly for the tumour vas-culature 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 viamore complex, multi-phase [73] or multi-scale approaches [74, 75, 76, 77].This would enable us to better capture the different time-scales on whichthe oxygen dynamics and cell adaptation velocity change. As shown in Fig-ure 17, variations in oxygen levels emerge naturally within spatially-resolved35odels. Here, cell killing leads to tissue re-oxygenation which, in turn, dis-rupts the CSC niche. Depending on the time scale over which the cells adaptto their new environmental conditions, this may increase the overall radio-sensitivity. Understanding and accounting for such phenomena is particularlyrelevant for predicting responses to RT and comparing alternative treatmentprotocols. Figure 17: Evolution of the population n ( x, z, t ) in the spatial and phenotypic dimensionsfollowing a cycle of fractionated radiotherapy (5 × R In the extinction scenario, or post administration of high RT doses, thenumber of cells in the population can become low and our continuum modelmay cease to be valid. In such conditions, stochastic effects which are ne-glected herein may become important. As in [78, 79, 80], stochastic and meanfield approaches may be combined with hybrid discrete-continuum techniquesto account for small population effects and to study their impact on the prob-ability of tumour extinction.In this paper, we considered only single dose and fractionated treatmentprotocols. In future work, we could investigate alternative strategies, suchas adaptive therapeutic protocols [81] and/or multi-drug treatments, whichhave been proposed as an effective way to overcome radio-resistance. Fromthis point of view, considerable efforts have been invested in designing treat-ments that exploit features of CSCs, such as their metabolic plasticity [10].Motivated by recent metabolically-structured models [24, 21, 25], a naturalextension of our model would be to include a “metabolic dimension” in or-der to investigate the interplay between stemness, metabolic switching andresistance. A biologically informed model that incorporates metabolic andphenotypic effects, together with the tumour micro-environment and vascular36emodelling lies at the heart of a mathematical program that would enablesystematic comparison with in vivo observations. The framework and resultsoutlined in this work represent a first step towards achieving this long-termgoal.
Appendix A. Spatial Model
We outline here the set up for the 1D simulations presented in Section 5. Asa full description of the spatial model goes beyond the scope of the presentwork, we focus on the main changes to (3)-(9). We now view the oxygenconcentration c as a dependent variable, rather than a prescribed function.We suppose that oxygen is supplied to the region by blood vessels on thedomain boundary ∂ Ω (see Figure A.18). Oxygen diffuses from the boundaryinto the tissue where it is consumed by the tumour cells at rates which dependon their phenotype and the local oxygen concentration. The evolution of thedimensionless cell density, n = n ( x , z, t ), is driven by a phenotypic flux ofthe same form as in Equation (3) but a spatial flux is included to accountfor random motion in the spatial dimension. Figure A.18: Schematic representation of the phenotypic and spatial model.
As shown in Figure A.18, we consider a fixed tissue slice where the oxygensupply (i.e. vasculature) is confined to one of the tissue boundaries. Given theassumed symmetry of the problem, we can consider a 1D Cartesian geometrywith x ∈ [0 , L ]. The spatial model is defined by the following system of37oupled PDEs: ∂n∂t = D N ∂ n∂x (cid:124) (cid:123)(cid:122) (cid:125) spatial flux + ∂∂z (cid:18) θ ∂n∂z − nv z ( z, c ) (cid:19) + F ( z, c, φ, t ) n, (A.1a) ∂c∂t = D C ∂ n∂x − Γ( t, x, c ) , (A.1b) θ ∂n∂z − nv z = 0 , z ∈ { , } , x ∈ [0 , L ] , t > , (A.1c) ∂n∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 = ∂n∂x (cid:12)(cid:12)(cid:12)(cid:12) x = L = 0 , z ∈ (0 , , t > , (A.1d) ∂c∂x (cid:12)(cid:12)(cid:12)(cid:12) x = L = 0 , c (0 , t ) = c ∞ , t > , (A.1e) n ( x, z,
0) = n ( x, z ) x ∈ [0 , L ] , z ∈ (0 , , (A.1f) c ( x,
0) = c ( x ) x ∈ [0 , L ] , (A.1g) φ ( x, t ) = (cid:90) n ( x, z, t ) dz, (A.1h)Γ( t, x, c ) = (cid:90) γ ( z, c ) n ( x, z, t ) dz, (A.1i) F ( z, c, φ, t ) = p ( z, c ) (1 − φ ) − f ( z ) − gH ( c N − c ) (cid:124) (cid:123)(cid:122) (cid:125) necrosis − N (cid:88) i log (cid:18) SF ( z, c ) (cid:19) δ ( t − t i ) . (A.1j)In Equation (A.1), D N and D C are the assumed constant spatial diffusioncoefficient for the cells and oxygen, respectively, while γ denotes the rate atwhich cells of phenotype z consume oxygen and Γ the net rate of oxygenconsumption at position x and time t . The advection velocity v z is as de-fined by Eq. (9), while the fitness function F is analogous to that definedin Section 2.1, with an additional term to account for necrosis. The latteris assumed to occur at a constant rate g ≥
0, independent of cell pheno-type, when the oxygen concentration falls below a threshold value, c N ≥ SF given in § oxygen-enhancement ratio (OER) [34, 33].According to the oxygen fixation hypothesis [82], part of the biological dam-age induced by radiation is indirect, being mediated by the presence of free38adicals. Thus, when oxygen is limited, radio-sensitivity is accordingly re-duced. Based on experiments, the range of oxygen concentrations at whichthis effect is relevant corresponds to more severe levels of hypoxia (where c ∼ .
5% or lower). We do not consider such situations for the well-mixedmodel, where we consider mild hypoxia. However, accounting for the OERwill be important for the spatially extended model. Recall from Section 3.2that hypoxia is a favourable niche for CSCs. Therefore the OER will endowthem with additional protection from radiation. Denoting by c RH the oxygenthreshold at which the OER becomes active, we use the following functionalform for the survival fraction when simulating the spatially-extended model: SF ( z, c ) = exp [ − α ( z ) d − β ( z ) d ] c > c RH exp (cid:20) − α ( z ) OER d − β ( z ) OER d (cid:21) c < c RH . (A.2)In Equation (A.2), α and β are defined by (8). We note that in the maintext, we consider c = 1 (normoxia) and c = 0 . γ which depends on their phenotype, z . As mentioned previously, stem cellsare known to have a glycolytic metabolism and, thus, we assume that theyconsume less oxygen than cancer cells. Consequently, we consider γ to be amonotonically increasing function of the phenotypic variable z which asymp-totes to its maximum value for z > . γ ( z, c ) = H ( c − c N ) (cid:104) γ max − γ max e − k γ z (cid:105) . (A.3)In Equation (A.3), H = H ( x ) is the Heaviside function (i.e. H ( x ) = 1 if x > H ( x ) = 0 if x ≤ Appendix A.1. Parameters
The model contains a large number of parameters, most of which willvary in value between tumours and patients. The main focus of this work39 arameter Value Units Reference LabelPhenotypic Diffusion θ × − hr − -Advectionvelocity v z Eq (9) V ± { , , } × − hr − - ξ ± { . , . , . } - - ω ± { , } - -Fitness F Eq (6)-(7) p max hr − [84] K H, .
05 - - g p max hr − [84] K H, . g d f { . , . } hr − - k f
10 - -Φ max cell/cm [85]SurvivalFraction SF Eq (8)/Eq (A.2) α min,max Table 2 Gy [12] β min,max Table 2 Gy − [12] ξ R n φ hr − - σ D N . × − mm hr − SDomain Size L 0.45 mm - SOxygen Diffusion D c . × − mm hr − - SConsumption γ Eq (A.3) γ max . × − g(cell hr) − [86] S k γ
10 - - SOxygenthresholds c ∞ c H c N is to study the role played by phenotypic advection (as it interacts with cellproliferation and apoptosis, as well as competition mechanisms). On thisbasis, we decided to perform a parameter sweep for parameters associatedwith the advection velocity, while holding all other model parameters fixed40t values previously reported in the literature, where such values exist. Themain challenge is to identify the phenotypically dependent parameters, suchas the growth rate in Equation (6b). As most data reported in the literaturerefer to processes, such as cell proliferation, at the population/cell colonyand do not account for phenotypic variation, it was difficult to estimateparameters that characterise the phenotypic variation in these processes.We based our estimates of the proliferation rate on the doubling timesreported by [84] for two breast cancer cell lines, MCF-7 and BT-549. Theformer belong to the class of laminal -like cells which are characterised bylow stemness levels [89] and high proliferation rates (doubling time 1 . .
016 hr − ). On the other hand, BT-549 belong to the classof triple-negative cells whose population is dominated by highly aggressivebut slowly proliferating stem-like cells [89] (doubling time 3 . .
008 hr − c ∞ ) to be at a pressure of 100 mmHg [87]. Given that atmospheric pressurecorresponds to 760 mmHg with 21% O , the oxygen tension correspondingto c ∞ is about 8% O . The hypoxic and necrotic thresholds ( c H and c N )are then equivalent to oxygen pressures of 2 . O and 0 . O in line with[90, 88]. These values can be converted into oxygen concentrations by use ofHenry’s law [87], see Table A.4. Appendix B. Linear Stability Analysis.
As mentioned in Section 3.3, in order to compute the largest eigenvalue λ numerically we rely on the Chebfun package for MATLAB [61]. In orderto solve the eigenvalue problem we first make the following substitution inEquation (11c): δn = y ( z ) exp (cid:20) θ (cid:90) z v z ( s ) ds (cid:21) . (B.1)It is straightforward to show that the function y satisfies the following eigen-value problem: θ d ydz + q ( z ; c, ¯ φ ) y − p ( z ; c )¯ n (cid:90) y ( s ) k ( s, z ) ds = λy (B.2)41here q ( z ; c, ¯ φ ) = p ( z ; c )(1 − ¯ φ ) − f ( z ) − dv z dz − v z θ , (B.3)and k ( s, z ) = exp (cid:20) θ (cid:90) zs v z ( p ) dp (cid:21) , (B.4) dydz = 0 at z = 0 , . (B.5)Note that the integral in Equation (B.2)is of the form of a Fredholmintegral which is built in the Chebfun package [61]. The above differentialequation for ¯ n = 0 corresponds to the standard form of a Schr¨odinger-type, Sturm-Liouville eigenvalue problem, where the
Hermiticity of the differentialoperator implies the existence of purely real eigenvalues. In the case of thenull steady state, the eigenvalue problem simplifies to: θ d ydz + ˜ q ( z ; c ) y = λy (B.6a) dydz = 0 at z = 0 , . (B.6b)where ˜ q ( z ; c ) = q ( z ; c,
0) as defined in (B.3). Therefore, applying the
SturmOscillation Theorem [91] to (B.6) we deduce that σ ( M ) has infinitely manysimple and real eigenvalues which can be enumerated in strictly decreasingorder: λ > λ > . . . , lim n →∞ λ n = −∞ . (B.7)We conclude that the trivial steady state is either a stable node (if λ < λ > λ , analytical approximations andbounds can be obtained via the so-called Rayleigh quotient R ( y ). If wemultiply Equation (B.6a) by y and integrate by parts, then we obtain: R ( y ) = 1 (cid:107) y (cid:107) L (cid:90) (cid:26) θy d ydz + ˜ q ( z ; c ) y (cid:27) dz, (B.8a)where y also satisfies the Neumann boundary conditions B.6b. We deducethat the following therefore holds: λ = sup y ∈ E, y (cid:54) =0 R ( y ) (B.8b)where E is the set of twice differentiable functions that satisfy condition (B.6b).42 emma 1. If the function ˜ q is such that max z ∈ (0 , ˜ q < then the null steadystate is stable.Proof. Consider the numerator of the quotient defining R ( y ): (cid:90) (cid:26) θy d ydz + ˜ q ( z ; c ) y (cid:27) dz = − (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:18) (cid:20) y dydz (cid:21) + (cid:90) ˜ q ( z ; c ) y − θ (cid:18) dydz (cid:19) dz ≤ (cid:90) ˜ q ( z ; c ) y dz. (B.9a)We deduce that R ( y ) ≤ (cid:90) ˜ q ( z ; c ) y (cid:107) y (cid:107) dz = R up ( y ) . (B.9b)It is therefore apparent that if the function q is negative throughout thedomain, then R up is negative for any choice of y ∈ E . In such a case, wehave that: λ = sup y ∈ E, y (cid:54) =0 R ( y ) < sup y ∈ E, y (cid:54) =0 R up ( y ) < . (B.9c)We now show that under normoxia, q < Lemma 2.
If the model proliferation rate, apoptosis rate and phenotypicadvection velocity and diffusion coefficient are chosen such that: (cid:90) (cid:26) p ( z, c ) − f ( z ) − v z θ (cid:27) dz > , (B.10) then the trivial steady state is unstable.Proof. Consider y ≡
1, then y ∈ V and λ = R ( y ) where: R ( y ) = (cid:90) (cid:26) p ( z, c ) − f ( z ) − v z θ (cid:27) dz > . Consequently, sup y ∈ V R ( y ) ≥ R ( y ) >
0, and our steady state is unstable.43 emark.
Note that for (B.10) to hold we require (cid:82) ( p − f ) dz > so thatcell proliferation dominates apoptosis. Based on the functional form definedin Section 2.1, we have that: I ( c ; d f ) = (cid:90) ( p − f ) dz = √ g p ( c ) (cid:20) Z (cid:18) z − . √ g (cid:19) + √ g p ( c ) Z (cid:18) z √ g (cid:19) + d f k f e − k f z (cid:21) z =1 z =0 ∼ √ g p ( c )2 + √ g p ( c ) − d f k f (B.11) where Z denotes the cumulative distribution function for the normal distri-bution. We note that I (1; d f ) > while I (0 . d f ) < for all values ofthe parameters listed in Table 1. We conclude that under normoxia thereis a threshold V + ( ξ + , ω ) such that the system is unstable for all choices of V + < V + ( ξ + , ω ) : V + = (cid:115) I (1; d f ) θI v ( ξ + , ω + ) , (B.12a) where I v ( ξ + , ω + ) = (cid:90) (cid:32) V ∗ + tanh (cid:32) z ω + ξ + (cid:33) tanh (cid:32) (1 − z ) ξ + (cid:33)(cid:33) dz. (B.12b) We note also that higher values of θ favour instability of the trivial solutionas 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 sufficientconditions holds, we rely on numerical estimates of the eigenvalue λ . Asshown in Figure B.19, and as expected based on the above findings, whenthe magnitude of the velocity V + is small, λ > ξ and the trivialsolution is unstable. By contrast, as the magnitude of the advection velocityincreases, its steepness, ξ , determines the stability of the trivial solution.Using this estimate, we can identify the region of stability of the trivial steadystate (see Figure 8 in Section 3.3). We remark that the boundary betweenthe regions of stability is non-smooth. This is because λ = λ ( ξ ) plateaus as ξ (cid:28) λ ( ξ ),we observe that the sharp change in the profile of λ as ξ decreases occurs44 a) (b) Figure B.19: Linear stability analysis of the trivial solution: plot of the two largesteigenvalues λ ( ξ ) and λ ( ξ ) for (a) V + = 6 × − , ω + = 1 and d f = 0 . V + = 8 × − , ω + = 1, d f = 0 . λ > ξ . In (b), λ changessign as ξ increases and we can identify a critical value of ξ at which the trivial solutionloses stability, favouring the emergence of a nontrivial, phenotypic cell distribution. where | λ − λ | attains its minimum value. It is possible to show that thetwo 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 hypoxicenvironment. We confirm the numerical simulations from § Lemma 3.
Under hypoxia (i.e. when c = 0 . ), and for the parameter valueslisted in Table A.4, the trivial steady state is always unstable.Proof. Let us consider as a trial function: y = 1( πκ ) / exp (cid:18) − z κ (cid:19) + Az , (B.13)where a small parabolic correction is added to the standard Gaussian, theconstant A being chosen to ensure that the boundary condition at z = 1 issatisfied: A = e − κ π / κ / ; (B.14)the derivative y (cid:48) at z = 0 vanishes, by construction. We now want to showthat the Rayleigh quotient is positive for such a choice of the test function y which implies that the trivial steady state is unstable.45iven that the denominator of R ( y ) is always positive, its sign will bedetermined by the numerator R n ( y ) that is: R n ( y ) = (cid:90) (cid:18) p − f − v z θ (cid:19) y dz − (cid:0)(cid:0)(cid:0)(cid:0)(cid:18) v z y (cid:12)(cid:12)(cid:12)(cid:12) + (cid:90) v z y dydz − θ (cid:18) dydz (cid:19) dz (B.15)Computing the derivative of y and denoting the Gaussian by y , we obtain: y = y + 2 Az y + A z , (B.16a) y (cid:48) = z κ y − Az κ y + 4 z A , (B.16b) yy (cid:48) = − zκ y + Az (cid:18) − z κ (cid:19) y + 2 A z . (B.16c)Recalling that the constant A is exponentially small in κ while y grows onlyas a power law of κ − , the terms multiplied by A will be negligible and thesign of R n ( y ) will be determined also by the leading term: R n ( y ) = I + O ( A ) , (B.16d) I = (cid:90) ( p − f ) y dz − θ (cid:90) m ( z ) y dz, (B.16e) where m ( z ) = v z θ + zκ . (B.16f)Proving instability therefore reduces to show that I is positive for the rangeof parameters and functional forms considered in hypoxic condition. We doso finding a lower bound on the value on I , exploiting the quick decay of thefunction y , whose mass is concentrated in a neighborhood of z = 0. Giventhat p (0) − f (0) > m (0) = 0, provided that m does not grow too fastnear z = 0, we can intuitively see that the major contribution to the integral I will be positive. We will now expand this intuitive argument with a morerigorous calculation.We first focus on I (1)0 = (cid:82) ( p − f ) y dz , the contribution in (B.16e) dueto cell proliferation. We can compute this integral exactly as the integrandcomprises products of exponentials, that can be re-written as the integral ofGaussian distribution: I (1)0 = (cid:34) p ( c ) √ ζ κ e − . g + c ζ Z (cid:18) z − c ζ (cid:19) (B.17a)46 p ( c ) √ ζ κ Z (cid:18) zζ (cid:19) − d f e − k f + c fκ Z (cid:32) √ z − c f ) κ (cid:33)(cid:35) , (B.17b)where 2 ζ , = ( κ g , ) / ( g , + κ ), c = 0 . ζ ) /g and c f = k f κ /
2, while Z is again the normal cumulative distribution function as in Lemma 2.We now focus on the term in (B.20) which depends on m . In this casethe integral can not be computed exactly and we will therefore find a lowerbound for its contribution instead. This is achieved by decomposing the fulldomain [0 ,
1] into two three sub-domains. This will allow us to balance therapid growth of the function m with the quicker decay of y away from z = 0: (cid:90) m y dz = (cid:90) z κ m y dz + (cid:90) z κz κ m y dz + (cid:90) z κ m y dz. (B.18)where z , are positive constants such that 0 < z < z < κ − . Note thatwe have the freedom of choosing their values with the aim of making thequantity in (B.18) as small as possible. It is straightforward to see that m attains its maximum value at z = 1 as both v z and z/κ attain maxima there.We now choose the value of κ so that the derivative of m at z = 0 vanishes: m (cid:48) ( z ) = (cid:18) v (cid:48) z ( z )2 θ + 1 κ (cid:19) ⇒ κ = (cid:115) θ | v (cid:48) z (0) | . (B.19a)However, by definition (see Equation (9c)), under hypoxia, the advectionvelocity v z ( z ) = v − z ( z ) is such that | v (cid:48) z ( z ) | ≤ | v (cid:48) z (0) | for all z ∈ (0 , ω − = 2. Consequently, we have that m ( z ) is a non-decreasingfunction of z , i.e. m (cid:48) ( z ) ≥ (cid:90) m y dz ≤ m ( z κ ) (cid:90) z κ y dz + m ( z κ ) (cid:90) z κz κ y dz + m (1) (cid:90) z κ y dz = m ( z κ ) [ Z ] √ z + m ( z κ ) [ Z ] √ z √ z + 1 κ [ Z ] √ κ √ z ≤ m ( z κ )2 + m ( z κ ) [ Z ] √ z √ z + 1 κ [ Z ] ∞√ z (B.19b)47et us reiterate that we want z and z to be such that m ( z κ ) and m ( z κ )are not too large while [ Z ] √ z and [ Z ] ∞√ z are sufficiently small. In thisway, the growth of m is balanced by the exponential decay of the Gaussianfunction y . In particular, we choose z = √ z = 5 / √
2. Combiningthe above with the estimate from Equation (B.17), we obtain: I > I (1)0 − θm ( z κ )2 − θm ( z κ ) [ Z ] √ z √ z − v (cid:48) z (0) θ [ Z ] ∞√ z = I low . (B.20)We can compute the values of κ and I low associated with the value of the Figure B.20: Plot of the lower bound I low and the standard deviation κ as defined by (B.20)and (B.19a) respectively for parameter regime considered in the paper (note that d f isfixed to its maximum values 0 .
015 as this gives the smaller bound I max ). magnitude V − and steepness ξ − considered in the paper (without loss ofgenerality, we only consider d f = 0 .
015 as I low decreases with d f ). As shownin Figure B.20, for all such values, we have that I low >
0. As I > I low , wetherefore have that generically I is also positive. We estimate A ≤ O (10 − )which justifies us dropping the O ( A ) in (B.16d). Consequently, we concludethat R n ( y ) is positive and so is the quotient R . Hence, in hypoxia, the trivialsteady state is always unstable. 48 cknowledgements The authors wish to thank Professor Philip K. Maini for helpful commentsand feedback on the manuscript. GC is supported by from EPSRC and MRCCentre for Doctoral Training in Systems Approaches to Biomedical Scienceand Cancer Research UK. C.Z. acknowledges Breast Cancer Research Foun-dation (BCRF). P.G.K. acknowledges support from the Leverhulme Trustvia a Visiting Fellowship and thanks the Mathematical Institute of the Uni-versity of Oxford for its hospitality during part of this work.
References [1] F. Bray, J. Ferlay, I. Soerjomataram, et al. , Global cancer statistics2018: GLOBOCAN estimates of incidence and mortality worldwide for36 cancers in 185 countries, CA: A Cancer Journal for Clinicians 68 (6)(2018) 394–424. doi:10.3322/caac.21492 .[2] D. Hanahan, R. A. Weinberg, The hallmarks of cancer, Cell 100 (1)(2000) 57–70. doi:10.1016/S0092-8674(00)81683-9 .[3] T. Reya, S. Morrison, M. a. Clarke, Stem cells, cancer, and cancer stemcells, Nature 414 (2001) 105–111. doi:10.1038/35102167 .[4] K. Rycaj, D. G. Tang, Cancer stem cells and radioresistance, Inter-national Journal of Radiation Biology 90 (8) (2014) 615–621. doi:10.3109/09553002.2014.892227 .[5] M. Baumann, M. Krause, R. Hill, Exploring the role of cancer stemcells in radioresistance, Nature Reviews Cancer 8 (7) (2008) 545–554. doi:10.1038/nrc2419 .[6] A. Schulz, M. F., A. Dubrovska, B. K., Cancer stem cells and ra-dioresistance: Dna repair and beyond, Cancers 11 (2019) 862. doi:10.3390/cancers11060862 .[7] H. Enderling, Cancer stem cells and tumor dormancy, in
Systems Biologyof Tumor Dormancy , H. Enderling, N. Almog and L. Hlatky (Eds.)(2013) 55–71. 498] D. Kong, C. J. Hughes, H. L. Ford, Cellular Plasticity in Breast CancerProgression and Therapy, Frontiers in Molecular Biosciences 7 (2020)72. doi:10.3389/fmolb.2020.00072 .[9] M. Shibata, M. Hoque, Targeting cancer stem cells: A strategy for ef-fective eradication of cancer, Cancers 11 (5) (2019) 732. doi:10.3390/cancers11050732 .[10] V. Snyder, T. C. Reed-Newman, L. Arnold, et al. , Cancer stem cellmetabolism and potential therapeutic targets, Frontiers in Oncology 8(2018) 203. doi:10.3389/fonc.2018.00203 .[11] B. Prager, Q. Xie, S. Bao, et al. , Cancer stem cells: the architects ofthe tumor ecosystem, Cell Stem Cell 24 (2019) 41–53. doi:10.1016/j.stem.2018.12.009 .[12] R. Saga, Y. Matsuya, R. et al. . Takahashi, Analysis of the high-dose-range radioresistance of prostate cancer cells, including cancer stem cells,based on a stochastic model, Journal of Radiation Research 60 (3) (2019)298–307. doi:10.1093/jrr/rrz011 .[13] J. G. Scott, A. Dhawan, A. Hjelmeland, et al. , Recasting the can-cer stem cell hypothesis: unification using a continuum model of mi-croenvironmental forces, Current Stem Cell Reports 5 (2019) 22–30. doi:10.1007/s40778-019-0153-0 .[14] A. Dirkse, A. Golebiewska, T. Buder, et al. , Stem cell-associated het-erogeneity in glioblastoma results from intrinsic tumor plasticity shapedby the microenvironment, Nature Communications 10 (2019) 1787. doi:10.1038/s41467-019-09853-z .[15] H. Soleymani Abyaneh, N. Gupta, A. Alshareef, et al. , Hypoxia In-duces the Acquisition of Cancer Stem-like Phenotype Via Upregulationand Activation of Signal Transducer and Activator of Transcription-3 (STAT3) in MDA-MB-231, a Triple Negative Breast Cancer CellLine, Cancer Microenvironment 11 (2-3) (2018) 141–152. doi:10.1007/s12307-018-0218-0 .[16] G. Fanelli, A. Naccarato, C. Scatena, Recent Advances in Cancer Plas-ticity: Cellular Mechanisms, Surveillance Strategies, and Therapeutic50ptimization, Frontiers in Oncology 10 (2020) 569. doi:10.3389/fonc.2020.00569 .[17] I. Tirosh, A. Venteicher, C. Hebert, et al. , Single cell rna-seq sup-ports a developmental hierarchy in human oligodendroglioma, Nature539 (2016) 309–313. doi:10.1038/nature20123 .[18] A. Venteicher, I. Tirosh, C. Hebert, et al. , Decoupling genetics, lineages,and microenvironment in idh-mutant gliomas by single-cell rna-seq, Sci-ence 355 (6332) (2017) 1391. doi:10.1126/science.aai8478 .[19] R. Chisholm, T. Lorenzi, J. Clairambault, Cell population heterogeneityand evolution towards drug resistance in cancer: Biological and math-ematical assessment, theoretical treatment optimisation, Biochimica etBiophysica Acta 1860 (11) (2016) 2627–2645. doi:10.1016/j.bbagen.2016.06.009 .[20] S. Shen, J. Clairambault, Cell plasticity in cancer cell populations,F1000Research 9 (jun 2020). doi:10.12688/f1000research.24803.1 .[21] A. Hodgkinson, L. Le Cam, D. Trucu, et al. , Spatio-genetic and pheno-typic modelling elucidates resistance and re-sensitisation to treatmentin heterogeneous melanoma, Journal of Theoretical Biology 466 (2019)84–105. doi:10.1016/j.jtbi.2018.11.037 .[22] T. Lorenzi, R. Chisholm, J. Clairambault, Tracking the evolution ofcancer cell populations through the mathematical lens of phenotype-structured equations, Biology Direct 11 (1) (2016) 43. doi:10.1186/s13062-016-0143-4 .[23] A. Lorz, T. Lorenzi, J. Clairambault, et al. , Modeling the Effects ofSpace Structure and Combination Therapies on Phenotypic Heterogene-ity and Drug Resistance in Solid Tumors, Bulletin of Mathematical Bi-ology 77 (1) (2015) 1–22. doi:10.1007/s11538-014-0046-4 .[24] A. Ardaˇseva, R. A. Gatenby, A. R. A. Anderson, et al. , A mathe-matical dissection of the adaptation of cell populations to fluctuat-ing oxygen levels, Bulletin of Mathematical Biology 82 (2020). doi:10.1007/s11538-020-00754-7 .5125] C. Villa, M. Chaplain, T. Lorenzi, Evolutionary dynamics in vascu-larised tumours under chemotherapy, Vietnam Journal of Mathematics(oct 2020). doi:10.1007/s10013-020-00445-9 .[26] F. Thomas, D. Fisher, P. Fort, et al. , Applying ecological and evolution-ary theory to cancer: a long and winding road, Evolutionary Applica-tions 6 (1) (2013) 1–10. doi:10.1111/eva.12021 .[27] T. Lorenzi, R. H. Chisholm, et al. , Dissecting the dynamics of epige-netic changes in phenotype-structured populations exposed to fluctuat-ing environments, Journal of Theoretical Biology 386 (2015) 166–176. doi:10.1016/j.jtbi.2015.08.031 .[28] R. E. A. Stace, T. Stiehl, M. A. Chaplain, et al. , Discrete and continuumphenotype-structured models for the evolution of cancer cell populationsunder chemotherapy, Math. Model. Nat. Phenom 15 (2020) 14. doi:10.1051/mmnp/2019027 .[29] D. Garnier, O. Renoult, M. Alves-Guerra, et al. , Glioblastoma stem-like cells, Metabolic strategy to kill a challenging target, Frontiers inOncology 9 (2019) 118. doi:10.3389/fonc.2019.00118 .[30] S. Liu, Y. Cong, D. Wang, et al. , Breast cancer stem cells transitionbetween epithelial and mesenchymal states reflective of their normalcounterparts, Stem Cell Reports 2 (1) (2013) 78–91. doi:doi:10.1016/j.stemcr.2013.11.009 .[31] F. Pistollato, S. Abbadi, E. Rampazzo, et al. , Intratumoral Hypoxic Gra-dient Drives Stem Cells Distribution and MGMT Expression in Glioblas-toma, Stem Cells 28 (5) (2010) 851–862. doi:10.1002/stem.415 .[32] F. Pistollato, H. Chen, B. et al. . Rood, Hypoxia and HIF1 α Repressthe Differentiative Effects of BMPs in High-Grade Glioma, Stem Cells27 (1) (2009) 7–17. doi:10.1634/stemcells.2008-0402 .[33] T. Lewin, H. Byrne, P. Maini, et al. , The importance of dead materialwithin a tumour on the dynamics in response to radiotherapy, Physicsin medicine and biology 65 (2020) 015007. doi:10.1088/1361-6560/ab4c27 . 5234] T. Lewin, P. K. Maini, E. G. Moros, et al. , A three phase model toinvestigate the effects of dead material on the growth of avascular tu-mours, Mathematical modeling of natural phenomena 15 (2020) 22. doi:10.1051/mmnp/2019039 .[35] E. Markert, A. Vazquez, Mathematical models of cancermetabolism, Cancer & Metabolism 3 (2015) 14. doi:doi:10.1186/s40170-015-0140-6. [36] M. Lee, J. S. Lee, Exploiting tumor cell senescence in anticancer therapy,BMB Reports 47 (2) (2014) 51–59. doi:10.5483/BMBRep.2014.47.2.005 .[37] A. O. Pisco, S. Huang, Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: ’What does not kill me strengthensme’, British Journal of Cancer 112 (11) (2015) 1725–1732. doi:10.1038/bjc.2015.146 .[38] M. H¨ockel, K. Schlenger, M. Mitze, et al. , Hypoxia and radiation re-sponse in human tumors, Seminars in Radiation Oncology 6 (1) (1996)3–9. doi:10.1016/S1053-4296(96)80031-2 .[39] B. S. Sørensen, M. R. Horsman, Tumor Hypoxia: Impact on RadiationTherapy and Molecular Pathways, Frontiers in Oncology 10 (2020) 562. doi:10.3389/fonc.2020.00562 .[40] S. R. McKeown, Defining normoxia, physoxia and hypoxia in tu-mours—implications for treatment response, The British Journal of Ra-diology 87 (1035) (2014) 20130676. doi:10.1259/bjr.20130676 .[41] A. Ayob, T. Ramasamy, Cancer stem cells as key drivers of tumourprogression, Journal of Biomedical Science 25 (1) (mar 2018). doi:10.1186/s12929-018-0426-4 .[42] S. Conley, E. Gheordunescu, P. et al. . Kakarala, Antiangiogenic agentsincrease breast cancer stem cells via the generation of tumor hypoxia,Proceedings of the National Academy of Sciences of the United Statesof America 109 (8) (2012) 2784–2789. doi:10.1073/pnas.1018866109 .[43] J. Lan, H. Lu, D. Samanta, et al. , Hypoxia-inducible factor 1-dependentexpression of adenosine receptor 2B promotes breast cancer stem cell53nrichment, Proceedings of the National Academy of Sciences of theUnited States of America 115 (41) (2018) E9640–E9648. doi:10.1073/pnas.1809695115 .[44] V. Snyder, T. Reed-Newman, L. Arnold, et al. , Cancer stem cellmetabolism and potential therapeutic targets, Frontiers in Oncology 8(jun 2018). doi:10.3389/fonc.2018.00203 .[45] G. Driessens, B. Beck, A. Caauwe, et al. , Defining the mode of tumourgrowth by clonal analysis, Nature 488 (7412) (2012) 527–530. doi:10.1038/nature11344 .[46] S. Bao, Q. Wu, R. McLendon, et al. , Glioma stem cells promote radiore-sistance by preferential activation of the dna damage response, Nature444 (2006) 756–760. doi:10.1038/nature05236 .[47] D. W. Clark, K. Palle, Aldehyde dehydrogenases in cancer stem cells:potential as therapeutic targets, Annals of Translational medicine 4 (24)(2016) 518. doi:10.21037/atm.2016.11.82 .[48] D. M., R. Cho, N. Lobo, et al. , Association of reactive oxygen specieslevels and radioresistance in cancer stem cells, Nature 458 (7239) (2009)780–783. doi:10.1038/nature07733 .[49] G. Vassalli, Aldehyde Dehydrogenases: Not Just Markers, but Func-tional Regulators of Stem Cells, Stem Cells International 2019 (2019). doi:10.1155/2019/3904645 .[50] M. Horsman, L. Mortensen, J. Petersen, et al. , Imaging hypoxia to im-prove radiotherapy outcome, Nature Reviews Clinical Oncology 9 (12)(2012) 674–687. doi:10.1038/nrclinonc.2012.171 .[51] J. Moulder, S. Rockwell, Tumor hypoxia: its impact on cancer therapy,Cancer and Metastasis Reviews 5 (4) (1987) 313–341. doi:10.1007/BF00055376 .[52] K. Leder, K. Pitter, Q. LaPlant, et al. , Mathematical modeling of pdgfdriven glioblastoma reveals optimized radiation dosing schedules, Cell156 (3) (2014) 603–616. doi:10.1016/j.cell.2013.12.029 .5453] S. J. McMahon, The linear quadratic model: usage, interpretation andchallenges, Physics in Medicine & Biology 64 (1) (dec 2018). doi:10.1088/1361-6560/aaf26a .[54] D. J. Brenner, M. Martel, E. Hall, Fractionated regimens for stereotacticradiotherapy of recurrent tumors in the brain, International Journalof Radiation Oncology, Biology, Physics 21 (3) (1991) 819–824. doi:10.1016/0360-3016(91)90703-7 .[55] R. G. Dale, The application of the linear-quadratic dose-effect equationto fractionated and protracted radiotherapy, British Journal of Radiol-ogy 58 (690) (1985) 515–528. doi:10.1259/0007-1285-58-690-515 .[56] H. Hammerlindl, H. Schaider, Tumor cell-intrinsic phenotypic plastic-ity facilitates adaptive cellular reprogramming driving acquired drugresistance, J Cell Commun Signal. 12 (2018) 133–141. doi:10.1007/s12079-017-0435-1 .[57] A. Shiraishi, K. Tachi, N. Essid, et al. , Hypoxia promotes the phenotypicchange of aldehyde dehydrogenase activity of breast cancer stem cells,Cancer Science 108 (3) (2017) 362–372. doi:10.1111/cas.13147 .[58] B. Perthame, Transport equations in biology, Frontiers in Mathematics,Birkh¨auser Basel, 2006.[59] G. F. Webb, Population Models Structured by Age, Size, and SpatialPosition, Springer Berlin Heidelberg, Berlin, Heidelberg, 2008, Ch. 1,pp. 1–49. doi:10.1007/978-3-540-78273-5\_1 .[60] A. Gerisch, M. A. Chaplain, Robust numerical methods for taxis-diffusion-reaction systems: Applications to biomedical problems, Math-ematical and Computer Modelling 43 (1-2) (2006) 49–75. doi:10.1016/j.mcm.2004.05.016 .[61] T. A. Driscoll, N. Hale, L. N. Trefethen, Chebfun Guide, Pafnuty Pub-lications, 2014.[62] C. Calabrese, H. Poppleton, M. Kocak, et al. , A Perivascular Nichefor Brain Tumor Stem Cells, Cancer Cell 11 (1) (2007) 69–82. doi:10.1016/j.ccr.2006.11.020 . 5563] M. Taylor, T. Kron, Consideration of the radiation dose delivered awayfrom the treatment field to patients in radiotherapy, Journal of medicalphysics 36 (2011) 59–71. doi:10.4103/0971-6203.79686 .[64] K. M. Arnold, N. J. Flynn, A. Raben, et al. , The Impact of Radia-tion on the Tumor Microenvironment: Effect of Dose and FractionationSchedules, Cancer Growth and Metastasis 11 (2018) 117906441876163. doi:10.1177/1179064418761639 .[65] B. Fenton, E. Lord, S. Paoni, Effects of Radiation on Tumor Intravascu-lar Oxygenation, Vascular Configuration, Development of Hypoxia, andClonogenic Survival, Radiation Research Society 155 (2) (2001) 360–368. doi:10.1667/0033-7587(2001)155[0360:eoroti]2.0.co;2 .[66] H. Kempf, M. Bleicher, M. Meyer-Hermann, Spatio-Temporal Dynamicsof Hypoxia during Radiotherapy, PLoS ONE 10 (8) (2015) e0133357. doi:10.1371/journal.pone.0133357 .[67] D. Hormuth, A. Jarrett, T. Yankeelov, Forecasting tumor and vascu-lature response dynamics to radiation therapy via image based mathe-matical modeling, Radiation Oncology 15 (2020) 1–14. doi:10.1186/s13014-019-1446-2 .[68] P. Carmeliet, R. K. Jain, Principles and mechanisms of vessel normal-ization for cancer and other angiogenic diseases, Nature Reviews DrugDiscovery 10 (2011) 417–427. doi:10.1038/nrd3455 .[69] R. K. Jain, Antiangiogenesis strategies revisited: From starving tumorsto alleviating hypoxia, Cancer Cell 26 (2014) 605–622. doi:10.1016/j.ccell.2014.10.006 .[70] I. Roeder, R. Lorenz, Asymmetry of stem cell fate and the potentialimpact of the niche observations, simulations, and interpretations, StemCell Reviews 2 (3) (2006) 171–180. doi:10.1007/s12015-006-0045-4 .[71] P. Hahnfeldt, D. Panigrahy, J. Folkman, L. Hlatky, Tumor developmentunder angiogenic signaling: A dynamical theory of tumor growth, treat-ment response, and postvascular dormancy, Cancer Research 59 (19)(1999) 4770–4775. 5672] I. J. Stamper, M. R. Owen, P. K. Maini, H. M. Byrne, Oscil-latory dynamics in a model of vascular tumour growth - implica-tions for chemotherapy, Biology Direct 5 (2010) 27. doi:10.1186/1745-6150-5-27 .[73] M. E. Hubbard, H. M. Byrne, Multiphase modelling of vascular tumourgrowth in two spatial dimensions, Journal of Theoretical Biology 316(2013) 70–89. doi:10.1016/j.jtbi.2012.09.031 .[74] H. M. Byrne, Dissecting cancer through mathematics: From the cellto the animal model, Nature Reviews Cancer 10 (3) (2010) 221–230. doi:10.1038/nrc2808 .[75] P. Macklin, S. McDougall, A. R. A. Anderson, et al. , Multiscalemodelling and nonlinear simulation of vascular tumour growth, Jour-nal of Mathematical Biology 58 (2009) 765–798. doi:10.1007/s00285-008-0216-9 .[76] V. Vavourakis, P. Wijeratne, R. Shipley, , et al. , A validated multiscalein-silico model for mechano-sensitive tumour angiogenesis and growth,PLoS Computational Biology 13 (01 2017). doi:10.1371/journal.pcbi.1005259 .[77] J. Walpole, J. A. Papin, S. M. Peirce, Multiscale computa-tional models of complex biological systems, Annual Review ofBiomedical Engineering 15 (1) (2013) 137–154. doi:10.1146/annurev-bioeng-071811-150104 .[78] A. Ardaˇseva, R. A. Anderson, A, R. A. Gatenby, et al. , Compara-tive study between discrete and continuum models for the evolution ofcompeting phenotype-structured cell populations in dynamical environ-ments, Physical Review E 102 (4) (2020) 042404. arXiv:2004.00914 , doi:10.1103/PhysRevE.102.042404 .[79] B. Franz, M. B. Flegg, et al. , Multiscale reaction-diffusion algorithms:PDE-assisted Brownian dynamics, SIAM Journal on Applied Math-ematics 73 (3) (2013) 1224–1247. arXiv:1206.5860 , doi:10.1137/120882469 . 5780] F. Spill, P. Guerrero, T. Alarcon, et al. , Hybrid approaches for multiple-species stochastic reaction-diffusion models, Journal of ComputationalPhysics 299 (2015) 429–445. doi:10.1016/j.jcp.2015.07.002 .[81] R. A. Gatenby, A. Silva, R. Gillies, et al. , Adaptive therapy, CancerResearch 69 (11) (2009) 4894–4903. doi:10.1158/0008-5472 .[82] E. J. Hall, A. J. Giaccia, Radiobiology for the radiologist [electronicresource], seventh edition Edition, Ebook central, Lippincott Williams& Wilkins, Philadelphia, 2012.[83] P. Sonveaux, F. V´ egran, T. Schroeder, et al. , Targeting lactate-fueledrespiration selectively kills hypoxic tumor cells in mice, The Journal ofClinical Investigation 118 (12) (2008) 3930–42. doi:10.1172/JCI36843 .[84] K. J. Sweeney, A. Swarbrick, et al. , Lack of relationship between CDKactivity and G1 cyclin expression in breast cancer cells, Oncogene 16 (22)(1998) 2865–2878. doi:10.1038/sj.onc.1201814 .[85] U. Del Monte, Does the cell number 109 still really fit one gram of tumortissue?, Cell Cycle 8 (3) (2009) 505–506. doi:10.4161/cc.8.3.7608 .[86] J. W. Boag, Cell respiration as a function of oxygen tension, Inter-national Journal of Radiation Biology 18 (5) (1970) 475–478. doi:10.1080/09553007014551361 .[87] T. D. Lewin, P. Maini, E. Moros, et al. , The Evolution of TumourComposition During Fractionated Radiotherapy: Implications for Out-come, Bulletin of Mathematical Biology 80 (5) (2018) 1207–1235. doi:10.1007/s11538-018-0391-9 .[88] I. M. Pires, Z. Bencokova, M. Milani, et al. , Effects of acute versuschronic hypoxia on DNA damage responses and genomic instability,Cancer Research 70 (3) (2010) 925–935. doi:10.1158/0008-5472.CAN-09-2715 .[89] S. Ricardo, A. Vieira, R. Gerhard, et al. , Breast cancer stem cell mark-ers CD44, CD24 and ALDH1: Expression distribution within intrinsicmolecular subtype, Journal of Clinical Pathology 64 (11) (2011) 937–944. doi:10.1136/jcp.2011.090456 .5890] N. Ng, K. Purshouse, I. Foskolou, , et al. , Challenges to dna replicationin hypoxic conditions, The FEBS Journal 285 (2017) 1563–1571. doi:10.1111/febs.14377doi:10.1111/febs.14377