A population model with non-neutral mutations using branching processes with immigration

We consider a stationary continuous model of random size population with non-neutral mutations using a continuous state branching process with non-homogeneous immigration. We assume the type (or mutation) of the immigrants is random given by a constant mutation rate measure. We determine some genealogical properties of this process such as: distribution of the time to the most recent common ancestor (MRCA), bottleneck effect at the time to the MRCA (which might be drastic for some mutation rate measures), favorable type for the MRCA, asymptotics of the number of ancestors.


1.1.
Motivations. Galton-Watson (GW) processes are branching processes modeling discrete populations in discrete time. Since (non-degenerate) GW processes either become extinct or blow up at infinity, one needs to consider a stationary version of the GW process, such as a sub-critical GW process with immigration to model the stationary population. It is well known that the rescaled limit in time and space of GW processes are continuous-state branching (CB) processes, see Lamperti [27], and that the rescaled limit of GW processes with immigration are CB processes with immigration (CBI), see Kawazu and Watanabe [26]. Sub-critical CB process becomes extinct a.s., and conditioning this process not to become extinct gives a CBI with a particular immigration which is a natural continuous model for populations with stationary random size. For the study of genealogical properties of CBI we are interested in, see Chen and Delmas [14] and Bi [9], for a more general immigration. The aim of this paper is to consider the simplest CB process (with quadratic branching mechanism) and an immigration taking into account non-neutral mutations. We shall prove the existence of this stationary continuous process and give some genealogical properties, such as a bottleneck effect at the time to the most recent common ancestor (TMRCA), favorable mutation for the most recent common ancestor (MRCA) of the process. as well as Nordborg and Joyce [19] for continuous models in which the type of the mutant does not depend on the type of the parent. In Fearnhead [21] and Taylor [33] the MRCA is studied. In particular it is shown in [21] that the expected fitness of MRCA is greater than that of a randomly chosen individual.
Notice that the non-neutral models studied by Donnelly and Kurtz [17] are nonstationary. Such models could be made stationary by conditioning on the non-extinction of all the types, see Foucart and Hénard [25] for a work in this direction. For nonstationary models see also Bianconi, Ferretti and Franz [10] for constant size discrete population in discrete time. In those latter models the non-neutral mutations are described using an immigration at constant rate but with various fitness.
1.3. Random size population models. Another large literature is devoted to random size population using branching processes. Neutral models are now well known, see Bertoin [7,8] for GW processes, Champagnat, Lambert and Richard [12] for Crump-Mode-Jagers branching processes and Abraham and Delmas [1,2] with a CB process presentation in [1] and a genealogical tree approach using continuous random tree presentation in [2]. Conditioning on non-extinction will provide a stationary model, see [14] in the CB process setting.
One can also use multi-type processes for non-neutral mutation models; in all those models the rate of mutation is proportional to the size of the current population. For the discrete setting, we refer to Athreay and Ney [5] on multi-type GW processes, see also Buiculescu [11] or Nakagawa [30] for sub-critical multi-type GW processes conditioned on non-extinction. Those latter processes provide natural stationary models. Similar results exist for multi-type CB process or non-homogeneous super-processes (which correspond to infinitely many types) conditioned on non-extinction see Champagnat and Roelly [13] for the former case and Delmas and Hénard [16] for the more general latter case. A non-stationary model with non-neutral mutation is also given in Abraham, Delmas and Hoscheit [3], with a model of immigration with (random) increasing fitness. In this model, the whole population is again a CB process.
1.4. The model. Our present model of a CB process with non-homogeneous immigration follows the approach of [22] and [10] as we consider non-neutral mutation provided by an immigration at constant rate.
For simplicity, we shall restrict our-self to the quadratic branching mechanism ψ(λ) = βλ 2 , with β > 0. For θ > 0, let Y θ = (Y θ t , t ≥ 0) denote a CB process with branching mechanism ψ θ (λ) = ψ(λ + θ) − ψ(θ). It is well known that Y θ is stochastically larger than Y q for θ ≤ q. In particular, we shall say that the type (or mutation) θ is more advantageous than the type q. We shall consider a stationary CBI (Z t , t ∈ R) with non-homogeneous immigration such that at rate 2β dt there is an immigration of a CB process starting with an infinitesimal mass and of type θ, with θ chosen according to a σ-finite measure µ(dθ). We call µ the mutation rate measure.
In [14], the immigration was homogeneous and the mutation rate measure was a Dirac mass measure; we shall call this model CBI with neutral mutations. In the framework of [14] the process Z is distributed as the initial CB process conditioned on non-extinction (or Qprocess) under its stationary measure. The immigration can also be seen as the descendants of an immortal individual. This description of the genealogy using an immortal individual is in the same spirit as the bottom individual in the modified look-down process in [18]. Even if this interpretation is no more valid in our setting, we might keep the corresponding vocabulary as MRCA or TMRCA.
The mutation rate measure allows to consider non-neutral mutations. In our model different CB processes with different branching mechanism coexist at the same time and all mutations eventually die out. Contrary to [3], who considered only advantageous mutations (that is advantageous immigration with rate proportional to the size of the population), the type of the immigrants in our model is random and there is no improvement of the type as time goes on. One of the advantages of our model is that it has a stationary version, which we shall consider. Notice also that the size of the population is random (and different from an homogeneous CBI unless the mutation rate measure is a constant time a Dirac mass).
1.5. Presentation of the results. After some preliminaries on CB processes in Section 2, we define precisely our model in Section 3. In particular, we give an integrability condition on the mutation rate measure µ for the process Z to be well defined (Theorem 3.1) and we check that Z is continuous (Theorem 3.3). We give the expectation of Z t (Corollary 3.4) which might be infinite and characterize the mutation rate measure for which the population size is always strictly positive (Proposition 3.5). Notice only this case is biologically meaningful. We also give (Lemma 3.2) the distribution of time to the most recent common ancestor (TMRCA) which is seen as the first immigration time of an ancestor of the current population living at time t = 0.
We study the type Θ of the MRCA in fact that the type of the first immigrant having descendants at time t = 0 in Section 4. In particular, we get that if µ is a probability measure, then Θ is stochastically more favorable than the type of a random immigrants given by µ (that is Θ is stochastically less than Θ ′ with probability measure µ).
Using arguments close to [14] we give in Section 5 the distribution of the size Z A of the population at the TMRCA (Proposition 5.3) and check that the size of the population at the TMRCA is stochastically smaller than the size of the population at fixed time (which is the stationary measure). This can be interpreted as a bottleneck effect.
In Section 6.1 we give (Lemma 6.1) the asymptotic number of immigrants who still have descendants in the current population at time t = 0. In Section 6.2, we give a precise description of the genealogical structure of the population relying on the tree structure of the Brownian excursion. We study in Section 6.3 the asymptotic number M s of ancestors s unit of time in the past of the current population. In particular we get (Proposition 6.5) that βM s is of order Z 0 /s, which is similar to the CBI case with neutral mutation see [14] or Berestycki, Berestycki and Limic [6] for Λ-coalescent models. We also give the fluctuations (Theorem 6.6) which are similar to the neutral case if the mutation rate behaves nicely. Section 7 is devoted to the stable mutation rate: with c > 0 and α ∈ (0, 1). In this case we have E[Z 0 ] = +∞ and E[Z A ] finite iff α ∈ (1/2, 1). In particular for α ∈ (1/2, 1) we have a drastic bottleneck effect as the ratio E[Z A ]/E[Z 0 ] is equal to 0. We also prove (Proposition 7.2) that the type of the MRCA is (stochastically) more advantageous than the type of an individual taken at random in the current population, see also [21] for similar behavior in a different model. We conjecture this result holds for any mutation rate measure. We get in Section 7.3 that the number of families at s unit of time in the past for the model with non-neutral mutations behave as s −α and that this result can not be compared to the neutral case in [14] even with stable branching mechanism where the number of families at s unit of time in the past is of order 1/| log(s)|. Concerning the fluctuations of the number of ancestors, M s , we get that results given in Theorem 6.6 holds iff α ∈ (0, 1/2) and we get a deterministic limit for α ∈ (1/2, 1). We interpret this latter phenomenon as a law of large number effect from the large number of small populations generated by immigrants with very disadvantageous mutations which is preponderant to the fluctuations of the number of ancestors in each of the immigrant populations.

