{'fraction_non_alphanumeric': 0.041652518392756084, 'fraction_numerical': 0.012365591397849462, 'mean_word_length': 4.641021548284118, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A simple mathematical model of gradual Darwinian evolution: Emergence of a Gaussian trait distribution in adaptation along a fitness gradient

We consider a simple mathematical model of gradual Darwinian evolution in continuous time and continuous trait space, due to intraspecific competition for common resource in an asexually reproducing population in constant environment, while far from evolutionary stable equilibrium. The model admits exact analytical solution. In particular, Gaussian distribution of the trait emerges from generic initial conditions.


Introduction
The question considered in this paper is: suppose a population evolves according to the Darwin's mechanism involving mutations and natural selection, and some of its quantitative traits change gradually, what is the rate of this gradual change? This question may be not the most interesting when applied to analysis of past evolution, say fossil records, where the epochs of such gradual changes are relatively short compared to much longer epoch when the species appear unchanged ("punctuated equilibrium", Eldredge and Gould, 1972). However, the speed of evolution is crucial in constant competition of taxa ("Red Queen" hypothesis, Liow et al 2011; "evolutional arms race", Dawkins and Krebs 1979), and is of considerable practical importance in relation to present day phenomena such as adaptation of pathogens to existing methods of treatment, or adaptation of endangered species to changing environmental conditions.
The literature dedicated to this subject is vast. Here we mention here only some cornerstone publications, most relevant to the present communication, to designate its context and motivation.
Quantitative approach to evolution dates back at least to Fisher's (1930) book, which contained his famous "Fundamental Theorem of Natural Selection", stating that the rate of increase of the mean fitness of a population at any moment of time, attributable to to natural selection, equals the genetic variance of fitness of that population at that moment of time. This result is as elegant and powerful as it is difficult to apply correctly, since its deceptively simple words encrypt complicated concepts, as it took nearly 40 years to figure out (Price, 1972). The next question is, of course, what determines this variance in the population fitness, and how to predict it.
Adaptive dynamics is mainly concerned with qualitative questions such as direction of evolution, stability of evolutionary steady states and speciation due to branching (Geritz et al, 1998;Bowers, 2011). On the quantitative level, the fundamental for adaptive dynamics is the "canonical equation". The influential paper by Dieckmann and Law (1996) gives this equation in the form (for a single selected trait) dµ dt = K(µ) where µ(t) is the average value of the trait at time t, r(x, µ) measures fitness of individuals with trait value x in the environment of resitent trait values µ and the coefficient K(µ) is described as a "non-negative coefficient, . . . that scales the rate of evolutionary change". Dieckmann and Law have endeavoured to derive this coefficient based on a certain miscroscopic model of mutation and selection processes, and come to the result that it is equal to the variance of the population with respect to the quantiative trait, which, when applied to the fitness, exactly reproduces the Fisher's result. In their derivation, Dieckmann and Law made an essential assumption, with reference to separation of times between mutations and selection in the limit of slow evolution and the competitive exclusion principle, that at each moment of time, the selection reduces the population to a certain type, which however changes in time due to random mutations ("quasi-monomorphic framework"). This framework, under the name of the "Trait Substitution Sequence" model, has been rigorously justified by Champagnat (2006), using a stochastic model, under certain asymptotic assumption about the mutation rate. Slightly simplifying, the key assumption is that mutations are so rare that for a given population size, there is sufficient time between consecutive mutations for the whole population to convert to the new trait value if it is fitter than the previous. They also consider the opposite limit which they call "large-population limit with accelerated births and deaths", in which the population is so big and mutations so frequent that at any time there the population consists of many different types. This leads to a deterministic model in the form of an integro-differential equation, which is akin to Fisher's reaction-diffusion equation, only for population distribution in the trait space, and is a time-dependent variation of the Kimura (1965) model. Such deterministic models have been studied in many works, e.g. Calsina and Perelló (1995); Gudelj et al (2006); Schuster (2011). These works typically concentrate on the analysis of stationary solutions corresponding to the evolutionary stable states, rather than quantifying the speed with which these states may be approached. Quantitative genetics has developed a number of its own approaches and investigated a great number of complicated problems related to quantitative description of evolution. One approach is through the method of moments. An example is the paper by Barton and Turelli (1987) which considered multilocus determination of a quantitative trait in a sexually reproducing population, and in particular presented an infinite chain of ordinary differential equations for the moments of allelic distribution. A more recent example is a paper by Sasaki and Dieckmann (2011), which looks at multiple-peak character distributions, called "olygomorphs", whereas each of the morphs is treated by the method of moments. The chain of equations for the moments is typically treated using a "closure" procedure, say by assuming that the distribution is Gaussian. Gaussian distributions naturally occurred in a number of studies, e.g. by Kimura (1965), Lande (1975 and Turelli (1984), as stationary solution at stabilizing selection. However, stabilizing selection is hardly relevant to description of intermediate states of continuing evolution; hence as far as we can see, the question whether and when Gaussian distribution may actually realize during such evolution remains open.
Here we aim to analyze the speed of evolution while far from any evolutionary stable state, based on the simplest possible meaningful model. This is a deterministic integro-partial differential equation, which is similar to various forms postulated or derived elsewhere. We also provide a simple derivation of this model from "first principles", avoiding to make non-verifiable assumptions, for fear that the ultimate results may become artefacts of any such assumptions. We consider single asexually reproducing species, stick with phenotypic description and use dynamic formalism, leaving all stochastics within the underlying population dynamics model of intraspecific competition. As a model of intraspecific competition, we consider a predator-prey system where various predator populations depend on a common prey and do not interact otherwise. All these assumptions are admittedly rather restrictive; we believe that ab initio approach not only helps to make clear what postulates the ultimate results depend on, but can also suggest a basis for subsequent generalizations for more realistic assumptions. The simplicity of the resulting mathematical model allows its exhaustive rigorous study. In particular, the Gaussian shape of the trait distribution does indeed emerge spontaneously. The practical utility of the model is illustrated by providing treatment of a more realistic model via asymptotic methods.
The structure of the paper is as follows. Section 2 introduces the main equation and its variations. Section 3 describes the "normal solution" of the simplified version of the model, which underlies subsequent analysis. Section 4 goes on to consider the general solution of the simplified model. These results are extended to the more generic version of the model by means of a perturbation theory in Section 5. Section 6 is dedicated to the discussion of the results. We conclude with Appendix A with the derivation of our model "from the first principles" and Appendix B with the proof of the theorem presented in Section 4.

The model
We consider diffusive, mutational spread of a phenotypic trait distribution during a transient phase of optimizing evolution, described by a deterministic model of the form where x is a continuous trait, u(x, t) ≥ 0 is population density in the trait space so that N (t) > 0 is the total population at time t, r(x) is the fitness of trait value x measured as the low-density reproduction rate of the given type and corresponding carrying capacity of the habitat, and C(x) and D(x) > 0 represent mutations, with C(x) for the directed component and D(x) for the diffusive component. The case C(x) = D ′ (x) of this equation can be obtained as a special case of deterministic integro-partial differential equation (4.5) derived in Champagnat et al (2006) as a weak limit of a stochastic model, with convolution kernels U = const and V = const, (the "mean-field" case), and with an appropriate choice of functions b and d. A simple non-rigorous derivation of (2) through a continuous-trait limit of a deterministic population dynamics model is given in Appendix A. The model is admittedly rather simplified and ignores many important evolutionary factors, e.g. frequency dependent selection. Some equivalent forms of special cases of this equation found in literature will be mentioned below. We shall first look at solutions of (2) in simplifying assumptions regarding its coefficients, and then relax those assumptions by means of a perturbation theory. So, in the simplified version, we take that D(x) = D = const (which is a significant limitation as normally one would expect that mutation rate is proportional to birth rate so should be varying together with the reproduction rate r(x)). Further, the coefficient C(x) represents possible mutations bias with respect to the selected trait. In many studies it is assumed to be zero, however it may correspond to nonselective evolution, discussed e.g. by Koonin (2009). For the sake of simplicity, we take C = const; then without loss of generality we take C = 0, as a nonzero constant C would simply add −C to the trait change rate. So, Function r(x) plays the role of the relative fitness of the subpopulation with trait value x. Considering this function a constant would not be appropriate as it would remove any selection; in our simplified version we take it to be a linear function r(x) = r 0 + r 1 x, where r 0 , r 1 = const. This results in the equation which is essentially identical e.g. to equation (2.1) postulated by Calsina and Perelló (1995). Substitution brings the evolution equation (3) to the well known (e.g. Taylor and Jonker, 1978;Page and Nowak, 2002;Karev, 2010) "replicator-mutator" form, Note that by virtue of (5) and the definition of N in (3), we automatically have at any time. Function p(x, t) is the probability density function (PDF) of the population in the trait space x at time t, andr(t) is the mean fitness of the population at that time. If solution of the replicator-mutator equation (6) is known, then the total population size can be found by solving equation which then recovers a solution to the original evolution equation (3) via (5).
brings equation (3) to the form ∂v ∂t provided that the factor P satisfieṡ Hence the integral part of the equation uncouples from the differential, the closed linear differential equation (10) can be solved first, and the linearizing factor can be found afterwards via (12) which, for the initial condition P (0) = 1, gives the explicit expression (this is continuous-trait variant of linearization used by Schuster (2011)). In ecological terms, the linear equation (10) corresponds to the case when the subpopulations with different traits x are in no direct competition with each other or even within themselves and, aside from mutations described by the term Dvxx, each subpopulations multiplies by a Malthusian law with its own reproduction coefficient r(x). For the linear fitness function r(x), we have the alternative forms of (4) respectively as ∂p ∂t and ∂v ∂t In the next two sections, we concentrate on the solution of equation (4) and its quivalent forms (13) and (14), before relaxing the simplifying assumptions in Section 5.

The normal solution
We look for solutions that describe the dynamic change of the population during its gradual adaptation, while far from any evolutionary steady state. We note that equation (13) admits a family of exact self-similar solutions, which are Gaussians, or PDFs of normal distributions, As is easily verified by direct substitution, function (15) is a solution of equation (13), provided that the parameters of the normal distribution obey the following system of ordinary differential equation: We shall call (15) a normal solution of (4). We stress here that functional form (15) is not an arbitrary assumption, but an exact consequence of the evolution equation (4), once appropriate initial conditions are supplied. These initial conditions should, of course, be Gaussian. However, numerical simulations shown in fig. 1 suggest that the general solution at arbitrary initial distributions asymptotically become normal as time increases, so the special class (15) should in fact be fully representative. This is indeed the case as we show below.

General solution
We formulate properties of the solution to (13) for a wide class of initial conditions, which generalize the properties of the normal solutions shown above. We do that in terms of the moments of the function p(x, t) considered as PDF of the random quantity x at a given moment of time t. Namely, we use the mean, and the variance Theorem 1 Let p(x, t) be a solution of (13) with initial condition p(x, 0) = p 0 (x).
We assume that that the initial condition satisfies Then this solution can be written in the form where w(·, t) is a PDF of a zero-mean, unit-variance distribution, and for all t ≥ 0, (C1) w(·, t) has moments of all orders, the mean of PDF p(·, t) varies according to (C4) and its variance according to where The proof is given in Appendix B. Conclusion (C3) states that equation (16a) remains exact in the general case. This is, of course, a special case of the generic Fisher-Price law (Fisher, 1930;Price, 1972) as applied to the current model.
On the contrary, conclusion (C4) states that equation (16b) is not exact and in the general case requires a correction associated with the instant value of skewness γ of the distribution. However, (16b) is "asymptotically valid" for large t, as the skewness, according to (C2) vanishes in the limit t → ∞. Moreover, from the proof we can see that the asymptotic order of γ is such that the second term in (C4) is asymptotically smaller than the first term in the limit t → ∞.

Perturbation theory
The previous results were for a simplified version of the model, where dependencies C(x), D(x) were replaced by constants and r(x) was replaced by a linear function, to allow a simple analytical solution. Now we would like to ensure that these solutions are not exceptional and small violation of the simplifying assumptions will not lead to completely different solutions. So, we now consider generic smooth dependencies for C(x), D(x) and r(x), but assume that the variation of x across the population at any given time is smaller than the typical scale at which these functions vary significantly. So, we develop a perturbation theory where the small parameter is σ, the standard deviation of the selected trait in the population. Admittedly σ is not a constant but a dynamic variable; this however does not affect the formal asymptotic expansions. Equivalently one could use the inital value σ 0 as the small parameter. This, however, complicates notation, so we do it only in one place where such explicit treatment is essential. We require that functions C(x), D(x) and r(x) and the necessary number of their derivatives are bounded and that D(x) is everywhere bigger than some positive constant, so that |D ′ (x)/D(x)| is bounded. We silently assume that all the moments are at most O (1), before establishing their actual asymptotic orders more accurately. All results in this section are obtained formally, without any attempts of rigorous justification.
We start from the integro-partial differential equation (2) and bring it to the form of replicator-mutator equation by substitution which gives Now we pass from the probability density p(x, t) to the normalized (zero-mean, unit-variance) probability density w, through p's mean and its variance via substitution (for chain rule differerentiation, it is convenient to distinguish the time variables in the old coordinates (x, t) and the new coordinates (η, τ ); for functions of one variable such as σ and µ this distinction is of course not important). Then (27) is subject to constraints Let us expand functions of x in Taylor series, and consider the formal asymptotic expansion of equation (27) in the small parameter σ.
It is straightforward to see that if w dη = 1 at τ = 0, it remains so for all τ > 0, so constraint (28) is always satisfied. The constraints (29) and (30) lead to asymptotic series forμ andσ. Further, multiplying both sides of equation (27) by η n , n ≥ 3, and integrating over η leads to asymptotic series for the moments Mn. In this way, we obtaiṅ So, in the selected asymptotic orders, of all the moments only the skewness γ = M 3 affects the evolution speed. To see how big is its effect, we now need to consider σ = σ(t) as a function of time rather than merely a small parameter. To estimate the upper bound for γ(t), we consider γ as a function of σ, with the initial condition where the integrand D 1 /D 0 depends on σ via µ = µ(t) and t = t(σ). Let us assume that |D 1 /D 0 | ≤ G = const for the whole solution under consideration. Then it follows that So, our upper bound for skewness γ consists of two components: one is related to the decaying contribution of the initial condition, |γ 0 (σ 0 /σ) 3 | < |γ 0 |, and γ 0 (σ 0 /σ) 3 = O t −3/2 → 0 as t → ∞, and the other is the contribution of the perturbation, which is O (σ). So if the effect of the initial skewness can be neglected, say γ 0 = O (σ), then we have γ = O (σ), and from (31) we have finally our asymptotic evolution equationṡ This (asymptotically) closed system of ordinary differential is a generalization of the previous result (16) and it asserts that that simple system of two equations remains approximately valid even if fitness gradient and mutation diffusivity are not constant, only subject to the drift term due to mutation bias −C(µ), as we discussed above. Moreover, it provides higher-order corrections to that simple system, if necessary. Further increase in asymptotic accuracy will require including into consideration higher-order moments. Let us consider higher-order asymptotic equations for a special case when D(x) = D = const, C(x) = C = const so only the fitness r(x) is trait-dependent. Then reasoning as before, we see that in this case γ = O σ 3 , M 4 = 3 + O σ 3 , and thereforė The third term in the right-hand side of (33)(b) describes the effect on the variance of the stabilizing (r ′′ < 0) or disruptive (r ′′ > 0) selection. Note that at a fitness maximum, r ′ = 0, r ′′ < 0, equation (33)(b) gives a stationary variance at σ 2 = −2D/r ′′ (µ) 1/2 , in agreement with the classical result by Kimura (1965).

Discussion
We have considered solutions of a simple model of gradual Darwinian evolution in continuous phenotypicl trait space and continuous time, while far away from evolutionary stable equilibrium. In biological terms, the main predictions of the model are: -Equation (8) states that total population size N is described by Verhulst dynamics, where the mean reproduction rater(t) is itself dynamically changing. -Equations (16a), (20), (32)(a), (33)(a) describe the dynamics of the mean trait value µ(t), which has a selective component, proportional to the local fitness gradient r ′ (µ) and instant distribution variance σ 2 (t), and a non-selective component, due to mutation bias C(µ). The selective component corresponds to the canonical equation (1) of the adaptive dynamics and to Fisher's Fundamental Theorem of Natural Selection. That is, not only population evolves faster in a stronger fitness gradient, but also populations which are more diverse with respect to the selected trait, evolve faster, and those that are very homogeneous in the selected trait, evolve slower. The Gaussian distribution is preserved even if the gradient of fitness and the mutation parameters vary with the change of the trait coordinate, provided that their variation is slight within the spread of the population.
Since the rate of evolution is proportional to the variance, and the variance in the long run is constantly growing, these result predict that the evolution has a propensity to accelerate, for as long as the assumptions of the model remain valid. This property is illustrated by a simulation presented in fig. 2. Our underlying assumptions are self-limiting. If coefficients C, D and r ′ are varying significantly on the trait space scale of xc, then our perturbation theory will remain valid only for as long as the variance is small enough, σ 2 ≪ x 2 c . On the other hand, this variance grows unbounded, σ 2 ≈ 2Dt, hence this theory may be only valid for a limited time interval, t ≪ x 2 c /(2D). However, this limitation is consistent with our goal which was to consider the transient phase of the evolution, before the population reached the optimum, as stationary distribution around the optimum is considered in other works.
To conclude, we discuss the role of mutations, specifically the value of coefficient D in the above results. According to equation (20), the instantaneous rate of change of the average value of the trait does not depend directly on D, which is entirely consistent with the general Fisher-Price law which claims that it is determined only by intensity of selection and variance of the instant trait distribution. As to the variance itself, then according to (21) in the long run it always grows with the rate directly proportional to D. The convergence of the distribution shape to Gaussian also depends on mutation. According to estimates obtained in Appendix B, γ(t) ≈ Q/(r 3 1 D 3/2 t 9/2 ), t → ∞, where Q is a constant depending on initial conditions. Then, accepting smallness of skewness γ as a measure of convergence to Gaussian distribution, we see that this convergence is due to both selection (r 1 ) and mutation (D) and so will be very slow in the limit of r 1 → 0 or D → 0. So, for a fixed time interval and very near to the fitness extremum, or for very small mutation rate, our results become inapplicable. This is of little surprise as mutation without selection as well as selection without mutation are completely different biological situations as well as completely different mathematical problems. The case of small mutation rate is considered in a recent paper by Sasaki and Dieckmann (2011), who start from selection without mutation and then add mutation as a perturbation, thus presenting an approach which is complementary to the one given here, and more appropriate for analysis near fitness extrema, as opposed to nearly-linear fitness landscapes considered here. Many if not all individual ideas of our derivation are found in literature; however we could not find such derivation as a whole, so present it here in its entirety, step by step, emphasizing all assumptions made on the way, in order to identify the limitations of the resulting model, which may lead to ways to overcome these limitations. We start from the population dynamics model of Lotka-Volterra-Gause type, where the predator population consists of a number of subpopulations z j , differing in their parameters, all consuming the same prey (resource) R, mutating into each other, and not interacting otherwise: where α is the low-density reproduction rate of the resource in absence of grazing pressure, β j describes the per capita grazing pressure by trait j on the common resource, the coefficient at the quadratic term in the first equation of (34) is unity due to the choice of the unit of measurement for R, (−γ j + δ j R) is the resource-dependent reproduction rate of trait j, and m j is the contribution of mutations, defined as Here we chose index j so it enumerates subpopulations monotonically with respect to the selected trait. The coefficients m ± j are the rates of mutation of subpopulation j that respectively increase or decrease the trait index j. So, in this model different types differ both by internal dynamics and by their interaction with the common resource, and their competition is indirect, only via dependence on the common resource. The population model with mutations (34,35) is admittedly very simple and ignores many biological aspects which can be very important in real life. For instance, it ignores frequency dependent selection, as the population dynamics of each type does not depend on its abundance relative to other types, but only on the "mean field", represented by the common resource. Also, we shall consider "smooth" dependence of the coefficients in this model on the one-dimensional index j, which will effectively impose one-dimensionality on the resulting model, which is also a significant restricting assumption, see Metz et al (2008). We stress however that our purpose here is not to derive a biologically realistic model, but only to present a simple example of how equations of evolution can be derived from population dynamics equations.
For mathematical simplicity, we take the set of subpopulations {j} to be infinite, which in practice only means that the overall diversity of the population observed during the stage of evolution under consideration is much smaller than all that are theoretically possible. Also for mathematical simplicity only, we take that mutations are so small they occur only between subpopulations adjacent in the trait index j; this assumption is not essential: admitting mutations to any small distance in j leads to the same form of dynamic equation in the continuous limit. What is important is that we neglect the probability of mutations that change the trait significantly, in comparison with the diversity of the population at any given moment of time; taking those into account would make the resulting equation a bit more complicated (mutations will be described by an integral rather than differential term, see e.g. Champagnat (2006)). There are of course also lethal mutations but they can be considered as contributing to the death rate and thus implicitly included in γ j .
Our next assumption is that the dynamics of resource (prey) are much faster than those of predators, so the resource concentration can be adiabatically eliminated. We find the quasistationary resource concentration as which gives the predators dynamics equations as We change the dynamic variables to measure the subpopulation sizes by their grazing pressure, u j = β j z j , and turn to the continuous limit via x = jh, h → 0, where x is a physical measurement of the trait. This leads to the integro-differential equation where and all limits are taken as h → 0, j = x/h, x = const. In particular, we assume that dependence of mutation rate m ± j on j and the asymmetry between beneficial and deleterious mutation rates transforms equation (37) to Finally, in this paper, we consider the case δ(x) = const, and we rescale population density so this constant is unity, which gives the model (2). In biological terms, this means that we assume that the competition between different types is in the speed of reproduction which is assumed strictly proportional to the carrying capacity, i.e. the equilibrium population density of trait x in the given habitat in absence of other types. This choice is arbitrary and is made from consideration of mathematical simplicity rather than motivated by any specific biological examples. We note however, that the effect of this, as well as many other simplifying assumptions made above, can be relaxed by perturbation theory, of the kind considered in Section 5.

B Proof of Theorem 1
First two brief preliminary remarks.
-We shall use the identity (7) in our proof. As already noted, this is ensured by construction of the equation, and it also can be easily verified that once this property is satisfied for the initial condition according to (A1), it is then preserved by the equation (6). -We shall, without loss of generality, assume that x + = 0; otherwise, we just make trans- Note also that some of the variables used here coincide by name with variables used in Appendix A but have different meaning (m, δ).
Proof 1 • We start by proving (C1). We do that using the equivalent linear equation (14). First, we construct the fundamental solution of that equation, that is, a generalized solution V (x, t ; ξ) with initial condition V (x, 0 ; ξ) = δ(x − ξ), where δ() is Dirac's delta-function. This can be obtained from the normal solution p * (x, t ; µ, σ) (15,16) with initial conditions µ(0) = ξ and σ(0) = 0, transformed by (9) and (12), which leads to Hence for the generic initial condition v( where K 1 = e r 0 t+Dr 2 1 t 3 /3 2 √ πDt ,σ 2 = 2Dt,μ = r 1 Dt 2 are some known functions of time. Normalization of this solution gives the PDF p(·, t), the moments of which are found as The corresponding double integrals are absolutely convergent under the assumed properties of p 0 . Furthermore,Ĩ 0 > 0. Hence existence of moments of all orders at all t for the PDF p follows. In particular, the mean µ(t) =M 1 (t) and variance σ 2 (t) =M 2 (t) are defined for all t ≥ 0. The normalized PDF w(η, t) is obtained from p(x, t) via substitution x = µ + ση, and is zero-mean and unit-variance by construction. Its moments of all orders exist as they are the standardized moments of PDF p. This is the conclusion (C1) of the Theorem. 2 • Now that the existence of µ(t) and σ(t) for all t ≥ 0 is established, we can proceed to prove the predictions (C3) and (C4) about their dynamics. To find the rate of change of the mean µ(t), let us first note that for a linear fitness function r(x), we havē Now let us multiply both sides of (13) by x and integrate them over x ∈ R. This gives where we used the definition of the mean (17) and the integral proprotional to D vanishes when integrated by parts (the limit lim x→±∞ [x∂xp(x, t)] = 0 can be verified using (40)). Now we do equivalent transformations in this equation, where we have used definitions (7), (17) and (18). According to (41), the last two terms cancel out, which delivers the conclusion (C3) of the theorem.
Transformation (19) ensures some (already mentioned) identities for w(η, t), which we will now need stated explicitly. Namely, follows from (7), follows from (17) and follows from (18). Substitution of (19) into the linear evolution equation (13), with account of (42-44) and the already proved identity (C3), leads to the following differential equation Considering the first three moments of both sides of this equation, we see that the subspace of functions defined by (42-44) is an invariant set of this equation if and only iḟ which is conclusion (C4) of the Theorem. 3 • It remains to prove (C2). We do it by the method of moments. An exact expression for the normalized PDF w(η, t) is obtained from (40) via substitution x = µ + ση, where µ = µ(t) and σ = σ(t) are the true mean and standard deviation of the v-distribution (as opposed to the "rough guess" values of the same,μ(t) andσ(t)), and an appropriate normalization. This gives where K 2 is a coefficient depending only on t but not on η, chosen so that ∞ −∞ w dη = 1.
The moments of PDF w are found as The already mentioned absolute convergence of double integrals (47) means that Fubini's theorem applies and we can change the order of integration. On doing so, and also introducing notations where φ is the initial PDF, "biased" by selection towards fitter trait values: φ(ξ, t) = p 0 (ξ)e r 1 tξ .
For brevity, we shall now omit dependence on time, until we start studying the t → ∞ asymptotics. By using the binomial formula, the moment integrals are rewritten as are the moments of the standard normal distribution, and where, in turn, Φ ℓ are the moments of the biased initial condition, Combining these together, we obtain where [x] denotes the integer part of x. In particular, hence for the moments we have These expressions can be used for determining the true mean and variance in terms of the initial PDF, via the identities (43) and (44). Namely, we have Then, With account of these, finally we have an exact formula for the moments of w in terms of the initial PDF p 0 : So, the problem of the t → ∞ asymptotics of the moments is reduced to asymptotics of Φn. In this limit, we haveσ = (2Dt) 1/2 → ∞ andμ = r 1 Dt 2 → ∞. We also introduce s = r 1 t for brevity, and s → ∞, too.
In terms of s, integrals Φn defined by (49) are, up to the signs, the bilateral Laplace image of the initial distribution p 0 for n = 0, and derivatives of that image for n > 0. Asymptotics of Laplace images are known to be very sensitive to analytical properties of the originals. So at this point our assumptions (A2) and (A3) become essential. In accordance with the second preliminary remark, we set x + = 0 and x − = −W , W > 0, without loss of generality. Then the initial PDF has the form: Then for the biased moments we have asymptotic series am(−1) m+ℓ (m + ℓ)! s m+ℓ+1 + e.s.t.
Summarising, we have for all n that where An are the moments of the standard normal distribution, (48). Since the normal distribution is uniquely characterized by its moments, convergence of w in moments here implies convergence of distributions, and we have the claimed (C2). ⊓ ⊔ Two final remarks.
-Equations (55) and (56), meaning that the "crude guesses"μ andσ happen to be asymptotically accurate for the true µ and σ, are a direct consequence of the choice x + = 0, as can be verified by repeating the calculations with generic x + . There is a simple interpretation of this fact in biological terms. In the initial PDF p 0 (x) there is only a finite range x ∈ [x − , x + ] of traits present, and the descendants of individuals with different traits make different contributions to the overall population at different times; but in the present model, the distribution of traits in the population in the long run is such as if descendants of the fittest ancestors, x = x + , dominated in it, regardless the effect of mutations. Mathematically, this is eventually down to linearity of (10). With initial conditions that are not finite-supported, the problem becomes significantly more complicated: the "dominant ancestor" trait will keep changing with time. -According to (59), we have γ = O t −9/2 , so the related correction 1 2 r 1 σ 2 γ in (C4) is O t −7/2 , which is asymptotically smaller than the main term D/σ, which is O t −1/2 .