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

Aleksandra Arda\v{s}eva; Robert A. Gatenby; Alexander R. A. Anderson,; Helen M. Byrne; Philip K. Maini; Tommaso Lorenzi

arXiv:1905.11712·q-bio.PE·April 3, 2020

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

Aleksandra Arda\v{s}eva, Robert A. Gatenby, Alexander R. A. Anderson,, Helen M. Byrne, Philip K. Maini, Tommaso Lorenzi

PDF

TL;DR

This paper models the evolutionary dynamics of two competing populations with phenotypic variations in environments with periodic nutrient fluctuations, revealing how oscillation rates influence phenotypic variation strategies.

Contribution

It introduces a novel non-local PDE model for competing phenotype-structured populations in fluctuating environments, analyzing long-term evolutionary behavior.

Findings

01

Slow nutrient oscillations favor rare phenotypic variations.

02

Fast oscillations favor higher phenotypic variation rates.

03

Results have implications for understanding cancer metabolism.

Abstract

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 spontaneous phenotypic variations at different rates. The phenotypic state of each individual is represented by a continuous variable, and the…

Tables1

Table 1. Table 1: Parameter values used to carry out numerical simulations.

Parameter	Description	Value / Value range
$γ$	Maximum proliferation rate	100
$d$	Death rate due to competition	0.01
$β_{L}$	Rate of phenotypic variation of population $L$	0.01
$β_{H}$	Rate of phenotypic variation of population $H$	[0.01, 0.1]

Equations197

ρ_{H} (t) = \int_{R} n_{H} (x, t) d x, ρ_{L} (t) = \int_{R} n_{L} (x, t) d x, ρ (t) = ρ_{H} (t) + ρ_{L} (t) .

ρ_{H} (t) = \int_{R} n_{H} (x, t) d x, ρ_{L} (t) = \int_{R} n_{L} (x, t) d x, ρ (t) = ρ_{H} (t) + ρ_{L} (t) .

μ_{i} (t) = \frac{1}{ρ _{i} ( t )} \int_{R} x n_{i} (x, t) d x, σ_{i}^{2} (t) = \frac{1}{ρ _{i} ( t )} \int_{R} x^{2} n_{i} (x, t) d x - μ_{i}^{2} (t) .

μ_{i} (t) = \frac{1}{ρ _{i} ( t )} \int_{R} x n_{i} (x, t) d x, σ_{i}^{2} (t) = \frac{1}{ρ _{i} ( t )} \int_{R} x^{2} n_{i} (x, t) d x - μ_{i}^{2} (t) .

