Evolutionary branching via replicatormutator equations
Abstract.
We consider a class of nonlocal reactiondiffusion problems, referred to as replicatormutator equations in evolutionary genetics. For a confining fitness function, we prove wellposedness and write the solution explicitly, via some underlying Schrödinger spectral elements (for which we provide new and nonstandard 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 eigenelements2010 Mathematics Subject Classification:
92B05, 92D15, 35K15, 45K051. Introduction
In this paper we first study the existence, uniqueness and long time behaviour of solutions , , , to the integrodifferential 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 unimodal to multimodal shape, thus reproducing evolutionary branching.
The above equation is referred to as a replicatormutator 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 meanfield theory for the evolution of RNA virus population. Without mutations, and under the constraint of constant mass , the dynamics is given by
(2) 
with in [40]. In this context, represents the density of a population (at time and per unit of phenotypic trait) on a onedimensional phenotypic trait space. The function represents the fitness of the phenotype and models the individual reproductive success; thus the nonlocal term
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 replicatormutator problem (1) from individual based models.
Equation (1) is supplemented with a nonnegative and bounded, initial data such that , so that, formally, for later times. Indeed, integrating formally (1) over , the total mass
solves the initial value problem
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 AvronHerbst 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 nonlocal 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 welladapted. 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 “cutoff 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 wellposedness, 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 nonstandard estimates on the eigenelements. Also, from a modelling perspective, this enables to reproduce evolutionary branching, consisting of the spontaneous splitting from unimodal to multimodal distribution of the trait.
Such splitting phenomena have long been discussed and analysed in different frameworks, see e.g. [29] via HamiltonJacobi technics, [42] within finite populations, or [27] for a LotkaVolterra system in a bounded domain. In a replicatormutator 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 nonlocal diffusion (mutation kernel), namely
Their approach [30], [18] is based on Cumulant Generating Functions (CGF): it turns out that the CGF satisfies a first order nonlocal 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 wellpossessedness 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
Proposition 2.1 (Spectral basis).
Let satisfy Assumption 1. Then the operator
(3) 
is essentially selfadjoint on , and has discrete spectrum: there exists an orthonormal basis of consisting of eigenfunctions of
with corresponding eigenvalues
of finite multiplicity.
Remark 1.
In the quantum mechanics terminology, is known as the ground state, corresponding to the boundstate 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
(4) 
where is the energy functional given by
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 as ).
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, nonnegative, and compactly supported in function with . We define the test function
From the variational formula (4), we have
The first integral in the right hand side is given by
as . The second integral gives
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 :
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 squareintegrable [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
(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 CourantFisher’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 socalled pseudopolynomials (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. , and weighted 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 nonstandard estimates on the eigenfunctions.
We define
the mass of the th eigenfunction of the Hamiltonian . In the sequel, by
we mean that there is such that, for all , .
Proposition 2.4 ( norm of eigenfunctions).
Let satisfy Assumption 2. Then we have
(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 ,
(7) 
Proof.
Remark 3.
The correct power in (7) can be retrieved by a standard homogeneity argument. Indeed, defining for , we get
so that
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
and integrating over , we get
Integrating by parts and recalling that eigenfunctions are normalized in , we obtain
so that
Next, it follows from Assumption 2 (and ) that there is such that for all , and thus
Now, by Lemma 2.5, we have
which, combined with (5), implies (6). The proposition is proved. ∎
Proposition 2.6 ( norm of eigenfunctions).
Let satisfy Assumption 2. Then we have
(9) 
Proof.
Proposition 2.7 (Weighted norm of eigenfunctions).
Let satisfy Assumption 2. Then we have
Proof.
From Assumption 2 and , we can find large enough so that the following facts hold for all : there are and such that
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
are respectively super and subsolutions of the eigenvalue equation
so that
(10) 
by the comparison principle. An analogous estimate holds on .
3. Wellposedness 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 nonlocal 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
Recall that the Sobolev space is defined as
where the derivative is understood in the distributional sense. We denote by the Hilbert space with inner product defined by
and with usual inner product
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 RieszFréchetKolmogorov theorem, see e.g. [8, Theorem 4.26]. Let a bounded sequence of functions of : there is such that, for all ,
We first need to show the uniform smallness of the tails of . Let . Select large enough so that for all . Then
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 ,
so that
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

[label=()]

For all , all ,
where the time derivative is understood in the distributional sense. Equivalently, for all , all ,
(11) 
is a continuous function on .

.
Here is our main mathematical result.
Theorem 3.4 (Solving replicatormutator problem).
Proof.
We proceed by necessary and sufficient condition. Let be a solution, in the sense of Definition 3.3. We define the function as
(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
since .
We now show that solves the linear Cauchy problem
(14) 
Indeed, formally for the moment,
so that
since solves (1). Those computations can be made rigorous in the distributional sense. Indeed, for a test function , set
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 wellposedness of the linear Cauchy problem (14) is postponed to the next subsection: from Proposition 3.6, we know that, for all ,
Now, the estimates on the eigenvalues and the norm of eigenfunctions, namely Proposition 2.3 and Proposition 2.6, allow to write
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
(15) 
On the other hand, we claim that, for all ,
(16) 
which follows formally by integrating (14) over . To prove (16) rigorously, notice first that by Proposition 2.3 and 2.4, the series
converges for all . Hence , the total mass of , is given by
Next, for any , thanks to Proposition 2.3 and 2.4, so that is differentiable on and
the last equality following by similar arguments based on Proposition 2.7. Hence (16) is proved. From (15), (16) and , we deduce that
for all . As a conclusion, (13) is inverted into
(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 ,