Evolutionary dynamics in vascularised tumours under chemotherapy

We consider a mathematical model for the evolutionary dynamics of tumour cells in vascularised tumours under chemotherapy. The model comprises a system of coupled partial integro-differential equations for the phenotypic distribution of tumour cells, the concentration of oxygen and the concentration of a chemotherapeutic agent. In order to disentangle the impact of different evolutionary parameters on the emergence of intra-tumour phenotypic heterogeneity and the development of resistance to chemotherapy, we construct explicit solutions to the equation for the phenotypic distribution of tumour cells and provide a detailed quantitative characterisation of the long-time asymptotic behaviour of such solutions. Analytical results are integrated with numerical simulations of a calibrated version of the model based on biologically consistent parameter values. The results obtained provide a theoretical explanation for the observation that the phenotypic properties of tumour cells in vascularised tumours vary with the distance from the blood vessels. Moreover, we demonstrate that lower oxygen levels may correlate with higher levels of phenotypic variability, which suggests that the presence of hypoxic regions supports intra-tumour phenotypic heterogeneity. Finally, the results of our analysis put on a rigorous mathematical basis the idea, previously suggested by formal asymptotic results and numerical simulations, that hypoxia favours the selection for chemoresistant phenotypic variants prior to treatment. Consequently, this facilitates the development of resistance following chemotherapy.


Introduction
Previous empirical and theoretical work has suggested that spatial variations in oxygen levels can foster the emergence of intra-tumour phenotypic heterogeneity [4,25,28,32,49,55,57,70]. In particular, it has been hypothesised that the nonlinear interplay between impaired oxygen delivery caused by structural abnormalities present in the tumour vasculature [21,24,41,42,60,72,73], limited oxygen diffusion and oxygen consumption by tumour cells may lead to the creation of distinct ecological niches in vascularised tumours, whereby tumour cells with different phenotypic characteristics can be selected [4,17,26,36,39,48]. This hypothesis is supported by a growing body of experimental and clinical studies indicating that: well-oxygenated parts of the tumour are densely populated by cells characterised by higher oxygen uptake and faster proliferation via aerobic pathways; hypoxic parts of the tumour (i.e. regions where oxygen levels are below normal physiological levels) are mainly occupied by cells that display higher levels of hypoxia-inducible factors, such as HIF-1 [21,29,45,60,65,66,68,71,77], which typically correlate with slower proliferation via anaerobic pathways and higher levels of resistance to chemotherapy.
In this paper, we use a mathematical model for the evolutionary dynamics of tumour cells in vascularised tumours under chemotherapy to gain a deeper understanding of the adaptive process that underpins the emergence of intra-tumour phenotypic heterogeneity and the development of resistance to chemotherapeutic agents.
The model comprises a system of coupled partial integro-differential equations for the phenotypic distribution of tumour cells, the concentration of oxygen and the concentration of a chemotherapeutic agent.
Non-local partial differential equations (PDEs) similar to the one that governs the evolution of the phenotypic distribution of tumour cells in our model and related integro-differential equations have recently received increasing attention from the mathematical community -see for instance [2,3,9,12,13,14,16,22,37,40,56]. In particular, our work follows earlier papers on the analysis and numerical simulation of integro-differential equations and non-local PDEs modelling the emergence of intra-tumour phenotypic heterogeneity [40,52,53,56,74]. The focus of these papers is on the cases where tumour cells do not change their phenotypic state or the rate of phenotypic variation is small. The main novelty of our work is that we allow tumour cells to undergo spontaneous epimutations (i.e. heritable phenotypic variations that occur randomly due to non-genetic instability and are not induced by any selective pressure [38]) and we do not impose any smallness assumptions on the rate at which such phenotypic changes occur. In this more general scenario, building upon the method of proof presented in [5,7,18,51], we carry out an analytical study of evolutionary dynamics. In particular, we construct explicit solutions to the equation for the phenotypic distribution of tumour cells and, considering the case where the concentrations of oxygen and of chemotherapeutic agent are stationary, we provide a detailed quantitative characterisation of the long-time asymptotic behaviour of such solutions. The analytical results obtained are integrated with numerical simulations of a calibrated version of the model based on biologically consistent parameter values, in order to further assess the impact of the dynamics of oxygen and of chemotherapeutic agent on the phenotypic evolution of tumour cells.
The paper is organised as follows. In Section 2, we introduce the equations of the model and the underlying modelling assumptions. In Section 3, we present the results of our analytical study of evolutionary dynamics. In Section 4, we report on a sample of numerical solutions that confirm and extend the analytical results obtained. Section 5 concludes the paper and provides a brief overview of possible research perspectives.

