Evolutionary branching

# Evolutionary branching via replicator-mutator equations

Matthieu Alfaro and Mario Veruete IMAG, Université de Montpellier, CC051, 34095, Montpellier, France
###### Abstract.

We consider a class of non-local reaction-diffusion problems, referred to as replicator-mutator equations in evolutionary genetics. For a confining fitness function, we prove well-posedness and write the solution explicitly, via some underlying Schrödinger spectral elements (for which we provide new and non-standard estimates). As a consequence, the long time behaviour is determined by the principal eigenfunction or ground state. Based on this, we discuss (rigorously and via numerical explorations) the conditions on the fitness function and the mutation rate for evolutionary branching to occur.

###### Key words and phrases:
Evolutionary genetics, dynamics of adaptation, branching phenomena, long time behaviour, Schrödinger eigenelements
###### 2010 Mathematics Subject Classification:
92B05, 92D15, 35K15, 45K05

## 1. Introduction

In this paper we first study the existence, uniqueness and long time behaviour of solutions , , , to the integro-differential Cauchy problem

 (1)

which serves as a model for the dynamics of adaptation, and where is a confining fitness function (see below for details). Next, we enquire on the possibility, depending on the function and the parameter , for a solution to split from uni-modal to multi-modal shape, thus reproducing evolutionary branching.

The above equation is referred to as a replicator-mutator model. This type of model has found applications in different fields such as economics and biology [25], [4]. In the field of evolutionary genetics, a free spatial version of equation (1) was introduced by Tsimring, Levine and Kessler in [40], where they propose a mean-field theory for the evolution of RNA virus population. Without mutations, and under the constraint of constant mass , the dynamics is given by

 ∂u∂t=u(W(x)−∫RW(y)u(t,y)dy), (2)

with in [40]. In this context, represents the density of a population (at time and per unit of phenotypic trait) on a one-dimensional phenotypic trait space. The function represents the fitness of the phenotype and models the individual reproductive success; thus the non-local term

 ¯¯¯u(t):=∫RW(y)u(t,y)dy

stands for the mean fitness at time .

As a first step to take into account evolutionary phenomena, mutations are modelled by the local diffusion operator , where is the mutation rate, so that equation (2) is transferred into (1). We refer to the recent paper [41] for a rigorous derivation of the replicator-mutator problem (1) from individual based models.

Equation (1) is supplemented with a non-negative and bounded, initial data such that , so that, formally, for later times. Indeed, integrating formally (1) over , the total mass

 m(t):=∫Ru(t,x)dx

solves the initial value problem

 ddtm(t)=(1−m(t))¯¯¯u(t), m(0)=1.

Hence, by Gronwall’s lemma, , as long as is meaningful.

The case of linear fitness function, , was the first introduced in [40], but little was known concerning existence and behaviours of solutions. Let us here mention the main result of Biktashev [5]: for compactly supported initial data, solutions converge, as goes to infinity, to a Gaussian profile, where the convergence is understood in terms of the moments of . In a recent paper [2], Alfaro and Carles proved that, thanks to a tricky change of unknown based on the Avron-Herbst formula (coming from quantum mechanics), equation (1) can be reduced to the heat equation. This enables to compute the solution explicitly and describe contrasted behaviours depending on the tails of the initial datum: either the solution is global and tends, as tends to infinity, to a Gaussian profile which is centred around (acceleration) and is flattening (extinction in infinite horizon), or the solution becomes extinct in finite time (or even immediately) thus contradicting the conservation of the mass, previously formally observed.

For quadratic fitness functions, , it turns out that the equation can again be reduced to the heat equation [3], up to an additional use of the generalized lens transform of the Schrödinger equation. In the case , for any initial data, there is extinction at a finite time which is always bounded from above by . Roughly speaking, both the right and left tails quickly enlarge, so that, in order to conserve the mass, the central part is quickly decreasing. Then the non-local mean fitness term becomes infinite very quickly and equation (1) becomes meaningless (extinction). On the other hand, when , for any initial data, the solution is global and tends, as tends to infinity, to an universal stationary Gaussian profile.