Preliminaries and notations
We consider a quadratic branching mechanism ψ(λ) = βλ 2 for some fixed β > 0. We will consider a family (ψ θ , θ ≥ 0) of (sub)-critical branching mechanism defined by: For every fixed θ ≥ 0, let P ψ θ x be the law of a CB process, Y θ = (Y θ t , t ≥ 0), started at mass x with branching mechanism ψ θ . Let E ψ θ x be the corresponding expectation and N ψ θ be the canonical measure (excursion measure) associated to Y θ . In particular N ψ θ is a σ-finite measure on the set D 0 of continuous functions from (0, ∞) to [0, ∞) having zero as a trap (for a function f , this means f (s) = 0 implies f (t) = 0 for all t ≥ s). According to [1], see also Abraham, Delmas and Voisin [4], it is possible to define the processes (Y θ , θ ≥ 0) on the same space so that a.s.
When there is no confusion we shall write Y for Y θ that is for example We recall some well known results on quadratic CB processes. For every t ≥ 0 and λ > −2θ/(1 − e −θt ), we have: Notice that u θ satisfies the backward and forward equations for λ ≥ 0 and t ≥ 0: with initial conditions u θ (λ, 0) = λ and u θ (0, t) = 0. It is easy to deduce that for t ≥ 0: where we set for t ≥ 0: It is easy to get that for t ≥ 0 and λ ≥ 0: be the lifetime of Y and set c θ (t) = N ψ θ [ζ > t]. Then we have: Notice that for t, s > 0, we have c θ (t + s) = u θ (c θ (t), s). We also have for t > 0: 3. Definition and properties of the total size process For a Borel measure µ and a measurable non-negative function f defined on the same space, we will write µ, f = f (x) µ(dx).
Let µ be a non-zero Borel σ-finite measure on (0, +∞), which we shall call a mutation rate measure. Consider under P a Poisson point measure . Let E be the expectation corresponding to the probability measure P. For i ∈ I, we shall call Y i a family, θ i its type (or mutation) and t i its birth time. Define the super-process Z = (Z t , t ∈ R) by: with the convention that Y i t = 0 for t < 0 and δ θ denotes the Dirac mass at θ. By construction Z is a stationary Markovian σ-finite measure-valued process. We shall consider the corresponding total size process Z = (Z t , t ∈ R) defined by: Notice that Z is stationary but it is not Markovian unless µ is a constant times a Dirac mass. The process Z is a CB process with a non-homogeneous immigration. It will represent the evolution of a random size population with non-neutral mutations in a stationary regime. The genealogy of Z will be defined in Section 6.2. First we will consider the condition on µ such that Z is well defined.
The random variable Z t is finite a.s. if and only if the following conditions are satisfied: The distribution of Z t is characterized by its Laplace transform, for λ ≥ 0: Proof. By the exponential formula, one obtains that, for F non-negative measurable, (4): Letting λ → 0 entails that The right hand side of (10) is equivalent to the existence of some λ > 0, such that: As log(1 + λ/2θ) is equivalent to | log θ| (resp. λ/2θ) as θ goes to 0+ (resp. +∞), we deduce that (11) holds if and only if (8) holds.
Before giving other properties of the process Z, we shall study the time A to the first immigration time of an ancestor (or equivalently the TMRCA) of the current population living at time 0 which is defined as: We have for all t ≥ 0: Under conditions (8), we get that A is a.s. finite.
Proof. The property of the Poisson random measure implies that for t ≥ 0: where we used (7) for the last equality. Under conditions (8), we get that is finite for any t > 0, which implies, thanks to dominated convergence, that lim t→+∞ P(A < t) = 1 that is A is a.s. finite. Theorem 3.3. Under conditions (8), the process Z is continuous.
In particular, we deduce that under conditions (8), Z is a stationary Markov process with values in the set of finite measures on R + .