Mathematical model
We model the evolution of tumour cells within a region of a vascularised tumour along with the dynamical interactions that occur between tumour cells and both oxygen and a chemotherapeutic agent, which are released from the intra-tumoural vascular network.
The tumour region is approximated as a bounded set Ω ⊂ R d , with smooth boundary ∂Ω, where d = 1, 2, 3 depending on the biological scenario under study. The spatial position of tumour cells is described by a vector x ∈ Ω and every cell is characterised by the expression levelŷ ∈ R >0 of a hypoxia-inducible factor. We assume that there is a sufficiently high level of expression of the hypoxia-inducible factorŷ H conferring both the highest rate of cellular division via anaerobic energy pathways and the highest level of resistance to chemotherapy, while there is a sufficiently low level of expression of the hypoxia-inducible factorŷ L <ŷ H providing the highest rate of cellular division via aerobic energy pathways. Under these assumptions, we model the phenotypic state of tumour cells by means of a scalar variable y ∈ R which represents a rescaled level of expression of the hypoxia-inducible factor so that values of y → 1 correspond to phenotypic variants with higher rates of cellular division via anaerobic energy pathways and higher levels of chemoresistance (i.e. anaerobic and chemoresistant phenotypic variants), whereas values of y → 0 correspond to phenotypic variants with higher rates of cellular division via aerobic energy pathways (i.e. aerobic phenotypic variants).
The phenotypic distribution of tumour cells at time t ∈ [0, ∞) and position x is described by the function n(t, x, y), while the functions s(t, x) and c(t, x) describe, respectively, the oxygen concentration and the concentration of the chemotherapeutic agent at time t and position x. Moreover, at each time t, we define the density of tumour cells at position x as while the local mean phenotypic state and the related variance are defined, respectively, as

