Evolutionary dynamics of competing phenotype-structured populations in periodically fluctuating environments

Living species, ranging from bacteria to animals, exist in environmental conditions that exhibit spatial and temporal heterogeneity which requires them to adapt. Risk-spreading through spontaneous phenotypic variations is a known concept in ecology, which is used to explain how species may survive when faced with the evolutionary risks associated with temporally varying environments. In order to support a deeper understanding of the adaptive role of spontaneous phenotypic variations in fluctuating environments, we consider a system of non-local partial differential equations modelling the evolutionary dynamics of two competing phenotype-structured populations in the presence of periodically oscillating nutrient levels. The two populations undergo heritable, spontaneous phenotypic variations at different rates. The phenotypic state of each individual is represented by a continuous variable, and the phenotypic landscape of the populations evolves in time due to variations in the nutrient level. Exploiting the analytical tractability of our model, we study the long-time behaviour of the solutions to obtain a detailed mathematical depiction of the evolutionary dynamics. The results suggest that when nutrient levels undergo small and slow oscillations, it is evolutionarily more convenient to rarely undergo spontaneous phenotypic variations. Conversely, under relatively large and fast periodic oscillations in the nutrient levels, which bring about alternating cycles of starvation and nutrient abundance, higher rates of spontaneous phenotypic variations confer a competitive advantage. We discuss the implications of our results in the context of cancer metabolism.


Introduction
Organisms of various scales, ranging from bacteria to animals, exist in fluctuating environments. For example, in order to cope with changes in nutrient availability, they are required to adapt. When the fluctuations are regular and the populations have sufficient time to sense the changes and react, a highly plastic phenotype, which allows individuals in the population to acquire different traits based on environmental cues, is an optimal strategy (Xue and Leibler 2018). An alternative strategy that is more suitable for dealing with irregular and unpredictable changes in the environment is risk spreading, which is also known as bet-hedging (Cohen 1966). Here, the population diversifies its phenotypes such that each sub-population is adapted to a specific environment. This ensures that at least some fraction of the population will survive in the face of sudden environmental changes (Philippi and Seger 1989). Phenotypic heterogeneity, a characteristic feature of a risk spreading strategy, is observed in many systems, including bacterial populations (Kussell and Leibler 2005) and solid tumours (Marusyk et al. 2012).
Bet-hedging is typically proposed to occur within the context of bacterial populations, where experimental support for stochastic phenotype switching is available (Kussell and Leibler 2005;Smits et al. 2006;Veening et al. 2008;Beaumont et al. 2009;Acar et al. 2008). The classic example of bet-hedging is bacterial persistence. During antibiotic treatment a small fraction of slowly growing bacteria, that are resistant to the antibiotic, is able to survive. After the treatment is over, the original population is restored, resulting in resistance to the antibiotic (Balaban et al. 2004). Schreiber et al. (2016) showed that fluctuations in nutrient levels alter the metabolism of bacteria and promote phenotypic heterogeneity. Risk spreading strategies have been observed in other organisms, such as fungi and slime moulds (Kussell and Leibler 2005). It is also hypothesized to be present in cancer where irregular vasculature can cause significant fluctuations within the tumour microenvironment (Gillies et al. 2018). Experimental and theoretical work suggest that intermittent lack of oxygen, i.e. cycling hypoxia, leads to clonal diversity, promotes metastasis and selects for more aggressive phenotypes (Cairns et al. 2001;Robertson-Tessi et al. 2015;Chen et al. 2018;Nichol et al. 2016).
Mathematically, competition between populations evolving in fluctuating environments has been studied using different modelling approaches, including deterministic predator-prey models and stochastic models (Cushing 1980;Chesson 1994;Roxburgh et al. 2004;Anderies and Beisner 2000). Previous work suggests that the likelihood of species coexistence is increased by temporal variations in the environment. More recent models have looked at adaptive strategies, including stochastic phenotype switching, emerging in stochastic environments (Fudenberg and Imhof 2012;Müller et al. 2013;Ashcroft et al. 2014;Wienand et al. 2017;Gravenmier et al. 2018). Müller et al. (2013) theoretically investigated the environmental conditions that would lead to the emergence of bet-hedging, noting the importance of the fluctuation timescales on the success of the adaptation strategy. A rapidly fluctuating environment selects the phenotype that is adapted to averaged conditions, whereas in a slowly varying environment, having two distinct specialists is beneficial. The bet-hedging population was shown to be most successful in an environment that fluctuates on an intermediate timescale.
Most of the experimental and theoretical models that have been developed to explore the dynamics of phenotypic changes in fluctuating environments consider the state of the environment and the phenotypic state of the individuals to be binary-i.e. the environment switches between two extreme conditions and individuals are allowed to jump between two antithetical phenotypic states that are each adapted to opposing environmental conditions (Acar et al. 2008;Müller et al. 2013;Wienand et al. 2017). However, in many cases of biological and ecological interest Natura non facit saltus, and it might therefore be relevant to consider the occurrence of intermediate environmental conditions and the existence of a spectrum of possible phenotypic states.
In light of these considerations, we present here a novel mathematical model for the evolutionary dynamics of two competing phenotype-structured populations in periodically fluctuating environments. The phenotypic state of each individual is represented by a continuous variable, and the phenotypic fitness landscape of the populations evolves in time due to variations in the concentration of a nutrient. In order to assess the evolutionary role that heritable, spontaneous phenotypic variations play in environmental adaptation, we focus on the case where the two populations undergo such phenotypic variations with different probabilities.
In our model, the phenotype distribution of the individuals within each population is described by a population density function that is governed by a parabolic partial differential equation (PDE), whereby a linear diffusion operator models the occurrence of spontaneous phenotypic variations, while a non-local reaction term takes into account the effects of asexual reproduction and intrapopulation competition. The two non-local parabolic PDEs for the population density functions are coupled through an additional non-local term modelling the effects of interpopulation competition. In such a mathematical framework, the fact that the two populations undergo phenotypic variations with different probabilities translates into the assumption that the two PDEs have different diffusion coefficients.
Mathematical models formulated in terms of integrodifferential equations and nonlocal parabolic PDEs like those considered here have been increasingly used to achieve a more in-depth theoretical understanding of the mechanisms underlying phenotypic adaptation in a variety of biological contexts (Alfaro and Veruete 2018;Alfaro et al. 2017;Almeida et al. 2019;Bouin and Calvez 2014;Bouin et al. 2012;Busse et al. 2016;Cuadrado 2004, 2007;Calsina et al. 2013;Chisholm et al. 2016Chisholm et al. , 2015Cuadrado 2009;Delitala and Lorenzi 2012a, b;Delitala et al. 2013;Desvillettes et al. 2008;Diekmann et al. 2005;Domschke et al. 2017;Iglesias and Mirrahimi 2018;Lam 2017a, b;Lorenzi et al. 2014Lorenzi et al. , 2016Lorenzi et al. , 2015aLorenzi et al. , b, 2018Lorz et al. 2013Lorz et al. , 2015Nordmann et al. 2018;Perthame 2006;Pouchol et al. 2018;Turanova 2015). In particular, our work follows earlier papers on non-local parabolic PDEs modelling the evolutionary dynamics of populations structured by continuous traits in periodically-fluctuating environments (Iglesias and Mirrahimi 2018;Lorenzi et al. 2015a). Compared to these previous studies, which considered scalar equations modelling the dynamics of single population, our model comprises a system of coupled equations modelling the dynamics of competing populations. This requires a novel extension of the methods developed in Lorenzi et al. (2015a) to characterise the qualitative and quantitative properties of the solutions.
Exploiting the analytical tractability of our model, we study the long-time behaviour of the solutions in order to obtain a detailed mathematical depiction of the evolutionary dynamics. Moreover, the asymptotic results are compared to numerical solutions of the model equations. Our analytical and numerical results clarify the role of heritable, spontaneous phenotypic variations as drivers of adaptation in periodically fluctuating environments.
The paper is organised as follows. In Sect. 2 we introduce the mathematical model. In Sect. 3 we carry out an analytical study of evolutionary dynamics. In Sect. 4 we integrate the analytical results with numerical simulations. In Sect. 5 we discuss the biological relevance of our theoretical findings in the context of cancer cell metabolism and tumour-microenvironment interactions. Section 6 concludes the paper and provides a brief overview of possible research perspectives.