Proof.
To prove the continuity of the process Z, we notice that by stationarity, we just need to prove the continuity of Z on [0, 1]. Let c > 0 be a finite constant and consider the truncated Notice that Z c and (Z t , t ∈ [0, 1]) coincide on {A ≤ c}. Since A is a.s. finite, to get the continuity of Z on [0, 1], we just need to prove that Z c is continuous. We shall check the Kolmogorov criterion for Z c . Let λ ≥ 0 and γ ≥ 0, 0 ≤ s ≤ t ≤ 1. We have: where we used (4) for the fourth equality, and the equality (3)) for the fifth. Notice that for fixed r > 0, there exists a constant C r > 0 such that for all θ > 0: Therefore, there exists a constant c 1 ≥ 1 such that for λ, γ ∈ R: We deduce that under conditions (8), the function is analytic in (λ, γ) in a neighborhood of 0 for example on {(λ, γ); |λ| + |γ| ≤ 1/4c 1 }. Taking γ = −λ, this implies that for |λ| ≤ 1/8c 1 , we have: Using (13), an easy computation yields that there exists a constant c 2 such that: Then using that ∆ θ t−s ≤ β(t − s), we get there exists a constant c 3 such that: This gives the Kolmogorov criterion for Z c . Thus Z c is continuous, which ends the proof.
We give the first moment of Z.
Proof. Using (9), we get: We give a criterion for Z to reach 0. See also Foucart and Bravo [24] for such a criterion for CBI.