Dynamics of tumour cells
The local phenotypic distribution of tumour cells n(t, x, y) is governed by the following non-local PDE In the reaction-diffusion equation (3), the diffusion term models the effect of spontaneous epimutations, which occur at rate β > 0 [19,50], while the non-local reaction term models the effect of cell division and death. The function R y, ρ(t, x), s(t, x), c(t, x) represents the fitness of tumour cells in the phenotypic state y at position x and time t under the local environmental conditions given by the cell density ρ(t, x), the oxygen concentration s(t, x) and the concentration of chemotherapeutic agent c(t, x) (i.e. the phenotypic fitness landscape of the tumour at position x and time t). In particular, we consider with p(y, s) := f (y) + g(y, s).
Here, ζ > 0, f (y) is a C 2 -function such that g(y, s) is a C 2 -function of y and a C 1 -function of s that satisfies the following assumptions arg max y∈R g(y, s) = 0, g(0, s) > 0, ∂ 2 yy g(·, s) < 0 ∀ s ∈ (0, ∞) and lim s→∞ g(0, s) > f (1) and k(y, c) is a C 2 -function of y and a C 1 -function of c that satisfies the following assumptions arg min Definition (4) along with assumptions (9) and (10) models a biological scenario whereby the background fitness of tumour cells in the phenotypic state y at position x and time t is given by a function p(y, s(t, x)), the value of which is reduced: • due to competition for limited space, by a certain amount which is the same for all phenotypic variants and is proportional to ρ(t, x), with a proportionality constant ζ that is related to the local carrying capacity of the tumour; • due to the cytotoxic action of the chemotherapeutic agent, by a certain amount k(y, c) which increases monotonically with the concentration of the chemotherapeutic agent c and is smaller for phenotypic variants with y → 1, which are characterised by higher levels of chemoresistance, and is null for the phenotypic variant corresponding to y = 1, since such a phenotypic variant is assumed to be completely resistant to the chemotherapeutic agent.
Definition (5) corresponds to the case where the background fitness p(y, s) is defined as a linear combination of the background fitness associated with anaerobic energy pathways f (y) and that associated with aerobic energy pathways g(y, s). In particular, assumptions (6)-(8) translate into mathematical terms the following biological ideas: • The state y = 1 corresponds to the phenotypic variant with the maximal background fitness associated with anaerobic energy pathways, whereas the state y = 0 corresponds to the phenotypic variant with the maximal background fitness associated with aerobic energy pathways.
• Due to the fact that less fit phenotypic variants are driven to extinction by natural selection, the background fitness associated with anaerobic (or aerobic) energy pathways can be negative for phenotypic variants with values of y sufficiently far from 1 (or 0).
• Because of the fitness cost associated with a less efficient anaerobic metabolism [11], the maximal background fitness of aerobic phenotypic variants in well-oxygenated environments is larger than the maximal background fitness of anaerobic phenotypic variants.
• In the absence of oxygen, the background fitness p(y, s) coincides with the background fitness associated with anaerobic energy pathways f (y).
• The larger is the oxygen concentration, the stronger is the impact of the background fitness associated with aerobic energy pathways g(y, s) on the background fitness p(y, s).
In particular, following the modelling strategies presented in [52], here we use the definitions where ϕ > 0 is the maximal background fitness of anaerobic phenotypic variants, γ s > ϕ is the maximal background fitness of aerobic phenotypic variants, α s > 0 and α c > 0 are the Michaelis-Menten constants of oxygen and of chemotherapeutic agent respectively, and γ c > 0 is the maximal reduction of the background fitness due to the cytotoxic action of the chemotherapeutic agent. Definitions (11) satisfy assumptions (6)-(10), ensure analytical tractability of the model and lead to a fitness function R y, ρ, s, c that is close to the approximate fitness landscapes which can be inferred from experimental data through regression techniques -see, for instance, equation (1) in [59]. In fact, with these definitions, after a little algebra, the difference p(y, s) − k(y, c) in (4) can be rewritten as where a(s, c) and Here, a(s(t, x), c(t, x)) is the maximum fitness, h(s(t, x), c(t, x)) is the fittest phenotypic state and b(s(t, x), c(t, x)) is a nonlinear selection gradient at position x and time t under the environmental conditions corresponding to the oxygen concentration s(t, x) and the concentration of chemotherapeutic agent c(t, x).

Dynamics of abiotic factors
The oxygen concentration s(t, x) and the concentration of chemotherapeutic agent c(t, x) are governed by the following partial integro-differential equations and coupled with (3) and subject to zero-flux boundary conditions, i.e. ∇ x s · u = 0 and ∇ x c · u = 0 on ∂Ω, where u is the unit normal to ∂Ω that points outward from Ω. In (16) and (17), the parameters D s > 0 and D c > 0 are the diffusion coefficients of oxygen and of chemotherapeutic agent, the functions r s (y, s) and r c (y, c) are the consumption rates of oxygen and of chemotherapeutic agent by tumour cells in the phenotypic state y, the parameters λ s > 0 and λ c > 0 are the natural decay rates of oxygen and of chemotherapeutic agent, and the source terms q s (t, x) and q c (t, x) model the influx of oxygen and of chemotherapeutic agent from the intra-tumoural blood vessels at position x ∈ Ω and at time t. We assume that oxygen is consumed only by phenotypic variants corresponding to values of y for which the fitness associated with aerobic energy pathways g(y, s) is positive and we let oxygen consumption occur at a rate proportional to g(y, s). Moreover, we assume that the chemotherapeutic agent is consumed by phenotypic variants corresponding to different y at different rates proportional to the amount k(y, c) by which their background fitness is reduced due to the cytotoxic action of the chemotherapeutic agent. In accordance with these assumptions, we use the following definitions r s (y, s) := η s (g(y, s)) + and r c (y, s) := η c k(y, c), where η s > 0 and η c > 0 are constants of proportionality and (·) + denotes the positive part of (·). As done in [74], we let ω ⊂ Ω be the set of points within the tumour tissue which are occupied by blood vessels and, since we do not consider the formation of new blood vessels, we assume ω to be given and remain constant in time. Therefore, we define the source terms q s and q c as where 1 ω is the indicator function of the set ω, and i s (t, x) and i c (t, x) are the rates of inflow of oxygen and of chemotherapeutic agent through intra-tumoural blood vessels at position x ∈ ω and time t.

