A stochastic model for the evolution of species with random fitness

We generalize the evolution model introduced by Guiol, Machado and Schinazi (2010). In our model at odd times a random number X of species is created. Each species is endowed with a random fitness with arbitrary distribution on $[0, 1]$. At even times a random number Y of species is removed, killing the species with lower fitness. We show that there is a critical fitness $f_c$ below which the number of species hits zero i.o. and above of which this number goes to infinity. We prove uniform convergence for the distribution of surviving species and describe the phenomena which could not be observed in previous works with uniformly distributed fitness.


Introduction
During the history of our planet, species have emerged and have become extinct, some have lasted a relatively brief period, others are still present in a more or less unchanged form after millions of years. It is widely accepted that the driving engine of evolution is natural selection or "survival of the fittest". It is therefore interesting to provide mathematical models for the evolution of species.
Guiol, Machado and Schinazi [6] proposed a model where creation and deletion of species is driven by chance in the sense that at each step with probability p one new species is created and its fitness is chosen uniformly in [0, 1], while with probability 1 − p the least fit species (if there are species alive at that time) is removed. One motivation for the study of this model is that its long-term behaviour is similar to the one which simulations show for the Bak-Sneppen model: there is a critical value for the fitness and species with smaller fitness disappear, while species will a larger fitness persist indefinitely. Bak and Sneppen [1] modelled a simple ecosystem where the population size is constant and at each step not only the least fit is removed, but also its neighbours are replaced by new species (proximity may be seen as representing ecological links between species). It has proven difficult to obtain rigorous results for this model (see for instance [8]) and this motivates the search for similar, more tractable models.
Several papers have studied the GMS model: [3] gives a law of the iterated logarithm and a central limit theorem for number of species with supercritical fitness which go extinct (this number is negligible with respect to n); [5] studies the maximal fitness ever appeared in the subcritical case.
The model has been generalized in [9] and [2]: there is still a toss of a coin to decide for creation or deletion, but instead of adding/removing one species at a time, increments are arbitrary random variables. Even with these assumptions, the same cut-off phenomenon of [6] appears.
In the original GMS, the lengths of subsequent births and deaths are geometrically distributed random variables (with parameters which sum up to 1) and in [2,9] they are geometrical convolutions of certain laws (where the parameters of these geometrically distributed number of convolutions, again, sum up to 1). In our model we group all subsequent creations and deletions: the length of subsequent creations {X n } n∈N and the length of subsequent annihilations {Y n } n∈N are such that {(X n , Y n )} n∈N is an i.i.d. sequence with arbitrary distribution. Whence our results apply to the models in [2,6,9] (see Section 2.1). Besides, in the older papers the fitness is assigned uniformly while we use a general distribution µ. If µ has atoms, a new phenomenon appears: there might be a 1 fitness which acts as a barrier eventually protecting all species with higher fitness (see Corollary 2.4 and the subsequent discussion for details).
Here is the outline of the paper. In Section 2 we give the formal construction of the process and the necessary definitions. We state our main result, Theorem 2.2, which describe the asymptotic expression of the proportion of species in a generic (Borel) range of fitness. The asymptotic behavior of a single fitness is described by Theorem 2.3. Corollary 2.4 and the subsequent discussion gives some details on the number of species which are killed. Section 2.1 is devoted to a detailed comparison with previous works; we explain why our work is a generalization of the previous models and which new phenomena arise.
In Section 3 we study some examples: the original GMS with an atomic measure µ (see Section 3.1), a Markov model which cannot be treated by using previously known results (see Section 3.2) and a model which is related to Branching Processes (see Section 3.3). We also give a counterexample to be compared with Theorem 2.2 (2).
All the proofs are in Section 4 which contains a couple of results which are worth mentioning: a Law of Large Numbers (Proposition 4.2) and Proposition 4.1 which identifies the set of fitness which become empty i.o. (and the total amount of time they are empty).