The aforementioned cases and share the property of being unbounded from above, meaning that some phenotypes are infinitely well-adapted. This unlimited growth rate of in (1) yields rich mathematical behaviours (acceleration, extinction) but is not admissible for biological applications. To deal with such a problem, for the linear fitness case, some works consider a “cut-off version” of (1) at large [40], [35], [37], or provide a proper stochastic treatment for large phenotypic trait region [34].

On the other hand, is referred to as a confining fitness function, typically preventing extinction phenomena. However, it does not suffice to take into account more realistic cases for which fitness functions are defined by a linear combination of two components (e.g. birth and death rates), each maximized by different optimal values of the underlying trait, a typical case being .

Our main goal is thus to provide a rigorous treatment of the Cauchy problem (1) when the fitness function is confining. For a relatively large class of such fitness functions, we prove well-posedness, and show that the solution of (1) converges to the principal eigenfunction (or ground state) of the underlying Schrödinger operator divided by its mass. This requires rather non-standard estimates on the eigenelements. Also, from a modelling perspective, this enables to reproduce evolutionary branching, consisting of the spontaneous splitting from uni-modal to multi-modal distribution of the trait.

Such splitting phenomena have long been discussed and analysed in different frameworks, see e.g. [29] via Hamilton-Jacobi technics, [42] within finite populations, or [27] for a Lotka-Volterra system in a bounded domain. In a replicator-mutator context, let us notice that, while branching in (1) is mainly induced by the fitness function, it was recently obtained in [18] through different means. Precisely, the authors study the case of linear fitness but non-local diffusion (mutation kernel), namely

 ∂tu=J∗u−u+u(x−∫Ryu(t,y)dy).

Their approach [30], [18] is based on Cumulant Generating Functions (CGF): it turns out that the CGF satisfies a first order non-local partial differential equation that can be explicitly solved, thus giving access to many informations such as mean trait, variance, position of the leading edge. When a purely deleterious mutation kernel balances the infinite growth rate of , they reveal some branching scenarios.

The paper is organized as follows. In Section 2 we present the underlying linear material. In Section 3 we prove the well-possessedness of the Cauchy problem associated to (1). We also provide an explicit expression of the solution and studies its long time behaviour. In Section 4 we discuss, through rigorous details or numerical explorations, the conditions on the shape of the fitness function and on the mutation parameter for branching phenomenon to occur. Finally, we briefly conclude in Section 5.

## 2. Some spectral properties

In this section, we present some linear material. We first quote some very classical results [39], [33], [1], [21], [22], [20], [15] for Schrödinger operators, and then prove less standard estimates on the eigenfunctions, which are crucial for later analysis.

### 2.1. Confining fitness functions and eigenvalues properties

Confining fitness functions tend to at infinity. In quantum mechanics, this corresponds to potentials describing the evolution of quantum particles subject to an external field force that prevents them from escaping to infinity, that is, particles have a high probability of presence in a bounded spatial region.

###### Assumption 1 (Confining fitness function).

The fitness function is continuous and confining, that is

 lim|x|→∞W(x)=−∞.
###### Proposition 2.1 (Spectral basis).

Let satisfy Assumption 1. Then the operator

 H:=−σ2d2dx2−W(x) (3)

is essentially self-adjoint on , and has discrete spectrum: there exists an orthonormal basis of consisting of eigenfunctions of

 Hϕk=λkϕk,∥ϕk∥L2(R)=1,

with corresponding eigenvalues

 λ0<λ1≤λ2≤⋯≤λk→+∞,

of finite multiplicity.

###### Remark 1.

Proposition 2.1 is a classical result [33], [38, Chapter 3, Theorem 1.4 and Theorem 1.6] for Schrödinger operators with confining potential . In this paper, the minus sign is due to the biological interpretation of the fitness function .

In the quantum mechanics terminology, is known as the ground state, corresponding to the bound-state of minimal energy . In this paper we refer to the couple indistinctly as ground state/ground state energy or as principal eigenfunction/principal eigenvalue.

The principal eigenvalue can be characterised by the variational formulation

 λ0=inf{E(u):u∈C∞c(R),∥u∥L2(R)=1}, (4)