Model description
We study the evolutionary dynamics of two competing phenotype-structured populations in a well-mixed system. Individuals within the two populations reproduce asexually, die and undergo heritable, spontaneous phenotypic variations. We assume the two populations differ only in the rate at which they undergo phenotypic variations. We label the population undergoing phenotypic variations at a higher rate by the letter H , while the other population is labelled by the letter L.
We represent the phenotypic state of each individual by a continuous variable x ∈ R, and we describe the phenotype distributions of the two populations at time t ∈ [0, ∞) by means of the population density functions n H (x, t) ≥ 0 and n L (x, t) ≥ 0. We define the size of population H , the size of population L and the total number of individuals inside the system at time t, respectively, as (1) Moreover, we define, respectively, the mean phenotypic state and the related variance of each population i ∈ {H , L} at time t as In the mathematical framework of our model, the function σ 2 i (t) provides a measure of the level of phenotypic heterogeneity in the i th population. Finally, we introduce a function S(t) ≥ 0 to model the concentration of a nutrient that is equally available to the two populations at time t, which we assume is given.
The evolution of the population density functions is governed by the following system of non-local parabolic PDEs for (x, t) ∈ R × (0, ∞). (3) In the system of PDEs (3), the diffusion terms model the effects of spontaneous phenotypic variations, which occur at rates β H > 0 and β L > 0, with The functional R x, S(t), ρ(t) models the fitness of individuals in the phenotypic state x at time t under the environmental conditions given by the nutrient concentration S(t) and the total number of individuals ρ(t)-i.e. the functional R x, S(t), ρ(t) can be seen as the phenotypic fitness landscape of the two populations at time t. Throughout the paper, we define this fitness functional as Definition (5) translates into mathematical terms the following biological ideas: (i) all else being equal, individuals die due to interpopulation and intrapopulation competition at rate dρ(t), with the parameter d > 0 being related to the carrying capacity of the system in which the two populations are contained; (ii) individuals in the phenotypic state x proliferate and die under natural selection at rate p(x, S(t)) (i.e. the function p(x, S) is a net proliferation rate). We focus on a scenario corresponding to the biological assumptions given hereafter. Under Assumptions 1 and 2, we define the net proliferation rate as with 0 < ζ ≤ γ . The parameters γ and ζ model, respectively, the maximum proliferation rate of the phenotypic variants best adapted to nutrient-rich and nutrient-scarce environments. Definition (6) ensures analytical tractability of the model and leads to a fitness functional that is close to the approximate fitness landscapes which can be inferred from experimental data through regression techniques-see, for instance, equation (1) in Otwinowski and Plotkin (2014). In fact, after a little algebra, definition (6) can be rewritten as and Under the environmental conditions defined by the nutrient concentration S, the function 0 ≤ ϕ(S) ≤ 1 represents the fittest phenotypic state, γ g(S) > 0 is the maximum fitness, and h(S) can be seen as a nonlinear selection gradient that quantifies the intensity of natural selection. Throughout the paper we will refer to g(S) as the rescaled maximum fitness.
In accordance with Assumptions 1 and 2, Eq. (7) shows that definition (6) is such that the fittest phenotypic state ϕ(S) belongs to the interval [0, 1] for any nutrient concentration S ≥ 0, i.e. ϕ : R ≥0 → [0, 1]. In particular, under starvation conditions (i.e. if S = 0) the fittest phenotypic state is ϕ(0) = 1, while increasing nutrient concentrations correspond to values of the fittest phenotypic state closer to 0, i.e. ϕ (S) < 0 for all S ≥ 0 and ϕ(S) → 0 as S → ∞. Furthermore, the fact that the function p(x, S) is negative for values of x sufficiently far from the fittest phenotypic state ϕ(S) captures the idea that less fit variants are driven to extinction by natural selection. These observations are illustrated by the plots in Fig. 1.

