A Mathematical Study of the Influence of Hypoxia and Acidity on the Evolutionary Dynamics of Cancer

Hypoxia and acidity act as environmental stressors promoting selection for cancer cells with a more aggressive phenotype. As a result, a deeper theoretical understanding of the spatio-temporal processes that drive the adaptation of tumour cells to hypoxic and acidic microenvironments may open up new avenues of research in oncology and cancer treatment. We present a mathematical model to study the influence of hypoxia and acidity on the evolutionary dynamics of cancer cells in vascularised tumours. The model is formulated as a system of partial integro-differential equations that describe the phenotypic evolution of cancer cells in response to dynamic variations in the spatial distribution of three abiotic factors that are key players in tumour metabolism: oxygen, glucose and lactate. The results of numerical simulations of a calibrated version of the model based on real data recapitulate the eco-evolutionary spatial dynamics of tumour cells and their adaptation to hypoxic and acidic microenvironments. Moreover, such results demonstrate how nonlinear interactions between tumour cells and abiotic factors can lead to the formation of environmental gradients which select for cells with phenotypic characteristics that vary with distance from intra-tumour blood vessels, thus promoting the emergence of intra-tumour phenotypic heterogeneity. Finally, our theoretical findings reconcile the conclusions of earlier studies by showing that the order in which resistance to hypoxia and resistance to acidity arise in tumours depend on the ways in which oxygen and lactate act as environmental stressors in the evolutionary dynamics of cancer cells.


Introduction
Cancer is a dynamic disease, the characteristics of which are constantly evolving. This is reflected in the fact that the genotypic and phenotypic properties of cancer cells may change across space and time within the same tumour, and the dynamics of tumours with the same histological features are still likely to vary across patients. Moreover, since the same cancer clones may arise through different evolutionary pathways, the fact that two tumours have a similar clonal composition at a 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 emergence and development of intra-and inter-tumour heterogeneity, which are major obstacles to cancer eradication [29,48].
Clinical evidence suggests that cancer cells and the tumour microenvironment mutually shape each other [24]. This supports the idea that tumours can be seen as evolving ecosystems where cancer cells with different phenotypic characteristics proliferate, die, undergo genotypic and phenotypic changes, and compete for space and resources under the selective pressure exerted by the various 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-evolutionary process in which spatial variability of the concentration of abiotic factors (i.e. substrates and metabolites) across the tumour supports the formation of distinct ecological niches whereby cells with different phenotypic characteristics may be selected [14,26,31].
In normal tissues, cells produce the energy required to sustain their proliferation via oxidative phosphorylation (i.e. they rely on oxygen as their primary source of energy) and turn to glycolysis only when oxygen is scarce. In tumours, the presence of hypoxic regions (i.e. regions where the oxygen levels are below the physiological ones) induces cells to transiently switch to a glycolytic metabolic phenotype (i.e. to rely on glucose as their primary source of energy) [63]. Cancer cells eventually acquire such a glycolytic phenotype and express it also in aerobic conditions, leading to the so-called Warburg effect [33]. The interplay between the high glycolytic rate of cancer cells and low perfusion in tumours brings about accumulation of lactate (i.e. a waste product of glycolysis), which causes acidity levels in the tumour microenvironment to rise (i.e. the pH level drops) [60].
Since hypoxia and acidity act as environmental stressors promoting selection for cancer cells with a more aggressive phenotype [57,60], an in-depth theoretical understanding of the spatiotemporal processes that drive the adaptation of tumour cells to hypoxic and acidic microenvironments 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 enable extrapolation beyond scenarios which can be investigated through experiments and may reveal emergent 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 explore how nonlinear interactions between cancer cells and abiotic components of the tumour microenvironment may shape tumour growth. The Gatenby-Gawlinski model has recently been extended in [58], in order to take into account the presence of cells with different phenotypic characteristics within the tumour. Hybrid cellular automaton models have been employed to study the impact of hypoxia and acidity on tumour growth and invasion [3,5,26,34,57]. A mechanical model of tumour growth whereby cells are allowed to switch between aerobic and anaerobic metabolism was presented in [9]. Integro-differential equations and partial integro-differential have been used in [7,15,42,44,64] to investigate the ecological role of hypoxia in the development of intra-tumour phenotypic heterogeneity.
In this paper, we complement these earlier studies by presenting a mathematical model to study the influence of hypoxia and acidity on the evolutionary dynamics of cancer cells in vascularised tumours. The model comprises a system of partial integro-differential equations that describe the phenotypic evolution of cancer cells in response to dynamic variations in the spatial distribution of 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 equations and the underlying modelling assumptions. In Section 3, we summarise the main results of numerical simulations of the model and discuss their biological implications. Section 4 concludes the paper and provides a brief overview of possible research perspectives.