where is the energy functional given by

 E(u)=σ2∫R∣∣∣∂u∂x∣∣∣2dx+∫R−W(x)|u(x)|2dx.

Using concentrated test functions, the above formula enables to understand the behaviour of the principal eigenvalue as the mutation rate tends to 0. The following will be used in Section 4 to prove some branching phenomena.

###### Proposition 2.2 (Asymptotics for λ0(σ) as σ→0).

Let satisfy Assumption 1. Assume that reaches a global maximum at . Then .

###### Proof.

For the convenience of the reader, we give the proof of this standard fact. Let be a smooth, non-negative, and compactly supported in function with . We define the test function

 pσ(x):=1σ1/4p(x−α√σ).

From the variational formula (4), we have

 −M≤λ0(σ)≤σ2∫R|∂xpσ(x)|2dx+∫R−W(x)|pσ(x)|2dx.

The first integral in the right hand side is given by

 σ2∫R|∂xpσ(x)|2dx=σ∥∥p′∥∥L2(R)→0,

as . The second integral gives

 ∫R−W(x)|pσ(x)|2dx=∫R−W(α+√σy)p(y)2dy

which, by the -dominated convergence theorem tends to as . ∎

In the subsequent sections, we will quote results on the spectral properties of Schrödinger operators, in particular an asymptotics for the eigenvalues as . As far as we know, the available results require to assume that the fitness is polynomial.

###### Assumption 2 (Polynomial confining fitness function).

The fitness function is a real polynomial of degree :

 W(x)=−x2s+2s−1∑k=0wkxk,

for some integer and some real numbers , .

Under Assumption 2, elliptic regularity theory insures that the eigenfunctions are infinitely differentiable. Furthermore, all the derivatives of each eigenfunction are square-integrable [17]. Notice that it is also known that all eigenfunctions actually belong to the Schwartz space .

###### Proposition 2.3 (Asymptotics for eigenvalues).

Let satisfy Assumption 2. Then all eigenvalues of are simple and

 λk∼Cs,σk2ss+1 as k→+∞, (5)

where , with being the gamma function.

We refer to [39], [15] and the references therein for more details on the above asymptotic formula. Furthermore, in the case of a symmetric fitness , the simplicity of eigenvalues enforce all eigenfunctions to be even or odd. In particular the principal eigenfunction (ground state) is even since it is known to have constant sign.

###### Remark 2.

Assume that is such that for some polynomial as in Assumption 2 and some constant . From Courant-Fisher’s theorem, that is the variational characterization of the eigenvalues, we deduce that , where are the eigenvalues of the Hamiltonian with potential . Hence, share with the asymptotics (5), which is the keystone for deriving the estimates on eigenfunctions in subsection 2.2, and thereafter our main results in Section 3. Hence, our results apply to such fitness functions, covering in particular the case of the so-called pseudo-polynomials (i.e. smooth functions which coincide, outside of a compact region, with a polynomial as in Assumption 2), which are relevant for numerical computations.

### 2.2. L1, L∞ and weighted L1 norms of the eigenfunctions

In the study of spectral properties of Schrödinger operators, efforts tend to concentrate around asymptotic estimates of eigenvalues or on the regularity and decay of eigenfunctions [39], [1], [10], [16]. Much less attention has been given to estimate the and norms of eigenfunctions. One reason is that the natural framework for eigenfunctions of the Hamiltonian , defined in (3), is . On the other hand, the biological nature of problem (1) suggests and as natural spaces for the solution . We therefore provide in this subsection rather non-standard estimates on the eigenfunctions.

We define

 mk:=∫Rϕk(x)dx,

the mass of the -th eigenfunction of the Hamiltonian . In the sequel, by

 Ak≲Bk

we mean that there is such that, for all , .

###### Proposition 2.4 (L1 norm of eigenfunctions).

Let satisfy Assumption 2. Then we have

 |mk|≤∥ϕk∥L1(R)≲k12(s+1). (6)

Before proving the above proposition, we need the following lemma which is of independent interest.

###### Lemma 2.5.

Let and be given. Then there is a constant such that, for all ,

 ∥f∥L1(Rd)≤C ∥f∥1−δL2(Rd)∥∥xN/2f∥∥δL2(Rd),δ:=d/N. (7)