Type of the MRCA
We assume that conditions (8) hold.
Because of the stationarity, we shall focus on the MRCA of the current population living at time 0. Recall the TMRCA is given by (12). We set i 0 ∈ I the (unique) index i such that A = −t i . We shall say that Y i 0 is the oldest family. We define the type of the MRCA that is of the oldest immigrant family as: Θ = θ i 0 . We give the joint distribution of the TMRCA and the type of the MRCA.
Lemma 4.1. We have for every t ∈ R, θ > 0, Proof. For f non-negative measurable, we get: We deduce that P(A ∈ dt, Θ ∈ dθ) = 2βc θ (t) P(A < t) dtµ(dθ). Then, using Lemma 3.2, it follows that: Using Lemma 3.2, we can derive the distribution µ MRCA t of the type of MRCA given the TMRCA being equal to t: Notice that the function θ → θ(e 2βθt − 1) −1 decreases to 0 as θ increases to +∞. Intuitively, the distribution of the type of the MRCA is more likely to focus on the favorable θ (that is θ small which corresponds to large population) than that on the θ large (which corresponds to small population). In particular if µ is a probability measure then µ MRCA t is stochastically smaller that µ. This means the type of the MRCA (given {A = t}) is stochastically less, which means stochastically more favorable, than the type of a random immigrant. Notice also that µ MRCA t is stochastically decreasing with t which means that the oldest family, Y i 0 , is stochastically increasing with |t i 0 |.

Bottleneck effect
We assume that conditions (8) hold. We consider Z −A , which we shall denote Z A , the size of the population at the TMRCA:  [14], it is easy to get the following result.
Lemma 5.1. The joint distribution of (Z A , Z O , Z I , A, Θ) is characterized by: for λ, γ, η ∈ [0, +∞), and t, θ ∈ (0, +∞), We deduce the following result. Now we concentrate on the population size at the MRCA. Recall ∆ θ t defined in (3). Proposition 5.3. Let t ∈ (0, +∞). We have for η ≥ 0: Furthermore conditionally on A (or not), Z A is stochastically smaller than Z 0 , that is for all z > 0: The fact that Z A is stochastically smaller than Z 0 corresponds to the bottleneck effect.
Proof. Using (16), we get: Then, using (4), (7) and (3), it is easy to get (17). This readily implies (18). We now prove the stochastic order. First notice that ∆ θ t ≤ 1/(2θ). We deduce that for all η ≥ 0, we have: ]. This means that Z A is smaller than Z 0 in the Laplace transform order. We will however prove the stronger result on the stochastic order.
We deduce from (17) that conditionally on . As recalled in Section 2, it is possible to define on the same space two CB processes Y 1 and Y 2 such that Y 1 ≤ Y 2 a.e. and Y 1 (resp. Y 2 ) is distributed under N ψ 1/(2∆ θ t ) [dY ] (resp. N ψ θ [dY ]) since θ ≤ 1/(2∆ θ t ). We deduce that Z 1 can be defined on a possible enlarged space so that there exists a PPM and such that a.s. for all i ∈ I, we deduce that Z A (conditionally on {A = t}) is stochastically less than Z 0 . This gives the first part of (19). By integrating the first part of (19) with respect to P(A ∈ dt), we deduce the second part of (19).
As a direct consequence of (18) we can compute the expectation of Z A , which will be used in Section 7.