Model description
We consider a one-dimensional region of vascularised tissue of length L > 0. We describe the spatial position of every tumour cell in the tissue region by a scalar variable x ∈ [0, L] and we assume a blood vessel to be present at x = 0 (cf. the schematic in Figure 1a). Moreover, building upon the modelling framework developed in [15,42,44,64], we model the phenotypic state of every cell by a vector y = (y 1 , y 2 ) ∈ [0, 1] 2 (cf. the schematics in Figure 1b). Here, y 1 ∈ [0, 1] represents the normalised level of expression of an acidity-resistant gene (e.g. the LAMP2 gene), while y 2 ∈ [0, 1] represents the normalised level of expression of a hypoxia-resistant gene (e.g. the GLUT-1 gene) [19,26].
We describe the phenotypic distribution of tumour cells at position x and time t ∈ [0, T], with T > 0, by means of the local population density function n(t, x, y) (i.e. the local phenotypic distribution of tumour cells). We define the cell density ρ(t, x), the local mean level of expression of the acidity-resistant gene µ 1 (t, x) and the local mean level of expression of the hypoxia-resistant gene µ 2 (t, x) as for i = 1, 2. Moreover, we define the phenotypic distribution of tumour cells across the whole tissue region f (t, y) as the mean value of n(t, x, y) on the interval [0, L], i.e.
Similarly, we define the levels of expression of the acidity-resistant gene and the hypoxia-resistant gene across the whole tissue region as the mean values of µ 1 (t, x) and µ 2 (t, x) on the interval [0, L], respectively, i.e.
The local concentrations of oxygen, glucose and lactate at position x and time t are denoted by 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 length L. A blood vessel is assumed to be present at x = 0. b) Schematics illustrating the relationships between the values of the variables y 2 and y 1 modelling the phenotypic state of tumour cells and the levels of resistance to hypoxia and acidosis.

Dynamics of tumour cells
We describe the dynamics of tumour cells through the following balance equation for the local population density function n(t, x, y) with (t, x, y) ∈ (0, T] × (0, L) × (0, 1) 2 , 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 leave the tissue region) and zero-flux boundary conditions on the boundary of the square [0, 1] 2 (i.e. we assume that cells cannot have normalised levels of gene expression smaller than 0 or larger than 1). The first term on the right-hand side of the partial integro-differential equation (4) models the effect of undirected, random movement, which is described through Fick's first law of diffusion with diffusivity β n > 0, while the second term models the effect of heritable, spontaneous phenotypic changes, which occur at rate θ > 0. The function R(S o , S g , S l , ρ, y) represents the fitness of tumour cells in the phenotypic state y at position x and time t under the environmental conditions given by 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 following definition proliferation and death due to oxygen-driven selection − D(S l , y 1 ) death due to lactate-driven selection − d(ρ). death due to competition for space (5) Here, the function P (S o , S g , y 2 ) is the rate at which cells with level of expression y 2 of the hypoxiaresistant gene proliferate via oxidative phosphorylation and glycolysis, and die due to oxygendriven selection (i.e. P (S o , S g , y 2 ) is a net proliferation rate). The function D(S l , y 1 ) is the rate at which cells with level of expression y 1 of the acidity-resistant gene die due to lactate-driven selection. The function d(ρ) models the rate of cell death due to competition for space associated with saturation of the cell density.