The process and its asymptotic behaviour
We start by giving a formal description of the process. Let {X n , Y n , f n,i } n,i∈N be a family of nonnegative random variables and, for all n ∈ N, denote by f n the sequence {f n,i } i∈N . Suppose that (1) for every n ∈ N, (X n , Y n ) and f n are independent, (3) all f n,i are distributed according to a measure µ on R.
Roughly speaking, X n counts the new species at time n, Y n counts the deaths and f n,i the fitness of a newly created species. In order to avoid trivial cases we suppose that E[X k ] and E[Y k ] are both in (0, +∞]; moreover we assume that at least one of these two expected values is finite. Note that in this case {X n } n∈N , {Y n } n∈N and {X n − Y n } n∈N are all i.i.d. families, but X i and Y i might be dependent. From now on, we will denote by (X, Y ) a couple with the same law as (X 1 , Y 1 ). For every fixed n ∈ N also {f n,i } i∈N might be dependent (for instance they can be generated by a Markov Chain or f n,1 = f n,i for all i ∈ N).
We will assume that µ([0, 1]) = 1; there is no loss of generality, since any measure on R can be mapped to a measure supported in [0, 1]. We denote its cumulative distribution function by F = F µ and we define F (f − ) := lim a→f − F (a).
Let Z n be the number of species alive at time n. We start at time 0 with Z 0 = 0 (Z 0 could be a random variable with an arbitrary distribution on N).
At time 1, X 1 species are generated and to each of them we assign a random fitness with law µ. More precisely the fitness of the i-th created species is f 1,i for all 1 ≤ i ≤ X 1 . Thus Z 1 = Z 0 + X 1 . The procedure is repeated at any odd time: Z 2k+1 = Z 2k + X k+1 , meaning that X k+1 species are created and their fitness f k+1,1 , . . . , f k+1,X k+1 are assigned. For any set A ⊆ [0, 1] we denote by Z n (A) the total number of species alive at time n and with fitness in A. The fitness of a species does not change during its entire lifetime, and species may disappear only at even times.
At time 2k + 2, a number Y k+1 ∧ Z 2k+1 of species are removed and removal starts from the least fit. This means that All species with fitness not larger than x − are removed and Z 2k+2 (A) = 0 for all A ⊆ [0, x + ). A number M k+1 := Y k+1 − Z 2k+1 ([0, x − ]) of species is removed from the set of species with fitness equal to x + : Given a Borel set A ⊆ [0, 1] such that µ(A) > 0, we define the number of species created in A as 1l A (f n,i ).
(2.1) By our assumptions, for any A, we have that {( X n , Y n )} n∈N are i.i.d. and E[ X n ] = µ(A)E[X n ] (where a · (+∞) = +∞, if a > 0 and 0 · ∞ = 0). Henceforth, an interval I ⊆ [0, 1] (either closed or not) such that 0 ∈ I is called a left interval. We note that, for a left interval I such that µ(I) > 0, {Z 2n (I)} n∈N is the queuing process (see [4,Chapter VI.9]) associated to the i.i.d. increments { X n − Y n } n∈N (see Section 4 for details). We will often make use of the expected value We define the critical parameter:  It is a consequence of the following theorem that, when A ⊆ [0, 1] is a Borel set, either there is extinction in A or there is survival. Indeed, if there is no extinction in A then P ∞ (A) > 0, thus Z n (A) → +∞ almost surely. By a standard argument this implies survival.
The example given in Section 3.4 shows that, if {f n,i } i∈N are just dependent, then the conclusion in Theorem 2.2(2) might be false.
The following theorem describes the long-term behaviour of a fixed fitness. Note that all f > f c belong to case (1), while all f < f c belong to (2). If f = f c , then case (2) applies if and only if ( where X n is the number of species created in the Borel set A (see equation (2.1)).
Here is a more explicit description. First of all, by (5) a.s. there are no more species killed in [f, 1] eventually as n → +∞ but by (2) the number of species killed in [f c , 1] diverges almost surely. If , so that by (4) K n ([f c , 1])/n goes to zero a.s. and τ n ([0, f c ])/n → 0 almost surely as n → +∞ (see Proposition 4.1(2)).
If µ({f c }) > 0 then we have the following possibilities: then, just as before, a.s. the species killed in [f c , 1] eventually will have fitness f c and the fraction of species alive with fitness f c converges to the same positive limit. This time K n ({f c })/n has a positive limit: −E[F (f − c )X − Y ]/2 and τ n ([0, f c ))/n converges to a positive limit almost surely as n → +∞ (see Proposition 4.1(3)). • then, by Theorem 2.2(2), every species with fitness f c is eventually killed a.s. and )/n tends to 0, a.s., thus outside a negligible proportion, the killed species all have fitness f c , whence K n ({f c })/n has the same positive limit as before. Finally, τ n ([0, f c ])/n → 0 almost surely as n → +∞ (see Proposition 4.1(2)).

