# On Discrete-Time Multiallelic Evolutionary Dynamics Driven by Selection

Authors: Thierry E. Huillet
Journal: Journal of Probability and Statistics
 (2010)

---

## Abstract

We revisit some problems arising in the context of multiallelic discrete-time evolutionary dynamics driven by fitness. We consider both the deterministic and the stochastic setups and for the latter both the Wright-Fisher and the Moran approaches. In the deterministic formulation, we construct a Markov process whose Master equation identifies with the nonlinear
deterministic evolutionary equation. Then, we draw the attention on a class of fitness matrices that plays some role in the important matter of polymorphism: the class of strictly ultrametric fitness matrices. In the random cases, we focus on fixation probabilities, on various conditionings on nonfixation, and on (quasi)stationary distributions.


---

## Body

## 1. Introduction

Population genetics aims at elucidating the fate of genotype frequencies undergoing the basic evolutionary processes when various driving “forces” such as fitness, mutation, or recombination are at stake in the gene pool. This requires to clarify the updating mechanisms of the gene frequency-distributions over time. Another important additional driving source is the genetic drift whose nature is exclusively random. The corresponding field of interest is the statistical theory arising from this aspect of the gene replacement processes and it requires some use of the Markov chain theory (see [1]).In this note, we revisit the basics of both the deterministic and stochastic dynamics arising in discrete-time asexualmultiallele evolutionary genetics driven only by fitness. We do not touch at all upon other important mechanisms such as mutation. We start with the haploid case with K alleles before switching to the more interesting diploid case.Let us summarize and comment the material developed in Section2. In the deterministic haploid case, a vector of fitness is attached to the alleles. The interest is on the evolution of the allelic frequency distribution over time. The updates of the allele frequency distributions are driven by the relative fitnesses of the alleles, ending up in a state where only the fittest monomorphic state will survive. This state is also an extremal point of the simplex over which the dynamics takes place. From this dynamics, it appears that the mean fitness increases as time passes by, the rate of increase being the variance in relative fitness (a well-known particular incarnation of the Fisher theorem of natural selection). We recall the background and the evolution equations.The haploid replicator dynamics, as a nonlinear updating mapping from the simplex to the simplex, may be viewed as the discrete-time nonlinear Master equation of some random Markov process. We supply a construction for this process thereby giving a stochastic interpretation to the original deterministic formulation of the dynamics.In the diploid case, there is a similar deterministic updating dynamics but now on the full array of the genotype frequencies. It involves the fitness matrix attached to the genotypes. When mating is random so that the Hardy-Weinberg law applies, we may look at the induced marginal allelic frequencies dynamics. The updating dynamics looks quite similar to the one occurring in the haploid case except that the mean fitness is now the mean fitness quadratic form in the current frequencies whereas marginal fitnesses are no longer constant but affine functions in these frequencies. As for the haploid case, it is possible to construct a Markov process whose nonlinear Master equation coincides with the diploid replicator dynamics. We supply this construction which we believe is new.In the diploid context, the Fisher theorem still holds true but, as a result of the fitness landscape being more complex, there is a possibility for a polymorphic equilibrium state to emerge. Due to its major evolutionary interest, our subsequent concern is to identify examples of diploid dynamics leading to a uniquepolymorphic state on the simplex, either unstableorstable. We start with the unstable case and draw the attention on a class of fitness matrices leading to a unique unstable polymorphic equilibrium state: the class of strictly potential matrices. Strictly ultrametric matrices are particular instances of strictly potential matrices [2] which therefore will display unstable polymorphism as well. There is a useful canonical representation of strictly ultrametric matrices due to [3] which we recall which helps giving specific examples of strictly ultrametric matrices. When dealing with the class of strictly potential fitness matrices, the mean fitness quadratic form is definite-positive; we derive a related class of fitness matrices leading to a definite-negative mean fitness quadratic form. For this class of matrices, there will also be a unique polymorphic equilibrium state for the diploid dynamics and it will be stable. We also draw the attention on a subclass of the latter: the so-called “antistrictly ultrametric matrices.” For such matrices, among other things, the fitness of each homozygote should not exceed the ones of all the heterozygotes carrying the allele of the homozygote. To the best of our knowledge, the large class of potential fitness matrices as natural candidates for polymorphic states to emerge was not discussed in the literature.Section3 is devoted to the stochastic version of these considerations when the transitions in the constitutive allelic population sizes are given by a K-dimensional Wright-Fisher model with total constant-size N (see [1, 4]). It takes into account an additional important driving source of evolution, namely, the genetic drift, whose nature is exclusively random. Under fitness only and in particular in the absence of mutations, the multiallelic Wright-Fisher model is a transient Markov chain on a (N+K-1K-1)-dimensional state-space whose absorbing states are the monomorphic states. We give an expression for the fixation probabilities of this process. Then, we develop four conditioning problems: conditioning on fixating in a given monomorphic state, conditioning on avoiding the extremal states before the current instant, conditioning on nonfixation at each transition time, and conditioning on avoiding the extremal states in the remote future. For the last three conditioned processes, the equilibrium structure of the fitness mechanism shows up.Finally, we run into similar considerations for theMoran model. When dealing with the fixation probabilities in this Moran context, we suggest a new mean-field approximation of these probabilities which is based on a well-known explicit formula for the 2-alleles case [1]. It concerns the case of multiplicative fitnesses only. Finally, we consider the Moran model conditioned on nonfixation at each transition time. We exploit the reversible character of this process to derive a new explicit product formula for its invariant probability measure.


## 2. Deterministic Evolutionary Dynamics

We start with the haploid case before moving to the more interesting diploid case.


### 2.1. Single Locus: Haploid Population withK Alleles

ConsiderK alleles Ak,k=1,…,K attached to a single locus. Suppose that the current time-t allelic frequency distribution is given by the column vector x:=xk,k=1,…,K. (Throughout, a boldface variable, say x, will represent a column-vector so that its transpose, say x*, will be a line-vector.) We therefore have x∈SK={x≥0:|x|:=∑k=1Kxk=1}, the K-simplex as a convex subset of ℝK with dimension K-1. Let w:=wk>0,k=1,…,K denote the absolute fitnesses of the alleles. Let
(2.1)w(x):=∑lwlxl=w*x
be the mean fitness of the population at time t. We will also need
(2.2)σ2(x)=∑k=1Kxk(wk-w(x))2,
the variance in absolute fitness, and(2.3)σ¯2(x)=∑k=1Kxk(wkw(x)-1)2=σ2(x)w(x)2,
the variance in relative fitness wk/w(x).


#### 2.1.1. Dynamics

From the deterministic evolutionary genetics point of view, the discrete-time update of the allele frequency distribution on the simplexSK is given by (2.4)xk′=pk(x):=xkwkw(x),k=1,…,K.
(the symbol ′ is a common and useful notation to denote the updated frequency.) The quantity wk/w(x)-1 therefore interprets as the frequency-dependent Malthus growth rate parameter of xk. As required, the vector p(x):=pk(x),k=1,…,K, maps SK into SK. In vector form, with Dx:=diag(xk,k=1,…,K), the nonlinear deterministic evolutionary dynamics reads
(2.5)x′=p(x)=1w(x)Dwx=1w(x)Dxw,
or, with Δx:=x′-x, the increment of x(2.6)Δx=(1w(x)Dw-I)x.
Avoiding the trivial case where fitnesses are all equal, without loss of generality, we can assume that either w1≥⋯≥wK=1 or w1≤⋯≤wK=1. Thus that allele A1 or AK has largest fitness.


#### 2.1.2. Mean Fitness Increase

According to the dynamical system (2.4), unless the equilibrium state is attained, the absolute mean fitness w(x) increases:
(2.7)Δw(x)=w(x′)-w(x)=∑kwkΔxk=∑kwkxk(wkw(x)-1)=∑kwk2xkw(x)-w(x)>0.
The mean fitness is maximal at equilibrium. The rate of increase of w(x) is(2.8)Δw(x)w(x)=∑kxk(wkw(x)-1)2=∑k(Δxk)2xk,
which is the variance in relative fitness σ¯2(x) defined in (2.3). These last two facts are sometimes termed the 1930s Fisher fundamental theorem of natural selection. The equilibria of (2.4) are the extremal states (0-faces) of the boundary of SK. To make it simple, if there is an allele whose fitness is strictly larger than the ones of the others, the deterministic evolutionary dynamics (2.4) will attain an equilibrium where only the fittest will survive; starting from any initial state of SK which is not an extremal (or monomorphic) point, the haploid trajectories will converge to this fittest state.


#### 2.1.3. A Stochastic Interpretation of the Deterministic Dynamics (2.4)

A vectorx of SK can be thought of as a probability vector. The dynamical equation (2.4), as a nonlinear update mapping from SK to SK, may be viewed as the discrete-time nonlinear Master equation of some Markov process whose construction we now give. Suppose that we have a population of N haploid individuals each of which can be of one among K types or colors (carrying one among the K possible alleles). We will need to introduce an extra color-state, say ∂={0}, which will be absorbing for the process we will now construct. Let K(t):=Kn(t),n=1,…,N be the random color distribution of the population at time t, therefore with enlarged state-space {0,1,…,K}. Assume that the individuals are indistinguishable leading to the exchangeability property: Kn(t)=dK1(t),n=2,…,N (equality in distribution). Let Ut,m,t=1,2,…;m=1,…,N be N i.i.d. driving sequences of uniformly distributed random variables on [0,1], independent of K(t). Let w¯k:=wk/∑kwk,k=1,…,K. To decide the allele Km(t+1) carried by the individual number m∈{1,…,N} at time t+1, with 1(A) being the indicator function of the event A, consider the random Markovian dynamical system:(2.9)1(Km(t+1)=k,τm>t+1)=1(w¯kN∑n=1N1(Kn(t)=k,τn>t)>Ut+1,m).
Here k∈{1,…,K} and τn is the first time that Kn(t) hits the absorbing state ∂. As a result of Kn(t)=dK1(t), we naturally assume τn=dτ1,n=2,…,N. For each n, our model therefore attributes a positive probability that Kn(t)=0 for all t≥τn. Although, in principle, there is a possibility that the type of the mth particle is the one of the fictitious unobservable allele A0, as a result of (2.9), the sample paths of K(t) leading to this A0 are ruled out because focus is on the observable states only.In words, for the dynamics (2.9), the observable event Km(t+1)=k is realized (together then with τm>t+1) if the proportion at t of type-k individuals, weighted by the corresponding scaled fitness w¯k, is large enough (compared to Ut+1) and of course if the whole process was not absorbed at {0} in the previous step. Taking first the expectation with respect to the driving noise Ut+1,m in (2.9), we get
(2.10)ℙ(Km(t+1)=k,τm>t+1∣K(t)>0)=w¯kN∑n=1N1(Kn(t)=k,τn>t).
Putting zk(t):=ℙ(K1(t)=k,τ1>t), recalling (Kn(t);τn)=d(K1(t);τ1),n=2,…,N, and taking the expectation with respect to K(t), we get an unnormalized version of (2.4):
(2.11)zk(t+1)=w¯kzk(t),k∈{1,…,K}.
We have 1>ℙ(τ1>t)=∑k=1Kzk(t)=∑k=1Kxk(0)w¯kt→0, geometrically fast. Defining the normalized conditional probabilities
(2.12)xk(t)=zk(t)∑k=1Kzk(t)=ℙ(K1(t)=k∣τ1>t),
we obtain the normalized haploid dynamics (2.4):
(2.13)xk′=wkxk∑k=1Kwkxk,k∈{1,…,K}.
It now may be viewed as the nonlinear Master equation of some stochastic Markov process. Let us make some miscellaneous remarks.(i) Clearly this construction makes also sense if N=1 (a single particle). (ii) When N is finite, we should stress that the initial condition can be chosen to be deterministic, say with xk(0)=zk(0)=ik/N, for some sequence of integers ik≥0 satisfying i1+⋯+iK=N (i0=0) and quantifying the initial population sizes. It could also be chosen to be random, with xk(0) defining the initial probability distribution of the alleles. This occurs in the large N limit if ik=⌊Nxk(0)⌋ so that ik/N→xk(0). The latter choice may therefore be interpreted as a large N limit of the former one. (iii) In the stochastic interpretation (2.9) of the deterministic dynamics (2.4), xk(t) can be interpreted either as the probability that the random allele carried by a typical individual is Ak or like the expected proportion of the individuals of type k within the whole population (a frequentist point of view). (iv) The appeal to the coffin state ∂ was a necessary step to understand the normalization zk→xk.(v) Even though K(t) is exchangeable, it is not true that, with n1≠n2 any two distinct individuals, their random labels Kn1(t) and Kn2(t) are independent. The random algorithm allowing to update the joint types of Kn1(t) and Kn2(t) could be written down but is much more involved.


### 2.2. Single Locus: Diploid Population withK Alleles

We now run into similar considerations but with diploid populations.


#### 2.2.1. Joint Evolutionary Dynamics

Letwk,l≥0,k,l=1,…,K stand for the absolute fitness of the genotypes AkAl attached to a single locus. Assume wk,l=wl,k (wk,l being proportional to the probability of an AkAl surviving to maturity, it is natural to take wk,l=wl,k). Let then W be the symmetric fitness matrix with k,l-entry wk,l.Assume that the current frequency distribution at timet of the genotypes AkAl is given by xk,l. Let X be the frequencies array with k,l-entry xk,l. The joint evolutionary dynamics in the diploid case is given by the updating:(2.14)xk,l′=xk,lwk,lω(X),
where the mean fitness ω is now given by ω(X)=∑k,lxk,lwk,l. The relative fitness of the genotype AkAl is wk,l/ω(X). The joint dynamics takes the matrix form:
(2.15)X′=1ω(X)X∘W=1ω(X)W∘X,
where ∘ stands for the (commutative) Hadamard product of matrices.LetJ be the flat K×K matrix whose entries are all 1. Then
(2.16)ΔX:=X′-X=1ω(X)(X-J)∘W=1ω(X)W∘(X-J).
We will also let(2.17)σ2(X)=∑k,l=1K‍xk,l(wk,l-ω(X))2
stand for the genotypic variance in absolute fitness and(2.18)σ¯2(X)=∑k,l=1Kxk,l(ωk,lω(X)-1)2=σ2(X)ω(X)2
will stand for the diploid variance in relative fitness.Consider the problem of evaluating the increase of the mean fitness. We have(2.19)Δω(X)=∑k,lΔxk,lwk,l=∑k,lxk,l(wk,l2ω(X)-wk,l)=ω(X)σ¯2(X)>0
with a relative rate of increase: Δw(X)/w(X)=σ¯2(X). This is the full diploid version of the Fisher theorem.


#### 2.2.2. Marginal Allelic Dynamics

Assuming a Hardy-Weinberg equilibrium, the frequency distribution at timet, say xk,l, of the genotypes AkAl is given by xk,l=xkxl where xk=∑lxk,l is the marginal frequency of allele Ak in the whole genotypic population . The whole frequency information is now enclosed within x:=xk,k=1,…,K. For instance, the mean fitness is now given by the quadratic form: ω(x):=∑k,lxkxlwk,l=x*Wx, with x* being the transposed line vector of the column vector x=X1 (1 the unit K-vector). We will also let(2.20)σ2(x)=∑k,l=1Kxkxl(wk,l-ω(x))2
stand for the genotypic variance in absolute fitness and(2.21)σ¯2(x)=∑k,l=1Kxkxl(wk,lω(x)-1)2=σ2(x)ω(x)2
will stand for the diploid variance in relative fitness.Consider now the update of the allelic marginal frequenciesx themselves. If we first define the frequency-dependent marginal fitness of Ak by wk(x)=(Wx)k:=∑lwk,lxl, the marginal dynamics is given as in (2.4) by:(2.22)xk′=xkwk(x)ω(x)=1ω(x)xk(Wx)k=:pk(x),k=1,…,K.
This dynamics involves a multiplicative interaction between xk and (Wx)k, the kth entry of the image Wx of x by W. In (2.22) there is a normalization by the quadratic form ω(x)=x*Wx. In vector form (2.22) reads
(2.23)x′=1ω(x)DxWx=1ω(x)DWxx=:p(x),
where p maps SK into SK. Iterating, the time-t frequency distribution is(2.24)x(t)=p(p(⋯ttimes⋯(p(x0)))).Except for the fact that the mean fitnessω in (2.22) is now a quadratic form in x and that the marginal fitness of Ak is now frequency-dependent, depending linearly on x, as far as the marginal frequencies are concerned, the updating formalism (2.22) in the diploid case looks very similar to the one in (2.4) describing the haploid case.In the diploid case, assuming fitnesses to bemultiplicative, say with wk,l=wkwl, then wk(x)/x*Wx=wk/∑lwlxl and the dynamics (2.22) boils down to (2.4). However, the mean fitness in this case is ω(x)=(∑lwlxl)2 and not w(x)=∑lwlxl as in the haploid case.


#### 2.2.3. A Stochastic Interpretation of the Deterministic Dynamics (2.22)

As for the haploid case, there is a Markov chain governed by the Master equation (2.22). Consider a population of diploid individuals. The number of alleles N in this population is therefore twice the number of genes. Each allele can be of one among K types or colors (carrying one among the K possible alleles). As before, we introduce an extra color-state, say ∂={0} which is absorbing for the process to be constructed. For c=1,2, let Kc(t):=Knc(t),n=1,…,N be two independent copies of the random color distribution of the allelic population at time t. Let K(t)=(K1(t),K2(t)). Assume that the alleles are indistinguishable within each sample, leading to Knc(t)=dK1c(t),n=2,…,N,c=1,2. For c=1,2, let Ut,mc,t=1,2,…;m=1,…,N be two mutually independent i.i.d. driving N-sequences of uniformly distributed random variables on [0,1] and independent of K(t). To decide the type of the random allele Km(t+1),m=1,…,N, at time t+1, consider now the Markovian dynamical system:(2.25)1(Km(t+1)=k,τm>t+1)=1(1N∑n=1N1(Kn1(t)=k,τn1>t)>Ut+1,m1)1(1N∑n=1NW¯k,Kn2(t)1(τn2>t)>Ut+1,m2).
Here k∈{1,…,K},τnc are the first hitting times of each Knc(t) of the absorbing state ∂,c=1,2 and W¯:=W/∥W∥ for any matrix norm ∥W∥, say, for example, ∥W∥=∑k,lwk,l. We assume τnc=dτ1c,n=2,…,N,c=1,2.In words, for this new dynamics, the observable eventKm(t+1)=k is seen to be realized together with τm>t+1 if two natural independent conditions are now satisfied which can be read from the two indicator functions in the right-hand side of (2.25):(i)
first, the proportion att of type-k individuals of the first copy should be large enough (compared to Ut+11);(ii)
second, for the second sample copyK2(t), the expected fitness of the genotypes AkAl,l=1,…,K, containing allele Ak at t should be large enough (compared to Ut+12).As for the haploid case, it is necessary that both processesK(t) should not be absorbed at {0} in the previous step. Taking first the expectation with respect to the independent driving noises Ut+1,mc in (2.25), we get
(2.26)ℙ(Km(t+1)=k,τm>t+1∣K(t)>0)=1N∑n=1N1(Kn1(t)=k,τn1>t)·1N∑n=1NW¯k,Kn2(t)1(τn2>t).
Putting zk(t)=ℙ(K1(t)=k,τ1>t), taking the expectation with respect to K(t) and using our independence and exchangeability hypotheses, we get
(2.27)zk(t+1)=zk(t)·𝔼[W¯k,K12(t)1(τ12>t)]=zk(t)·∑l=1kW¯k,lzl(t),
corresponding to an unnormalized version of (2.22):
(2.28)zk(t+1)=zk(t)(W¯z(t))k,k∈{1,…,K}.
Defining the normalized conditional probabilities
(2.29)xk(t)=zk(t)∑k=1Kzk(t)=ℙ(K1(t)=k∣τ1>t),
we obtain the normalized nonlinear Markovian diploid dynamics (2.22):
(2.30)xk′=xk(Wx)k∑k=1Kxk(Wx)k,k∈{1,…,K}.
Note that this construction makes sense if N=2 (a single individual and its 2 alleles of K possible types). The need for two copies of K(t) was governed by the quadratic character of the interaction appearing in the numerator of (2.22).


#### 2.2.4. Increase of Mean Fitness