Modelling oxygen-driven selection
Based on the theoretical results and experimental data presented in [36,63], we focus on a scenario corresponding to the following biological assumptions. Assumption 2. Cells proliferate at a rate that depends on the concentrations of oxygen and glucose. Moreover, the trade-off between increase in cell death associated with sensitivity to hypoxia and decrease in cell proliferation associated with acquisition of resistance to hypoxia results in the existence of a level of expression of the hypoxia-resistant gene which is the fittest in that: a lower level of gene expression would correlate with a lower resistance to hypoxia, and thus a higher death rate; a higher level of gene expression would correlate with a larger fitness cost, and thus a lower proliferation rate. Cells with levels of gene expression that are closer to the fittest one are more likely to survive than the others. Hence, the farther the gene expression level of a cell is from the fittest one, the more likely is that the cell will die due to a form of oxygen-driven selection.  Under Assumptions 1 and 2, we define the net proliferation rate P (S o , S g , y 2 ) as In (6), the function p o (S o ) models the rate of cell proliferation via oxidative phosphorylation, while the function p g (S o , S g ) models the rate of cell proliferation via glycolysis. Furthermore, the third term in the definition given by (6)  Remark 1. In (6), the distance between [7,41], it leads to a smoother fitness function which is closer to the approximate fitness landscapes which can be inferred from experimental data through regression techniques.

Under Assumptions 3 and 4, we use the definitions of the functions
with and In (7), the parameters γ o > 0 and γ g > 0 model the maximum rates of cell proliferation via oxidative phosphorylation and glycolysis, respectively. The parameters α o > 0 and α g > 0 are the Michaelis-Menten constants of oxygen and glucose. The weight function w(S o ) defined via (8) ensures that Assumption 3 is satisfied, while definition (9) of ϕ o (S o ) is such that Assumption 4 is satisfied (cf. the plot in Figure 2a).

Modelling lactate-driven selection
Based on theoretical results and experimental data presented in [57], we focus on a scenario corresponding to the following biological assumptions.
Assumption 5. There exist two threshold levels of the lactate concentration L M > L m > 0 such that 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, the trade-off between increase in cell death associated with sensitivity to acidity and decrease in cell proliferation associated with acquisition of resistance to acidity results in the existence of a level of expression of the acidity-resistant gene which is the fittest in that: a lower level of gene expression would correlate with a lower resistance to acidity, and thus a higher death rate; a higher level of gene expression would correlate with a larger fitness cost, and thus a lower proliferation rate. Cells with levels of gene expression that are closer to the fittest one are more likely to survive than the others. Hence, the farther the gene expression level of a cell is from the fittest one, the more likely is 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 associated with the variable y 1 ) may vary with the lactate concentration. In particular: in mildly-acidic environments (i.e. when S l ≤ L m ), the fittest level of gene expression is the minimal one (i.e. y 1 = 0); in highly-acidic environments (i.e. when S l ≥ L M ) the fittest level of gene expression is the maximal one (i.e. y 1 = 1); in moderately-acidic environments (i.e. when L m < S l < L M ), the fittest level of gene expression is a monotonically increasing function of the lactate concentration (i.e. it increases from y 1 = 0 to y 1 = 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 1 ) as In (10), the parameter η l > 0 is a selection gradient that quantifies the intensity of lactate-driven selection and the function ϕ l (S l ) is the fittest level of expression of the acidity-resistant gene under the environmental conditions given by the lactate concentration S l . Considerations analogous to those made in Remark 1 on the term (6) apply to the term y 1 − ϕ l (S l ) 2 in (10).
Finally, we use the definition of the function ϕ l (S l ) given hereafter (cf. the plot in Figure 2b), so that Assumption 7 is satisfied:

Modelling competition for space
We define the rate of cell death due to competition for space associated with saturation of the cell density as d(ρ) := κ ρ.
Here, the proportionality constant κ > 0 is related to the local carrying capacity of the tumour, which may vary depending on the tumour type.