Comparison with previous works.
Our process extends those appeared in [2,6,9]. Aside from our general choice for the fitness law, the birth-and-death mechanism that we study is more general than those adopted in these papers.
One way to see the original GMS (see [6]) as a particular case of our process is by observing that the random sequences of consecutive births X k and consecutive deaths Y k have right-shifted Geometric distribution with parameter 1 − p and p respectively.
In general, consider a process {Z n } n∈N where at each step either a species is created (along with its fitness) or the least-fit species, if any, is removed. Denote by X 1 > 0 the length of the first stretch of "creations", followed by a stretch of "annihilations" of length Y 1 > 0, then another stretch of "creations" of length X 2 followed by a stretch of Y 2 "annihilations" and so on. Suppose that {X n } n∈N and {Y n } n∈N are two i.i.d. sequences. It is clear that there is a connection between our process and this one, namely for every set Therefore, the long-term behaviour of {Z n (I)} n∈N can be derived simply by studying {Z n (I)} n∈N . Our work can also be considered as a generalization of [2] and [9] whose models are essentially equivalent. Indeed, in [2], a single family of Z-valued variables {U n } n∈N is considered. In this process, U n > 0 means that U n species are created, while U n < 0 means that −U n species are killed. In this case the laws of length of a "creation" stretch X i and "annihilation" stretch Y i are necessarily geometric random convolutions of the law of U n conditioned on {U n > 0} and {U n < 0} respectively. Moreover, the sum of the parameters of these geometric convolutions must be 1 − P(U 1 = 0). Therefore, a model constructed from the variables {U n } n∈N can be considered as a particular case of our model: take for instance X n := U n 1l {Un>0} , Y n := −U n 1l {Un<0} and consider the process {Z 2n (A)} n∈N . In Section 3.2 we consider a particular case of our process which cannot be obtained with a single family of variables describing simultaneously creations and annihilations.
Observe that in Theorem 2.2 we used Z n as a normalizing factor for Z n (A) but there are two other natural choices: n (to compare with [2,9]) and N n (to compare with [3,6]).
as n → +∞. Hence Theorem 2.2(3) can be equivalently written in terms of the timescale n or N n (in this last case we obtain a generalization of Proposition 4.2(1) to Borel sets).
then Z n ∼ N n almost surely as n → +∞. Indeed one can use the same kind of arguments used in the proof of Theorem 2.2(2), to prove that Z n and ⌊(n+1)/2⌋ i=1 X i are asymptotic and the remaining terms are negligible. Roughly speaking, changing timescale turns out to be just a linear rescaling.
We note that for the GMS model and its generalizations, with µ ∼ U([0, 1]) (where U(I) is the uniform distribution on I), the fraction of surviving species in any I ⊆ [f c , 1] is proportional to µ(I). This is still true in our case when I ⊆ (f c , 1], but it does not hold for 1])/n → 0 (the exact rate of convergence for the GMS is studied in [3]), while again this needs not to be true if f c is an atom for µ.