Analysis of evolutionary dynamics
To obtain an analytical description of the evolutionary dynamics, we focus on a biological scenario whereby the initial phenotype distributions of the two populations are Gaussians, that is, we study the behaviour of the solution to the system of non-local parabolic equations (3) subject to the initial condition given by the pair n H (x, 0) and n L (x, 0) with

Remark 1
The choice of initial condition (12) is consistent with much of the previous work on the mathematical analysis of the evolutionary dynamics of continuous traits, which relies on the prima facie assumption that population densities are Gaussians (Rice 2004).
Before turning to the case of periodically fluctuating environments in Sect. 3.2, we consider the case of constant environments in Sect. 3.1. The proofs of the results presented in Sect. 3.1 and Sect. 3.2 rely on the results established by Proposition 1.

Proposition 1
Under assumptions (5), (7) and (11), the system of non-local PDEs (3) subject to the initial condition (12) admits the exact solution with the population size, ρ i (t), the mean phenotypic state, μ i (t), and the inverse of the related variance, v i (t) = 1/σ 2 i (t), being solutions of the Cauchy problem where Proof Substituting definitions (5), (7) and (11) into the non-local PDE (3) for n i (x, t) yields Building upon the results presented in Almeida et al. (2019), Chisholm et al. (2016) and Lorenzi et al. (2015a), we make the ansatz (13) and substituting this ansatz into Eq. (16) we find Equating the coefficients of the zero-order, first-order and second-order terms in x in (17) produces a system of differential equations. Namely, the second-order terms in x yield the following differential equation Moreover, equating the coefficients of the first-order terms in x, and eliminating v i from the resulting equation, yields Lastly, choosing x = μ i in Eq. (17) gives and eliminating v i from Eq. (20) we find with the function F i (t) being defined according to (15). Under the initial condition (12), we have Imposing these initial conditions on the system of differential equations (18)-(21), we arrive at the Cauchy problem (14) for the functions v i (t), μ i (t) and ρ i (t).

Evolutionary dynamics in constant environments
Focussing on the case of constant environments, we let the nutrient concentration be constant and thus we make the assumption which implies that In this case, our main results are summarised by Theorem 3.
Theorem 3 Under assumptions (4), (5), (7), (11) and the additional assumption (22), the solution of the system of PDEs (3) subject to the initial condition (12) is of the Gaussian form (13) and satisfies the following: and Proof Under the additional assumption (22), Proposition 1 ensures that the population density function n i (x, t) is of the Gaussian form (13) with the population size, ρ i (t), the mean phenotypic state, μ i (t), and the inverse of the related variance, v i (t) = 1/σ 2 i (t), being governed by the Cauchy problem (14) with g(t) ≡ g and ϕ(t) ≡ ϕ.
In this framework, we divide the proof of Theorem 3 into four steps. We study the asymptotic behaviour of v i (t), μ i (t) and F i (t) for t → ∞ (Step 1). We show that ρ i (t) is non-negative and uniformly bounded (Step 2). Finally, we prove claim (25) (Step 3), and we conclude with the proof of claims (27) and (28) (Step 4).
Step 1: asymptotic behaviour of v i (t), μ i (t) and F i (t) for t → ∞. Solving the separable first-order differential equation (14) 1 for v i (t) and imposing the initial condition (14) 4 gives Moreover, solving the differential equation (14) 2 for μ i (t) by the integrating factor method and imposing the initial condition (14) 4 yields from which, using the positivity of v i (t), we find that Lastly, noting that, under the additional assumption (22), the function F i (t) defined by (15) reads as the asymptotic results (30) and (32) allow us to conclude that Step 2: non-negativity and boundedness of ρ i (t). Solving the differential equation (14) 3 for ρ i and imposing the initial condition (14) 4 yields This result, along with the positivity of ρ 0 i , implies that Moreover, substituting (33) into the differential equation (14) 3 for ρ i yields Estimating from above the right-hand side of the latter differential equation by using the non-negativity of ρ j (t) [cf. the uniform lower bound (36)], the positivity of v i (t) [cf. expression (29)] and the fact that g < 2 [cf. definition (11)], we obtain the differential inequality which gives the uniform upper bound Step 3: proof of claim (25). Combining the asymptotic result (34) with the expression (35) for ρ i we find that for some positive constant C. Since the function ρ(t) is non-negative [cf. the uniform lower bound (36)], the asymptotic relation (38) ensures that Under assumption (4) and the additional assumption (24), claim (25) follows from the asymptotic result (39).
Step 4: proof of claims (27) and (28). As long as ρ H (t) > 0, we can compute the quotient of ρ L (t) and ρ H (t) through (35). In so doing we find Using the limit (34) for F i , we then have for some positive constant C. Under assumption (4), the asymptotic relation (41) gives and, since ρ L is uniformly bounded from above [cf. the uniform upper bound (37)], from (41) we conclude that We can rewrite the differential equation (14) 3 for ρ L as where the function η(t) is defined as Using the asymptotic results (30), (32) and (42), we see that Solving the differential equation (43) complemented with the initial condition ρ L (0) = ρ 0 L yields ) The result (44) ensures that in the asymptotic regime t → ∞ we have and, under the additional assumption (26), we also have for some positive constant C. These asymptotic relations, along with the expression (45) for ρ L , allow us to conclude that Claims (27) and (28) follow from the asymptotic results (42) and (46), and the asymptotic results (32) and (30) with i = L.
The asymptotic results established by Theorem 3 provide a mathematical formalisation of the idea that in constant environments: 1. populations undergoing spontaneous phenotypic variation at a rate that is too large compared to the maximum fitness will ultimately go extinct [cf. point (i)]; 2. ceteris paribus, if at least one population undergoes spontaneous phenotypic variation at a rate sufficiently small compared to the maximum fitness [cf. point (ii)] then: 2(a). the population with the lower rate of phenotypic variation will outcompete the other population; 2(b). the equilibrium phenotype distribution of the surviving population will be unimodal with the mean phenotype corresponding to the fittest phenotypic state and the related variance being directly proportional to the rate of phenotypic variations.