###### Proof.

Let denote the open -dimensional ball of radius and center . We write

 ∥f∥L1(Rd)=∫BR|f(x)|dx+∫Rd∖BR1|x|N/2|x|N/2|f(x)|dx=:I1+I2.

By the Cauchy-Schwarz inequality we have

and

 I2 ≤ ≤ C(∫+∞R1rNrd−1dr)1/2∥∥xN/2f∥∥L2(Rd) ≤ C2R(d−N)/2∥∥xN/2f∥∥L2(Rd),

for some and . Summarizing,

 (8)

Now, we select

 R=⎛⎜⎝C2(N−d)C1d∥∥xN/2f∥∥L2(Rd)∥f∥L2(Rd)⎞⎟⎠2/N

which minimizes the right hand side of (8) and yields (7). ∎

###### Remark 3.

The correct power in (7) can be retrieved by a standard homogeneity argument. Indeed, defining for , we get

 ∥fλ∥L1(Rd)=1λd∥f∥L1(Rd),
 ∥fλ∥L2(Rd)=1λd/2∥f∥L2(Rd),
 ∥∥xN/2fλ∥∥L2(Rd)=1λ(N+d)/2∥∥xN/2f∥∥L2(Rd),

so that

 1λd≲(1λd/2)1−δ(1λ(N+d)/2)δ.

Powers of in both sides must coincide, which enforces .

We can now estimate the mass of the eigenfunctions.

###### Proof of Proposition 2.4.

Up to subtracting a constant to , we can assume without loss of generality that . Multiplying by the eigenvalue equation

 −σ2ϕ′′k−Wϕk=λkϕk

and integrating over , we get

 ∫R−σ2ϕ′′k(x)ϕk(x)dx+∫R−W(x)ϕk(x)2dx=λk∫Rϕ2k(x)dx.

Integrating by parts and recalling that eigenfunctions are normalized in , we obtain

 σ2∫Rϕ′k(x)2dx+∫R−W(x)ϕk(x)2dx=λk,

so that

 ∫R−W(x)ϕk(x)2dx≤λk.

Next, it follows from Assumption 2 (and ) that there is such that for all , and thus

 ∥xsϕk∥2L2(R)≤λkγ.

Now, by Lemma 2.5, we have

 ∥ϕk∥L1(R)≲∥xsϕk∥1/(2s)L2(R)≲λ1/4sk,

which, combined with (5), implies (6). The proposition is proved. ∎

###### Proposition 2.6 (L∞ norm of eigenfunctions).

Let satisfy Assumption 2. Then we have

 ∥ϕk∥L∞(R)≲ks2(s+1). (9)
###### Proof.

Since , we have

 ∥ϕk∥2L∞(R)≲∥ϕk∥L2(R)∥∥ϕ′k∥∥L2(R)≲λ1/2k,

and the conclusion follows from (5). ∎

###### Proposition 2.7 (Weighted L1 norm of eigenfunctions).

Let satisfy Assumption 2. Then we have

 ∥Wϕk∥L1(R)≲k5s+22s+2.
###### Proof.

From Assumption 2 and , we can find large enough so that the following facts hold for all : there are and such that

 W(x)+λk≥0 ∀x∈(−yk,xk) W(x)+λk=0 ∀x∈{−yk,xk} W(x)+λk≤0 ∀x∈R∖(−yk,xk),

and is decreasing on . Assumption 2 implies that and thus, from Proposition 2.3, . Next, up to enlarging if necessary, it follows from Assumption 2 that for all . As a result, functions

 ϕ±(x):=±∥ϕk∥∞e−(x−2xk)

are respectively super and sub-solutions of the eigenvalue equation

 −ϕ′′k−(W(x)+λk)ϕk=0,

so that

 |ϕk(x)|≤∥ϕk∥∞e−(x−2xk)∀x∈(2xk,+∞), (10)

by the comparison principle. An analogous estimate holds on .

In order to estimate , we split the domain of integration into three parts: , and . Setting

 Ii:=∫Ωi|W(x)ϕk(x)|dx=∫Ωi−W(x)|ϕk(x)|dx,