\left\{\begin{array}[]{ll}\displaystyle{\frac{\partial n_{H}}{\partial t}=\beta_{H}\frac{\partial^{2}n_{H}}{\partial x^{2}}\;+\;R\big{(}x,S(t),\rho(t)\big{)}\,n_{H}},\\ \\ \displaystyle{\frac{\partial n_{L}}{\partial t}=\beta_{L}\frac{\partial^{2}n_{L}}{\partial x^{2}}\;\;\;+\;R\big{(}x,S(t),\rho(t)\big{)}\,n_{L}},\end{array}\right.\quad\text{for }(x,t)\in\mathbb{R}\times(0,\infty).

\left\{\begin{array}[]{ll}\displaystyle{\frac{\partial n_{H}}{\partial t}=\beta_{H}\frac{\partial^{2}n_{H}}{\partial x^{2}}\;+\;R\big{(}x,S(t),\rho(t)\big{)}\,n_{H}},\\ \\ \displaystyle{\frac{\partial n_{L}}{\partial t}=\beta_{L}\frac{\partial^{2}n_{L}}{\partial x^{2}}\;\;\;+\;R\big{(}x,S(t),\rho(t)\big{)}\,n_{L}},\end{array}\right.\quad\text{for }(x,t)\in\mathbb{R}\times(0,\infty).

β_{H} > β_{L} .

β_{H} > β_{L} .

R\big{(}x,S,\rho\big{)}=p(x,S)\;-\;d\rho.

R\big{(}x,S,\rho\big{)}=p(x,S)\;-\;d\rho.

p\big{(}x,S(t)\big{)}=\gamma\,\frac{S(t)}{1+S(t)}\,\left(1-x^{2}\right)\;+\;\zeta\,\left(1-\frac{S(t)}{1+S(t)}\right)\left[1-\left(1-x\right)^{2}\right],

p\big{(}x,S(t)\big{)}=\gamma\,\frac{S(t)}{1+S(t)}\,\left(1-x^{2}\right)\;+\;\zeta\,\left(1-\frac{S(t)}{1+S(t)}\right)\left[1-\left(1-x\right)^{2}\right],

p\big{(}x,S\big{)}=\gamma\,g(S)-h(S)\left(x-\varphi(S)\right)^{2}

p\big{(}x,S\big{)}=\gamma\,g(S)-h(S)\left(x-\varphi(S)\right)^{2}

g (S) = \frac{1}{1 + S} (S + \frac{ζ}{γ} \frac{ζ}{ζ + γ S}), φ (S) = \frac{ζ}{ζ + γ S}

g (S) = \frac{1}{1 + S} (S + \frac{ζ}{γ} \frac{ζ}{ζ + γ S}), φ (S) = \frac{ζ}{ζ + γ S}

h (S) = ζ + (γ - ζ) \frac{S}{1 + S} .

h (S) = ζ + (γ - ζ) \frac{S}{1 + S} .

ζ = γ .

ζ = γ .

g (S) = \frac{1}{1 + S} (S + \frac{1}{1 + S}), φ (S) = \frac{1}{1 + S} and h (S) \equiv γ .

g (S) = \frac{1}{1 + S} (S + \frac{1}{1 + S}), φ (S) = \frac{1}{1 + S} and h (S) \equiv γ .

g (t) \equiv g (S (t)) and φ (t) \equiv φ (S (t)) .

g (t) \equiv g (S (t)) and φ (t) \equiv φ (S (t)) .

n_{i} (x, 0) = ρ_{i}^{0} \frac{v _{i}^{0}}{2 π} exp [- \frac{v _{i}^{0}}{2} (x - μ_{i}^{0})^{2}] for i \in {H, L},

n_{i} (x, 0) = ρ_{i}^{0} \frac{v _{i}^{0}}{2 π} exp [- \frac{v _{i}^{0}}{2} (x - μ_{i}^{0})^{2}] for i \in {H, L},

n_{i} (x, t) = ρ_{i} (t) \frac{v _{i} ( t )}{2 π} exp [- \frac{v _{i} ( t )}{2} (x - μ_{i} (t))^{2}] for i \in {H, L},

n_{i} (x, t) = ρ_{i} (t) \frac{v _{i} ( t )}{2 π} exp [- \frac{v _{i} ( t )}{2} (x - μ_{i} (t))^{2}] for i \in {H, L},

\left\{\begin{array}[]{ll}\displaystyle{v_{i}^{\prime}(t)=2\left(\gamma-\beta_{i}v_{i}^{2}(t)\right)},\\ \\ \displaystyle{\mu_{i}^{\prime}(t)=\frac{2\gamma}{v_{i}(t)}\left(\varphi(t)-\mu_{i}(t)\right)},\\ \\ \displaystyle{\rho_{i}^{\prime}(t)=\left(F_{i}(t)-d\rho(t)\right)\rho_{i}(t)},\\ \\ v_{i}(0)=v^{0}_{i},\quad\mu_{i}(0)=\mu^{0}_{i},\quad\rho_{i}(0)=\rho^{0}_{i},\\ \\ \rho(t)=\rho_{H}(t)+\rho_{L}(t),\end{array}\right.\qquad\text{for }\;i\in\left\{H,L\right\},

\left\{\begin{array}[]{ll}\displaystyle{v_{i}^{\prime}(t)=2\left(\gamma-\beta_{i}v_{i}^{2}(t)\right)},\\ \\ \displaystyle{\mu_{i}^{\prime}(t)=\frac{2\gamma}{v_{i}(t)}\left(\varphi(t)-\mu_{i}(t)\right)},\\ \\ \displaystyle{\rho_{i}^{\prime}(t)=\left(F_{i}(t)-d\rho(t)\right)\rho_{i}(t)},\\ \\ v_{i}(0)=v^{0}_{i},\quad\mu_{i}(0)=\mu^{0}_{i},\quad\rho_{i}(0)=\rho^{0}_{i},\\ \\ \rho(t)=\rho_{H}(t)+\rho_{L}(t),\end{array}\right.\qquad\text{for }\;i\in\left\{H,L\right\},

F_{i} (t) \equiv F_{i} (t, v_{i} (t), μ_{i} (t)) = γ g (t) - \frac{γ}{v _{i} ( t )} - γ (μ_{i} (t) - φ (t))^{2} .

F_{i} (t) \equiv F_{i} (t, v_{i} (t), μ_{i} (t)) = γ g (t) - \frac{γ}{v _{i} ( t )} - γ (μ_{i} (t) - φ (t))^{2} .

\frac{\partial n _{i}}{\partial t} = β_{i} \frac{\partial ^{2} n _{i}}{\partial x ^{2}} + [γ g (t) - γ (x - φ (t))^{2} - d ρ (t)] n_{i}, n_{i} \equiv n_{i} (x, t) .

\frac{\partial n _{i}}{\partial t} = β_{i} \frac{\partial ^{2} n _{i}}{\partial x ^{2}} + [γ g (t) - γ (x - φ (t))^{2} - d ρ (t)] n_{i}, n_{i} \equiv n_{i} (x, t) .

\frac{ρ _{i}^{'}}{ρ _{i}} + \frac{v _{i}^{'}}{2 v _{i}}

\frac{ρ _{i}^{'}}{ρ _{i}} + \frac{v _{i}^{'}}{2 v _{i}}

v_{i}^{'} + 2 β_{i} v_{i}^{2} = 2 γ .

v_{i}^{'} + 2 β_{i} v_{i}^{2} = 2 γ .

μ_{i}^{'} = \frac{2 γ ( φ - μ _{i} )}{v _{i}} .

μ_{i}^{'} = \frac{2 γ ( φ - μ _{i} )}{v _{i}} .

\frac{ρ _{i}^{'}}{ρ _{i}} + \frac{v _{i}^{'}}{2 v _{i}} = - β_{i} v_{i} + γ g - γ (μ_{i} - φ)^{2} - d ρ

\frac{ρ _{i}^{'}}{ρ _{i}} + \frac{v _{i}^{'}}{2 v _{i}} = - β_{i} v_{i} + γ g - γ (μ_{i} - φ)^{2} - d ρ

ρ_{i}^{'} = (F_{i} - d ρ) ρ_{i},

ρ_{i}^{'} = (F_{i} - d ρ) ρ_{i},

v_{i} (0) = v_{i}^{0}, μ_{i} (0) = μ_{i}^{0} and ρ_{i} (0) = ρ_{i}^{0} .

v_{i} (0) = v_{i}^{0}, μ_{i} (0) = μ_{i}^{0} and ρ_{i} (0) = ρ_{i}^{0} .

S (t) \equiv \overline{S} \geq 0,

S (t) \equiv \overline{S} \geq 0,

g (t) \equiv \overline{g} and φ (t) \equiv \overline{φ} .

g (t) \equiv \overline{g} and φ (t) \equiv \overline{φ} .

β_{L} \geq γ \overline{g}

β_{L} \geq γ \overline{g}

t \to \infty lim ρ_{H} (t) = 0 and t \to \infty lim ρ_{L} (t) = 0;

t \to \infty lim ρ_{H} (t) = 0 and t \to \infty lim ρ_{L} (t) = 0;

β_{L} < γ \overline{g}

β_{L} < γ \overline{g}

t \to \infty lim ρ_{H} (t) = 0, t \to \infty lim ρ_{L} (t) = \frac{γ}{d} (γ \overline{g} - β_{L})

t \to \infty lim ρ_{H} (t) = 0, t \to \infty lim ρ_{L} (t) = \frac{γ}{d} (γ \overline{g} - β_{L})

t \to \infty lim μ_{L} (t) = \overline{φ}, t \to \infty lim σ_{L}^{2} (t) = \frac{β _{L}}{γ} .

t \to \infty lim μ_{L} (t) = \overline{φ}, t \to \infty lim σ_{L}^{2} (t) = \frac{β _{L}}{γ} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

∎

11institutetext: Aleksandra Ardaševa 22institutetext: Wolfson Centre for Mathematical Biology, Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, UK

22email: [email protected] 33institutetext: Robert A Gatenby44institutetext: Department of Integrated Mathematical Oncology, H. Lee Moffitt Cancer Center, Tampa, Florida, USA

44email: [email protected] 55institutetext: Alexander R A Anderson66institutetext: Department of Integrated Mathematical Oncology, H. Lee Moffitt Cancer Center, Tampa, Florida, USA

66email: [email protected] 77institutetext: Helen M Byrne 88institutetext: Wolfson Centre for Mathematical Biology, Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, UK

88email: [email protected] 99institutetext: Philip K Maini 1010institutetext: Wolfson Centre for Mathematical Biology, Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, UK

1010email: [email protected] 1111institutetext: Tommaso Lorenzi 1212institutetext: School of Mathematics and Statistics, University of St Andrews, St Andrews, KY16 9SS, UK

1212email: [email protected]

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

Aleksandra Ardaševa

Robert A Gatenby

Alexander R A Anderson

Helen M Byrne

Philip K Maini

Tommaso Lorenzi

(Received: date / Accepted: date)

Abstract

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.

Keywords:

Periodically fluctuating environments Evolutionary dynamics Spontaneous phenotypic variation Bet-hedging Non-local partial differential equations

1 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 xue2018benefits . 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 cohen1966optimizing . 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 philippi1989hedging . Phenotypic heterogeneity, a characteristic feature of a risk spreading strategy, is observed in many systems, including bacterial populations kussell2005phenotypic and solid tumours marusyk2012intra .

Bet-hedging is typically proposed to occur within the context of bacterial populations, where experimental support for stochastic phenotype switching is available kussell2005phenotypic ; smits2006phenotypic ; veening2008bistability ; beaumont2009experimental ; acar2008stochastic . 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 balaban2004bacterial . Schreiber et al. schreiber2016phenotypic 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 kussell2005phenotypic . It is also hypothesized to be present in cancer where irregular vasculature can cause significant fluctuations within the tumour microenvironment gillies2018eco . 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 cairns2001acute ; robertson2015impact ; chen2018intermittent ; nichol2016stochasticity .

Mathematically, competition between populations evolving in fluctuating environments has been studied using different modelling approaches, including deterministic predator-prey models and stochastic models cushing1980two ; chesson1994multispecies ; roxburgh2004intermediate ; anderies2000fluctuating . 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 fudenberg2012phenotype ; muller2013bet ; ashcroft2014fixation ; wienand2017evolution ; gravenmier2018adaptation . Muller et al. muller2013bet 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 acar2008stochastic ; muller2013bet ; wienand2017evolution . 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 non-local 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 alfaro2018evolutionary ; alfaro2017effect ; almeida2019evolution ; bouin2014travelling ; bouin2012invasion ; busse2016mass ; calsina2004small ; calsina2007asymptotic ; calsina2013asymptotics ; chisholm2016evolutionary ; chisholm2015emergence ; cuadrado2009equilibria ; delitala2012mathematical ; delitala2012asymptotic ; delitala2013mathematical ; desvillettes2008selection ; diekmann2005dynamics ; domschke2017structured ; iglesias2018long ; lam2017dirac ; lam2017stability ; lorenzi2014asymptotic ; lorenzi2016tracking ; lorenzi2015dissecting ; lorenzi2015mathematical ; lorenzi2018role ; lorz2015modeling ; lorz2013populational ; mirrahimi2015time ; Nordmann2018 ; perthame2006transport ; pouchol2018asymptotic ; turanova2015model . 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 iglesias2018long ; lorenzi2015dissecting . 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 lorenzi2015dissecting 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 Section 2 we introduce the mathematical model. In Section 3 we carry out an analytical study of evolutionary dynamics. In Section 4 we integrate the analytical results with numerical simulations. In Section 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.

2 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\in\mathbb{R}$ , and we describe the phenotype distributions of the two populations at time $t\in[0,\infty)$ by means of the population density functions $n_{H}(x,t)~{}\geq~{}0$ and $n_{L}(x,t)~{}\geq~{}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

[TABLE]

Moreover, we define, respectively, the mean phenotypic state and the related variance of each population $i\in\{H,L\}$ at time $t$ as

[TABLE]

In the mathematical framework of our model, the function $\sigma^{2}_{i}(t)$ provides a measure of the level of phenotypic heterogeneity in the $i^{th}$ population. Finally, we introduce a function $S(t)\geq 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

[TABLE]

In the system of PDEs (3), the diffusion terms model the effects of spontaneous phenotypic variations, which occur at rates $\beta_{H}>0$ and $\beta_{L}>0$ , with

[TABLE]

The functional $R\big{(}x,S(t),\rho(t)\big{)}$ 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 $\rho(t)$ – i.e. the functional $R\big{(}x,S(t),\rho(t)\big{)}$ can be seen as the phenotypic fitness landscape of the two populations at time $t$ . Throughout the paper, we define this fitness functional as

[TABLE]

Definition (5) translates into mathematical terms the following biological ideas: (i) all else being equal, individuals die due to nterpopulation and intrapopulation competition at rate $d\rho(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.

Assumption 2.1

Phenotypic variants with $x\to 0$ have a competitive advantage over the other phenotypic variants when the nutrient concentration is high (i.e. if $S(t)\gg 1$ ).

Assumption 2.2

Phenotypic variants with $x\to 1$ are favoured over the other phenotypic variants when the nutrient concentration is low (i.e. if $S(t)\ll 1$ ).

Under Assumptions 2.1 and 2.2, we define the net proliferation rate as

[TABLE]

with $0<\zeta\leq\gamma$ . The parameters $\gamma$ and $\zeta$ 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 otwinowski2014inferring . In fact, after a little algebra, definition (6) can be rewritten as

[TABLE]

with

[TABLE]

and

[TABLE]

Under the environmental conditions defined by the nutrient concentration $S$ , the function $0\leq\varphi(S)\leq 1$ represents the fittest phenotypic state, $\gamma\,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 2.1 and 2.2, equation (7) shows that definition (6) is such that the fittest phenotypic state $\varphi(S)$ belongs to the interval $[0,1]$ for any nutrient concentration $S\geq 0$ , i.e. $\varphi:\mathbb{R}_{\geq 0}\to[0,1]$ . In particular, under starvation conditions (i.e. if $S=0$ ) the fittest phenotypic state is $\varphi(0)=1$ , while increasing nutrient concentrations correspond to values of the fittest phenotypic state closer to [math], i.e. $\varphi^{\prime}(S)<0$ for all $S\geq 0$ and $\varphi(S)\to 0$ as $S\to\infty$ . Furthermore, the fact that the function $p(x,S)$ is negative for values of $x$ sufficiently far from the fittest phenotypic state $\varphi(S)$ captures the idea that less fit variants are driven to extinction by natural selection. These observations are illustrated by the plots in Figure 1.

Henceforth for simplicity we assume

[TABLE]

Under assumption (10), definitions (8) and (9) become, respectively,

[TABLE]

Moreover, since we assume the function $S(t)$ to be given, we use the notation

[TABLE]

3 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

[TABLE]

where $\rho^{0}_{i}\in\mathbb{R}_{>0}$ , $v^{0}_{i}\in\mathbb{R}_{>0}$ and $\mu^{0}_{i}\in\mathbb{R}$ .

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 rice2004evolutionary .

Before turning to the case of periodically fluctuating environments in Section 3.2, we consider the case of constant environments in Section 3.1. The proofs of the results presented in Section 3.1 and Section 3.2 rely on the results established by the 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

[TABLE]

with the population size, $\rho_{i}(t)$ , the mean phenotypic state, $\mu_{i}(t)$ , and the inverse of the related variance, $v_{i}(t)=1/\sigma^{2}_{i}(t)$ , being solutions of the Cauchy problem

[TABLE]

where

[TABLE]

Proof

Substituting definitions (5), (7) and (11) into the non-local PDE (3) for $n_{i}(x,t)$ yields

[TABLE]

Building upon the results presented in almeida2019evolution ; chisholm2016evolutionary ; lorenzi2015dissecting , we make the ansatz (13) and substituting this ansatz into equation (16) we find

[TABLE]

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 for $v_{i}$ alone

[TABLE]

Moreover, equating the coefficients of the first-order terms in $x$ , and eliminating $v^{\prime}_{i}$ from the resulting equation, yields

[TABLE]

Lastly, choosing $x=\mu_{i}$ in equation (17) gives

[TABLE]

and eliminating $v^{\prime}_{i}$ from equation (20) we find

[TABLE]

with the function $F_{i}(t)$ being defined according to (15). Under the initial condition (12), we have

[TABLE]

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)$ , $\mu_{i}(t)$ and $\rho_{i}(t)$ . ∎

3.1 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

[TABLE]

which implies that

[TABLE]

In this case, our main results are summarised by Theorem 3.1.

Theorem 3.1

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:

(i)

if

[TABLE]

then

[TABLE]

(ii)

if

[TABLE]

then

[TABLE]

and

[TABLE]

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, $\rho_{i}(t)$ , the mean phenotypic state, $\mu_{i}(t)$ , and the inverse of the related variance, $v_{i}(t)=1/\sigma_{i}^{2}(t)$ , being governed by the Cauchy problem (14) with $g(t)\equiv\overline{g}$ and $\varphi(t)\equiv\overline{\varphi}$ .

In this framework, we divide the proof of Theorem 3.1 into four steps. We study the asymptotic behaviour of $v_{i}(t)$ , $\mu_{i}(t)$ and $F_{i}(t)$ for $t\to\infty$ (Step 1). We show that $\rho_{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)$ , $\mu_{i}(t)$ and $F_{i}(t)$ for $t\to\infty$ .

Solving the separable first-order differential equation (14)1 for $v_{i}(t)$ and imposing the initial condition (14)4 gives

[TABLE]

which implies that

[TABLE]

Moreover, solving the differential equation (14)2 for $\mu_{i}(t)$ by the integrating factor method and imposing the initial condition (14)4 yields

[TABLE]

from which, using the positivity of $v_{i}(t)$ , we find that

[TABLE]

Lastly, noting that, under the additional assumption (22), the function $F_{i}(t)$ defined by (15) reads as

[TABLE]

the asymptotic results (30) and (32) allow us to conclude that

[TABLE]

Step 2: non-negativity and boundedness of $\rho_{i}(t)$ .

Solving the differential equation (14)3 for $\rho_{i}$ and imposing the initial condition (14)4 yields

[TABLE]

This result, along with the positivity of $\rho^{0}_{i}$ , implies that

[TABLE]

Moreover, substituting (33) into the differential equation (14)3 for $\rho_{i}$ yields

[TABLE]

with $j=L$ if $i=H$ and $j=H$ if $i=L$ . Estimating from above the right-hand side of the latter differential equation by using the non-negativity of $\rho_{j}(t)$ [cf. the uniform lower bound (36)], the positivity of $v_{i}(t)$ [cf. expression (29)] and the fact that $\overline{g}<2$ [cf. definition (11)], we obtain the differential inequality

[TABLE]

which gives the uniform upper bound

[TABLE]

Step 3: proof of claim (25).

Combining the asymptotic result (34) with the expression (35) for $\rho_{i}$ we find that

[TABLE]

for some positive constant $C$ . Since the function $\rho(t)$ is non-negative [cf. the uniform lower bound (36)], the asymptotic relation (38) ensures that

[TABLE]

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 $\rho_{H}(t)>0$ , we can compute the quotient of $\rho_{L}(t)$ and $\rho_{H}(t)$ through (35). In so doing we find

[TABLE]

Using the limit (34) for $F_{i}$ , we then have

[TABLE]

for some positive constant $C$ . Under assumption (4), the asymptotic relation (41) gives

[TABLE]

and, since $\rho_{L}$ is uniformly bounded from above [cf. the uniform upper bound (37)], from (41) we conclude that

[TABLE]

We can rewrite the differential equation (14)3 for $\rho_{L}$ as

[TABLE]

where the function $\eta(t)$ is defined as

[TABLE]

Using the asymptotic results (30), (32) and (42), we see that

[TABLE]

Solving the differential equation (43) complemented with the initial condition $\rho_{L}(0)=\rho_{L}^{0}$ yields chisholm2016evolutionary

[TABLE]

The result (44) ensures that in the asymptotic regime $t\to\infty$ we have

[TABLE]

and, under the additional assumption (26), we also have

[TABLE]

for some positive constant $C$ . These asymptotic relations, along with the expression (45) for $\rho_{L}$ , allow us to conclude that

[TABLE]

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.1 provide a mathematical formalisation of the idea that in constant environments:

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. 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. 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.

3.2 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,\infty)\to\mathbb{R}_{\geq 0}$ satisfy the assumptions

[TABLE]

which implies that the functions $g(t)$ and $\varphi(t)$ satisfy the assumptions

[TABLE]

Our main results are summarised by Theorem 3.2, the proof of which relies on the results established by Lemma 1 and Lemma 2.

Lemma 1

Under assumptions (5), (7), (11) and (47), the unique real $T$ -periodic solution of the problem

[TABLE]

is

[TABLE]

and satisfies the integral identity

[TABLE]

Lemma 2

Let

[TABLE]

and

[TABLE]

with $u_{i}(t)$ given by (1). Under assumptions (5), (7), (11), (47) and the additional assumption

[TABLE]

the unique real non-negative $T$ -periodic solution of the problem

[TABLE]

is

[TABLE]

and satisfies the integral identity

[TABLE]

We refer the interested reader to lorenzi2015dissecting for the proofs of Lemma 1 and Lemma 2.

Theorem 3.2

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:

(i)

if

[TABLE]

then

[TABLE]

(ii)

if

[TABLE]

and

[TABLE]

then

[TABLE]

and

[TABLE]

with $w_{i}(t)$ and $u_{i}(t)$ given by (55) and (1), respectively.

Proof

Proposition 1 ensures that the population density function $n_{i}(x,t)$ is of the Gaussian form (13) with the population size, $\rho_{i}(t)$ , the mean phenotypic state, $\mu_{i}(t)$ , and the inverse of the related variance, $v_{i}(t)=1/\sigma_{i}^{2}(t)$ , being governed by the Cauchy problem (14). In this framework, we prove Theorem 3.2 in 4 steps. In Step 1, we study the asymptotic behaviour of $v_{i}(t)$ , $\mu_{i}(t)$ and $F_{i}(t)$ for $t\to\infty$ . In Step 2, we show that $\rho_{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)$ , $\mu_{i}(t)$ and $F_{i}(t)$ for $t\to\infty$ .

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.1 still holds and

[TABLE]

Moreover, using the asymptotic result (62) along with the linear differential equation (14)2 for $\mu_{i}(t)$ one can easily show that

[TABLE]

where $u_{i}(t)$ is a $T$ -periodic solution of the differential equation (49). Lemma 1 ensures that $u_{i}(t)$ is given by (1). Lastly, for $F_{i}(t)$ defined according to (15), the asymptotic results (62) and (63) allow us to conclude that

[TABLE]

Step 2: non-negativity and boundedness of $\rho_{i}(t)$ .

Proceeding in a similar way as in the proof of Theorem 3.1 (cf. Step 2 in the proof of Theorem 3.1), one can prove that

[TABLE]

Step 3: proof of claim (58).

Solving the differential equation (14)3 for $\rho_{i}$ and imposing the initial condition (14)4 yields

[TABLE]

with $F_{i}(t)$ defined according to (15). Combining the asymptotic result (64) with the expression (66) for $\rho_{i}(t)$ gives

[TABLE]

for some positive constant $C$ . Hence, using the fact that the functions $g(t)$ , $\varphi(t)$ and $u_{i}(t)$ are $T$ -periodic and considering ${\rm m}\to\infty$ , we find

[TABLE]

for some positive constant $C$ . Since the function $\rho(t)$ is non-negative [cf. the uniform lower bound (65)], the asymptotic relation (68) ensures that if

[TABLE]

then

[TABLE]

This proves that if assumption (57) is satisfied then claim (58) is verified.

Step 4: proof of claims (60) and (61).

Let

[TABLE]

As long as $\rho_{j}(t)>0$ , we can compute the quotient of $\rho_{i}(t)$ and $\rho_{j}(t)$ through (66). In so doing, using the asymptotic relation (67) for $\rho_{i}(t)$ and $\rho_{j}(t)$ and considering ${\rm m}\to\infty$ we obtain

[TABLE]

for some positive constant $C$ , with $\Lambda_{i}$ and $\Lambda_{j}$ defined according to (52). The asymptotic relation (71) allows us to conclude that

[TABLE]

Since $\rho_{i}$ is uniformly bounded from above [cf. the uniform upper bound (65)], the asymptotic result (72) implies that

[TABLE]

We can rewrite the differential equation (14)3 for $\rho_{i}$ as

[TABLE]

where the function $\eta(t)$ is defined as

[TABLE]

Using the asymptotic results (62), (63) and (73) we see that $\eta(t)\to 0$ exponentially fast as $t\to\infty$ . Hence, $\rho_{i}(t)\to\tilde{\rho}_{i}(t)$ as $t\to\infty$ , with $\tilde{\rho}_{i}(t)$ being the solution of the differential equation

[TABLE]

with

[TABLE]

subject to an initial condition $0<\tilde{\rho}_{i}(0)<\infty$ . We note that: (i) the function $\tilde{\rho}_{i}(t)$ is uniformly bounded as it satisfies the upper and lower bounds

[TABLE]

(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$ , $\varphi$ and $u_{i}$ are $T$ -periodic Lipschitz continuous functions of $t$ [cf. assumptions (48) and expression (1)]. Therefore, the conditions of Massera’s Convergence Theorem massera1950existence ; smith1986massera are satisfied and this allows us to conclude that

[TABLE]

with $w_{i}(t)$ being a non-negative $T$ -periodic solution of the differential equation (54). Under the additional assumption (59), that is,

[TABLE]

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) with $i=\underset{k\in\left\{H,L\right\}}{\operatorname{arg}\,\operatorname{min}}\;\Lambda_{k}$ . ∎

In summary, Theorem 3.2 gives an explicit characterisation of the long-term limit of $v_{i}(t)$ , $\mu_{i}(t)$ and $\rho_{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

[TABLE]

where

[TABLE]

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 3.2 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

[TABLE]

Integrating both sides of the above equation with respect to $t$ between [math] and $T$ , and using the fact that $u_{i}(T)-\varphi(T)=u_{i}(0)-\varphi(0)$ , yields

[TABLE]

Therefore, definition (52) can be rewritten as

[TABLE]

If $S\equiv\bar{S}$ then $g(t)\equiv\overline{g}$ and $\varphi(t)\equiv\overline{\varphi}$ (i.e. $\varphi^{\prime}\equiv 0$ ). In this case,

[TABLE]

and (78) allows one to see that $\Lambda_{i}=\sqrt{\beta_{i}}$ . Hence, if $S$ is constant then the results of Theorem 3.2 reduce to the results of Theorem 3.1. Moreover, the first term in the expression (78) for $\Lambda_{i}$ is clearly a monotonically increasing function of $\beta_{i}$ , whereas the factor in front of the integral in the second term is a monotonically decreasing function of $\beta_{i}$ . Hence, if the mean value of the $T-$ periodic function $\left(\varphi(t)-u_{i}(t)\right)\varphi^{\prime}(t)$ is sufficiently small then $\Lambda_{H}>\Lambda_{L}$ , while if such a mean value is sufficiently large then $\Lambda_{H}<\Lambda_{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 $\varphi(t)$ .

The asymptotic results established by Theorem 3.2 formalise mathematically the idea that in periodically fluctuating environments:

populations undergoing spontaneous phenotypic variations at a rate too large compared to the mean value of the maximum fitness will ultimately go extinct [cf. point (i)]; 2. 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. 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]; 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]; 4. 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.

4 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 Section 4.1, we describe the numerical methods employed and the set-up of numerical simulations. In Section 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.

4.1 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\in[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 $\Delta 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 leveque2007finite . On the other hand, we use the Matlab built-in solver ode45 to solve numerically the Cauchy problem (14) for $v_{i}(t)$ , $\mu_{i}(t)$ and $\rho_{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 $\beta_{i}\ll\gamma$ for $i\in\left\{H,L\right\}$ . Furthermore, given the values of $\gamma$ and $\beta_{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 voorde2019improving .

We remark that the value of the parameter $\beta_{L}$ and the range of values of the parameter $\beta_{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_{i}^{0}=20$ , $\mu_{i}^{0}=0$ and $\rho_{i}^{0}\approx 800$ for $i\in\left\{H,L\right\}$ .

4.2 Main results

We consider the following definition of the nutrient concentration

[TABLE]

In definition (79), the parameter $M>0$ represents the mean nutrient concentration, while the parameter $A\geq 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)\geq 0$ , i.e. $0\leq A\leq 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 $\varphi(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)$ , $\mu_{i}(t)$ and $\rho_{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 computing numerically the integrals of the components of the numerical solution of the system of PDEs (3) (cf. solid lines in the left column of Figure 3) and the population sizes obtained by solving numerically the Cauchy problem (14) (cf. dashed and dotted lines in the left column of Figure 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 Figure 3) and the population density functions (13) with $v_{i}(t)$ , $\mu_{i}(t)$ and $\rho_{i}(t)$ given by the numerical solutions of the Cauchy problem (14) (cf. dashed and dotted lines in the right column of Figure 3).

In accord with the results of Theorem 3.1, when the nutrient concentration is constant (i.e. $S(t)\equiv M$ ), the population with the lower rate of phenotypic variations (i.e. population $L$ ) outcompetes the other population (vid. Figure 3A). The size of the surviving population $\rho_{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 3.2 (vid. Remark 3), a similar outcome is observed in the presence of a low nutrient variability (vid. Figure 3B). In fact, in this case $\Lambda_{L}<\Lambda_{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. Figure 3C). This is due to the fact that in this case the condition $\Lambda_{H}<\Lambda_{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 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 Section 3, we can further assess the range of environmental conditions under which higher rates of spontaneous phenotypic variations will represent a source of competitive advantage. In more detail, as shown by the asymptotic results established by Theorem 3.2, 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 $\Lambda_{H}$ and $\Lambda_{L}$ given by (52). In particular, the population characterised by the lower value of this quantity will ultimately be selected. Therefore, we computed $\Lambda_{H}$ and $\Lambda_{L}$ for different values of the period of the nutrient oscillations $T$ and different values of the rate of spontaneous phenotypic variations $\beta_{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. Figure 4A); an environment whereby the nutrient is scarce and undergoes small-amplitude periodic oscillations, i.e. $M$ and $A$ are both relatively small (vid. Figure 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. Figure 4C).

The results obtained are summarised by the plots in Figure 4. As we would expect (cf. Remark 3), if the nutrient concentration undergoes low-amplitude periodic oscillations then $\Lambda_{L}<\Lambda_{H}$ for all values of $T$ and $\beta_{H}$ considered (vid. Figure 4A and Figure 4B). On the other hand, when periodic oscillations can bring about sufficiently high levels of nutrient variability, there is a region of the $\beta_{H}-T$ plane where $\Lambda_{H}<\Lambda_{L}$ (vid. Figure 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 $\Lambda_{H}<\Lambda_{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 $\beta_{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 Figure 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 $\rho_{H}(x,t)$ and $\rho_{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. Figure 5A) population $L$ outcompetes population $H$ . This is due to the fact that, as shown by Figure 8A in Appendix A, the fittest phenotypic state $\varphi(t)$ undergoes very small periodic oscillations and its value remains close to [math] (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. Figure 5D) population $L$ is selected against population $H$ . This is due to the fact that, as shown by Figure 8D in Appendix A, both the fittest phenotypic state $\varphi(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. Figure 5B and Figure 5C) population $H$ outcompetes population $L$ . This is due to the fact that, as shown by Figure 8B and Figure 8C in Appendix A, the fittest phenotypic state $\varphi(t)$ fluctuates periodically between $1$ and a positive value close to [math]. 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 and Figure 5C), the mean phenotypic state of the surviving population (i.e. population $H$ ) undergoes rapid transitions between $1$ and a positive value close to [math]. This can be biologically seen as the emergence of a bet-hedging behaviour.

5 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 tumour-microenvironment 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 H2O and CO2). Known as the Warburg effect warburg1925metabolism , 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 gillies2018eco . 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 kimura1996fluctuations . 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 michiels2016cycling . 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\rightarrow 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\rightarrow 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 Figure 6.

Regions of normoxia and chronic hypoxia are analogous to cases A and B in Figure 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 Figure 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 chen2018intermittent , and to observe cells adopting bet-hedging as a survival strategy in response to drastic variation in oxygen levels. These conclusions are summarised in the table provided in Figure 7.

6 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 $\gamma$ ) and the maximum proliferation rate of the phenotypic variants best adapted to nutrient-scarce environments (i.e. the parameter $\zeta$ ) 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 vander2009understanding , 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 $\zeta<\gamma$ . 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 okuyama2015demographic . 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 bouin2014travelling ; lorenzi2018role ; lorz2015modeling ; alfaro2013travelling ; mirrahimi2015asymptotic ; lam2014evolution , 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 hassell1994species ; kneitel2004trade . 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 tomaszewski2017oxygen ; michiels2016cycling , and localised regions of cycling hypoxia and chronic hypoxia have been identified in tumour xenografts matsuo2014magnetic . 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.

Acknowledgements.

AA is supported by funding from the Engineering and Physical Sciences Research Council (EPSRC) and the Medical Research Council (MRC) (grant no. EP/L016044/1) and in part by the Moffitt Cancer Center PSOC, NIH/NCI (grant no. U54CA193489). RG and ARAA are supported by Physical Sciences Oncology Network (PSON) grant from the National Cancer Institute (grant no. U54CA193489) as well as the Cancer Systems Biology Consortium grant from the National Cancer Institute (grant no. U01CA23238). ARAA and RG would also like to acknowledge support from the Moffitt Cancer Center of Excellence for Evolutionary Therapy.

Appendix A Supplementary figures

Bibliography76

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) B. Xue, S. Leibler, Proceedings of the National Academy of Sciences 115 (50), 12745 (2018)
2(2) D. Cohen, Journal of Theoretical Biology 12 (1), 119 (1966)
3(3) T. Philippi, J. Seger, Trends in Ecology & Evolution 4 (2), 41 (1989)
4(4) E. Kussell, S. Leibler, Science 309 (5743), 2075 (2005)
5(5) A. Marusyk, V. Almendro, K. Polyak, Nature Reviews Cancer 12 (5), 323 (2012)
6(6) W.K. Smits, O.P. Kuipers, J.W. Veening, Nature Reviews Microbiology 4 (4), 259 (2006)
7(7) J.W. Veening, W.K. Smits, O.P. Kuipers, Annual Review of Microbiology 62 , 193 (2008)
8(8) H.J. Beaumont, J. Gallie, C. Kost, G.C. Ferguson, P.B. Rainey, Nature 462 (7269), 90 (2009)

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

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

Abstract

Keywords:

1 Introduction

2 Model description

Assumption 2.1

Assumption 2.2

3 Analysis of evolutionary dynamics

Remark 1

Proposition 1

Proof

3.1 Evolutionary dynamics in constant environments

Theorem 3.1

Proof

Step 1: asymptotic behaviour of vi(t)v_{i}(t)vi​(t), μi(t)\mu_{i}(t)μi​(t) and Fi(t)F_{i}(t)Fi​(t) for t→∞t\to\inftyt→∞.

Step 2: non-negativity and boundedness of ρi(t)\rho_{i}(t)ρi​(t).

Step 3: proof of claim (25).

Step 4: proof of claims (27) and (28).

3.2 Evolutionary dynamics in periodically fluctuating environments

Lemma 1

Lemma 2

Theorem 3.2

Proof

Step 1: asymptotic behaviour of vi(t)v_{i}(t)vi​(t), μi(t)\mu_{i}(t)μi​(t) and Fi(t)F_{i}(t)Fi​(t) for t→∞t\to\inftyt→∞.

Step 2: non-negativity and boundedness of ρi(t)\rho_{i}(t)ρi​(t).

Step 3: proof of claim (58).

Step 4: proof of claims (60) and (61).

Remark 2

Remark 3

4 Numerical simulations

4.1 Numerical methods and set-up of numerical simulations

4.2 Main results

5 Interpretation of the results in the context of cancer metabolism

6 Conclusions

Acknowledgements.

Appendix A Supplementary figures

Step 1: asymptotic behaviour of $v_{i}(t)$ , $\mu_{i}(t)$ and $F_{i}(t)$ for $t\to\infty$ .

Step 2: non-negativity and boundedness of $\rho_{i}(t)$ .

Step 1: asymptotic behaviour of $v_{i}(t)$ , $\mu_{i}(t)$ and $F_{i}(t)$ for $t\to\infty$ .

Step 2: non-negativity and boundedness of $\rho_{i}(t)$ .