Evolutionary dynamics in periodically fluctuating environments
We now focus on the case of environments that undergo fluctuations with period T > 0, and we assume the nutrient concentration to be Lipschitz continuous and T -periodic, i.e. we let S : [0, ∞) → R ≥0 satisfy the assumptions which implies that the functions g(t) and ϕ(t) satisfy the assumptions (48) Our main results are summarised by Theorem 4, the proof of which relies on the results established by Lemmas 1 and 2. (5), (7), (11) and (47), the unique real T -periodic solution of the problem

Lemma 1 Under assumptions
is and satisfies the integral identity

the unique real non-negative T -periodic solution of the problem
is and satisfies the integral identity We refer the interested reader to Lorenzi et al. (2015a) for the proofs of Lemmas 1 and 2.
Theorem 4 Under assumptions (4), (5), (7), (11) and the additional assumptions (47), the solution of the system of PDEs (3) subject to the initial condition (12) is of the Gaussian form (13) and satisfies the following: and and with w i (t) and u i (t) given by (55) and (50), respectively.
Proof Proposition 1 ensures that the population density function n i (x, t) is of the Gaussian form (13) with the population size, ρ i (t), the mean phenotypic state, μ i (t), and the inverse of the related variance, v i (t) = 1/σ 2 i (t), being governed by the Cauchy problem (14). In this framework, we prove Theorem 4 in 4 steps. In Step 1, we study the asymptotic behaviour of v i (t), μ i (t) and F i (t) for t → ∞. In Step 2, we show that ρ i (t) is non-negative and uniformly bounded. In Step 3, we prove claim (58). Finally, we prove claims (60) and (61) in Step 4.
Step 1: asymptotic behaviour of v i (t), μ i (t) and F i (t) for t → ∞. Since the differential equation (14) 1 does not depend on S(t), the expression (29) of v i (t) obtained in the proof of Theorem 3 still holds and Moreover, using the asymptotic result (62) along with the linear differential equation (14) 2 for μ i (t) one can easily show that where u i (t) is a T -periodic solution of the differential equation (49). Lemma 1 ensures that u i (t) is given by (50). Lastly, for F i (t) defined according to (15), the asymptotic results (62) and (63) allow us to conclude that Step 2: non-negativity and boundedness of ρ i (t). Proceeding in a similar way as in the proof of Theorem 3 (cf. Step 2 in the proof of Theorem 3), one can prove that Step 3: proof of claim (58). Solving the differential equation (14) 3 for ρ i and imposing the initial condition (14) 4 yields with F i (t) defined according to (15). Combining the asymptotic result (64) with the expression (66) for ρ i (t) gives for some positive constant C. Hence, using the fact that the functions g(t), ϕ(t) and u i (t) are T -periodic and considering m → ∞, we find for some positive constant C. Since the function ρ(t) is non-negative [cf. the uniform lower bound (65)], the asymptotic relation (68) ensures that if This proves that if assumption (57) is satisfied then claim (58) is verified.
Step 4: proof of claims (60) and (61). Let As long as ρ j (t) > 0, we can compute the quotient of ρ i (t) and ρ j (t) through (66). In so doing, using the asymptotic relation (67) for ρ i (t) and ρ j (t) and considering m → ∞ we obtain for some positive constant C, with Λ i and Λ j defined according to (52). The asymptotic relation (71) allows us to conclude that Since ρ i is uniformly bounded from above [cf. the uniform upper bound (65)], the asymptotic result (72) implies that ρ j (t) → 0 exponentially fast as t → ∞.
We can rewrite the differential equation (14) 3 for ρ i as where the function η(t) is defined as Using the asymptotic results (62), (63) and (73) we see that η(t) → 0 exponentially fast as t → ∞. Hence, ρ i (t) →ρ i (t) as t → ∞, withρ i (t) being the solution of the differential equationρ subject to an initial condition 0 <ρ i (0) < ∞. We note that: (i) the functionρ i (t) is uniformly bounded as it satisfies the upper and lower bounds (ii) the function f is Lipschitz continuous in the first variable; (iii) the function f is Lipschitz continuous and T -periodic in the second variable, since g, ϕ and u i are Tperiodic Lipschitz continuous functions of t [cf. assumptions (48) and expression (50)]. Therefore, the conditions of Massera's Convergence Theorem (Massera et al. 1950;Smith 1986) are satisfied and this allows us to conclude that with w i (t) being a non-negative T -periodic solution of the differential equation (54). Under the additional assumption (59), that is, Lemma 2 ensures that w i (t) is given by (55). Claims (60) and (61) follow from the asymptotic results (73) and (76) along with the asymptotic result (63) and (62) In summary, Theorem 4 gives an explicit characterisation of the long-term limit of v i (t), μ i (t) and ρ i (t) for the surviving population i and shows that the surviving population is the one characterised by the larger positive value of the quantity Remark 2 Since the functions u i (t) and w i (t) are T -periodic and satisfy the integral identities (51) and (56), respectively, the results established by Theorem 4 show that the long-term limits of the size and the mean phenotypic state of the surviving population are periodic functions of time with period T and mean values given by (51) and (56), respectively.
Remark 3 Using the differential equation (49) for u i (t) one can easily obtain Integrating both sides of the above equation with respect to t between 0 and T , and using the fact that u Therefore, definition (52) can be rewritten as If S ≡S then g(t) ≡ g and ϕ(t) ≡ ϕ (i.e. ϕ ≡ 0). In this case,