Dynamics of abiotic factors
Oxygen and glucose are consumed by tumour cells, while lactate is produced by tumour cells as a waste product of glycolysis. Moreover, oxygen, glucose and lactate diffuse in space and decay over time. Hence, their dynamics are governed by the following balance equations for the functions consumption by tumour cells consumption by tumour cells (14) and production by tumour cells (15) with (t, x) ∈ (0, T] × (0, L), subject to suitable boundary conditions (see considerations below) and initial conditions. In (13)-(15), the parameters β o > 0, β g > 0 and β l > 0 are the diffusion coefficients of oxygen, glucose and lactate, respectively, while the parameters λ o > 0, λ g > 0 and λ l > 0 are the natural decay rates of the three abiotic factors. The third term on the right-hand side of (13) is the consumption rate of oxygen by tumour cells, which is proportional to the product between the cell density ρ and the rate of cell proliferation via oxidative phosphorylation p o (S o ), which is defined via (7). The parameter ζ o > 0 is a conversion factor linked to oxygen consumption by the cells. Analogous considerations hold for the third term on the right-hand side of (14), which models the consumption rate of glucose by tumour cells. Furthermore, the third term on the right-hand side of (15) is the production rate of lactate by tumour cells, which is assumed to be proportional to the product between the cell density ρ 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 > 0 linked to lactate production by the cells.
We assume that oxygen and glucose enter the spatial domain through the blood vessel only, while lactate is flushed out through the blood vessel only. Hence, focussing on the case where the inflow rate of oxygen and glucose and the outflow rate of lactate are constant, we complement (13)- (15) with the following boundary conditions at x = 0 where S o > 0 and S g > 0 are related to the average physiological levels of oxygen and glucose in proximity to blood vessels, while S l > 0 is a small parameter of value close to zero. In particular, we have S o > O M and S l < L m . Moreover, we assume that far from the blood vessel the concentrations 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 regions distant from blood vessels. In particular, we have S o < O m . Furthermore, we model abnormal accumulation of lactate, which is expected to occur far from blood vessels, imposing zero-flux boundary conditions. Therefore, we complement (13)-(15) with the following boundary conditions

Main results
In this section, we present the results of numerical simulations of the mathematical model defined by (4) coupled with (13)-(15) and we discuss their biological relevance. First, we describe the setup of numerical simulations (see Section 3.1). Next, we present a sample of numerical solutions that 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 tumour cells and the emergence of phenotypic heterogeneity (see Section 3.3). Finally, we present the results of numerical simulations that reveal the existence of alternative evolutionary pathways that may lead to the development of resistance to hypoxia and acidity in vascularised tumours (see Section 3.4).

Set-up of numerical simulations
Numerical simulations are carried out assuming L = 400 µm, which is chosen coherently with experimental data reported in [54], and t ∈ [0, T], where the final time T > 0 is such that the solutions of the model equations are at numerical equilibrium for t = T.
Initial conditions. We consider (13), (14) and (15) subject, respectively, to the following initial conditions Here, the functions S 0 o (x) and S 0 g (x) (see Figure 3) are defined in such a way as to match the experimental 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) correspond to a situation in which the initial distributions of oxygen and glucose match with experimental equilibrium distributions of such abiotic factors and lactate is present at a uniform level which is below the threshold level L m , that is, the level below which the environment surrounding the cells is 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 0 (x, y) : which corresponds to a biological scenario in which at the initial time t = 0 most tumour cells are concentrated near the blood vessel and are characterised by the minimal expression level of both the hypoxia-resistant gene and the acidity-resistant gene.
The values of S o and S g correspond to the average physiological levels of oxygen and glucose in proximity to blood vessels reported in [54]. Moreover, the values of S o and S g correspond to the 0.1% of S o and the 1% of S g , respectively. This is because, based on experimental data reported in [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. Parameter values. Unless otherwise explicitly stated, we use the values of the model parameters listed in Table 1, which are chosen to be consistent with the existing literature, except for the values of the parameters η o , η l , λ l and ζ l that are model specific in that we could not find them in the literature and are defined on the basis of the following considerations. The value of the conversion factor for lactate production ζ l is chosen to be the same as the value of the conversion factor 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)) is similar 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 parameters is performed in Section 3.4.
Numerical methods. Numerical solutions are constructed using a uniform discretisation of the interval [0, L] as the computational domain of the independent variable x and a uniform discretisation of the square [0, 1] 2 as the computational domain of the independent variable y. We also discretise the interval [0, T] with a uniform step. The method for constructing numerical solutions is based on an explicit finite difference scheme in which a three-point and a five-point stencils are used to approximate the diffusion terms in x and y, respectively, and an implicitexplicit finite difference scheme is used for the reaction terms [37,45]. All numerical computations are performed in Matlab.