Analysis of evolutionary dynamics
In order to obtain an analytical description of the evolutionary dynamics of tumour cells, in this section we construct explicit solutions of (3) and study the long-time asymptotic behaviour of such solutions in the case where the concentrations of oxygen and of chemotherapeutic agent are stationary, i.e. we let with S(x) and C(x) being given functions such that Moreover, in agreement with much of the previous work on the mathematical analysis of the evolutionary dynamics of continuously-structured populations [61,63], we consider the case where at time t = 0 the local phenotypic distribution of tumour cells is of the following Gaussian form where Our main results are summarised by Theorem 1, the proof of which relies on the results established by Proposition 1. (4) and (12)- (15), the non-local PDE (3) complemented with (23) and (24) admits the exact solution

Proposition 1. Under assumptions
and Proof. Substituting (4) and (12) into (3) yields Building upon the results presented in [5,18,51], we make the ansatz (25). Substituting this ansatz into (29) and introducing the notation v(t, Equating the second-order terms in y gives the following differential equation for v alone Moreover, equating the coefficients of the first-order terms in y, and eliminating ∂ t v from the resulting equation, yields Lastly, choosing y = µ in (30) gives and eliminating ∂ t v from (33) we obtain with the function F (t, x) being defined according to (27). Under the initial condition (23), we have . Imposing these initial conditions on (31), (32) and (34), we arrive at the Cauchy problem (26) for the functions v(t, x), µ(t, x) and ρ(t, x). Theorem 1. Under assumptions (4), (12)-(15), (21) and (22), the solution of the Cauchy problem defined by complementing (3) with (23) and (24) is of the Gaussian form (25) and such that for every x ∈ Ω with and where the functions a(S, C), b(S, C) and h(S, C) are defined according to (13)- (15).
Proof. Proposition 1 ensures that for any (t, x) ∈ [0, ∞) × Ω the solution of the Cauchy problem defined by complementing (3) with (23) and (24) is of the Gaussian form (25). Therefore, building upon the method of proof presented in [7,18], we prove Theorem 1 by studying the behaviour of the components of the solution to the Cauchy problem (26) for t → ∞. Throughout the proof we use the abridged notation Step 1: asymptotic behaviour of v(t, x) for t → ∞. Solving (26) for Step 2: asymptotic behaviour of which implies that Step 3: asymptotic behaviour of ρ(t, x) for t → ∞. Under assumptions (21) and (22), using the abridged notation (39), we rewrite (27) as Substituting into (26) 3 and defining Solving (44) subject to the initial condition ρ(0, x) = ρ 0 (x) yields The asymptotic results (41) and (43) ensure that and, therefore, (45) allows us to conclude that On the other hand, the asymptotic result (46) implies that in the asymptotic regime t → ∞ we have and also that, under the additional assumption b(x) β < a(x), for some positive function A(x). These asymptotic relations, along with (45), allow us to conclude that Taken together, the asymptotic results (47) and (48) ensure that Claims (35)- (38) follow from the asymptotic results (41), (43) and (49).
The asymptotic results established by Theorem 1 provide a mathematical formalisation of the idea that, when the concentrations of oxygen and of chemotherapeutic agent are stationary (i.e. s(t, x) ≡ S(x) and c(t, x) ≡ C(x)), the tumour cell density ρ(t, x), the local mean phenotypic state µ(t, x) and the related variance σ 2 (t, x) converge to some equilibrium values ρ ∞ (x), µ ∞ (x) and σ 2 ∞ (x), respectively, which are determined by the concentration of oxygen and the concentration of chemotherapeutic agent - ). This is illustrated by the heat maps in Figure 1, which show how, for the biologically consistent parameter values listed in Table 1, the values of ρ ∞ , µ ∞ and σ 2 ∞ given by (36)-(38) vary as functions of S and C. Notice that the parameter values in Table 1 are such that ρ ∞ > 0.
These results demonstrate that spatial variations of the oxygen concentration determine spatial variations of the tumour cell density, of the local mean phenotypic state and of the related variance. Specifically, under the parameter values listed in Table 1, the tumour cell density ρ ∞ is an increasing function of the oxygen concentration. Moreover, the local mean phenotypic state µ ∞ decreases from values close to y = 1 (i.e. the state corresponding to the phenotypic variant with the highest rate of cellular division via anaerobic energy pathways) to values close to y = 0 (i.e. the state corresponding to the phenotypic variant with the highest rate of cellular division via aerobic energy pathways) for increasing values of the oxygen concentration. This suggests that aerobic phenotypic variants are to be expected to colonise oxygenated regions of the tumour, while anaerobic phenotypic variants are likely to populate poorly-oxygenated regions. Finally, the local phenotypic variance σ 2 ∞ is a decreasing function of the oxygen concentration, which supports the idea that higher levels of phenotypic variability may occur in hypoxic regions of the tumour.  ) and c(t, x) that converge to some limits s ∞ (x) and c ∞ (x) as t → ∞, we expect the long-time asymptotic limit of the local phenotypic distribution of tumour cells n(t, x, y) to be of the Gaussian form and .