ϕ(t) − u i (t)) ϕ (t)
is sufficiently small then Λ H > Λ L , while if such a mean value is sufficiently large then Λ H < Λ L . We expect the latter scenario to occur when the variability and the rate of change of S(t) are sufficiently high so as to cause substantial and sufficiently fast variations in the value of ϕ(t).
The asymptotic results established by Theorem 4 formalise mathematically the idea that in periodically fluctuating environments: 2. ceteris paribus, if at least one population undergoes spontaneous phenotypic variations at a rate sufficiently small compared to the mean value of the maximum fitness, then the following behaviours are possible: 2(a). when environmental conditions are relatively stable, the population with the lower rate of phenotypic variations will outcompete the other population [cf. point (ii) and Remark 3]; 2(b). when environmental conditions undergo drastic changes, either both populations go extinct [cf. point (i) and Remark 3] or the population with the higher rate of phenotypic variations will outcompete the other population [cf. point (ii) and Remark 3]; 2(c). the phenotype distribution of the surviving population will be unimodal, and both the population size and the mean phenotype will become periodic [cf. point (ii) and Remark 2]; 2(d). ultimately, the population size and the mean phenotype will both oscillate with the same period as the fluctuating environment, and the mean value (with respect to time) of the mean phenotype will be the same as the mean value of the fittest phenotypic state with the related variance being directly proportional to the rate of phenotypic variations [cf. point (ii) and Remark 2].
These biological implications are reinforced by the numerical solutions presented in the next section.

Numerical simulations
In this section, we construct numerical solutions to the system of non-local parabolic PDEs (3) subject to the initial condition (12). In Sect. 4.1, we describe the numerical methods employed and the set-up of numerical simulations. In Sect. 4.2, we present a sample of numerical solutions that confirm the results of our analysis of evolutionary dynamics, both in the case where S(t) is constant and when S(t) oscillates periodically.

Numerical methods and set-up of numerical simulations
We select a uniform discretisation consisting of 2000 points on the interval [− 5, 5] as the computational domain of the independent variable x and impose no flux boundary conditions. Moreover, we assume t ∈ [0, t f ], with t f > 0 being the final time of simulations, and we discretise the interval [0, t f ] with the uniform step Δt = 0.0001. The method for constructing numerical solutions to the system of non-local parabolic PDEs (3) is based on an explicit finite difference scheme in which a three-point stencil is used to approximate the diffusion terms and an explicit finite difference scheme is used for the reaction term (LeVeque 2007). On the other hand, we use the Matlab built-in solver ode45 to solve numerically the Cauchy problem (14) for v i (t), μ i (t) and ρ i (t).
The parameter values used to carry out numerical simulations are listed in Table 1. In summary, to capture the fact that rates of spontaneous phenotypic variation are small, in general, and much smaller than maximum proliferation rates, in particular, we assume β i γ for i ∈ {H , L}. Furthermore, given the values of γ and β i , we fix the value of d to be such that the long-term limit (27) of the size of the population L is approximatively 10 4 , which is consistent with biological data from the existing literature regarding in vitro cell populations (Voorde et al. 2019).
We remark that the value of the parameter β L and the range of values of the parameter β H reported in Table 1 are such that neither condition (24) nor condition (57) are met in all cases on which we report here. This ensures that the two populations do not simultaneously go extinct.
We consider both populations to have the same initial phenotypic distribution (12) with v 0 i = 20, μ 0 i = 0 and ρ 0 i ≈ 800 for i ∈ {H , L}.