Again, the mean fitnessω(x), as a Lyapunov function, increases as time passes by. We indeed have
(2.31)Δω(x)=ω(x′)-ω(x)=1ω(x)2∑k,lxkwk(x)wk,lxlwl(x)-∑k,lxkwk,lxl>0,
because, defining 0<X(x):=∑k,lxk(1-wk(x)/ω(x))wk,l(1-wl(x)/ω(x))xl, we have
(2.32)Δω(x)=X(x)+2ω(x)(∑kxkwk(x)2-ω(x)2)>0.
Its partial rate of increase due to frequency shifts only is δω(x):=∑kΔxkwk(x). It satisfies(2.33)δω(x)ω(x)=∑kxk(wk(x)ω(x)-1)2=∑k(Δxk)2xk=σ¯A2(x)2,
where σ¯A2(x) is the allelic variance in relative fitness:(2.34)σ¯A2(x):=2∑k=1Kxk(wk(x)ω(x)-1)2.


#### 2.2.5. An Alternative Representation of the Allelic Dynamics

There is an alternative vectorial representation of the dynamics (2.22). Define the symmetric positive-definite matrix G(x)=Dx(I-1x*) with quadratic entries in the frequencies:
(2.35)G(x)k,l=xk(δk,l-xl).
Introduce the quantity VW(x)=(1/2)logω(x), which is half the logarithm of the mean fitness. Then, (2.22) may be recast as the gradient-like dynamics:(2.36)Δx=1ω(x)G(x)Wx=G(x)∇VW(x),
with |Δx|=1*Δx=0 as a result of 1*G(x)=0*. Note
(2.37)∇VW(x)*Δx=∇VW(x)*G(x)∇VW(x)≥0.
The dynamics (2.36) is of pure gradient-type with respect to the Svirezhev-Shashahani distance metric dG(x,x′); see [5, 6]. For this metric, the distance between x and x′=x+Δx of SK is
(2.38)dG(x,x′)=(Δx*G-1Δx)1/2=(∑k=1Kxk-1(Δxk)2)1/2.
From (2.33) and (2.34), this quantity, which is the length of Δx satisfying |Δx|=0, is also the square-root of half the allelic variance (the standard deviation) in relative fitness.


### 2.3. Equilibria (Diploid Case)

The mean fitness increase phenomenon occurs till the evolutionary dynamics reaches an equilibrium state. We wish to briefly discuss the questions relative to equilibria in the diploid case.


#### 2.3.1. Preliminaries

In contrast with the haploid case, in the diploid situation, the dynamics (2.22) can have more complex equilibrium points, satisfying wk(xeq)=w1(xeq),k=2,…,K, and ∑lxeq,l=1. To avoid linear manifolds of equilibria, we first assume that all principal minors of W are nonsingular and also that the fitnesses of all homozygotes wk,k are positive. In this case, from the Bézout theorem, the number of equilibria is finite and less or equal than the number 2K-1 of faces of SK. Note that the K extremal endpoints of SK (0-faces) are always monomorphic fixed points of (2.22).An instructive example fulfilling these preliminary conditions isW=I. There are 2K-1 equilibrium points (the barycenters of the (Kk+1)k-dimensional faces, k=0,…,K-1), but only one polymorphic equilibrium which is the barycenter xB of SK. This point is the one with minimal fitness and it is unstable. The 0-faces are stable fixed points whereas the barycenters of the k faces with k∈{1,…,K-2} are saddle-points. The simplex SK could be partitioned into pieces each of which being the attraction basins of the stable 0-face states: in contrast with the haploid case, the type of the survivor is not necessarily the one of the fittest; it will depend on the initial condition.Similar conclusions can be drawn if instead ofW=I, we start with W=(I-Dλ)-1 where λ:=(λk,k=1,…,K) satisfies 0≤λ<1 (meaning 0≤λk<1, for all k). In this case again, there is only one unstable polymorphic equilibrium which is easily seen to be xeq=(1-λ)/(K-|λ|)∈SK.Due to its evolutionary interest, we would like now to discuss the opportunity for a polymorphic state to be stable. Under the above assumptions onW, a unique stable internal (polymorphic) equilibrium state can exist, necessary and sufficient conditions being that (i) there is a unique z>0 for which Wz=1 and (ii)W has exactly one strictly positive dominant eigenvalue and at least one strictly negative eigenvalue or else the sequence of principal minors of W alternates in sign (see Kingman [7]). If this is the case, the equilibrium polymorphic state is xeq=z/|z|. It is stable in the sense that it is a local maximum of the mean fitness ω(x)=x*Wx, with ω(xeq)=1/|z|. Since |W|≠0, the linearization of p(x) at xeq has no eigenvalue of modulus 1 and so xeq is hyperbolic and/or isolated (see [8]). A stable isolated polymorphic state is asymptotically Lyapunov stable.Under these additional assumptions therefore onW, starting from any initial condition in the interior of SK, all trajectories will be attracted by xeq=z/|z| which is an isolated global maximum of ω(x).When there is no such unique globally stable polymorphic equilibrium, all trajectories will still converge but perhaps to a local equilibrium state where some alleles get extinct. Which allele and how many alleles are concerned seems to be an unsolved problem in its full generality.


#### 2.3.2. Special Classes of Fitness Matrices Leading to a Polymorphic State

We now draw the attention on a particular class of fitness matrices that lead to a polymorphic state, either unstable or stable. We start with the unstable case, extending the above special diagonal caseW=(I-Dλ)-1 leading to a unique unstable polymorphic state.(i) The Unstable Case
LetΛ≥0 be a symmetric irreducible strictly substochastic matrix satisfying Λ1:=q<1: the positive mass-defect vector of Λ is 1-q. Let λ>0. Define the symmetric strictly potential matrix: W=λ-1(I-Λ)-1≥0, with W-1=λ(I-Λ) defining a strictly row-diagonally dominant Stieltjes matrix with the properties [3]: (W-1)k,k>0,(W-1)k,l≤0 for k≠l and (W-1)k,k+∑l≠k(W-1)k,l>0 for all k. Then
(2.39)λ1>z=W-11=λ(I-Λ)1=λ(1-q)>0.
The vector z=λ(1-q) is called the equilibrium potential of W. We have |z|=λ(K-|q|). For this class of W therefore, Wz=1 admits a positive solution z.
Conversely, given a nonsingular matrixW≥0 satisfying Wz=1 for some z≥0, the matrix Λ=I-λ-1W-1 defines a substochastic matrix if and only W-1 satisfies (W-1)k,k>0,(W-1)k,l≤0 for k≠l and λ≥maxk(W-1)k,k. Then, W-1 is row-diagonally dominant and W=λ-1(I-Λ)-1 is a potential matrix.Strictly ultrametric (sUm) matrices are special classes of positive-definite and symmetric strictly potential matrices [2]. An sUm matrix W is symmetric with nonnegative entries, satisfying (i) wk,l≥min{wk,j,wj,l}, for all j,k,l and (ii) wk,k>maxl≠k{wk,l}, for all k (If in condition (ii),≥ is substituted for >, then W is simply an ultrametric matrix and this new condition is implied by (i)). If W is an sUm matrix, the fitness dynamics will admit an unstable polymorphic equilibrium state, as a result of W being positive-definite.Remark 2.1.
Suppose thatΛ is substochastic and primitive. Then W=λ-1(I-Λ)-1>0 is an ultrametric matrix. If V is the Hadamard reciprocal of W with entries vk,l=1/wk,l, it satisfies that vk,l≤max{vk,j,vj,l}. Therefore V is an ultrametric distance associated to the ultrametric potential W. Tree matrices are ultrametric matrices that are not sUm.(ii) The Stable Case
To produce a stable equilibrium state from the sUm matrix construction, letW=λ-1(I-Λ)-1 define a symmetric strictly potential matrix as before. Then, there is a z=λ(1-q)>0 for which Wz=1. With α>1, define
(2.40)W̃:=α|z|J-W.
With W̃=[w̃k,l], we have minw̃k,l=α/|z|-maxwk,l and we can choose α>1 so that W̃≥0. We have (α-1)W̃z=1 and δ*W̃δ=-δ*Wδ<0 for all δ≠0 satisfying |δ|=0 showing that xeq:=z/|z|=(1-q)/(K-|q|) is now a stable polymorphic state for W̃. If W is an sUm matrix, then clearly W̃ satisfies the “anti-sUm” property expressing a fitness domination of the heterozygotes AkAl over the homozygotes:
(2.41)w̃k,l≤max{w̃k,j,w̃j,l},∀j,k,l,w̃k,k<minl≠k{w̃k,l},∀k.
Nonnegative symmetric negative-definite anti-sUm fitness matrices W̃ will therefore display a stable polymorphic equilibrium state xeq. Note that xeq*W̃xeq=(α-1)/|z| is now the maximal value of the mean fitness.Example 2.2.
WhenK=2, with s>-1, h>0, and sh>-1, let
(2.42)W=[1+s1+sh1+sh1]>0
define the fitness matrix with selection parameter s and dominance h. This W is sUm if and only if s<0 and h>1. The equilibrium state is xeq*:=(h/(2h-1);(h-1)/(2h-1)) and it is unstable. This W is anti-sUm if and only if s>0 and h>1. The equilibrium state is the same but it is now stable.Note that a singular multiplicative fitness matrix of the formW=ww* cannot be a strictly potential matrix because its determinant |W| is zero.


#### 2.3.3. A General Construction to Produce sUm and Anti-sUm Matrices

Consider the problem consisting in splitting binarily and recursively the set{1,…,K} till complete reduction to singletons (leaves) which are left behind in the process. For instance, consider the refinement sequence with K=6 of {1,…,6}≡(123456): (123456)→((23)(1456))→(((2)(3))((16)(45)))→(((1)(6))((4)(5))).Starting from the left, there are2K-1=11 blocks of symbols (the total number of nodes in the splitting binary tree with K leaves). To each encountered block, numbered from l=1 to 2K-1, starting from the left, attach a vector ul of size K with ith entry ul(i)=1 if symbol i is in the block string,0 otherwise. For instance, from the above sequence, u1*=(1,1,1,1,1,1),u2*=(0,1,1,0,0,0),u3*=(1,0,0,1,1,1),…, and u10*=(0,0,0,1,0,0),u11*=(0,0,0,0,1,0). To each such ul, attach a number sl which is >0 if |ul|=1 (the leaves) and ≥0 if |ul|>1 (the internal nodes, including the root). Then (see [3]), for any choice of sl respecting these constraints
(2.43)W=∑l=12K-1slulul*≥0
is an sUm matrix and any sUm matrix can be represented in this way. Because for the ul corresponding to the leaves sl>0, the diagonal terms of W are necessarily >0.Since for each set{1,…,K}, there are bK splitting tree sequences where bK satisfies bK=∑k=1K-1(Kk)bkbK-k,k≥2,b1=1, there are many ways to generate an sUm matrix.Clearly, for each splitting procedure, withλ-1:=∑l=12K-1sl>0,W may be written as (2.44)W=λ-1(J-∑l=22K-1‍λsl(J-ulul*))≥0,
where the matrices J-ulul* take values in {0,1}.Now, withγ-1>0, any matrix of the form (2.45)W̃=(λ-1+γ-1)J-W=γ-1(J+∑l=22K-1γsl(J-ulul*))≥0
is an anti-sUm matrix. Assuming γ=1 and sl=shl,W̃ may be written under the form: W̃=J+sA where A:=∑l=22K-1hl(J-ulul*). It has at most 2K-1 independent parameters, namely, the hl,l=2,…,2K-1, and s. If the hl are known, then W̃ is a one-parameter family of fitness matrices. (Fitness matrices of the form J+sA were considered in [9] in the context of the estimation of s problem).Examples 2.3.
Assume also thatsl=s/(2(K-1)), for all l, with s>0 being a selection parameter. Then, defining the (0,1]-valued matrix
(2.46)A:=12(K-1)∑l=22K-1(J-ulul*)>0,
an anti-sUm matrix of the form W̃=J+sA will admit a stable polymorphic equilibrium. Clearly, A itself is anti-sUm. Because of this, there is a zA>0 such that AzA=1. Thus (J+sA)zA=(|zA|+s)1 showing that, with z=zA/(|zA|+s),W̃z=(J+sA)z=1. We thus have |z|=|zA|/(|zA|+s) and so xeq=z/|z|=zA/|zA|. Furthermore, the equilibrium mean fitness for such models is xeq*W̃xeq=1/|z|=(|zA|+s)/|zA|>1.
For the following simple sequence example withK=4 (four alleles, say A,C,T,G), (1234)→((1)(234))→((24)(3))→((2)(4)) we find
(2.47)A=16[5666635465456453]
which itself clearly is a symmetric anti-sUm matrix, together with W̃=J+sA. For this example, xeq*=(1/13)(8,1,3,1) and the equilibrium mean fitness is 1+s/13.
Note that takingsl=s/(2(K-1)) just for the indices l that were initially chosen to satisfy sl>0 would also lead to anti-sUm matrices A>0 and W̃=J+sA>0. In this case, the sum (2.46) defining A should then be restricted to the indices l from 2 to 2K-1 for which sl>0. Proceeding in this extreme way for the above simple example, we get the borderline anti-sUm shapes:
(2.48)A=16[3444434444344443],W̃=[s+223+2s33+2s33+2s33+2s3s+223+2s33+2s33+2s33+2s3s+223+2s33+2s33+2s33+2s3s+22].
Here, xeq* is simply the barycenter of the 4-simplex.Remark 2.4.
The ultrametric conditions (2.41) should be compared with the so-called “triangle inequality” conditions (leading to stable polymorphism) pointed out in [10], which read
(2.49)w̃k,l<w̃k,j+w̃l,j,∀k≠landatleastonej≠k,j≠l,w̃k,l>w̃k,k+w̃l,l2,∀k≠l.
It is clear that the class of anti-sUm matrices is a particular subclass of the Lewontin one.


## 2.1. Single Locus: Haploid Population withK Alleles

ConsiderK alleles Ak,k=1,…,K attached to a single locus. Suppose that the current time-t allelic frequency distribution is given by the column vector x:=xk,k=1,…,K. (Throughout, a boldface variable, say x, will represent a column-vector so that its transpose, say x*, will be a line-vector.) We therefore have x∈SK={x≥0:|x|:=∑k=1Kxk=1}, the K-simplex as a convex subset of ℝK with dimension K-1. Let w:=wk>0,k=1,…,K denote the absolute fitnesses of the alleles. Let
(2.1)w(x):=∑lwlxl=w*x
be the mean fitness of the population at time t. We will also need
(2.2)σ2(x)=∑k=1Kxk(wk-w(x))2,
the variance in absolute fitness, and(2.3)σ¯2(x)=∑k=1Kxk(wkw(x)-1)2=σ2(x)w(x)2,
the variance in relative fitness wk/w(x).


### 2.1.1. Dynamics

From the deterministic evolutionary genetics point of view, the discrete-time update of the allele frequency distribution on the simplexSK is given by (2.4)xk′=pk(x):=xkwkw(x),k=1,…,K.
(the symbol ′ is a common and useful notation to denote the updated frequency.) The quantity wk/w(x)-1 therefore interprets as the frequency-dependent Malthus growth rate parameter of xk. As required, the vector p(x):=pk(x),k=1,…,K, maps SK into SK. In vector form, with Dx:=diag(xk,k=1,…,K), the nonlinear deterministic evolutionary dynamics reads
(2.5)x′=p(x)=1w(x)Dwx=1w(x)Dxw,
or, with Δx:=x′-x, the increment of x(2.6)Δx=(1w(x)Dw-I)x.
Avoiding the trivial case where fitnesses are all equal, without loss of generality, we can assume that either w1≥⋯≥wK=1 or w1≤⋯≤wK=1. Thus that allele A1 or AK has largest fitness.


### 2.1.2. Mean Fitness Increase

According to the dynamical system (2.4), unless the equilibrium state is attained, the absolute mean fitness w(x) increases:
(2.7)Δw(x)=w(x′)-w(x)=∑kwkΔxk=∑kwkxk(wkw(x)-1)=∑kwk2xkw(x)-w(x)>0.
The mean fitness is maximal at equilibrium. The rate of increase of w(x) is(2.8)Δw(x)w(x)=∑kxk(wkw(x)-1)2=∑k(Δxk)2xk,
which is the variance in relative fitness σ¯2(x) defined in (2.3). These last two facts are sometimes termed the 1930s Fisher fundamental theorem of natural selection. The equilibria of (2.4) are the extremal states (0-faces) of the boundary of SK. To make it simple, if there is an allele whose fitness is strictly larger than the ones of the others, the deterministic evolutionary dynamics (2.4) will attain an equilibrium where only the fittest will survive; starting from any initial state of SK which is not an extremal (or monomorphic) point, the haploid trajectories will converge to this fittest state.


### 2.1.3. A Stochastic Interpretation of the Deterministic Dynamics (2.4)

A vectorx of SK can be thought of as a probability vector. The dynamical equation (2.4), as a nonlinear update mapping from SK to SK, may be viewed as the discrete-time nonlinear Master equation of some Markov process whose construction we now give. Suppose that we have a population of N haploid individuals each of which can be of one among K types or colors (carrying one among the K possible alleles). We will need to introduce an extra color-state, say ∂={0}, which will be absorbing for the process we will now construct. Let K(t):=Kn(t),n=1,…,N be the random color distribution of the population at time t, therefore with enlarged state-space {0,1,…,K}. Assume that the individuals are indistinguishable leading to the exchangeability property: Kn(t)=dK1(t),n=2,…,N (equality in distribution). Let Ut,m,t=1,2,…;m=1,…,N be N i.i.d. driving sequences of uniformly distributed random variables on [0,1], independent of K(t). Let w¯k:=wk/∑kwk,k=1,…,K. To decide the allele Km(t+1) carried by the individual number m∈{1,…,N} at time t+1, with 1(A) being the indicator function of the event A, consider the random Markovian dynamical system:(2.9)1(Km(t+1)=k,τm>t+1)=1(w¯kN∑n=1N1(Kn(t)=k,τn>t)>Ut+1,m).
Here k∈{1,…,K} and τn is the first time that Kn(t) hits the absorbing state ∂. As a result of Kn(t)=dK1(t), we naturally assume τn=dτ1,n=2,…,N. For each n, our model therefore attributes a positive probability that Kn(t)=0 for all t≥τn. Although, in principle, there is a possibility that the type of the mth particle is the one of the fictitious unobservable allele A0, as a result of (2.9), the sample paths of K(t) leading to this A0 are ruled out because focus is on the observable states only.In words, for the dynamics (2.9), the observable event Km(t+1)=k is realized (together then with τm>t+1) if the proportion at t of type-k individuals, weighted by the corresponding scaled fitness w¯k, is large enough (compared to Ut+1) and of course if the whole process was not absorbed at {0} in the previous step. Taking first the expectation with respect to the driving noise Ut+1,m in (2.9), we get
(2.10)ℙ(Km(t+1)=k,τm>t+1∣K(t)>0)=w¯kN∑n=1N1(Kn(t)=k,τn>t).
Putting zk(t):=ℙ(K1(t)=k,τ1>t), recalling (Kn(t);τn)=d(K1(t);τ1),n=2,…,N, and taking the expectation with respect to K(t), we get an unnormalized version of (2.4):
(2.11)zk(t+1)=w¯kzk(t),k∈{1,…,K}.
We have 1>ℙ(τ1>t)=∑k=1Kzk(t)=∑k=1Kxk(0)w¯kt→0, geometrically fast. Defining the normalized conditional probabilities
(2.12)xk(t)=zk(t)∑k=1Kzk(t)=ℙ(K1(t)=k∣τ1>t),
we obtain the normalized haploid dynamics (2.4):
(2.13)xk′=wkxk∑k=1Kwkxk,k∈{1,…,K}.
It now may be viewed as the nonlinear Master equation of some stochastic Markov process. Let us make some miscellaneous remarks.(i) Clearly this construction makes also sense if N=1 (a single particle). (ii) When N is finite, we should stress that the initial condition can be chosen to be deterministic, say with xk(0)=zk(0)=ik/N, for some sequence of integers ik≥0 satisfying i1+⋯+iK=N (i0=0) and quantifying the initial population sizes. It could also be chosen to be random, with xk(0) defining the initial probability distribution of the alleles. This occurs in the large N limit if ik=⌊Nxk(0)⌋ so that ik/N→xk(0). The latter choice may therefore be interpreted as a large N limit of the former one. (iii) In the stochastic interpretation (2.9) of the deterministic dynamics (2.4), xk(t) can be interpreted either as the probability that the random allele carried by a typical individual is Ak or like the expected proportion of the individuals of type k within the whole population (a frequentist point of view). (iv) The appeal to the coffin state ∂ was a necessary step to understand the normalization zk→xk.(v) Even though K(t) is exchangeable, it is not true that, with n1≠n2 any two distinct individuals, their random labels Kn1(t) and Kn2(t) are independent. The random algorithm allowing to update the joint types of Kn1(t) and Kn2(t) could be written down but is much more involved.