Numerical simulations
We complement the analytical results of evolutionary dynamics established in Section 3 with numerical solutions of the model equations. In Section 4.1, we describe the set-up of numerical simulations and the methods employed to construct numerical solutions. In Section 4.2, we consider the case of a one-dimensional spatial domain whereby the concentrations of oxygen and of chemotherapeutic agent are stationary. In Section 4.3, we focus on the case of a two-dimensional spatial domain and let the dynamics of oxygen and of chemotherapeutic agent be governed by (16) and (17). All simulations are carried out using the parameter values listed in Table 1, which are chosen to be consistent with the existing literature.

Set-up of numerical simulations and numerical methods
Set-up of numerical simulations of Section 4.2. For the numerical simulations we report on in Section 4.2, we define Ω := (0, 0.05) and assume that increasing values of x ≡ x correspond to increasing values of the distance from a blood vessel located in x = 0. Under the parameter values listed in Table 1, the values of x are in units of cm. Coherently with assumptions (21) and (22), we let the concentrations of oxygen and of chemotherapeutic agent be stationary and given by with the functions S(x) and C(x) defined as shown by the plots in Figure 2. Here, the oxygen concentration S(x) is defined in such a way as to match the experimental oxygen distribution presented in [34, Fig. 3]. Furthermore, the concentration of chemotherapeutic agent C(x) is defined in such a way as to have a behaviour qualitatively similar to that of S(x) and the value of C(0) is chosen in agreement with experimental data presented in [34].  Fig. 3]. The conversion from mmHg of pO 2 to g cm −3 of oxygen concentration was performed using the conversion factor 1 mmHg = 4.6 × 10 −8 g cm −3 , which was estimated using the ideal gas law. The concentration of chemotherapeutic agent C(x) is defined in such a way as to have a behaviour which is qualitatively similar to that of S(x) and the value of C(0) is chosen in agreement with experimental data presented in [34].
We complement (3) with the initial condition (23) and assume Assumptions (54) Table 1, the values of x ∈ Ω are in units of cm. We let the dynamics of oxygen and of chemotherapeutic agent be governed by (16)- (18). Moreover, we assume the rate of inflow of oxygen and of chemotherapeutic agent through intra-tumoural blood vessels to be constant in time and the same for all vessels, i.e. we define the functions i s (t, x) and i c (t, x) in (20) as with the values of I s and I c being those given in Table 1.
We complement (3) with the initial condition defined via (23) and (54), while (16) and (17) are complemented with the following initial conditions with the values of S 0 and C 0 being those given in Table 1. These initial conditions correspond to a biological scenario whereby at the initial time t = 0 tumour cells are uniformly distributed across the spatial domain Ω and are mainly found in the phenotypic state y = 0.5, while the oxygen and the chemotherapeutic agent are concentrated in correspondence of the blood vessels.  is based on an explicit finite difference scheme in which a three-point stencil is used to approximate the diffusion term in y and an explicit finite difference scheme is used for the non-local reaction term [46]. Furthermore, the method for solving numerically (16) and (17) subject to the zero-flux boundary conditions (18) is based on an explicit finite difference scheme whereby a five-point stencil is used to approximate the diffusion terms and an explicit finite difference scheme is used for the other terms [46]. Finally, numerical solutions to the Cauchy problem (26) are constructed using the explicit Euler method. All numerical computations are performed in Matlab.