Main results
We consider the following definition of the nutrient concentration In definition (79), the parameter M > 0 represents the mean nutrient concentration, while the parameter A ≥ 0 models the semi-amplitude of the oscillations of the nutrient concentration, which have period T > 0. Clearly, we consider only values of M and A such that S(t) ≥ 0, i.e. 0 ≤ A ≤ M. We start by exploring three prototypical scenarios exemplified by different values of the parameter A. In particular, we choose M = 1 and compare the numerical solutions obtained for A = 0 (i.e. constant nutrient concentration), A = 0.5 (i.e. lower nutrient variability) and A = 1 (i.e. higher nutrient variability). Figure 2 displays plots of the nutrient concentration S(t), the rescaled maximum fitness g(t) and the fittest phenotypic state ϕ(t) corresponding to such choices of the parameter A. These plots show that, as one would expect, higher nutrient variability brings about more pronounced variations in the rescaled maximum fitness and the fittest phenotypic state. Figure 3 shows a comparison between the exact solutions (13)-with v i (t), μ i (t) and ρ i (t) obtained by solving numerically the Cauchy problem (14)-and the numerical solutions of the system of non-local parabolic PDEs (3) subject to the initial condition (12). In agreement with the results established by Proposition 1, for all values of A considered, there is a perfect match between the population sizes obtained by  Fig. 3). Similarly, there is excellent agreement between the population density functions obtained by solving numerically the system of PDEs (3) (cf. solid lines in the right column of Fig. 3) and the population density functions (13) with v i (t), μ i (t) and ρ i (t) given by the numerical solutions of the Cauchy problem (14) (cf. dashed and dotted lines in the right column of Fig. 3).
In accord with the results of Theorem 3, when the nutrient concentration is constant (i.e. S(t) ≡ M), the population with the lower rate of phenotypic variations (i.e. population L) outcompetes the other population (vid. Fig. 3a). The size of the surviving population ρ L (t) reaches the asymptotic value (27) and the phenotype distribution at the end of the simulations n L (x, t f ) is Gaussian with mean and variance equal to the asymptotic values (28).
In agreement with the results established by Theorem 4 (vid. Remark 3), a similar outcome is observed in the presence of a low nutrient variability (vid. Fig. 3b). In fact, in this case Λ L < Λ H . On the contrary, the population with the higher rate of phenotypic variations (i.e. population H ) outcompetes the other population when the nutrient variability is sufficiently high (vid. Fig. 3c). This is due to the fact that in this case the condition Λ H < Λ L is met. As expected (cf. Remark 2), since A > 0 both the size and the mean phenotype of the surviving population become T -periodic, with mean values given by (51) and (56), respectively. Moreover, the phenotype distribution  Table 1 with β H = 0.025. Right column. Plots of the corresponding phenotype distribution of the surviving population obtained by solving numerically the system of PDEs (3) subject to the initial condition (12). In particular, the plot of n L (x, t f ) is shown in row a, while the plots in rows b and c are, respectively, those of n L (x, t) and n H (x, t) at t = t 1 and t = t 2 , with t 1 and t 2 being highlighted in the corresponding plots in the left column. The dotted and dashed lines correspond to the exact phenotype distributions (13) with v i (t), μ i (t) and ρ i (t) given by the numerical solutions of the Cauchy problem (14) (color figure online) of the surviving population remains Gaussian with variance given by (61). This implies that if population H outcompetes population L then the variance of the phenotype distribution (i.e. the level of phenotypic heterogeneity) will be ultimately larger than in the case where population L is selected. Taken together, these results demonstrate that when the nutrient concentration is constant, or in the presence of a low level of nutrient variability, it is evolutionarily more desirable to rarely undergo spontaneous phenotypic variations, since environmental conditions are stable. Conversely, when nutrient variability is high (i.e. alternating cycles of starvation and nutrient abundance occur), higher rates of spontaneous phenotypic variations constitute a competitive advantage, as they allow for a quicker adaptation to changeable environmental conditions, and higher level of phenotypic heterogeneity emerge.
Exploiting the results of the analysis of evolutionary dynamics developed in Sect. 3, we can further assess the range of environmental conditions under which higher rates of spontaneous phenotypic variations will represent a source of competitive advan- Low mean, small amplitude High mean, small amplitude tage. In more detail, as shown by the asymptotic results established by Theorem 4, provided that condition (57) is not satisfied (i.e. at least one population survives), the outcome of competition between population H and population L in periodically fluctuating environments can be predicted by computing the value of the quantities Λ H and Λ L given by (52). In particular, the population characterised by the lower value of this quantity will ultimately be selected. Therefore, we computed Λ H and Λ L for different values of the period of the nutrient oscillations T and different values of the rate of spontaneous phenotypic variations β H . We used the values of the other evolutionary parameters reported in Table 1 and considered possible values of the environmental parameters M and A corresponding to three different scenarios: an environment whereby the nutrient is abundant and undergoes small-amplitude periodic oscillations, i.e. M is relatively large and A is relatively small (vid. Fig. 4a); an environment whereby the nutrient is scarce and undergoes small-amplitude periodic oscillations, i.e. M and A are both relatively small (vid. Fig. 4b); and an environment whereby periodic oscillations can induce a sufficiently high variability of nutrient concentration, i.e. M = A and different values of A are allowed (vid. Fig. 4c).
The results obtained are summarised by the plots in Fig. 4. As we would expect (cf. Remark 3), if the nutrient concentration undergoes low-amplitude periodic oscillations then Λ L < Λ H for all values of T and β H considered (vid. Fig. 4a, b). On the other hand, when periodic oscillations can bring about sufficiently high levels of nutrient variability, there is a region of the β H − T plane where Λ H < Λ L (vid. Fig. 4c). When the value of A is either low or high, this region is small and concentrated in the bottom left corner of the plane. For intermediate values of A the region where Λ H < Λ L is wider and such that the smaller the value T (i.e. the higher the frequency of the nutrient oscillations) the wider the range of values of β H that belong to it. Furthermore, we can investigate how the fluctuations in the nutrient level affect the phenotypic distribution of each population at any given time point by constructing numerical solutions to the system of non-local PDEs (1) subject to initial condition (12) with S(t) defined according to (79). This is demonstrated by the plots in Fig. 5, which show sample dynamics of the nutrient concentration S(t), the phenotype distributions n H (x, t) and n L (x, t), and the population sizes ρ H (x, t) and ρ L (x, t), for different values of semi-amplitude A and mean M of the fluctuations.
When the nutrient is abundant and experiences fluctuations of relatively low level (vid. Fig. 5a) population L outcompetes population H . This is due to the fact that, as shown by Fig. 8a in "Appendix", the fittest phenotypic state ϕ(t) undergoes very small periodic oscillations and its value remains close to 0 (i.e. the value of the phenotypic variable x corresponding to the fittest phenotypic state when nutrient is abundant). Moreover, the rescaled maximum fitness g(t) undergoes very small periodic oscillations and its value remains close to 1.
When the nutrient level is uniformly low and undergoes relatively small oscillations (vid. Fig. 5d) population L is selected against population H . This is due to the fact that, as shown by Fig. 8d in "Appendix", both the fittest phenotypic state ϕ(t) and the rescaled maximum fitness g(t) undergo small periodic oscillations and their values remain close to 1. We recall that x = 1 is the value of the phenotypic variable corresponding to the fittest phenotypic state when nutrient is scarce.
For the cases when the populations experience fluctuations of relatively high level (vid. Fig. 5b, c) population H outcompetes population L. This is due to the fact that, as shown by Fig. 8b, c in "Appendix", the fittest phenotypic state ϕ(t) fluctuates periodically between 1 and a positive value close to 0. Moreover, the rescaled maximum fitness g(t) undergoes small periodic oscillations and its value remains close to 1.
Note that in all cases the phenotype distribution of the surviving population remains unimodal with maximum at the mean phenotypic state. Ultimately, both the size and the mean phenotypic state of the surviving population oscillate periodically with the same period as the nutrient concentration S(t). Furthermore, when the populations experience fluctuations of relatively high level (i.e. Figure 5b, c), the mean phenotypic state of the surviving population (i.e. population H ) undergoes rapid transitions between 1 and a positive value close to 0. This can be biologically seen as the emergence of a bet-hedging behaviour.

