A mathematical study of the influence of hypoxia and acidity on the evolutionary dynamics of cancer
AA mathematical study of the influence of hypoxia and acidityon the evolutionary dynamics of cancer ∗ Giada Fiandaca † Marcello Delitala ‡ Tommaso Lorenzi § Abstract
Hypoxia and acidity act as environmental stressors promoting selection for cancer cellswith a more aggressive phenotype. As a result, a deeper theoretical understanding of thespatio-temporal processes that drive the adaptation of tumour cells to hypoxic and acidic mi-croenvironments may open up new avenues of research in oncology and cancer treatment. Wepresent a mathematical model to study the influence of hypoxia and acidity on the evolution-ary dynamics of cancer cells in vascularised tumours. The model is formulated as a system ofpartial integro-differential equations that describe the phenotypic evolution of cancer cells inresponse to dynamic variations in the spatial distribution of three abiotic factors that are keyplayers in tumour metabolism: oxygen, glucose and lactate. The results of numerical simula-tions of a calibrated version of the model based on real data recapitulate the eco-evolutionaryspatial dynamics of tumour cells and their adaptation to hypoxic and acidic microenviron-ments. Moreover, such results demonstrate how nonlinear interactions between tumour cellsand abiotic factors can lead to the formation of environmental gradients which select for cellswith phenotypic characteristics that vary with distance from intra-tumour blood vessels, thuspromoting the emergence of intra-tumour phenotypic heterogeneity. Finally, our theoreticalfindings reconcile the conclusions of earlier studies by showing that the order in which re-sistance to hypoxia and resistance to acidity arise in tumours depend on the ways in whichoxygen and lactate act as environmental stressors in the evolutionary dynamics of cancer cells.
Keywords:
Mathematical oncology, Intra-tumour heterogeneity, Eco-evolutionary dynamics,Vascularised tumours, Partial integro-differential equations
Cancer is a dynamic disease, the characteristics of which are constantly evolving. This is reflectedin the fact that the genotypic and phenotypic properties of cancer cells may change across spaceand time within the same tumour, and the dynamics of tumours with the same histological featuresare still likely to vary across patients. Moreover, since the same cancer clones may arise throughdifferent evolutionary pathways, the fact that two tumours have a similar clonal composition ata given point in time does not necessarily indicate that they share similar evolutionary histories,and does not rule out the possibility that their future evolution will diverge significantly [46].These sources of variability within and between tumours provide the substrate for the emergenceand development of intra- and inter-tumour heterogeneity, which are major obstacles to cancereradication [29, 48]. ∗ This work was supported by the MIUR grant “Dipartimenti di Eccellenza 2018-2022”. † Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Corso Duca degli Abruzzi, 24,Torino, Italy ‡ Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Corso Duca degli Abruzzi, 24,Torino, Italy § Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, Corso Duca degli Abruzzi, 24,Torino, Italy a r X i v : . [ q - b i o . T O ] S e p linical evidence suggests that cancer cells and the tumour microenvironment mutually shapeeach other [24]. This supports the idea that tumours can be seen as evolving ecosystems wherecancer cells with different phenotypic characteristics proliferate, die, undergo genotypic and pheno-typic changes, and compete for space and resources under the selective pressure exerted by the var-ious components of the tumour microenvironment [23, 28, 32, 36, 39, 40, 50, 52, 62]. In this light,intra-tumour phenotypic heterogeneity can be regarded as the outcome of an eco-evolutionaryprocess in which spatial variability of the concentration of abiotic factors ( i.e. substrates andmetabolites) across the tumour supports the formation of distinct ecological niches whereby cellswith different phenotypic characteristics may be selected [14, 26, 31].In normal tissues, cells produce the energy required to sustain their proliferation via oxidativephosphorylation ( i.e. they rely on oxygen as their primary source of energy) and turn to glycolysisonly when oxygen is scarce. In tumours, the presence of hypoxic regions ( i.e. regions where theoxygen levels are below the physiological ones) induces cells to transiently switch to a glycolyticmetabolic phenotype ( i.e. to rely on glucose as their primary source of energy) [63]. Cancer cellseventually acquire such a glycolytic phenotype and express it also in aerobic conditions, leadingto the so-called Warburg effect [33]. The interplay between the high glycolytic rate of cancercells and low perfusion in tumours brings about accumulation of lactate ( i.e. a waste product ofglycolysis), which causes acidity levels in the tumour microenvironment to rise ( i.e. the pH leveldrops) [60].Since hypoxia and acidity act as environmental stressors promoting selection for cancer cellswith a more aggressive phenotype [57, 60], an in-depth theoretical understanding of the spatio-temporal processes that drive the adaptation of tumour cells to hypoxic and acidic microenviron-ments may open up new avenues of research in oncology and cancer treatment [47]. In this regard,mathematical models can be an important source of support to cancer research, as they enableextrapolation beyond scenarios which can be investigated through experiments and may revealemergent phenomena that would otherwise remain unobserved [4, 13, 16, 17, 30]. For instance,in their pioneering paper [25], Gatenby and Gawlinski used a reaction-diffusion system to explorehow nonlinear interactions between cancer cells and abiotic components of the tumour microen-vironment may shape tumour growth. The Gatenby-Gawlinski model has recently been extendedin [58], in order to take into account the presence of cells with different phenotypic characteristicswithin the tumour. Hybrid cellular automaton models have been employed to study the impactof hypoxia and acidity on tumour growth and invasion [3, 5, 26, 34, 57]. A mechanical model oftumour growth whereby cells are allowed to switch between aerobic and anaerobic metabolismwas presented in [9]. Integro-differential equations and partial integro-differential have been usedin [7, 15, 42, 44, 64] to investigate the ecological role of hypoxia in the development of intra-tumourphenotypic heterogeneity.In this paper, we complement these earlier studies by presenting a mathematical model to studythe influence of hypoxia and acidity on the evolutionary dynamics of cancer cells in vascularisedtumours. The model comprises a system of partial integro-differential equations that describe thephenotypic evolution of cancer cells in response to dynamic variations in the spatial distributionof three abiotic factors that are key players in tumour metabolism: oxygen, glucose and lactate.The remainder of the paper is organised as follows. In Section 2, we present the model equa-tions and the underlying modelling assumptions. In Section 3, we summarise the main results ofnumerical simulations of the model and discuss their biological implications. Section 4 concludesthe paper and provides a brief overview of possible research perspectives. We consider a one-dimensional region of vascularised tissue of length L >
0. We describe thespatial position of every tumour cell in the tissue region by a scalar variable x ∈ [0 , L] and weassume a blood vessel to be present at x = 0 ( cf. the schematic in Figure 1a). Moreover, buildingupon the modelling framework developed in [15, 42, 44, 64], we model the phenotypic state ofevery cell by a vector y = ( y , y ) ∈ [0 , ( cf. the schematics in Figure 1b). Here, y ∈ [0 , e.g. the LAMP2 gene),while y ∈ [0 ,
1] represents the normalised level of expression of a hypoxia-resistant gene ( e.g. theGLUT-1 gene) [19, 26].We describe the phenotypic distribution of tumour cells at position x and time t ∈ [0 , T], withT >
0, by means of the local population density function n ( t, x, y ) ( i.e. the local phenotypicdistribution of tumour cells). We define the cell density ρ ( t, x ), the local mean level of expressionof the acidity-resistant gene µ ( t, x ) and the local mean level of expression of the hypoxia-resistantgene µ ( t, x ) as ρ ( t, x ) := (cid:90) [0 , n ( t, x, y ) d y , µ i ( t, x ) := 1 ρ ( t, x ) (cid:90) [0 , y i n ( t, x, y ) d y (1)for i = 1 ,
2. Moreover, we define the phenotypic distribution of tumour cells across the wholetissue region f ( t, y ) as the mean value of n ( t, x, y ) on the interval [0 , L], i.e. f ( t, y ) := 1L (cid:90) L0 n ( t, x, y ) d x. (2)Similarly, we define the levels of expression of the acidity-resistant gene and the hypoxia-resistantgene across the whole tissue region as the mean values of µ ( t, x ) and µ ( t, x ) on the interval [0 , L],respectively, i.e. ν ( t ) := 1L (cid:90) L0 µ ( t, x ) d x and ν ( t ) := 1L (cid:90) L0 µ ( t, x ) d x. (3)The local concentrations of oxygen, glucose and lactate at position x and time t are denotedby S o ( t, x ), S g ( t, x ) and S l ( t, x ), respectively.Figure 1: a) Schematic of the spatial domain defined as a one-dimensional region of vascularised tissue of lengthL. A blood vessel is assumed to be present at x = 0. b) Schematics illustrating the relationships between the valuesof the variables y and y modelling the phenotypic state of tumour cells and the levels of resistance to hypoxiaand acidosis. We describe the dynamics of tumour cells through the following balance equation for the localpopulation density function n ( t, x, y ) ∂n∂t = β n ∂ n∂x (cid:124) (cid:123)(cid:122) (cid:125) undirected, randomcell movement + θ ∆ y n (cid:124) (cid:123)(cid:122) (cid:125) spontaneousphenotypic changes + R ( S o , S g , S l , ρ, y ) n (cid:124) (cid:123)(cid:122) (cid:125) proliferation and death , (4)3ith ( t, x, y ) ∈ (0 , T] × (0 , L) × (0 , , subject to suitable initial conditions. We complement (4)with zero-flux boundary conditions at x = 0 and x = L ( i.e. we assume that cells cannot leavethe tissue region) and zero-flux boundary conditions on the boundary of the square [0 , ( i.e. weassume that cells cannot have normalised levels of gene expression smaller than 0 or larger than1). The first term on the right-hand side of the partial integro-differential equation (4) models theeffect of undirected, random movement, which is described through Fick’s first law of diffusion withdiffusivity β n >
0, while the second term models the effect of heritable, spontaneous phenotypicchanges, which occur at rate θ >
0. The function R ( S o , S g , S l , ρ, y ) represents the fitness of tumourcells in the phenotypic state y at position x and time t under the environmental conditions givenby the concentrations of abiotic factors S o ( t, x ), S g ( t, x ) and S l ( t, x ), and the cell density ρ ( t, x )( i.e. R ( S o , S g , S l , ρ, y ) is the phenotypic fitness landscape of the tumour). We use the followingdefinition R ( S o , S g , S l , ρ, y ) := P ( S o , S g , y ) (cid:124) (cid:123)(cid:122) (cid:125) proliferation anddeath due tooxygen-driven selection − D ( S l , y ) (cid:124) (cid:123)(cid:122) (cid:125) death due tolactate-driven selection − d ( ρ ) . (cid:124) (cid:123)(cid:122) (cid:125) death due tocompetitionfor space (5)Here, the function P ( S o , S g , y ) is the rate at which cells with level of expression y of the hypoxia-resistant gene proliferate via oxidative phosphorylation and glycolysis, and die due to oxygen-driven selection ( i.e. P ( S o , S g , y ) is a net proliferation rate). The function D ( S l , y ) is the rateat which cells with level of expression y of the acidity-resistant gene die due to lactate-drivenselection. The function d ( ρ ) models the rate of cell death due to competition for space associatedwith saturation of the cell density. Based on the theoretical results and experimental data presented in [36, 63], we focus on a scenariocorresponding to the following biological assumptions.
Assumption 1.
There exist two threshold levels of the oxygen concentration O M > O m > such that the environment surrounding the cells is: hypoxic if S o ≤ O m ; moderately oxygenated if O m < S o < O M ; normoxic (i.e. well oxygenated) if S o ≥ O M . Assumption 2.
Cells proliferate at a rate that depends on the concentrations of oxygen andglucose. Moreover, the trade-off between increase in cell death associated with sensitivity to hypoxiaand decrease in cell proliferation associated with acquisition of resistance to hypoxia results in theexistence of a level of expression of the hypoxia-resistant gene which is the fittest in that: a lowerlevel of gene expression would correlate with a lower resistance to hypoxia, and thus a higher deathrate; a higher level of gene expression would correlate with a larger fitness cost, and thus a lowerproliferation rate. Cells with levels of gene expression that are closer to the fittest one are morelikely to survive than the others. Hence, the farther the gene expression level of a cell is from thefittest one, the more likely is that the cell will die due to a form of oxygen-driven selection.
Assumption 3.
In normoxic environments (i.e. when S o ≥ O M ), the energy required for cell pro-liferation is produced via oxidative phosphorylation and the cell proliferation rate is a monotonicallyincreasing function of the concentration of oxygen. In hypoxic environments (i.e. when S o ≤ O m ),the energy required for cell proliferation is produced via glycolysis and the cell proliferation rateis a monotonically increasing function of the concentration of glucose. In moderately-oxygenatedenvironments (i.e. when O m < S o < O M ), the energy required for cell proliferation is producedvia both oxidative phosphorylation and glycolysis. Moreover, the cell proliferation rate is a mono-tonically increasing function of the concentrations of oxygen and glucose, and lower values of theoxygen concentration correlate with a greater tendency of the cells to proliferate via glycolysis. Assumption 4.
The fittest level of expression of the hypoxia-resistant gene (i.e. the gene asso-ciated with the variable y ) may vary with the oxygen concentration. In particular: in normoxic nvironments (i.e. when S o ≥ O M ), the fittest level of gene expression is the minimal one (i.e. y = 0 ); in hypoxic environments (i.e. when S o ≤ O m ) the fittest level of gene expression is themaximal one (i.e. y = 1 ); in moderately-oxygenated environments (i.e. when O m < S o < O M ),the fittest level of gene expression is a monotonically decreasing function of the oxygen concentra-tion (i.e. it decreases from y = 1 to y = 0 when the oxygen concentration increases). Under Assumptions 1 and 2, we define the net proliferation rate P ( S o , S g , y ) as P ( S o , S g , y ) := p o ( S o ) (cid:124) (cid:123)(cid:122) (cid:125) proliferation viaoxidative phosphorylation + p g ( S o , S g ) (cid:124) (cid:123)(cid:122) (cid:125) proliferation viaglycolysis − η o (cid:0) y − ϕ o ( S o ) (cid:1) (cid:124) (cid:123)(cid:122) (cid:125) death due tooxygen-driven selection . (6)In (6), the function p o ( S o ) models the rate of cell proliferation via oxidative phosphorylation, whilethe function p g ( S o , S g ) models the rate of cell proliferation via glycolysis. Furthermore, the thirdterm in the definition given by (6) is the rate of death induced by oxygen-driven selection. Here,the parameter η o > ϕ o ( S o ) is the fittest level of expression of the hypoxia-resistant gene under theenvironmental conditions given by the oxygen concentration S o . Remark 1. In (6) , the distance between y and ϕ o ( S o ) is computed as (cid:0) y − ϕ o ( S o ) (cid:1) . Alterna-tively, one could compute such a distance as (cid:12)(cid:12) y − ϕ o ( S o ) (cid:12)(cid:12) . However, we have chosen (cid:0) y − ϕ o ( S o ) (cid:1) over (cid:12)(cid:12) y − ϕ o ( S o ) (cid:12)(cid:12) because, as discussed in [7, 41], it leads to a smoother fitness function which iscloser to the approximate fitness landscapes which can be inferred from experimental data throughregression techniques. Under Assumptions 3 and 4, we use the definitions of the functions p o ( S o ), p g ( S o , S g ) and ϕ o ( S o ) given hereafter p o ( S o ) := γ o S o α o + S o w ( S o ) , p g ( S o , S g ) := γ g S g α g + S g (cid:0) − w ( S o ) (cid:1) , (7)with w ( S o ) := S o ≥ O M − O M − S o O M − O m O m < S o < O M S o ≤ O m (8)and ϕ o ( S o ) := S o ≥ O M O M − S o O M − O m O m < S o < O M S o ≤ O m . (9)In (7), the parameters γ o > γ g > α o > α g > w ( S o ) defined via (8)ensures that Assumption 3 is satisfied, while definition (9) of ϕ o ( S o ) is such that Assumption 4 issatisfied ( cf. the plot in Figure 2a). Based on theoretical results and experimental data presented in [57], we focus on a scenariocorresponding to the following biological assumptions.5 ssumption 5.
There exist two threshold levels of the lactate concentration L M > L m > suchthat the environment surrounding the cells is: mildly acidic if S l ≤ L m ; moderately acidic if L m < S l < L M ; highly acidic if S l ≥ L M . Assumption 6.
Cells die at a rate that depends on the concentration of lactate. Moreover, thetrade-off between increase in cell death associated with sensitivity to acidity and decrease in cellproliferation associated with acquisition of resistance to acidity results in the existence of a level ofexpression of the acidity-resistant gene which is the fittest in that: a lower level of gene expressionwould correlate with a lower resistance to acidity, and thus a higher death rate; a higher level ofgene expression would correlate with a larger fitness cost, and thus a lower proliferation rate. Cellswith levels of gene expression that are closer to the fittest one are more likely to survive than theothers. Hence, the farther the gene expression level of a cell is from the fittest one, the more likelyis that the cell will die due to a form of lactate-driven selection.
Assumption 7.
The fittest level of expression of the acidity-resistant gene (i.e. the gene associatedwith the variable y ) may vary with the lactate concentration. In particular: in mildly-acidicenvironments (i.e. when S l ≤ L m ), the fittest level of gene expression is the minimal one (i.e. y = 0 ); in highly-acidic environments (i.e. when S l ≥ L M ) the fittest level of gene expression isthe maximal one (i.e. y = 1 ); in moderately-acidic environments (i.e. when L m < S l < L M ), thefittest level of gene expression is a monotonically increasing function of the lactate concentration(i.e. it increases from y = 0 to y = 1 when the lactate concentration increases). Under Assumptions 5 and 6, we define the rate of cell death due to lactate-driven selection D ( S l , y ) as D ( S l , y ) := η l (cid:0) y − ϕ l ( S l ) (cid:1) . (10)In (10), the parameter η l > ϕ l ( S l ) is the fittest level of expression of the acidity-resistant gene underthe environmental conditions given by the lactate concentration S l . Considerations analogous tothose made in Remark 1 on the term (cid:0) y − ϕ o ( S o ) (cid:1) in (6) apply to the term (cid:0) y − ϕ l ( S l ) (cid:1) in (10).Finally, we use the definition of the function ϕ l ( S l ) given hereafter ( cf. the plot in Figure 2b), sothat Assumption 7 is satisfied: ϕ l ( S l ) := S l ≤ L m S l − L m L M − L m L m < S l < L M S l ≥ L M . (11) We define the rate of cell death due to competition for space associated with saturation of the celldensity as d ( ρ ) := κ ρ. (12)Here, the proportionality constant κ > Oxygen and glucose are consumed by tumour cells, while lactate is produced by tumour cells as awaste product of glycolysis. Moreover, oxygen, glucose and lactate diffuse in space and decay over6igure 2: a) Plot of the fittest level of expression of the hypoxia-resistant gene ϕ o ( S o ) defined via (9). The vertical,dashed lines highlight the oxygen levels O m and O M . Hence, the white region corresponds to normoxic conditions,the pale-blue region corresponds to moderately-oxygenated environments and the blue region corresponds to hypoxicconditions. b) Plot of the fittest level of expression of the acidity-resistant gene ϕ l ( S l ) defined via (11). Thevertical, dashed lines highlight the lactate levels L m and L M . Hence, the white region corresponds to mildly-acidicconditions, the pale-green region corresponds to moderately-acidic conditions and the green region corresponds tohighly-acidic conditions. time. Hence, their dynamics are governed by the following balance equations for the functions S o ( t, x ), S g ( t, x ) and S l ( t, x ), respectively, ∂S o ∂t = β o ∂ S o ∂x (cid:124) (cid:123)(cid:122) (cid:125) diffusion − λ o S o (cid:124) (cid:123)(cid:122) (cid:125) naturaldecay − ζ o p o ( S o ) ρ, (cid:124) (cid:123)(cid:122) (cid:125) consumption by tumour cells (13) ∂S g ∂t = β g ∂ S g ∂x (cid:124) (cid:123)(cid:122) (cid:125) diffusion − λ g S g (cid:124) (cid:123)(cid:122) (cid:125) naturaldecay − ζ g p g ( S o , S g ) ρ (cid:124) (cid:123)(cid:122) (cid:125) consumption by tumour cells (14)and ∂S l ∂t = β l ∂ S l ∂x (cid:124) (cid:123)(cid:122) (cid:125) diffusion − λ l S l (cid:124)(cid:123)(cid:122)(cid:125) naturaldecay + ζ l p g ( S o , S g ) ρ, (cid:124) (cid:123)(cid:122) (cid:125) production by tumour cells (15)with ( t, x ) ∈ (0 , T] × (0 , L), subject to suitable boundary conditions (see considerations below) andinitial conditions.In (13)-(15), the parameters β o > β g > β l > λ o > λ g > λ l > ρ and the rate of cell proliferation via oxidative phosphorylation p o ( S o ), which is definedvia (7). The parameter ζ o > ρ and the rate of cell proliferation via glycolysis p g ( S o , S g )defined via (7). The constant of proportionality is the conversion factor ζ l > x = 0 S o ( t,
0) = S o , S g ( t,
0) = S g , S l ( t,
0) = S l for all t > , (16)where S o > S g > S l > S o > O M and S l < L m . Moreover, we assume that far from the blood vessel theconcentrations of oxygen and glucose drop to some lower values 0 < S o < S o and 0 < S g < S g ,which correspond to the levels of oxygen and glucose which are typically observed in regionsdistant from blood vessels. In particular, we have S o < O m . Furthermore, we model abnormalaccumulation of lactate, which is expected to occur far from blood vessels, imposing zero-fluxboundary conditions. Therefore, we complement (13)-(15) with the following boundary conditionsat x = L S o ( t, L) = S o , S g ( t, L) = S g , ∂S l ( t, L) ∂x = 0 for all t > . (17) In this section, we present the results of numerical simulations of the mathematical model definedby (4) coupled with (13)-(15) and we discuss their biological relevance. First, we describe the set-up of numerical simulations (see Section 3.1). Next, we present a sample of numerical solutionsthat summarise the spatial dynamics of tumour cells and abiotic factors (see Section 3.2). Then,we report on the results of numerical simulations showing the evolutionary dynamics of tumourcells and the emergence of phenotypic heterogeneity (see Section 3.3). Finally, we present theresults of numerical simulations that reveal the existence of alternative evolutionary pathwaysthat may lead to the development of resistance to hypoxia and acidity in vascularised tumours(see Section 3.4).
Numerical simulations are carried out assuming L = 400 µm , which is chosen coherently withexperimental data reported in [54], and t ∈ [0 , T], where the final time T > t = T. Initial conditions.
We consider (13), (14) and (15) subject, respectively, to the following initialconditions S o (0 , x ) = S o ( x ) , S g (0 , x ) = S g ( x ) and S l (0 , x ) = S l ( x ) ≡ S l . (18)Here, the functions S o ( x ) and S g ( x ) (see Figure 3) are defined in such a way as to match theexperimental equilibrium distributions of oxygen and glucose presented in [54, Fig. 2], while S l is the same small parameter used in (16), i.e. S l < L m . Initial conditions (18) correspondto a situation in which the initial distributions of oxygen and glucose match with experimentalequilibrium distributions of such abiotic factors and lactate is present at a uniform level which isbelow the threshold level L m , that is, the level below which the environment surrounding the cellsis mildly acidic and the fittest level of expression of the acidity-resistant gene is the minimal one.Moreover, we complement (4) with the following initial condition n (0 , x, y ) = n ( x, y ) := 200 exp (cid:18) − ( x − . − | y − | . (cid:19) , (19)which corresponds to a biological scenario in which at the initial time t = 0 most tumour cells areconcentrated near the blood vessel and are characterised by the minimal expression level of boththe hypoxia-resistant gene and the acidity-resistant gene.8 oundary conditions. We use the following values of the parameters S o , S g and S l in (16) S o = 2 . × − g/cm , S g = 1 . × − g/cm , S l = 10 − g/cm (20)and the following values of the parameters S o and S g in (17) S o = 2 × − g/cm , S g = 1 . × − g/cm . (21)The values of S o and S g correspond to the average physiological levels of oxygen and glucose inproximity to blood vessels reported in [54]. Moreover, the values of S o and S g correspond to the0 .
1% of S o and the 1% of S g , respectively. This is because, based on experimental data reportedin [54], we expect the concentrations of oxygen and glucose at 400 µm from the blood vessel ( i.e. at x = L) to drop, respectively, below the 0 .
1% and the 1% of their value near the blood vessel.Figure 3:
Plots of the initial distribution of oxygen S o ( x ) (blue curve, axis on the left) and the initial distributionof glucose S g ( x ) (red curve, axis on the right) which are defined in such a way as to match the experimentalequilibrium distributions of oxygen and glucose presented in [54, Fig. 2]. The vertical, dashed lines highlight thepoints x O M and x O m such that S o ( x O M ) = O M and S o ( x O m ) = O m . Hence, the white region corresponds tonormoxic conditions, the pale-blue region corresponds to moderately-oxygenated environments and the blue regioncorresponds to hypoxic conditions. Parameter values.
Unless otherwise explicitly stated, we use the values of the model parame-ters listed in Table 1, which are chosen to be consistent with the existing literature, except for thevalues of the parameters η o , η l , λ l and ζ l that are model specific in that we could not find themin the literature and are defined on the basis of the following considerations. The value of theconversion factor for lactate production ζ l is chosen to be the same as the value of the conversionfactor for glucose consumption ζ g . Furthermore, the value of the rate of natural decay of lactate λ l is such that the distribution of lactate at numerical equilibrium ( i.e. the graph of S l (T , x )) issimilar to the lactate distributions reported in [54]. Finally, the values of the selection gradients η o and η l are chosen with exploratory aim and a systematic sensitivity analysis of these parametersis performed in Section 3.4. Numerical methods.
Numerical solutions are constructed using a uniform discretisation ofthe interval [0 , L] as the computational domain of the independent variable x and a uniformdiscretisation of the square [0 , as the computational domain of the independent variable y .We also discretise the interval [0 , T] with a uniform step. The method for constructing numericalsolutions is based on an explicit finite difference scheme in which a three-point and a five-pointstencils are used to approximate the diffusion terms in x and y , respectively, and an implicit-explicit finite difference scheme is used for the reaction terms [37, 45]. All numerical computationsare performed in Matlab . 9able 1:
Parameter values used in numerical simulations.Par. Biological meaning Value Units Ref. β o Diffusion coefficient of oxygen 1 . × − cm s − [54] β g Diffusion coefficient of glucose 1 . × − cm s − [54] β l Diffusion coefficient of lactate 1 . × − cm s − [54] β n Cell motility 10 − cm s − [64] θ Rate of phenotypic changes 10 − s − [64] λ o Rate of natural decay of oxygen 1 . × − s − [54] λ g Rate of natural decay of glucose 2 . × − s − [54] λ l Rate of natural decay of lactate 5 × − s − ad hoc γ o Max. prolif. rate via oxidative phosphorylation 3 . × − s − [54] α o Michaelis-Menten constant of oxygen 6 . × − g cm − [54] γ g Max. prolif. rate via glycolysis 3 . × − s − [54] α g Michaelis-Menten constant of glucose 9 × − g cm − [54] η o Selection gradient related to oxygen 3 . × − s − ad hoc η l Selection gradient related to lactate 10 − s − ad hoc κ Rate of cell death due to competition for space 2 × − cm s − cells − [42] ζ o Conversion factor for oxygen consumption 10 − g cells − [42] ζ g Conversion factor for glucose consumption 10 − g cells − [42] ζ l Conversion factor for lactate production 10 − g cells − ad hoc O m Threshold level of oxygen for hypoxic env. 8 . × − g cm − [63] O M Threshold level of oxygen for normoxic env. 4 . × − g cm − [63] L m Threshold level of lactate for mildly-acidic env. 2 × − g cm − [57] L M Threshold level of lactate for highly-acidic env. 7 . × − g cm − [57]L Max. distance from blood vessel 400 µm [54] The dynamics of the density of tumour cells and the concentrations of abiotic factors are illustratedby the plots in Figure 4. In summary, the cell density and the concentration of lactate at successivetimes ( i.e. the graphs of ρ ( t, x ) for four different values of t and S l ( t, x ) for three different valuesof t ) are displayed in Figure 4a and Figure 4c, respectively, while the concentrations of oxygenand glucose at time T ( i.e. the graphs of S o (T , x ) and S g (T , x )) are displayed in Figure 4b.The curves in Figure 4a show that cell movement in concert with cell proliferation and deathleads tumour cells, which are initially concentrated near the blood vessel at x = 0 ( cf. the initialcondition n ( x, y ) defined via (19)), to spread into the surrounding tissue ( i.e. ρ ( t, x ) behaves likean invasion front). At every position, the cell density grows until a saturation value, which varieswith the distance from the blood vessel, is reached ( i.e. ρ ( t, x ) appears to converge to a form ofgeneralised transition wave [10, 11]). Moreover, the curves in Figure 4b show that the reaction-diffusion dynamics of oxygen and glucose, along with the inflow through the blood vessel and theconsumption by tumour cells, lead the concentrations of such abiotic factors to converge to somestable values which decrease as the distance from the blood vessel increases ( i.e. at t = T, S o ( t, x )and S g ( t, x ) appear to be at numerical equilibrium and are monotonically decreasing functionsof x ). Furthermore, the curves in Figure 4c show that the interplay of production by tumourcells, outflow through the blood vessel and reaction-diffusion dynamics leads the concentration oflactate to grow at every position until a saturation value, which varies with the distance from theblood vessel, is reached ( cf. the graph of S l (T , x )).Notice that the distributions of oxygen and glucose at the final time T are close to the initialdistributions S o ( x ) and S g ( x ) displayed in Figure 3. This is to be expected. In fact, since S o ( x ) and S g ( x ) are defined in such a way as to match experimental equilibrium distributions of oxygen andglucose, under the biologically informed parameter values ( cf. Table 1) and boundary conditions( cf. (16), (17) and (20), (21)) used here, the concentrations S o ( t, x ) and S g ( t, x ) reach quicklynumerical equilibrium. We verified via additional numerical simulations (results not shown) that,as one would expect, the concentrations of oxygen and glucose at numerical equilibrium do notdepend on the choice of the initial conditions. 10he dashed lines in Figure 4 highlight the spatial positions x O M and x O m at which theoxygen concentration at time T crosses, respectively, the threshold values O M and O m ( i.e. S o (T , x O M ) = O M and S o (T , x O m ) = O m ). Hence, the white region ( i.e. the interval [0 , x O M ]), thepale-blue region ( i.e. the interval ( x O M , x O m )) and the blue region ( i.e. the interval [ x O m , L]) cor-respond to normoxic, moderately-oxygenated and hypoxic environmental conditions, respectively.Under normoxic conditions (white regions), cells proliferate via oxidative phosphorylation at rate γ o S o α o + S o , while cell proliferation occurs via glycolysis at rate γ g S g α g + S g in hypoxic conditions (blueregions). Moreover, in moderately-oxygenated environments (pale-blue regions), cells proliferatevia oxidative phosphorylation at rate γ o S o α o + S o w ( S o ) and via glycolysis at rate γ g S g α g + S g (cid:0) − w ( S o ) (cid:1) ,where w ( S o ) decreases from 1 to 0 when S o decreases from O M to O m ( cf. the definition of w ( S o )given by (8)). Coherently with experimental evidence indicating that using glycolysis to produceenergy required for cell proliferation is less efficient than employing oxidative phosphorylation, thebiologically informed parameter values considered here are such that the cell proliferation rate innormoxic conditions is higher than the cell proliferation rate in moderately-oxygenated environ-ments, which in turn is higher than the cell proliferation rate in hypoxic conditions ( cf. the plotin Figure 5, blue curve). As a result, the saturation value of the cell density decreases with thedistance from the blood vessel ( i.e. ρ (T , x ) is a monotonically decreasing function of x ). Further-more, the production rate of lactate by tumour cells is proportional to the rate of proliferation viaglycolysis and, therefore, it increases with the distance from the blood vessel at x = 0 ( cf. the plotin Figure 5, red curve). As a result, the saturation value of the lactate concentration increaseswith the distance from the blood vessel and, in agreement with the lactate distributions reportedin [54], S l (T , x ) is a monotonically increasing function of x .Figure 4: a) Plots of the cell density ρ ( t, x ) at four successive time instants (yellow, orange, red and burgundylines). The burgundy line highlights ρ (T , x ). b) Plots of the concentrations of oxygen S o (T , x ) (blue line) andglucose S g (T , x ) (red line). c) Plots of the concentration of lactate S l ( t, x ) at three successive time instants (light-green, green and dark-green lines). The dark-green line highlights S l (T , x ). In every panel, the vertical, dashed lineshighlight the points x O M and x O m such that S o (T , x O M ) = O M and S o (T , x O m ) = O m . Hence, the white regioncorresponds to normoxic conditions, the pale-blue region corresponds to moderately-oxygenated environments andthe blue region corresponds to hypoxic conditions. As discussed in the previous section, reaction-diffusion dynamics of abiotic factors and mutualinteractions between abiotic factors and tumour cells lead to the emergence of spatial variations inthe concentrations of oxygen and lactate – i.e. the oxygen concentration S o (T , x ) is a monotoni-cally decreasing function of x while the lactate concentration S l (T , x ) is a monotonically increasingfunction of x ( cf. plots in Figure 4b and Figure 4c). Spatial variability of oxygen and lactate con-11igure 5: Plots of the cell proliferation rate p o ( S o (T , x )) + p g ( S o (T , x ) , S g (T , x )) (blue curve, axis on the left)and the production rate of lactate ζ l p g ( S o (T , x ) , S g (T , x )) ρ (T , x ) (red curve, axis on the right) defined via (7).The vertical, dashed lines highlight the points x O M and x O m such that S o (T , x O M ) = O M and S o (T , x O m ) = O m .Hence, the white region corresponds to normoxic conditions, the pale-blue region corresponds to moderately-oxygenated environments and the blue region corresponds to hypoxic conditions. centrations can lead to the formation of environmental gradients resulting in the selection for cellswith phenotypic characteristics that vary with distance from the blood vessel, thus promotingthe emergence of intra-tumour phenotypic heterogeneity. This is demonstrated by the numericalresults in Figure 6.The dashed lines in Figure 6a and Figure 6b highlight the fittest levels of expression of thehypoxia-resistant gene (see Figure 6a) and the acidity-resistant gene (see Figure 6b) at time T[ i.e. the graphs of ϕ o ( S o (T , x )) and ϕ l ( S l (T , x ))], while the solid lines display the local meanlevels of expression of the hypoxia-resistant gene (see Figure 6a) and the acidity-resistant gene(see Figure 6b) at four successive times ( i.e. the graphs of µ ( t, x ) and µ ( t, x ) for four differentvalues of t ). In analogy with Figure 4, the vertical, dashed lines in Figure 6a highlight the spatialpositions x O M and x O m at which the oxygen concentration at time T crosses, respectively, thethreshold values O M and O m ( i.e. S o (T , x O M ) = O M and S o (T , x O m ) = O m ). Hence, thewhite region corresponds to normoxic conditions, the pale-blue region corresponds to moderately-oxygenated environments and the blue region corresponds to hypoxic conditions. Moreover, thevertical, dashed lines in Figure 6b highlight the spatial positions x L m and x L M at which the lactateconcentration at time T crosses, respectively, the threshold values L m and L M ( i.e. S l (T , x L m ) = L m and S l (T , x L M ) = L M ). Hence, the white region corresponds to mildly-acidic conditions, thepale-green region corresponds to moderately-acidic conditions and the green region correspondsto highly-acidic conditions.As shown by the plots in Figure 6a and Figure 6b, the local mean levels of expression of thehypoxia- and acidity-resistant genes converge, as time passes, to the fittest ones ( i.e. µ (T , x )matches with ϕ o ( S o (T , x )) and µ (T , x ) matches with ϕ l ( S l (T , x )), the values of which vary withthe distance from the blood vessel depending on the local concentrations of oxygen and lactate.In more detail, the local mean level of expression of the hypoxia-resistant gene at time T isthe minimal one in normoxic conditions ( i.e. µ (T , x ) ≡ x ∈ [0 , x O M ]), the maximal one inhypoxic conditions ( i.e. µ (T , x ) ≡ x ∈ [ x O m , L]) and increases with the oxygen concentrationin moderately-oxygenated environments ( i.e. µ (T , x ) increases monotonically from 0 to 1 for x ∈ ( x O M , x O m )). Furthermore, the local mean level of expression of the acidity-resistant geneat time T is the minimal one in the mildly-acidic region ( i.e. µ (T , x ) ≡ x ∈ [0 , x L m ]), themaximal one in highly-acidic conditions ( i.e. µ (T , x ) ≡ x ∈ [ x L M , L]) and increases withthe lactate concentration in moderately-acidic environments ( i.e. µ (T , x ) increases monotonicallyfrom 0 to 1 for x ∈ ( x L m , x L M )).Finally, the plots in Figures 6c-e show that, at every position x ∈ [0 , L], the local phenotypicdistribution of tumour cells at time T ( i.e. the local population density function n (T , x, y )) isunimodal and attains its maximum at the fittest phenotypic state y = (cid:0) ϕ o ( S o (T , x )) , ϕ l ( S l (T , x ) (cid:1) .The numerical results of Figure 6 are complemented by the numerical results displayed in Figure 7,12hich summarise the time-evolution of the phenotypic distribution of tumour cells across the wholetissue region ( i.e. the function f ( t, y ) defined via (2)) and show that the maximum point of thedistribution departs from the point y = (0 ,
0) ( i.e. the point corresponding to the minimalexpression level of both the acidity-resistant gene and the hypoxia-resistant gene) – cf. the initialcondition n ( x, y ) defined via (19) – and moves toward the point y = (1 , i.e. the degree of malignancy of the tumour increases over time).Figure 6: a) Plots of the normalised level of expression of the hypoxia-resistant gene µ ( t, x ) at four successivetime instants (yellow, orange, red and light-green lines). The light-green line highlights µ (T , x ) and the burgundy,dashed line highlights the fittest level of expression of the hypoxia-resistant gene ϕ o ( S o (T , x )) defined via (9). Thevertical, dashed lines highlight the points x O M and x O m such that S o (T , x O M ) = O M and S o (T , x O m ) = O m .Hence, the white region corresponds to normoxic conditions, the pale-blue region corresponds to moderately-oxygenated environments and the blue region corresponds to hypoxic conditions. b) Plots of the normalised levelof expression of the acidity-resistant gene µ ( t, x ) at four successive time instants (yellow, orange, red and light-bluelines). The light-blue line highlights µ (T , x ) and the burgundy, dashed line highlights the fittest level of expressionof the acidity-resistant gene ϕ l ( S l (T , x )) defined via (11). The vertical, dashed lines highlight the points x L m and x L M such that S l (T , x L m ) = L m and S l (T , x L M ) = L M . Hence, the white region corresponds to mildly-acidicconditions, the pale-green region corresponds to moderately-acidic conditions and the green region corresponds tohighly-acidic conditions. c) - e) Plots of the local phenotypic distribution of tumour cells n (T , x, y ) at the points x = x a , x = x b and x = x c highlighted, respectively, by the circle, square and triangle markers shown in panels a)and b). As discussed in the previous section, the degree of malignancy of the tumour increases over timeand the majority of cancer cells are ultimately characterised by high levels of expression of thehypoxia-resistant gene and the acidity-resistant gene. Furthermore, the results we report on inthis section indicate that the dynamics of the levels of expression of the acidity-resistant gene andthe hypoxia-resistant gene across the whole tissue region ( i.e. the functions ν ( t ) and ν ( t ) definedvia (3)) are strongly affected by the values of the selection gradient related to oxygen η o and the13igure 7: Plots of the phenotypic distribution of tumour cells across the whole tissue region f ( t, y ) definedaccording to (2) at six successive time instants t < t < . . . < t . selection gradient related to lactate η l . In fact, as demonstrated by the plot in Figure 8, whichdisplays the curve φ ( t ) = (cid:0) ν ( t ) , ν ( t ) (cid:1) for different values of the ratio η o /η l , depending on thevalues of ratio between these parameters one has that resistance to hypoxia and resistance to acidityarise via alternative evolutionary pathways. Coherently with the results presented in the previoussection, φ ( t ) departs from the point (0 . , .
15) ( i.e. the point corresponding to the initial levelsof expression of the acidity-resistant gene and the hypoxia-resistant gene across the whole tissueregion) and, for all values of η o /η l considered, ultimately converges to a point corresponding to ahigh expression level of both the acidity-resistant gene and the hypoxia-resistant gene. However,larger values of the ratio η o /η l lead the level of expression of the hypoxia-resistant gene ν ( t ) toincrease faster than the level of expression of the acidity-resistant gene ν ( t ), while smaller valuesof η o /η l correlate with a faster increase of ν ( t ) and a slower increase of ν ( t ). Furthermore, forintermediate values of the ratio η o /η l we observe a simultaneous increase of the values of ν ( t ) and ν ( t ), whereas sufficiently large and sufficiently small values of η o /η l correlate with a decouplingbetween the increase of ν ( t ) and ν ( t ). In more detail, if the ratio η o /η l is sufficiently high, first ν ( t ) increases while ν ( t ) remains almost constant and then, when ν ( t ) is sufficiently high, ν ( t )starts increasing as well. On the other hand, in the case where η o /η l is sufficiently low, we havethat ν ( t ) increases first and then ν ( t ) starts increasing as soon as ν ( t ) becomes sufficiently high.These results communicate the biological notion that: the strength of the selective pressuresexerted by oxygen and lactate on tumour cells, which are quantified by the values of the selectiongradients η o and η l , may shape the emergence of hypoxic resistance and acidic resistance in tu-mours; the order in which such forms of resistance develop depend on the intensity of oxygen-drivenselection in relation to the intensity of lactate-driven selection. In this work, we have developed a mathematical modelling approach to investigate the influenceof hypoxia and acidity on the evolutionary dynamics of cancer cells in vascularised tumours.The results of numerical simulations of a calibrated version of the model based on real datarecapitulate the eco-evolutionary spatial dynamics of tumour cells and their adaptation to hypoxicand acidic microenvironments. In particular, the results obtained indicate that tumour cells char-acterised by lower levels of expression of hypoxia- and acidity-resistant genes are to be expected tocolonise well-oxygenated and mildly-acidic regions of vascularised tissues, whereas cells expressing14igure 8:
Plot of the curve φ ( t ) = (cid:0) ν ( t ) , ν ( t ) (cid:1) for different values of the ratio between the selection gradientrelated to oxygen η o and the selection gradient related to lactate η l . The functions ν ( t ) and ν ( t ) defined via (3)model the levels of expression of the acidity-resistant gene and the hypoxia-resistant gene across the whole tissueregion, respectively. a more aggressive phenotype characterised by higher levels of resistance to hypoxia and aciditywill ultimately populate tissue regions corresponding to hypoxic and acidic microenvironments.Such theoretical findings recapitulate histological data on ductal carcinoma in situ showing thatthe levels at which the acidity-resistant gene LAMP2 and the hypoxia-resistant gene GLUT-1are expressed by cancer cells increase moving from the walls to the centre of the milk duct ( i.e. moving from more oxygenated and less acidic regions to regions that are less oxygenated and moreacidic) [19, 26].Moreover, our theoretical findings reconcile the conclusions of [26], suggesting that tumourcells acquire first resistance to hypoxia and then resistance to acidity, and the conclusions of [57],supporting the idea that the two forms of resistance are acquired in reverse order, by showing thatthe order in which resistance to hypoxia and resistance to acidity arise depend on the ways inwhich oxygen and lactate act as environmental stressors in the evolutionary dynamics of tumourcells, which are known to vary between tissue types and between patients [46].We conclude with a brief overview of possible research perspectives. Along the lines of [43],the modelling framework presented here could be extended to incorporate additional details of cellmovement and mechanical interactions between cells [2, 8, 12], which would make it possible toinvestigate the interplay between phenotypic evolution of cancer cells and tumour growth.Moreover, building on [7], it would be interesting to generalise the model to study the effects offluctuations in the inflow rate of oxygen and glucose and the outflow rate of lactate in the evolutionof cancer cells. In fact, when in hypoxic conditions, cancer cells are known to produce and secreteproangiogenic factors which induce the formation of new blood vessels departing from existing ones.Such an angiogenic process results in the formation of a disordered tumour vasculature wherebythe rates at which oxygen and glucose enter the tumour and the rate at which lactate is flushedout through intra-tumour blood vessels fluctuate over time, which impacts on the evolutionarydynamics of cancer cells [21, 35, 49, 51].It would also be interesting to extend the model in order to investigate the role of pheno-typic transitions triggered by hypoxia and acidity – such as the epithelial-mesenchymal transitioninduced by hypoxic environmental conditions [53, 59, 65] and the acquisition of the metastaticphenotype promoted by acidic microenvironments [19, 20, 22] – in the phenotypic adaptation ofcancer cells and tumour growth.Finally, since resistance to hypoxia is known to correlate with resistance to chemotherapyand radiotherapy [18, 20, 38, 56, 61], building on [15, 42, 44], it would be relevant for anti-cancer therapy to address numerical optimal control for an extended version of the model that15akes into account the effect of chemotherapy and/or radiotherapy [6, 55], which could informthe development of optimised cancer treatment protocols that exploit evolutionary and ecologicalprinciples [1, 27, 36, 50]. Conflict of interest
The authors declare that they have no conflict of interest.
References [1] Ahmet Acar, Daniel Nichol, Javier Fernandez-Mateos, George D Cresswell, Iros Barozzi, SungPil Hong, Nicholas Trahearn, Inmaculada Spiteri, Mark Stubbs, Rosemary Burke and others,Exploiting evolutionary steering to induce collateral drug sensitivity in cancer, Nature Commu-nications, 11(1), 1–14 (2020)[2] Davide Ambrosi and Luigi Preziosi, On the closure of mass balance models for tumor growth,Mathematical Models and Methods in Applied Sciences, 12(5), 737–754 (2002)[3] Alexander RA Anderson, Alissa M Weaver, Peter T Cummings and Vito Quaranta, Tumormorphology and phenotypic evolution driven by selective pressure from the microenvironment,Cell, 127(5), 905–915 (2006)[4] Alexander RA Anderson and Vito Quaranta, Integrative mathematical oncology, Nature Re-views Cancer, 8(3), 227–234 (2008)[5] Alexander RA Anderson, Katarzyna A Rejniak, Philip Gerlee and Vito Quaranta, Microen-vironment driven invasion: a multiscale multimodel investigation, Journal of MathematicalBiology, 58(4-5), 579 (2009)[6] Lu´ıs Almeida, Patrizia Bagnerini, Giulia Fabrini, Barry D Hughes and Tommaso Lorenzi,Evolution of cancer cell populations under cytotoxic therapy and treatment optimisation: insightfrom a phenotype-structured model, ESAIM: Mathematical Modelling and Numerical Analysis,53(4), 1157–1190 (2019)[7] Aleksandra Ardaˇseva, Robert A Gatenby, Alexander RA Anderson, Helen M Byrne, Philip KMaini and Tommaso Lorenzi, A mathematical dissection of the adaptation of cell populationsto fluctuating oxygen levels, Bulletin of Mathematical Biology, 82(6), 81 (2020)[8] Sergey Astanin and Luigi Preziosi, Multiphase models of tumour growth. In: N. Bellomo, E.de Angelis (eds) Selected Topics in Cancer Modeling, 1–31. Birkh¨auser Boston (2008)[9] Sergey Astanin and Luigi Preziosi, Mathematical modelling of the Warburg effect in tumourcords, Journal of Theoretical Biology, 258(4), 578 – 590 (2009)[10] Henri Berestycki and Fran¸cois Hamel, Generalized transition waves and their properties,Communications on Pure and Applied Mathematics, 65(5), 592–648 (2012)[11] Henri Berestycki and Gr´egoire Nadin, Asymptotic spreading for general heterogeneous Fisher-KPP type equations, Memoirs of the American Mathematical Society, In press (2020)[12] Helen M Byrne and Luigi Preziosi, Modelling solid tumour growth using the theory of mix-tures, Mathematical Medicine and Biology: a journal of the IMA, 20(4), 341–366 (2003)[13] Helen M Byrne, Dissecting cancer through mathematics: from the cell to the animal model,Nature Reviews Cancer, 10(3), 221–230 (2010)1614] Joseph J Casciari, Stratis V Sotirchos and Robert M Sutherland, Variations in tumor cellgrowth rates and metabolism with oxygen concentration, glucose concentration, and extracel-lular pH, Journal of Cellular Physiology, 151(2), 386–394 (1992)[15] Mark AJ Chaplain, Tommaso Lorenzi and Chiara Villa, Evolutionary dynamics in vascu-larised tumours under chemotherapy, Vietnam Journal of Mathematics, In press (2020)[16] Mark AJ Chaplain, Multiscale Modelling of Cancer: Micro-, Meso- and Macro-scales ofGrowth and Spread. In: M. Bizzarri (eds) Approaching Complex Diseases, 149–168. Springer(2020)[17] Rebecca H Chisholm, Tommaso Lorenzi and Jean Clairambault, Cell population heterogene-ity and evolution towards drug resistance in cancer: biological and mathematical assessment,theoretical treatment optimisation, Biochimica et Biophysica Acta (BBA)-General Subjects,1860(11), 2627–2645 (2016)[18] Jean-Philippe Cosse and Carine Michiels, Tumour hypoxia affects the responsiveness of can-cer cells to chemotherapy and promotes cancer progression, Anti-Cancer Agents in MedicinalChemistry (Formerly Current Medicinal Chemistry-Anti-Cancer Agents, 8(7), 790–797 (2008)[19] Mehdi Damaghi, Narges K Tafreshi, Mark C Lloyd, Robert W Sprung, Veronica C Estrella,Jonathan W Wojtkowiak, David L Morse, John M Koomen, Marilyn M Bui, Robert A Gatenbyand Robert J Gillies, Chronic acidosis in the tumour microenvironment selects for overexpressionof LAMP2 in the plasma membrane, Nature Communications, 6, 8752 (2015)[20] Katie DeClerck and Randolph C Elble, The role of hypoxia and acidosis in promoting metas-tasis and resistance to chemotherapy, Frontiers in Bioscience, 15, 213–225 (2010)[21] Mark W Dewhirst, Relationships between cycling hypoxia, HIF-1, angiogenesis and oxidativestress, Radiation Research, 172(6), 653–665 (2009)[22] Stefano Fais, Giulietta Venturi and Bob Gatenby, Microenvironmental acidosis in carcinogen-esis and metastases: new strategies in prevention and therapy, Cancer and Metastasis Reviews,33(4), 1095–1108 (2014)[23] Jill Gallaher and Alexander RA Anderson, Evolution of intratumoral phenotypic heterogene-ity: the role of trait inheritance, Interface Focus, 3(4), 20130016 (2016)[24] Jill A Gallaher, Joel Brown and Alexander RA Anderson, The impact of proliferation-migration tradeoffs on phenotypic evolution in cancer, Scientific Reports, 9, 2425 (2019)[25] Robert A Gatenby and Edward T Gawlinski, A reaction-diffusion model of cancer invasion,Cancer Research, 56(24), 5745–5753 (1996)[26] Robert A Gatenby, Kieran Smallbone, Philip K Maini, Fabrice Rose, James G Averill, Ray-mond B Nagle, Liam J Worrall and Robert J Gillies, Cellular adaptations to hypoxia andacidosis during somatic evolution of breast cancer, British Journal of Cancer, 97(5), 646–653(2007)[27] Robert A Gatenby, Ariosto S Silva, Robert J Gillies and B Roy Frieden, Adaptive therapy,Cancer Research, 69(11), 4894–4903 (2009)[28] Laura Gay, Ann-Marie Baker and Trevor A Graham, Tumour cell heterogeneity, Faculty of1000 Ltd, 5, (2016)[29] Robert J Gillies, Daniel Verduzco and Robert A Gatenby, Evolutionary dynamics of car-cinogenesis and why targeted therapy does not work, Nature Reviews Cancer, 12(7), 487–493(2012) 1730] Sara J Hamis, Gibin G Powathil, Can we crack cancer?. In: F. Matth¨aus, S. Matth¨aus, S.Harris, T. Hillen (eds) The Art of Theoretical Biology, 50–51. Springer (2020)[31] Michael Hockel and Peter Vaupel, Tumor hypoxia: definitions and current clinical, biologic,and molecular aspects, Journal of the National Cancer Institute, 93(4), 266–276 (2001)[32] Arig Ibrahim-Hashim, Mark Robertson-Tessi, Pedro M Enriquez-Navas, Mehdi Damaghi,Yoganand Balagurunathan, Jonathan W Wojtkowiak, Shonagh Russell, Kam Yoonseok, MarkC Lloyd, Marilyn M Bui and others, Defining cancer subpopulations by adaptive strategiesrather than molecular properties provides novel insights into intratumoral evolution, CancerResearch, 77(9), 2242–2254 (2017)[33] Jung-whan Kim and Chi V Dang, Cancer’s molecular sweet tooth and the Warburg effect,Cancer Research, 66(18), 8927–8930 (2006)[34] Yangjin Kim, Hyunji Kang, Gibin Powathil, Hyeongi Kim, Dumitru Trucu, Wanho Lee,Sean E Lawler and Mark AJ Chaplain, Role of extracellular matrix and microenvironmentin regulation of tumor growth and LAR-mediated invasion in glioblastoma, PloS One, 13(10)(2018)[35] Hiroyuki Kimura, Rod D Braun, Edgardo T Ong, Richard Hsu, Timothy W Secomb,Demetrios Papahadjopoulos, Keelung Hong and Mark W Dewhirst, Fluctuations in red cell fluxin tumor microvessels can lead to transient hypoxia and reoxygenation in tumor parenchyma,Cancer Research, 56(23), 5522–5528 (1996)[36] Kirill S Korolev, Joao B Xavier and Jeff Gore, Turning ecology and evolution against cancer,Nature Reviews Cancer, 14(5), 371–380 (2014)[37] Randall J LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equa-tions: steady-state and time-dependent problems, Society for Industrial and Applied Mathe-matics (SIAM), Philadelphia (2007)[38] Thomas D Lewin, Philip K Maini, Eduardo G Moros, Heiko Enderling and Helen M Byrne,The evolution of tumour composition during fractionated radiotherapy: implications for out-come, Bulletin of Mathematical Biology, 80(5), 1207–1235 (2018)[39] Mark C Lloyd, Jessica J Cunningham, Marilyn M Bui, Robert J Gillies, Joel S Brown, Joel Sand Robert A Gatenby, Darwinian dynamics of intratumoral heterogeneity: not solely randommutations but also variable environmental selection forces, Cancer Research, 76(11), 3136–3144(2016)[40] Lawrence A Loeb, A mutator phenotype in cancer, Cancer Research, 61(8), 3230–3239 (2001)[41] Tommaso Lorenzi, Rebecca H Chisholm and Jean Clairambault, Tracking the evolution ofcancer cell populations through the mathematical lens of phenotype-structured equations, Bi-ology Direct, 11(1), 43 (2016)[42] Tommaso Lorenzi, Chandrasekhar Venkataraman, Alexander Lorz and Mark AJ Chaplain,The role of spatial variations of abiotic factors in mediating intratumour phenotypic hetero-geneity, Journal of Theoretical Biology, 451(14), 101 –110 (2018)[43] Tommaso Lorenzi, Benoˆıt Perthame and Xinran Ruan, Invasion fronts and adaptivedynamics in a model for the growth of cell populations with heterogeneous mobility,https://arxiv.org/abs/2007.13084 (2020)[44] Alexander Lorz, Tommaso Lorenzi, Jean Clairambault, Alexandre Escargueil and BenoˆıtPerthame, Modeling the effects of space structure and combination therapies on phenotypicheterogeneity and drug resistance in solid tumors, Bulletin of Mathematical Biology, 77(1),1–22 (2015) 1845] Alexander Lorz, Tommaso Lorenzi, Michael Hochberg, Jean Clairambault and BenoˆıtPerthame, Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies, ESAIM: Mathematical Modelling and Numerical Analysis, 47(2), 377–399(2013)[46] Carlo C Maley, Athena Aktipis, Trevor A Graham, Andrea Sottoriva, Amy M Boddy,Michalina Janiszewska, Ariosto S Silva, Marco Gerlinger, Yinyin Yuan, Kenneth J Pienta etal., Classifying the evolutionary and ecological features of neoplasms, Nature Reviews Cancer,17(10), 605–619 (2017)[47] Ubaldo E Martinez-Outschoorn, Maria Peiris-Pag´es, Richard G Pestell, Federica Sotgia,Michael P Lisanti, Cancer metabolism: a therapeutic perspective, Nature Reviews ClinicalOncology, 14(1), 11–31 (2016)[48] Andriy Marusyk, Vanessa Almendro and Kornelia Polyak, Intra-tumour heterogeneity: alooking glass for cancer?, Nature Reviews Cancer, 12(5), 323–334 (2012)[49] Shingo Matsumoto, Hironobu Yasui, James B Mitchell and Murali C Krishna, Imaging cyclingtumor hypoxia, Cancer Research, 70(24), 10019–10023 (2010)[50] Lauren MF Merlo, John W Pepper, Brian J Reid, and Carlo C Maley, Cancer as an evolu-tionary and ecological process, Nature Reviews Cancer, 6(12), 924–935 (2006)[51] Carine Michiels, C´eline Tellier and Olivier Feron, Cycling hypoxia: A key feature of the tumormicroenvironment, Biochimica et Biophysica Acta (BBA)-Reviews on Cancer, 1866(1), 76–86(2016)[52] Franziska Michor and Kornelia Polyak, The origins and implications of intratumor hetero-geneity, Cancer Prevention Research, 3(11), 1361–1364 (2010)[53] Ashish Misra, Chhiti Pandey, Siu Kwan Sze and Thirumaran Thanabalu, Hypoxia activatedEGFR signaling induces epithelial to mesenchymal transition (EMT), PloS One, 7(11), e49766(2012)[54] Hamid R Molavian, Mohammad Kohandel, Michael Milosevic and Sivabal Sivaloganathan,Fingerprint of cell metabolism in the experimentally observed interstitial pH and pO2 in solidtmors, Cancer Research, 69(23), 9141–9147 (2009)[55] Camille Pouchol, Jean Clairambault, Alexander Lorz and Emmanuel Tr´elat, Asymptotic anal-ysis and optimal control of an integro-differential system modelling healthy and cancer cellsexposed to chemotherapy, Journal de Math´ematiques Pures et Appliqu´ees, 116, 268–308 (2018)[56] Sotiris Prokopiou, Eduardo G Moros, Jan Poleszczuk, Jimmy Caudell, Javier F Torres-Roca,Kujtim Latifi, Jae K Lee, Robert Myerson, Louis B Harrison and Heiko Enderling, A prolifera-tion saturation index to predict radiation response and personalize radiotherapy fractionation,Radiation Oncology, 10(1), 1–8 (2015)[57] Mark Robertson-Tessi, Robert J Gillies, Robert A Gatenby and Alexander RA Anderson,Impact of metabolic heterogeneity on tumor Growth, invasion, and treatment outcomes, CancerResearch, 75(8), 1567–1579 (2015)[58] Maximilian AR Strobl, Andrew L Krause, Mehdi Damaghi, Robert J Gillies, Alexander RAAnderson and Philip K Maini, Mix and match: phenotypic coexistence as a key facilitator ofcancer invasion, Bulletin of Mathematical Biology, 82(1), 1–26 (2020)[59] Shing Y Tam, Vincent WC Wu and Helen KW Law, Hypoxia-induced epithelial-mesenchymaltransition in cancers: HIF-1 α and Beyond, Frontiers in Oncology, 10, 486 (2020)1960] Xiaohu Tang, Joseph E Lucas, Julia Ling-Yu Chen, Gregory LaMonte, Jianli Wu, MichaelChangsheng Wang, Constantinos Koumenis and Jen-Tsan Chi, Functional interaction betweenresponses to lactic acidosis and hypoxia regulates genomic transcriptional outputs, Cancer Re-search, 72(2), 491–502 (2012)[61] Beverly A Teicher, Hypoxia and drug resistance, Cancer and Metastasis Reviews, 13(2), 139–168 (1994)[62] Catherine Vander Linden and Cyril Corbet, Reconciling environment-mediated metabolicheterogeneity with the oncogene-driven cancer paradigm in precision oncology, Seminars in Celland Developmental Biology, 98, 202–210 (2019)[63] Peter Vaupel, Oliver Thews and Michael Hoeckel, Treatment resistance of solid tumors: roleof hypoxia and anemia, Medical Oncology, 18, 243 (2001)[64] Chiara Villa, Mark AJ Chaplain and Tommaso Lorenzi, Modelling the emergence of pheno-typic heterogeneity in vascularised tumours, https://arxiv.org/abs/1910.08566 (2019)[65] Wenjing Zhang, Xinpeng Shi, Ying Peng, Meiyan Wu, Pei Zhang, Ruyi Xie, Yao Wu, QingqingYan, Side Liu Jide Wang, HIF-1 αα