One-dimensional numerical results under stationary concentrations of oxygen and of chemotherapeutic agent
The sample of numerical results presented in Figure 3 refer to the case where the oxygen concentration s(t, x) ≡ S(x) and the concentration of cytotoxic agent c(t, x) ≡ 0, while the results presented in Figure 4 refer to the case where s(t, x) ≡ S(x) and c(t, x) ≡ C(x), with S(x) and C(x) being defined as illustrated by the plots in Figure 2.
Agreement between analytical and numerical results. In agreement with the results established by Proposition 1, the numerical results displayed in the top rows of Figure 3 and Figure 4 show that there is a perfect match between the cell density ρ(t, x), the local mean phenotypic state µ(t, x) and the related variance σ 2 (t, x) computed via numerical integration of the local cell phenotypic distribution n(t, x, y), which is obtained by solving numerically (3) subject to the initial condition defined via (23) and (54), and the corresponding quantities obtained by solving numerically the Cauchy problem (26) complemented with (54). Similarly, the sample of numerical results presented in the bottom rows of Figure 3 and Figure 4 show that the local cell phenotypic distribution n(t, x, y) matches the exact local cell phenotypic distribution (25). Moreover, in accordance with the asymptotic results established by Theorem 1, the cell density, the local mean phenotypic state and the related variance converge, respectively, to the equilibrium values ρ ∞ (x), µ ∞ (x) and σ 2 ∞ (x) given by (36)- (38).
Tumour cell dynamics in the absence of chemotherapeutic agent. The numerical results of Figure 3 show that, in the absence of chemotherapeutic agent, since the stationary oxygen concentration S(x) decreases monotonically with the distance from the blood vessel located at x = 0 (vid. Figure 2), the cell density ρ(t, x) at equilibrium is maximal in the vicinity of the blood vessel (cf. red line), where the oxygen concentration is higher, and decreases monotonically as the distance from the vessel increases (cf. blue and green lines). Accordingly, the local mean phenotypic state at equilibrium increases from values close to y = 0 (i.e. the state corresponding to the phenotypic variant with the highest rate of cellular division via aerobic energy pathways) to values close to y = 1 (i.e. the state corresponding to the phenotypic variant with the highest rate of cellular division via anaerobic energy pathways) moving away from the blood vessel. Moreover, the local phenotypic variance σ 2 (t, x) at equilibrium is a monotonically increasing function of the distance from the blood vessel (i.e. local phenotypic variability increases with the distance from the blood vessel).  Table 1.
Tumour cell dynamics in presence of chemotherapeutic agent. A comparison of the numerical results of Figure 3 and Figure 4 reveals that in the regions in close proximity of the blood vessel (cf. red lines), where its concentration is higher, the chemotherapeutic agent leads to the occurrence of a population bottleneck in tumour cells, which results into: a reduction of the equilibrium value of the cell density ρ(t, x); a selective sweep toward more resistant phenotypic variants, as demonstrated by the fact that the equilibrium value of the local mean phenotypic state µ(t, x) shifts from values closer to y = 0 (i.e. the state corresponding to the phenotypic variant with the highest rate of cellular division via aerobic energy pathways) to values closer to y = 1 (i.e. the state corresponding to the anaerobic phenotypic variant with the highest level of resistance to chemotherapy); a reduction of the equilibrium value of the local phenotypic variance σ 2 (t, x). Moreover, moving away from the blood vessel, since its concentration decreases, the chemotherapeutic agent has a weaker impact on the dynamics of tumour cells (cf. blue lines). As a result, the evolution of tumour cells in regions distal to the blood vessel is hardly affected by the chemotherapeutic agent (cf. green lines).  Table 1.