Examples and counterexamples
3.1. The original GMS model. The original GMS process can be seen as the particular case where X has a geometric law G(1−p) while Y has a geometric law G(p). Thus E[ /p is the relevant term for computing f c according to equation (2.2).
In this example we choose µ := αδ 1/2 + (1 − α)ν (where ν ∼ U([0, 1])); the case α = 0 is discussed in [6]. Roughly speaking, every time a new species is born we toss a (possibly biased) coin: with probability α we assign to the new species a fitness 1/2 and with probability 1 − α the fitness is drawn uniformly and independently in [0, 1]. Clearly where f c is given by equation (2.2) and it is the unique solution in . More interesting is the cumulative limit distribution (see equation (2.3)) To avoid useless complications, we discuss just the "fair coin" case α = 1/2. In this case we have There are five typical situations that we can explore and they are represented by the following table where we choose (1 − p)/p = 7/8, 3/4, 1/2, 1/4, 1/8.
starting from a birth. Thus the probability of a birth after the birth P ++ = p, the probability of death after the birth is P +− = 1 − p and so on. This can be seen as a particular case of our process where X has a geometric law G(1 − p) while Y has a geometric law G(1 − q). We assume that p > q; As before we choose µ := αδ 1/2 + (1 − α)ν (where ν ∼ U([0, 1])); thus the cumulative distribution function is still given by equation (3.6). q) and, according to equation (2.3), . As before, we discuss just the "fair coin" case α = 1/2. In this case we have We retrieve the same typical cases as before by choosing (1 − p)/(1 − q) = 7/8, 3/4, 1/2, 1/4, 1/8. 3.3. The Branching Process case. In this example we consider the case where P(Y n = 1) = 1 while X k has a generic discrete distribution on N. In order to avoid a trivial situation we assume that P(X k = 1) < 1. In the following, we make use of the generating function of the variables {X n } n∈N , that is, Φ(z) := ∞ n=0 P(X 0 = n)z n for all z ∈ [0, 1]. Similarly, the generating function of the random number of species whose fitness belongs to The peculiarity of this case is the fact that the process can studied by means of a branching process and the probability of survival of a fitness can be computed in terms of the probability of survival of the branching process. 3.4. Counterexample of Theorem 2.2(2) for dependent {f n,i } i∈N . We define Y n := 1 a.s. and f n,i := f n,1 for all i > 1 and n ∈ N (where {f n,1 } n∈N is an i.i.d. sequence distributed according to µ). We construct the sequence {X n } n∈N as X n := g(H n ), for a suitable choice of an i.i.d. sequence {H n } n∈N and a function g.
Let us define T k := min{i : H i > n k } for all k ∈ N (clearly T 0 = 1); note that T k ∼ G(1/2 n k ). Since lim m→+∞ P(T k ≤ m, H T k ≤ n k+1 ) = P(H T k ≤ n k+1 ) = 1 − 1/2 k+1 then there exists τ k ∈ N such that P(T k ≤ τ k , H T k ≤ n k+1 ) ≥ 1 − 1/2 k . The sequence {τ k } k≥1 can be always constructed iteratively as a nondecreasing sequence. It is not difficult to prove, by using the Borel-Cantelli Lemma, that the event Ω 0 := k∈N {H T k ≤ n k+1 , T k ≤ τ k } has positive probability.
We are now ready to define g(i) := (k + 1)! k j=0 τ j for all i = n k + 1, . . . , n k+1 (for all k ∈ N); X n := g(H n ) for all n ∈ N. On Ω 0 we have Roughly speaking, this means that, on Ω 0 , for every ε > 0 infinitely often the last generation represents at least a fraction 1 − ε of the entire population. Whence due to our choice of {f n,i } i∈N , on Ω 0 , for every Borel set A ⊆ [0, 1] such that µ(A) > 0 and for every ε > 0, we have Z n (A)/Z n ≥ 1−ε infinitely often.