## 2.1.1. Dynamics

From the deterministic evolutionary genetics point of view, the discrete-time update of the allele frequency distribution on the simplexSK is given by (2.4)xk′=pk(x):=xkwkw(x),k=1,…,K.
(the symbol ′ is a common and useful notation to denote the updated frequency.) The quantity wk/w(x)-1 therefore interprets as the frequency-dependent Malthus growth rate parameter of xk. As required, the vector p(x):=pk(x),k=1,…,K, maps SK into SK. In vector form, with Dx:=diag(xk,k=1,…,K), the nonlinear deterministic evolutionary dynamics reads
(2.5)x′=p(x)=1w(x)Dwx=1w(x)Dxw,
or, with Δx:=x′-x, the increment of x(2.6)Δx=(1w(x)Dw-I)x.
Avoiding the trivial case where fitnesses are all equal, without loss of generality, we can assume that either w1≥⋯≥wK=1 or w1≤⋯≤wK=1. Thus that allele A1 or AK has largest fitness.


## 2.1.2. Mean Fitness Increase

According to the dynamical system (2.4), unless the equilibrium state is attained, the absolute mean fitness w(x) increases:
(2.7)Δw(x)=w(x′)-w(x)=∑kwkΔxk=∑kwkxk(wkw(x)-1)=∑kwk2xkw(x)-w(x)>0.
The mean fitness is maximal at equilibrium. The rate of increase of w(x) is(2.8)Δw(x)w(x)=∑kxk(wkw(x)-1)2=∑k(Δxk)2xk,
which is the variance in relative fitness σ¯2(x) defined in (2.3). These last two facts are sometimes termed the 1930s Fisher fundamental theorem of natural selection. The equilibria of (2.4) are the extremal states (0-faces) of the boundary of SK. To make it simple, if there is an allele whose fitness is strictly larger than the ones of the others, the deterministic evolutionary dynamics (2.4) will attain an equilibrium where only the fittest will survive; starting from any initial state of SK which is not an extremal (or monomorphic) point, the haploid trajectories will converge to this fittest state.


## 2.1.3. A Stochastic Interpretation of the Deterministic Dynamics (2.4)

A vectorx of SK can be thought of as a probability vector. The dynamical equation (2.4), as a nonlinear update mapping from SK to SK, may be viewed as the discrete-time nonlinear Master equation of some Markov process whose construction we now give. Suppose that we have a population of N haploid individuals each of which can be of one among K types or colors (carrying one among the K possible alleles). We will need to introduce an extra color-state, say ∂={0}, which will be absorbing for the process we will now construct. Let K(t):=Kn(t),n=1,…,N be the random color distribution of the population at time t, therefore with enlarged state-space {0,1,…,K}. Assume that the individuals are indistinguishable leading to the exchangeability property: Kn(t)=dK1(t),n=2,…,N (equality in distribution). Let Ut,m,t=1,2,…;m=1,…,N be N i.i.d. driving sequences of uniformly distributed random variables on [0,1], independent of K(t). Let w¯k:=wk/∑kwk,k=1,…,K. To decide the allele Km(t+1) carried by the individual number m∈{1,…,N} at time t+1, with 1(A) being the indicator function of the event A, consider the random Markovian dynamical system:(2.9)1(Km(t+1)=k,τm>t+1)=1(w¯kN∑n=1N1(Kn(t)=k,τn>t)>Ut+1,m).
Here k∈{1,…,K} and τn is the first time that Kn(t) hits the absorbing state ∂. As a result of Kn(t)=dK1(t), we naturally assume τn=dτ1,n=2,…,N. For each n, our model therefore attributes a positive probability that Kn(t)=0 for all t≥τn. Although, in principle, there is a possibility that the type of the mth particle is the one of the fictitious unobservable allele A0, as a result of (2.9), the sample paths of K(t) leading to this A0 are ruled out because focus is on the observable states only.In words, for the dynamics (2.9), the observable event Km(t+1)=k is realized (together then with τm>t+1) if the proportion at t of type-k individuals, weighted by the corresponding scaled fitness w¯k, is large enough (compared to Ut+1) and of course if the whole process was not absorbed at {0} in the previous step. Taking first the expectation with respect to the driving noise Ut+1,m in (2.9), we get
(2.10)ℙ(Km(t+1)=k,τm>t+1∣K(t)>0)=w¯kN∑n=1N1(Kn(t)=k,τn>t).
Putting zk(t):=ℙ(K1(t)=k,τ1>t), recalling (Kn(t);τn)=d(K1(t);τ1),n=2,…,N, and taking the expectation with respect to K(t), we get an unnormalized version of (2.4):
(2.11)zk(t+1)=w¯kzk(t),k∈{1,…,K}.
We have 1>ℙ(τ1>t)=∑k=1Kzk(t)=∑k=1Kxk(0)w¯kt→0, geometrically fast. Defining the normalized conditional probabilities
(2.12)xk(t)=zk(t)∑k=1Kzk(t)=ℙ(K1(t)=k∣τ1>t),
we obtain the normalized haploid dynamics (2.4):
(2.13)xk′=wkxk∑k=1Kwkxk,k∈{1,…,K}.
It now may be viewed as the nonlinear Master equation of some stochastic Markov process. Let us make some miscellaneous remarks.(i) Clearly this construction makes also sense if N=1 (a single particle). (ii) When N is finite, we should stress that the initial condition can be chosen to be deterministic, say with xk(0)=zk(0)=ik/N, for some sequence of integers ik≥0 satisfying i1+⋯+iK=N (i0=0) and quantifying the initial population sizes. It could also be chosen to be random, with xk(0) defining the initial probability distribution of the alleles. This occurs in the large N limit if ik=⌊Nxk(0)⌋ so that ik/N→xk(0). The latter choice may therefore be interpreted as a large N limit of the former one. (iii) In the stochastic interpretation (2.9) of the deterministic dynamics (2.4), xk(t) can be interpreted either as the probability that the random allele carried by a typical individual is Ak or like the expected proportion of the individuals of type k within the whole population (a frequentist point of view). (iv) The appeal to the coffin state ∂ was a necessary step to understand the normalization zk→xk.(v) Even though K(t) is exchangeable, it is not true that, with n1≠n2 any two distinct individuals, their random labels Kn1(t) and Kn2(t) are independent. The random algorithm allowing to update the joint types of Kn1(t) and Kn2(t) could be written down but is much more involved.


## 2.2. Single Locus: Diploid Population withK Alleles

We now run into similar considerations but with diploid populations.


### 2.2.1. Joint Evolutionary Dynamics

Letwk,l≥0,k,l=1,…,K stand for the absolute fitness of the genotypes AkAl attached to a single locus. Assume wk,l=wl,k (wk,l being proportional to the probability of an AkAl surviving to maturity, it is natural to take wk,l=wl,k). Let then W be the symmetric fitness matrix with k,l-entry wk,l.Assume that the current frequency distribution at timet of the genotypes AkAl is given by xk,l. Let X be the frequencies array with k,l-entry xk,l. The joint evolutionary dynamics in the diploid case is given by the updating:(2.14)xk,l′=xk,lwk,lω(X),
where the mean fitness ω is now given by ω(X)=∑k,lxk,lwk,l. The relative fitness of the genotype AkAl is wk,l/ω(X). The joint dynamics takes the matrix form:
(2.15)X′=1ω(X)X∘W=1ω(X)W∘X,
where ∘ stands for the (commutative) Hadamard product of matrices.LetJ be the flat K×K matrix whose entries are all 1. Then
(2.16)ΔX:=X′-X=1ω(X)(X-J)∘W=1ω(X)W∘(X-J).
We will also let(2.17)σ2(X)=∑k,l=1K‍xk,l(wk,l-ω(X))2
stand for the genotypic variance in absolute fitness and(2.18)σ¯2(X)=∑k,l=1Kxk,l(ωk,lω(X)-1)2=σ2(X)ω(X)2
will stand for the diploid variance in relative fitness.Consider the problem of evaluating the increase of the mean fitness. We have(2.19)Δω(X)=∑k,lΔxk,lwk,l=∑k,lxk,l(wk,l2ω(X)-wk,l)=ω(X)σ¯2(X)>0
with a relative rate of increase: Δw(X)/w(X)=σ¯2(X). This is the full diploid version of the Fisher theorem.


### 2.2.2. Marginal Allelic Dynamics

Assuming a Hardy-Weinberg equilibrium, the frequency distribution at timet, say xk,l, of the genotypes AkAl is given by xk,l=xkxl where xk=∑lxk,l is the marginal frequency of allele Ak in the whole genotypic population . The whole frequency information is now enclosed within x:=xk,k=1,…,K. For instance, the mean fitness is now given by the quadratic form: ω(x):=∑k,lxkxlwk,l=x*Wx, with x* being the transposed line vector of the column vector x=X1 (1 the unit K-vector). We will also let(2.20)σ2(x)=∑k,l=1Kxkxl(wk,l-ω(x))2
stand for the genotypic variance in absolute fitness and(2.21)σ¯2(x)=∑k,l=1Kxkxl(wk,lω(x)-1)2=σ2(x)ω(x)2
will stand for the diploid variance in relative fitness.Consider now the update of the allelic marginal frequenciesx themselves. If we first define the frequency-dependent marginal fitness of Ak by wk(x)=(Wx)k:=∑lwk,lxl, the marginal dynamics is given as in (2.4) by:(2.22)xk′=xkwk(x)ω(x)=1ω(x)xk(Wx)k=:pk(x),k=1,…,K.
This dynamics involves a multiplicative interaction between xk and (Wx)k, the kth entry of the image Wx of x by W. In (2.22) there is a normalization by the quadratic form ω(x)=x*Wx. In vector form (2.22) reads
(2.23)x′=1ω(x)DxWx=1ω(x)DWxx=:p(x),
where p maps SK into SK. Iterating, the time-t frequency distribution is(2.24)x(t)=p(p(⋯ttimes⋯(p(x0)))).Except for the fact that the mean fitnessω in (2.22) is now a quadratic form in x and that the marginal fitness of Ak is now frequency-dependent, depending linearly on x, as far as the marginal frequencies are concerned, the updating formalism (2.22) in the diploid case looks very similar to the one in (2.4) describing the haploid case.In the diploid case, assuming fitnesses to bemultiplicative, say with wk,l=wkwl, then wk(x)/x*Wx=wk/∑lwlxl and the dynamics (2.22) boils down to (2.4). However, the mean fitness in this case is ω(x)=(∑lwlxl)2 and not w(x)=∑lwlxl as in the haploid case.


### 2.2.3. A Stochastic Interpretation of the Deterministic Dynamics (2.22)

As for the haploid case, there is a Markov chain governed by the Master equation (2.22). Consider a population of diploid individuals. The number of alleles N in this population is therefore twice the number of genes. Each allele can be of one among K types or colors (carrying one among the K possible alleles). As before, we introduce an extra color-state, say ∂={0} which is absorbing for the process to be constructed. For c=1,2, let Kc(t):=Knc(t),n=1,…,N be two independent copies of the random color distribution of the allelic population at time t. Let K(t)=(K1(t),K2(t)). Assume that the alleles are indistinguishable within each sample, leading to Knc(t)=dK1c(t),n=2,…,N,c=1,2. For c=1,2, let Ut,mc,t=1,2,…;m=1,…,N be two mutually independent i.i.d. driving N-sequences of uniformly distributed random variables on [0,1] and independent of K(t). To decide the type of the random allele Km(t+1),m=1,…,N, at time t+1, consider now the Markovian dynamical system:(2.25)1(Km(t+1)=k,τm>t+1)=1(1N∑n=1N1(Kn1(t)=k,τn1>t)>Ut+1,m1)1(1N∑n=1NW¯k,Kn2(t)1(τn2>t)>Ut+1,m2).
Here k∈{1,…,K},τnc are the first hitting times of each Knc(t) of the absorbing state ∂,c=1,2 and W¯:=W/∥W∥ for any matrix norm ∥W∥, say, for example, ∥W∥=∑k,lwk,l. We assume τnc=dτ1c,n=2,…,N,c=1,2.In words, for this new dynamics, the observable eventKm(t+1)=k is seen to be realized together with τm>t+1 if two natural independent conditions are now satisfied which can be read from the two indicator functions in the right-hand side of (2.25):(i)
first, the proportion att of type-k individuals of the first copy should be large enough (compared to Ut+11);(ii)
second, for the second sample copyK2(t), the expected fitness of the genotypes AkAl,l=1,…,K, containing allele Ak at t should be large enough (compared to Ut+12).As for the haploid case, it is necessary that both processesK(t) should not be absorbed at {0} in the previous step. Taking first the expectation with respect to the independent driving noises Ut+1,mc in (2.25), we get
(2.26)ℙ(Km(t+1)=k,τm>t+1∣K(t)>0)=1N∑n=1N1(Kn1(t)=k,τn1>t)·1N∑n=1NW¯k,Kn2(t)1(τn2>t).
Putting zk(t)=ℙ(K1(t)=k,τ1>t), taking the expectation with respect to K(t) and using our independence and exchangeability hypotheses, we get
(2.27)zk(t+1)=zk(t)·𝔼[W¯k,K12(t)1(τ12>t)]=zk(t)·∑l=1kW¯k,lzl(t),
corresponding to an unnormalized version of (2.22):
(2.28)zk(t+1)=zk(t)(W¯z(t))k,k∈{1,…,K}.
Defining the normalized conditional probabilities
(2.29)xk(t)=zk(t)∑k=1Kzk(t)=ℙ(K1(t)=k∣τ1>t),
we obtain the normalized nonlinear Markovian diploid dynamics (2.22):
(2.30)xk′=xk(Wx)k∑k=1Kxk(Wx)k,k∈{1,…,K}.
Note that this construction makes sense if N=2 (a single individual and its 2 alleles of K possible types). The need for two copies of K(t) was governed by the quadratic character of the interaction appearing in the numerator of (2.22).


### 2.2.4. Increase of Mean Fitness

Again, the mean fitnessω(x), as a Lyapunov function, increases as time passes by. We indeed have
(2.31)Δω(x)=ω(x′)-ω(x)=1ω(x)2∑k,lxkwk(x)wk,lxlwl(x)-∑k,lxkwk,lxl>0,
because, defining 0<X(x):=∑k,lxk(1-wk(x)/ω(x))wk,l(1-wl(x)/ω(x))xl, we have
(2.32)Δω(x)=X(x)+2ω(x)(∑kxkwk(x)2-ω(x)2)>0.
Its partial rate of increase due to frequency shifts only is δω(x):=∑kΔxkwk(x). It satisfies(2.33)δω(x)ω(x)=∑kxk(wk(x)ω(x)-1)2=∑k(Δxk)2xk=σ¯A2(x)2,
where σ¯A2(x) is the allelic variance in relative fitness:(2.34)σ¯A2(x):=2∑k=1Kxk(wk(x)ω(x)-1)2.


### 2.2.5. An Alternative Representation of the Allelic Dynamics

There is an alternative vectorial representation of the dynamics (2.22). Define the symmetric positive-definite matrix G(x)=Dx(I-1x*) with quadratic entries in the frequencies:
(2.35)G(x)k,l=xk(δk,l-xl).
Introduce the quantity VW(x)=(1/2)logω(x), which is half the logarithm of the mean fitness. Then, (2.22) may be recast as the gradient-like dynamics:(2.36)Δx=1ω(x)G(x)Wx=G(x)∇VW(x),
with |Δx|=1*Δx=0 as a result of 1*G(x)=0*. Note
(2.37)∇VW(x)*Δx=∇VW(x)*G(x)∇VW(x)≥0.
The dynamics (2.36) is of pure gradient-type with respect to the Svirezhev-Shashahani distance metric dG(x,x′); see [5, 6]. For this metric, the distance between x and x′=x+Δx of SK is
(2.38)dG(x,x′)=(Δx*G-1Δx)1/2=(∑k=1Kxk-1(Δxk)2)1/2.
From (2.33) and (2.34), this quantity, which is the length of Δx satisfying |Δx|=0, is also the square-root of half the allelic variance (the standard deviation) in relative fitness.


## 2.2.1. Joint Evolutionary Dynamics

Letwk,l≥0,k,l=1,…,K stand for the absolute fitness of the genotypes AkAl attached to a single locus. Assume wk,l=wl,k (wk,l being proportional to the probability of an AkAl surviving to maturity, it is natural to take wk,l=wl,k). Let then W be the symmetric fitness matrix with k,l-entry wk,l.Assume that the current frequency distribution at timet of the genotypes AkAl is given by xk,l. Let X be the frequencies array with k,l-entry xk,l. The joint evolutionary dynamics in the diploid case is given by the updating:(2.14)xk,l′=xk,lwk,lω(X),
where the mean fitness ω is now given by ω(X)=∑k,lxk,lwk,l. The relative fitness of the genotype AkAl is wk,l/ω(X). The joint dynamics takes the matrix form:
(2.15)X′=1ω(X)X∘W=1ω(X)W∘X,
where ∘ stands for the (commutative) Hadamard product of matrices.LetJ be the flat K×K matrix whose entries are all 1. Then
(2.16)ΔX:=X′-X=1ω(X)(X-J)∘W=1ω(X)W∘(X-J).
We will also let(2.17)σ2(X)=∑k,l=1K‍xk,l(wk,l-ω(X))2
stand for the genotypic variance in absolute fitness and(2.18)σ¯2(X)=∑k,l=1Kxk,l(ωk,lω(X)-1)2=σ2(X)ω(X)2
will stand for the diploid variance in relative fitness.Consider the problem of evaluating the increase of the mean fitness. We have(2.19)Δω(X)=∑k,lΔxk,lwk,l=∑k,lxk,l(wk,l2ω(X)-wk,l)=ω(X)σ¯2(X)>0
with a relative rate of increase: Δw(X)/w(X)=σ¯2(X). This is the full diploid version of the Fisher theorem.


## 2.2.2. Marginal Allelic Dynamics

Assuming a Hardy-Weinberg equilibrium, the frequency distribution at timet, say xk,l, of the genotypes AkAl is given by xk,l=xkxl where xk=∑lxk,l is the marginal frequency of allele Ak in the whole genotypic population . The whole frequency information is now enclosed within x:=xk,k=1,…,K. For instance, the mean fitness is now given by the quadratic form: ω(x):=∑k,lxkxlwk,l=x*Wx, with x* being the transposed line vector of the column vector x=X1 (1 the unit K-vector). We will also let(2.20)σ2(x)=∑k,l=1Kxkxl(wk,l-ω(x))2
stand for the genotypic variance in absolute fitness and(2.21)σ¯2(x)=∑k,l=1Kxkxl(wk,lω(x)-1)2=σ2(x)ω(x)2
will stand for the diploid variance in relative fitness.Consider now the update of the allelic marginal frequenciesx themselves. If we first define the frequency-dependent marginal fitness of Ak by wk(x)=(Wx)k:=∑lwk,lxl, the marginal dynamics is given as in (2.4) by:(2.22)xk′=xkwk(x)ω(x)=1ω(x)xk(Wx)k=:pk(x),k=1,…,K.
This dynamics involves a multiplicative interaction between xk and (Wx)k, the kth entry of the image Wx of x by W. In (2.22) there is a normalization by the quadratic form ω(x)=x*Wx. In vector form (2.22) reads
(2.23)x′=1ω(x)DxWx=1ω(x)DWxx=:p(x),
where p maps SK into SK. Iterating, the time-t frequency distribution is(2.24)x(t)=p(p(⋯ttimes⋯(p(x0)))).Except for the fact that the mean fitnessω in (2.22) is now a quadratic form in x and that the marginal fitness of Ak is now frequency-dependent, depending linearly on x, as far as the marginal frequencies are concerned, the updating formalism (2.22) in the diploid case looks very similar to the one in (2.4) describing the haploid case.In the diploid case, assuming fitnesses to bemultiplicative, say with wk,l=wkwl, then wk(x)/x*Wx=wk/∑lwlxl and the dynamics (2.22) boils down to (2.4). However, the mean fitness in this case is ω(x)=(∑lwlxl)2 and not w(x)=∑lwlxl as in the haploid case.