Two-dimensional numerical results under dynamical concentrations of oxygen and of chemotherapeutic agent
In the remainder of this section, we use the notation x ≡ (x 1 , x 2 ). The sample of numerical results presented in Figure 5 and Figure 6 refer to the case where the oxygen concentration s(t, x 1 , x 2 ) is governed by (16), subject to the initial condition (56) and the boundary condition (18), while the concentration of chemotherapeutic agent c(t, x 1 , x 2 ) ≡ 0. On the other hand, the results presented in Figure 7 and Figure 8 refer to the case where s(t, x 1 , x 2 ) and c(t, x 1 , x 2 ) are governed by (16) and (17), respectively, subject to the initial conditions (56) and the boundary conditions (18). In both cases, the set of points within the tumour tissue which are occupied by blood vessels (i.e. the set ω) is defined as illustrated by the plots in the first panels of Figure 5 and Figure 7.
Agreement between analytical and numerical results. The sample of numerical results presented in Figure 5 and Figure 7 show that, in the case of constant influx from intra-tumoural blood vessels, the concentration of oxygen s(t, x 1 , x 2 ) and the concentration of chemotherapeutic agent c(t, x 1 , x 2 ) obtained by solving numerically (16) and (17), subject to the initial conditions (56) and the boundary conditions (18), converge to some equilibria s ∞ (x 1 , x 2 ) and c ∞ (x 1 , x 2 ). As a result, in agreement with our expectation based on the results established by Theorem 1 (cf. Remark 1), the cell density ρ(t, x 1 , x 2 ) and the local mean phenotypic state µ(t, x 1 , x 2 ) computed via numerical integration of the local cell phenotypic distribution n(t, x 1 , x 2 , y), which is obtained by solving numerically (3) subject to the initial condition defined via (23) and (54), converge to the equilibrium values ρ ∞ (x 1 , x 2 ) and µ ∞ (x 1 , x 2 ) given by (51) and (52). Moreover, the sample of numerical results presented in Figure 6 and Figure 8 show that the local phenotypic distribution of tumour cells n(t, x 1 , x 2 , y) converges to the equilibrium phenotypic distribution n ∞ (x 1 , x 2 , y) given by (50).
Emergence of spatial gradients of oxygen and of chemotherapeutic agent. The numerical results of Figure 5 and Figure 7 show that, as one would expect based on the experimental results presented by [34], the equilibrium concentration of oxygen s(T, x 1 , x 2 ) and the equilibrium concentration of chemotherapeutic agent c(T, x 1 , x 2 ) are maximal in the vicinity of the blood vessels and decrease monotonically with the distance from the blood vessels. Moreover, these results demonstrate that the nonlinear interplay between the spatial distribution of the blood vessels, the reaction-diffusion dynamics of oxygen and of chemotherapeutic agent, and their consumption by tumour cells leads naturally to the emergence of spatial inhomogeneities in the equilibrium concentrations of such abiotic factors.
Tumour cell dynamics. The plots in Figures 5-8 demonstrate that the qualitative behaviour of the numerical results obtained under stationary concentrations of oxygen and of chemotherapeutic agents displayed in Figure 3 and Figure 4 remains unchanged when dynamical concentrations of oxygen and of chemotherapeutic agent are considered. Specifically, in the absence of chemotherapy, when moving away from the blood vessels, the equilibrium value of the cell density ρ(t, x 1 , x 2 ) decreases, the local mean phenotypic state µ(t, x 1 , x 2 ) at equilibrium increases from values close to y = 0 to values close to y = 1, and the equilibrium value of the related variance σ 2 (t, x 1 , x 2 ) increases (vid. Figure 5 and Figure 6). When chemotherapy is administered, its effect is more pronounced in the proximity of the blood vessels and consists in a reduction of the equilibrium value of ρ(t, x 1 , x 2 ), a shift of the equilibrium value of µ(t, x 1 , x 2 ) toward y = 1 and a reduction of the equilibrium value of σ 2 (t, x 1 , x 2 ) compared to the case where the chemotherapeutic agent is not present. Moreover, the evolutionary dynamics of tumour cells is weakly affected by chemotherapy in regions far from the blood vessels, where the concentration of chemotherapeutic agent is lower (vid. Figure 7 and Figure 8).  Table 1.