Dynamics of the cell density and the concentrations of abiotic factors
The dynamics of the density of tumour cells and the concentrations of abiotic factors are illustrated by the plots in Figure 4. In summary, the cell density and the concentration of lactate at successive times (i.e. the graphs of ρ(t, x) for four different values of t and S l (t, x) for three different values of t) are displayed in Figure 4a and Figure 4c, respectively, while the concentrations of oxygen and 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 death leads tumour cells, which are initially concentrated near the blood vessel at x = 0 (cf. the initial condition n 0 (x, y) defined via (19)), to spread into the surrounding tissue (i.e. ρ(t, x) behaves like an invasion front). At every position, the cell density grows until a saturation value, which varies with the distance from the blood vessel, is reached (i.e. ρ(t, x) appears to converge to a form of generalised transition wave [10,11]). Moreover, the curves in Figure 4b show that the reactiondiffusion dynamics of oxygen and glucose, along with the inflow through the blood vessel and the consumption by tumour cells, lead the concentrations of such abiotic factors to converge to some stable 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 functions of x). Furthermore, the curves in Figure 4c show that the interplay of production by tumour cells, outflow through the blood vessel and reaction-diffusion dynamics leads the concentration of lactate to grow at every position until a saturation value, which varies with the distance from the blood 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 initial distributions S 0 o (x) and S 0 g (x) displayed in Figure 3. This is to be expected. In fact, since S 0 o (x) and S 0 g (x) are defined in such a way as to match experimental equilibrium distributions of oxygen and glucose, 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 quickly numerical 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 not depend on the choice of the initial conditions.  (8)). Coherently with experimental evidence indicating that using glycolysis to produce energy required for cell proliferation is less efficient than employing oxidative phosphorylation, the biologically informed parameter values considered here are such that the cell proliferation rate in normoxic conditions is higher than the cell proliferation rate in moderately-oxygenated environments, which in turn is higher than the cell proliferation rate in hypoxic conditions (cf. the plot in Figure 5, blue curve). As a result, the saturation value of the cell density decreases with the distance from the blood vessel (i.e. ρ(T, x) is a monotonically decreasing function of x). Furthermore, the production rate of lactate by tumour cells is proportional to the rate of proliferation via glycolysis and, therefore, it increases with the distance from the blood vessel at x = 0 (cf. the plot in Figure 5, red curve). As a result, the saturation value of the lactate concentration increases with the distance from the blood vessel and, in agreement with the lactate distributions reported in [54], S l (T, x) is a monotonically increasing function of x. 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.