## 2.2.3. A Stochastic Interpretation of the Deterministic Dynamics (2.22)

As for the haploid case, there is a Markov chain governed by the Master equation (2.22). Consider a population of diploid individuals. The number of alleles N in this population is therefore twice the number of genes. Each allele can be of one among K types or colors (carrying one among the K possible alleles). As before, we introduce an extra color-state, say ∂={0} which is absorbing for the process to be constructed. For c=1,2, let Kc(t):=Knc(t),n=1,…,N be two independent copies of the random color distribution of the allelic population at time t. Let K(t)=(K1(t),K2(t)). Assume that the alleles are indistinguishable within each sample, leading to Knc(t)=dK1c(t),n=2,…,N,c=1,2. For c=1,2, let Ut,mc,t=1,2,…;m=1,…,N be two mutually independent i.i.d. driving N-sequences of uniformly distributed random variables on [0,1] and independent of K(t). To decide the type of the random allele Km(t+1),m=1,…,N, at time t+1, consider now the Markovian dynamical system:(2.25)1(Km(t+1)=k,τm>t+1)=1(1N∑n=1N1(Kn1(t)=k,τn1>t)>Ut+1,m1)1(1N∑n=1NW¯k,Kn2(t)1(τn2>t)>Ut+1,m2).
Here k∈{1,…,K},τnc are the first hitting times of each Knc(t) of the absorbing state ∂,c=1,2 and W¯:=W/∥W∥ for any matrix norm ∥W∥, say, for example, ∥W∥=∑k,lwk,l. We assume τnc=dτ1c,n=2,…,N,c=1,2.In words, for this new dynamics, the observable eventKm(t+1)=k is seen to be realized together with τm>t+1 if two natural independent conditions are now satisfied which can be read from the two indicator functions in the right-hand side of (2.25):(i)
first, the proportion att of type-k individuals of the first copy should be large enough (compared to Ut+11);(ii)
second, for the second sample copyK2(t), the expected fitness of the genotypes AkAl,l=1,…,K, containing allele Ak at t should be large enough (compared to Ut+12).As for the haploid case, it is necessary that both processesK(t) should not be absorbed at {0} in the previous step. Taking first the expectation with respect to the independent driving noises Ut+1,mc in (2.25), we get
(2.26)ℙ(Km(t+1)=k,τm>t+1∣K(t)>0)=1N∑n=1N1(Kn1(t)=k,τn1>t)·1N∑n=1NW¯k,Kn2(t)1(τn2>t).
Putting zk(t)=ℙ(K1(t)=k,τ1>t), taking the expectation with respect to K(t) and using our independence and exchangeability hypotheses, we get
(2.27)zk(t+1)=zk(t)·𝔼[W¯k,K12(t)1(τ12>t)]=zk(t)·∑l=1kW¯k,lzl(t),
corresponding to an unnormalized version of (2.22):
(2.28)zk(t+1)=zk(t)(W¯z(t))k,k∈{1,…,K}.
Defining the normalized conditional probabilities
(2.29)xk(t)=zk(t)∑k=1Kzk(t)=ℙ(K1(t)=k∣τ1>t),
we obtain the normalized nonlinear Markovian diploid dynamics (2.22):
(2.30)xk′=xk(Wx)k∑k=1Kxk(Wx)k,k∈{1,…,K}.
Note that this construction makes sense if N=2 (a single individual and its 2 alleles of K possible types). The need for two copies of K(t) was governed by the quadratic character of the interaction appearing in the numerator of (2.22).


## 2.2.4. Increase of Mean Fitness

Again, the mean fitnessω(x), as a Lyapunov function, increases as time passes by. We indeed have
(2.31)Δω(x)=ω(x′)-ω(x)=1ω(x)2∑k,lxkwk(x)wk,lxlwl(x)-∑k,lxkwk,lxl>0,
because, defining 0<X(x):=∑k,lxk(1-wk(x)/ω(x))wk,l(1-wl(x)/ω(x))xl, we have
(2.32)Δω(x)=X(x)+2ω(x)(∑kxkwk(x)2-ω(x)2)>0.
Its partial rate of increase due to frequency shifts only is δω(x):=∑kΔxkwk(x). It satisfies(2.33)δω(x)ω(x)=∑kxk(wk(x)ω(x)-1)2=∑k(Δxk)2xk=σ¯A2(x)2,
where σ¯A2(x) is the allelic variance in relative fitness:(2.34)σ¯A2(x):=2∑k=1Kxk(wk(x)ω(x)-1)2.


## 2.2.5. An Alternative Representation of the Allelic Dynamics

There is an alternative vectorial representation of the dynamics (2.22). Define the symmetric positive-definite matrix G(x)=Dx(I-1x*) with quadratic entries in the frequencies:
(2.35)G(x)k,l=xk(δk,l-xl).
Introduce the quantity VW(x)=(1/2)logω(x), which is half the logarithm of the mean fitness. Then, (2.22) may be recast as the gradient-like dynamics:(2.36)Δx=1ω(x)G(x)Wx=G(x)∇VW(x),
with |Δx|=1*Δx=0 as a result of 1*G(x)=0*. Note
(2.37)∇VW(x)*Δx=∇VW(x)*G(x)∇VW(x)≥0.
The dynamics (2.36) is of pure gradient-type with respect to the Svirezhev-Shashahani distance metric dG(x,x′); see [5, 6]. For this metric, the distance between x and x′=x+Δx of SK is
(2.38)dG(x,x′)=(Δx*G-1Δx)1/2=(∑k=1Kxk-1(Δxk)2)1/2.
From (2.33) and (2.34), this quantity, which is the length of Δx satisfying |Δx|=0, is also the square-root of half the allelic variance (the standard deviation) in relative fitness.


## 2.3. Equilibria (Diploid Case)

The mean fitness increase phenomenon occurs till the evolutionary dynamics reaches an equilibrium state. We wish to briefly discuss the questions relative to equilibria in the diploid case.


### 2.3.1. Preliminaries

In contrast with the haploid case, in the diploid situation, the dynamics (2.22) can have more complex equilibrium points, satisfying wk(xeq)=w1(xeq),k=2,…,K, and ∑lxeq,l=1. To avoid linear manifolds of equilibria, we first assume that all principal minors of W are nonsingular and also that the fitnesses of all homozygotes wk,k are positive. In this case, from the Bézout theorem, the number of equilibria is finite and less or equal than the number 2K-1 of faces of SK. Note that the K extremal endpoints of SK (0-faces) are always monomorphic fixed points of (2.22).An instructive example fulfilling these preliminary conditions isW=I. There are 2K-1 equilibrium points (the barycenters of the (Kk+1)k-dimensional faces, k=0,…,K-1), but only one polymorphic equilibrium which is the barycenter xB of SK. This point is the one with minimal fitness and it is unstable. The 0-faces are stable fixed points whereas the barycenters of the k faces with k∈{1,…,K-2} are saddle-points. The simplex SK could be partitioned into pieces each of which being the attraction basins of the stable 0-face states: in contrast with the haploid case, the type of the survivor is not necessarily the one of the fittest; it will depend on the initial condition.Similar conclusions can be drawn if instead ofW=I, we start with W=(I-Dλ)-1 where λ:=(λk,k=1,…,K) satisfies 0≤λ<1 (meaning 0≤λk<1, for all k). In this case again, there is only one unstable polymorphic equilibrium which is easily seen to be xeq=(1-λ)/(K-|λ|)∈SK.Due to its evolutionary interest, we would like now to discuss the opportunity for a polymorphic state to be stable. Under the above assumptions onW, a unique stable internal (polymorphic) equilibrium state can exist, necessary and sufficient conditions being that (i) there is a unique z>0 for which Wz=1 and (ii)W has exactly one strictly positive dominant eigenvalue and at least one strictly negative eigenvalue or else the sequence of principal minors of W alternates in sign (see Kingman [7]). If this is the case, the equilibrium polymorphic state is xeq=z/|z|. It is stable in the sense that it is a local maximum of the mean fitness ω(x)=x*Wx, with ω(xeq)=1/|z|. Since |W|≠0, the linearization of p(x) at xeq has no eigenvalue of modulus 1 and so xeq is hyperbolic and/or isolated (see [8]). A stable isolated polymorphic state is asymptotically Lyapunov stable.Under these additional assumptions therefore onW, starting from any initial condition in the interior of SK, all trajectories will be attracted by xeq=z/|z| which is an isolated global maximum of ω(x).When there is no such unique globally stable polymorphic equilibrium, all trajectories will still converge but perhaps to a local equilibrium state where some alleles get extinct. Which allele and how many alleles are concerned seems to be an unsolved problem in its full generality.


### 2.3.2. Special Classes of Fitness Matrices Leading to a Polymorphic State

We now draw the attention on a particular class of fitness matrices that lead to a polymorphic state, either unstable or stable. We start with the unstable case, extending the above special diagonal caseW=(I-Dλ)-1 leading to a unique unstable polymorphic state.(i) The Unstable Case
LetΛ≥0 be a symmetric irreducible strictly substochastic matrix satisfying Λ1:=q<1: the positive mass-defect vector of Λ is 1-q. Let λ>0. Define the symmetric strictly potential matrix: W=λ-1(I-Λ)-1≥0, with W-1=λ(I-Λ) defining a strictly row-diagonally dominant Stieltjes matrix with the properties [3]: (W-1)k,k>0,(W-1)k,l≤0 for k≠l and (W-1)k,k+∑l≠k(W-1)k,l>0 for all k. Then
(2.39)λ1>z=W-11=λ(I-Λ)1=λ(1-q)>0.
The vector z=λ(1-q) is called the equilibrium potential of W. We have |z|=λ(K-|q|). For this class of W therefore, Wz=1 admits a positive solution z.
Conversely, given a nonsingular matrixW≥0 satisfying Wz=1 for some z≥0, the matrix Λ=I-λ-1W-1 defines a substochastic matrix if and only W-1 satisfies (W-1)k,k>0,(W-1)k,l≤0 for k≠l and λ≥maxk(W-1)k,k. Then, W-1 is row-diagonally dominant and W=λ-1(I-Λ)-1 is a potential matrix.Strictly ultrametric (sUm) matrices are special classes of positive-definite and symmetric strictly potential matrices [2]. An sUm matrix W is symmetric with nonnegative entries, satisfying (i) wk,l≥min{wk,j,wj,l}, for all j,k,l and (ii) wk,k>maxl≠k{wk,l}, for all k (If in condition (ii),≥ is substituted for >, then W is simply an ultrametric matrix and this new condition is implied by (i)). If W is an sUm matrix, the fitness dynamics will admit an unstable polymorphic equilibrium state, as a result of W being positive-definite.Remark 2.1.
Suppose thatΛ is substochastic and primitive. Then W=λ-1(I-Λ)-1>0 is an ultrametric matrix. If V is the Hadamard reciprocal of W with entries vk,l=1/wk,l, it satisfies that vk,l≤max{vk,j,vj,l}. Therefore V is an ultrametric distance associated to the ultrametric potential W. Tree matrices are ultrametric matrices that are not sUm.(ii) The Stable Case
To produce a stable equilibrium state from the sUm matrix construction, letW=λ-1(I-Λ)-1 define a symmetric strictly potential matrix as before. Then, there is a z=λ(1-q)>0 for which Wz=1. With α>1, define
(2.40)W̃:=α|z|J-W.
With W̃=[w̃k,l], we have minw̃k,l=α/|z|-maxwk,l and we can choose α>1 so that W̃≥0. We have (α-1)W̃z=1 and δ*W̃δ=-δ*Wδ<0 for all δ≠0 satisfying |δ|=0 showing that xeq:=z/|z|=(1-q)/(K-|q|) is now a stable polymorphic state for W̃. If W is an sUm matrix, then clearly W̃ satisfies the “anti-sUm” property expressing a fitness domination of the heterozygotes AkAl over the homozygotes:
(2.41)w̃k,l≤max{w̃k,j,w̃j,l},∀j,k,l,w̃k,k<minl≠k{w̃k,l},∀k.
Nonnegative symmetric negative-definite anti-sUm fitness matrices W̃ will therefore display a stable polymorphic equilibrium state xeq. Note that xeq*W̃xeq=(α-1)/|z| is now the maximal value of the mean fitness.Example 2.2.
WhenK=2, with s>-1, h>0, and sh>-1, let
(2.42)W=[1+s1+sh1+sh1]>0
define the fitness matrix with selection parameter s and dominance h. This W is sUm if and only if s<0 and h>1. The equilibrium state is xeq*:=(h/(2h-1);(h-1)/(2h-1)) and it is unstable. This W is anti-sUm if and only if s>0 and h>1. The equilibrium state is the same but it is now stable.Note that a singular multiplicative fitness matrix of the formW=ww* cannot be a strictly potential matrix because its determinant |W| is zero.


### 2.3.3. A General Construction to Produce sUm and Anti-sUm Matrices

Consider the problem consisting in splitting binarily and recursively the set{1,…,K} till complete reduction to singletons (leaves) which are left behind in the process. For instance, consider the refinement sequence with K=6 of {1,…,6}≡(123456): (123456)→((23)(1456))→(((2)(3))((16)(45)))→(((1)(6))((4)(5))).Starting from the left, there are2K-1=11 blocks of symbols (the total number of nodes in the splitting binary tree with K leaves). To each encountered block, numbered from l=1 to 2K-1, starting from the left, attach a vector ul of size K with ith entry ul(i)=1 if symbol i is in the block string,0 otherwise. For instance, from the above sequence, u1*=(1,1,1,1,1,1),u2*=(0,1,1,0,0,0),u3*=(1,0,0,1,1,1),…, and u10*=(0,0,0,1,0,0),u11*=(0,0,0,0,1,0). To each such ul, attach a number sl which is >0 if |ul|=1 (the leaves) and ≥0 if |ul|>1 (the internal nodes, including the root). Then (see [3]), for any choice of sl respecting these constraints
(2.43)W=∑l=12K-1slulul*≥0
is an sUm matrix and any sUm matrix can be represented in this way. Because for the ul corresponding to the leaves sl>0, the diagonal terms of W are necessarily >0.Since for each set{1,…,K}, there are bK splitting tree sequences where bK satisfies bK=∑k=1K-1(Kk)bkbK-k,k≥2,b1=1, there are many ways to generate an sUm matrix.Clearly, for each splitting procedure, withλ-1:=∑l=12K-1sl>0,W may be written as (2.44)W=λ-1(J-∑l=22K-1‍λsl(J-ulul*))≥0,
where the matrices J-ulul* take values in {0,1}.Now, withγ-1>0, any matrix of the form (2.45)W̃=(λ-1+γ-1)J-W=γ-1(J+∑l=22K-1γsl(J-ulul*))≥0
is an anti-sUm matrix. Assuming γ=1 and sl=shl,W̃ may be written under the form: W̃=J+sA where A:=∑l=22K-1hl(J-ulul*). It has at most 2K-1 independent parameters, namely, the hl,l=2,…,2K-1, and s. If the hl are known, then W̃ is a one-parameter family of fitness matrices. (Fitness matrices of the form J+sA were considered in [9] in the context of the estimation of s problem).Examples 2.3.
Assume also thatsl=s/(2(K-1)), for all l, with s>0 being a selection parameter. Then, defining the (0,1]-valued matrix
(2.46)A:=12(K-1)∑l=22K-1(J-ulul*)>0,
an anti-sUm matrix of the form W̃=J+sA will admit a stable polymorphic equilibrium. Clearly, A itself is anti-sUm. Because of this, there is a zA>0 such that AzA=1. Thus (J+sA)zA=(|zA|+s)1 showing that, with z=zA/(|zA|+s),W̃z=(J+sA)z=1. We thus have |z|=|zA|/(|zA|+s) and so xeq=z/|z|=zA/|zA|. Furthermore, the equilibrium mean fitness for such models is xeq*W̃xeq=1/|z|=(|zA|+s)/|zA|>1.
For the following simple sequence example withK=4 (four alleles, say A,C,T,G), (1234)→((1)(234))→((24)(3))→((2)(4)) we find
(2.47)A=16[5666635465456453]
which itself clearly is a symmetric anti-sUm matrix, together with W̃=J+sA. For this example, xeq*=(1/13)(8,1,3,1) and the equilibrium mean fitness is 1+s/13.
Note that takingsl=s/(2(K-1)) just for the indices l that were initially chosen to satisfy sl>0 would also lead to anti-sUm matrices A>0 and W̃=J+sA>0. In this case, the sum (2.46) defining A should then be restricted to the indices l from 2 to 2K-1 for which sl>0. Proceeding in this extreme way for the above simple example, we get the borderline anti-sUm shapes:
(2.48)A=16[3444434444344443],W̃=[s+223+2s33+2s33+2s33+2s3s+223+2s33+2s33+2s33+2s3s+223+2s33+2s33+2s33+2s3s+22].
Here, xeq* is simply the barycenter of the 4-simplex.Remark 2.4.
The ultrametric conditions (2.41) should be compared with the so-called “triangle inequality” conditions (leading to stable polymorphism) pointed out in [10], which read
(2.49)w̃k,l<w̃k,j+w̃l,j,∀k≠landatleastonej≠k,j≠l,w̃k,l>w̃k,k+w̃l,l2,∀k≠l.
It is clear that the class of anti-sUm matrices is a particular subclass of the Lewontin one.


## 2.3.1. Preliminaries

In contrast with the haploid case, in the diploid situation, the dynamics (2.22) can have more complex equilibrium points, satisfying wk(xeq)=w1(xeq),k=2,…,K, and ∑lxeq,l=1. To avoid linear manifolds of equilibria, we first assume that all principal minors of W are nonsingular and also that the fitnesses of all homozygotes wk,k are positive. In this case, from the Bézout theorem, the number of equilibria is finite and less or equal than the number 2K-1 of faces of SK. Note that the K extremal endpoints of SK (0-faces) are always monomorphic fixed points of (2.22).An instructive example fulfilling these preliminary conditions isW=I. There are 2K-1 equilibrium points (the barycenters of the (Kk+1)k-dimensional faces, k=0,…,K-1), but only one polymorphic equilibrium which is the barycenter xB of SK. This point is the one with minimal fitness and it is unstable. The 0-faces are stable fixed points whereas the barycenters of the k faces with k∈{1,…,K-2} are saddle-points. The simplex SK could be partitioned into pieces each of which being the attraction basins of the stable 0-face states: in contrast with the haploid case, the type of the survivor is not necessarily the one of the fittest; it will depend on the initial condition.Similar conclusions can be drawn if instead ofW=I, we start with W=(I-Dλ)-1 where λ:=(λk,k=1,…,K) satisfies 0≤λ<1 (meaning 0≤λk<1, for all k). In this case again, there is only one unstable polymorphic equilibrium which is easily seen to be xeq=(1-λ)/(K-|λ|)∈SK.Due to its evolutionary interest, we would like now to discuss the opportunity for a polymorphic state to be stable. Under the above assumptions onW, a unique stable internal (polymorphic) equilibrium state can exist, necessary and sufficient conditions being that (i) there is a unique z>0 for which Wz=1 and (ii)W has exactly one strictly positive dominant eigenvalue and at least one strictly negative eigenvalue or else the sequence of principal minors of W alternates in sign (see Kingman [7]). If this is the case, the equilibrium polymorphic state is xeq=z/|z|. It is stable in the sense that it is a local maximum of the mean fitness ω(x)=x*Wx, with ω(xeq)=1/|z|. Since |W|≠0, the linearization of p(x) at xeq has no eigenvalue of modulus 1 and so xeq is hyperbolic and/or isolated (see [8]). A stable isolated polymorphic state is asymptotically Lyapunov stable.Under these additional assumptions therefore onW, starting from any initial condition in the interior of SK, all trajectories will be attracted by xeq=z/|z| which is an isolated global maximum of ω(x).When there is no such unique globally stable polymorphic equilibrium, all trajectories will still converge but perhaps to a local equilibrium state where some alleles get extinct. Which allele and how many alleles are concerned seems to be an unsolved problem in its full generality.


## 2.3.2. Special Classes of Fitness Matrices Leading to a Polymorphic State