Conclusions and research perspectives
Conclusions. The results of our analysis of evolutionary dynamics recapitulate previous theoretical results [4,6,8,25,32,39,43,52,53,55,70,74] and experimental data [60,66,70,71] by demonstrating that spatial inhomogeneities in the concentration of oxygen promote the selection of different phenotypic variants at  (3) and (16) imposing the initial conditions defined via (23), (54) and (56), the boundary condition (18) and assuming c(t, x 1 , x 2 ) ≡ 0 (i.e. in the absence of chemotherapeutic agent). The set ω in (20) consists of the parts of Ω highlighted in red in the first panel of Figure 5. In the second, third and fourth panels, the bullets on the axis of abscissas highlight the value of the local mean phenotypic state µ(T, x 1 , x 2 ) and the black, dashed lines highlight the asymptotic limit (50) (51), (52) and (53) Table 1. different positions within the tumour. More specifically, our analytical results indicate that the tumour tissue in the vicinity of blood vessels is to be expected to be densely populated by aerobic phenotypic variants, while poorly oxygenated regions of the tumour are more likely to be sparsely populated by anaerobic phenotypic variants. Furthermore, the analytical results obtained support the idea that higher levels of phenotypic variability may occur in hypoxic regions of the tumour, which provides a theoretical basis for experimental results such as those presented by [10].
Coherently with observations made in previous theoretical and experimental studies [1,15,62,69,76], our analytical results also suggest that hypoxia favours the selection for chemoresistant phenotypic variants prior to treatment. Consequently, this facilitates the development of resistance following chemotherapy. Moreover, these results put on a rigorous mathematical basis the idea, previously suggested by formal analysis and numerical simulations [52,64], that chemotherapy removes the selective barrier limiting the growth of chemoresistant phenotypic variants by killing aerobic phenotypic variants in well-oxygenated regions of the tumour.
The results of our analysis of evolutionary dynamics are confirmed by the sample of numerical results presented. Such numerical results also indicate that gradients of oxygen and of chemotherapeutic agents, which are released from the intra-tumoural vascular network, naturally emerge in vascularised tumours due to the nonlinear interplay between the spatial distribution of the blood vessels, the reaction-diffusion dynamics of oxygen and of chemotherapeutic agents, and their consumption by tumour cells.
Research perspectives. We plan to extend our analytical results to the case where spatial movement of tumour cells is incorporated into the model. Based on the formal asymptotic results that we presented in [74], in the case where cell movement is modelled through Fick's first law, we expect the qualitative behaviour of  Table 1. the results obtained in this paper to remain unchanged in the asymptotic regime where the rate of spontaneous phenotypic variation and the cell diffusivity tend to zero. However, further developments of the method of proof employed here are required in order to carry out a similar analysis of evolutionary dynamics in more general scenarios.
It would also be useful to consider a discrete version of the continuum model presented here, whereby the dynamics of tumour cells would be described in terms of a branching random walk, while the concentrations of oxygen and of chemotherapeutic agent would be governed by discrete balance equations. This would make it possible to study stochastic effects which are relevant at low tumour cell densities and cannot be easily captured  (3), (16) and (17) imposing the initial conditions defined via (23), (54) and (56), and the boundary conditions (18). The set ω in (20) consists of the parts of Ω highlighted in red in the first panel of Figure 7. The white, dashed lines in the first and second panels highlight the 1D cross-section corresponding to x 2 = 0.4 and the bullets highlight the points (0.15, 0.4), (0.16, 0.4) and (0.3, 0.4). In the third, fourth and fifth panels, the filled bullets on the axis of abscissas highlight the value of the mean phenotypic state µ(T, x 1 , x 2 ), while the empty bullets highlight the corresponding values obtained in the case where c(t, x 1 , x 2 ) ≡ 0 (i.e. in the absence of chemotherapeutic agent). Moreover, the black, dashed lines highlight the asymptotic limit (50) with ρ ∞ (x), µ ∞ (x) and σ 2 ∞ (x) computed through (51), (52) and (53) Table 1. by deterministic models formulated in terms of differential equations. In this regard, we plan to extend the formal methods that we developed in [67] in order to identify the corresponding discrete counterpart of our continuum model. Finally, building upon the ideas presented in [7,8], it would be interesting to study the effect on the evolutionary dynamics of tumour cells of fluctuations in the rate of oxygen inflow, which are known to influence intra-tumour phenotypic heterogeneity [30,55,64]. It would also be interesting to include the effect of temporal variations in the spatial distribution of intra-tumoural blood vessels, which would make it possible to explore the influence of angiogenesis on the evolutionary dynamics of tumour cells in vascularised tumours. In addition, the model considered here could be further developed to incorporate a more comprehensive description of cell metabolism that captures acidosis and enhanced tumour invasiveness caused by the presence of hypoxic cells [23,26,27,31,43,54,57,65,77].