Asymptotics for the number of families and ancestors
We assume that conditions (8) hold.
We set: We have the following result.
where we used (6) and (7) for the first equality. As s → 0+, Λ(s) goes to infinity. Then notice that (N Λ −1 (s) , s ≥ 0) is a Poisson process with parameter 1. We deduce the result from the strong law of large numbers for Lévy processes.
6.2. The genealogical tree. In order to consider the number of ancestors M s at time −s of the current population living at time 0, we need to introduce the genealogical tree for a CB process, see Le Gall [28] or Duquesne and Le Gall [20]. Since the branching mechanism is quadratic, we will code the genealogical tree using Brownian excursion. Let W = (W t , t ∈ R + ) be a Brownian motion. We consider the Brownian motion W θ = (W θ t , t ∈ R + ) with negative drift and the corresponding reflected process above its minimum H θ = (H θ (t), t ∈ R + ): We deduce from equation (1.7) in [20] that H θ is the height process associated to the branching mechanism ψ θ . For a function H, we set: max(H) = max(H(t), t ∈ R + ).
Let N ψ θ [dH θ ] be the excursion measure of H θ above 0 normalized such that N ψ θ [max(H θ ) ≥ r] = c θ (r). Let (L x t (H θ ), t ∈ R + , x ∈ R + ) be the local time of H θ at time t and level x. Let ζ = inf{t > 0; H θ (t) = 0} be the duration of the excursion H θ under N ψ θ [dH θ ]. We recall that (L r ζ (H θ ), r ∈ R + ) under N ψ θ is distributed as Y θ under N ψ θ . From now on we shall identify Y θ with (L r ζ (H θ ), r ∈ R + ) and write N ψ θ for N ψ θ . When there is no confusion, we shall write H for H θ and Y for Y θ . We now recall the construction of the genealogical tree of the CB process Y from H.
Let f be a continuous non-negative function defined on [0, +∞), such that f (0) = 0, with compact support. We set ζ f = sup{t; f (t) > 0}, with the convention sup ∅ = 0. Let d f be the non-negative function defined by: To emphasize the dependence of Y and R in H, we may write Y a (H) and R a,b (H). We give the joint distribution of (Y a , Y b , R a,b ).
Lemma 6.2. Let 0 < a < b. For λ, ρ, η ≥ 0, we have: r). Proof. We have: (λ, η), a), where we used the property of the PPM k∈Ka δ H k for the second equality, and as Y s = 0 on {ζ < s}, for the third equality.
In order to describe the genealogical structure of Z, following the beginning of Section 3, we shall consider under P the PPM: Let s > 0. We define the super-process for the number of ancestors for the population at time 0 of each family: Then the number of ancestors at time −s of the population at time 0, is: We will first consider the joint distribution of (Z 0 , Z −s , M s ). Proposition 6.3. Let ρ, λ and η be non-negative measurable functions defined on R + . We have: with for r > 0: In particular, we deduce that a.s.: Proof. We have: Using Lemma 6.2, we get: Then use (4) to end the proof.
In fact, we shall see that the convergence in the previous Remark is an a.s. convergence.
Proposition 6.5. We have a.s.: Notice the order of M s do not depend on µ: the non-neutral mutations do not change the asymptotics of the number of ancestors.
Proof. Let λ = 0 and η(θ) = η 0 . We deduce from (23) that: Therefore, conditionally on Z −s , the number of ancestors M s is a Poisson random variable with mean W s = t i ≤−s c θ i (s)Y i −s−t i . We first prove that a.s.: Notice that the function x → x/(e x −1) is decreasing on (0, +∞) and is less than 1. We deduce that βsc θ (s) ≤ 1 and therefore βsW s ≤ Z −s .
Recall that M s is increasing, is conditionally on Z −s a Poisson random variable with mean W s . Then use (24) and properties of Poisson distributions to get that a.s.: lim s→0+ βsM s /βsW s = 1.