Proofs
We note that, for any fixed left interval I, {Z 2n (I)} n∈N is a random walk on N (with increments depending on the position). More precisely it is the queuing process associated to the i.i.d. increments { X n − Y n } n∈N , as defined by equation (2.1); indeed We denote by {S n (I)} n∈N the random walk with independent increments, where S n (I) : At time 0 we have S 0 = Z 0 (I) = 0 and for all n  Note that the case where P( X = Y ) = 1 is trivial, since it means that µ(I) = 1 and X = Y = c ∈ (0, +∞) a.s.; thus, Z n (I) equals c when n is odd and 0 when n is even.
Proof. Recall the relation between the random walks {S n (I)} n∈N and {Z 2n (I)} n∈N given by equation (4.9). In particular the return times to 0 of the second process are the weak descending ladder times of the first one, that is, the times n such that S n (I) ≤ S i (I) for all i ≤ n. We denote by {T i } i∈N the sequence of intervals between two consecutive weak descending ladder times of {S n (I)} n∈N (that is, the times between two consecutive returns at 0 of {Z 2n (I)} n∈N ). Note that {T i } i∈N are i.i.d random variables and T 1 = min{n ≥ 1 : S n ≤ 0}.
(1) If E[µ(I)X − Y ] > 0 (either finite or infinite) then, by the SLLN, S n (I) → +∞ a.s., hence the same happens to the process {Z 2n (I)} n∈N since Z 2n (I) ≥ S n (I) for all n ∈ N (see equation (4.9) and the remark afterwards). This implies that inf n≥0 S n (I) =: S −∞ > −∞ a.s. and the Markov chain {Z 2n (I)} n∈N is transient. As a consequence P(sup n∈N τ n < +∞) = 1. (2) When the distribution of X n − Y n is not degenerate (that is, it is not δ 0 ) then according to [4,Theorem 4  (1) For every interval I ⊆ [0, 1], When 0 ≥ ∆ > −∞ then for every ε > 0 there exists n 0 such that for every n ≥ n 0 we have |S n (I)/n − ∆| < ε/2. Consider a weak descending ladder time n 1 ≥ n 0 ; it is clear that, for every n ≥ n 1 then min k≤n S k (I) = S kn (I) for some k n such that n ≥ k n ≥ n 1 ≥ n 0 .
it suffices to prove that a.s.
We recall that X i is a sum of X i Bernoulli random variables of parameter µ(A) and that the family of these Bernoulli variables is independent of the family {X n } n∈N . Thus and we are therefore left to prove that a.s.
so that after taking expectation (over N n ) and the statement follows from the Borel-Cantelli's Lemma. (3) The a.s. convergence Z n /n → E[X − Y ]/2 as n → +∞ comes from Proposition 4.2.
As for the second part, if A is an interval then the claim follows trivially by applying Proposition 4.2 to Z n (I) and Z n . Let so that equation (4.16) implies equation (2.4). By the above arguments, it suffices to show the following: if P n , P are probability measures on R, so that for every Borel set A, P n (A) → P (A), then sup t |F n (t) − F (t)| =: F n − F ∞ → 0, where F n and F are the corresponding distribution functions. Let {x i } be the set of atoms on P , let {p i } be their masses and let H(t) = ∞ i=1 p i I (−∞,t] (x i ) be the distribution function of the subprobability measure ∞ i=1 p i δ xi . Let, for every i, p n i = P n ({x i }), by assumption p n i → p i . Let H n be the distribution function of the subprobability measure ∞ i=1 p n i δ xi . Since P n ({x 1 , x 2 , . . .}) → P ({x 1 , x 2 , . . .}), from Scheffe's theorem, it follows that H n − H ∞ → 0. Since for every t, F n (t) → F (t), we have that (F n (t) − H n (t)) → (F (t) − H(t)) for every t. The function F (t) − H(t) is a continuous distribution function of a subprobability measure, so the pointwise convergence implies uniform convergence. So, The following Lemma is well known and we include it for the sake of completeness. Lemma 4.3. Consider two N-valued random variables X and Y on N with laws ρ X and ρ Y respectively; let Φ X (z) = E[z X ] = i=0 ρ X (i)z i and Φ Y (z) = E[z Y ] = i=0 ρ Y (i)z i the corresponding generating functions. Let {X i } i∈N be a i.i.d. sequence of random variables with law ρ X . Finally let {Z i } i∈N be a generic sequence of N-valued random variables with laws {ρ Zn } n∈N and generating functions {Φ Zn } n∈N .

12
(3) If Z = Y i=0 X i then the law of Z is ρ Z = i∈N ρ Y (i)( * i ρ X ) (where * i ρ X is the convolution of i copies of ρ X ) and Φ Z = Φ Y • Φ X .

Proof.
(1) It is straightforward. (2) The explicit expression of the law is trivial and where the last equality comes from the independence. The case of the process {Z n ([0, f ))} n∈N is completely analogous. In this case we just need to note that F (f − c )Φ ′ (1) ≤ 1.