We now draw the attention on a particular class of fitness matrices that lead to a polymorphic state, either unstable or stable. We start with the unstable case, extending the above special diagonal caseW=(I-Dλ)-1 leading to a unique unstable polymorphic state.(i) The Unstable Case
LetΛ≥0 be a symmetric irreducible strictly substochastic matrix satisfying Λ1:=q<1: the positive mass-defect vector of Λ is 1-q. Let λ>0. Define the symmetric strictly potential matrix: W=λ-1(I-Λ)-1≥0, with W-1=λ(I-Λ) defining a strictly row-diagonally dominant Stieltjes matrix with the properties [3]: (W-1)k,k>0,(W-1)k,l≤0 for k≠l and (W-1)k,k+∑l≠k(W-1)k,l>0 for all k. Then
(2.39)λ1>z=W-11=λ(I-Λ)1=λ(1-q)>0.
The vector z=λ(1-q) is called the equilibrium potential of W. We have |z|=λ(K-|q|). For this class of W therefore, Wz=1 admits a positive solution z.
Conversely, given a nonsingular matrixW≥0 satisfying Wz=1 for some z≥0, the matrix Λ=I-λ-1W-1 defines a substochastic matrix if and only W-1 satisfies (W-1)k,k>0,(W-1)k,l≤0 for k≠l and λ≥maxk(W-1)k,k. Then, W-1 is row-diagonally dominant and W=λ-1(I-Λ)-1 is a potential matrix.Strictly ultrametric (sUm) matrices are special classes of positive-definite and symmetric strictly potential matrices [2]. An sUm matrix W is symmetric with nonnegative entries, satisfying (i) wk,l≥min{wk,j,wj,l}, for all j,k,l and (ii) wk,k>maxl≠k{wk,l}, for all k (If in condition (ii),≥ is substituted for >, then W is simply an ultrametric matrix and this new condition is implied by (i)). If W is an sUm matrix, the fitness dynamics will admit an unstable polymorphic equilibrium state, as a result of W being positive-definite.Remark 2.1.
Suppose thatΛ is substochastic and primitive. Then W=λ-1(I-Λ)-1>0 is an ultrametric matrix. If V is the Hadamard reciprocal of W with entries vk,l=1/wk,l, it satisfies that vk,l≤max{vk,j,vj,l}. Therefore V is an ultrametric distance associated to the ultrametric potential W. Tree matrices are ultrametric matrices that are not sUm.(ii) The Stable Case
To produce a stable equilibrium state from the sUm matrix construction, letW=λ-1(I-Λ)-1 define a symmetric strictly potential matrix as before. Then, there is a z=λ(1-q)>0 for which Wz=1. With α>1, define
(2.40)W̃:=α|z|J-W.
With W̃=[w̃k,l], we have minw̃k,l=α/|z|-maxwk,l and we can choose α>1 so that W̃≥0. We have (α-1)W̃z=1 and δ*W̃δ=-δ*Wδ<0 for all δ≠0 satisfying |δ|=0 showing that xeq:=z/|z|=(1-q)/(K-|q|) is now a stable polymorphic state for W̃. If W is an sUm matrix, then clearly W̃ satisfies the “anti-sUm” property expressing a fitness domination of the heterozygotes AkAl over the homozygotes:
(2.41)w̃k,l≤max{w̃k,j,w̃j,l},∀j,k,l,w̃k,k<minl≠k{w̃k,l},∀k.
Nonnegative symmetric negative-definite anti-sUm fitness matrices W̃ will therefore display a stable polymorphic equilibrium state xeq. Note that xeq*W̃xeq=(α-1)/|z| is now the maximal value of the mean fitness.Example 2.2.
WhenK=2, with s>-1, h>0, and sh>-1, let
(2.42)W=[1+s1+sh1+sh1]>0
define the fitness matrix with selection parameter s and dominance h. This W is sUm if and only if s<0 and h>1. The equilibrium state is xeq*:=(h/(2h-1);(h-1)/(2h-1)) and it is unstable. This W is anti-sUm if and only if s>0 and h>1. The equilibrium state is the same but it is now stable.Note that a singular multiplicative fitness matrix of the formW=ww* cannot be a strictly potential matrix because its determinant |W| is zero.


## 2.3.3. A General Construction to Produce sUm and Anti-sUm Matrices

Consider the problem consisting in splitting binarily and recursively the set{1,…,K} till complete reduction to singletons (leaves) which are left behind in the process. For instance, consider the refinement sequence with K=6 of {1,…,6}≡(123456): (123456)→((23)(1456))→(((2)(3))((16)(45)))→(((1)(6))((4)(5))).Starting from the left, there are2K-1=11 blocks of symbols (the total number of nodes in the splitting binary tree with K leaves). To each encountered block, numbered from l=1 to 2K-1, starting from the left, attach a vector ul of size K with ith entry ul(i)=1 if symbol i is in the block string,0 otherwise. For instance, from the above sequence, u1*=(1,1,1,1,1,1),u2*=(0,1,1,0,0,0),u3*=(1,0,0,1,1,1),…, and u10*=(0,0,0,1,0,0),u11*=(0,0,0,0,1,0). To each such ul, attach a number sl which is >0 if |ul|=1 (the leaves) and ≥0 if |ul|>1 (the internal nodes, including the root). Then (see [3]), for any choice of sl respecting these constraints
(2.43)W=∑l=12K-1slulul*≥0
is an sUm matrix and any sUm matrix can be represented in this way. Because for the ul corresponding to the leaves sl>0, the diagonal terms of W are necessarily >0.Since for each set{1,…,K}, there are bK splitting tree sequences where bK satisfies bK=∑k=1K-1(Kk)bkbK-k,k≥2,b1=1, there are many ways to generate an sUm matrix.Clearly, for each splitting procedure, withλ-1:=∑l=12K-1sl>0,W may be written as (2.44)W=λ-1(J-∑l=22K-1‍λsl(J-ulul*))≥0,
where the matrices J-ulul* take values in {0,1}.Now, withγ-1>0, any matrix of the form (2.45)W̃=(λ-1+γ-1)J-W=γ-1(J+∑l=22K-1γsl(J-ulul*))≥0
is an anti-sUm matrix. Assuming γ=1 and sl=shl,W̃ may be written under the form: W̃=J+sA where A:=∑l=22K-1hl(J-ulul*). It has at most 2K-1 independent parameters, namely, the hl,l=2,…,2K-1, and s. If the hl are known, then W̃ is a one-parameter family of fitness matrices. (Fitness matrices of the form J+sA were considered in [9] in the context of the estimation of s problem).Examples 2.3.
Assume also thatsl=s/(2(K-1)), for all l, with s>0 being a selection parameter. Then, defining the (0,1]-valued matrix
(2.46)A:=12(K-1)∑l=22K-1(J-ulul*)>0,
an anti-sUm matrix of the form W̃=J+sA will admit a stable polymorphic equilibrium. Clearly, A itself is anti-sUm. Because of this, there is a zA>0 such that AzA=1. Thus (J+sA)zA=(|zA|+s)1 showing that, with z=zA/(|zA|+s),W̃z=(J+sA)z=1. We thus have |z|=|zA|/(|zA|+s) and so xeq=z/|z|=zA/|zA|. Furthermore, the equilibrium mean fitness for such models is xeq*W̃xeq=1/|z|=(|zA|+s)/|zA|>1.
For the following simple sequence example withK=4 (four alleles, say A,C,T,G), (1234)→((1)(234))→((24)(3))→((2)(4)) we find
(2.47)A=16[5666635465456453]
which itself clearly is a symmetric anti-sUm matrix, together with W̃=J+sA. For this example, xeq*=(1/13)(8,1,3,1) and the equilibrium mean fitness is 1+s/13.
Note that takingsl=s/(2(K-1)) just for the indices l that were initially chosen to satisfy sl>0 would also lead to anti-sUm matrices A>0 and W̃=J+sA>0. In this case, the sum (2.46) defining A should then be restricted to the indices l from 2 to 2K-1 for which sl>0. Proceeding in this extreme way for the above simple example, we get the borderline anti-sUm shapes:
(2.48)A=16[3444434444344443],W̃=[s+223+2s33+2s33+2s33+2s3s+223+2s33+2s33+2s33+2s3s+223+2s33+2s33+2s33+2s3s+22].
Here, xeq* is simply the barycenter of the 4-simplex.Remark 2.4.
The ultrametric conditions (2.41) should be compared with the so-called “triangle inequality” conditions (leading to stable polymorphism) pointed out in [10], which read
(2.49)w̃k,l<w̃k,j+w̃l,j,∀k≠landatleastonej≠k,j≠l,w̃k,l>w̃k,k+w̃l,l2,∀k≠l.
It is clear that the class of anti-sUm matrices is a particular subclass of the Lewontin one.


## 3. Stochastic Evolutionary Dynamics

We now switch to the random point of view of multiallelic evolutionary dynamics driven by selection. There are two models of interest: the Wright-Fisher and the Moran models. We start with Wright-Fisher.


### 3.1. The Wright-Fisher Model

The Wright-Fisher model is a discrete space-time model which takes into account another important driving source of evolution, namely, the genetic drift whose nature is exclusively random.


#### 3.1.1. The Model and Its First Properties

Consider an allelic population with constant sizeN. In the haploid (diploid) case, N is (twice) the number of real individuals. Let i:=ik and i′:=ik′,k=1,…,K be two vectors of integers quantifying the size of the allelic populations at two consecutive generations t and t+1. We will let SK,N={iintegers:|i|=∑k=1Kik=N}. Suppose that the stochastic evolutionary dynamics is now given by a Markov chain whose one-step transition matrix P from states I=i to I′=i′ is given by the multinomial Wright-Fisher model:(3.1)ℙ(It+1′=i′∣It=i)=:P(i,i′)=(Ni1′⋯iK′)∏k=1Kpk(iN)ik′.
Therefore, given It=i, to form the next generation, each allele chooses its type independently of the others with probability p, where p=(pk,k=1,…,K) is given either by (2.4) in the haploid case or by (2.22) in the diploid case. In the diploid case, the mechanism p is assumed to present a unique polymorphic state, either stable or unstable.The state-space dimension of the Markov chain governed by (3.1) is n=(N+K-1K-1) (the number of compositions of integer N into K nonnegative parts which is also the number of ways to assign N indistinguishable balls into K distinguishable boxes). To view P(i,i′) as a standard transition matrix of some Markov chain, we need first to order the states i and i′ in (3.1). Starting from the bottom right corner of P states should be arranged in decreasing order when listing the entries of P moving up and left along the lines and columns, respectively; or equivalently, starting from the top left corner of P states should be arranged in increasing order when moving down and right.For example, we can order the statesi∈SK,N in decreasing order from n to 1, as follows. Let H(i)=(1/N)∑k=1K(N+1)-(k-1)ik be a base-(N+1) code of the state i that will serve as a ranking function. The largest state in for which H(i) is maximal (equal to 1) is in*:=(N,0,…,0). Given a state i, define the subsequent state in decreasing order as
(3.2)σ-(i)=argminj:H(j)<H(i)(H(i)-H(j)).
Then in*=(N,0,…,0),in-1*=(N-1,1,0,…,0),…,i1*=(0,…,0,N) and to pass from state i to the next state σ-(i) in this decreasing sequence, there is a unique δi with entries in ℤ satisfying |δi|=0 and such that σ-(i)=i-δi. This way to order the n states is consistent with the reverse lexicographic order. For instance, if N=3,K=4, there are (63)=20 states ordered in decreasing order as follows:3000>2100>2010>2001>1200>1110>1101>1020>1011>1002>0300>0210>0201>0120>0111>0102>0030>0021>0012>0003.The way the digits are propagated from left to right is clear and the consecutiveδi can easily be obtained. Proceeding in this way to order states, P(i,i′) is a well-defined conventional object (matrix).From (3.1), the marginal transition matrix from i to Ik′=ik′ is binomial bin(N,pk(i/N)) with:
(3.3)P(i,ik′)=(Nik′)pk(iN)(1-pk(iN))N-ik′.
Given I=i, the kth component Ik′ of the updated state is random with
(3.4)𝔼i(Ik′N)=pk(iN),σi2(Ik′N)=pk(i/N)(1-pk(i/N))N,
suggesting that, in the large N population limit, the deterministic evolutionary dynamics should be recovered. Indeed, from (3.1), the Laplace-Stieltjes transform of the joint law of It+1′∣It=i reads
(3.5)𝔼i(e-∑k=1KλkIk′/N)=(∑k=1Kpk(iN)e-λk/N)N~N↑∞e-∑k=1Kλkpk(i/N).
From the strong law of large numbers therefore, if ik:=⌊Nxk⌋,k=1,…,K, then, given x=xk,k=1,…,K,
(3.6)Ik′N→N↑∞a.s.xk′=pk(x)
which is (2.22): when the population under study is very large, random fluctuations as modelled by (3.1) can be ignored so that the gene frequencies evolve deterministically.


#### 3.1.2. Fixation Probabilities

Letel be the K-null vector except for its lth entry which is 1. The extremal pure states SK,Nex:={ilex:=Nel,l=1,…,K} are all absorbing for this Markov chain because pk(ilex/N)=δk,l and, from (3.1), any additional fixed point which p could have on the boundary-faces of SK which are not points does not give rise to an absorbing state for P. Under our assumptions on p, the chain is not recurrent, rather it is transient. Depending on the initial condition, say i0, the chain will necessarily end up in one of the extremal states ilex, with some fixation probability, say πl(i0), which can be computed as follows. Let πl:=πl(i),i∈SK,N be the harmonic function of the Wright-Fisher Markov chain which is the smallest solution to the boundary problem:(3.7)(I-P)πl=0ifi∈SK,N∖SK,Nex,πl=1(i=ilex)ifi∈SK,Nex.
We also have
(3.8)ℙ(Iτ=Nel∣I0=i0)=πl(i0),
where τ (<∞ almost surely) is the random hitting time of SK,Nex for It and the πl(i0)s are normalized so as ∑lπl(i0)=1. Thus πl(i0) are the searched probabilities to end up in state ilex starting from state i0.From (3.7), πl(i0) is known if i0∈SK,Nex. The remaining unknown restriction, say πlQ, of πl=πl(i0) to the nonextremal states is easily seen to be(3.9)πlQ=(I-Q)-1pilex,i0∈SK,N∖SK,Nex.
In (3.9), Q is obtained from P after erasing the lines and columns corresponding to all the K extremal states and pilex is the ilex-column of P where the entries corresponding to the extremal states have been deleted. When dealing with the Wright-Fisher model, Q is a positive matrix and also pilex>0, therefore πl(i0)>0 for all i0∈SK,N∖SK,Nex and this for each l: starting from any state i0 which is not an extremal state, there is a positive probability to hit any of the extremal states. Fixation of the state ilex means extinction of the remaining monomorphic states. It would therefore be of interest to understand the structure of the set 𝒜l={i0:πl(i0)>∑k≠lπk(i0)}={i0:πl(i0)>1/2}, for each l, which is the stochastic version of the attraction basin of ilex, especially when l is the label of the extremal state with largest fitness wk,k. If i0∈𝒜l indeed, the probability to end up in ilex is larger than the probability to end up in any other extremal state.Unfortunately, the development of the inverse ofI-Q in terms of its adjugate matrix in (3.9) shows that these fixation probabilities have a very complex determinantal alternating structure and the question of identifying 𝒜l is very complex.Withπl being the solution (3.9) to the Dirichlet problem with boundary conditions (3.7), the equilibrium measure of the chain therefore is
(3.10)πeq(i0):=∑l=1Kπl(i0)δilex,
which depends on i0. Necessarily, one allele will fixate and there is no polymorphic equilibrium state even when dealing with diploid populations. Which allele and with what probability will depend on the initial condition. Thanks to fluctuations, the picture therefore looks very different from the one pertaining to the deterministic theory.Of importance also is the time it takes to get extinct. It relies on similar techniques. For instance, the expected overall fixation timeα(i0):=𝔼i0(τ) solves the boundary problem:
(3.11)(I-P)α=1,i0∈SK,N∖SK,Nex,α=0,i0∈SK,Nex,
where α:=α(i0),i0∈SK,N. The restriction αQ of α to the nonextremal states therefore is (3.12)αQ=(I-Q)-11,i0∈SK,N∖SK,Nex.


#### 3.1.3. ConditioningIt on Nonfixation

There are four places where questions relative to conditioning on fixation are relevant in this context.  (Similar conditioning problems were considered in [11] in the context of the 2-alleles Wright-Fisher diffusion).(i) Consider the full fixation vector πl:=πl(i0). Remove from πl the states i0 for which πl(i0)=0 only, keeping the one, ilex, for which πl(ilex)=1. The size of this vector, say πlR, is nR=(N+K-1K-1)-K+1. Let R be the corresponding nR×nR-reduced transition matrix obtained from P after erasing the lines and columns corresponding to all the K monomorphic states except ilex. The Markov chain It conditioned to exit in the extremal state ilex only admits the stochastic transition matrix:
(3.13)Rl:=DπlR-1RDπlR.
It is obtained from R after a diagonal Doob transform based on πlR. The chain governed by Rl admits a unique absorbing state which is ilex. The entries of Rl are
(3.14)Rl(i,i′)=πlR(i′)πlR(i)R(i,i′),
and for this new conditioned Markov chain, transitions to states i′ for which πlR(i′)>πlR(i) are favored.(ii) Let us now consider again the fully reduced transition matrix Q obtained from P after erasing the lines and columns corresponding to all the K monomorphic states. The matrix Q is substochastic and irreducible, with Q1<1. The law of the process corresponding to It conditioned on avoiding the monomorphic states before t evolves as follows: with τ:=⋀l=1Kτl being the time needed for first hitting one of the extremal states for It, let πt(i)=ℙ(It=i∣τ>t). Then, with πt=πt(i),i∈SK,N∖SK,Nex,
(3.15)πt+1*=πt*Qπt*Q1
is the nonlinear Master Equation governing its evolution [12]. The reduced state-space dimension of this Markov chain is nQ=(N+K-1K-1)-K and πt→t↑∞π∞ where π∞ is the left Perron probability eigenvector of Q associated to the dominant Perron eigenvalue ρQ<1, namely, ρQπ∞*=π∞*Q.  (π∞ is called a Yaglom limit (see [13]) or a quasistationary distribution). If the process is started using this limiting quasistationary distribution, it remains in the same state over time and the fixation time τ is geometrically distributed with success probability ρQ.(iii) One can define another stochastic process I¯t which admits the stochastic transition matrix: Q¯:=DQ1-1Q again defined on the reduced state-space. For each t, we have
(3.16)Q¯(i,i′)=ℙ(It+1=i′∣It=i,τ>1),
and the conditioning on nonfixation occurs at each transition time. This process is an ergodic Markov chain with invariant probability measure solving π¯eq*=π¯eq*Q¯. It has the following closed-form determinantal expression (see [14, Section  6] and [15, page 1559]):(3.17)π¯eq(i)=|(I-Q¯)[i,i]|∑i|(I-Q¯)[i,i]|,i∈SK,N∖SK,Nex,
where (I-Q¯)[i,i] is the submatrix resulting from the deletion of row i and column i of I-Q¯. The question as to whether the process governed by Q¯ is reversible or not arises. Defining the transition matrix of the time-reversed process Q⟵ by
(3.18)Q⟵*=Dπ¯eqQ¯Dπ¯eq-1,
it does not hold that Q⟵=Q¯ and so detailed balance does not hold. Indeed, Q¯ is similar to the transition matrix of an ergodic Wright-Fisher model and Wright-Fisher chains are not reversible.(iv) If we condition on nonfixation in the remote future (see [16] for additional details), we get a Markov chain whose stochastic transition matrix is
(3.19)Q̃=ρQ-1Dψ∞-1QDψ∞.
Here ψ∞ is the positive right Perron eigenvector of Q associated to the Perron eigenvalue ρQ<1 satisfying ρQψ∞=Qψ∞. This vector can be chosen so that ∑kπ∞,kψ∞,k=1, where π∞ is again the left Perron probability eigenvector of Q associated to ρQ<1 (see [17]). With π̃t(·)=lims↑∞ℙ(It=·∣τ>t+s), we have π̃t+1*=π̃t*Q̃. The process governed by Q̃ is an ergodic Markov chain whose invariant probability measure is π̃∞=π∞∘ψ∞, the Schur product of π∞ and ψ∞ with kth entry π̃∞,k=π∞,kψ∞,k.For the last three conditionings, it is difficult to extract some information on the limiting distribution, eitherπ∞ or π¯eq or π̃∞, respectively. This is because it would suppose to solve the eigenvalue problems explicitly which is out of reach, at least theoretically. However, assuming a diploid population with a polymorphic equilibrium state xeq for p, we expect that these distributions will present a global (local) maximum (minimum) near xeq if xeq is stable (unstable). These limiting distributions should be more sharply peaked around the extremum if we consider the conditioning (iv) compared to (ii) because, the latter conditioning being more stringent than the former, it should charge more heavily the sample paths that stay away from the monomorphic states.Finally, we would like to stress that all these considerations are also relevant in the context of another fundamental stochastic model arising in the context of evolutionary genetics. We will give some elements of how to proceed with this model presenting very different properties.