we decompose . Notice that,

 ∫+∞2xk−W(x)e−(x−2xk)dx=[−P(x)e−(x−2xk)]+∞2xk=P(2xk),

where is a polynomial of same degree as . Hence, from (10) and Proposition 2.6, we get

 I1 ≲ ∥ϕk∥∞∫+∞2xk|W(x)|e−(x−2xk)dx≲ks2(s+1)P(2xk) ≲ ks2(s+1)x2sk≲ks2(s+1)k2ss+1=k5s2s+2.

By monotonicity of on ,

 I2≲xk|W(2xk)|∥ϕk∥∞≲x1+2skks2(s+1)≲k5s+22s+2.

Remember that in so that

 I3≲xkλk∥ϕk∥∞≲xkk2ss+1ks2(s+1)≲k5s+22s+2.

Finally, . ∎

## 3. Well-posedness and long time behaviour

In this section we show that the Cauchy problem (1) has a unique smooth solution which is global in time. Keystones are the change of variable (13) that links the non-local equation (1) to a linear parabolic problem, and our previous estimates on the underlying eigenelements. Equipped with the representation (12) of the solution, we then prove convergence in any , , to the principal eigenfunction normalized by its mass.

Up to subtracting a constant to the confining fitness function , we can assume without loss of generality (recall the mass conservation property) that .

### 3.1. Functional framework

For a negative confining fitness function (see Assumption 1), we set

 L2−W(R):={v:R→R,∥v∥2L2−W(R):=∫R−W(x)v2(x)dx<+∞}.

Recall that the Sobolev space is defined as

 W1,2(R)={f:R→R,f∈L2(R),f′∈L2(R)},

where the derivative is understood in the distributional sense. We denote by the Hilbert space with inner product defined by

 (v,v)V:=∫Rdvdx(x)dvdx(x)dx+∫R−W(x)v(x)v(x)dx,

and with usual inner product

 (v,v)H=∫Rv(x)v(x)dx.

By Assumption 1, it is straightforward that , so that . Moreover, the following holds.

###### Lemma 3.1.

The embedding is dense, continuous and compact.

###### Proof.

This is very classical but, for the convenience of the reader, we present the details. Since and is dense in , it follows that is dense in . Next, since for ,

the embedding is continuous.

The proof of compactness follows by the Riesz-Fréchet-Kolmogorov theorem, see e.g. [8, Theorem 4.26]. Let a bounded sequence of functions of : there is such that, for all ,

 ∥vn∥2V=∫Rv′n(x)2dx+∫R−W(x)vn(x)2dx

We first need to show the uniform smallness of the tails of . Let . Select large enough so that for all . Then

 ∥vn∥2L2(R∖[−a,a]) =∫−a−∞vn(x)2dx+∫+∞avn(x)2dx ≤ε∫−a−∞−W(x)vn(x)2dx+ε∫+∞a−W(x)vn(x)2dx ≤Mε.

Next, for a compact set , we need to show the uniform smallness of the norm of . Let . By Morrey’s theorem, there is such that, for all and ,

 |vn(x+h)−vn(x)|≤C|h|1/2∥v′n∥L2(R)≤C|h|1/2M1/2,

so that

 ∥vn(⋅+h)−vn∥2L2(K)≤C2M|K||h|≤ε,

for all , if is sufficiently small. The lemma is proved. ∎

### 3.2. Main results

We first define the notion of solution to the Cauchy problem (1).

###### Definition 3.2 (Admissible initial data).

We say that a function is an admissible initial data if , and .

###### Definition 3.3 (Solution of the Cauchy problem (1)).

Let be an admissible initial data. We say that is a (global) solution of the Cauchy problem (1) if, for any , , , and

1. [label=()]

2. For all , all ,

 ddt∫Ru(t,x)v(x)dx+σ2∫R∂u∂x(t,x)dvdx(x)dx+∫R−W(x)u(t,x)v(x)dx+¯¯¯u(t)∫Ru(t,x)v(x)dx=0,

where the time derivative is understood in the distributional sense. Equivalently, for all , all ,

 (11)
3. is a continuous function on .

4. .

Here is our main mathematical result.

###### Theorem 3.4 (Solving replicator-mutator problem).