Evolutionary dynamics of tumour cells and emergence of phenotypic heterogeneity
As discussed in the previous section, reaction-diffusion dynamics of abiotic factors and mutual interactions between abiotic factors and tumour cells lead to the emergence of spatial variations in the concentrations of oxygen and lactate -i.e. the oxygen concentration S o (T, x) is a monotonically decreasing function of x while the lactate concentration S l (T, x) is a monotonically increasing function of x (cf. plots in Figure 4b and Figure 4c). Spatial variability of oxygen and lactate con- Hence, the white region corresponds to normoxic conditions, the pale-blue region corresponds to moderatelyoxygenated environments and the blue region corresponds to hypoxic conditions.
centrations can lead to the formation of environmental gradients resulting in the selection for cells with phenotypic characteristics that vary with distance from the blood vessel, thus promoting the emergence of intra-tumour phenotypic heterogeneity. This is demonstrated by the numerical results in Figure 6. The dashed lines in Figure 6a and Figure 6b highlight the fittest levels of expression of the hypoxia-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 mean levels 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 µ 2 (t, x) and µ 1 (t, x) for four different values of t). In analogy with Figure 4, the vertical, dashed lines in Figure 6a  . Hence, the white region corresponds to normoxic conditions, the pale-blue region corresponds to moderatelyoxygenated environments and the blue region corresponds to hypoxic conditions. Moreover, the vertical, dashed lines in Figure 6b highlight the spatial positions x Lm and x L M at which the lactate concentration at time T crosses, respectively, the threshold values L m and L M (i.e. S l (T, x Lm ) = L m and S l (T, x L M ) = L M ). Hence, the white region corresponds to mildly-acidic conditions, the pale-green region corresponds to moderately-acidic conditions and the green region corresponds to highly-acidic conditions.
As shown by the plots in Figure 6a and Figure 6b, the local mean levels of expression of the hypoxia-and acidity-resistant genes converge, as time passes, to the fittest ones (i.e. µ 2 (T, x) matches with ϕ o (S o (T, x)) and µ 1 (T, x) matches with ϕ l (S l (T, x)), the values of which vary with the 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 is the minimal one in normoxic conditions (i.e. µ 2 (T, x) ≡ 0 for x ∈ [0, x O M ]), the maximal one in hypoxic conditions (i.e. µ 2 (T, x) ≡ 1 for x ∈ [x Om , L]) and increases with the oxygen concentration in moderately-oxygenated environments (i.e. µ 2 (T, x) increases monotonically from 0 to 1 for x ∈ (x O M , x Om )). Furthermore, the local mean level of expression of the acidity-resistant gene at time T is the minimal one in the mildly-acidic region (i.e. µ 1 (T, x) ≡ 0 for x ∈ [0, x Lm ]), the maximal one in highly-acidic conditions (i.e. µ 1 (T, x) ≡ 1 for x ∈ [x L M , L]) and increases with the lactate concentration in moderately-acidic environments (i.e. µ 1 (T, x) increases monotonically from 0 to 1 for x ∈ (x Lm , x L M )).
Finally, the plots in Figures 6c-e show that, at every position x ∈ [0, L], the local phenotypic distribution of tumour cells at time T (i.e. the local population density function n(T, x, y)) is unimodal and attains its maximum at the fittest phenotypic state y = ϕ o (S o (T, x)), ϕ l (S l (T, x) . The numerical results of Figure 6 are complemented by the numerical results displayed in Figure 7, which summarise the time-evolution of the phenotypic distribution of tumour cells across the whole tissue region (i.e. the function f (t, y) defined via (2)) and show that the maximum point of the distribution departs from the point y = (0, 0) (i.e. the point corresponding to the minimal expression level of both the acidity-resistant gene and the hypoxia-resistant gene) -cf. the initial condition n 0 (x, y) defined via (19) -and moves toward the point y = (1, 1), which corresponds to the maximal expression level of both the hypoxia-resistant gene and the acidity-resistant gene (i.e. the degree of malignancy of the tumour increases over time). Hence, the white region corresponds to normoxic conditions, the pale-blue region corresponds to moderatelyoxygenated environments and the blue region corresponds to hypoxic conditions. b) Plots of the normalised level of expression of the acidity-resistant gene µ 1 (t, x) at four successive time instants (yellow, orange, red and light-blue lines). The light-blue line highlights µ 1 (T, x) and the burgundy, dashed line highlights the fittest level of expression of the acidity-resistant gene ϕ l (S l (T, x)) defined via (11). The vertical, dashed lines highlight the points x Lm and x L M such that S l (T, x Lm ) = Lm and S l (T, x L M ) = L M . Hence, the white region corresponds to mildly-acidic conditions, the pale-green region corresponds to moderately-acidic conditions and the green region corresponds to highly-acidic conditions. c) -e) Plots of the local phenotypic distribution of tumour cells n(T, x, y) at the points x = xa, x = x b and x = xc highlighted, respectively, by the circle, square and triangle markers shown in panels a) and b).