### 3.2. TheK-Alleles Moran Model

We now focus on the Moran model.


#### 3.2.1. The Multiallelic Moran Model

Letα,β∈{1,…,K}. In the Moran version of the stochastic evolution,  given It=I=i, the only accessible values of I′ are the neighboring states: iα,β′:=i+dα,β where dα,β*:=(0,…,0,-1,0,…,1,0,…,0). Here -1 is in position α and 1 in position β≠α corresponding to the transfer of an individual from the box α (if nonempty) to the box β. With n(i)=#{k:ik>0} being the number of nonempty entries of i, there are n(i)(K-1)≤K(K-1) accessible states from i. The Moran stochastic evolutionary dynamics is now given by a Markov chain whose one-step transition matrix P from states I=i to I′=i′ is(3.20)ℙ(It+1=i′∣It=i)=0ifi′≠iα,β′,ℙ(It+1=iα,β′∣It=i)=:P(i,iα,β′)=iαNpβ(iN),
where p=(pβ(i/N),β=1,…,K) is given either by (2.4) in the haploid case or by (2.22) in the diploid case.SummingP(i,iα,β′) over α,β, β≠α in (3.20), we get the holding probability
(3.21)ℙ(It+1=i∣It=i)=1-∑α,β:β≠αiαNpβ(iN)=∑αiαNpα(iN),
completing the characterization of the K-alleles Moran model. Under our assumptions on p, the holding probabilities are equal to 1 only for the extremal states i∈SK,Nex which are therefore the only absorbing states of the Moran chain, just like for the Wright-Fisher model. The drift at state i is(3.22)𝔼(It+1-It∣It=i)=∑αiαN∑β≠αpβ(iN)dα,β.Let us compute the Laplace-Stieltjes Transform ofI′ in the context of a Moran model. Omitting the argument i/N in pβ(i/N), we get the factorized form: (3.23)𝔼i(e-∑kλkIk′)=∑α,β:α≠βe-∑kλkiα,β′(k)P(i,iα,β′)+e-∑kλkik∑βiβNpβ=e-∑kλkik(∑α,β:α≠βe-∑kλkdα,β(k)P(i,iα,β′)+∑βiβNpβ)=e-∑kλkik(∑α,β:α≠βe-(λβ-λα)iαNpβ+∑βiβNpβ)=e-∑kλkik(∑βe-λβpβ∑α≠βiαNeλα+∑βiβNpβ)=e-∑kλkik(∑βe-λβpβ(∑αiαNeλα-iβNeλβ)+∑βiβNpβ)=(e-∑kλkik)(∑αiαNeλα)(∑βe-λβpβ).
Putting λl=0 if l≠k, the kth marginal reads(3.24)𝔼i(e-λkIk′)=e-λkik(1-ikN+eλkikN)(1-pk+e-λkpk)
showing that Ik(t) is of the random walk type. Indeed, we get ℙi(Ik′=ik′)=0 if ik′≠ik±1 or ik′≠ik and (3.25)ℙi(Ik′=ik)=(1-ikN)(1-pk)+ikNpk,ℙi(Ik′=ik+1)=(1-ikN)pk=∑l≠kP(i,il,k′),ℙi(Ik′=ik-1)=ikN(1-pk)=∑l≠kP(i,ik,l′).
We have (3.26)𝔼i(Ik′N)=ikN+1N(pk-ikN);σi2(Ik′N)=1N2(ikN(1-ikN)+pk(1-pk)).
There is only a small correction (of order N-1) of the updated mean to its current value and fluctuations around the mean are small too (of order N-1). The evolution process is very slow.


#### 3.2.2. Fixation Probabilities

As a random walk model, the Moran model has a much simpler transition matrixP of the Jacobi type. The equilibrium measure of the chain again is(3.27)πeq:=∑l=1Kπl(i0)δilex,
where πl again solves the Dirichlet problem (3.7) but with this new simpler Jacobi P. For the Moran model, the explicit expression (3.9) of the fixation probability simplifies a little bit because pilex(i)≠0 only for the K-1 neighboring states of ilex that is {i:i+dα,β=ilex} for some dα,β.WhenK=2 (2 alleles), the random walk transition probabilities (p1+p2=1) (3.28)ℙ(It+1=i1,2′∣It=i)=:P(i,i+(-11))=i1Np2(iN),ℙ(It+1=i2,1′∣It=i)=:P(i,i+(1-1))=i2Np1(iN)
are the probabilities that the first component It,1 of It moves down and up by one unit, respectively. In this case, the Dirichlet problem giving the fixation probabilities solves explicitly. With ϕ(i0)=1+∑i=1i0-1∏i1=1i(i1p2(i/N)/i2p1(i/N)) (ϕ(0)=0) being the harmonic function of the 2-alleles chain, we easily get that (3.29)π1(i0,N-i0)=ϕ(i0)ϕ(N)
is the probability that the extremal state i1ex=(N,0) is reached given i0=(i0,N-i0). Assuming a model with multiplicative fitnesses: pα(i/N)=(iα/N)(wα/ω(i/N)), then ϕ takes the simple form (i1+i2=N) (3.30)ϕ(i0)=1+∑i=1i0-1∏i1=1ii1p2(i/N)i2p1(i/N)=1+∑i=1i0-1(w2w1)i
showing (see [1, page 109]) that(3.31)π1(i0,N-i0)=1-(w2/w1)i01-(w2/w1)N.
Assuming w1=1+s/N and w2=1, putting i0=[Nx0], for large N, we get [18] (3.32)π1(Nx0,N-i0)~1-e-sx01-e-s.In the general fitness case(3.33)ϕ(i0)=1+∑i=1i0-1‍∏i1=1ii1p2(i/N)i2p1(i/N)=1+∑i=1i0-1‍∏i1=1i(Wi/N)2(Wi/N)1,
where, as usual, (Wi/N)k=∑l=12wk,lil/N,k,l=1,2, leading to (3.34)(WiN)1=w1,2+(w1,1-w1,2)i1N,(WiN)2=w2,2+(w2,1-w2,2)i1N.This 2-alleles exact solution can be used in the full diploidK-alleles Moran case with multiplicative fitnesses. Indeed, from this, with i0=(i1,…,iK), the fixation probability of Al can be conjectured to be approximated qualitatively by(3.35)πl(ilex)=1,πl(i0)=1-((∑k≠likwk/(N-il))/wl)il1-((∑k≠likwk/(N-il))/wl)Nifi0≠ilex.
This can be justified as follows: mark one particular box with size il, corresponding to the allele Al with fitness wl. Then clump the K-1 remaining boxes into a single box with size N-il, corresponding to a fictitious allele with average fitness ∑k≠likwk/(N-il). We are left with a 2-alleles Moran multiplicative fitness model for which (from the 2-alleles exact solution) the fixation probability of Al is given by (3.35). This formula constitutes sort of a mean field approximation to the full Dirichlet problem associated to the Moran model.Assumingwk~1+sk/N, a Kimura-like approximation of (3.35) would lead for large N to (3.36)πl(Nx0)~1-e(xl/(1-xl))∑k≠lxksk1-e(1/(1-xl))∑k≠lxksk,
where x0=(x1,…,xK) is now a point of the continuous K-simplex different from the l-unit vector el:=(0,…,0,1,0,…,0).From (3.35), when the fitness of allele Al is large (small), compared to the average fitness of the remaining alleles, then the fixation probability of Al gets close to 1 (resp., to 0). Note also that the larger il is, the larger the fixation probability is.As required also, for allk≠l,πl(ikex)=(1-(wk/wl)0)/(1-(wk/wl)N)=0. As another particular initial configuration case, suppose that we start from the 2-alleles type state: i0=:i0(k,l)=(0,…,0,N-1,0,…,0,1,0,…,0) where the 1 is in position l (i.e., il=1) and the entry N-1 in position k≠l (i.e., ik=N-1). Although for this choice of the initial state, the fixation of Ak is very likely, there still is a positive probability that allele Al gets fixed which is seen to be from (3.35): (3.37)πl(i0(k,l))=1-wk/wl1-(wk/wl)N,
depending only on the relative fitnesses of Al and Ak (see [19, 20] for a similar expression). As N gets large, this probability gets close to 1-wk/wl if wk<wl and close to (wl/wk)N-1~0 if wk>wl.Conversely, assuming the 2-alleles type statei0 to be defined by il=N-1 and ik=1, (3.38)πl(i0)=1-(wk/wl)N-11-(wk/wl)N,
which, as required, gets close to 1 as N gets large if wk<wl and close to wl/wk if wk>wl.


#### 3.2.3. ConditioningIt on Nonfixation

The conditioning developments discussed for the Wright-Fisher model are also relevant in the Moran model context substituting theP of Moran for the P of Wright-Fisher. Let us revisit condition (iii).(iii) With Q being now the reduced substochastic matrix of the full Moran transition matrix defined in (3.20), consider the stochastic process I¯t with stochastic transition matrix: Q¯:=DQ1-1Q defined on the reduced state-space with dimension nQ. The process I¯t is again an ergodic Markov chain with invariant probability measure solving π¯eq*=π¯eq*Q¯. With Q⟵ being the transition matrix of the time-reversed process it now holds that Q⟵=Q¯ and so detailed balance holds when dealing with the Moran case (see also [21]). This will be proved if we can exhibit an equilibrium probability measure π¯eq such that(3.39)π¯eq(i)Q¯(i,i′)=Q¯(i′,i)π¯eq(i′),
for all neighboring state i′=i+di of i.Withj=(j1,…,jK) any terminal state, suppose that we want to use (3.39) to compute π¯eq(j) starting from the smallest available state in the system which is j0=(0,…,0,1,N-1). This is possible because π¯eq(j) may be represented as(3.40)π¯eq(j)=π¯eq(j0)∏i=j0j-dK,1Q¯(i,i′)Q¯(i′,i),
where π¯eq(j0) can be chosen so that ∑jπ¯eq(j)=1. Let us give some details.(i) Note first, by reversing path, that for two consecutive states(i,i′), the ratio Q¯(i,i′)/Q¯(i′,i) can be computed. We have(3.41)Q¯(i,i′)Q¯(i′,i)=(Q1)i′(Q1)iQ(i,i′)Q(i′,i),
where i′=i+di for some di of the form dk,l and therefore(3.42)Q(i,i′)=ikNpl(iN),Q(i′,i)=Q(i′,i′+dl,k)=il+1Npk(i+dk,lN).
Clearly, from such a structure of the entries of Q¯, for each possible transition i→i′, the ratio Q¯(i,i′)/Q¯(i′,i) will only depend separately on a ratio involving the terminal and the initial states i′ and i (the detailed balance condition holds).(ii) Second, the sequence ofi-s in (3.40) is governed by the following path starting from j0 and ending up in the target state j: we can use the following sequence of dk,ls:(3.43)(dK,K-1⋯jK-1-1dK,K-1)(dK,K-2⋯jK-2dK,K-1)⋯(dK,1⋯j1dK,1),
filling up successively the entries of j to the left of the last entry of j0 by using the N-1 individuals of the reservoir state j0=(0,…,0,1,N-1). By doing so, each intermediate state i is separated from the next i′ by some clearly identified di, and after evaluating the probability ratio Q¯(i,i′)/Q¯(i′,i) for each consecutive states of this sequence, we are done. Because there exists a probability distribution π¯eq(j) such that the reversibility identity (3.39) holds, then this Moran process is reversible with (3.40) as its stationary distribution. Using (3.42), this constitutes an exact explicit product-form formula.


## 3.1. The Wright-Fisher Model

The Wright-Fisher model is a discrete space-time model which takes into account another important driving source of evolution, namely, the genetic drift whose nature is exclusively random.


### 3.1.1. The Model and Its First Properties

Consider an allelic population with constant sizeN. In the haploid (diploid) case, N is (twice) the number of real individuals. Let i:=ik and i′:=ik′,k=1,…,K be two vectors of integers quantifying the size of the allelic populations at two consecutive generations t and t+1. We will let SK,N={iintegers:|i|=∑k=1Kik=N}. Suppose that the stochastic evolutionary dynamics is now given by a Markov chain whose one-step transition matrix P from states I=i to I′=i′ is given by the multinomial Wright-Fisher model:(3.1)ℙ(It+1′=i′∣It=i)=:P(i,i′)=(Ni1′⋯iK′)∏k=1Kpk(iN)ik′.
Therefore, given It=i, to form the next generation, each allele chooses its type independently of the others with probability p, where p=(pk,k=1,…,K) is given either by (2.4) in the haploid case or by (2.22) in the diploid case. In the diploid case, the mechanism p is assumed to present a unique polymorphic state, either stable or unstable.The state-space dimension of the Markov chain governed by (3.1) is n=(N+K-1K-1) (the number of compositions of integer N into K nonnegative parts which is also the number of ways to assign N indistinguishable balls into K distinguishable boxes). To view P(i,i′) as a standard transition matrix of some Markov chain, we need first to order the states i and i′ in (3.1). Starting from the bottom right corner of P states should be arranged in decreasing order when listing the entries of P moving up and left along the lines and columns, respectively; or equivalently, starting from the top left corner of P states should be arranged in increasing order when moving down and right.For example, we can order the statesi∈SK,N in decreasing order from n to 1, as follows. Let H(i)=(1/N)∑k=1K(N+1)-(k-1)ik be a base-(N+1) code of the state i that will serve as a ranking function. The largest state in for which H(i) is maximal (equal to 1) is in*:=(N,0,…,0). Given a state i, define the subsequent state in decreasing order as
(3.2)σ-(i)=argminj:H(j)<H(i)(H(i)-H(j)).
Then in*=(N,0,…,0),in-1*=(N-1,1,0,…,0),…,i1*=(0,…,0,N) and to pass from state i to the next state σ-(i) in this decreasing sequence, there is a unique δi with entries in ℤ satisfying |δi|=0 and such that σ-(i)=i-δi. This way to order the n states is consistent with the reverse lexicographic order. For instance, if N=3,K=4, there are (63)=20 states ordered in decreasing order as follows:3000>2100>2010>2001>1200>1110>1101>1020>1011>1002>0300>0210>0201>0120>0111>0102>0030>0021>0012>0003.The way the digits are propagated from left to right is clear and the consecutiveδi can easily be obtained. Proceeding in this way to order states, P(i,i′) is a well-defined conventional object (matrix).From (3.1), the marginal transition matrix from i to Ik′=ik′ is binomial bin(N,pk(i/N)) with:
(3.3)P(i,ik′)=(Nik′)pk(iN)(1-pk(iN))N-ik′.
Given I=i, the kth component Ik′ of the updated state is random with
(3.4)𝔼i(Ik′N)=pk(iN),σi2(Ik′N)=pk(i/N)(1-pk(i/N))N,
suggesting that, in the large N population limit, the deterministic evolutionary dynamics should be recovered. Indeed, from (3.1), the Laplace-Stieltjes transform of the joint law of It+1′∣It=i reads
(3.5)𝔼i(e-∑k=1KλkIk′/N)=(∑k=1Kpk(iN)e-λk/N)N~N↑∞e-∑k=1Kλkpk(i/N).
From the strong law of large numbers therefore, if ik:=⌊Nxk⌋,k=1,…,K, then, given x=xk,k=1,…,K,
(3.6)Ik′N→N↑∞a.s.xk′=pk(x)
which is (2.22): when the population under study is very large, random fluctuations as modelled by (3.1) can be ignored so that the gene frequencies evolve deterministically.


### 3.1.2. Fixation Probabilities

Letel be the K-null vector except for its lth entry which is 1. The extremal pure states SK,Nex:={ilex:=Nel,l=1,…,K} are all absorbing for this Markov chain because pk(ilex/N)=δk,l and, from (3.1), any additional fixed point which p could have on the boundary-faces of SK which are not points does not give rise to an absorbing state for P. Under our assumptions on p, the chain is not recurrent, rather it is transient. Depending on the initial condition, say i0, the chain will necessarily end up in one of the extremal states ilex, with some fixation probability, say πl(i0), which can be computed as follows. Let πl:=πl(i),i∈SK,N be the harmonic function of the Wright-Fisher Markov chain which is the smallest solution to the boundary problem:(3.7)(I-P)πl=0ifi∈SK,N∖SK,Nex,πl=1(i=ilex)ifi∈SK,Nex.
We also have
(3.8)ℙ(Iτ=Nel∣I0=i0)=πl(i0),
where τ (<∞ almost surely) is the random hitting time of SK,Nex for It and the πl(i0)s are normalized so as ∑lπl(i0)=1. Thus πl(i0) are the searched probabilities to end up in state ilex starting from state i0.From (3.7), πl(i0) is known if i0∈SK,Nex. The remaining unknown restriction, say πlQ, of πl=πl(i0) to the nonextremal states is easily seen to be(3.9)πlQ=(I-Q)-1pilex,i0∈SK,N∖SK,Nex.
In (3.9), Q is obtained from P after erasing the lines and columns corresponding to all the K extremal states and pilex is the ilex-column of P where the entries corresponding to the extremal states have been deleted. When dealing with the Wright-Fisher model, Q is a positive matrix and also pilex>0, therefore πl(i0)>0 for all i0∈SK,N∖SK,Nex and this for each l: starting from any state i0 which is not an extremal state, there is a positive probability to hit any of the extremal states. Fixation of the state ilex means extinction of the remaining monomorphic states. It would therefore be of interest to understand the structure of the set 𝒜l={i0:πl(i0)>∑k≠lπk(i0)}={i0:πl(i0)>1/2}, for each l, which is the stochastic version of the attraction basin of ilex, especially when l is the label of the extremal state with largest fitness wk,k. If i0∈𝒜l indeed, the probability to end up in ilex is larger than the probability to end up in any other extremal state.Unfortunately, the development of the inverse ofI-Q in terms of its adjugate matrix in (3.9) shows that these fixation probabilities have a very complex determinantal alternating structure and the question of identifying 𝒜l is very complex.Withπl being the solution (3.9) to the Dirichlet problem with boundary conditions (3.7), the equilibrium measure of the chain therefore is
(3.10)πeq(i0):=∑l=1Kπl(i0)δilex,
which depends on i0. Necessarily, one allele will fixate and there is no polymorphic equilibrium state even when dealing with diploid populations. Which allele and with what probability will depend on the initial condition. Thanks to fluctuations, the picture therefore looks very different from the one pertaining to the deterministic theory.Of importance also is the time it takes to get extinct. It relies on similar techniques. For instance, the expected overall fixation timeα(i0):=𝔼i0(τ) solves the boundary problem:
(3.11)(I-P)α=1,i0∈SK,N∖SK,Nex,α=0,i0∈SK,Nex,
where α:=α(i0),i0∈SK,N. The restriction αQ of α to the nonextremal states therefore is (3.12)αQ=(I-Q)-11,i0∈SK,N∖SK,Nex.


### 3.1.3. ConditioningIt on Nonfixation