Interpretation of the results in the context of cancer metabolism
The generality of the model and the robustness of the results make our conclusions applicable to a broad range of asexual populations evolving in fluctuating environments. As an example, in this section we discuss the biological implications of our mathematical results in the context of cancer cell metabolism and tumourmicroenvironment interactions.
Cancers begin from single cells that grow to form organ-like masses within multicellular organisms. A fundamental property of cancer cells is a self-defined fitness function such that their proliferation is determined by their heritable phenotypic properties and the local environmental selection forces. That is, individual cancer cells have the capacity to evolve novel phenotypes and to adapt to the often harsh intratumoural environment. In contrast, while normal epithelial cells have the capacity to change their phenotype to some degree, e.g. they can only do so within normal physiological constraints in response to stress. In other words, because the proliferation of normal cells is entirely governed by local tissue constraints, they cannot, unlike cancer cells, evolve adaptations to many non-physiological conditions. This difference in evolutionary capacity (or reaction norm) can confer a significant adaptive advantage to cancer cells. For example, cancer cells often metabolise glucose using glycolytic (converting glucose to lactic acid) pathways even when the oxygen concentration is sufficient for the aerobic mechanism (converting glucose to H 2 O and CO 2 ). Known as the Warburg effect (Warburg 1925), this can be understood in a Darwinian context as a form of niche construction because inefficiency of ATP (energy currency of cells) production is offset by the evolutionary advantage of generating a locally acidic environment. Cancer cells can evolve adaptive strategies to survive and proliferate in such an environment but normal cells cannot.
Angiogenesis is another form of niche construction as cancer cells, acting as a loosely organised group, produce and excrete pro-vascular proteins such as VEGF. Importantly, angiogenesis in cancers will occur entirely through these local interactions so that new vascular sprouts will emerge from the nearest vessel regardless of its flow capacity. Furthermore, cancer cells, once receiving blood flow, have no evolutionary imperative to invest resources in vascular maturation that permits blood flow to other regions of the tumour (Gillies et al. 2018). Such an unregulated vascular network will be highly unstable with cycles of growth and regression and blood flow dynamics that are inevitably disordered.
A number of studies have demonstrated that this disordered process of angiogenesis produces stochastic variations in blood flow leading to cycles of perfusion, cessation of flow, and then re-perfusion (Kimura et al. 1996). This produces corresponding fluctuations in local environmental conditions that are dependent on blood flow, including oxygen and glucose, retention of metabolites such as acid, and important signalling molecules such as testosterone or oestrogen. Particularly, regions of normoxia (normal levels of oxygen), chronic hypoxia (low oxygen level) and cycling hypoxia have been distinguished in experimental and clinical studies (Michiels et al. 2016). If we assume that S(t) in our model represents the oxygen level at time t, then the phenotypic variants best adapted to oxygen-rich environments, and thus displaying a regular metabolism, are those with x → 0, while the phenotypic variants best adapted to oxygen-low environments-i.e. the phenotypic variants that proliferate through the consumption of glucose, which is usually abundant-correspond to x → 1. For simplicity, we ignore any costs associated with the choice of metabolic preference, so that assumption (10) is satisfied. This is illustrated by the schemes in Fig. 6.
Regions of normoxia and chronic hypoxia are analogous to cases a and b in Fig. 4 where low levels of environmental variability are observed. Our results support the idea that these regions will be mainly populated by phenotypic variants best adapted to either oxygen-rich or oxygen-low environments. Moreover, since our results indicate that a higher rate of phenotypic variation does not constitute a competitive advantage in the presence of small environmental fluctuations, we expect relatively low levels of phenotypic heterogeneity to be observed in regions of either normoxia or chronic hypoxia.
On the contrary, regions of cyclic hypoxia are characterised by high variability in the oxygen levels. This can be related to case c in Fig. 4 where, under nutrient fluctuations leading to drastic environmental changes, having a higher rate of phenotypic variation represents a competitive advantage, and the mean phenotype switches between the two extreme phenotypic states (i.e. x = 0 and x = 1). Thus, in these regions we would expect to have higher levels of phenotypic heterogeneity, consistent with previous experimental findings (Chen et al. 2018), and to observe cells adopting bet-hedging as

