Density dependent replicator-mutator models in directed evolution
DDENSITY DEPENDENT REPLICATOR-MUTATOR MODELS INDIRECTED EVOLUTION
MATTHIEU ALFARO
AND
MARIO VERUETE imag , Universit´e de Montpellier , cnrs , Montpellier, France. Abstract.
We analyze a replicator-mutator model arising in the context of directed evolu-tion [23], where the selection term is modulated over time by the mean-fitness. We combinea Cumulant Generating Function approach [13] and a spatio-temporal rescaling related tothe Avron-Herbst formula [1] to give of a complete picture of the Cauchy problem. Besidesits well-posedness, we provide an implicit/explicit expression of the solution, and analyzeits large time behaviour. As a by product, we also solve a replicator-mutator model wherethe mutation coefficient is socially determined, in the sense that it is modulated by themean-fitness. The latter model reveals concentration or anti diffusion/diffusion phenomena.
Contents
1. Introduction 22. Main results 53. Special solutions 63.1. Steady states 73.2. Gaussian solutions 74. The
CGF approach 85. The solution implicitly/explicitly 116. When the mutation coefficient is socially determined 136.1. Gaussian solutions: concentration vs. extinction 136.2. General solutions 16References 18
E-mail address : [email protected], [email protected] . Date : January 24, 2019. a r X i v : . [ m a t h . A P ] J a n DENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION Introduction
In this paper, we analyze a mathematical model of a directed evolution process which isdensity-dependent. The model in consideration is derived in [23] (see below for details) andis given by the following nonlocal equation(1.1) ∂u∂t ( t, x ) = σ ∂ u∂x ( t, x ) + u ( t, x ) x − x ( t ) x ( t ) , t > , x ∈ R , where the nonlocal term is given by(1.2) x ( t ) := (cid:90) R x u ( t, x ) dx, and will also be denoted u ( t ), to remind its dependence on the solution itself. We will provethe well-posedness of the Cauchy problem associated with (1.1), obtain the solution throughan implicit/explicit expression, and analyze its long time behaviour.As a by product of our analysis of (1.1), we will collect results on the well-posedness andthe long time behaviour of the solution to the Cauchy problem associated with equation(1.3) ∂v∂t ( t, x ) = σ x ( t ) ∂ v∂x ( t, x ) + v ( t, x )( x − x ( t )) , t > , x ∈ R , where x ( t ) := (cid:82) R xv ( t, x ) dx will also be denoted v ( t ).Throughout this work, we assume that the initial data u and v , are non-negative andhave unitary mass(1.4) (cid:90) R u ( x ) dx = (cid:90) R v ( x ) dx = 1 , so that, formally , (cid:82) R u ( t, x ) dx = (cid:82) R v ( t, x ) dx = 1 for later times. Indeed, if we formallyintegrate (1.1) over x ∈ R , we see that the total mass M ( t ) := (cid:82) R u ( t, x ) dx solves the Cauchyproblem(1.5) M (cid:48) ( t ) = 1 − M ( t ) , M (0) = 1 , so that, by the Cauchy-Lipschitz theorem, M ( t ) ≡ t >
0. Hence, u ( t, · ) is a probabilitydensity on R and x ( t ) = u ( t ) is its mean. As far as (1.3) is concerned, we reach(1.6) M (cid:48) ( t ) = v ( t )(1 − M ( t )) , M (0) = 1 , so that, by Gronwall’s lemma, M ( t ) = 1 as long as “ v ( t ) is meaningful ”. We refer to [1] forsituations where, in some sense, v ( t ) blows up in finite time, thus leading to extinction of thesolution, which contradicts the formal conservation of the mass.Let us now comment on the emergence of equations (1.1) and (1.3), termed as replicator-mutator models, in evolutionary biology and their relevance in biotechnology.Replicator-mutator models aim at describing Darwinian evolutionary processes, whosefundamental elements are mutations and selection. Originally introduced by Taylor andJonker [22], the replicator dynamics was developed for symmetric games in order to describethe evolution (determined by the payoff matrix) of the frequencies of strategies in a popula-tion, see [16] for a complete derivation. Nevertheless, such an approach neglects the effect ofmutations. As an attempt to fill this gap, replicator-mutator models have been developed,starting with the work of Kimura [17] and followed by models considering different types ofmutations, both in the local [6], [1, 2], [4] and nonlocal [17], [12], [8, 9, 10], [13, 14] cases. ENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION It is important to mention the diversity of research areas where this type of model is applied:population genetics [15], game theory [7], auto-catalytic reaction networks [21] and languageevolution [18]. As pointed out by Schuster and Sigmund [20], in the ordinary differentialequation case, several evolutionary models in different biological fields lead independentlyto the same class of replicator dynamics, showing some sort of universal structure; this ideais also developed in [19], where authors show how apparently very different formulations ofevolutionary dynamics are part of a single unified framework given by the replicator-mutatorequation.When mutations are modeled by the local diffusion operator, and under the constrain ofmass, the replicator-mutator equation typically takes the form(1.7) ∂u∂t = σ ∂ u∂x (cid:124) (cid:123)(cid:122) (cid:125) mutations + u (cid:18) W ( x ) − (cid:90) R W ( y ) u ( t, y ) dy (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) selection . In the context of evolutionary biology, u ( t, x ) represents the density of population, at time t ,per unit of phenotypic trait on a one-dimensional phenotypic trait space. The function W ( x )stands for the fitness of the phenotype x and models the individual reproductive success.Hence the nonlocal term u ( t ) = W ( t ) = (cid:82) R W ( y ) u ( t, y ) dy is the mean fitness at time t , andcan be seen as a Lagrange multiplier for the mass to be conserved, thus yielding an equationon the frequencies.Due to their nonlocal structure, the analysis of replicator-mutator equations often requiresnew approaches compared with local- pde (e.g. comparison principle does not hold), both fromthe analytic and numerical viewpoint. In the framework of (1.7), we mention the works [2],[4] for the cases of quadratic or confining fitness functions, but now stick to the case wherethe fitness linearly depends on the phenotypic trait, say W ( x ) = x .In this setting, the authors of [1] proved that the solution of the replicator mutator Cauchyproblem (1.7) with linear fitness can be written explicitly in terms of the solution of the Heatequation. On the other hand, when a probability density (or kernel) J models mutations, theequation is recast ∂u∂t = J ∗ u − u + u (cid:18) x − (cid:90) R y u ( t, y ) dy (cid:19) , for which a Cumulant Generating Functions (CGF) approach has been developed in [13]: itturns out that the CGF satisfies a first order non-local partial differential equation that canbe explicitly solved, thus giving access to many information such as mean fitness (or meantrait since W ( x ) = x ), variance, position of the leading edge. In the present paper, we shallcombine these two techniques to first reach an implicit expression of the mean fitness u ( t ) ofthe solution to (1.1), and then to obtain a so called implicit/explicit expression of the solution u ( t, x ) itself. Additionally, this procedure allow us to develop a simple numerical strategy forsolving the Cauchy problem associated to (1.1).In this work, and in contrast with (1.7), we consider the case when the fitness function ismodulated by u ( t ) = x ( t ), the mean-fitness at time t , as can be seen in (1.1). Let us nowmake a short description of the emergence of model (1.1) in [23] to which we refer for furtherdetails.The original mutator model, yielding the aforementioned replicator-mutator model (1.7),is the continuous time model(1.8) ∂u∂t ( t, x ) = ( x − x ( t )) u ( t, x ) , DENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION and considers the so-called
Malthusian fitness , that is the rate of growth of a particulargenotype.On the other hand, in the case of non-overlapping generations, one may start from thediscrete time model(1.9) u ( t + 1 , x ) = xu ( t, x ) x ( t )for the change in the so-called Darwinian fitness (success in propagating genes to the nextgeneration) distribution. The Darwinian fitness being nonnegative, equation (1.9) is supple-mented with u ( t = 0 , · ) supported in [0 , + ∞ ). Model (1.9) is immediately recast(1.10) u ( t + 1 , x ) − u ( t, x ) = x − x ( t ) x ( t ) u ( t, x ) . Formally, at least for small times and narrow distributions, the above is approached by thecontinuous in time selection model(1.11) ∂u∂t ( t, x ) = x − x ( t ) x ( t ) u ( t, x ) , which becomes (1.1) after inserting the effect of mutations through a diffusion term. Noticehowever that, in this derivation of (1.1), the fitness of an individual with trait x depends onlyon x and thus selection is not density dependent.In [23], the authors consider directed evolution, that is a laboratory technique that mimicsnatural evolution and can be used for example to acquire proteins with new or improved prop-erties. More precisely, they consider a high-throughput compartmentalized directed evolutionprocess. Genotypes inside the compartments have different phenotypes. They not only pooltheir activity but also share the total number of produced copies, which makes the selectiondensity dependent. In this process, the analogous of (1.9) is given by(1.12) u ( t + 1 , x ) = [ g ( l ) x + (1 − g ( l )) x ( t )] u ( t, x ) x ( t )where the constant g (+ ∞ ) = 0 < g ( l ) < g (0) depends on l , a Poisson parametermeasuring the occupancy of compartments. The analogous of (1.10) is then(1.13) u ( t + 1 , x ) − u ( t, x ) = g ( l ) x − x ( t ) x ( t ) u ( t, x ) . Compared to (1.10), the presence of the coefficient g ( l ) < g ( l ) is small (meaning l large), the process is sloweddown and the validity of (1.11) as a continuous in time approximation should hold in muchmore cases than previously, meaning for larger time periods but also for more various shapesof distribution.Hence, the compartmentalization that takes place in directed evolution (but also in naturalsystems like viruses with multiple infections) and the associated frequency dependence arethe key tools to reach the continuous time model (1.1), starting from a discrete time modelwritten in terms of the Darwinian fitness, see [23].Last, our second main focus, namely equation (1.3), corresponds to the replicator-mutatormodel (1.7), but with the additional effect that the mutations are frequency-dependent: thediffusion coefficient is modulated by the mean trait u ( t ) = x ( t ) and is thus “socially deter-mined”. ENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION Main results
We here state our main results on (1.1), those on (1.3) will be gathered in the final section,where we transfer our developments on the solution u of (1.1) to the solution v of (1.3).We first need to define an admissible set where to look after solutions. We denote by A ,the set of non-negative functions f ∈ L ( R ) such that (cid:90) R f ( x ) dx = 1 , which decrease faster than any exponential function, that is(2.1) lim x →±∞ f ( x ) e α | x | = 0 , ∀ α > , and, last, such that m := (cid:90) R xf ( x ) dx > . Remark . Notice that the case m < m > m = 0 is singular, as clear from the equation and the Gaussian case studiedin subsection 3.2. Notice also that assumption (2.1) could be relaxed by only assuming limit x → + ∞ when m > x → −∞ when m < iii ) in thebelow definition by adding some integrability properties as x → −∞ . For simplicity, in thiswork, we stick to (2.1). Definition 2.2 (Notion of solution) . Let u ∈ A be given. We say that u = u ( t, x ) is a globalsolution of (1.1) starting from u if ( i ) u ∈ C ∞ ((0 , + ∞ ) × R ) , ( ii ) for all t ≥ , u ( t, · ) ∈ A , ( iii ) for all t > , ∂ t u ( t, · ) , ∂ x u ( t, · ) and ∂ xx u ( t, · ) decrease faster than any exponential func-tion, in the sense of (2.1) , ( iv ) u solves (1.1) in the classical sense, ( v ) u ( t, · ) → u in L ( R ) , as t → . Our main result consists in the well-posedness of the Cauchy problem (1.1) with, moreover,an implicit/explicit expression of the solution.
Theorem 2.3 (The solution of the Cauchy problem) . Let u ∈ A be given. Then there is aunique global solution of (1.1) starting from u (in the sense of Definition 2.2). Moreover,its mean u ( t ) is implicitly (and uniquely) determined by u ( t ) = (cid:115) m + 2 σ t + 2 (cid:90) t C (cid:48)(cid:48) (cid:18)(cid:90) s dyu ( y ) (cid:19) ds where m = (cid:82) R xu ( x ) dx > and C ( z ) := ln (cid:0)(cid:82) R u ( x ) e zx dx (cid:1) , z ≥ , is the cumulantgenerating function ( cgf ) of u . Last, u ( t, x ) is given by u ( t, x ) = w (cid:18) t, x + 2 (cid:90) t (cid:90) s σ u ( τ ) dτ ds (cid:19) exp (cid:32) − t + x (cid:90) t dsu ( s ) + (cid:90) t σ (cid:18)(cid:90) s dτu ( τ ) (cid:19) ds (cid:33) , where w = w ( t, y ) is the solution of the Heat equation ∂ t w = σ ∂ yy w starting from u . DENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION
The proof is based in a combination of the approaches of [13] and [1]. In the course of theproof, we collect the expressions and some estimates on the mean ¯ u ( t ) and the variance V ( t ) := (cid:90) R ( x − ¯ u ( t )) u ( t, x ) dx = (cid:90) R x u ( t, x ) dx − (cid:18)(cid:90) R xu ( t, x ) dx (cid:19) of the solution at time t . Corollary 2.4 (Long time behaviour) . Let u ∈ A be given. Then the mean u ( t ) of the globalsolution given by Theorem 2.3 satisfies u ( t ) (cid:40) ≥ √ σ t in any case ∼ √ σ t if sup supp( u ) < + ∞ , and, in any case, (2.2) (cid:90) + ∞ dyu ( y ) = + ∞ . The variance V ( t ) of the global solution satisfies V ( t ) (cid:40) ≥ σ t in any case ∼ σ t if sup supp( u ) < + ∞ , where supp( u ) denotes the support of u . The above shows that, as t → + ∞ , the solution moves to the right since u ( t ) → + ∞ , and isflattening since V ( t ) → + ∞ . In particular, for one side compactly supported initial data, thepropagation (asymptotically) occurs at constant speed √ σ , which is in contrast with theacceleration phenomenon proved in [1]. This is also true for Gaussian data as will be noticedin subsection 3.2. The role of (2.2) is to provide “a kind of upper bound” on u ( t ) when, forexample, the initial tails are still lighter than any exponential but heavier than Gaussian, e.g.of the magnitude e − x α (1 < α <
2) or even e − x ln x as x → + ∞ .The paper is organized as follows. In Section 3, we investigate special solutions, in particularGaussian ones which provide a preliminary understanding of equation (1.1). In Section 4, webegin the proof of our main result, Theorem 2.3, by using the cgf approach introducedin [13]. The proof of Theorem 2.3 and its corollary are completed in Section 5 thanks tothe methodology of [1]. As a by product of our analysis, we collect a numerical strategyfor solving the Cauchy problem. Last, Section 6 is devoted to obtain the companion resultson (1.3) by expressing v ( t, x ) in terms of u ( t, x ).3. Special solutions
In this section, we investigate special solutions to (1.1), in particular non-negative, inte-grable steady states and Gaussian solutions. In the first case we prove the non existence ofsuch steady states. This is due to the particular form of the fitness function and is an analo-gous result to that one obtained by Alfaro an Carles in [1]. The situation is different in thecase of a confining fitness function [4, 2] where the existence and uniqueness of a non-negativestationary solution is ensured and corresponds to the ground-state of the Schr¨odinger Hamil-tonian where the potential is the opposite of fitness function. In the second case, Gaussian
ENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION solutions are computed explicitly by solving the differential equations describing the evolutionof the mean and the inverse of the variance.3.1. Steady states.Proposition 3.1 (Steady state) . There is no nontrivial non-negative steady state φ = φ ( x ) solving (1.1) and satisfying φ ( ±∞ ) = 0 .Proof. A steady state with x (cid:54) = 0 must solve φ (cid:48)(cid:48) ( x ) + σ x − xx φ ( x ) = 0 for all x ∈ R . Letting ψ ( x ) := φ ( x − ( σ x ) / x ), this is recast ψ (cid:48)(cid:48) ( x ) − xψ ( x ) = 0 , x ∈ R . Hence ψ ( x ) is a linear combination of the Airy functions Ai( x ) and Bi( x ). From ψ (+ ∞ ) = 0,we deduce that ψ ( x ) is a multiple of Ai( x ). Hence either it is trivial, or it changes sign. (cid:3) Gaussian solutions.
We investigate the propagation of a Gaussian initial data, whichis relevant for biological applications. In the context of evolutionary genetics, families ofGaussian solutions for nonlinear and nonlocal equations can be found in [6], [1, 2]. In a differ-ent context involving logarithmic non-linearities, we also refer to [5], [11] for the Schr¨odingerequation and to [3] for the Heat equation.
Proposition 3.2 (Propagation of Gaussian initial data) . For a > and m ∈ R , let usdefine (3.1) a ( t ) := a a σ t , m ( t ) := + (cid:113) σ t + ta + m when m > ± (cid:113) σ t + ta when m = 0 − (cid:113) σ t + ta + m when m < . Then (3.2) u ( t, x ) := (cid:114) a ( t )2 π e − a ( t ) ( x − m ( t ))22 solves (1.1) in (0 , + ∞ ) × R and starts from u ( x ) = (cid:112) a π e − a x − m .Proof. From the ansatz (3.2), we compute ∂ t u ( t, x ) = u ( t, x ) (cid:20) a (cid:48) ( t )2 a ( t ) − a (cid:48) ( t ) ( x − m ( t )) a ( t ) m (cid:48) ( t )( x − m ( t )) (cid:21) ∂ xx u ( t, x ) = u ( t, x ) (cid:2) − a ( t ) + a ( t )( x − m ( t )) (cid:3) u ( t, x ) x − x ( t ) x ( t ) = u ( t, x ) x − m ( t ) m ( t ) . We plug the above into equation (1.1), and identify the x , the x and the x coefficients toobtain three differential equations. The first one is(3.3) − a (cid:48) ( t )2 = σ a ( t ) , whose solution, starting from a , is given by the first part in (3.1). Using (3.3), we see thatthe two other equations reduce to a ( t ) m ( t ) m (cid:48) ( t ) = 1, which is solved as m ( t ) = 2 σ t + 2 ta + m , DENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION and thus the second part in (3.1). (cid:3)
Notice that functions m ( t ) and V ( t ) = a ( t ) are respectively the mean and the varianceof the density u ( t, · ). Hence Proposition 3.2 shows that, starting from a Gaussian profile,the solution remains a Gaussian function, is asymptotically propagating at constant speedand flattening since m ( t ) ∼ ±√ σ t , V ( t ) ∼ σ t , as t → + ∞ , see Figure 1. Notice thatthe direction of propagation is determined by the initial value m : towards the right when m >
0, the left when m <
0. The case m = 0 is singular because of the multiplicity ofsolutions to (3.1), two Gaussian solutions emerge, propagating to the left and to the right,see Figure 1. Propagation - -
20 0 20 400.5 m = - PropagationPropagation - -
20 0 20 400.5 m = Propagation - -
20 0 20 400.5 m = Figure 1.
Evolution of Gaussian solutions for σ = 1, a = 1 and (from leftto right) m = − m = 0 and m = 4.4. The CGF approach
In this section, we assume that we are equipped with a solution of (1.1) starting from u ∈ A . We define the Cumulative Generating Function ( cgf )(4.1) C ( t, z ) := ln (cid:18)(cid:90) R u ( t, x ) e zx dx (cid:19) , t ≥ , z ≥ , of such a solution, in the spirit of [13]. From Definition 2.2, C ( t, z ) is well defined and smoothon (0 , + ∞ ) × [0 , + ∞ ). We shall derive an implicit expression for C ( t, z ), and then for thevalue u ( t ) of the mean. This crucial step will enable us to complete the analysis in the nextsection, which is much in the spirit of [1].The following computations are validated by the notion of solution adopted in Definition 2.2and, possibly, the dominated convergence theorem. Observe that ∂ t C ( t, z ) = (cid:82) R ∂ t u ( t, x ) e zx dx (cid:82) R u ( t, x ) e zx dx , and that ∂ z C ( t, z ) = (cid:82) R xu ( t, x ) e zx dx (cid:82) R u ( t, x ) e zx dx , ∂ zz C ( t, z ) = (cid:82) R x u ( t, x ) e zx dx (cid:82) R u ( t, x ) e zx dx − (cid:18) (cid:82) R xu ( t, x ) e zx dx (cid:82) R u ( t, x ) e zx dx (cid:19) . Notice in particular that the mean is reached through ∂ z C ( t, z = 0) = (cid:90) R xu ( t, x ) dx = u ( t ) , ENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION whereas the variance is reached through ∂ zz C ( t, z = 0) = (cid:90) R x u ( t, x ) dx − (cid:18)(cid:90) R xu ( t, x ) dx (cid:19) = V ( t ) . Hence, multiplying equation (1.1) by e zx , integrating over x ∈ R and finally dividing by (cid:82) R u ( t, x ) e zx , we obtain the following nonlocal first order pde (4.2) ∂ t C ( t, z ) = σ z − ∂ z C ( t, z ) ∂ z C ( t, , where we have used integration by parts to get (cid:90) R ∂ xx u ( t, x ) e zx dx = z (cid:90) R u ( t, x ) e zx dx. Furthermore, since (cid:82) R u ( t, x ) dx = 1, the condition C ( t,
0) = 0 must hold for any t ≥
0. As aresult we are facing the problem(4.3) ∂ t C ( t, z ) = σ z − ∂ z C ( t, z ) u ( t ) t ≥ , z ≥ ,C (0 , z ) = C ( z ) z ≥ ,C ( t,
0) = 0 t ≥ , where C ( z ) = ln (cid:0)(cid:82) R u ( x ) e zx dx (cid:1) is the cgf of the initial data u . Notice that, from astraightforward computation, C (cid:48)(cid:48) ( z ) = (cid:0)(cid:82) R x u ( x ) e zx dx (cid:1) (cid:0)(cid:82) R u ( x ) e zx dx (cid:1) − (cid:0)(cid:82) R xu ( x ) e zx dx (cid:1) (cid:0)(cid:82) R u ( x ) e zx dx (cid:1) . Hence, from Cauchy Schwarz inequality, we see that C (cid:48)(cid:48) ( z ) ≥ z ≥
0, and even C (cid:48)(cid:48) ( z ) > z ≥ u (cid:54)≡
0. This convexity property of cumulant generating functionsis well-known in probability, and will be used in the following.Fix t > z >
0. For s such that − t ≤ s and (cid:82) s dτu ( t + τ ) ≤ z , set ψ ( s ) := C (cid:18) t + s, z − (cid:90) s dτu ( t + τ ) (cid:19) , which we differentiate to get ψ (cid:48) ( s ) = ∂ t C (cid:18) t + s, z − (cid:90) s dτu ( t + τ ) (cid:19) − u ( t + s ) ∂ z C (cid:18) t + s, z − (cid:90) s dτu ( t + τ ) (cid:19) = σ (cid:18) z − (cid:90) s dτu ( t + τ ) (cid:19) − . (4.4)As a result, C ( t, z ) − C (cid:18) z + (cid:90) t dτu ( τ ) (cid:19) = C ( t, z ) − C (cid:18) z + (cid:90) − t dτu ( t + τ ) (cid:19) = ψ (0) − ψ ( − t )= (cid:90) − t ψ (cid:48) ( s ) ds. DENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION
From (4.4), we deduce that C ( t, z ) − C (cid:18) z + (cid:90) t dτu ( τ ) (cid:19) = σ (cid:90) − t (cid:18) z − (cid:90) s dτu ( t + τ ) (cid:19) ds − t = σ tz − σ z (cid:90) − t (cid:90) s dτu ( t + τ ) ds + σ (cid:90) − t (cid:18)(cid:90) s dτu ( t + τ ) (cid:19) ds − t = σ tz + 2 σ z (cid:90) − t (cid:90) τ − t u ( t + τ ) dsdτ + σ (cid:90) − t (cid:18)(cid:90) s dτu ( t + τ ) (cid:19) ds − t, by Fubini-Tonelli theorem. We thus conclude that C ( t, z ) = C (cid:18) z + (cid:90) t dτu ( τ ) (cid:19) + σ tz + 2 σ z (cid:90) t yu ( y ) dy + σ (cid:90) − t (cid:18)(cid:90) s dτu ( t + τ ) (cid:19) ds − t. This is an implicit expression for C ( t, z ) since it still involves u ( t ) = ∂ z C ( t, z and evaluating at z = 0, we reach another implicit formula for the meanvalue of the solution(4.5) u ( t ) = C (cid:48) (cid:18)(cid:90) t dyu ( y ) (cid:19) + 2 σ (cid:90) t yu ( y ) dy, whereas differentiating twice with respect to z and evaluating at z = 0, we reach the varianceof the solution(4.6) V ( t ) = C (cid:48)(cid:48) (cid:18)(cid:90) t dyu ( y ) (cid:19) + 2 σ t (cid:40) ≥ σ t in any case ∼ σ t if sup supp( u ) < + ∞ , where the last estimate follows from [13, Lemma 4.5].Next, differentiating (4.5) with respect to t , we get u (cid:48) ( t ) = 1 u ( t ) C (cid:48)(cid:48) (cid:18)(cid:90) t dyu ( y ) (cid:19) + 2 σ tu ( t ) , so that ddt u ( t ) = 2 V ( t ) and thus(4.7) u ( t ) (cid:40) ≥ √ σ t in any case ∼ √ σ t if sup supp( u ) < + ∞ . Also, integrating ddt u ( t ) = 2 V ( t ), we collect the very useful expression(4.8) u ( t ) = (cid:115) m + 2 σ t + 2 (cid:90) t C (cid:48)(cid:48) (cid:18)(cid:90) s dyu ( y ) (cid:19) ds . Remark . Notice that if u ( x ) = (cid:112) a π e − a x − m is a Gaussian initial data as in subsec-tion 3.2, we know that its cgf is C ( z ) = m z + 12 a z . ENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION Hence u ( t ) = (cid:113) m + 2 σ t + a t and V ( t ) = a + 2 σ t , which agrees with the values of m ( t )and a ( t ) in Proposition 3.2.5. The solution implicitly/explicitly
In this section, we mainly complete the proof of Theorem 2.3 and of Corollary 2.4. As aby product, we present some numerical strategies for solving (1.1).First, assume that we are equipped with a solution u ( t, x ) of (1.1) starting from u ∈ A .In particular we know from Section 4 that u ( t ) is given by (4.8). Following [1], we write(5.1) u ( t, x ) = w (cid:18) t, x + 2 (cid:90) t (cid:90) s σ u ( τ ) dτ ds (cid:19) e − t + x (cid:82) t dsu ( s ) + (cid:82) t σ (cid:16)(cid:82) s dτu ( τ ) (cid:17) ds . This clearly defines a unique w ( t, y ), t ≥ y ∈ R , which has to satisfy(5.2) (cid:40) ∂ t w = σ ∂ yy ww | t =0 = u . We refer to [1] for more details on such a change of unknown function, which enables to reducesome replicator-mutator equations to the Heat equation.Now, our main task is to show that the problem (4.8), is globally well-posed.
Lemma 5.1.
For any given
T > , there is a unique solution u ∈ C ([0 , T ]) to (5.3) u ( t ) = (cid:115) m + 2 σ t + 2 (cid:90) t C (cid:48)(cid:48) (cid:18)(cid:90) s dyu ( y ) (cid:19) ds . Proof.
Recall that m >
0. We define the following subspace of C ([0 , T ]), equipped with the L ∞ norm, X := { h ∈ C ([0 , T ]) , h ( t ) ≥ h (0) = m for all 0 ≤ t ≤ T } , which is a Banach space. Now, we define the operator L : h ∈ X (cid:55)→ q ∈ X , and q is given by q ( t ) := (cid:115) m + 2 σ t + 2 (cid:90) t C (cid:48)(cid:48) (cid:18)(cid:90) s dτh ( τ ) (cid:19) ds . The proof consists in mimicking that of the Banach fixed point theorem. In order to proveexistence, we introduce the sequence ( q n ) ∈ X N defined inductively by(5.4) q ( t ) = m q n +1 ( t ) = (cid:114) m + 2 σ t + 2 (cid:82) t C (cid:48)(cid:48) (cid:16)(cid:82) s dτq n ( τ ) (cid:17) ds . DENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION
We have | q n +1 ( t ) − q n ( t ) | = (cid:12)(cid:12) q n +1 ( t ) − q n ( t ) (cid:12)(cid:12) q n +1 ( t ) + q n ( t ) ≤ m (cid:90) t (cid:12)(cid:12)(cid:12)(cid:12) C (cid:48)(cid:48) (cid:18)(cid:90) s dτq n ( τ ) (cid:19) − C (cid:48)(cid:48) (cid:18)(cid:90) s dτq n − ( τ ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ds ≤ m (cid:107) C (cid:107) L ∞ (0 ,T ) (cid:90) t (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) s (cid:18) q n ( τ ) − q n − ( τ ) (cid:19) dτ (cid:12)(cid:12)(cid:12)(cid:12) ds ≤ m (cid:107) C (cid:107) L ∞ (0 ,T ) (cid:90) t (cid:90) s | q n ( τ ) − q n − ( τ ) | m dτ ds ≤ m (cid:107) C (cid:107) L ∞ (0 ,T ) (cid:90) t (cid:90) s (cid:107) q n − q n − (cid:107) L ∞ (0 ,s ) dτ ds ≤ m (cid:107) C (cid:107) L ∞ (0 ,T ) (cid:90) t s (cid:107) q n − q n − (cid:107) L ∞ (0 ,s ) ds ≤ k T (cid:90) t (cid:107) q n − q n − (cid:107) L ∞ (0 ,s ) ds, where k := m (cid:107) C (cid:107) L ∞ (0 ,T ) . Then we straightforwardly prove by induction that, for any n ≥
0, any 0 ≤ s ≤ T ,(5.5) (cid:107) q n +1 − q n (cid:107) L ∞ (0 ,s ) ≤ (cid:107) q − q (cid:107) L ∞ (0 ,T ) ( k T s ) n n ! . In particular, the series (cid:80) (cid:107) q n +1 − q n (cid:107) L ∞ (0 ,T ) converges and therefore ( q n ) converges uni-formly in C (0 , T ) to some q which is a fixed point of L .As far as uniqueness is concerned, if q and ˜ q are two fixed points of L then the argumentto reach (5.5) also yields (cid:107) q − ˜ q (cid:107) L ∞ (0 ,T ) ≤ (cid:107) q − ˜ q (cid:107) L ∞ (0 ,T ) ( kT ) n n ! . Letting n → + ∞ enforces q ≡ ˜ q . (cid:3) Completion of the proof of Theorem 2.3 and Corollary 2.4.
Hence, any solution in the senseof Definition 2.2 has to be given by (5.1), where u ( t ) is given by Lemma 5.1. Notice that,from (5.3), u has to be smooth.Conversely, it is now a matter of straightforward computations — based on (5.1),(5.2)and (5.3)— to check that this does provide the solution. Notice that, in particular, the initialdata is understood in the sense of Definition 2.2 ( v ) since, in the Cauchy problem for the heatequation (5.2), the initial data is understood in the sense w ( t, · ) → u in L ( R ) as t → (cid:3) Numerical implications.
Three major difficulties are encountered when setting up a nu-merical strategy for problem (1.1). The first one is the nonlocal nature of the equation. Atthis stage, two natural approaches can be considered: i ) use finite differences and approximatethe nonlocal term by a Riemann sum; ii ) apply a splitting method, computing alternativelythe solution of the non-local term and the resulting local reaction-diffusion equation. ENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION The second difficulty is the unboundedness of the domain. To manage it, one usuallyimposes artificial boundary conditions and solves numerically the equation on a sufficientlylarge domain, approximating the true solution by the latter.The third complexity comes from the propagation of the solution, at least linear in viewof (4.7), making it difficult to track it over time.Previously, we gave an implicit/explicit construction of the solution u = u ( t, x ) of theCauchy problem associated with equation (1.1), which actually provides a strategy for solvingnumerically the problem.The idea is to initially find an approximation of the nonlocal term u ( t ) on a time interval[0 , T ] , T >
0. This step is performed using the fixed-point iteration sequence (5.4), forwhich the estimate (5.5) provides a convergence rate. This requires in particular to compute(analytically or numerically) the cumulant generating function C ( z ) of the initial datum u ( x ).The next step consists in calculating the solution w = w ( t, x ) of the heat equation withinitial datum u ( x ) and evaluating it in the values indicated by formula (5.1). Obviously, inthe case where an explicit solution w = w ( t, x ) of the heat equation can be found, artificialboundary conditions for solving the heat equation are not needed, leading to a better numeri-cal solution. On the other hand, in the case where no closed form for w = w ( t, x ) is available,it is important to previously solve numerically the heat equation in a sufficiently large spatialdomain. For instance, if we want to solve numerically (1.1) in [ x − , x + ] × [0 , T ] ⊂ R × R ,assuming m > m < w needs to be evaluated in space up to x = x + + 2 (cid:82) t (cid:82) s σ u ( τ ) dτ ds and since u ( t ) ≥ m >
0, then we should solve the heat equation at most in [ x − , x + + σ T m ] × [0 , T ].In figure 2b, we have represented the numerical solution obtained by this method for theinitial condition u ( x ) = [ 12 ,
32 ] ( x ), so that m = 1 and C ( z ) = ln (cid:16) e z sinh( z/ z (cid:17) , and σ = 1;in this case, w is known to be given by w ( t, x ) = 12 (cid:18) erf (cid:18) − x √ t (cid:19) − erf (cid:18) − x √ t (cid:19)(cid:19) . When the mutation coefficient is socially determined
This section is devoted to the analysis of equation (1.3) supplemented with a non-negativeinitial data v ∈ L ( R ) satisfying (cid:82) R v ( x ) dx = 1, and decreasing faster than any exponentialdata in the sense of (2.1). We denote m := (cid:82) R xv ( x ) dx . Interestingly, one can avoid theassumption m > Gaussian solutions: concentration vs. extinction.Proposition 6.1 (Propagation of Gaussian initial data) . For a > and m ∈ R . Thenthere is a unique Gaussian (6.1) v ( t, x ) := (cid:114) a ( t )2 π e − a ( t ) ( x − m ( t ))22 which solves, at least locally in time, (1.3) and starts from u ( x ) = (cid:112) a π e − a x − m . Itsbehaviour depends strongly on the parameters a > , m ∈ R and σ > . DENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION (a) ×××××××××× × × - - x u (b) Figure 2. (a): The first four approximations, for 0 ≤ t ≤
3, of the nonlocalterm u ( t ) computed via the fixed point iteration (5.4). (b): Numerical solutionobtained by the method described in Section 5, starting from u = [1 / , / ,with σ = 1. The red points are the points on the graph u ( t, · ) with ab-scissa x = u ( t ). The green points are the maxima of u ( t, · ). This reveals thedissymmetry of the solution.( i ) Assume m < − a √ σ . Then a ( t ) blows up in finite time: there is < T (cid:63) < + ∞ suchthat m ( t ) → m ( T (cid:63) ) < and V ( t ) := 1 a ( t ) → , as t (cid:37) T (cid:63) . ( ii ) Assume m = − a √ σ . Then a ( t ) and m ( t ) are global and m ( t ) (cid:37) and V ( t ) := 1 a ( t ) → , as t → + ∞ . ( iii ) Assume m > − a √ σ . Then a ( t ) and m ( t ) are global and m ( t ) ∼ C m e √ σ t → + ∞ and V ( t ) := 1 a ( t ) ∼ C V e √ σ t → + ∞ , as t → + ∞ , where C m := a m √ σ + √ a √ σ > and C V := a m √ σ + √ √ a > . Remark . In Figure 3, we have subdivided the phase plane ( m , a ) into regions corre-sponding to the cases ( i ) , ( ii ) and ( iii ).Case ( i ) corresponds to a concentration phenomena in finite time: the Gaussian solutionconverges, in finite time, to a Dirac mass centered at x = m ( T (cid:63) ) <
0. See Figure 4.Case ( ii ) corresponds to a concentration phenomena in infinite time: the Gaussian solutionconverges, at large times, to a Dirac mass centered at x = 0. See Figure 5.In case ( iii ), if m < i ) and ( ii ), the mean of the solu-tion manages to “cross” the value zero. In case ( iii ) convergence to a accelerating ( m ( t ) ∼ C m e √ σ t ) and flattening ( V ( t ) ∼ C V e √ σ t ) profile always occurs at large times. Moreover,the variance of the solution is decreasing before the mean reaches zero, and then starts toincrease after the mean crosses zero, which reveals an anti-diffusion/diffusion phenomenon,see Figure 6. On the other hand, if m ≥
0, the solution does not anti-diffuse and onlydiffuses while the mean tends to infinity, see Figure 7.
ENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION Figure 3.
Vector field defined by the differential system (6.2) with σ = 1,describing the dynamics of Gaussian solutions. In yellow, the set of initialconditions for which a blows up in finite time T (cid:63) and in red, dark blue and lightblue those for which both a and m are globally defined. The red dashed curveis the set of values defined by m = − / ( a √ σ ), for which a tends to infinityand m tends to zero as time goes to infinity. The dark blue region correspondsto the values leading to an anti-diffusion/diffusion behaviour. The light blueregion corresponds to the values leading to a pure diffusion behaviour. Proof.
As in the proof of Proposition 3.2, we plug the ansatz (6.1) into (1.3) and arrive at(6.2) (cid:40) m (cid:48) ( t ) = a ( t ) a (cid:48) ( t ) = − σ a ( t ) m ( t ) , so that m (cid:48)(cid:48) ( t ) = 2 σ m ( t ), m (0) = m , m (cid:48) (0) = a which is globally solved as(6.3) m ( t ) = √ a m √ σ a √ σ e √ σ t − √ − a m √ σ a √ σ e −√ σ t . From the first equation in (6.2) we deduce(6.4) a ( t ) = 2 √ a ( √ a m √ σ ) e √ σ t + ( √ − a m √ σ ) e −√ σ t , DENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION
Propagation - - - x Figure 4.
Case ( i ) concentration in finite time. The values of the parametersare a = 5 / m = − / σ = 1. It follows that T (cid:63) ≈ .
878 and m ( T (cid:63) ) ≈ − . < Propagation - - - - - - x Figure 5.
Case ( ii ) concentration in infinite time: the solution converges toa Dirac mass at zero. The values of the parameters are a = 3 / m = − √ / σ = 1.as long as blow-up does not occur. We now easily see that in case ( i ) blow up occurs at time T (cid:63) = 12 √ σ ln (cid:32) a m √ σ − √ a m √ σ + √ (cid:33) ∈ (0 , + ∞ ) . In case ( ii ) blow up does not occur, m ( t ) = m e −√ σ t and a ( t ) = a e √ σ t . In case ( iii )blow up does not occur and conclusion follows from (6.3) and (6.4). (cid:3) General solutions.
Here we consider a non-negative initial data v ∈ L ( R ) satisfying (cid:82) R v ( x ) dx = 1, (2.1) and m := (cid:82) R xv ( x ) dx >
0, and investigate possible solutions v = v ( t, x ) of (1.3) starting from v . From Theorem 2.3, we are equipped with the unique solution u = u ( t, x ) of equation (1.1) starting from v . Roughly speaking, the understanding of (1.3) ENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION Propagation -
10 10 20 x Figure 6.
Case ( iii ) with m <
0, anti-diffusion/diffusion behaviour. Thevalues of the parameters are a = 2 / m = − / σ = 1. Propagation x Figure 7.
Case ( iii ) with m ≥
0, the solution is flattening and accelerating.The values of the parameters are a = 3 / m = 7 / σ = 1.now reduces to that of the o.d.e. Cauchy problem(6.5) (cid:40) ϕ (cid:48) ( t ) = u ( ϕ ( t )) for t > ,ϕ (0) = 0 . Indeed, defining(6.6) v ( t, x ) := u ( ϕ ( t ) , x ) , we have v ( t ) = u ( ϕ ( t )) and one checks that, since u solves (1.1) (and starts from v ), v = v ( t, x ) solves (1.3) (and starts from v ) as long as the solution ϕ ( t ) of (6.5) exists. Lemma 6.3.
The solution ϕ ( t ) of the o.d.e. Cauchy problem (6.5) is global.Proof. If ϕ blows up in finite time, say at time T >
0, then we obtain from the o.d.e. that (cid:90) + ∞ dyu ( y ) = T. Hence we deduce from (4.5) that u ( t ) ≤ (cid:107) C (cid:48) (cid:107) L ∞ (0 ,T ) + 2 σ (cid:90) t yu ( y ) dy. But, from (4.8), we know that yu ( y ) ≤ √ σ and thus u ( t ) ≤ (cid:107) C (cid:48) (cid:107) L ∞ (0 ,T ) + √ σ t, DENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION which contradicts u ∈ L (0 , + ∞ ). This concludes the proof of the lemma and, clearly, ofestimate (2.2). (cid:3) Since u is smooth, so is ϕ . Also, ϕ (cid:48) ( t ) ≥ m > ϕ ( t ) tends to + ∞ as t → + ∞ . Hence ϕ : (0 , + ∞ ) → (0 , + ∞ ) is a smooth diffeomorphism. In other words, (6.6)shows that there is a one-to-one correspondence between solutions to (1.3) and that of (1.1).We omit the full details but this clearly concludes the analysis of the Cauchy problem as-sociated with (1.3), whose surprising qualitative behaviours have already been underlined insubsection 6.1. Acknowledgements.
The authors are very grateful to Anton Zadorin for sharing his exper-tise on the models appearing in [23], and for a precise relecture of the preprint, in particularthe introduction. M. Alfaro is supported by the ANR i-site muse , project michel ◦ ANR idex -0006). M. Veruete is supported by the National Council for Science andTechnology of Mexico.
References [1]
M. Alfaro and R. Carles , Explicit solutions for replicator-mutator equations: extinction versus accel-eration , SIAM J. Appl. Math., 74 (2014), pp. 1919–1934.[2]
M. Alfaro and R. Carles , Replicator-mutator equations with quadratic fitness , Proc. Amer. Math. Soc.,145 (2017), pp. 5315–5327.[3]
M. Alfaro and R. Carles , Superexponential growth or decay in the heat equation with a logarithmicnonlinearity , Dyn. Partial Differ. Equ., 14 (2017), pp. 343–358.[4]
M. Alfaro and M. Veruete , Evolutionary branching via replicator–mutator equations , J. Dynam. Dif-ferential Equations, (2018).[5]
I. Bia(cid:32)lynicki-Birula and J. Mycielski , Nonlinear wave mechanics , Ann. Physics, 100 (1976), pp. 62–93.[6]
V. N. Biktashev , A simple mathematical model of gradual darwinian evolution: emergence of a gaussiantrait distribution in adaptation along a fitness gradient , J. Math. Biol., 68 (2014), pp. 1225–1248.[7]
I. Bomze and R. Burger , Stability by mutation in evolutionary games , Games Econom. Behav., 11(1995), pp. 146–172.[8]
R. B¨urger , On the maintenance of genetic variation: global analysis of Kimura’s continuum-of-allelesmodel , J. Math. Biol., 24 (1986), pp. 341–351.[9] ,
Mutation-selection balance and continuum-of-alleles models , Math. Biosci., 91 (1988), pp. 67–83.[10] ,
Perturbations of positive semigroups and applications to population genetics , Math. Z., 197 (1988),pp. 259–272.[11]
R. Carles and I. Gallagher , Universal dynamics for the defocusing logarithmic schrdinger equation ,Duke Math. J., 167 (2018), pp. 1761–1801.[12]
W. H. Fleming , Equilibrium distributions of continuous polygenic traits , SIAM J. Appl. Math., 36 (1979),pp. 148–168.[13]
M.-E. Gil, F. Hamel, G. Martin, and L. Roques , Mathematical properties of a class of integro-differential models from population genetics , SIAM J. Appl. Math., 77 (2017), pp. 1536–1561.[14]
M.-E. Gil, F. Hamel, G. Martin, and L. Roques , Dynamics of fitness distributions in the presenceof a phenotypic optimum: an integro-differential approach , HAL preprint, (2018).[15]
K. Hadeler , Stable polymorphisms in a selection model with mutation , SIAM J. Appl. Math., 41 (1981),pp. 1–7.[16]
J. Hofbauer and K. Sigmund , The Theory of Evolution and Dynamical Systems: Mathematical Aspectsof Selection , London Mathematical Society Student Texts, Cambridge University Press, 1988.[17]
M. Kimura , A stochastic model concerning the maintenance of genetic variability in quantitative charac-ters. , Proc. Natl. Acad. Sci. USA, 54 (1965), pp. 731–736.[18]
M. Nowak, N. Komarova, and P. Niyogi , Evolution of universal grammar , Science, 291 (2001),pp. 114–118.[19]
K. Page and M. Nowak , Unifying evolutionary dynamics , J. Theoret. Biol., 219 (2002), pp. 93–98.
ENSITY DEPENDENT REPLICATOR-MUTATOR MODELS IN DIRECTED EVOLUTION [20] P. Schuster and K. Sigmund , Replicator dynamics , J. Theoret. Biol., 100 (1983), pp. 533–538.[21]
P. Stadler and P. Schuster , Mutation in autocatalytic reaction networks , J. Math. Biol., 30 (1992),pp. 597–632.[22]
P. Taylor and L. Jonker , Evolutionary stable strategies and game dynamics , Math. Biosci., 40 (1978),pp. 145 – 156.[23]