There are four places where questions relative to conditioning on fixation are relevant in this context.  (Similar conditioning problems were considered in [11] in the context of the 2-alleles Wright-Fisher diffusion).(i) Consider the full fixation vector πl:=πl(i0). Remove from πl the states i0 for which πl(i0)=0 only, keeping the one, ilex, for which πl(ilex)=1. The size of this vector, say πlR, is nR=(N+K-1K-1)-K+1. Let R be the corresponding nR×nR-reduced transition matrix obtained from P after erasing the lines and columns corresponding to all the K monomorphic states except ilex. The Markov chain It conditioned to exit in the extremal state ilex only admits the stochastic transition matrix:
(3.13)Rl:=DπlR-1RDπlR.
It is obtained from R after a diagonal Doob transform based on πlR. The chain governed by Rl admits a unique absorbing state which is ilex. The entries of Rl are
(3.14)Rl(i,i′)=πlR(i′)πlR(i)R(i,i′),
and for this new conditioned Markov chain, transitions to states i′ for which πlR(i′)>πlR(i) are favored.(ii) Let us now consider again the fully reduced transition matrix Q obtained from P after erasing the lines and columns corresponding to all the K monomorphic states. The matrix Q is substochastic and irreducible, with Q1<1. The law of the process corresponding to It conditioned on avoiding the monomorphic states before t evolves as follows: with τ:=⋀l=1Kτl being the time needed for first hitting one of the extremal states for It, let πt(i)=ℙ(It=i∣τ>t). Then, with πt=πt(i),i∈SK,N∖SK,Nex,
(3.15)πt+1*=πt*Qπt*Q1
is the nonlinear Master Equation governing its evolution [12]. The reduced state-space dimension of this Markov chain is nQ=(N+K-1K-1)-K and πt→t↑∞π∞ where π∞ is the left Perron probability eigenvector of Q associated to the dominant Perron eigenvalue ρQ<1, namely, ρQπ∞*=π∞*Q.  (π∞ is called a Yaglom limit (see [13]) or a quasistationary distribution). If the process is started using this limiting quasistationary distribution, it remains in the same state over time and the fixation time τ is geometrically distributed with success probability ρQ.(iii) One can define another stochastic process I¯t which admits the stochastic transition matrix: Q¯:=DQ1-1Q again defined on the reduced state-space. For each t, we have
(3.16)Q¯(i,i′)=ℙ(It+1=i′∣It=i,τ>1),
and the conditioning on nonfixation occurs at each transition time. This process is an ergodic Markov chain with invariant probability measure solving π¯eq*=π¯eq*Q¯. It has the following closed-form determinantal expression (see [14, Section  6] and [15, page 1559]):(3.17)π¯eq(i)=|(I-Q¯)[i,i]|∑i|(I-Q¯)[i,i]|,i∈SK,N∖SK,Nex,
where (I-Q¯)[i,i] is the submatrix resulting from the deletion of row i and column i of I-Q¯. The question as to whether the process governed by Q¯ is reversible or not arises. Defining the transition matrix of the time-reversed process Q⟵ by
(3.18)Q⟵*=Dπ¯eqQ¯Dπ¯eq-1,
it does not hold that Q⟵=Q¯ and so detailed balance does not hold. Indeed, Q¯ is similar to the transition matrix of an ergodic Wright-Fisher model and Wright-Fisher chains are not reversible.(iv) If we condition on nonfixation in the remote future (see [16] for additional details), we get a Markov chain whose stochastic transition matrix is
(3.19)Q̃=ρQ-1Dψ∞-1QDψ∞.
Here ψ∞ is the positive right Perron eigenvector of Q associated to the Perron eigenvalue ρQ<1 satisfying ρQψ∞=Qψ∞. This vector can be chosen so that ∑kπ∞,kψ∞,k=1, where π∞ is again the left Perron probability eigenvector of Q associated to ρQ<1 (see [17]). With π̃t(·)=lims↑∞ℙ(It=·∣τ>t+s), we have π̃t+1*=π̃t*Q̃. The process governed by Q̃ is an ergodic Markov chain whose invariant probability measure is π̃∞=π∞∘ψ∞, the Schur product of π∞ and ψ∞ with kth entry π̃∞,k=π∞,kψ∞,k.For the last three conditionings, it is difficult to extract some information on the limiting distribution, eitherπ∞ or π¯eq or π̃∞, respectively. This is because it would suppose to solve the eigenvalue problems explicitly which is out of reach, at least theoretically. However, assuming a diploid population with a polymorphic equilibrium state xeq for p, we expect that these distributions will present a global (local) maximum (minimum) near xeq if xeq is stable (unstable). These limiting distributions should be more sharply peaked around the extremum if we consider the conditioning (iv) compared to (ii) because, the latter conditioning being more stringent than the former, it should charge more heavily the sample paths that stay away from the monomorphic states.Finally, we would like to stress that all these considerations are also relevant in the context of another fundamental stochastic model arising in the context of evolutionary genetics. We will give some elements of how to proceed with this model presenting very different properties.


## 3.1.1. The Model and Its First Properties

Consider an allelic population with constant sizeN. In the haploid (diploid) case, N is (twice) the number of real individuals. Let i:=ik and i′:=ik′,k=1,…,K be two vectors of integers quantifying the size of the allelic populations at two consecutive generations t and t+1. We will let SK,N={iintegers:|i|=∑k=1Kik=N}. Suppose that the stochastic evolutionary dynamics is now given by a Markov chain whose one-step transition matrix P from states I=i to I′=i′ is given by the multinomial Wright-Fisher model:(3.1)ℙ(It+1′=i′∣It=i)=:P(i,i′)=(Ni1′⋯iK′)∏k=1Kpk(iN)ik′.
Therefore, given It=i, to form the next generation, each allele chooses its type independently of the others with probability p, where p=(pk,k=1,…,K) is given either by (2.4) in the haploid case or by (2.22) in the diploid case. In the diploid case, the mechanism p is assumed to present a unique polymorphic state, either stable or unstable.The state-space dimension of the Markov chain governed by (3.1) is n=(N+K-1K-1) (the number of compositions of integer N into K nonnegative parts which is also the number of ways to assign N indistinguishable balls into K distinguishable boxes). To view P(i,i′) as a standard transition matrix of some Markov chain, we need first to order the states i and i′ in (3.1). Starting from the bottom right corner of P states should be arranged in decreasing order when listing the entries of P moving up and left along the lines and columns, respectively; or equivalently, starting from the top left corner of P states should be arranged in increasing order when moving down and right.For example, we can order the statesi∈SK,N in decreasing order from n to 1, as follows. Let H(i)=(1/N)∑k=1K(N+1)-(k-1)ik be a base-(N+1) code of the state i that will serve as a ranking function. The largest state in for which H(i) is maximal (equal to 1) is in*:=(N,0,…,0). Given a state i, define the subsequent state in decreasing order as
(3.2)σ-(i)=argminj:H(j)<H(i)(H(i)-H(j)).
Then in*=(N,0,…,0),in-1*=(N-1,1,0,…,0),…,i1*=(0,…,0,N) and to pass from state i to the next state σ-(i) in this decreasing sequence, there is a unique δi with entries in ℤ satisfying |δi|=0 and such that σ-(i)=i-δi. This way to order the n states is consistent with the reverse lexicographic order. For instance, if N=3,K=4, there are (63)=20 states ordered in decreasing order as follows:3000>2100>2010>2001>1200>1110>1101>1020>1011>1002>0300>0210>0201>0120>0111>0102>0030>0021>0012>0003.The way the digits are propagated from left to right is clear and the consecutiveδi can easily be obtained. Proceeding in this way to order states, P(i,i′) is a well-defined conventional object (matrix).From (3.1), the marginal transition matrix from i to Ik′=ik′ is binomial bin(N,pk(i/N)) with:
(3.3)P(i,ik′)=(Nik′)pk(iN)(1-pk(iN))N-ik′.
Given I=i, the kth component Ik′ of the updated state is random with
(3.4)𝔼i(Ik′N)=pk(iN),σi2(Ik′N)=pk(i/N)(1-pk(i/N))N,
suggesting that, in the large N population limit, the deterministic evolutionary dynamics should be recovered. Indeed, from (3.1), the Laplace-Stieltjes transform of the joint law of It+1′∣It=i reads
(3.5)𝔼i(e-∑k=1KλkIk′/N)=(∑k=1Kpk(iN)e-λk/N)N~N↑∞e-∑k=1Kλkpk(i/N).
From the strong law of large numbers therefore, if ik:=⌊Nxk⌋,k=1,…,K, then, given x=xk,k=1,…,K,
(3.6)Ik′N→N↑∞a.s.xk′=pk(x)
which is (2.22): when the population under study is very large, random fluctuations as modelled by (3.1) can be ignored so that the gene frequencies evolve deterministically.


## 3.1.2. Fixation Probabilities

Letel be the K-null vector except for its lth entry which is 1. The extremal pure states SK,Nex:={ilex:=Nel,l=1,…,K} are all absorbing for this Markov chain because pk(ilex/N)=δk,l and, from (3.1), any additional fixed point which p could have on the boundary-faces of SK which are not points does not give rise to an absorbing state for P. Under our assumptions on p, the chain is not recurrent, rather it is transient. Depending on the initial condition, say i0, the chain will necessarily end up in one of the extremal states ilex, with some fixation probability, say πl(i0), which can be computed as follows. Let πl:=πl(i),i∈SK,N be the harmonic function of the Wright-Fisher Markov chain which is the smallest solution to the boundary problem:(3.7)(I-P)πl=0ifi∈SK,N∖SK,Nex,πl=1(i=ilex)ifi∈SK,Nex.
We also have
(3.8)ℙ(Iτ=Nel∣I0=i0)=πl(i0),
where τ (<∞ almost surely) is the random hitting time of SK,Nex for It and the πl(i0)s are normalized so as ∑lπl(i0)=1. Thus πl(i0) are the searched probabilities to end up in state ilex starting from state i0.From (3.7), πl(i0) is known if i0∈SK,Nex. The remaining unknown restriction, say πlQ, of πl=πl(i0) to the nonextremal states is easily seen to be(3.9)πlQ=(I-Q)-1pilex,i0∈SK,N∖SK,Nex.
In (3.9), Q is obtained from P after erasing the lines and columns corresponding to all the K extremal states and pilex is the ilex-column of P where the entries corresponding to the extremal states have been deleted. When dealing with the Wright-Fisher model, Q is a positive matrix and also pilex>0, therefore πl(i0)>0 for all i0∈SK,N∖SK,Nex and this for each l: starting from any state i0 which is not an extremal state, there is a positive probability to hit any of the extremal states. Fixation of the state ilex means extinction of the remaining monomorphic states. It would therefore be of interest to understand the structure of the set 𝒜l={i0:πl(i0)>∑k≠lπk(i0)}={i0:πl(i0)>1/2}, for each l, which is the stochastic version of the attraction basin of ilex, especially when l is the label of the extremal state with largest fitness wk,k. If i0∈𝒜l indeed, the probability to end up in ilex is larger than the probability to end up in any other extremal state.Unfortunately, the development of the inverse ofI-Q in terms of its adjugate matrix in (3.9) shows that these fixation probabilities have a very complex determinantal alternating structure and the question of identifying 𝒜l is very complex.Withπl being the solution (3.9) to the Dirichlet problem with boundary conditions (3.7), the equilibrium measure of the chain therefore is
(3.10)πeq(i0):=∑l=1Kπl(i0)δilex,
which depends on i0. Necessarily, one allele will fixate and there is no polymorphic equilibrium state even when dealing with diploid populations. Which allele and with what probability will depend on the initial condition. Thanks to fluctuations, the picture therefore looks very different from the one pertaining to the deterministic theory.Of importance also is the time it takes to get extinct. It relies on similar techniques. For instance, the expected overall fixation timeα(i0):=𝔼i0(τ) solves the boundary problem:
(3.11)(I-P)α=1,i0∈SK,N∖SK,Nex,α=0,i0∈SK,Nex,
where α:=α(i0),i0∈SK,N. The restriction αQ of α to the nonextremal states therefore is (3.12)αQ=(I-Q)-11,i0∈SK,N∖SK,Nex.


## 3.1.3. ConditioningIt on Nonfixation

There are four places where questions relative to conditioning on fixation are relevant in this context.  (Similar conditioning problems were considered in [11] in the context of the 2-alleles Wright-Fisher diffusion).(i) Consider the full fixation vector πl:=πl(i0). Remove from πl the states i0 for which πl(i0)=0 only, keeping the one, ilex, for which πl(ilex)=1. The size of this vector, say πlR, is nR=(N+K-1K-1)-K+1. Let R be the corresponding nR×nR-reduced transition matrix obtained from P after erasing the lines and columns corresponding to all the K monomorphic states except ilex. The Markov chain It conditioned to exit in the extremal state ilex only admits the stochastic transition matrix:
(3.13)Rl:=DπlR-1RDπlR.
It is obtained from R after a diagonal Doob transform based on πlR. The chain governed by Rl admits a unique absorbing state which is ilex. The entries of Rl are
(3.14)Rl(i,i′)=πlR(i′)πlR(i)R(i,i′),
and for this new conditioned Markov chain, transitions to states i′ for which πlR(i′)>πlR(i) are favored.(ii) Let us now consider again the fully reduced transition matrix Q obtained from P after erasing the lines and columns corresponding to all the K monomorphic states. The matrix Q is substochastic and irreducible, with Q1<1. The law of the process corresponding to It conditioned on avoiding the monomorphic states before t evolves as follows: with τ:=⋀l=1Kτl being the time needed for first hitting one of the extremal states for It, let πt(i)=ℙ(It=i∣τ>t). Then, with πt=πt(i),i∈SK,N∖SK,Nex,
(3.15)πt+1*=πt*Qπt*Q1
is the nonlinear Master Equation governing its evolution [12]. The reduced state-space dimension of this Markov chain is nQ=(N+K-1K-1)-K and πt→t↑∞π∞ where π∞ is the left Perron probability eigenvector of Q associated to the dominant Perron eigenvalue ρQ<1, namely, ρQπ∞*=π∞*Q.  (π∞ is called a Yaglom limit (see [13]) or a quasistationary distribution). If the process is started using this limiting quasistationary distribution, it remains in the same state over time and the fixation time τ is geometrically distributed with success probability ρQ.(iii) One can define another stochastic process I¯t which admits the stochastic transition matrix: Q¯:=DQ1-1Q again defined on the reduced state-space. For each t, we have
(3.16)Q¯(i,i′)=ℙ(It+1=i′∣It=i,τ>1),
and the conditioning on nonfixation occurs at each transition time. This process is an ergodic Markov chain with invariant probability measure solving π¯eq*=π¯eq*Q¯. It has the following closed-form determinantal expression (see [14, Section  6] and [15, page 1559]):(3.17)π¯eq(i)=|(I-Q¯)[i,i]|∑i|(I-Q¯)[i,i]|,i∈SK,N∖SK,Nex,
where (I-Q¯)[i,i] is the submatrix resulting from the deletion of row i and column i of I-Q¯. The question as to whether the process governed by Q¯ is reversible or not arises. Defining the transition matrix of the time-reversed process Q⟵ by
(3.18)Q⟵*=Dπ¯eqQ¯Dπ¯eq-1,
it does not hold that Q⟵=Q¯ and so detailed balance does not hold. Indeed, Q¯ is similar to the transition matrix of an ergodic Wright-Fisher model and Wright-Fisher chains are not reversible.(iv) If we condition on nonfixation in the remote future (see [16] for additional details), we get a Markov chain whose stochastic transition matrix is
(3.19)Q̃=ρQ-1Dψ∞-1QDψ∞.
Here ψ∞ is the positive right Perron eigenvector of Q associated to the Perron eigenvalue ρQ<1 satisfying ρQψ∞=Qψ∞. This vector can be chosen so that ∑kπ∞,kψ∞,k=1, where π∞ is again the left Perron probability eigenvector of Q associated to ρQ<1 (see [17]). With π̃t(·)=lims↑∞ℙ(It=·∣τ>t+s), we have π̃t+1*=π̃t*Q̃. The process governed by Q̃ is an ergodic Markov chain whose invariant probability measure is π̃∞=π∞∘ψ∞, the Schur product of π∞ and ψ∞ with kth entry π̃∞,k=π∞,kψ∞,k.For the last three conditionings, it is difficult to extract some information on the limiting distribution, eitherπ∞ or π¯eq or π̃∞, respectively. This is because it would suppose to solve the eigenvalue problems explicitly which is out of reach, at least theoretically. However, assuming a diploid population with a polymorphic equilibrium state xeq for p, we expect that these distributions will present a global (local) maximum (minimum) near xeq if xeq is stable (unstable). These limiting distributions should be more sharply peaked around the extremum if we consider the conditioning (iv) compared to (ii) because, the latter conditioning being more stringent than the former, it should charge more heavily the sample paths that stay away from the monomorphic states.Finally, we would like to stress that all these considerations are also relevant in the context of another fundamental stochastic model arising in the context of evolutionary genetics. We will give some elements of how to proceed with this model presenting very different properties.


## 3.2. TheK-Alleles Moran Model

We now focus on the Moran model.


### 3.2.1. The Multiallelic Moran Model

Letα,β∈{1,…,K}. In the Moran version of the stochastic evolution,  given It=I=i, the only accessible values of I′ are the neighboring states: iα,β′:=i+dα,β where dα,β*:=(0,…,0,-1,0,…,1,0,…,0). Here -1 is in position α and 1 in position β≠α corresponding to the transfer of an individual from the box α (if nonempty) to the box β. With n(i)=#{k:ik>0} being the number of nonempty entries of i, there are n(i)(K-1)≤K(K-1) accessible states from i. The Moran stochastic evolutionary dynamics is now given by a Markov chain whose one-step transition matrix P from states I=i to I′=i′ is(3.20)ℙ(It+1=i′∣It=i)=0ifi′≠iα,β′,ℙ(It+1=iα,β′∣It=i)=:P(i,iα,β′)=iαNpβ(iN),
where p=(pβ(i/N),β=1,…,K) is given either by (2.4) in the haploid case or by (2.22) in the diploid case.SummingP(i,iα,β′) over α,β, β≠α in (3.20), we get the holding probability
(3.21)ℙ(It+1=i∣It=i)=1-∑α,β:β≠αiαNpβ(iN)=∑αiαNpα(iN),
completing the characterization of the K-alleles Moran model. Under our assumptions on p, the holding probabilities are equal to 1 only for the extremal states i∈SK,Nex which are therefore the only absorbing states of the Moran chain, just like for the Wright-Fisher model. The drift at state i is(3.22)𝔼(It+1-It∣It=i)=∑αiαN∑β≠αpβ(iN)dα,β.Let us compute the Laplace-Stieltjes Transform ofI′ in the context of a Moran model. Omitting the argument i/N in pβ(i/N), we get the factorized form: (3.23)𝔼i(e-∑kλkIk′)=∑α,β:α≠βe-∑kλkiα,β′(k)P(i,iα,β′)+e-∑kλkik∑βiβNpβ=e-∑kλkik(∑α,β:α≠βe-∑kλkdα,β(k)P(i,iα,β′)+∑βiβNpβ)=e-∑kλkik(∑α,β:α≠βe-(λβ-λα)iαNpβ+∑βiβNpβ)=e-∑kλkik(∑βe-λβpβ∑α≠βiαNeλα+∑βiβNpβ)=e-∑kλkik(∑βe-λβpβ(∑αiαNeλα-iβNeλβ)+∑βiβNpβ)=(e-∑kλkik)(∑αiαNeλα)(∑βe-λβpβ).
Putting λl=0 if l≠k, the kth marginal reads(3.24)𝔼i(e-λkIk′)=e-λkik(1-ikN+eλkikN)(1-pk+e-λkpk)
showing that Ik(t) is of the random walk type. Indeed, we get ℙi(Ik′=ik′)=0 if ik′≠ik±1 or ik′≠ik and (3.25)ℙi(Ik′=ik)=(1-ikN)(1-pk)+ikNpk,ℙi(Ik′=ik+1)=(1-ikN)pk=∑l≠kP(i,il,k′),ℙi(Ik′=ik-1)=ikN(1-pk)=∑l≠kP(i,ik,l′).
We have (3.26)𝔼i(Ik′N)=ikN+1N(pk-ikN);σi2(Ik′N)=1N2(ikN(1-ikN)+pk(1-pk)).
There is only a small correction (of order N-1) of the updated mean to its current value and fluctuations around the mean are small too (of order N-1). The evolution process is very slow.


### 3.2.2. Fixation Probabilities