Conclusions
We have presented a mathematical model for the evolutionary dynamics of two asexual phenotype-structured populations competing in periodically oscillating environments. The two populations undergo heritable, spontaneous phenotypic variations at different rates and their fitness landscape is dynamically sculpted by the occurrence of fluctuations in the concentration of a nutrient.
Our analytical results formalise the idea that when nutrient levels experience small and slow periodic oscillations, and thus environmental conditions are relatively stable, it is evolutionarily more efficient to rarely undergo spontaneous phenotypic variations. Conversely, under relatively large and fast periodic oscillations in the nutrient levels, which lead to alternating cycles of starvation and nutrient abundance, higher rates of spontaneous phenotypic variations can confer a competitive advantage, as they may allow for a quicker adaptation to changeable environmental conditions. In the latter case, our results indicate that higher levels of phenotypic heterogeneity are to be expected compared to those observed in slowly fluctuating environments. Finally, our results suggest that bet-hedging evolutionary strategies, whereby individuals switch between antithetical phenotypic states, can naturally emerge in the presence of relatively large and fast nutrient fluctuations leading to drastic environmental changes.
We conclude with an outlook on possible extensions of the present work. The focus of this paper has been on the case where the maximum proliferation rate of the phenotypic variants best adapted to nutrient-rich environments (i.e. the parameter γ ) and the maximum proliferation rate of the phenotypic variants best adapted to nutrient-scarce environments (i.e. the parameter ζ ) are the same. However, there are biological scenarios whereby the ability to survive in harsh environments comes with a fitness cost. For instance, cells are known to turn to glucose for their energy production when oxygen is in short supply, i.e. they produce energy through anaerobic glycolysis instead of using oxidative phosphorylation that requires aerobic conditions. Since anaerobic glycolysis is far less efficient in terms of produced energy than oxidative phosphorylation (Vander Heiden et al. 2009), the proliferation rate of glucose-dependent phenotypic variants might be lower than that of phenotypic variants best adapted to oxygenated environments. Therefore, it would be interesting to extend our analytical results to the case where ζ < γ . It would also be interesting to consider cases where an integral operator is used, in place of a differential operator, to describe the effect of phenotypic variations. In this case, we expect the qualitative behaviour of the results presented here to remain unchanged when the transition between phenotypic states is modelled via Gaussian kernels of sufficiently small variance. From a mathematical point of view, this would require further development of the methods of proof presented here in order to carry out a similar asymptotic analysis of evolutionary dynamics.
Another natural extension of the model is consideration of the feedback from the populations on the nutrient level. In fact, most existing models of evolutionary dynamics in a fluctuating environment do not account for this feedback and potentially nonlinear dynamical interactions between individuals and the surrounding environment. However, changes in the population dynamics are known to affect the outcome of interspecies competition in the presence of time variations in the availability of nutrients (Okuyama 2015). For instance, in the context of solid tumours, one could consider both negative feedback mechanisms that regulate population growth, such as nutrient consumption, and positive feedback mechanisms that promote the supply of nutrient, such as angiogenesis, i.e. hypoxia-induced formation of blood vessels. Consideration of these mechanisms is expected to affect the advantages gained by each population, and therefore the phenotypic composition of the tumour, in a given environment.
An additional development of our study would be to incorporate into the model a spatial structure, as done for instance in Bouin and Calvez (2014), Lorenzi et al. (2018), Lorz et al. (2015), Alfaro et al. (2013),  and Lam and Lou (2014), and let multiple nutrient sources with different inflows be distributed across the spatial domain. This would lead individuals to experience nutrient fluctuations of variable amplitudes and frequencies depending on their spatial position, thus leading to the emergence of multiple local niches whereby different phenotypic variants could be selected. Such an extension of our study would be relevant in several biological contexts. In fact, spatial niche partitioning is known to have an important impact on interspecies competition, as it promotes the coexistence between species best adapted to different local environmental conditions (Hassell et al. 1994;Kneitel and Chase 2004). Moreover, in the context of cancer research, clinical images and histological data have revealed the existence of considerable levels of spatial heterogeneity in oxygen distribution within tumours (Tomaszewski et al. 2017;Michiels et al. 2016), and localised regions of cycling hypoxia and chronic hypoxia have been identified in tumour xenografts (Matsuo et al. 2014). In this regard, a spatially-stuctured version of our model could shed new light on the ways in which the interplay between spatial and temporal variability of oxygen levels may dictate the phenotypic composition and the level of phenotypic heterogeneity of malignant tumours.