Random evolutionary dynamics driven by fitness and house-of-cards mutations. Sampling formulae
Thierry Huillet (LPTM)
TL;DR
This paper analyzes the impact of house-of-cards mutations on evolutionary dynamics, deriving sampling formulas and Ewens sampling distributions for multi-allelic Wright-Fisher models under weak selection and mutation regimes.
Contribution
It introduces new sampling formulae for evolutionary models with house-of-cards mutations and extends Ewens sampling formulas to these complex scenarios.
Findings
Derived normalizing partition functions for invariant distributions.
Established generalized Ewens sampling formulas for house-of-cards mutation models.
Addressed sampling issues in infinite-alleles weak limits.
Abstract
We first revisit the multi-allelic mutation-fitness balance problem, especially when mutations obey a house of cards condition, where the discrete-time deterministic evolutionary dynamics of the allelic frequencies derives from a Shahshahani potential. We then consider multi-allelic Wright-Fisher stochastic models whose deviation to neutrality is from the Shahsha-hani mutation/selection potential. We next focus on the weak selection, weak mutation cases and, making use of a Gamma calculus, we compute the normalizing partition functions of the invariant probability densities appearing in their Wright-Fisher diffusive approximations. Using these results, Generalized Ewens sampling formulae (ESF) from the equilibrium distributions are derived. We start treating the ESF in the mixed mutation/selection potential case and then we restrict ourselves to the ESF in the simpler house-of-cards…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Random evolutionary dynamics driven by fitness and house-of-cards mutations.
Sampling formulae
Thierry E. Huillet
Laboratoire de Physique Théorique et Modélisation
CNRS-UMR 8089 et Université de Cergy-Pontoise, 2 Avenue Adolphe Chauvin, 95302, Cergy-Pontoise, FRANCE
E-mail: [email protected]
Abstract.
We first revisit the multi-allelic mutation-fitness balance problem, especially when mutations obey a house of cards condition, where the discrete-time deterministic evolutionary dynamics of the allelic frequencies derives from a Shahshahani potential. We then consider multi-allelic Wright-Fisher stochastic models whose deviation to neutrality is from the Shahshahani mutation/selection potential. We next focus on the weak selection, weak mutation cases and, making use of a Gamma calculus, we compute the normalizing partition functions of the invariant probability densities appearing in their Wright-Fisher diffusive approximations. Using these results, Generalized Ewens sampling formulae (ESF) from the equilibrium distributions are derived. We start treating the ESF in the mixed mutation/selection potential case and then we restrict ourselves to the ESF in the simpler house-of-cards mutations only situation. We also address some issues concerning sampling problems from infinitely-many alleles weak limits.
Keywords: Evolutionary genetics, fitness landscape, house-of-cards mutations, Shahshahani gradient, Wright-Fisher random genetic drift, Gamma calculus, generalized Ewens sampling formulae.
1. Introduction
Population genetics is concerned with the fate of multi-allelic population frequencies when various driving ‘forces’ such as selection or mutation are involved. We will briefly revisit the basics of the deterministic dynamics arising in discrete-time asexual evolutionary genetics when the origin of motion is either fitness or mutations or both. We will mostly consider the multi-allelic diploid case under the Hardy-Weinberg hypothesis. Firstly, we will consider evolution under general fitness mechanisms (Section ), then we deal with general mutation mechanisms (Section ). Some particular fitness/mutation patterns are discussed in the process. In some cases, the dynamics driven by fitness only takes the form of a Shahshahani gradient dynamics, as deriving from a Shahshahani selection potential, [30], [5]. In Section , we give one way to combine the fitness and the mutation effects. All these issues are of course part of the standard models discussed for example in [8], [2], [21], [22], [26] and [20]. We then focus on a reversible mutation pattern in which mutation probabilities between any two states only depend on the target state (the house of cards condition for mutations). When combined with selection effects, the dynamics driven both by selection and house-of-cards mutations takes the form of a Shahshahani gradient-like dynamics with drift deriving from a Shahshahani potential mixing additively the mutation and selection potentials, [15].
On top of such deterministic dynamics, we add (Section ) a random genetic drift while considering a multi-allelic Wright-Fisher Markov chain whose deviation to neutrality appears as a drift deriving from the Shahshahani mutation/selection potential. Some scalings brings this Markov chain into multi-allelic Wright-Fisher diffusion processes with a unique invariant probability density describing the joint allelic frequencies at equilibrium, [32]. We consider the weak selection, weak mutation cases and we compute the normalizing partition functions of the corresponding invariant densities.
Generalized Ewens sampling formulae (ESF) from such equilibrium allelic frequencies are then desirable and these are obtained in Section ; they make use of the explicit expressions of the partition functions just introduced. Such ancient and long-standing questions go back to [37]; see [34] and [25]. We start treating the ESF in the mixed mutation/selection potential case (Section ) and then we restrict ourselves to the ESF in the simpler house-of-cards mutations only context (Section ). We show that such problems are amenable to the evaluation of some functionals of skew-symmetric Dirichlet distributed random variables (rvs) and for this purpose, we make extensive use of a ‘Gamma-calculus’ precisely designed to evaluate such functionals. We also address some issues concerning the availability of sampling formulas in some infinitely-many alleles weak limits.
2. Evolution driven by fitness: the deterministic point of view
2.1. Single locus: the case of a diploid population with alleles
We briefly describe the frequency distribution dynamics when fitness (or selection) only drives the process. We consider diploid populations.
2.1.1. Joint evolutionary dynamics
Consider alleles , and let be the genotypes attached to a single locus. Let , stand for the absolute fitness of the genotypes . We shall assume ( being the probability of an surviving to maturity, it is natural to take ). Let then be the symmetric fitness matrix with entry .
Assume the current frequency distribution at generation of the genotypes is given by Let be the frequencies array with entry , obeying . The joint evolutionary dynamics in the diploid case is given by the updating111The symbol + is a common and useful notation to denote the updated frequency from generation to .:
[TABLE]
where the relative fitness of the genotype is . The joint dynamics takes on the matrix form:
[TABLE]
where stands for the (commutative) Hadamard product of matrices.
With a column-vector of ones, its transpose and the matrix whose entries are all (the identity for )222In the sequel, a boldface variable, say , will represent a column-vector so that its transpose, say , will be a row-vector. Similarly, will stand for the transpose of some matrix ., then
[TABLE]
Let
[TABLE]
stand respectively for the genotypic variance in absolute fitness and the diploid variance in relative fitness. Note that, owing to and ,
[TABLE]
The partial increase of the mean fitness (where the change of mean fitness only comes through changes in the frequencies ) is given by
[TABLE]
with a relative rate of increase: . The latter result (3) constitutes the diploid version of the Fisher fundamental theorem of natural selection for asexual populations; see ([8], Sections , [28]) and [27]) for a deeper insight on its meaning.
2.1.2. Marginal multi-allelic dynamics
Assuming a Hardy-Weinberg equilibrium, the frequency distribution at generation say , of the genotypes is given by: where is the marginal frequency of type- allele in a genotypic population. The marginal allelic frequency vector is ( is the unit -column-vector) and the mean allelic fitness is given by the quadratic form: . The mean fitness is homogeneous of degree in the variables . Let
[TABLE]
be respectively the genotypic variance in absolute fitness and the diploid variance in relative fitness.
If we first define the frequency-dependent marginal fitness of by , the marginal dynamics is given by:
[TABLE]
Letting diag if the allelic frequency distribution is summarized in the column-vector , , (5) reads in vector form333 is the Schur product of vector and vector See [21] page 238 for a similar notational convenience.
[TABLE]
where maps the dimensional simplex
[TABLE]
into itself.
**Two particular cases:
**
The haploid case. Let denote the absolute fitnesses of the alleles and suppose, without loss of generality, that (so that allele has largest fitness ). Let denote the mean fitness of the population. Plugging in (6), the dynamics (6) boils down to
[TABLE]
giving the update of the frequency distribution of alleles in an haploid population where alleles (and not pairs of alleles) are attached to some locus. Along (7), the absolute mean fitness increases. Indeed, with :
[TABLE]
and it is except when is such that , the set of alleles with maximal fitness. Such s are equilibrium states of (7). In particular, an extremal vector with is a pure (or monomorphic) equilibrium state.
The mean fitness is maximal at equilibrium. The relative rate of increase of is:
[TABLE]
where is the variance in relative fitness given by
[TABLE]
Thus the population mean fitness is non-decreasing. As a consequence, if there is an unique allele whose fitness strictly dominates the ones of the others, starting from any initial state which is not an extremal point of , the haploid trajectories will ultimately converge to this fittest state (survival of the fittest allele).
The diploid case with multiplicative fitnesses. Suppose that , or in matrix form that Then selection acts on the gametes rather than on the genotypes. Observing , the dynamics (6) boils down to the one (7) of haploid populations. However, the mean fitness in this case is and not as in the haploid case.
2.1.3. Increase of mean fitness for diploid populations
Similarly, defining , in the diploid case we have
[TABLE]
The mean fitness for diploid populations, as a Lyapunov function, increases as time passes by, vanishing only when the process has reached equilibrium. Equation (10) constitutes the mean fitness increase theorem.
The partial rate of increase of the mean fitness due to frequency shifts only (see [8]) is It satisfies
[TABLE]
where is the allelic variance in relative fitness
[TABLE]
Equation (11) constitutes the diploid version of the Fisher fundamental theorem of natural selection under the Hardy-Weinberg condition involving random mating. See [28], [8] and [3].
2.1.4. **Shahshahani gradient-like representation of the
allelic dynamics **(6)
There is an alternative vectorial representation of the dynamics (6) emphasizing its gradient-like character. With , define the matrix . It is symmetric, positive semi-definite whose quadratic form vanishes only for the constants and** **for all . is partially invertible on the space orthogonal to the constants with left-inverse
[TABLE]
so with** ** for all , obeying** .** Note and for all** **. Looking for a left-inverse in the weaker sense of the quadratic form, that is satisfying
[TABLE]
for all ** with every would do for any real number . In particular leading to **.
Introduce the quantity . Then, (5) may be recast as the gradient-like dynamics:
[TABLE]
with as a result of observing** **. Note
[TABLE]
Based on [30], [31], the dynamics (13) is of gradient-type with respect to the Shahshahani-Svirezhev metric . Its piecewise trajectories are perpendicular to the level surfaces of with respect to the scalar product given by
[TABLE]
We also have
[TABLE]
From (11) and (12), , which is the length of , is also the square-root of half the allelic variance (the standard deviation) in relative fitness.
2.2. Frequency-dependent fitness
Consider the dynamics (6) in vector form. Let the frequency-dependent marginal fitness of be defined by , , not necessarily of the linear form for some fitness matrix With a new column vector of frequency-dependent marginal fitnesses, we can first think of defining the dynamics on by
[TABLE]
where is the new mean fitness of the allelic population. Unless as before, such dynamics cannot be of Shahshahani-gradient type. This suggests to consider the alternative gradient-like dynamics to (14), still on the simplex:
[TABLE]
where and A particular case is where is now frequency-dependent and symmetric for each and . Fitness landscapes can be more general than quadratic forms.
Examples:
-
Suppose where for all Then the mean fitness is homogeneous of degree in the variables .
-
With , let and suppose is diagonal, so that . Such selection models were considered in [34], [10], [13], [17] and [4], in relation to heterozygozity.
3. Evolution driven by mutation
We now briefly describe the frequency distribution dynamics when mutation is the only driving source of motion.
Assume alleles mutate according to the scheme: with probability satisfying and for all Let be the mutation pattern matrix; we shall assume that the non-negative matrix is irreducible.
3.1. Frequency dynamics under mutation only
Considering first an updating mechanism of the frequencies where only mutations operate, we get
[TABLE]
whose meaning is the one of a “master equation”. In vector form, with the transpose of
[TABLE]
and the update of the frequencies with mutations is given by the linear transformation
[TABLE]
The vector is called the mutation load. We have and if and only if is stochastic: meaning
[TABLE]
Also and then maps into because if , then . The matrix is stochastic and irreducible and so, by Perron-Frobenius theorem, it has a unique strictly positive probability left-eigenvector associated to the real dominant eigenvalue Let be this row-vector, so obeying , or . Under the irreducibility assumption on , the frequencies dynamics involving only mutations has a unique polymorphic equilibrium fixed point . When is primitive then This shows that, at generation ,
[TABLE]
regardless of the initial condition belonging to . The equilibrium vector is asymptotically stable. These considerations are the same as limit (ergodic) theorems for Markov chains.
Note finally that from (17):
[TABLE]
where is the quadratic mutation potential. The probability right-eigenvector of uniquely solves with maximal.
3.2. Special mutation patterns
There are many special mutation scenarii which deserve interest.
Reversible mutations: Let solve Define
[TABLE]
We have , so is the stochastic matrix of the time-reversed mutation process at equilibrium with invariant measure . If , then the mutation pattern is said to be time-reversible (detailed balance holds). In this case
[TABLE]
with . When reversible, the matrix is diagonally similar to the matrix which is symmetric with real eigenvalues, so and have real eigenvalues.
If , then and is symmetric itself and thus doubly stochastic. In that case, (the barycenter of ) and such mutation patterns are reversible. Let us give some Examples:
- A model with symmetric mutations is obtained for instance while assuming multiplicative mutations: In this case, with the column vector of the s,
[TABLE]
This mutation pattern is reversible with .
-
Alternatively, assuming for all leads to which is also symmetric.
-
Alternatively, while considering additive mutations: with ,
[TABLE]
and . This mutation pattern is reversible.
It is not necessary that in order to have doubly stochastic. It suffices to impose . In that case although the overall input-output mutation probabilities attached to any state coincide and the equilibrium state again matches with the barycenter of . But since , such mutation patterns are not reversible.
Assume the mutation probabilities only depend on the initial state, that is: for all . Then
[TABLE]
This mutation model is not reversible because
[TABLE]
If for all , the equilibrium state is
[TABLE]
One-way (irreversible) mutations: assume This model is clearly not reversible and, when the associated mutation matrix is irreducible, it has a non-trivial This model includes the cyclic mutation pattern for which , and with
[TABLE]
as in the previous example .
(Kingman house-of-cards mutations, [22]).
Assume the mutation probabilities now only depend on the terminal state, that is: for all , still with Throughout, we assume . Let . Then, , where , and
[TABLE]
The equilibrium state is and it is stable. Note that . This mutation model is reversible because
[TABLE]
In this model the coordinates are decoupled: , depends only on . We shall come back at length to this mutation pattern in the sequel.
- Note that if , namely if , then , (else ) and
[TABLE]
-
In the latter case, in particular. If and , in view of , goes fast to .
-
Another very special case is when . Here and . Starting from any initial condition, the dynamics moves in one-step to (inside the simplex ) and stays there for ever.
4. Evolution driven by combined fitness and mutation forces
Let us now consider the dynamics driven both by fitness and mutation.
4.1. The combined fitness/mutation frequency dynamics; [15], [26], [14], [20]
Combining the fitness and mutation effects consists in applying first the fitness mapping and then let mutation act on the result. Proceeding in this way, we get the ‘fitness-first’ dynamics [15]:
[TABLE]
defining a new nonlinear transformation. Alternatively, where is the new mapping from to to consider. Component-wise, this is also as required
[TABLE]
We have: and so the extremal states are not invariant under and from the fixed-point theorem, there exists some equilibrium state in . Using the representation (13), (17) and (18):
[TABLE]
This is not a Shahshahani gradient-like dynamics in general. [15] also considers a continuous-time version of (22).
When (no mutation) (22) boils down into (13) and when (no selection), (22) boils down into (18).
When both (no mutation) and (no selection), with corresponding neutral
4.2. Fitness/mutation frequency dynamics in the house of cards
condition for mutations
This is a remarkable case where the allelic dynamics driven both by fitness and mutation forces has a Shahshahani gradient-like structure, [15]. Indeed, from (19), and (20) boils down to
[TABLE]
We now have
[TABLE]
where . Indeed, and
[TABLE]
Furthermore,
[TABLE]
Thus (22) can alternatively be written as
[TABLE]
which is of Shahshahani gradient-type with combined mutation/selection potential
[TABLE]
Note that is homogeneous of degree .
4.3. House of cards and polymorphism
If some polymorphic state exists in the interior of the simplex, then where solves
[TABLE]
for some arbitrary Lagrangian parameter ( is an extremum of under the constraint fixed).
-
If (no mutation), the searched is the one solving, as required, up to a multiplicative constant, else if
-
If (no selection), and the searched is the one solving, as required, up to a multiplicative constant so that if .
So where is some fixed point of the map
[TABLE]
If there is a stable fixed point of , such that (the purely selection dynamics has a stable polymorphic equilibrium state in the simplex), then clearly, as long as , , as a fixed point of , exists as well and is unique, as a maximum of the concave potential . If is stable indeed, then so is if , because the new mean fitness function inherits the concavity of if and is concave as well if . If , then . Note however that entails that a stable can be switched to an unstable : strong mutations () can destroy a stable polymorphic state of the selection equation.
5. Evolutionary dynamics: the stochastic point of view
We now consider a stochastic version of evolutionary dynamics biased by selection/mutation effects, thereby adding a “random genetic drift” to the deterministic dynamics.
5.1. Random genetic drift and the Wright-Fisher model
We will consider the multi-allelic Wright-Fisher model with bias, (see [8], section ).
5.1.1. The discrete-time model and first properties
Consider a multi-allelic population with constant size In the haploid (diploid) case, is (twice) the number of effective loci. Let and be two vectors of integers quantifying the size of the allelic populations at two consecutive generations and where is some small parameter fixing the time elapsed between two consecutive generations. With , therefore and both and belong to . Suppose the stochastic evolutionary dynamics is now given by a Markov chain whose one-step probability transition matrix from state to state is given by the multinomial Wright-Fisher model
[TABLE]
Here the s are the coordinates of some mapping translating some bias from neutrality (where is simply the identity). The state-space dimension of this Markov chain is (the number of compositions of integer into non-negative parts). In other words, with ,
[TABLE]
5.1.2. Scalings
If we assume that for some drift function (meaning is the frequency shift per generation of duration ), from the mean and covariance structure of multinomial distributions, with , , we get
[TABLE]
With and to the first-order in ,
[TABLE]
since the deviation of from induce second-order effects in on the variance. This means
[TABLE]
where is an iid dimensional Gaussian sequence with zero mean and covariance matrix For any , using the multinomial structure of the Wright-Fisher transition matrix, we have
[TABLE]
The stochastic dynamics (26) is a discretized version of the diffusion process with continuous sample paths on (making )
[TABLE]
where is a dimensional standard Brownian motion. Such results on the diffusive approximation could be made rigorous using theorems from [6].
5.1.3. Speed (invariant) densities of the scaled process
The speed density of the diffusion process (27), cancelling the probability flux of its Kolmogorov forward (or Fokker-Planck) equation, is (up to a multiplicative constant)
[TABLE]
possibly not normalizable into a probability density function. If, for some potential , then, by applying the gradient theorem to the line integral over some path in ending in
[TABLE]
It may happen that, for some suitable choice of , some normalizing constant turns to a probability density . In such cases, is the limiting equilibrium probability distribution of some random variable describing the asymptotic allelic frequencies. Note that a non-integrable case is the neutral case with no drift (or ), with speed measure on the simplex, but also as in (13), with and non-summable associated speed measure . We shall now consider situations including mutations where is integrable and thus normalizable.
5.1.4. **The mutation/selection potential under the house of
cards condition for mutations444In the sequel, we consider a mean fitness but a frequency-dependent mean fitness would do as well provided is bounded as .**
For instance, if now is the Shahshahani mutation/selection potential under the house of cards condition for mutations (24)
[TABLE]
then, with
[TABLE]
[TABLE]
If , we recognize the skew Dirichlet distribution on the simplex with parameters ,
[TABLE]
with normalizing partition function . We shall call this distribution so that here . Except in the homogeneous case when all for each Dirichlet distributed random variables with distribution (30) on the simplex are not exchangeable and is not invariant under a permutation of its coordinates, [23]. The same holds a fortiori for (29) when .
With as in (28), the deterministic dynamics corresponding to this particular potential is thus , as a version of (23) when the time elapsed between two consecutive generations no longer is but .
5.1.5. Weak mutation probabilities
If (small mutation probabilities, with all now defining mutation rates), with and
[TABLE]
where
[TABLE]
In the asymmetric mutation case when for some common , the probability distribution of on the simplex is not exchangeable either. For similar shapes of the equilibrium distribution in this weak mutation setup, see [37], [24], [25] and [38].
With , , (31) is also
[TABLE]
where
[TABLE]
and is Dirichlet distributed.
5.1.6. Weak mutation and weak selection probabilities
If in addition for some new symmetric fitness differential matrix (involving an order correction to the neutral model of selection), then and
[TABLE]
so that
[TABLE]
for some normalizing constant
[TABLE]
With (the multidimensional beta function) and , we have
[TABLE]
where
[TABLE]
with The latter is the unique invariant probability density of the diffusion process on the simplex, [32], [7]
[TABLE]
where the drift is
[TABLE]
Namely,
[TABLE]
A unit time in the latter diffusion process represents generations of the discrete-time model. See [37], [25], [12] and [1].
5.2. Computing the partition functions and
5.2.1. A constructive formula for computing with Dirichlet distribution
The following ‘Gamma-calculus’ result will be useful (see [18]):
Theorem 1**.**
Consider an asymmetric Dirichlet distributed random variable on the simplex, viz: and let , , with
* Let be any Borel-measurable function for which*
[TABLE]
Then, with , independent random variables defined by , , where gamma, we have
[TABLE]
* If is homogeneous of degree , i.e. if , , , and if then, with ,*
[TABLE]
This shows that computing the expected value of some functional with respect to the asymmetric Dirichlet distribution on the simplex can be achieved while averaging over an identically distributed gamma distributed sample; a much simpler task.
5.2.2. Estimating in (32)
Theorem 2**.**
When ,
[TABLE]
where the summation runs over all integer ordered partitions of as and .
Proof: Consider the statement of Theorem 1 with . This function is homogeneous of degree . Thus
[TABLE]
where the random vector has independent components with gamma Thus, with
[TABLE]
When is large enough, is close to , with
[TABLE]
In the latter expression, the summation runs over all square arrays whose integral non-negative entries sum to (the set )666The number of ordered partitions of the integer into non-negative integral summands, , is . There are thus ways to realize , with non-negative integers., and , integers. Recalling the moment structure of gamma random variables (rvs), , we thus get
[TABLE]
5.2.3. Computing in (34)
We now show that
Theorem 3**.**
[TABLE]
Proof: Consider the statement of Theorem 1. From (38), the right hand-side quantity, namely
[TABLE]
may be interpreted as the Laplace transform in the variable of appearing in the left-hand side. Inverting this Laplace transform and putting yields . This can be used to compute with We have
[TABLE]
where,
[TABLE]
Recalling , we get
[TABLE]
Owing to , this is also
[TABLE]
The right hand-side quantity of statement (38) is thus
[TABLE]
[TABLE]
It is the Laplace-Stieltjes transform of .
But is the Laplace-Stieltjes transform of and the one of . Inverting the Laplace transform and evaluating the result at yields
[TABLE]
Recalling , we thus obtain
[TABLE]
Corollary 4**.**
The alternating series expansion of in (42) is convergent.
Proof: Let
[TABLE]
so that . Let be the largest of the s. Then and
[TABLE]
where is now a vector with nonnegative integral entries summing to . With we now have
[TABLE]
which are the integral moments of order of a Dirichlet random vector (see (58) below). Thus, and is absolutely convergent hence convergent.
Remarks: As gets large, the contribution of to becomes rapidly smaller and smaller suggesting that only a few first terms of the series-expansion of should lead to a satisfactory approximation. When , the first-order term is
[TABLE]
Note finally that, with and , the factors
[TABLE]
appearing in the series expansion of are the normalizing constants of a Dirichlet distribution .
6. Generalized Ewens sampling formulae
6.1. The weak mutation/selection potential case
Let
[TABLE]
whose support is the simplex We wish to consider sampling problems within (as a random partition of the unit interval ) describing the equilibrium distribution of the allelic frequencies, [29]. Such generalized Ewens sampling formulae (ESF) were also considered in [34]-[36], in this context.
6.1.1. Moments
We first consider the simpler problem of computing the moments of . With non-negative integers, the moments of this stationary distribution are:
[TABLE]
requiring the previous computation of . With , we find from Equation (42)
[TABLE]
Note that, from the identity
[TABLE]
[TABLE]
where is homogeneous of degree In particular
[TABLE]
The marginal moments are obtained from in (45), while considering except for
6.1.2. Generalized Ewens Sampling Formula
Take a random (uniform) sequential sample without replacement from describing the random equilibrium distribution of the allelic frequencies. Suppose there are distinct visited types of alleles in the process; suppose also that are the types of the visited alleles and that are the number of visits to (or hits of) alleles number , entailing and . We let ; shifting only the entries of by and leaving the other ones unchanged. From (44) and (42), we obtain:
Theorem 5**.**
With summing to and , the probability of such an occupancy event is
[TABLE]
[TABLE]
where is given by (42).
Considering the marginal event “”, it holds that
[TABLE]
where the sum runs over all positive integers summing to . There are terms in this sum. And summing the latter probability over gives
6.1.3. Generalized ESF in the ‘Kingman limit’, [23]
Let for each and while We call it the limit and we wish to consider the generalized Ewens sampling formula in this limiting situation corresponding to a case with infinitely many alleles per locus, [33], [11], [39].
For instance, take , where , with , and
With (48) and (49) mentioned below, we have
Theorem 6**.**
In the limit, the probability of the occupancy event: “”, with summing to exists and is given by
[TABLE]
where and are given by (49) and (48) respectively.
*Proof: *In the limit,
Consider first * * and split the sum * * into and where corresponds to those arrays summing to and with th row and th column equal to [math]. We have
\begin{array}[]{l}\mathcal{Z}\left(\mathbf{\theta}\right)=1+\frac{\Gamma\left(\left|\mathbf{\theta}\right|\right)}{\prod_{k=1}^{K}\Gamma\left(\theta_{k}\right)}\sum_{N\geq 1}\left(-1\right)^{N}\sum_{\left|\left[\mathbf{n}\right]\right|=N}\frac{\prod_{k=1}^{K}\Gamma\left(\theta_{k}+\alpha_{k}\left(\left[\mathbf{n}\right]\right)\right)}{\Gamma\left(\left|\mathbf{\theta}\right|+2N\right)}\prod_{k,l=1}^{K}\frac{\overline{W}_{k,l}^{n_{k,l}}}{n_{k,l}!}\\ =1+\frac{\Gamma\left(\left|\mathbf{\theta}\right|\right)}{\prod_{k=1}^{K}\Gamma\left(\theta_{k}\right)}\sum_{N\geq 1}\left(-1\right)^{N}(\sum_{\left|\left[\mathbf{n}\right]\right|=N}^{1}\frac{\prod_{k=1}^{K}\Gamma\left(\theta_{k}\right)}{\Gamma\left(\left|\mathbf{\theta}\right|+2N\right)}\prod_{k,l=1}^{K}\frac{\overline{W}_{k,l}^{n_{k,l}}}{n_{k,l}!}+\\ +\sum_{\left|\left[\mathbf{n}\right]\right|=N}^{2}\frac{\prod_{k=1}^{K}\Gamma\left(\theta_{k}+\alpha_{k}\left(\left[\mathbf{n}\right]\right)\right)}{\Gamma\left(\left|\mathbf{\theta}\right|+2N\right)}\prod_{k,l=1}^{K}\frac{\overline{W}_{k,l}^{n_{k,l}}}{n_{k,l}!})\end{array}
[TABLE]
Consider now * * and split the sum * * into the two parts and where corresponds to those arrays whose entries sum to and with th row and th column equal to [math], but only for those . We have
\begin{array}[]{l}\mathcal{Z}\left(\mathbf{\theta}+\mathbf{m}_{p}\right)=1+\frac{\Gamma\left(\left|\mathbf{\theta}\right|+m\right)}{\prod_{q=1}^{p}\Gamma\left(\theta_{k_{q}}+m_{q}\right)\prod_{k\neq\left\{k_{1},...,k_{p}\right\}}\Gamma\left(\theta_{k}\right)}\times\\ \sum_{N\geq 1}\left(-1\right)^{N}(\sum_{\left|\left[\mathbf{n}\right]\right|=N}^{1}\frac{\prod_{q=1}^{p}\Gamma\left(\theta_{k_{q}}+m_{q}+\alpha_{k_{q}}\left(\left[\mathbf{n}\right]\right)\right)\prod_{k\neq\left\{k_{1},...,k_{p}\right\}}\Gamma\left(\theta_{k}\right)}{\Gamma\left(\left|\mathbf{\theta}\right|+m+2N\right)}\prod_{k,l=1}^{K}\frac{\overline{W}_{k,l}^{n_{k,l}}}{n_{k,l}!}\\ +\sum_{\left|\left[\mathbf{n}\right]\right|=N}^{2}\frac{\prod_{q=1}^{p}\Gamma\left(\theta_{k_{q}}+m_{q}+\alpha_{k_{q}}\left(\left[\mathbf{n}\right]\right)\right)\prod_{k\neq\left\{k_{1},...,k_{p}\right\}}\Gamma\left(\theta_{k}+\alpha_{k}\left(\left[\mathbf{n}\right]\right)\right)}{\Gamma\left(\left|\mathbf{\theta}\right|+m+2N\right)}\prod_{k,l=1}^{K}\frac{\overline{W}_{k,l}^{n_{k,l}}}{n_{k,l}!})\end{array}
[TABLE]
The limit expression of the occupancy event generalizes the Ewens partition distribution in the infinitely-many-alleles population genetics model with symmetric selection/mutation, as studied in [13] and [17].
6.1.4. Marginal distributions and frequency spectrum
Let
[TABLE]
be the th marginal of as given by (43). Then
[TABLE]
as the frequency spectrum, is the density of alleles at equilibrium in a neighborhood of , see [9], [10]. It is a uniform mixture of the marginals whose precise computation is possible but involved. Note that, from (45), the marginal moments of are known.
6.2. The house-of-cards mutations potential case** **
We finally consider sampling formulae when the multi-allelic population is at equilibrium and subject only to mutation driving forces.
If is the mutation potential under the house of cards condition
[TABLE]
then, with
[TABLE]
[TABLE]
the Dirichlet distribution on the simplex with parameters ,
When dealing with weak mutations, the allelic equilibrium distribution is
[TABLE]
the asymmetric Dirichlet distribution on the simplex with parameters , (unless for some common in which case boils down to the symmetric Dirichlet distribution ). When dealing with the potential solely arising from mutations (i.e. when avoiding selection), the allelic equilibrium probability distribution is thus the one of , corresponding to Dirichlet spacings. A much simpler situation than the previous one involving both mutation and selection but still not so obvious ([9], section ). We shall derive some sampling formulae in this context.
6.2.1. The frequency of a typical allele at equilibrium
Pick at random an allele from . It has frequency with probability and so the frequency of a size-biased picked allele at equilibrium is
[TABLE]
an homogeneous functional of degree . Applying of Theorem 1
[TABLE]
The variance of requires the computation of the expected value of , as an homogeneous functional of degree , which is
[TABLE]
The full law of could be obtained while considering This functional is no longer homogeneous and of Theorem 1 should then be applied, with some combinatorics involved which we skip.
6.2.2. Smallest and largest allelic frequencies from
We shall now use of Theorem 1 to compute the joint distribution of the largest and smallest allelic frequencies in a distributed population at equilibrium. Suppose and consider the spacings’ functional
[TABLE]
Then is the required probability, assuming to be the order statistics of . The case () gives the probability respectively .
From statement of Theorem 1 indeed, the quantity
[TABLE]
interprets as the Laplace transform of Inverting this Laplace transform and putting yields . From this, we obtain directly
[TABLE]
where is the -fold convolution of the functions
[TABLE]
evaluated at . If , (54) gives the tail probability distribution of whereas gives the probability distribution of
We now briefly show in outline that these tools are also useful in the computation of simple sampling formulae.
6.2.3. Sampling and the Dirichlet multinomial distribution
Let be iid uniform throws on . Let be an integral-valued random vector which counts the number of visits to the different types of alleles in a sample. Hence, if is the allele type which the th trial meets, then ,
With and , conditionally given , we have the multinomial distribution:
[TABLE]
Averaging over and applying of Theorem 1 to compute , we find
[TABLE]
where , This distribution is known as the Dirichlet multinomial distribution.
Applying Bayes formula, the posterior distribution of given is determined by its density at a point on the simplex as
[TABLE]
This shows, as is well-known, that , where is obtained by shifting . In particular
[TABLE]
6.2.4. Pòlya urn sequence
This suggests the following recursive approach to the sampling formula where successive samples are now drawn from the corresponding iterative posterior distributions. More specifically, let be the types of the successive alleles thus drawn. Then,
[TABLE]
[TABLE]
[TABLE]
The joint distribution of reads
[TABLE]
The sequence is a Pòlya urn sequence.
6.2.5. Asymptotics of the occupancy vector
The joint conditional generating function of the full occupancy vector reads
[TABLE]
which is homogeneous with degree allowing to compute Further, with , gamma, , as above, using independence between and gamma and recalling
\begin{array}[]{l}\mathbf{E}\left(\prod_{k=1}^{K}u_{k}^{B_{k}/m}\right)=\frac{\Gamma\left(\left|\mathbf{\theta}\right|\right)}{\Gamma\left(\left|\mathbf{\theta}\right|+m\right)}\mathbf{E}\left[\left(\sum_{k=1}^{K}u_{k}^{1/m}T_{k}\right)^{m}\right]\\ \underset{m\uparrow\infty}{\sim}\frac{\Gamma\left(\left|\mathbf{\theta}\right|\right)}{\Gamma\left(\left|\mathbf{\theta}\right|+m\right)}\mathbf{E}\left[\left(\sum_{k=1}^{K}T_{k}\right)^{m}\left(1+\frac{1}{m}\sum_{k=1}^{K}\overline{T}_{k}\log u_{k}\right)^{m}\right]\\ \underset{m\uparrow\infty}{\sim}\frac{\Gamma\left(\left|\mathbf{\theta}\right|\right)}{\Gamma\left(\left|\mathbf{\theta}\right|+m\right)}\mathbf{E}\left(\sum_{k=1}^{K}T_{k}\right)^{m}\mathbf{E}\left(\prod_{k=1}^{K}u_{k}^{\overline{T}_{k}}\right)=\mathbf{E}\left(\prod_{k=1}^{K}u_{k}^{\overline{T}_{k}}\right)=\mathbf{E}\left(\prod_{k=1}^{K}u_{k}^{S_{k}}\right),\end{array}
because . This shows that
[TABLE]
Note that, applying the strong law of large numbers (conditionally given ), the above convergence in law also holds almost surely: the normalized occupancy vector from an sample converges to itself.
6.2.6. **ESF from equilibrium distribution driven solely by
mutations**
Firstly, with , as we already observed, it holds
[TABLE]
where if is an integer. Therefore, considering the probability of an occupancy event as in (46):
Proposition 7**.**
With summing to , the Ewens sampling formula under asymmetric mutations only is
[TABLE]
where
If we now consider the limit for mutations,
[TABLE]
showing that there is no proper limit of the occupancy probability (59).
Remark (standard Ewens sampling formula): Suppose , (the symmetric Dirichlet mutation model). Owing to exchangeability of the alleles
[TABLE]
where the occupancy vector is, say, over the first alleles of .
Consider the limit where , while , [23]. Owing to et , we get
[TABLE]
which is formula of [16], for example. In contrast with the asymmetric Dirichlet model, the symmetric Dirichlet model admits a proper limit occupancy probability. This is one of the many facets777In [16], the ESF (61) is called the first ESF. A second ESF rather deals with the occupancy vector , , which counts the number of alleles in the sample with representatives. of the standard ESF, see [16] and [11].
Coming back to (59) in the asymmetric Dirichlet case, this suggests to consider the following limit: choose in such a way that as , with none of the . Such a limiting model for the mutation rates was considered in [19].
Examples are:
, , for some with
, with
In the limit, therefore
[TABLE]
The asymmetric Dirichlet model has a proper limit occupancy probability (59), as given by (62).
To summarize:
Theorem 8**.**
Consider the sampling problem from . With summing to and , the probability of the occupancy event
[TABLE]
is given by (59). From (61) there is no non-degenerate limit of this probability in the limit but there is one in the limit, given by (62).
6.2.7. Moments, marginal distributions and frequency spectrum
Recalling (58), the marginal moments of are . So beta with marginal density equal to
[TABLE]
This gives the empirical average of in a neighborhood of as
[TABLE]
Else, , as the frequency spectrum, is the density of alleles at equilibrium in a neighborhood of It is a uniform mixture of beta distributed rvs, . Near , and , where and
Acknowledgments: T. Huillet acknowledges partial support both from the “Chaire Modélisation mathématique et biodiversité” and the labex MME-DII Center of Excellence (Modèles mathématiques et économiques de la dynamique, de l’incertitude et des interactions, ANR-11-LABX-0023-01 project).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] Barbour, A.D., Ethier, S.N. and Griffiths R.C. A transition function expansion for a diffusion model with selection. The Annals of Appl. Prob. Vol. 10, No. 1, 123-162, 2000.
- 2[2] Bürger, R. “The mathematical theory of selection, recombination, and mutation,” Wiley Series in Mathematical and Computational Biology. John Wiley & Sons, Ltd., Chichester, 2000.
- 3[3] Castilloux, A.M. and Lessard, S. The fundamental theorem of natural selection in Ewens’ sense (case of many loci). Theor. Popul. Biol. 48(3), 306-315, 1995.
- 4[4] Charlesworth, B. and Jain, K. Purifying selection, drift, and reversible mutation with arbitrarily high mutation rates. Genetics, 198(4), 1587-1602, 2014.
- 5[5] Edalat, A. Shahshahani Gradients. http://sharif.edu/~shahyad/edalat.pdf.
- 6[6] Ethier, S.N. and Kurtz, T.G. “Markov processes. Characterization and convergence”. John Wiley & Sons, Inc., New York, 1986.
- 7[7] Evans, S.N. Diffusions on the simplex from Brownian motions on hypersurfaces. Lecture Notes-Monograph Series. Institute of Mathematical Statistics. 35-48, 2003.
- 8[8] Ewens, W.J. “Mathematical population genetics. I. Theoretical introduction,” Second edition. Interdisciplinary Applied Mathematics, 27. Springer-Verlag, New York, 2004.