As a random walk model, the Moran model has a much simpler transition matrixP of the Jacobi type. The equilibrium measure of the chain again is(3.27)πeq:=∑l=1Kπl(i0)δilex,
where πl again solves the Dirichlet problem (3.7) but with this new simpler Jacobi P. For the Moran model, the explicit expression (3.9) of the fixation probability simplifies a little bit because pilex(i)≠0 only for the K-1 neighboring states of ilex that is {i:i+dα,β=ilex} for some dα,β.WhenK=2 (2 alleles), the random walk transition probabilities (p1+p2=1) (3.28)ℙ(It+1=i1,2′∣It=i)=:P(i,i+(-11))=i1Np2(iN),ℙ(It+1=i2,1′∣It=i)=:P(i,i+(1-1))=i2Np1(iN)
are the probabilities that the first component It,1 of It moves down and up by one unit, respectively. In this case, the Dirichlet problem giving the fixation probabilities solves explicitly. With ϕ(i0)=1+∑i=1i0-1∏i1=1i(i1p2(i/N)/i2p1(i/N)) (ϕ(0)=0) being the harmonic function of the 2-alleles chain, we easily get that (3.29)π1(i0,N-i0)=ϕ(i0)ϕ(N)
is the probability that the extremal state i1ex=(N,0) is reached given i0=(i0,N-i0). Assuming a model with multiplicative fitnesses: pα(i/N)=(iα/N)(wα/ω(i/N)), then ϕ takes the simple form (i1+i2=N) (3.30)ϕ(i0)=1+∑i=1i0-1∏i1=1ii1p2(i/N)i2p1(i/N)=1+∑i=1i0-1(w2w1)i
showing (see [1, page 109]) that(3.31)π1(i0,N-i0)=1-(w2/w1)i01-(w2/w1)N.
Assuming w1=1+s/N and w2=1, putting i0=[Nx0], for large N, we get [18] (3.32)π1(Nx0,N-i0)~1-e-sx01-e-s.In the general fitness case(3.33)ϕ(i0)=1+∑i=1i0-1‍∏i1=1ii1p2(i/N)i2p1(i/N)=1+∑i=1i0-1‍∏i1=1i(Wi/N)2(Wi/N)1,
where, as usual, (Wi/N)k=∑l=12wk,lil/N,k,l=1,2, leading to (3.34)(WiN)1=w1,2+(w1,1-w1,2)i1N,(WiN)2=w2,2+(w2,1-w2,2)i1N.This 2-alleles exact solution can be used in the full diploidK-alleles Moran case with multiplicative fitnesses. Indeed, from this, with i0=(i1,…,iK), the fixation probability of Al can be conjectured to be approximated qualitatively by(3.35)πl(ilex)=1,πl(i0)=1-((∑k≠likwk/(N-il))/wl)il1-((∑k≠likwk/(N-il))/wl)Nifi0≠ilex.
This can be justified as follows: mark one particular box with size il, corresponding to the allele Al with fitness wl. Then clump the K-1 remaining boxes into a single box with size N-il, corresponding to a fictitious allele with average fitness ∑k≠likwk/(N-il). We are left with a 2-alleles Moran multiplicative fitness model for which (from the 2-alleles exact solution) the fixation probability of Al is given by (3.35). This formula constitutes sort of a mean field approximation to the full Dirichlet problem associated to the Moran model.Assumingwk~1+sk/N, a Kimura-like approximation of (3.35) would lead for large N to (3.36)πl(Nx0)~1-e(xl/(1-xl))∑k≠lxksk1-e(1/(1-xl))∑k≠lxksk,
where x0=(x1,…,xK) is now a point of the continuous K-simplex different from the l-unit vector el:=(0,…,0,1,0,…,0).From (3.35), when the fitness of allele Al is large (small), compared to the average fitness of the remaining alleles, then the fixation probability of Al gets close to 1 (resp., to 0). Note also that the larger il is, the larger the fixation probability is.As required also, for allk≠l,πl(ikex)=(1-(wk/wl)0)/(1-(wk/wl)N)=0. As another particular initial configuration case, suppose that we start from the 2-alleles type state: i0=:i0(k,l)=(0,…,0,N-1,0,…,0,1,0,…,0) where the 1 is in position l (i.e., il=1) and the entry N-1 in position k≠l (i.e., ik=N-1). Although for this choice of the initial state, the fixation of Ak is very likely, there still is a positive probability that allele Al gets fixed which is seen to be from (3.35): (3.37)πl(i0(k,l))=1-wk/wl1-(wk/wl)N,
depending only on the relative fitnesses of Al and Ak (see [19, 20] for a similar expression). As N gets large, this probability gets close to 1-wk/wl if wk<wl and close to (wl/wk)N-1~0 if wk>wl.Conversely, assuming the 2-alleles type statei0 to be defined by il=N-1 and ik=1, (3.38)πl(i0)=1-(wk/wl)N-11-(wk/wl)N,
which, as required, gets close to 1 as N gets large if wk<wl and close to wl/wk if wk>wl.


### 3.2.3. ConditioningIt on Nonfixation

The conditioning developments discussed for the Wright-Fisher model are also relevant in the Moran model context substituting theP of Moran for the P of Wright-Fisher. Let us revisit condition (iii).(iii) With Q being now the reduced substochastic matrix of the full Moran transition matrix defined in (3.20), consider the stochastic process I¯t with stochastic transition matrix: Q¯:=DQ1-1Q defined on the reduced state-space with dimension nQ. The process I¯t is again an ergodic Markov chain with invariant probability measure solving π¯eq*=π¯eq*Q¯. With Q⟵ being the transition matrix of the time-reversed process it now holds that Q⟵=Q¯ and so detailed balance holds when dealing with the Moran case (see also [21]). This will be proved if we can exhibit an equilibrium probability measure π¯eq such that(3.39)π¯eq(i)Q¯(i,i′)=Q¯(i′,i)π¯eq(i′),
for all neighboring state i′=i+di of i.Withj=(j1,…,jK) any terminal state, suppose that we want to use (3.39) to compute π¯eq(j) starting from the smallest available state in the system which is j0=(0,…,0,1,N-1). This is possible because π¯eq(j) may be represented as(3.40)π¯eq(j)=π¯eq(j0)∏i=j0j-dK,1Q¯(i,i′)Q¯(i′,i),
where π¯eq(j0) can be chosen so that ∑jπ¯eq(j)=1. Let us give some details.(i) Note first, by reversing path, that for two consecutive states(i,i′), the ratio Q¯(i,i′)/Q¯(i′,i) can be computed. We have(3.41)Q¯(i,i′)Q¯(i′,i)=(Q1)i′(Q1)iQ(i,i′)Q(i′,i),
where i′=i+di for some di of the form dk,l and therefore(3.42)Q(i,i′)=ikNpl(iN),Q(i′,i)=Q(i′,i′+dl,k)=il+1Npk(i+dk,lN).
Clearly, from such a structure of the entries of Q¯, for each possible transition i→i′, the ratio Q¯(i,i′)/Q¯(i′,i) will only depend separately on a ratio involving the terminal and the initial states i′ and i (the detailed balance condition holds).(ii) Second, the sequence ofi-s in (3.40) is governed by the following path starting from j0 and ending up in the target state j: we can use the following sequence of dk,ls:(3.43)(dK,K-1⋯jK-1-1dK,K-1)(dK,K-2⋯jK-2dK,K-1)⋯(dK,1⋯j1dK,1),
filling up successively the entries of j to the left of the last entry of j0 by using the N-1 individuals of the reservoir state j0=(0,…,0,1,N-1). By doing so, each intermediate state i is separated from the next i′ by some clearly identified di, and after evaluating the probability ratio Q¯(i,i′)/Q¯(i′,i) for each consecutive states of this sequence, we are done. Because there exists a probability distribution π¯eq(j) such that the reversibility identity (3.39) holds, then this Moran process is reversible with (3.40) as its stationary distribution. Using (3.42), this constitutes an exact explicit product-form formula.


## 3.2.1. The Multiallelic Moran Model

Letα,β∈{1,…,K}. In the Moran version of the stochastic evolution,  given It=I=i, the only accessible values of I′ are the neighboring states: iα,β′:=i+dα,β where dα,β*:=(0,…,0,-1,0,…,1,0,…,0). Here -1 is in position α and 1 in position β≠α corresponding to the transfer of an individual from the box α (if nonempty) to the box β. With n(i)=#{k:ik>0} being the number of nonempty entries of i, there are n(i)(K-1)≤K(K-1) accessible states from i. The Moran stochastic evolutionary dynamics is now given by a Markov chain whose one-step transition matrix P from states I=i to I′=i′ is(3.20)ℙ(It+1=i′∣It=i)=0ifi′≠iα,β′,ℙ(It+1=iα,β′∣It=i)=:P(i,iα,β′)=iαNpβ(iN),
where p=(pβ(i/N),β=1,…,K) is given either by (2.4) in the haploid case or by (2.22) in the diploid case.SummingP(i,iα,β′) over α,β, β≠α in (3.20), we get the holding probability
(3.21)ℙ(It+1=i∣It=i)=1-∑α,β:β≠αiαNpβ(iN)=∑αiαNpα(iN),
completing the characterization of the K-alleles Moran model. Under our assumptions on p, the holding probabilities are equal to 1 only for the extremal states i∈SK,Nex which are therefore the only absorbing states of the Moran chain, just like for the Wright-Fisher model. The drift at state i is(3.22)𝔼(It+1-It∣It=i)=∑αiαN∑β≠αpβ(iN)dα,β.Let us compute the Laplace-Stieltjes Transform ofI′ in the context of a Moran model. Omitting the argument i/N in pβ(i/N), we get the factorized form: (3.23)𝔼i(e-∑kλkIk′)=∑α,β:α≠βe-∑kλkiα,β′(k)P(i,iα,β′)+e-∑kλkik∑βiβNpβ=e-∑kλkik(∑α,β:α≠βe-∑kλkdα,β(k)P(i,iα,β′)+∑βiβNpβ)=e-∑kλkik(∑α,β:α≠βe-(λβ-λα)iαNpβ+∑βiβNpβ)=e-∑kλkik(∑βe-λβpβ∑α≠βiαNeλα+∑βiβNpβ)=e-∑kλkik(∑βe-λβpβ(∑αiαNeλα-iβNeλβ)+∑βiβNpβ)=(e-∑kλkik)(∑αiαNeλα)(∑βe-λβpβ).
Putting λl=0 if l≠k, the kth marginal reads(3.24)𝔼i(e-λkIk′)=e-λkik(1-ikN+eλkikN)(1-pk+e-λkpk)
showing that Ik(t) is of the random walk type. Indeed, we get ℙi(Ik′=ik′)=0 if ik′≠ik±1 or ik′≠ik and (3.25)ℙi(Ik′=ik)=(1-ikN)(1-pk)+ikNpk,ℙi(Ik′=ik+1)=(1-ikN)pk=∑l≠kP(i,il,k′),ℙi(Ik′=ik-1)=ikN(1-pk)=∑l≠kP(i,ik,l′).
We have (3.26)𝔼i(Ik′N)=ikN+1N(pk-ikN);σi2(Ik′N)=1N2(ikN(1-ikN)+pk(1-pk)).
There is only a small correction (of order N-1) of the updated mean to its current value and fluctuations around the mean are small too (of order N-1). The evolution process is very slow.


## 3.2.2. Fixation Probabilities

As a random walk model, the Moran model has a much simpler transition matrixP of the Jacobi type. The equilibrium measure of the chain again is(3.27)πeq:=∑l=1Kπl(i0)δilex,
where πl again solves the Dirichlet problem (3.7) but with this new simpler Jacobi P. For the Moran model, the explicit expression (3.9) of the fixation probability simplifies a little bit because pilex(i)≠0 only for the K-1 neighboring states of ilex that is {i:i+dα,β=ilex} for some dα,β.WhenK=2 (2 alleles), the random walk transition probabilities (p1+p2=1) (3.28)ℙ(It+1=i1,2′∣It=i)=:P(i,i+(-11))=i1Np2(iN),ℙ(It+1=i2,1′∣It=i)=:P(i,i+(1-1))=i2Np1(iN)
are the probabilities that the first component It,1 of It moves down and up by one unit, respectively. In this case, the Dirichlet problem giving the fixation probabilities solves explicitly. With ϕ(i0)=1+∑i=1i0-1∏i1=1i(i1p2(i/N)/i2p1(i/N)) (ϕ(0)=0) being the harmonic function of the 2-alleles chain, we easily get that (3.29)π1(i0,N-i0)=ϕ(i0)ϕ(N)
is the probability that the extremal state i1ex=(N,0) is reached given i0=(i0,N-i0). Assuming a model with multiplicative fitnesses: pα(i/N)=(iα/N)(wα/ω(i/N)), then ϕ takes the simple form (i1+i2=N) (3.30)ϕ(i0)=1+∑i=1i0-1∏i1=1ii1p2(i/N)i2p1(i/N)=1+∑i=1i0-1(w2w1)i
showing (see [1, page 109]) that(3.31)π1(i0,N-i0)=1-(w2/w1)i01-(w2/w1)N.
Assuming w1=1+s/N and w2=1, putting i0=[Nx0], for large N, we get [18] (3.32)π1(Nx0,N-i0)~1-e-sx01-e-s.In the general fitness case(3.33)ϕ(i0)=1+∑i=1i0-1‍∏i1=1ii1p2(i/N)i2p1(i/N)=1+∑i=1i0-1‍∏i1=1i(Wi/N)2(Wi/N)1,
where, as usual, (Wi/N)k=∑l=12wk,lil/N,k,l=1,2, leading to (3.34)(WiN)1=w1,2+(w1,1-w1,2)i1N,(WiN)2=w2,2+(w2,1-w2,2)i1N.This 2-alleles exact solution can be used in the full diploidK-alleles Moran case with multiplicative fitnesses. Indeed, from this, with i0=(i1,…,iK), the fixation probability of Al can be conjectured to be approximated qualitatively by(3.35)πl(ilex)=1,πl(i0)=1-((∑k≠likwk/(N-il))/wl)il1-((∑k≠likwk/(N-il))/wl)Nifi0≠ilex.
This can be justified as follows: mark one particular box with size il, corresponding to the allele Al with fitness wl. Then clump the K-1 remaining boxes into a single box with size N-il, corresponding to a fictitious allele with average fitness ∑k≠likwk/(N-il). We are left with a 2-alleles Moran multiplicative fitness model for which (from the 2-alleles exact solution) the fixation probability of Al is given by (3.35). This formula constitutes sort of a mean field approximation to the full Dirichlet problem associated to the Moran model.Assumingwk~1+sk/N, a Kimura-like approximation of (3.35) would lead for large N to (3.36)πl(Nx0)~1-e(xl/(1-xl))∑k≠lxksk1-e(1/(1-xl))∑k≠lxksk,
where x0=(x1,…,xK) is now a point of the continuous K-simplex different from the l-unit vector el:=(0,…,0,1,0,…,0).From (3.35), when the fitness of allele Al is large (small), compared to the average fitness of the remaining alleles, then the fixation probability of Al gets close to 1 (resp., to 0). Note also that the larger il is, the larger the fixation probability is.As required also, for allk≠l,πl(ikex)=(1-(wk/wl)0)/(1-(wk/wl)N)=0. As another particular initial configuration case, suppose that we start from the 2-alleles type state: i0=:i0(k,l)=(0,…,0,N-1,0,…,0,1,0,…,0) where the 1 is in position l (i.e., il=1) and the entry N-1 in position k≠l (i.e., ik=N-1). Although for this choice of the initial state, the fixation of Ak is very likely, there still is a positive probability that allele Al gets fixed which is seen to be from (3.35): (3.37)πl(i0(k,l))=1-wk/wl1-(wk/wl)N,
depending only on the relative fitnesses of Al and Ak (see [19, 20] for a similar expression). As N gets large, this probability gets close to 1-wk/wl if wk<wl and close to (wl/wk)N-1~0 if wk>wl.Conversely, assuming the 2-alleles type statei0 to be defined by il=N-1 and ik=1, (3.38)πl(i0)=1-(wk/wl)N-11-(wk/wl)N,
which, as required, gets close to 1 as N gets large if wk<wl and close to wl/wk if wk>wl.


## 3.2.3. ConditioningIt on Nonfixation

The conditioning developments discussed for the Wright-Fisher model are also relevant in the Moran model context substituting theP of Moran for the P of Wright-Fisher. Let us revisit condition (iii).(iii) With Q being now the reduced substochastic matrix of the full Moran transition matrix defined in (3.20), consider the stochastic process I¯t with stochastic transition matrix: Q¯:=DQ1-1Q defined on the reduced state-space with dimension nQ. The process I¯t is again an ergodic Markov chain with invariant probability measure solving π¯eq*=π¯eq*Q¯. With Q⟵ being the transition matrix of the time-reversed process it now holds that Q⟵=Q¯ and so detailed balance holds when dealing with the Moran case (see also [21]). This will be proved if we can exhibit an equilibrium probability measure π¯eq such that(3.39)π¯eq(i)Q¯(i,i′)=Q¯(i′,i)π¯eq(i′),
for all neighboring state i′=i+di of i.Withj=(j1,…,jK) any terminal state, suppose that we want to use (3.39) to compute π¯eq(j) starting from the smallest available state in the system which is j0=(0,…,0,1,N-1). This is possible because π¯eq(j) may be represented as(3.40)π¯eq(j)=π¯eq(j0)∏i=j0j-dK,1Q¯(i,i′)Q¯(i′,i),
where π¯eq(j0) can be chosen so that ∑jπ¯eq(j)=1. Let us give some details.(i) Note first, by reversing path, that for two consecutive states(i,i′), the ratio Q¯(i,i′)/Q¯(i′,i) can be computed. We have(3.41)Q¯(i,i′)Q¯(i′,i)=(Q1)i′(Q1)iQ(i,i′)Q(i′,i),
where i′=i+di for some di of the form dk,l and therefore(3.42)Q(i,i′)=ikNpl(iN),Q(i′,i)=Q(i′,i′+dl,k)=il+1Npk(i+dk,lN).
Clearly, from such a structure of the entries of Q¯, for each possible transition i→i′, the ratio Q¯(i,i′)/Q¯(i′,i) will only depend separately on a ratio involving the terminal and the initial states i′ and i (the detailed balance condition holds).(ii) Second, the sequence ofi-s in (3.40) is governed by the following path starting from j0 and ending up in the target state j: we can use the following sequence of dk,ls:(3.43)(dK,K-1⋯jK-1-1dK,K-1)(dK,K-2⋯jK-2dK,K-1)⋯(dK,1⋯j1dK,1),
filling up successively the entries of j to the left of the last entry of j0 by using the N-1 individuals of the reservoir state j0=(0,…,0,1,N-1). By doing so, each intermediate state i is separated from the next i′ by some clearly identified di, and after evaluating the probability ratio Q¯(i,i′)/Q¯(i′,i) for each consecutive states of this sequence, we are done. Because there exists a probability distribution π¯eq(j) such that the reversibility identity (3.39) holds, then this Moran process is reversible with (3.40) as its stationary distribution. Using (3.42), this constitutes an exact explicit product-form formula.


## 4. Concluding Remarks

This paper studies a classical population genetic model describing a one-locus multiallelic population subject to natural selection, random mating and then random genetic drifts as from the Wright-Fisher and Moran models. Although the two-alleles version of this model is fairly well studied and understood, this is not so much the case of the multiallelic one, especially in the discrete-time context which we adopt here. Let us summarize our results emphasizing that the ones which we believe are new.Considering first the deterministic updating mechanisms driven by selection, we underline that it has the form of a nonlinear Master equation suggesting that it is possible to construct an underlying Markov process governed by this Master equation. We briefly and intuitively supply such a construction.In the diploid context, we pay attention on a class of fitness matrices that leads to polymorphism. Would the equilibrium polymorphic state be unstable, we suggest that the class of potential matrices constitute a large such admissible class of fitness matrices. It contains the class of strictly ultrametric matrices which therefore deserves some interest. To the best of our knowledge, there is no discussion of such fitness models in the population genetics context. Would the polymorphic state be stable, we derive a related class of fitness matrices leading to a definite-negative mean fitness quadratic form. Some simple examples are supplied and detailed.The last section is devoted to the stochastic version of these considerations taking into account an additional important driving source of evolution, namely, the random genetic drift. When driven by selection only and in particular in the absence of mutations, the multiallelic Wright-Fisher model is a transient Markov chain whose absorbing states are the monomorphic states. We give an expression for the fixation probabilities for this process. Then, we develop four conditioning problems: conditioning on fixating in a given monomorphic state, conditioning on avoiding the extremal states before the current instant, conditioning on nonfixation at each transition time, and conditioning on avoiding the extremal states in the remote future. Finally we run into similar considerations but for the Moran model. When dealing with the fixation probabilities in this Moran context, we suggest a mean-field approximation of these probabilities which is based on a well-known explicit formula for the 2-alleles case. It concerns the case of multiplicative fitnesses only. Finally, we consider the Moran model conditioned on nonfixation at each transition time. We exploit the reversible character of this process to derive an explicit product formula for its invariant probability measure.



---

*Source: 580762-2010-08-01.xml*