Let satisfy Assumption 2. For any admissible initial condition , there is a unique solution to the Cauchy problem (1), in the sense of Definition 3.3. Moreover the solution is smooth on and is given by

 u(t,x)=+∞∑k=0(u0,ϕk)L2(R)ϕk(x)e−λkt+∞∑k=0(u0,ϕk)L2(R)mke−λkt,t>0,x∈R, (12)

where are the eigenelements defined in Proposition 2.1, and

 mk=∫Rϕk(x)dx.
###### Proof.

We proceed by necessary and sufficient condition. Let be a solution, in the sense of Definition 3.3. We define the function as

 v(t,x):=u(t,x)exp(∫t0¯¯¯u(τ)dτ),0≤t≤T,x∈R. (13)

This function is well defined since by Definition 3.3 , the integral in the exponential is finite for all . Since and , it is straightforward to see that, for all , . Additionally, from and , one can see that . Last, due to

 ∫T0∥v(t,⋅)∥2Vdt

since .

We now show that solves the linear Cauchy problem

 {∂tv=σ2∂2xv+W(x)vv(0,x)=u0(x). (14)

Indeed, formally for the moment,

 ∂tv=(∂tu)e∫t0¯¯¯u(τ)dτ+u¯¯¯u∂2xv=(∂2xu)e∫t0¯¯¯u(τ)dτ,

so that

 ∂tv−σ2∂2xv−W(x)v=(∂tu+u¯¯¯u−σ2∂2xu−W(x)u)e∫t0¯¯¯u(τ)dτ=0

since solves (1). Those computations can be made rigorous in the distributional sense. Indeed, for a test function , set

 φ(t):=ψ(t)e∫t0¯¯¯u(τ)dτ,

and by Definition 3.3 , belongs to . Writing (11) with as test function yields the weak formulation of (14) with as test function, that is

for all .

The well-posedness of the linear Cauchy problem (14) is postponed to the next subsection: from Proposition 3.6, we know that, for all ,

 v(t)=+∞∑k=0(u0,ϕk)L2(R)ϕke−λkt in L2(R).

Now, the estimates on the eigenvalues and the norm of eigenfunctions, namely Proposition 2.3 and Proposition 2.6, allow to write

 v(t,x)=+∞∑k=0(u0,ϕk)L2(R)ϕk(x)e−λkt,0

Also, we know from the parabolic regularity theory and the comparison principle, that and that for all , .

Now, we show that the change of variable (13) can be inverted. For , multiplying (13) by and integrating over , we get

 ¯¯¯v(t):=∫RW(x)v(t,x)dx =¯¯¯u(t)exp(∫t0¯¯¯u(τ)dτ) =ddt(exp(∫t0¯¯¯u(τ)dτ)). (15)

On the other hand, we claim that, for all ,

 ddtmv(t)=¯¯¯v(t), (16)

which follows formally by integrating (14) over . To prove (16) rigorously, notice first that by Proposition 2.3 and 2.4, the series

 +∞∑k=0∣∣(u0,ϕk)L2(R)∣∣e−λkt∫R|ϕk(x)|dx

converges for all . Hence , the total mass of , is given by

 mv(t)=+∞∑k=0(u0,ϕk)L2(R)mke−λkt.

Next, for any , thanks to Proposition 2.3 and 2.4, so that is differentiable on and

 ddtmv(t)=+∞∑k=0(u0,ϕk)L2(R)mk(−λk)e−λkt=¯¯¯v(t),

the last equality following by similar arguments based on Proposition 2.7. Hence (16) is proved. From (15), (16) and , we deduce that

 exp(∫t0¯¯¯u(τ)dτ)=mv(t),

for all . As a conclusion, (13) is inverted into

 u(t,x)=v(t,x)mv(t),mv(t):=∫Rv(t,x)dx>0, (17)

for all , .

Conversely, we need to show that the function given by (17) is the solution of (1) in the sense of Definition 3.3. Let .

Since , the function is continuous on , which shows item of Definition 3.3.

Next, since and ,

 ∫T0∣∣¯¯¯u(t)∣∣dt=−∫T0ddtmv(t)mv(t)dt=−ln