This and (24) end the proof.
We have the following partial result on the fluctuations. Theorem 6.6. We have the convergence in distribution of ((Z 0 , with G a standard Gaussian random variable independent of Z 0 . Under the conditions: we have the convergence in distribution of ((Z 0 , , √ Z 0 G) with G and G ′ two independent standard Gaussian random variables independent of Z 0 .
If conditions (25) do not hold, we may have a very different behavior, see the stable case in Section 7.
We have the following representation of the limit in Theorem 6.6.

Stable mutation rate measure
Let c > 0 and α ∈ (0, 1). In this Section, we will consider the stable mutation rate: Notice that µ satisfies conditions (8). In addition, notice that E[Z 0 ] = +∞ and that µ, 1 = +∞ which in turn implies that {t; Z t = 0} = ∅ a.s. thanks to Proposition 3.5. 7.1. Bottleneck effect. We present a drastic bottleneck effect which was not observed in [14].

7.2.
Type of the MRCA and type of a random individual. Recall from Section 2 that a CB process Y θ has type θ and that Y θ is stochastically larger than Y q if θ ≤ q. We shall say that the type (or mutation) θ is more advantageous than the type q. The following proposition states that the type of the MRCA is (stochastically) more advantageous than the type of an individual taken at random at the current time.
Recall from Section 3 that the type of an individual in family Y i is θ i . We define the type Θ * of an individual taken at random at time 0 as follows: conditionally on Z, Θ * is equal to θ i with probability Y i −t i /Z 0 . Proposition 7.2. We have that Θ is stochastically smaller than Θ * : for all q ≥ 0, P(Θ * ≤ q) ≤ P(Θ ≤ q).
Proof. Firstly, we give the distribution of Θ. We have for θ > 0: where we used (16) for the first equality, the change of variables r = 2βqt (t fixed) and s = 2βθt as well as for the second equality. Set Q = 2c Θ α so that for q > 0: Then we deduce that Q is distributed as ES α /a 1 , where E is an exponential random variable with mean 1 independent of the random variable S whose density is: Secondly, we give the distribution of Θ * . Let F be a non-negative measurable function, we have: where we used the definition of Θ * for the first equality, the PPM properties for the third equality, (2) and (9) for the fourth equality. Set and use the change of variables q = λr (λ fixed), u = 1 − e −2βθt (θ fixed) and s = aθ/λ (θ fixed) with a α = a 1 /a 2 and b = 2/a to get: Set Q * = 2c Θ α * and we deduce that: Then similarly Q * is distributed as ES α * /a 1 with S * a random variable independent of E and with density: . We deduce that P(S * ≤ s) ≤ P(S ≤ s) for all s ≥ 0, that is S * is stochastically larger than S. This implies that Q * (resp. Θ * ) is stochastically larger than Q (resp. Θ).

Number of families.
We compare the number of families with the neutral case (quadratic and stable branching mechanism). Recall Λ(s) defined in (22). We deduce from (30) that Λ(s) = C 3 s −α . Lemma 6.1 implies that a.s.: We can compare (31) with the stationary stable case with immigration, see [14], that is ψ(λ) = λ a + bλ with 1 < a ≤ 2 and b > 0. According to Section 6 in [14], we have that the number N * s of families alive at time −s which are still alive at time 0 is a.s. equivalent, as s goes down to 0, to Λ * defined by (31) in [14]. Using that: Notice that Λ * (s) is equivalent to (a/a − 1)| log(s)| as s goes down to 0. This implies that a.s. lim s→0+ | log(s)| −1 N * s = a − 1 a · Therefore the number of families for stable case with neutral mutation is much smaller than that of CB process with non-neutral mutations (that is with stable rate of mutation).
We explain the contribution of −h(α) as a law of large number effect produced by the large number of (small) populations with large parameter θ. This effect is negligible if conditions (25) hold that is α ∈ (0, 1/2) but is significant for α ∈ [1/2, 1).