Alternative evolutionary pathways leading to the development of resistance to hypoxia and acidity
As discussed in the previous section, the degree of malignancy of the tumour increases over time and the majority of cancer cells are ultimately characterised by high levels of expression of the hypoxia-resistant gene and the acidity-resistant gene. Furthermore, the results we report on in this section indicate that the dynamics of the levels of expression of the acidity-resistant gene and the hypoxia-resistant gene across the whole tissue region (i.e. the functions ν 1 (t) and ν 2 (t) defined via (3)) are strongly affected by the values of the selection gradient related to oxygen η o and the selection gradient related to lactate η l . In fact, as demonstrated by the plot in Figure 8, which displays the curve φ(t) = ν 1 (t), ν 2 (t) for different values of the ratio η o /η l , depending on the values of ratio between these parameters one has that resistance to hypoxia and resistance to acidity arise via alternative evolutionary pathways. Coherently with the results presented in the previous section, φ(t) departs from the point (0.15, 0.15) (i.e. the point corresponding to the initial levels of expression of the acidity-resistant gene and the hypoxia-resistant gene across the whole tissue region) and, for all values of η o /η l considered, ultimately converges to a point corresponding to a high 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 ν 2 (t) to increase faster than the level of expression of the acidity-resistant gene ν 1 (t), while smaller values of η o /η l correlate with a faster increase of ν 1 (t) and a slower increase of ν 2 (t). Furthermore, for intermediate values of the ratio η o /η l we observe a simultaneous increase of the values of ν 1 (t) and ν 2 (t), whereas sufficiently large and sufficiently small values of η o /η l correlate with a decoupling between the increase of ν 1 (t) and ν 2 (t). In more detail, if the ratio η o /η l is sufficiently high, first ν 2 (t) increases while ν 1 (t) remains almost constant and then, when ν 2 (t) is sufficiently high, ν 1 (t) starts increasing as well. On the other hand, in the case where η o /η l is sufficiently low, we have that ν 1 (t) increases first and then ν 2 (t) starts increasing as soon as ν 1 (t) becomes sufficiently high. These results communicate the biological notion that: the strength of the selective pressures exerted by oxygen and lactate on tumour cells, which are quantified by the values of the selection gradients η o and η l , may shape the emergence of hypoxic resistance and acidic resistance in tumours; the order in which such forms of resistance develop depend on the intensity of oxygen-driven selection in relation to the intensity of lactate-driven selection.

Conclusions and research perspectives
In this work, we have developed a mathematical modelling approach to investigate the influence of 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 data recapitulate the eco-evolutionary spatial dynamics of tumour cells and their adaptation to hypoxic and acidic microenvironments. In particular, the results obtained indicate that tumour cells characterised by lower levels of expression of hypoxia-and acidity-resistant genes are to be expected to colonise well-oxygenated and mildly-acidic regions of vascularised tissues, whereas cells expressing a more aggressive phenotype characterised by higher levels of resistance to hypoxia and acidity will ultimately populate tissue regions corresponding to hypoxic and acidic microenvironments. Such theoretical findings recapitulate histological data on ductal carcinoma in situ showing that the levels at which the acidity-resistant gene LAMP2 and the hypoxia-resistant gene GLUT-1 are 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 more acidic) [19,26].
Moreover, our theoretical findings reconcile the conclusions of [26], suggesting that tumour cells 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 that the order in which resistance to hypoxia and resistance to acidity arise depend on the ways in which oxygen and lactate act as environmental stressors in the evolutionary dynamics of tumour cells, 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 cell movement and mechanical interactions between cells [2,8,12], which would make it possible to investigate 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 of fluctuations in the inflow rate of oxygen and glucose and the outflow rate of lactate in the evolution of cancer cells. In fact, when in hypoxic conditions, cancer cells are known to produce and secrete proangiogenic 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 whereby the rates at which oxygen and glucose enter the tumour and the rate at which lactate is flushed out through intra-tumour blood vessels fluctuate over time, which impacts on the evolutionary dynamics of cancer cells [21,35,49,51].
It would also be interesting to extend the model in order to investigate the role of phenotypic transitions triggered by hypoxia and acidity -such as the epithelial-mesenchymal transition induced by hypoxic environmental conditions [53,59,65] and the acquisition of the metastatic phenotype promoted by acidic microenvironments [19,20,22] -in the phenotypic adaptation of cancer cells and tumour growth.
Finally, since resistance to hypoxia is known to correlate with resistance to chemotherapy and radiotherapy [18,20,38,56,61], building on [15,42,44], it would be relevant for anticancer therapy to address numerical optimal control for an extended version of the model that takes into account the effect of chemotherapy and/or radiotherapy [6,55], which could inform the development of optimised cancer treatment protocols that exploit evolutionary and ecological principles [1,27,36,50].