# Fluctuations Analysis of Finite Discrete Birth and Death Chains with Emphasis on Moran Models with Mutations

Authors: Thierry E. Huillet
Journal: ISRN Biomathematics
 (2013)

---

## Abstract

The Moran model is a discrete-time birth and death Markov chain describing the evolution of the number of type 1 alleles in a haploid population with two alleles whose total sizeN is preserved during the course of evolution.
Bias mechanisms such as mutations or selection can affect its neutral dynamics. For the ergodic Moran model with mutations, we get interested in the fixation probabilities of a mutant, the growth rate of fluctuations, the first hitting time of the equilibrium state starting from state {0}, the first return time to the equilibrium state, and the first hitting time of {N} starting from {0}, together with the time needed for the walker to reach its invariant measure, again starting from {0}. For the last point, an appeal to the notion of Siegmund duality is necessary, and a cutoff phenomenon will be made explicit. We are interested in these problems in the large population size limit N→∞. The Moran model with mutations includes the heat exchange models of Ehrenfest and Bernoulli-Laplace as particular cases; these were studied from the point of view of the controversy concerning irreversibility (H-theorem) and the recurrence of states.


---

## Body

## 1. Introduction

In Section2, we start with generalities on discrete-time birth and death Markov chains with finite state space, whose transition matrix is of Jacobi type: spectral representation, reversibility versus H-theorem, large deviations of the empirical average, invariant measure, conditions of transience/recurrence, spectral positivity and stochastic monotonicity, distribution of hitting and first return times, and Greens function. The scale or harmonic function is shown to be of interest to determine the excursions heights would the chain be ergodic. After having fixed the background, we then proceed with the study of specific Moran chains presenting various bias mechanisms.The Moran chain model is a particular instance of a discrete-time birth and death Markov chain, describing the evolution of the number of type1 alleles in an haploid population with two alleles whose total size N is preserved during the course of evolution. Bias mechanisms such as mutations or selection can affect its neutral dynamics (there are alternative domains of science where this model makes sense, namely, economy and behavioral sciences. Here, the Markov state Xn takes on the interpretation of the number of agents that belong to the first strategy 1 (say right- or left-wing voters) at time n. See [1, 2] and references therein for a recent study with this physical image in mind.).For the ergodic Moran model with mutations [3], we get interested (in Sections 3 and 5) into the fixation probabilities of a mutant, the growth of fluctuations, the first hitting time of the equilibrium state starting from {0}, the first return time to the equilibrium state, and  the first hitting time of {N} starting from {0}, together with the time needed for the Moran walker to reach its invariant measure, again starting from {0}. The detailed growth rates of these quantities are shown to be exactly determined in terms of the mutation probabilities.(i)
For the first fixation probability point, a sharp threshold phenomenon is exhibited in the following sense: would the mutation probability favoring allele1 exceed the one favoring allele 2, the fixation probability of a type 1 mutant tends to 1 as N grows large with a 𝒪(1/N) correction. On the contrary, this fixation probability tends to 0 geometrically fast as N grows large. In the critical case when the two mutation probabilities are egal, this fixation probability goes to 1/2 from below.Under the additional presence of selection, the fixation condition can be achieved with probability1 even if the mutation probability favoring allele 2 exceeds the one favoring allele 1. The condition is that the relative fitness of allele 1 should be large enough with bound given in terms of a ratio of mutation probabilities, translating a mutation-selection balance.(ii)
For the last time to stationarity point, an appeal to the notion of Siegmund duality (requiring stochastic monotonicity) is necessary, and a cut-off phenomenon will be made explicit in Section4, following the general stream of ideas developed in [4–7]. The analysis is made simple because the spectral representation of the mutation Moran model is exactly known.We are interested into these problems in the large population size limitN→∞ (i.e., in the thermodynamic limit). From our analysis, it follows that the mutation Moran walker started at {0} reaches stationarity in about 𝒪(NlogN) time, far below the time needed to first return to the origin (or excursion’s length) which is geometrically large with N. This is because the mutational drift pushes the walker inside the bulk of the state space, making the reverse move to any of the two boundaries extremely costly.


## 2. Discrete-Time Birth and Death Processes

We start with generalities on discrete birth and death (BD) processes with finite state space, before particularizing the study to the Moran model under concern which is in this class.


### 2.1. Generalities on BD Processes with Finite State Space

LetP be a (N+1)2 tridiagonal (Jacobi) irreducible stochastic matrix
(1)P=[r0p0q1r1p1⋱⋱⋱qN-1rN-1pN-1qNrN],
with p0>0,qx,px>0, x=1,…,N-1, and qN>0, corresponding to the transition probability matrix of some ergodic discrete-time Markov chain Xn on the state space {0,…,N}. With 1 a column vector of 1s, P is stochastic if P1=1.The left (row-) eigenvectorπ′∶=(π0,π1,…,πN) associated to the eigenvalue 1 satisfies π′=π′P, where  ′ denotes transposition. Its components are given by π0 and πy=π0∏z=0y-1(pz/qz+1)>0, y=1,…,N, with π0 chosen so as ∑y=0Nπy=1. It corresponds to the invariant probability measure of the chain.The matrixP is diagonally equivalent to the symmetric matrix
(2)PS=[r0p0q1p0q1r1p1q2⋱⋱⋱pN-2qN-1rN-1pN-1qNpN-1qNrN],
so with real eigenvalues. Indeed, with Dπ∶=diag(π0,…,πN), we have P=Dπ-1/2PSDπ1/2. As a symmetric matrix, PS is diagonalizable by an orthogonal transformation, and so, P is diagonalizable.Such nearest-neighbors random walks models are well known to be reversible (detailed balance holds). Indeed,P=P←, where P← is the transition matrix of the time-reversed process, given by P←′=DπPDπ-1 (with P′ denoting the transpose of P). The transition matrix of reversible Markov chains has real eigenvalues.Under our assumptions, irreducible chains governed byP are positive recurrent, meaning that every state recurs with probability 1, the mean time between recurrence to state x being finite and equal to 1/πx.


### 2.2. The Karlin-McGregor (KMG) Spectral Theory

Consider the polynomials(𝔮y(t);y=0,…,N) in the variable t∈[-1,1], determined by 𝔮0(t)=1 and the 3-term recurrences:
(3)t𝔮0(t)=p0𝔮1(t)+r0𝔮0(t),t𝔮y(t)=py𝔮y+1(t)+ry𝔮y(t)+qy𝔮y-1(t),y∈{1,N-1},
therefore, with 𝔮0(t)=1,𝔮1(t)=(t-r0)/p0, 𝔮2(t)=[(t-r1)𝔮1(t)-q1]/p1,…. The polynomials 𝔮y(t) satisfy 𝔮y(1)=1, y≥0 and are of degree y in t. They are often called the random walk polynomials. Let tk, k=0,…,N be the zeroes of the polynomial of degree N+1:t→PN+1(t)=t𝔮N(t)-rN𝔮N(t)-qN𝔮N-1(t), namely,
(4)S:={tk:PN+1(tk)=0},
with 1=t0>t1>⋯>tN≥-1. The set S constitutes the spectrum of P (see [8, page 78]). The quantity 1-t1 is the spectral gap.Letμ(dt):=∑k=0Nμkδtk be the spectral probability measure on [-1,1] with respect to which (𝔮y(t),y≥1) are orthogonal, namely,
(5)γy∫-11𝔮x(t)𝔮y(t)μ(dt)=γy∑k=0Nμk𝔮x(tk)𝔮y(tk)=δx,y,
where γy:=πy/π0=1/∫-11𝔮y(t)2μ(dt)=∏z=0y-1(pz/qz+1), y≥0 are the potential coefficients.These polynomials are important in view of the KMG spectral representation theorem [8]. Indeed, with q(t)=(𝔮0(t),𝔮1(t),…,𝔮N(t))′, we have
(6)tq(t)=Pq(t),
and therefore, tnq(t)=Pnq(t). Selecting the row x, multiplying by 𝔮y(t), integrating with respect to μ, and applying the orthogonality relations (5), we get
(7)Pn(x,y)∶=Px(Xn=y)=γy∫-11tn𝔮x(t)𝔮y(t)μ(dt)=γy∑k=0Nμktkn𝔮x(tk)𝔮y(tk).Necessarily,μ0=π0 because by the ergodic theorem, for allx, Pn(x,y)→μ0γy=(μ0/π0)πy=πy as n↑∞. From (7), because 𝔮0(t)=1, we have P0(Xn=0)=∑k=0Nμktkn, and so the generating function of P0(Xn=0) is
(8)u0(z)∶=∑n≥1znP0(Xn=0)=∑k=0Nμk(1-ztk)-1=∫-11μ(dt)1-zt,
the Stieltjes transform of the spectral measure, satisfying u0(1)=∞ as a result of t0=1. On the other side, let ϕ0(z)=∑k≥1zkP(τ0,0*=k) be the probability generating function (pgf) of the first return time to zero τ0,0* probability. As it can easily be checked by renewal arguments (see [9, page  330]), u0(z)=1+u0(z)ϕ0(z), showing that
(9)u0(z)=11-ϕ0(z),ϕ0(z)=1-1u0(z)=1-1∑k=0Nμk(1-ztk)-1.We haveP(τ0,0*<∞)=ϕ0(1)=1 and ϕ0′(1)=E(τ0,0*)=μ0-1=π0-1.More generally, with|z|≤1, let
(10)𝔤z(x,y)=∑n≥0znPn(x,y)=(I-zP)-1(x,y)
be the Green potential function of the chain. We let 𝔤(x,y):=𝔤1(x,y), the expected sojourn time at y, starting from x. Note that 𝔤z(0,0)=u0(z) and 𝔤0(x,y)=1.(i)
Withτx,y=inf(n≥0:Xn=y∣X0=x) the first hitting time of y starting from x, by renewal arguments, the pgf of τx,y reads
(11)E(zτx,y)=𝔤z(x,y)𝔤z(x,x),withτx,x=0. We have
(12)E(zτy,x)=𝔤z(y,x)𝔤z(y,y)=γx𝔤z(x,x)γy𝔤z(y,y)𝔤z(x,y)𝔤z(x,x)=γx𝔤z(x,x)γy𝔤z(y,y)E(zτx,y).Withx<y, we have
(13)P(τx,y>n)=P(Xn<y∣X0=x)=∑0≤y′<yP(Xn=y′∣X0=x)=∑0≤y′<yPn(x,y′),P(τy,x>n)=P(Xn>x∣X0=y)=∑N≥y′>xP(Xn=y′∣X0=y)=∑N≥y′>xPn(y,y′),E(τx,y)=∑n≥0P(τx,y>n)=∑0≤y′<y𝔤(x,y′),E(τy,x)=∑n≥0P(τy,x>n)=∑N≥y′>x𝔤(y,y′).(ii)
Withτx,x*=inf(n≥1:Xn=x∣X0=x) the first return time to x, we have
(14)ϕx(z)∶=E(zτx,x*)=𝔤z(x,x)-1𝔤z(x,x)=1-1𝔤z(x,x)=:1-1ux(z),with(15)P(τx,x*<∞)=ϕx(1)=𝔤(x,x)-1𝔤(x,x).The probability that the first return time tox is finite is 1 if and only if 𝔤(x,x)=∞ (a recurrence condition), and this happens if and only if P (stochastic) is irreducible, in which case all states are in the same recurrence class. Else, if for some x, 𝔤(x,x)<∞ (a transience condition for x), this probability is smaller than 1, and one possible reason is that absorption occurs in at least one state ≠x; all states then but the absorbing ones (typically here, either {0} or {N}, or both endpoints) are in the same transience class, and P(τx,x*<∞)∈(0,1) because the walker may or not return to x before hitting one of the absorbing states. So P here is no longer irreducible, and this is a particular reducibility case where transience pops in. In the latter transient case however, by removing the rows and columns associated to the absorbing states, we are left with an irreducible substochastic transition matrix of smaller dimension to which the Karlin-McGregor spectral representation applies as well (see [8]). The dominant eigenvalue of the restricted matrix is now strictly smaller than 1.In view ofγx:=πx/π0=1/∫-11𝔮(t)2μ(dt) and in case of recurrence,
(16)E(τx,x*)=1πx.Equations (6) and (7) suggest that, with
(17)rk(x)=μkμ0𝔮x(tk),lk(y)=πy𝔮y(tk),
the vectors rk=(rk(0),…,rk(x),…,rk(N))′ and lk′=(lk(0),…,lk(y),…,lk(N)) can be chosen to be the right (column) and left (row) eigenvectors of P associated to the eigenvalue tk for which:(18)Prk=tkrk,lk′P=tklk′.In particular,l0(y)=πy𝔮y(1)=πy and r0(x)=1, as required. We have
(19)lj′rk=δj,k,
because lj′Prk=tklj′rk=tjlj′rk, with the right normalization. Equation (19) yields another orthogonality relation which is dual to the one in (5), namely,
(20)μkμ0∑x=0Nπx𝔮x(tj)𝔮x(tk)=δj,k.WithPk:=rklk′, we now get the decomposition of P into orthogonal projectors which is (7):
(21)P=∑k=0NtkPk,Pn=∑k=0NtknPk.Note the biorthogonality relations:(22)lj′Prk=tkδj,k,
from which it follows that
(23)∑x,y=0Nlj(x)Pn(x,y)rk(y)=∑x,y=0Nπx𝔮x(tj)Pn(x,y)μkμ0𝔮y(tk)=tknδj,k.There is a matrix version of all this. CallingR (resp., L) the matrix whose kth column (row) is rk (resp., lk′), (18) reads PR=RDt and LP=DtL where Dt:=   diag(t), t=(tk,k=0,…,N). Since by (19), LR=I:=   diag(1), we have LPR=DtLR=Dt, which is (22). Since RL=I, (LPR)n=Dtn=LPnR. If we let Q be the matrix whose entry (x,k) is  Q(x,k)=𝔮x(tk), then R=(1/μ0)QDμ and L=Q′Dπ is the matrix formulation of (17). Let us finally check that L and R (i.e., (17)) are indeed admissible left and right eigen matrices. The relation PR=RDt is also PQDμ=QDμDt=QDtDμ which is PQ=QDt which holds from (6). The relation LP=DtL is also Q′DπP=DtQ′Dπ or Q′P←′=DtQ′ which is P←Q=QDt or again PQ=QDt, because by the reversibility property, P←=P. If μ0R2=I (R is involutive), then μ0R=L which means Q′=QDμDπ-1. If in addition μ=π, then Q′=Q: the duality relations (5) and (20) coincide, and this is an exceptional case of self duality.The KMG spectral theory therefore gives a representation of the right and left eigenvectors ofP associated to its spectrum t, in terms of the orthogonal polynomials q(t) evaluated at t (the Q-matrix).The interpretation of the right eigenvectors ofP is as follows. We have Pnrk=tknrk. Fixing the row x, the left-hand-side of the last equality is Ex[rk(Xn)]. Therefore,
(24)Ex[rk(Xn)]=tknrk(x).Suppose that we want to computeEx[ϕ(Xn)] for some test function ϕ. Let then ϕ(y)=∑k=0Nck(ϕ)rk(y) be the decomposition of ϕ in the basis {rk}. We have 〈rj,rk〉π:=∑y=0Nπyrj(y)rk(y)=(μj/μ0)lj′rk=(μj/μ0)δj,k.Thus,〈rj,ϕ〉π=∑k=0Nck(ϕ)〈rj,rk〉π=(μj/μ0)cj(ϕ) and
(25)Ex[ϕ(Xn)]=∑k=0Ntkn·μ0μk〈rk,ϕ〉π·rk(x)=1πx∑k=0Ntkn·〈rk,ϕ〉π·lk(x).Whenever the spectral structure of the chain is known, the solutionEx[ϕ(Xn)] to the backward Kolmogorov equations are thus computable exactly. Let us now investigate some special classes of random walks (RWs).


### 2.3. Symmetric RW

Suppose thatN is even with N=2N0. Assume that a random walk Xn on {0,…,2N0} is given so that
(26)P=[01q10p1⋱⋱⋱q2N0-10p2N0-110]
is the transition matrix of some symmetric random walk reflected at the boundaries {0,2N0}. The KMG theory is also relevant for such transition matrices Pn. Let tk, k=0,…,2N0 be the zeroes of the polynomial of degree 2N0+1:t→P2N0+1(t)=t𝔮2N0(t)-𝔮2N0𝔮2N0-1(t), namely,
(27)S={tk:P2N0+1(tk)=0},
with 1=t0>t1≥⋯≥t2N0=-1. Then, μ(dt):=∑k=02N0μkδtk is the probability measure on [-1,1] with respect to which the polynomial (𝔮y(t),y≥1) associated to P is orthogonal. This measure is symmetric on [-1,1]. In particular, tN0=0 and t2N0=-1,t0=1 are eigenvalues of such Ps. When N is odd, say with N=2N0+1, the spectral measure μ(dt):=∑k=02N0+1μkδtk is again symmetric, but {0} no longer is an eigenvalue, and tN0>0.Example 1.
The simplest such example is the simple gambler random walkXn for which px=p, qx=q, x=1,…,N-1 (p+q=1) with a pure reflection at the endpoints (p0=qN=1). The invariant measure π is a truncated geometric distribution. It follows from [9, page 438] that tk=2pqcos(kπ/N),k=1,…,N-1,t0=1,tN=-1. The spectral gap tends to 1-2pq≠0 as N↑∞ if p≠1/2. The mass of the spectral measure can be seen to be equal to μN=μ0=π0=(1-(p/q))/[2(1-(p/q)N)], μk=(1-2μ0)/(N-1),  k=1,…,N-1 (with μ0=1/(2N) when p=1/2). It exhibits boundary effects causing deviation from the uniform measure on {0,…,N}. The orthogonal polynomials involve two sine functions. A closed form expression of the generating functions of P0(Xn=0) and P(τ0,0*=k) follows.
Additional Comment.Let P be the tridiagonal transition matrix associated to this RW Xn. Let P¯ be the transition matrix associated to X¯n:=N-Xn. Then, P¯(x,y)=P(N-x,N-y) has the same structure as P except that p and q are interchanged (a drift reversal property). P¯ is obtained from P by applying a permutation reversing the order of the rows and then of the columns. So P¯ and P share the same spectrum, and this is why tk is a symmetric function of p and q. However, the mass of the spectral measure of P¯ is different with: μ¯N=μ¯0=(1-(q/p))/[2(1-(q/p)N)].


### 2.4. Special Cases: Spectral Positivity

Let(28)S′=[p0q1p1⋱⋱qN-1pN-1qNpN]
be a subdiagonal matrix. It is the transpose of a superdiagonal matrix S. Then, with Dr:= diag(r0,…,rN)(29)S′S=[p0p0q1p0q1p1+q1p1q2⋱⋱⋱pN-2qN-1pN-1+qN-1pN-1qNpN-1qNpN+qN]=I-Dr+[0p0q1p0q10p1q2⋱⋱⋱pN-2qN-10pN-1qNpN-1qN0].Thus,(30)PS=2Dr-I+S′S
is the sum of a diagonal matrix and a symmetric positive definite matrix. From this, we conclude that if the holding probabilities rx≥1/2, for all x=0,…,N, then for all z∈ℝN+1∖{0},
(31)z′PSz=∑x=0N(2rx-1)|zx|2+|Sz|2>0,
and so PS and then P are positive definite with real positive eigenvalues. The Jacobi matrix P then has all its principal minors nonnegative and therefore is oscillatory (i.e., totally nonnegative and such that PN is totally positive) with all its minors nonnegative (see e.g., [10, page  99]). Oscillatory stochastic matrices have distinct positive eigenvalues, with: 1=t0>t1>⋯>tN>0. As a sequence of numbers, the right row eigenvector rk=(rk(0),…,rk(x),…,rk(N))′ associated to tk has exactly k sign changes (see [10, page  101]). And so do the numbers: (𝔮x(tk);x=0,…,N). The above condition rx≥1/2 is a simple sufficient condition allowing to exhibit examples of P that are oscillatory but of course, it is not necessary; the full condition under our assumptions being z′PSz>0 for all z∈ℝN+1∖{0} or
(32)z′Pz=∑x=0Nrxzx2+∑x=0N-1(px+qx+1)zxzx+1>0.Equivalently, ifsk are the positive eigenvalues of S′S, we need to check 2rk-1+sk>0 for all k.For a simple random walkXn whose (N+1)2-transition matrix P is spectrally nonnegative (resp., spectrally positive), there exists a symmetric random walk Ym on {0,…,2N} (resp., on {0,…,2N+1}), reflected at the endpoints, started at an even integer {0,2,4,…,2N} (resp., odd integer {1,3,…,2N+1}), such that {Xn}=d{Y2n/2}, (resp., {Xn}=d{(Y2n-1)/2}). This follows simply from adapting [11, Theorem  2.1] to our finite-dimensional case. As noted just before, the spectral measure of the random walk Ym is symmetric on [-1,1], and by passing to Xn, the spectrum is being folded: If ∑k=02Nμkδtk (resp., ∑k=02N+1μkδtk) is the symmetric spectral measure of Ym with tN=0 (resp., tN>0), then 2∑k=0Nμkδtk2 is the spectral measure of Xn. Let αy and βy be the up and down probabilities that Ym→Ym+1=Ym±1 in one step given Ym is in state y different from the endpoints, αy+βy=1, then
(33)qx=β2xβ2x-1,rx=β2xα2x-1+α2xβ2x+1,px=α2xα2x+1.This, together withp0=α1 and qN=β2N-1 (resp., qN=β2N), allows to determine recursively the transition matrix of Ym from the one of Xn. From these facts, we get the following proposition.Proposition 2.
If a BD chainXn is spectrally nonnegative, then it is stochastically monotone.Proof.
Under our hypothesis indeed,(34)px+qx+1=α2xα2x+1+β2x+2β2x+1<α2xα2x+1+β2x+1<1,
because this is also α2x<1. This condition means that Xn is stochastically monotone (see below the construction of the Siegmund dual of an RW where stochastic monotonicity is needed). See also [6, Lemma 2.4] for a different proof.Remark 3.
Suppose a BD chainXn with transition matrix P is not spectrally nonnegative, so P has strictly negative eigenvalues (tN<0). With ρ∈(0,1), consider the lazy modified stochastic matrix
(35)Pρ=(1-ρ)I+ρP.The eigenvalues ofPρ are tk(ρ)=1-ρ+ρtk, k=0,…,N with the same left and right eigenvectors as the ones of P. One can always choose ρ≤(1-tN)-1 small enough so that tN(ρ)≥0. For this ρ, Pρ is now spectrally nonnegative, and so Pρ governs a stochastically monotone chain Xn(ρ). The RW Xn(ρ), with shifted spectrum, clearly has the same invariant measure as the one of Xn.For example, solving (34) with the boundary conditions p0=α1 and qN=β2N gives the forward and backward recurrences of the unknown odd and even probability terms in the bulk:
(36)Forx=1,…,N-1:α2x+1=px(1-qx/(1-α2x-1)),α1=p0,Forx=N-1,…,1:α2x=px(1-qx+1/(1-α2(x+1))),α2N=1-qN.As a conclusion, it is easy to construct spectrally positive RWs by “squaring” two symmetric random walks, but given a nonsymmetric RW, it may be difficult to decide whether it is or not spectrally positive because one needs conversely to check whether all the aboveαs are probabilities. Spectral positivity determines stochastic monotonicity.


### 2.5. Scale Function and Excursion Height of an RW

Assume the heightH=h∈{1,…,N-1} for some excursion (sample paths of Xn between two consecutive visits to state 0). This event will be realized if and only if (i) downward paths started from h hit state {0} before hitting state h+1 and (ii) upward paths started at {0} first reach 1 (with probability p0) and then, paths started at 1 hit h without returning to {0} again in the intervening time. These two events are independent. Therefore, with τx,y, the first hitting time of y starting from x:(37)P(H=h)=p0P(τ1,h<τ1,0)P(τh,0<τh,h+1),h=1,…,N.AssumeX0=x. Let Xn∧τx,0 denote the random walk stopped when it first hits 0. Its transition matrix is P except for its first row which is (1,0,0,…,0), translating that state {0} is now absorbing. Define the scale (or harmonic) function φ of this random walk as the function which makes Mn:=φ(Xn∧τx,0) a martingale. The function φ is important because, as is well-known, for all 0<x<h≤N, with τx=τx,0∧τx,h, the first hitting time of {0,h} starting from x(38)P(Xτx=h)=φ(x)φ(h).Using this remark, we get(39)P(H=h)=p0P(τ1,h<τ1,0)P(τh,0<τh,h+1)=p0φ(1)φ(h)(1-φ(h)φ(h+1)),h∈{1,…,N-1}.Note thatP(H≥h)=p0/φ(h) if φ(1)=1. The event H=N may also occur, with probability P(H=N)=p0/φ(N).The interpretation the eventH=N takes on is twofold: either it is the maximal height of an excursion if the state {N} remains reflecting, either, if state {N} is itself assumed absorbing, P(Xτx=N)=φ(x)/φ(N)=(φ(x)/p0)P(H=N) constitutes the probability of a fixation event at {N}, starting from state x. We can thus compute the probability of a fixation before extinction at {0} occurs. It remains to compute φ (defined up to an arbitrary multiplicative constant). We wish to have Ex(Mn+1∣Mn=y)=y, leading to φ(0)=0 and
(40)φ(x)=qxφ(x-1)+rxφ(x)+pxφ(x+1),x≥1,
where px+qx+rx=1, x≥1. Imposing φ(1)=1 (fixing the arbitrary multiplicative constant), the searched “harmonic” function is φ(x)=:1+∑y=1x-1ψ(y) where ψ(y) satisfies qyψ(y-1)=pyψ(y). Thus, ψ(y)=∏z=1y(qz/pz) and
(41)φ(x)=1+∑y=1x-1∏z=1yqzpz,x≥1,φ(0):=0.Equations (39) and (41) characterize the law of the excursion height of the random walk Xn. We note that φ(x) increases with x. From (38) and (41), we get an explicit expression of the probability of the event τx,h<τx,0:
(42)P(τx,h<τx,0)=P(Xτx=h)=φ(x)φ(h),
where τx,h is the first hitting time of h starting from x, with 0<x<h≤N. In particular, if h=N with {N} also absorbing, P(τx,N<τx,0)=φ(x)/φ(N) is the fixation probability starting from x, whereas P(τx,0<τx,N)=1-P(τx,N<τx,0) is the extinction probability.


## 2.1. Generalities on BD Processes with Finite State Space

LetP be a (N+1)2 tridiagonal (Jacobi) irreducible stochastic matrix
(1)P=[r0p0q1r1p1⋱⋱⋱qN-1rN-1pN-1qNrN],
with p0>0,qx,px>0, x=1,…,N-1, and qN>0, corresponding to the transition probability matrix of some ergodic discrete-time Markov chain Xn on the state space {0,…,N}. With 1 a column vector of 1s, P is stochastic if P1=1.The left (row-) eigenvectorπ′∶=(π0,π1,…,πN) associated to the eigenvalue 1 satisfies π′=π′P, where  ′ denotes transposition. Its components are given by π0 and πy=π0∏z=0y-1(pz/qz+1)>0, y=1,…,N, with π0 chosen so as ∑y=0Nπy=1. It corresponds to the invariant probability measure of the chain.The matrixP is diagonally equivalent to the symmetric matrix
(2)PS=[r0p0q1p0q1r1p1q2⋱⋱⋱pN-2qN-1rN-1pN-1qNpN-1qNrN],
so with real eigenvalues. Indeed, with Dπ∶=diag(π0,…,πN), we have P=Dπ-1/2PSDπ1/2. As a symmetric matrix, PS is diagonalizable by an orthogonal transformation, and so, P is diagonalizable.Such nearest-neighbors random walks models are well known to be reversible (detailed balance holds). Indeed,P=P←, where P← is the transition matrix of the time-reversed process, given by P←′=DπPDπ-1 (with P′ denoting the transpose of P). The transition matrix of reversible Markov chains has real eigenvalues.Under our assumptions, irreducible chains governed byP are positive recurrent, meaning that every state recurs with probability 1, the mean time between recurrence to state x being finite and equal to 1/πx.


## 2.2. The Karlin-McGregor (KMG) Spectral Theory

Consider the polynomials(𝔮y(t);y=0,…,N) in the variable t∈[-1,1], determined by 𝔮0(t)=1 and the 3-term recurrences:
(3)t𝔮0(t)=p0𝔮1(t)+r0𝔮0(t),t𝔮y(t)=py𝔮y+1(t)+ry𝔮y(t)+qy𝔮y-1(t),y∈{1,N-1},
therefore, with 𝔮0(t)=1,𝔮1(t)=(t-r0)/p0, 𝔮2(t)=[(t-r1)𝔮1(t)-q1]/p1,…. The polynomials 𝔮y(t) satisfy 𝔮y(1)=1, y≥0 and are of degree y in t. They are often called the random walk polynomials. Let tk, k=0,…,N be the zeroes of the polynomial of degree N+1:t→PN+1(t)=t𝔮N(t)-rN𝔮N(t)-qN𝔮N-1(t), namely,
(4)S:={tk:PN+1(tk)=0},
with 1=t0>t1>⋯>tN≥-1. The set S constitutes the spectrum of P (see [8, page 78]). The quantity 1-t1 is the spectral gap.Letμ(dt):=∑k=0Nμkδtk be the spectral probability measure on [-1,1] with respect to which (𝔮y(t),y≥1) are orthogonal, namely,
(5)γy∫-11𝔮x(t)𝔮y(t)μ(dt)=γy∑k=0Nμk𝔮x(tk)𝔮y(tk)=δx,y,
where γy:=πy/π0=1/∫-11𝔮y(t)2μ(dt)=∏z=0y-1(pz/qz+1), y≥0 are the potential coefficients.These polynomials are important in view of the KMG spectral representation theorem [8]. Indeed, with q(t)=(𝔮0(t),𝔮1(t),…,𝔮N(t))′, we have
(6)tq(t)=Pq(t),
and therefore, tnq(t)=Pnq(t). Selecting the row x, multiplying by 𝔮y(t), integrating with respect to μ, and applying the orthogonality relations (5), we get
(7)Pn(x,y)∶=Px(Xn=y)=γy∫-11tn𝔮x(t)𝔮y(t)μ(dt)=γy∑k=0Nμktkn𝔮x(tk)𝔮y(tk).Necessarily,μ0=π0 because by the ergodic theorem, for allx, Pn(x,y)→μ0γy=(μ0/π0)πy=πy as n↑∞. From (7), because 𝔮0(t)=1, we have P0(Xn=0)=∑k=0Nμktkn, and so the generating function of P0(Xn=0) is
(8)u0(z)∶=∑n≥1znP0(Xn=0)=∑k=0Nμk(1-ztk)-1=∫-11μ(dt)1-zt,
the Stieltjes transform of the spectral measure, satisfying u0(1)=∞ as a result of t0=1. On the other side, let ϕ0(z)=∑k≥1zkP(τ0,0*=k) be the probability generating function (pgf) of the first return time to zero τ0,0* probability. As it can easily be checked by renewal arguments (see [9, page  330]), u0(z)=1+u0(z)ϕ0(z), showing that
(9)u0(z)=11-ϕ0(z),ϕ0(z)=1-1u0(z)=1-1∑k=0Nμk(1-ztk)-1.We haveP(τ0,0*<∞)=ϕ0(1)=1 and ϕ0′(1)=E(τ0,0*)=μ0-1=π0-1.More generally, with|z|≤1, let
(10)𝔤z(x,y)=∑n≥0znPn(x,y)=(I-zP)-1(x,y)
be the Green potential function of the chain. We let 𝔤(x,y):=𝔤1(x,y), the expected sojourn time at y, starting from x. Note that 𝔤z(0,0)=u0(z) and 𝔤0(x,y)=1.(i)
Withτx,y=inf(n≥0:Xn=y∣X0=x) the first hitting time of y starting from x, by renewal arguments, the pgf of τx,y reads
(11)E(zτx,y)=𝔤z(x,y)𝔤z(x,x),withτx,x=0. We have
(12)E(zτy,x)=𝔤z(y,x)𝔤z(y,y)=γx𝔤z(x,x)γy𝔤z(y,y)𝔤z(x,y)𝔤z(x,x)=γx𝔤z(x,x)γy𝔤z(y,y)E(zτx,y).Withx<y, we have
(13)P(τx,y>n)=P(Xn<y∣X0=x)=∑0≤y′<yP(Xn=y′∣X0=x)=∑0≤y′<yPn(x,y′),P(τy,x>n)=P(Xn>x∣X0=y)=∑N≥y′>xP(Xn=y′∣X0=y)=∑N≥y′>xPn(y,y′),E(τx,y)=∑n≥0P(τx,y>n)=∑0≤y′<y𝔤(x,y′),E(τy,x)=∑n≥0P(τy,x>n)=∑N≥y′>x𝔤(y,y′).(ii)
Withτx,x*=inf(n≥1:Xn=x∣X0=x) the first return time to x, we have
(14)ϕx(z)∶=E(zτx,x*)=𝔤z(x,x)-1𝔤z(x,x)=1-1𝔤z(x,x)=:1-1ux(z),with(15)P(τx,x*<∞)=ϕx(1)=𝔤(x,x)-1𝔤(x,x).The probability that the first return time tox is finite is 1 if and only if 𝔤(x,x)=∞ (a recurrence condition), and this happens if and only if P (stochastic) is irreducible, in which case all states are in the same recurrence class. Else, if for some x, 𝔤(x,x)<∞ (a transience condition for x), this probability is smaller than 1, and one possible reason is that absorption occurs in at least one state ≠x; all states then but the absorbing ones (typically here, either {0} or {N}, or both endpoints) are in the same transience class, and P(τx,x*<∞)∈(0,1) because the walker may or not return to x before hitting one of the absorbing states. So P here is no longer irreducible, and this is a particular reducibility case where transience pops in. In the latter transient case however, by removing the rows and columns associated to the absorbing states, we are left with an irreducible substochastic transition matrix of smaller dimension to which the Karlin-McGregor spectral representation applies as well (see [8]). The dominant eigenvalue of the restricted matrix is now strictly smaller than 1.In view ofγx:=πx/π0=1/∫-11𝔮(t)2μ(dt) and in case of recurrence,
(16)E(τx,x*)=1πx.Equations (6) and (7) suggest that, with
(17)rk(x)=μkμ0𝔮x(tk),lk(y)=πy𝔮y(tk),
the vectors rk=(rk(0),…,rk(x),…,rk(N))′ and lk′=(lk(0),…,lk(y),…,lk(N)) can be chosen to be the right (column) and left (row) eigenvectors of P associated to the eigenvalue tk for which:(18)Prk=tkrk,lk′P=tklk′.In particular,l0(y)=πy𝔮y(1)=πy and r0(x)=1, as required. We have
(19)lj′rk=δj,k,
because lj′Prk=tklj′rk=tjlj′rk, with the right normalization. Equation (19) yields another orthogonality relation which is dual to the one in (5), namely,
(20)μkμ0∑x=0Nπx𝔮x(tj)𝔮x(tk)=δj,k.WithPk:=rklk′, we now get the decomposition of P into orthogonal projectors which is (7):
(21)P=∑k=0NtkPk,Pn=∑k=0NtknPk.Note the biorthogonality relations:(22)lj′Prk=tkδj,k,
from which it follows that
(23)∑x,y=0Nlj(x)Pn(x,y)rk(y)=∑x,y=0Nπx𝔮x(tj)Pn(x,y)μkμ0𝔮y(tk)=tknδj,k.There is a matrix version of all this. CallingR (resp., L) the matrix whose kth column (row) is rk (resp., lk′), (18) reads PR=RDt and LP=DtL where Dt:=   diag(t), t=(tk,k=0,…,N). Since by (19), LR=I:=   diag(1), we have LPR=DtLR=Dt, which is (22). Since RL=I, (LPR)n=Dtn=LPnR. If we let Q be the matrix whose entry (x,k) is  Q(x,k)=𝔮x(tk), then R=(1/μ0)QDμ and L=Q′Dπ is the matrix formulation of (17). Let us finally check that L and R (i.e., (17)) are indeed admissible left and right eigen matrices. The relation PR=RDt is also PQDμ=QDμDt=QDtDμ which is PQ=QDt which holds from (6). The relation LP=DtL is also Q′DπP=DtQ′Dπ or Q′P←′=DtQ′ which is P←Q=QDt or again PQ=QDt, because by the reversibility property, P←=P. If μ0R2=I (R is involutive), then μ0R=L which means Q′=QDμDπ-1. If in addition μ=π, then Q′=Q: the duality relations (5) and (20) coincide, and this is an exceptional case of self duality.The KMG spectral theory therefore gives a representation of the right and left eigenvectors ofP associated to its spectrum t, in terms of the orthogonal polynomials q(t) evaluated at t (the Q-matrix).The interpretation of the right eigenvectors ofP is as follows. We have Pnrk=tknrk. Fixing the row x, the left-hand-side of the last equality is Ex[rk(Xn)]. Therefore,
(24)Ex[rk(Xn)]=tknrk(x).Suppose that we want to computeEx[ϕ(Xn)] for some test function ϕ. Let then ϕ(y)=∑k=0Nck(ϕ)rk(y) be the decomposition of ϕ in the basis {rk}. We have 〈rj,rk〉π:=∑y=0Nπyrj(y)rk(y)=(μj/μ0)lj′rk=(μj/μ0)δj,k.Thus,〈rj,ϕ〉π=∑k=0Nck(ϕ)〈rj,rk〉π=(μj/μ0)cj(ϕ) and
(25)Ex[ϕ(Xn)]=∑k=0Ntkn·μ0μk〈rk,ϕ〉π·rk(x)=1πx∑k=0Ntkn·〈rk,ϕ〉π·lk(x).Whenever the spectral structure of the chain is known, the solutionEx[ϕ(Xn)] to the backward Kolmogorov equations are thus computable exactly. Let us now investigate some special classes of random walks (RWs).


## 2.3. Symmetric RW

Suppose thatN is even with N=2N0. Assume that a random walk Xn on {0,…,2N0} is given so that
(26)P=[01q10p1⋱⋱⋱q2N0-10p2N0-110]
is the transition matrix of some symmetric random walk reflected at the boundaries {0,2N0}. The KMG theory is also relevant for such transition matrices Pn. Let tk, k=0,…,2N0 be the zeroes of the polynomial of degree 2N0+1:t→P2N0+1(t)=t𝔮2N0(t)-𝔮2N0𝔮2N0-1(t), namely,
(27)S={tk:P2N0+1(tk)=0},
with 1=t0>t1≥⋯≥t2N0=-1. Then, μ(dt):=∑k=02N0μkδtk is the probability measure on [-1,1] with respect to which the polynomial (𝔮y(t),y≥1) associated to P is orthogonal. This measure is symmetric on [-1,1]. In particular, tN0=0 and t2N0=-1,t0=1 are eigenvalues of such Ps. When N is odd, say with N=2N0+1, the spectral measure μ(dt):=∑k=02N0+1μkδtk is again symmetric, but {0} no longer is an eigenvalue, and tN0>0.Example 1.
The simplest such example is the simple gambler random walkXn for which px=p, qx=q, x=1,…,N-1 (p+q=1) with a pure reflection at the endpoints (p0=qN=1). The invariant measure π is a truncated geometric distribution. It follows from [9, page 438] that tk=2pqcos(kπ/N),k=1,…,N-1,t0=1,tN=-1. The spectral gap tends to 1-2pq≠0 as N↑∞ if p≠1/2. The mass of the spectral measure can be seen to be equal to μN=μ0=π0=(1-(p/q))/[2(1-(p/q)N)], μk=(1-2μ0)/(N-1),  k=1,…,N-1 (with μ0=1/(2N) when p=1/2). It exhibits boundary effects causing deviation from the uniform measure on {0,…,N}. The orthogonal polynomials involve two sine functions. A closed form expression of the generating functions of P0(Xn=0) and P(τ0,0*=k) follows.
Additional Comment.Let P be the tridiagonal transition matrix associated to this RW Xn. Let P¯ be the transition matrix associated to X¯n:=N-Xn. Then, P¯(x,y)=P(N-x,N-y) has the same structure as P except that p and q are interchanged (a drift reversal property). P¯ is obtained from P by applying a permutation reversing the order of the rows and then of the columns. So P¯ and P share the same spectrum, and this is why tk is a symmetric function of p and q. However, the mass of the spectral measure of P¯ is different with: μ¯N=μ¯0=(1-(q/p))/[2(1-(q/p)N)].


## 2.4. Special Cases: Spectral Positivity

Let(28)S′=[p0q1p1⋱⋱qN-1pN-1qNpN]
be a subdiagonal matrix. It is the transpose of a superdiagonal matrix S. Then, with Dr:= diag(r0,…,rN)(29)S′S=[p0p0q1p0q1p1+q1p1q2⋱⋱⋱pN-2qN-1pN-1+qN-1pN-1qNpN-1qNpN+qN]=I-Dr+[0p0q1p0q10p1q2⋱⋱⋱pN-2qN-10pN-1qNpN-1qN0].Thus,(30)PS=2Dr-I+S′S
is the sum of a diagonal matrix and a symmetric positive definite matrix. From this, we conclude that if the holding probabilities rx≥1/2, for all x=0,…,N, then for all z∈ℝN+1∖{0},
(31)z′PSz=∑x=0N(2rx-1)|zx|2+|Sz|2>0,
and so PS and then P are positive definite with real positive eigenvalues. The Jacobi matrix P then has all its principal minors nonnegative and therefore is oscillatory (i.e., totally nonnegative and such that PN is totally positive) with all its minors nonnegative (see e.g., [10, page  99]). Oscillatory stochastic matrices have distinct positive eigenvalues, with: 1=t0>t1>⋯>tN>0. As a sequence of numbers, the right row eigenvector rk=(rk(0),…,rk(x),…,rk(N))′ associated to tk has exactly k sign changes (see [10, page  101]). And so do the numbers: (𝔮x(tk);x=0,…,N). The above condition rx≥1/2 is a simple sufficient condition allowing to exhibit examples of P that are oscillatory but of course, it is not necessary; the full condition under our assumptions being z′PSz>0 for all z∈ℝN+1∖{0} or
(32)z′Pz=∑x=0Nrxzx2+∑x=0N-1(px+qx+1)zxzx+1>0.Equivalently, ifsk are the positive eigenvalues of S′S, we need to check 2rk-1+sk>0 for all k.For a simple random walkXn whose (N+1)2-transition matrix P is spectrally nonnegative (resp., spectrally positive), there exists a symmetric random walk Ym on {0,…,2N} (resp., on {0,…,2N+1}), reflected at the endpoints, started at an even integer {0,2,4,…,2N} (resp., odd integer {1,3,…,2N+1}), such that {Xn}=d{Y2n/2}, (resp., {Xn}=d{(Y2n-1)/2}). This follows simply from adapting [11, Theorem  2.1] to our finite-dimensional case. As noted just before, the spectral measure of the random walk Ym is symmetric on [-1,1], and by passing to Xn, the spectrum is being folded: If ∑k=02Nμkδtk (resp., ∑k=02N+1μkδtk) is the symmetric spectral measure of Ym with tN=0 (resp., tN>0), then 2∑k=0Nμkδtk2 is the spectral measure of Xn. Let αy and βy be the up and down probabilities that Ym→Ym+1=Ym±1 in one step given Ym is in state y different from the endpoints, αy+βy=1, then
(33)qx=β2xβ2x-1,rx=β2xα2x-1+α2xβ2x+1,px=α2xα2x+1.This, together withp0=α1 and qN=β2N-1 (resp., qN=β2N), allows to determine recursively the transition matrix of Ym from the one of Xn. From these facts, we get the following proposition.Proposition 2.
If a BD chainXn is spectrally nonnegative, then it is stochastically monotone.Proof.
Under our hypothesis indeed,(34)px+qx+1=α2xα2x+1+β2x+2β2x+1<α2xα2x+1+β2x+1<1,
because this is also α2x<1. This condition means that Xn is stochastically monotone (see below the construction of the Siegmund dual of an RW where stochastic monotonicity is needed). See also [6, Lemma 2.4] for a different proof.Remark 3.
Suppose a BD chainXn with transition matrix P is not spectrally nonnegative, so P has strictly negative eigenvalues (tN<0). With ρ∈(0,1), consider the lazy modified stochastic matrix
(35)Pρ=(1-ρ)I+ρP.The eigenvalues ofPρ are tk(ρ)=1-ρ+ρtk, k=0,…,N with the same left and right eigenvectors as the ones of P. One can always choose ρ≤(1-tN)-1 small enough so that tN(ρ)≥0. For this ρ, Pρ is now spectrally nonnegative, and so Pρ governs a stochastically monotone chain Xn(ρ). The RW Xn(ρ), with shifted spectrum, clearly has the same invariant measure as the one of Xn.For example, solving (34) with the boundary conditions p0=α1 and qN=β2N gives the forward and backward recurrences of the unknown odd and even probability terms in the bulk:
(36)Forx=1,…,N-1:α2x+1=px(1-qx/(1-α2x-1)),α1=p0,Forx=N-1,…,1:α2x=px(1-qx+1/(1-α2(x+1))),α2N=1-qN.As a conclusion, it is easy to construct spectrally positive RWs by “squaring” two symmetric random walks, but given a nonsymmetric RW, it may be difficult to decide whether it is or not spectrally positive because one needs conversely to check whether all the aboveαs are probabilities. Spectral positivity determines stochastic monotonicity.


## 2.5. Scale Function and Excursion Height of an RW

Assume the heightH=h∈{1,…,N-1} for some excursion (sample paths of Xn between two consecutive visits to state 0). This event will be realized if and only if (i) downward paths started from h hit state {0} before hitting state h+1 and (ii) upward paths started at {0} first reach 1 (with probability p0) and then, paths started at 1 hit h without returning to {0} again in the intervening time. These two events are independent. Therefore, with τx,y, the first hitting time of y starting from x:(37)P(H=h)=p0P(τ1,h<τ1,0)P(τh,0<τh,h+1),h=1,…,N.AssumeX0=x. Let Xn∧τx,0 denote the random walk stopped when it first hits 0. Its transition matrix is P except for its first row which is (1,0,0,…,0), translating that state {0} is now absorbing. Define the scale (or harmonic) function φ of this random walk as the function which makes Mn:=φ(Xn∧τx,0) a martingale. The function φ is important because, as is well-known, for all 0<x<h≤N, with τx=τx,0∧τx,h, the first hitting time of {0,h} starting from x(38)P(Xτx=h)=φ(x)φ(h).Using this remark, we get(39)P(H=h)=p0P(τ1,h<τ1,0)P(τh,0<τh,h+1)=p0φ(1)φ(h)(1-φ(h)φ(h+1)),h∈{1,…,N-1}.Note thatP(H≥h)=p0/φ(h) if φ(1)=1. The event H=N may also occur, with probability P(H=N)=p0/φ(N).The interpretation the eventH=N takes on is twofold: either it is the maximal height of an excursion if the state {N} remains reflecting, either, if state {N} is itself assumed absorbing, P(Xτx=N)=φ(x)/φ(N)=(φ(x)/p0)P(H=N) constitutes the probability of a fixation event at {N}, starting from state x. We can thus compute the probability of a fixation before extinction at {0} occurs. It remains to compute φ (defined up to an arbitrary multiplicative constant). We wish to have Ex(Mn+1∣Mn=y)=y, leading to φ(0)=0 and
(40)φ(x)=qxφ(x-1)+rxφ(x)+pxφ(x+1),x≥1,
where px+qx+rx=1, x≥1. Imposing φ(1)=1 (fixing the arbitrary multiplicative constant), the searched “harmonic” function is φ(x)=:1+∑y=1x-1ψ(y) where ψ(y) satisfies qyψ(y-1)=pyψ(y). Thus, ψ(y)=∏z=1y(qz/pz) and
(41)φ(x)=1+∑y=1x-1∏z=1yqzpz,x≥1,φ(0):=0.Equations (39) and (41) characterize the law of the excursion height of the random walk Xn. We note that φ(x) increases with x. From (38) and (41), we get an explicit expression of the probability of the event τx,h<τx,0:
(42)P(τx,h<τx,0)=P(Xτx=h)=φ(x)φ(h),
where τx,h is the first hitting time of h starting from x, with 0<x<h≤N. In particular, if h=N with {N} also absorbing, P(τx,N<τx,0)=φ(x)/φ(N) is the fixation probability starting from x, whereas P(τx,0<τx,N)=1-P(τx,N<τx,0) is the extinction probability.


## 3. A BD Example: The Moran Model

The example we have in mind is the2-allele Moran model with bias mechanism p. Let
(43)p(u):u∈[0,1]→[0,1],with0≤p(0),p(1)≤1
be continuous and p¯(u):=1-p(u).The Moran modelXn is characterized by the following transition probabilities, given the walker is in state x∈{0,…,N}:
(44)qx=xNp-(xN),rx=xNp(xN)+(1-xN)p-(xN),px=(1-xN)p(xN).For a Moran modelXn is the amount of say allele or type 1 individuals in a population of size N, and the reproduction law is obtained while choosing two individuals at random among N, one dying and the other one passing to the next generation while giving birth to an additional individual so as to keep constant the full population size N over the generations. Given then Xn=x, the transition Xn→Xn+1 (Xn→Xn-1) occurs with probability px (qx) translating the fact that the individual bound to die is chosen among type 2 (type 1) individuals with probability 1-x/N    (x/N), and the one which splits into two new born is chosen among the other type, with probability p(x/N) (resp., p-(x/N)), after deforming the neutral frequency x/N by the bias mechanism p. These bias p take into account various evolutionary forces which are the causes of deviation to neutrality (see [12–14] for a discussion on various models of utmost interest in population genetics, among which selection and mutations). The classical neutral Moran model is obtained when p(u)=u (no external deformation). Extended Moran models are studied in [15].Assumingp0=p(0)>0 and qN=1-p(1)>0, with y∈{1,…,N}, the corresponding Markov chain is ergodic with invariant distribution:
(45)πyπ0=∏z=0y-1pzqz+1=(Ny)∏z=0y-1p(z/N)p-((z+1)/N),
where
(46)π0=11+∑y=1N(Ny)∏z=0y-1(p(z/N)/p¯((z+1)/N)).We have the following obvious facts.Proposition 4.
LetXn be a Moran process with transition probabilities (44) for some bias p. Then, X-n:=N-Xn is again a Moran process with transition probabilities of the type (44) but with the new bias p~(u):=1-p(1-u), namely, q-x=pN-x=(x/N)(1-p~(x/N)), p-x=qN-x=(1-x/N)p~(x/N). The spectrum of the transition matrix P- of X-n is the same as the one P of Xn. If the bias is such that p~(u)=p(u), then the transition matrix of X-n coincides with the one of Xn.Proof.
We haveP-(x,y)=P(N-x,N-y), and P- is obtained from P applying the permutation σ:(0,1,…,N)→(N,…,1,0) reversing the order of the rows and then of the columns.


### 3.1. Fixation Probability

Forcing the states{0,N} to be absorbing, the scale function reads
(47)φ(0)=0,φ(x)=1+∑y=1x-11(N-1y)∏z=1yp¯(z/N)p(z/N),
with φ(x)/φ(N) being the fixation probability, whereas 1-φ(x)/φ(N) is the extinction probability, starting from x (for the neutral model with p(u)=u, the scale function is φ(x)=x, and the fixation probability is trivially x/N). The probability of a fixation before extinction occurs in an excursion reads
(48)P(H=N)=p0φ(N)=p(0)1+∑y=1N-1(1/(N-1y))∏z=1y(p¯(z/N)/p(z/N)).We have(49)P(H=N)=p0P(τ1,N<τ1,0),
where P(τ1,N<τ1,0)=1/φ(N) is the fixation probability of a mutant-type allele A1, for example, the probability that, starting from X0=1, the RW hits state {N} (fixation of the mutant-type allele 1) before hitting state {0} (extinction of allele 1).


### 3.2. Examples of Bias Mechanisms

We now give three examples of bias that we will deal with in the sequel.Moran Model with Mutations. A basic bias example is the mutation mechanism
(50)p(u)=(1-μ2)u+μ1(1-u),
where (μ1,μ2) are mutation probabilities in (0,1]. We  let  μ¯1:=1-μ1,   μ¯2:=1-μ2, μ=μ1+μ2, and  μ¯=1-μ.Schematically, type1 alleles are being produced whenever type 1 alleles do not mutate to type 2 (with probability μ¯2) or type 2 alleles mutate to type 1 (with probability μ1)(51)↻μ¯2A1⇌μ2μ1A2↺μ¯1.The transition matrixP of this model has therefore the particular BD structure
(52)qx=xN(μ¯1-μ¯xN),rx=1-(px+qx),px=(1-xN)(μ1+μ¯xN),
with p0=μ1, qN=μ2. When μ≠1, we introduce the auxiliary parameters
(53)ν1=μ1μ¯,ν2=μ2μ¯,ν=ν1+ν2=μμ¯,
with ν1+1=μ¯2/μ¯  andν2+1=μ¯1/μ¯.The drift of this RW is(54)E(Xn+1-Xn∣Xn=x)=px-qx=p(xN)-xN=μ1-μxN.The deterministic part of the dynamics ofXn has a stable point at x=Nμ1/μ, cancelling the drift. When p is nondecreasing, we have μ≤1. We note that in the expression of p~(u), the roles of μ1 and μ2 are interchanged.   So, p~(u)=p(u) if and only if μ1=μ2 (symmetric mutations).Ifμ1=μ2=1, we recover if N is even the heat-exchange Bernoulli-Laplace model [9] as a borderline example, but p(u)=1-u is strictly decreasing in this case. If μ1=μ2=1/2, then p(u)=1/2, and we obtain an aperiodic model amenable (through a suitable time substitution) to the Ehrenfest urn model provided N is even.For one-way mutations,(μ1,μ2)=(μ1,0) or (0,μ2) leading to the simpler forms μ1+uμ¯1 and μ¯2uofp(u), but with p(1)=1 and p(0)=0, respectively, corresponding to {N} (resp., {0}) being absorbing. We avoid this situation because it deserves another treatment due to transience.Moran Model with Mutations and Selection. Another example could be
(55)p(u)=μ1+u((1+s)μ¯2-μ1)1+su,
which is the composition of the selection bias mechanism
(56)p(u)=(1+s)u1+su(s>-1)
with the mutation mechanism (50). When considering the fixation probability P(H=N) for the Moran model with this particular mechanism, we obtain a closed-form formula for the fixation probability of a mutant which of course is an interesting issue in genetics.Cases with Positive Eigenvalues. We may look for conditions on the mechanism p leading to rx≥1/2 in which case the RW is spectrally positive. Assume p(u):u∈(0,1)→(0,1) is continuous and nondecreasing, with 0<p(0)≤p(1)<1.Then, as can easily be checked:rx≥1/2 for all x if and only if p(1/2)=1/2 with p(0)≤1/2≤p(1). Indeed, imposing rx≥1/2 for all x leads to p(u)≥1/2 if u≥1/2 and p(u)≤1/2 if u≤1/2. So, if p is nondecreasing with p(1/2)=1/2, then rx≥1/2.(i)
Take, for example, the mutation mechanism now satisfying0<μ1≤μ¯2<1. The condition p(1/2)=1/2 leads to μ1=μ2 (fair mutation). In this case, p(u)=μ1+(1-2μ1)u with 0<μ1≤1/2. If μ1=1/2, then p(u)=1/2 and qx=x/2N, rx=1/2, and px=(1/2)(1-x/N).For example, ifN=2,
(57)P=[μ¯1μ101412140μ1μ¯1]
has three distinct real nonnegative eigenvalues {1,μ¯1,(1-2μ1)/2} ordered in decreasing size. We note however that it is not necessary that μ1=μ2 for P to be spectrally positive; indeed, for the full Moran model with mutations μ1≠μ2, under the sole condition μ1+μ2≤1 (p is nondecreasing), the eigenvalues tk can be seen to be nonnegative. Still, we were not able to find a sharp general condition on p that would lead to a spectrally nonnegative Moran chain.(ii)
For the mechanismp of (55) combining mutations with selection, assuming 0<μ1≤μ¯2<1 to guarantee that p is nondecreasing, the condition p(1/2)=1/2 relates (s,μ1,μ2) to give s=(μ2-μ1)/(1/2-μ2), where both μ1,μ2<1/2. In this case, the model is spectrally positive.


### 3.3. Fixation Probabilities for the Moran Model with Mutations and/or Selection

We now study the particular Moran model with mutations and/or selection, starting with fixation probabilities.Moran with Mutations (sharp threshold phenomenon). We need to evaluate φ(N)=1+∑y=1N-1((N-1y)1/(N-1y))∏z=1y(p¯(z/N)/p(z/N)) for large N and then P(τ1,N<τ1,0)=1/φ(N) will follow.(i) Suppose first thatμ=1 and let ρ:=μ¯1/μ1.Proposition 5.
The fixation probability satisfies the sharp threshold property:(i)
ifμ1=1/2 (ρ=1), thenP(τ1,N<τ1,0)=1/φ(N)→1/2;(ii)
ifμ1<1/2 (ρ>1), thenP(τ1,N<τ1,0)=1/φ(N)~ρ-N→0;(iii)
ifμ1>1/2 (ρ<1), thenP(τ1,N<τ1,0)=1/φ(N)→1.Proof.
Withρ=μ¯1/μ1, we first have the identity for sums of inverse binomial coefficients (see e.g., [16, page  197])
(58)φ(N)∶=∑y=0N-1ρy(N-1y)=N(ρ+1)(1/ρ+1)N∑k=1N(ρk+1)(1/ρ+1)kk.
There is therefore a closed-form expression of the fixation probabilityP(τ1,N<τ1,0)=1/φ(N). Observing now
(59)φ(N)∶=∑y=0N-1ρy(N-1y)=1+ρN-1+ρN-2N-1+ρN-1+∑y=2N-3ρy(N-1y),
we get(i)
μ¯1=μ1  (ρ=1): φ(N)→2+. So: P(τ1,N<τ1,0)=1/φ(N)→1/2;(ii)
μ¯1>μ1 (ρ>1): φ(N)~ρN-1. So: P(τ1,N<τ1,0)=1/φ(N)~ρ-N→0;(iii)
μ¯1<μ1 (ρ<1): φ(N)→1+. So: P(τ1,N<τ1,0)=1/φ(N)→1.
Ifρ=1, because 1/(Ny)≤1/(N2) if y∈{2,N-2}: 2+2/N<φ(N+1)<2+2/N+2(N-3)/(N(N-1)) and φ1(N+1) is very close to 2+2/N. The fixation probability tends to 1/2 from below.
Ifρ≠1, 1+ρN+ρ/N+ρN-1/N<φ(N+1)<1+ρN + ρ/N+ρN-1/N + 2/(N(N-1))ρ2((1-ρN-3)/(1-ρ)), and φ(N+1) is very close to ρN if ρ>1 and very close to 1 up to a ρ/N factor if ρ<1.
Whenμ1<1/2, the probability, starting from {0} to hit the state {N} (fixation) starting from 1 before returning to {0}, tends to 0 geometrically fast. On the contrary, would μ1>1/2, this probability would converge to 1; mutations producing type 1 individuals are sufficiently large to enhance fixation of a mutant allele with probability 1. If μ1=1/2, the fixation probability approaches 1/2.(ii)
Suppose now thatμ≠1. With AN=N(μ¯1/μ¯)=N(ν2+1), BN=N(μ1/μ¯)=Nν1,
(60)φ(N)=1+∑y=1N-11(N-1y)∏z=1yμ¯1-μ¯(z/N)μ1+μ¯(z/N)=1+∑y=1N-11(N-1y)∏z=1yN(μ¯1/μ¯)-zN(μ1/μ¯)+z=1+∑y=1N-11(N-1y)Γ(AN)Γ(AN-y)Γ(BN+1)Γ(BN+1+y)=:1+∑y=1N-11(N-1y)ψN(y).Proceeding as before, we have
(61)1+ψN(1)N-1+ψN(N-1)<φ(N)<1+ψN(1)N-1+ψN(N-1)+O(N-2),
showing that φ(N)~Nlarge1+ψN(1)/(N-1)+ψN(N-1). The dominating term in φ(N) is thus either the one corresponding to y=1 or the one corresponding to y=N-1.Withr:=((μ¯1μ¯1/μ¯2μ¯2)(μ1μ1/μ2μ2))1/μ¯, by Stirling formula, we get
(62)φ(N)~Nlarge1+1N-1AN-1BN+1+Γ(AN)Γ(AN-N+1)Γ(BN+1)Γ(BN+N)~1+1N-1ν2+1ν1+Γ(N(ν2+1))Γ(Nν2+1)Γ(Nν1+1)Γ(N(ν1+1))=1+1N-1ν2+1ν1+ν1+1ν2+1(N(ν2+1)Nν2)(N(ν1+1)Nν1)~Nlarge1+1N-1ν2+1ν1+ν1+1ν2+1((ν2+1)ν2+1ν2ν2ν1ν1(ν1+1)ν1+1)N=:1+1N-1ν2+1ν1+ν1+1ν2+1rN.Now, wheneverμ1<μ2, regardless of whether μ>1 or μ<1,
(63)r∶=(ν2+1)ν2+1ν2ν2ν1ν1(ν1+1)ν1+1=(μ¯1μ¯1μ¯2μ¯2μ1μ1μ2μ2)1/μ¯>1.Indeed, letH(μ1)=-(μ1logμ1+μ¯1logμ¯1). The condition r>1 is also
(64)1μ¯(H(μ2)-H(μ1))>0.Ifμ1<μ2, then H(μ2)-H(μ1)<0 if and only if μ2>μ¯1 (or μ¯<0) and then (1/μ¯)(H(μ2)-H(μ1))>0.Ifμ1<μ2, then H(μ2)-H(μ1)>0 if and only if μ2<μ¯1 (or μ¯>0) and then (1/μ¯)(H(μ2)-H(μ1))>0.We obtained the following proposition.Proposition 6.
When
μ
≠
1,
the fixation probability also satisfies the following sharp threshold property.(i)
Ifμ1<μ2, the fixation probability tends to0geometrically fast becauseφ(N)~Nlarge(ν1+1)/(ν2+1)rN=μ¯2/μ¯1rN, leading toP(τ1,N<τ1,0)=1/φ(N)~Nlargeμ¯1/μ¯2r-N. The fixation probability of a mutant tends to0,algebraically fast at ratergiven by(63).(ii)
Ifμ1>μ2, regardless of whetherμ>1orμ<1,r<1. Thus,φ(N)~Nlarge1+(1/N)((ν2+1)/ν1)=1+(1/N)(μ¯1/μ1)leading toP(τ1,N<τ1,0)=1/φ(N)~Nlarge1/(1+(1/N)(μ¯1/μ1))→1from below.(iii)
Ifμ1=μ2(balanced mutations),ν1=ν2,r=1,
and therefore,(N)~Nlarge2+(1/N)(μ¯1/μ1)→2+: the fixation probability tends to1/2from below.(i) For theMoran model with selection mechanism p(u)=(1+s)u/(1+su) as in (56), gene A1 (resp., A2) has fitness w1=1+s (resp., w2=1) so that, with w:=w1/w2=1+s, the mutant relative fitness
(65)φ(x)=1+∑y=1x-1‍1(N-1y)∏z=1y‍p¯(z/N)p(z/N)=1+∑y=1x-1‍w-y=1-w-x1-w-1,
showing that:P(τ1,N<τ1,0)=1/φ(N)=(1-w-1)/(1-w-N). Whenever s>0 (s<0), the mutant gene A1 is selectively advantageous (deleterious) compared to A2. When w-1<1 (s>0), P(τ1,N<τ1,0) tends to 1-w-1, whereas when w-1>1 (s<0), this probability tends to 0 exponentially fast like (w-1-1)wN. When s=0 (neutral case), P(τ1,N<τ1,0)=1/N tends to 0 much slower.In this pure selection case, these results are well known (see, e.g., (3.66) page 109 of [14] and also [1] in [17]).(ii) For theMoran model with bias mechanism (55) involving both mutation and selection,
(66)φ(x)=1+∑y=1x-11(N-1y)∏z=1yp¯(z/N)p(z/N)=1+∑y=1x-11(N-1y)ψN(y),
where, with μ~=(1+s)μ¯2-μ1=μ¯+sμ¯2,AN=N(μ¯1/μ~), BN=N(μ1/μ~), ψN(y)=(Γ(AN)/Γ(AN-y))(Γ(BN+1)/Γ(BN+1+y)).From the preceding analysis, we conclude that with(67)r∶=(μ¯1μ¯1(μ¯2(1+s))μ¯2(1+s)μ1μ1(1-μ¯2(1+s))1-μ¯2(1+s))1/[μ¯+sμ¯2],φ(N)~Nlarge1+1N-1μ¯1μ1+(1+s)μ¯2μ¯1rN.We thus obtained the following proposition.Proposition 7.
For the Moran model with bias mechanism (55):(i)
ifw=1+s<μ¯1/μ¯2, the fixation probability tends to0  geometrically fast becauseφ(N)~Nlargew(μ¯2/μ¯1)rN,   leading toP(τ1,N<τ1,0)=1/φ(N)~Nlarge(1/w(μ¯2/μ¯1))r-N.   The fixation probability of a mutant tends to0, algebraically fast at rater  now given by (67);(ii)
if  w>μ¯1/μ¯2,  φ(N)~Nlarge1+(1/N)(μ¯1/μ1)  leading toP(τ1,N<τ1,0)=1/φ(N)~Nlarge1/(1+(1/N)(μ¯1/μ1))→1  from below;(iii)
in the critical regime when  w=μ¯1/μ¯2  and then  φ(N)~Nlarge2+(1/N)(μ¯1/μ1), the fixation probability tends again to  1/2.Proof.
Regardless of whetherμ<1 or μ>1, the extinction condition r>1 is achieved whenever
(68)1μ¯+sμ¯2(H(1-μ¯2(1+s))-H(μ1))>0.
By the above arguments, this leads to the following mutation-selection balance condition:μ1<1-μ¯2(1+s) which is also
(69)w∶=1+s<μ¯1μ¯2.
The selective advantage of alleleA1 (its relative fitness) should not be too large, and the upper bound is given by the ratio of the mutation probabilities μ¯1/μ¯2. Note that μ1<μ2 (resp., μ1>μ2) ⇒μ¯1/μ¯2>1 (resp., μ¯1/μ¯2<1). In the presence of selection, the extinction condition can be achieved even if μ1>μ2, provided that s<0 is small enough. When μ1=μ2, the extinction of the mutant allele occurs if it is selectively deleterious: s<0.To the best of the author’s knowledge, the results displayed in Propositions5–7 seem to be new.


### 3.4. Spectral Representation and Invariant Measure of the Moran Model with Mutations

Except for some exceptional special cases, the spectral measure associated to the Moran model with general bias mechanismp as in (44) is not known. One of the notable exceptions is the Moran model with positive mutation probabilities μ1,μ2; see [18–20]. The transition matrix P of this model has therefore the BD structure (52).For the Moran model with mutations, the orthogonal polynomials coincide with the dual Hahn polynomialsRx(λk):=Hk(x), where, with λk=k(k+Nν-1),
(70)Hk(x)=3F2((-k,k+Nν-1,-xNν1,-N);1),k=0,…,N
are the classical Hahn or Hahn-Eberlein polynomials. We note that (Hahn case) μ<1⇒ν1,ν2>0 with ν1,ν2→+∞ (μ→1-) and ν1,ν2→0 (μ→0+), whereas (Hahn-Eberlein case) μ>1⇒ν1,ν2<-1 with ν1,ν2→-∞ (μ→1+) and ν1,ν2→-1 (μ→2-).The caseμ<1 (μ>1) corresponds to a weak (strong) mutation regime.As a truncated hypergeometric series, the polynomialsRx(λk) are of degree k in x. In this case, thanks to the 3-term recurrence satisfied by the Hk(x), the eigenvalues of the transition matrix (52) are
(71)tk=1-λkN(N+Nν)=1-kN(μ+k-1Nμ¯),
depending only on the total mutation pressure μ. In particular, 1-t1=μ/N is the spectral gap, and tN=μ¯/N is the smallest eigenvalue (with tN<0 if μ>1).(i)
In the weak mutation regimeμ<1, the RW is spectrally positive with
(72)0<tN=μ¯N<⋯<t1=1-μN<t0=1.(ii)
In the strong mutation regime1<μ<2, the RW is not spectrally positive with
(73)tN=μ¯N<⋯<0<⋯<t1=1-μN<t0=1.WhenN gets large, there is a fixed (independent of N) number of negative eigenvalues. Only when μ=2 do we have a number of negative eigenvalues of order N/2.Invariant Measure.Furthermore, with (a)x:=Γ(a+x)/Γ(a) the rising factorial of a, the invariant probability measure is the generalized bivariate hypergeometric distribution:
(74)πx=(Nx)(Nν1)x(Nν2)N-x(Nν)N,x=0,…,N
which may also be recast as
(75)πx=(-Nν1x)(-Nν2N-x)(-NνN),x=0,…,N,with the convention(-ax)={-a}x/x!, {-a}x=-a(-a-1)⋯(-a-x+1), the falling factorials of -a with {-a}x=(-1)x(a)x.The mean and variance ofπx are
(76)m∞:=Nν1ν=Nμ1μ,σ∞2:=N2ν1ν2ν21+νNν+1=N2Nμ+μ¯μ1μ2μ2.Upon scaling, the law ofU∞:=X∞/N is very much peaked around its mean μ1/μ, with density
(77)fN(u)=NπNu=eNV(u)∫01dueNV(u),u∈(0,1),
where V(u)=logW(u) and, by Stirling formula,
(78)W(u)=(μ<1)(ν1+u)ν1+uν1ν1uu(ν2+1-u)ν2+1-uν2ν2(1-u)1-u×νν(ν+1)ν+1,ifν1,ν2>0,(79)W(u)=(μ>1)(-ν1)-ν1uu(-ν1-u)-ν1-u×(-ν2)-ν2(1-u)1-u(-ν2-1+u)-ν2-1+u×(-ν-1)-ν-1(-ν)-νifν1,ν2<-1,(80)W(u)=(μ=1)μ1uμ¯11-uuu(1-u)1-u.We haveV′(u)=0 at u*=ν1/ν=μ1/μ, and V′′(u*)=(-ν2/(1+ν))(1/ν1+1/ν2)=-μ3/μ1μ2<0.A saddle point estimate gives(81)∫01dueNV(u)=eNV(u*)2πN|V′′(u*)|,andsofN(u)~N|V′′(u*)|2πe-N(V(u*)-V(u)).(i) When bothν1,ν2>0 (μ<1), the distribution (74) is the negative hypergeometric distribution which can be obtained as a beta(Nν1,Nν2) mixture of the binomial distribution (Nx)px(1-p)N-x with parameter p∈[0,1], namely,
(82)πx=(Nx)∫01px(1-p)N-xΓ(Nν)Γ(Nν1)Γ(Nν2)×pNν1-1(1-p)Nν2-1dp,
with (by Stirling formula)
(83)π0=Γ(Nν)Γ(Nν1)Γ(Nν2)∫01pNν1-1(1-p)N(ν2+1)-1dp=Γ(Nν)Γ(Nν2)Γ(N(ν2+1))Γ(N(ν+1))~ν2νν+1ν2+1((ν2+1)ν2+1ν2ν2νν(ν+1)ν+1)N=:μ2μμ¯1ρN,
displaying exponential decay (ρ=(μ¯1μ¯1μμ/μ2μ2)1/μ¯<1). Note that, with
(84)ρ¯:=(ν1+1)ν1+1ν1ν1νν(ν+1)ν+1=(μ¯2μ¯2μμμ1μ1)1/μ¯<1,
it holds that
(85)πN=Γ(Nν)Γ(Nν1)Γ(Nν2)∫01pN(ν1+1)-1(1-p)Nν2-1dp=Γ(Nν)Γ(Nν1)Γ(N(ν1+1))Γ(N(ν+1))~μ1μμ¯2ρ¯N.Note that bothρ and ρ¯<1 because, with h(x):=-xlogx, μ<1⇒h(μ)>h(μ2)-h(μ¯1) and h(μ)>h(μ1)-h(μ¯2).Furthermore, the rater (63) appearing in the fixation probability is r=ρ/ρ¯.The case of symmetric mutationsμ1=μ2 (with μ<1) is obtained while setting ν1=ν2=ν/2.(ii) When bothν1,ν2<-1 (2>μ>1), the distribution πx in (75) constitutes the classical (generalized) hypergeometric distribution arising in sampling without replacement.By Stirling formula, for allb>a>0 (exponential growth),
(86)(NbNa)~12π(b-a)N/b[bbaa(b-a)b-a]N.This shows that, withρ<1 and ρ¯<1,
(87)π0=(-Nν2N)(-NνN)~ν2(ν+1)ν(ν2+1)[(-ν2)-ν2(-ν2-1)-ν2-1(-ν-1)-ν-1(-ν)-ν]N=:μ2μμ¯1ρN,πN=(-Nν1N)(-NνN)~ν1(ν+1)ν(ν1+1)[(-ν1)-ν1(-ν1-1)-ν1-1(-ν-1)-ν-1(-ν)-ν]N=:μ1μμ¯2ρ¯N.Note that again(88)ρ=(μ¯1μ¯1μμμ2μ2)1/μ¯,ρ¯=(μ¯2μ¯2μμμ1μ1)1/μ¯,
with ρ and ρ¯<1 because, with h(x):=-xlogx, μ>1⇒h(μ)<h(μ2)-h(μ¯1) and h(μ)<h(μ1)-h(μ¯2).(iii) When bothν1,ν2=-1 (μ1=μ2=1), we get the standard hypergeometric distribution πx=(Nx)(NN-x)/(2NN), with π0=πN=1/(2NN)~πN4-N.The case of symmetric mutationsμ1=μ2 (with μ>1) is obtained while setting ν1=ν2=ν/2.(iv) Whenμ1=μ2=1 (Bernoulli-Laplace), the transition probabilities read qx=(x/N)2,  rx=2(x/N)(1-x/N),px=(1-x/N)2. Here, πx=(Nx)(NN-x)/(2NN) (the standard hypergeometric distribution),   μk=((2N+1-2k)/(2N+1-k))(Nk)/(2N-kN), and tk=1-(k/N2)(2N+1-k).The expected return time to{0} is 22N/πN, which is huge when N is large, whereas the expected return time to {N/2} is of order πN/2, much smaller.(v) On the critical lineμ1+μ2=1, p(u)=μ1 is constant and so the transition probabilities become affine functions of the state: qx=μ¯1(x/N),rx=μ¯1+(2μ1-1)(x/N),andpx=μ1(1-x/N). Here, πx=(Nx)μ1xμ¯1N-x, μk=(Nk)μ1kμ¯1N-k are binomial bin(N,μ1)-distributed and self-dual, and tk=1-k/N, independent of μ1 (this RW is spectrally nonnegative). The Hahn polynomials boil down to Hk(x)=2F1((-k,-x-N);μ1-1), k=0,…,N, a class of Krawtchouk polynomials. When μ1=1/2 (the lazy Ehrenfest urn), the holding probabilities are rx=1/2, and both πx and μk are symmetric bin(N,1/2)-distributed.Spectral Measure When
μ
≠
1. The spectral probability measure is
(89)μk=2k+Nν-1k+Nν-1(Nk)(Nν2)N(k+Nν)N(Nν1)k(Nν2)k,k=0,…,N,
and the orthogonality relation reads
(90)∑k=0NμkRx(λk)Ry(λk)=(Nν2)NN!x!(N-x)!(Nν1)x(Nν2)N-xδx,y=π0πxδx,y.Notice that, in our notations:𝔮x(tk)=Rx(N(N+Nν)(1-tk))=Rx(N2(1-tk)/(μ¯)). In particular,
(91)𝔮x(t1)=:H1(x)=1-νNν1x,𝔮x(t2)=:H2(x)=1-Ax+Bx2,
where
(92)A=Nν+1Nν12Nν1+2+ν2-ν1(N-1)(Nν1+1),B=(Nν+1)(Nν+2)N2(Nν1+1)(N-1)ν1.Evolution of the Conditional Mean and Variance of
X
n
given
X
0
=
x.   𝔮x(t1) and 𝔮x(t2) being right eigenvectors of P,
(93)ExH1(Xn)=t1nH1(x),ExH2(Xn)=t2nH2(x)
leading, by virtue of (25), to the following peoposition.Proposition 8.
For the Moran model with mutations(94)Ex(Xn)=m∞(1-t1nH1(x))→n→∞m∞=Nν1ν=Nμ1μ,σx2(Xn)=Am∞B(1-t1nH1(x))-m∞2(1-t1nH1(x))2-1B+1Bt2nH2(x)→n→∞σ∞2.This shows how the mean and variance ofXn vary, during the course of evolution, as a function of n and the initial condition X0=x before reaching the mean and variance of the limit law. Concerning the mean value
(95)Ex(Xn)-m∞=t1n(x-m∞)
and depending on x>m∞ or x<m∞, Ex(Xn) approaches monotonically the limiting value m∞  from above or from below. Note however that for all initial conditions x,
(96)|Ex(Xn-m∞)|isadecreasing(negentropic)functionofn.Let us now consider the variance. Noting thatA/B-2m∞=(μ1-μ2)(2-μ)/(μ-1), rearranging the terms,
(97)σx2(Xn)=σ∞2-m∞H1(x)(μ2-μ1)(2-μ)μ¯t1n-m∞2H1(x)2t12n+1BH2(x)t2n
with (μ2-μ1)(2-μ)/μ¯>0 if (μ2>μ1 and μ<1) or (μ2<μ1 and μ>1) and t2<t12<t1<1. The corrective term to the limiting variance is a linear combination of three terms with geometric decay at different rates. For large n, the discrepancy between σx2(Xn) and σ∞2 is indeed controlled by the slower t1n term. In general, the variance does not increase monotonically towards its limiting value, starting from σx2(X0=x)=0.Proposition 9.
If for instance starting fromx<N μ1/μ=m∞ (which is H1(x)>0), depending on (μ2-μ1)(2-μ)/μ¯>0 or not, σx2(Xn) approaches its limiting value σ∞2 from below or from above, respectively. When (μ2-μ1)(2-μ)/μ¯<0, there is a time n=n*(x) at which the fluctuation σx2(Xn) is maximum, and this maximum exceeds σ∞2. The conclusions should be reversed would x>N μ1/μ=m∞ (or H1(x)<0).It can also be shown [21, page  47] that, for m≥0,
(98)covx(Xn,Xn+m)=t1mσx2(Xn),
meaning a geometric asymptotic decay of the autocovariance also controlled by the second largest eigenvalue t1.Evolution of Heterozygosity. Let
(99)𝔥(Xn)=Xn(N-Xn)=:XnX¯n
denote the heterozygosity of the chain. This function is the largest when Xn=N/2, that is, when both type population sizes Xn and X¯n are the closest. The exact evolution of Ex𝔥(Xn) and σx2(𝔥(Xn)) can be derived in a similar way using (25). For the variance σx2(𝔥(Xn)), 𝔮x(t4)=:H4(x), as a degree-4 polynomial in x, is needed. It can be shown that Ex𝔥(Xn) tends monotonically to σ∞2μ(N-1)~N2(μ1μ2/μ2), geometrically fast at rate t1, as n→∞. The variance tends to 0 geometrically fast as n→∞, and it is maximal at some intermediate time n=n*(x).Note that, with𝔤(Xn)=Xn-X¯n=2Xn-N denoting the population gap,
(100)Ex𝔤(Xn)=2m∞(1-t1nH1(x))-N,
and this function, just like Ex(Xn)-m∞, is not always (i.e., independently of x) decreasing with n.Increase of the Conditional Entropy Given
X
0
=
x  (H-Theorem). Let P(Xn∣X0=x)=Pn(x,Xn), n≥1. For any initial distribution P(X0=x), with
(101)Hn∶=-ElogP(Xn∣X0)=-∑x=0NP(X0=x)ElogP(Xn∣X0=x),n≥1,H0=0
the Shannon conditional entropy of the chain, we have [22]
(102)0≤Hn+1-Hn≤Hn-Hn-1.The first inequality shows thatHn is nondecreasing. It converges to its limiting maximal value, which is the Shannon entropy of the invariant measure. The second inequality shows that Hn is n-concave: the speed of increase of entropy decreases as time passes by. In particular, letting Hn(x):=-ElogP (Xn∣X0=x), for all fixed starting point x,H0(x)=0 and
(103)0≤Hn+1(x)-Hn(x)≤Hn(x)-Hn-1(x).The increase of entropy is a manifestation of irreversibility apparently in conflict with the reversibility property of all ergodic BD chains and their recurrence of states.Empirical Average and Large Deviation [23].Given X0=x0, let
(104)MnU:=1nSnU:=1n∑m=1nU(Xm)
denote the empirical average of some bounded measurable output function U of the (ergodic) chain. Examples of interesting Us are U(x)=x in which case MnU is directly the empirical average of the Xs, but also U(x)=1(x=y) in which case SnU is the fraction of time spent by the walker at state y (its local time). With ex0′ the row vector which is 0 everywhere but at position x0 where it is 1 and 1 a column vector of 1s, we have
(105)E(e-βSnU∣X0=x0)=∑x1,…,xne-β∑m=1nU(xm)∏m=1nP(xm-1,xm)=∑x1,…,xn,∏m=1nP(xm-1,xm)e-βU(xm)=ex0′Pβn1,
where, with wβ a column weight-vector with x entry wβ(x)=e-βU(x), Pβ:=PDwβ. By Perron-Frobenius theorem, we have
(106)-1nlogE(e-βSnU∣X0=x0)→n→∞-logλ(β)=:p(β),
where λ(β) is the dominant principal eigenvalue of Pβ and p(β) a concave function of β∈ℝ, with p(0)=0. Thus,
(107)1nlogP(1n∑m=1nU(Xm)∈n→∞dα)→n→∞f(α)dα,
where f(α)=infβ(αβ-p(β))≤0 is the concave Legendre transform of the pressure function p(β), giving the large deviation rate function of MnU. Clearly,
(108)1n∑m=1nU(Xm)→n→∞a.s.α*:=E(U(X∞)),
which is the value where f(α) is maximum (β*=f′(α*)=0) equal to 0.


## 3.1. Fixation Probability

Forcing the states{0,N} to be absorbing, the scale function reads
(47)φ(0)=0,φ(x)=1+∑y=1x-11(N-1y)∏z=1yp¯(z/N)p(z/N),
with φ(x)/φ(N) being the fixation probability, whereas 1-φ(x)/φ(N) is the extinction probability, starting from x (for the neutral model with p(u)=u, the scale function is φ(x)=x, and the fixation probability is trivially x/N). The probability of a fixation before extinction occurs in an excursion reads
(48)P(H=N)=p0φ(N)=p(0)1+∑y=1N-1(1/(N-1y))∏z=1y(p¯(z/N)/p(z/N)).We have(49)P(H=N)=p0P(τ1,N<τ1,0),
where P(τ1,N<τ1,0)=1/φ(N) is the fixation probability of a mutant-type allele A1, for example, the probability that, starting from X0=1, the RW hits state {N} (fixation of the mutant-type allele 1) before hitting state {0} (extinction of allele 1).


## 3.2. Examples of Bias Mechanisms

We now give three examples of bias that we will deal with in the sequel.Moran Model with Mutations. A basic bias example is the mutation mechanism
(50)p(u)=(1-μ2)u+μ1(1-u),
where (μ1,μ2) are mutation probabilities in (0,1]. We  let  μ¯1:=1-μ1,   μ¯2:=1-μ2, μ=μ1+μ2, and  μ¯=1-μ.Schematically, type1 alleles are being produced whenever type 1 alleles do not mutate to type 2 (with probability μ¯2) or type 2 alleles mutate to type 1 (with probability μ1)(51)↻μ¯2A1⇌μ2μ1A2↺μ¯1.The transition matrixP of this model has therefore the particular BD structure
(52)qx=xN(μ¯1-μ¯xN),rx=1-(px+qx),px=(1-xN)(μ1+μ¯xN),
with p0=μ1, qN=μ2. When μ≠1, we introduce the auxiliary parameters
(53)ν1=μ1μ¯,ν2=μ2μ¯,ν=ν1+ν2=μμ¯,
with ν1+1=μ¯2/μ¯  andν2+1=μ¯1/μ¯.The drift of this RW is(54)E(Xn+1-Xn∣Xn=x)=px-qx=p(xN)-xN=μ1-μxN.The deterministic part of the dynamics ofXn has a stable point at x=Nμ1/μ, cancelling the drift. When p is nondecreasing, we have μ≤1. We note that in the expression of p~(u), the roles of μ1 and μ2 are interchanged.   So, p~(u)=p(u) if and only if μ1=μ2 (symmetric mutations).Ifμ1=μ2=1, we recover if N is even the heat-exchange Bernoulli-Laplace model [9] as a borderline example, but p(u)=1-u is strictly decreasing in this case. If μ1=μ2=1/2, then p(u)=1/2, and we obtain an aperiodic model amenable (through a suitable time substitution) to the Ehrenfest urn model provided N is even.For one-way mutations,(μ1,μ2)=(μ1,0) or (0,μ2) leading to the simpler forms μ1+uμ¯1 and μ¯2uofp(u), but with p(1)=1 and p(0)=0, respectively, corresponding to {N} (resp., {0}) being absorbing. We avoid this situation because it deserves another treatment due to transience.Moran Model with Mutations and Selection. Another example could be
(55)p(u)=μ1+u((1+s)μ¯2-μ1)1+su,
which is the composition of the selection bias mechanism
(56)p(u)=(1+s)u1+su(s>-1)
with the mutation mechanism (50). When considering the fixation probability P(H=N) for the Moran model with this particular mechanism, we obtain a closed-form formula for the fixation probability of a mutant which of course is an interesting issue in genetics.Cases with Positive Eigenvalues. We may look for conditions on the mechanism p leading to rx≥1/2 in which case the RW is spectrally positive. Assume p(u):u∈(0,1)→(0,1) is continuous and nondecreasing, with 0<p(0)≤p(1)<1.Then, as can easily be checked:rx≥1/2 for all x if and only if p(1/2)=1/2 with p(0)≤1/2≤p(1). Indeed, imposing rx≥1/2 for all x leads to p(u)≥1/2 if u≥1/2 and p(u)≤1/2 if u≤1/2. So, if p is nondecreasing with p(1/2)=1/2, then rx≥1/2.(i)
Take, for example, the mutation mechanism now satisfying0<μ1≤μ¯2<1. The condition p(1/2)=1/2 leads to μ1=μ2 (fair mutation). In this case, p(u)=μ1+(1-2μ1)u with 0<μ1≤1/2. If μ1=1/2, then p(u)=1/2 and qx=x/2N, rx=1/2, and px=(1/2)(1-x/N).For example, ifN=2,
(57)P=[μ¯1μ101412140μ1μ¯1]
has three distinct real nonnegative eigenvalues {1,μ¯1,(1-2μ1)/2} ordered in decreasing size. We note however that it is not necessary that μ1=μ2 for P to be spectrally positive; indeed, for the full Moran model with mutations μ1≠μ2, under the sole condition μ1+μ2≤1 (p is nondecreasing), the eigenvalues tk can be seen to be nonnegative. Still, we were not able to find a sharp general condition on p that would lead to a spectrally nonnegative Moran chain.(ii)
For the mechanismp of (55) combining mutations with selection, assuming 0<μ1≤μ¯2<1 to guarantee that p is nondecreasing, the condition p(1/2)=1/2 relates (s,μ1,μ2) to give s=(μ2-μ1)/(1/2-μ2), where both μ1,μ2<1/2. In this case, the model is spectrally positive.


## 3.3. Fixation Probabilities for the Moran Model with Mutations and/or Selection

We now study the particular Moran model with mutations and/or selection, starting with fixation probabilities.Moran with Mutations (sharp threshold phenomenon). We need to evaluate φ(N)=1+∑y=1N-1((N-1y)1/(N-1y))∏z=1y(p¯(z/N)/p(z/N)) for large N and then P(τ1,N<τ1,0)=1/φ(N) will follow.(i) Suppose first thatμ=1 and let ρ:=μ¯1/μ1.Proposition 5.
The fixation probability satisfies the sharp threshold property:(i)
ifμ1=1/2 (ρ=1), thenP(τ1,N<τ1,0)=1/φ(N)→1/2;(ii)
ifμ1<1/2 (ρ>1), thenP(τ1,N<τ1,0)=1/φ(N)~ρ-N→0;(iii)
ifμ1>1/2 (ρ<1), thenP(τ1,N<τ1,0)=1/φ(N)→1.Proof.
Withρ=μ¯1/μ1, we first have the identity for sums of inverse binomial coefficients (see e.g., [16, page  197])
(58)φ(N)∶=∑y=0N-1ρy(N-1y)=N(ρ+1)(1/ρ+1)N∑k=1N(ρk+1)(1/ρ+1)kk.
There is therefore a closed-form expression of the fixation probabilityP(τ1,N<τ1,0)=1/φ(N). Observing now
(59)φ(N)∶=∑y=0N-1ρy(N-1y)=1+ρN-1+ρN-2N-1+ρN-1+∑y=2N-3ρy(N-1y),
we get(i)
μ¯1=μ1  (ρ=1): φ(N)→2+. So: P(τ1,N<τ1,0)=1/φ(N)→1/2;(ii)
μ¯1>μ1 (ρ>1): φ(N)~ρN-1. So: P(τ1,N<τ1,0)=1/φ(N)~ρ-N→0;(iii)
μ¯1<μ1 (ρ<1): φ(N)→1+. So: P(τ1,N<τ1,0)=1/φ(N)→1.
Ifρ=1, because 1/(Ny)≤1/(N2) if y∈{2,N-2}: 2+2/N<φ(N+1)<2+2/N+2(N-3)/(N(N-1)) and φ1(N+1) is very close to 2+2/N. The fixation probability tends to 1/2 from below.
Ifρ≠1, 1+ρN+ρ/N+ρN-1/N<φ(N+1)<1+ρN + ρ/N+ρN-1/N + 2/(N(N-1))ρ2((1-ρN-3)/(1-ρ)), and φ(N+1) is very close to ρN if ρ>1 and very close to 1 up to a ρ/N factor if ρ<1.
Whenμ1<1/2, the probability, starting from {0} to hit the state {N} (fixation) starting from 1 before returning to {0}, tends to 0 geometrically fast. On the contrary, would μ1>1/2, this probability would converge to 1; mutations producing type 1 individuals are sufficiently large to enhance fixation of a mutant allele with probability 1. If μ1=1/2, the fixation probability approaches 1/2.(ii)
Suppose now thatμ≠1. With AN=N(μ¯1/μ¯)=N(ν2+1), BN=N(μ1/μ¯)=Nν1,
(60)φ(N)=1+∑y=1N-11(N-1y)∏z=1yμ¯1-μ¯(z/N)μ1+μ¯(z/N)=1+∑y=1N-11(N-1y)∏z=1yN(μ¯1/μ¯)-zN(μ1/μ¯)+z=1+∑y=1N-11(N-1y)Γ(AN)Γ(AN-y)Γ(BN+1)Γ(BN+1+y)=:1+∑y=1N-11(N-1y)ψN(y).Proceeding as before, we have
(61)1+ψN(1)N-1+ψN(N-1)<φ(N)<1+ψN(1)N-1+ψN(N-1)+O(N-2),
showing that φ(N)~Nlarge1+ψN(1)/(N-1)+ψN(N-1). The dominating term in φ(N) is thus either the one corresponding to y=1 or the one corresponding to y=N-1.Withr:=((μ¯1μ¯1/μ¯2μ¯2)(μ1μ1/μ2μ2))1/μ¯, by Stirling formula, we get
(62)φ(N)~Nlarge1+1N-1AN-1BN+1+Γ(AN)Γ(AN-N+1)Γ(BN+1)Γ(BN+N)~1+1N-1ν2+1ν1+Γ(N(ν2+1))Γ(Nν2+1)Γ(Nν1+1)Γ(N(ν1+1))=1+1N-1ν2+1ν1+ν1+1ν2+1(N(ν2+1)Nν2)(N(ν1+1)Nν1)~Nlarge1+1N-1ν2+1ν1+ν1+1ν2+1((ν2+1)ν2+1ν2ν2ν1ν1(ν1+1)ν1+1)N=:1+1N-1ν2+1ν1+ν1+1ν2+1rN.Now, wheneverμ1<μ2, regardless of whether μ>1 or μ<1,
(63)r∶=(ν2+1)ν2+1ν2ν2ν1ν1(ν1+1)ν1+1=(μ¯1μ¯1μ¯2μ¯2μ1μ1μ2μ2)1/μ¯>1.Indeed, letH(μ1)=-(μ1logμ1+μ¯1logμ¯1). The condition r>1 is also
(64)1μ¯(H(μ2)-H(μ1))>0.Ifμ1<μ2, then H(μ2)-H(μ1)<0 if and only if μ2>μ¯1 (or μ¯<0) and then (1/μ¯)(H(μ2)-H(μ1))>0.Ifμ1<μ2, then H(μ2)-H(μ1)>0 if and only if μ2<μ¯1 (or μ¯>0) and then (1/μ¯)(H(μ2)-H(μ1))>0.We obtained the following proposition.Proposition 6.
When
μ
≠
1,
the fixation probability also satisfies the following sharp threshold property.(i)
Ifμ1<μ2, the fixation probability tends to0geometrically fast becauseφ(N)~Nlarge(ν1+1)/(ν2+1)rN=μ¯2/μ¯1rN, leading toP(τ1,N<τ1,0)=1/φ(N)~Nlargeμ¯1/μ¯2r-N. The fixation probability of a mutant tends to0,algebraically fast at ratergiven by(63).(ii)
Ifμ1>μ2, regardless of whetherμ>1orμ<1,r<1. Thus,φ(N)~Nlarge1+(1/N)((ν2+1)/ν1)=1+(1/N)(μ¯1/μ1)leading toP(τ1,N<τ1,0)=1/φ(N)~Nlarge1/(1+(1/N)(μ¯1/μ1))→1from below.(iii)
Ifμ1=μ2(balanced mutations),ν1=ν2,r=1,
and therefore,(N)~Nlarge2+(1/N)(μ¯1/μ1)→2+: the fixation probability tends to1/2from below.(i) For theMoran model with selection mechanism p(u)=(1+s)u/(1+su) as in (56), gene A1 (resp., A2) has fitness w1=1+s (resp., w2=1) so that, with w:=w1/w2=1+s, the mutant relative fitness
(65)φ(x)=1+∑y=1x-1‍1(N-1y)∏z=1y‍p¯(z/N)p(z/N)=1+∑y=1x-1‍w-y=1-w-x1-w-1,
showing that:P(τ1,N<τ1,0)=1/φ(N)=(1-w-1)/(1-w-N). Whenever s>0 (s<0), the mutant gene A1 is selectively advantageous (deleterious) compared to A2. When w-1<1 (s>0), P(τ1,N<τ1,0) tends to 1-w-1, whereas when w-1>1 (s<0), this probability tends to 0 exponentially fast like (w-1-1)wN. When s=0 (neutral case), P(τ1,N<τ1,0)=1/N tends to 0 much slower.In this pure selection case, these results are well known (see, e.g., (3.66) page 109 of [14] and also [1] in [17]).(ii) For theMoran model with bias mechanism (55) involving both mutation and selection,
(66)φ(x)=1+∑y=1x-11(N-1y)∏z=1yp¯(z/N)p(z/N)=1+∑y=1x-11(N-1y)ψN(y),
where, with μ~=(1+s)μ¯2-μ1=μ¯+sμ¯2,AN=N(μ¯1/μ~), BN=N(μ1/μ~), ψN(y)=(Γ(AN)/Γ(AN-y))(Γ(BN+1)/Γ(BN+1+y)).From the preceding analysis, we conclude that with(67)r∶=(μ¯1μ¯1(μ¯2(1+s))μ¯2(1+s)μ1μ1(1-μ¯2(1+s))1-μ¯2(1+s))1/[μ¯+sμ¯2],φ(N)~Nlarge1+1N-1μ¯1μ1+(1+s)μ¯2μ¯1rN.We thus obtained the following proposition.Proposition 7.
For the Moran model with bias mechanism (55):(i)
ifw=1+s<μ¯1/μ¯2, the fixation probability tends to0  geometrically fast becauseφ(N)~Nlargew(μ¯2/μ¯1)rN,   leading toP(τ1,N<τ1,0)=1/φ(N)~Nlarge(1/w(μ¯2/μ¯1))r-N.   The fixation probability of a mutant tends to0, algebraically fast at rater  now given by (67);(ii)
if  w>μ¯1/μ¯2,  φ(N)~Nlarge1+(1/N)(μ¯1/μ1)  leading toP(τ1,N<τ1,0)=1/φ(N)~Nlarge1/(1+(1/N)(μ¯1/μ1))→1  from below;(iii)
in the critical regime when  w=μ¯1/μ¯2  and then  φ(N)~Nlarge2+(1/N)(μ¯1/μ1), the fixation probability tends again to  1/2.Proof.
Regardless of whetherμ<1 or μ>1, the extinction condition r>1 is achieved whenever
(68)1μ¯+sμ¯2(H(1-μ¯2(1+s))-H(μ1))>0.
By the above arguments, this leads to the following mutation-selection balance condition:μ1<1-μ¯2(1+s) which is also
(69)w∶=1+s<μ¯1μ¯2.
The selective advantage of alleleA1 (its relative fitness) should not be too large, and the upper bound is given by the ratio of the mutation probabilities μ¯1/μ¯2. Note that μ1<μ2 (resp., μ1>μ2) ⇒μ¯1/μ¯2>1 (resp., μ¯1/μ¯2<1). In the presence of selection, the extinction condition can be achieved even if μ1>μ2, provided that s<0 is small enough. When μ1=μ2, the extinction of the mutant allele occurs if it is selectively deleterious: s<0.To the best of the author’s knowledge, the results displayed in Propositions5–7 seem to be new.


## 3.4. Spectral Representation and Invariant Measure of the Moran Model with Mutations

Except for some exceptional special cases, the spectral measure associated to the Moran model with general bias mechanismp as in (44) is not known. One of the notable exceptions is the Moran model with positive mutation probabilities μ1,μ2; see [18–20]. The transition matrix P of this model has therefore the BD structure (52).For the Moran model with mutations, the orthogonal polynomials coincide with the dual Hahn polynomialsRx(λk):=Hk(x), where, with λk=k(k+Nν-1),
(70)Hk(x)=3F2((-k,k+Nν-1,-xNν1,-N);1),k=0,…,N
are the classical Hahn or Hahn-Eberlein polynomials. We note that (Hahn case) μ<1⇒ν1,ν2>0 with ν1,ν2→+∞ (μ→1-) and ν1,ν2→0 (μ→0+), whereas (Hahn-Eberlein case) μ>1⇒ν1,ν2<-1 with ν1,ν2→-∞ (μ→1+) and ν1,ν2→-1 (μ→2-).The caseμ<1 (μ>1) corresponds to a weak (strong) mutation regime.As a truncated hypergeometric series, the polynomialsRx(λk) are of degree k in x. In this case, thanks to the 3-term recurrence satisfied by the Hk(x), the eigenvalues of the transition matrix (52) are
(71)tk=1-λkN(N+Nν)=1-kN(μ+k-1Nμ¯),
depending only on the total mutation pressure μ. In particular, 1-t1=μ/N is the spectral gap, and tN=μ¯/N is the smallest eigenvalue (with tN<0 if μ>1).(i)
In the weak mutation regimeμ<1, the RW is spectrally positive with
(72)0<tN=μ¯N<⋯<t1=1-μN<t0=1.(ii)
In the strong mutation regime1<μ<2, the RW is not spectrally positive with
(73)tN=μ¯N<⋯<0<⋯<t1=1-μN<t0=1.WhenN gets large, there is a fixed (independent of N) number of negative eigenvalues. Only when μ=2 do we have a number of negative eigenvalues of order N/2.Invariant Measure.Furthermore, with (a)x:=Γ(a+x)/Γ(a) the rising factorial of a, the invariant probability measure is the generalized bivariate hypergeometric distribution:
(74)πx=(Nx)(Nν1)x(Nν2)N-x(Nν)N,x=0,…,N
which may also be recast as
(75)πx=(-Nν1x)(-Nν2N-x)(-NνN),x=0,…,N,with the convention(-ax)={-a}x/x!, {-a}x=-a(-a-1)⋯(-a-x+1), the falling factorials of -a with {-a}x=(-1)x(a)x.The mean and variance ofπx are
(76)m∞:=Nν1ν=Nμ1μ,σ∞2:=N2ν1ν2ν21+νNν+1=N2Nμ+μ¯μ1μ2μ2.Upon scaling, the law ofU∞:=X∞/N is very much peaked around its mean μ1/μ, with density
(77)fN(u)=NπNu=eNV(u)∫01dueNV(u),u∈(0,1),
where V(u)=logW(u) and, by Stirling formula,
(78)W(u)=(μ<1)(ν1+u)ν1+uν1ν1uu(ν2+1-u)ν2+1-uν2ν2(1-u)1-u×νν(ν+1)ν+1,ifν1,ν2>0,(79)W(u)=(μ>1)(-ν1)-ν1uu(-ν1-u)-ν1-u×(-ν2)-ν2(1-u)1-u(-ν2-1+u)-ν2-1+u×(-ν-1)-ν-1(-ν)-νifν1,ν2<-1,(80)W(u)=(μ=1)μ1uμ¯11-uuu(1-u)1-u.We haveV′(u)=0 at u*=ν1/ν=μ1/μ, and V′′(u*)=(-ν2/(1+ν))(1/ν1+1/ν2)=-μ3/μ1μ2<0.A saddle point estimate gives(81)∫01dueNV(u)=eNV(u*)2πN|V′′(u*)|,andsofN(u)~N|V′′(u*)|2πe-N(V(u*)-V(u)).(i) When bothν1,ν2>0 (μ<1), the distribution (74) is the negative hypergeometric distribution which can be obtained as a beta(Nν1,Nν2) mixture of the binomial distribution (Nx)px(1-p)N-x with parameter p∈[0,1], namely,
(82)πx=(Nx)∫01px(1-p)N-xΓ(Nν)Γ(Nν1)Γ(Nν2)×pNν1-1(1-p)Nν2-1dp,
with (by Stirling formula)
(83)π0=Γ(Nν)Γ(Nν1)Γ(Nν2)∫01pNν1-1(1-p)N(ν2+1)-1dp=Γ(Nν)Γ(Nν2)Γ(N(ν2+1))Γ(N(ν+1))~ν2νν+1ν2+1((ν2+1)ν2+1ν2ν2νν(ν+1)ν+1)N=:μ2μμ¯1ρN,
displaying exponential decay (ρ=(μ¯1μ¯1μμ/μ2μ2)1/μ¯<1). Note that, with
(84)ρ¯:=(ν1+1)ν1+1ν1ν1νν(ν+1)ν+1=(μ¯2μ¯2μμμ1μ1)1/μ¯<1,
it holds that
(85)πN=Γ(Nν)Γ(Nν1)Γ(Nν2)∫01pN(ν1+1)-1(1-p)Nν2-1dp=Γ(Nν)Γ(Nν1)Γ(N(ν1+1))Γ(N(ν+1))~μ1μμ¯2ρ¯N.Note that bothρ and ρ¯<1 because, with h(x):=-xlogx, μ<1⇒h(μ)>h(μ2)-h(μ¯1) and h(μ)>h(μ1)-h(μ¯2).Furthermore, the rater (63) appearing in the fixation probability is r=ρ/ρ¯.The case of symmetric mutationsμ1=μ2 (with μ<1) is obtained while setting ν1=ν2=ν/2.(ii) When bothν1,ν2<-1 (2>μ>1), the distribution πx in (75) constitutes the classical (generalized) hypergeometric distribution arising in sampling without replacement.By Stirling formula, for allb>a>0 (exponential growth),
(86)(NbNa)~12π(b-a)N/b[bbaa(b-a)b-a]N.This shows that, withρ<1 and ρ¯<1,
(87)π0=(-Nν2N)(-NνN)~ν2(ν+1)ν(ν2+1)[(-ν2)-ν2(-ν2-1)-ν2-1(-ν-1)-ν-1(-ν)-ν]N=:μ2μμ¯1ρN,πN=(-Nν1N)(-NνN)~ν1(ν+1)ν(ν1+1)[(-ν1)-ν1(-ν1-1)-ν1-1(-ν-1)-ν-1(-ν)-ν]N=:μ1μμ¯2ρ¯N.Note that again(88)ρ=(μ¯1μ¯1μμμ2μ2)1/μ¯,ρ¯=(μ¯2μ¯2μμμ1μ1)1/μ¯,
with ρ and ρ¯<1 because, with h(x):=-xlogx, μ>1⇒h(μ)<h(μ2)-h(μ¯1) and h(μ)<h(μ1)-h(μ¯2).(iii) When bothν1,ν2=-1 (μ1=μ2=1), we get the standard hypergeometric distribution πx=(Nx)(NN-x)/(2NN), with π0=πN=1/(2NN)~πN4-N.The case of symmetric mutationsμ1=μ2 (with μ>1) is obtained while setting ν1=ν2=ν/2.(iv) Whenμ1=μ2=1 (Bernoulli-Laplace), the transition probabilities read qx=(x/N)2,  rx=2(x/N)(1-x/N),px=(1-x/N)2. Here, πx=(Nx)(NN-x)/(2NN) (the standard hypergeometric distribution),   μk=((2N+1-2k)/(2N+1-k))(Nk)/(2N-kN), and tk=1-(k/N2)(2N+1-k).The expected return time to{0} is 22N/πN, which is huge when N is large, whereas the expected return time to {N/2} is of order πN/2, much smaller.(v) On the critical lineμ1+μ2=1, p(u)=μ1 is constant and so the transition probabilities become affine functions of the state: qx=μ¯1(x/N),rx=μ¯1+(2μ1-1)(x/N),andpx=μ1(1-x/N). Here, πx=(Nx)μ1xμ¯1N-x, μk=(Nk)μ1kμ¯1N-k are binomial bin(N,μ1)-distributed and self-dual, and tk=1-k/N, independent of μ1 (this RW is spectrally nonnegative). The Hahn polynomials boil down to Hk(x)=2F1((-k,-x-N);μ1-1), k=0,…,N, a class of Krawtchouk polynomials. When μ1=1/2 (the lazy Ehrenfest urn), the holding probabilities are rx=1/2, and both πx and μk are symmetric bin(N,1/2)-distributed.Spectral Measure When
μ
≠
1. The spectral probability measure is
(89)μk=2k+Nν-1k+Nν-1(Nk)(Nν2)N(k+Nν)N(Nν1)k(Nν2)k,k=0,…,N,
and the orthogonality relation reads
(90)∑k=0NμkRx(λk)Ry(λk)=(Nν2)NN!x!(N-x)!(Nν1)x(Nν2)N-xδx,y=π0πxδx,y.Notice that, in our notations:𝔮x(tk)=Rx(N(N+Nν)(1-tk))=Rx(N2(1-tk)/(μ¯)). In particular,
(91)𝔮x(t1)=:H1(x)=1-νNν1x,𝔮x(t2)=:H2(x)=1-Ax+Bx2,
where
(92)A=Nν+1Nν12Nν1+2+ν2-ν1(N-1)(Nν1+1),B=(Nν+1)(Nν+2)N2(Nν1+1)(N-1)ν1.Evolution of the Conditional Mean and Variance of
X
n
given
X
0
=
x.   𝔮x(t1) and 𝔮x(t2) being right eigenvectors of P,
(93)ExH1(Xn)=t1nH1(x),ExH2(Xn)=t2nH2(x)
leading, by virtue of (25), to the following peoposition.Proposition 8.
For the Moran model with mutations(94)Ex(Xn)=m∞(1-t1nH1(x))→n→∞m∞=Nν1ν=Nμ1μ,σx2(Xn)=Am∞B(1-t1nH1(x))-m∞2(1-t1nH1(x))2-1B+1Bt2nH2(x)→n→∞σ∞2.This shows how the mean and variance ofXn vary, during the course of evolution, as a function of n and the initial condition X0=x before reaching the mean and variance of the limit law. Concerning the mean value
(95)Ex(Xn)-m∞=t1n(x-m∞)
and depending on x>m∞ or x<m∞, Ex(Xn) approaches monotonically the limiting value m∞  from above or from below. Note however that for all initial conditions x,
(96)|Ex(Xn-m∞)|isadecreasing(negentropic)functionofn.Let us now consider the variance. Noting thatA/B-2m∞=(μ1-μ2)(2-μ)/(μ-1), rearranging the terms,
(97)σx2(Xn)=σ∞2-m∞H1(x)(μ2-μ1)(2-μ)μ¯t1n-m∞2H1(x)2t12n+1BH2(x)t2n
with (μ2-μ1)(2-μ)/μ¯>0 if (μ2>μ1 and μ<1) or (μ2<μ1 and μ>1) and t2<t12<t1<1. The corrective term to the limiting variance is a linear combination of three terms with geometric decay at different rates. For large n, the discrepancy between σx2(Xn) and σ∞2 is indeed controlled by the slower t1n term. In general, the variance does not increase monotonically towards its limiting value, starting from σx2(X0=x)=0.Proposition 9.
If for instance starting fromx<N μ1/μ=m∞ (which is H1(x)>0), depending on (μ2-μ1)(2-μ)/μ¯>0 or not, σx2(Xn) approaches its limiting value σ∞2 from below or from above, respectively. When (μ2-μ1)(2-μ)/μ¯<0, there is a time n=n*(x) at which the fluctuation σx2(Xn) is maximum, and this maximum exceeds σ∞2. The conclusions should be reversed would x>N μ1/μ=m∞ (or H1(x)<0).It can also be shown [21, page  47] that, for m≥0,
(98)covx(Xn,Xn+m)=t1mσx2(Xn),
meaning a geometric asymptotic decay of the autocovariance also controlled by the second largest eigenvalue t1.Evolution of Heterozygosity. Let
(99)𝔥(Xn)=Xn(N-Xn)=:XnX¯n
denote the heterozygosity of the chain. This function is the largest when Xn=N/2, that is, when both type population sizes Xn and X¯n are the closest. The exact evolution of Ex𝔥(Xn) and σx2(𝔥(Xn)) can be derived in a similar way using (25). For the variance σx2(𝔥(Xn)), 𝔮x(t4)=:H4(x), as a degree-4 polynomial in x, is needed. It can be shown that Ex𝔥(Xn) tends monotonically to σ∞2μ(N-1)~N2(μ1μ2/μ2), geometrically fast at rate t1, as n→∞. The variance tends to 0 geometrically fast as n→∞, and it is maximal at some intermediate time n=n*(x).Note that, with𝔤(Xn)=Xn-X¯n=2Xn-N denoting the population gap,
(100)Ex𝔤(Xn)=2m∞(1-t1nH1(x))-N,
and this function, just like Ex(Xn)-m∞, is not always (i.e., independently of x) decreasing with n.Increase of the Conditional Entropy Given
X
0
=
x  (H-Theorem). Let P(Xn∣X0=x)=Pn(x,Xn), n≥1. For any initial distribution P(X0=x), with
(101)Hn∶=-ElogP(Xn∣X0)=-∑x=0NP(X0=x)ElogP(Xn∣X0=x),n≥1,H0=0
the Shannon conditional entropy of the chain, we have [22]
(102)0≤Hn+1-Hn≤Hn-Hn-1.The first inequality shows thatHn is nondecreasing. It converges to its limiting maximal value, which is the Shannon entropy of the invariant measure. The second inequality shows that Hn is n-concave: the speed of increase of entropy decreases as time passes by. In particular, letting Hn(x):=-ElogP (Xn∣X0=x), for all fixed starting point x,H0(x)=0 and
(103)0≤Hn+1(x)-Hn(x)≤Hn(x)-Hn-1(x).The increase of entropy is a manifestation of irreversibility apparently in conflict with the reversibility property of all ergodic BD chains and their recurrence of states.Empirical Average and Large Deviation [23].Given X0=x0, let
(104)MnU:=1nSnU:=1n∑m=1nU(Xm)
denote the empirical average of some bounded measurable output function U of the (ergodic) chain. Examples of interesting Us are U(x)=x in which case MnU is directly the empirical average of the Xs, but also U(x)=1(x=y) in which case SnU is the fraction of time spent by the walker at state y (its local time). With ex0′ the row vector which is 0 everywhere but at position x0 where it is 1 and 1 a column vector of 1s, we have
(105)E(e-βSnU∣X0=x0)=∑x1,…,xne-β∑m=1nU(xm)∏m=1nP(xm-1,xm)=∑x1,…,xn,∏m=1nP(xm-1,xm)e-βU(xm)=ex0′Pβn1,
where, with wβ a column weight-vector with x entry wβ(x)=e-βU(x), Pβ:=PDwβ. By Perron-Frobenius theorem, we have
(106)-1nlogE(e-βSnU∣X0=x0)→n→∞-logλ(β)=:p(β),
where λ(β) is the dominant principal eigenvalue of Pβ and p(β) a concave function of β∈ℝ, with p(0)=0. Thus,
(107)1nlogP(1n∑m=1nU(Xm)∈n→∞dα)→n→∞f(α)dα,
where f(α)=infβ(αβ-p(β))≤0 is the concave Legendre transform of the pressure function p(β), giving the large deviation rate function of MnU. Clearly,
(108)1n∑m=1nU(Xm)→n→∞a.s.α*:=E(U(X∞)),
which is the value where f(α) is maximum (β*=f′(α*)=0) equal to 0.


## 4. A Detailed Study of the Siegmund Dual of BD Chains with an Application to the Mutation Moran Model

In this Section, we will illustrate the power of the duality/intertwining relationship by considering the simplest Siegmund dual of a birth and death chain, with the mutation Moran birth and death example in mind. We also address here the question of computing the strong stationary time distribution that helps quantifying the “distance” to equilibrium of the original positive recurrent BD process.Definition 10 (see [24]).
Two discrete-time Markov processes(Xn,X^n;n≥0), with state spaces (𝒳,𝒴), possibly with substochastic transition kernels, are said to be dual with respect to some nonsingular duality kernel H≥0 on the product space (𝒳,𝒴) if ∀x∈𝒳,∀y∈𝒴,∀n∈ℕ:(109)ExH(Xn,y)=EyH(x,X^n).When state spaces(𝒳,𝒴)={0,…,N}2 are finite and identical, the duality kernel is a square-matrix, and the transition matrix of the dual process X^, say P,^ is obtained from the one P of the direct process by
(110)P^′=H-1PH,
where P^′ stands for the transpose of P^. Note that if P^ is an H-dual to P, then P is an H′-dual to P^. If H=H′, P^ is an H-dual to P, but also P is an H-dual to P^.The Siegmund Duality Kernel. (For nonneutral population genetics models, other duality kernels have recently been shown to be of interest, see [25, 26]). The Siegmund duality kernel is H(x,y)=1(x≤y). If, for a given process Xn, a process X^n exists satisfying the above condition, X^n is called the Siegmund dual of Xn, see [27]. Clearly, in the BD case for Xn, the condition for the Siegmund dual to exist is that Xn should be stochastically monotone in that, for all y≥0 and n≥0,Px(Xn>y) should be increasing with x.For positive recurrent birth and death processes, and for the Siegmund kernel, the transition matrixP^ of the dual process X^n reads
(111)P^=[r0-q1q1p1r^1q2⋱⋱⋱pN-1r^N-1qN01],
where r^y:=1-(py+qy+1),   y∈{1,…,N-1} (and q^y=py,   y=1,…,N-1,p^y=qy+1, y=0,…,N-1). It is again the one of a BD process (but not of a Moran BD process if P is a Moran transition matrix). Indeed, in this case, H is an upper-right triangular matrix with nonzero entries 1, whereas the nonnull entries of H-1 are the diagonal (with entries 1) and the upper diagonal with entries -1. The structure of P^ follows from this, and the duality relation: P^′=H-1PH.For this dual to exist, we need to ensure thatpy+qy+1≤1 for y∈{0,…,N-1} which is a necessary and sufficient condition to guarantee the stochastic monotonicity of Xn. This condition also reads: py+1-py≥-ry+1. The drift of X^n at y is
(112)f^(y)∶=p^y-q^y=qy+1-py=-f(y+1)+(py+1-py),
where f(y)∶=py-qy is the drift of Xn at y. We have
(113)-f(y+1)-ry+1≤f^(y)≤-f(y)+ry.We already know that ifP is a spectrally nonnegative BD matrix, the chain is stochastically monotone. Here is another sufficient condition relative to the specific ergodic Moran case.Proposition 11.
Consider the Moran modelXn with bias p satisfying 0<p(0), p(1)<1. If p(u) is nondecreasing, the condition px+qx+1≤1 is fulfilled (the Moran chain is stochastically monotone), and so the Siegmund dual exists.Proof.
We first need to guarantee that1-p0-q1≥0 corresponding to the row x=0 of P^, but this is true because p(1/N)≥p(0)≥Np(0)-(N-1). Next, when x=N, the inequality clearly holds because pN=qN+1=0. Finally, for x∈{1,…,N-1},
(114)px+qx+1=p(x/N)-xp(x/N)N+(x+1)N-(x+1)p((x+1)/N)N≤p(xN)(1-x/N-(x+1)N)+(x+1)N=[p(x/N)(N-2x-1)+x+1]N≤1
because
(115)p(xN)≤(1-x/(N-1))(1-2x/(N-1)),ifx∈{0,…,⌊N-12⌋}p(xN)≥(1-x/(N-1))(1-2x/(N-1)),ifx∈{⌈N-12⌉,…,N-1}.For the mutation Moran model in the weak mutation regime0<μ<1, p(u) is nondecreasing, and therefore, this chain is stochastically monotone. Although the condition that p(u) is nondecreasing is a sufficient condition to guarantee that the chain is stochastically monotone, it is not necessary.Proposition 12.
Consider the Moran modelXn with mutation bias p in the strong mutation regime 2>μ≥1. Although p(u) is not nondecreasing, the condition px+qx+1≤1 is still fulfilled would N be large enough (in which case the Moran chain with strong mutations is still stochastically monotone). And then, the Siegmund dual exists.Proof.
When2>μ>1, p(u)=μ1+μ¯u is monotone decreasing. We need to check that in that case alsopx+qx+1≤1 when x varies from 0 to N-1, as soon as N is large enough. After some elementary algebra,
(116)px+qx+1-1=-2μ¯(xN)2+(2μ¯1-μ-2μ¯N)(xN)-(μ¯1-μ¯1N+μ¯N2).
The minimum if this convex degree-2 polynomial is attained at the point x=N(2μ¯1-μ-2(μ¯/N))/(4μ¯) so, when N is large, to the left (right) of x=N/2 if μ2>μ1 (if μ1>μ2) with a negative value of this minimum. The maximum of px+qx+1-1 is thus attained either at x=N-1 (at x=0), with pN-1+qN-1=-μ¯2+(1/N)(μ¯2-μ¯/N) (resp., p0+q1-1=-μ¯1+(1/N)(μ¯1-μ¯/N)). These maxima are nonpositive for all N sufficiently large so that N2-N>max(-μ¯/μ¯2,-μ¯/μ¯1). Then, given a Moran model with 2>μ>1, for all N large enough, both pN-1+qN-1<0 and p0+q1-1<0.
Whenμ=1, p(u)=μ1 and px+qx+1 is a linear function of x. The maximum of px+qx+1-1 is thus attained either at x=N-1 (at x=0), with pN-1+qN-1=-(μ1/N)(N-1) (resp., p0+q1-1=-(μ¯1/N)(N-1)), both nonpositive. So px+qx+1≤1 for all x and the Siegmund dual also exists when μ=1, whatever the value of N.Remark 13.
Whenever0<μ<2, therefore, the Moran model with mutations is or can be made stochastically monotone. When μ=2, px+qx+1-1<0 except at least at both x=0 and x=N-1 where px+qx+1-1=1/N2>0. The Bernoulli-Laplace chain fails to be stochastically monotone.
From the structure ofP^, it is apparent that the dual process loses mass at y=0 and is absorbed at y=N. Let us therefore add a coffin state ∂:={-1} and consider the enlarged stochastic matrix which we will call P^∂:
(117)P^∂=[101-r0r0-q1q1p1r^1q2⋱⋱⋱pN-1r^N-1qN01].
The corresponding proper BD chain, call it∂X^n, now has two absorbing states, one at {-1} and one at {N}. Let now φ^(y), y=-1,0,1,…,N be the scale function of ∂X^n, solving P^∂φ^=φ^, forcing φ^(-1)=0. We can easily check that
(118)φ^(-1)=0,φ^(0)=1,φ^(y)=γyc∶=∑z=0yγz=1π0∑z=0yπz
so that the scale function of ∂X^n expresses in terms of the cumulative distribution πyc∶=∑z=0y‍πz of the invariant measure of the original process.
Letτ^y:=τ^y,-1∧τ^y,N be the infimum of the first hitting time of {-1} and {N} starting from y∈{0,…,N-1}. We have
(119)Py(∂X^τ^y=N)=φ^(y)φ^(N)=γycγNc=πyc.
Doob
h
-
Transform. Define a new transition matrix P~∂ by
(120)P~∂(x,y)=πycπxcP^∂(x,y),x,y∈{-1,0,…,N}2.
We have(121)P~∂=[100r0-q1p0+q1π0cπ1cp1r^1π2cπ1cq2⋱⋱⋱πN-2cπN-1cpN-1r^N-1πNcπN-1cqN01],
where state {-1} becomes isolated and disconnected. Deleting the row and column associated to the entry {-1} of P~∂, we get a stochastic matrix, call it P~, of a process X~n on {0,…,N} which corresponds to ∂X^n conditioned to first hit state {N} before state {-1}. The state {0} of this conditioned BD process now is partially reflecting whereas the remaining absorbing state, say a, is a={N}.
We will also need subsequently to introduceφ~(x), x=0,1,…,N, which is the new scale function of X~n. We have
(122)φ~(0)=0,φ^(1)=1,φ^(x)=1+∑y=1x-1∏z=1yq~zp~z=1+∑y=1x-1∏z=1yπz+1cπz-1cqz+1pz
because q~z=(πz+1c/πzc)qz+1, p~z=(πz-1c/πzc)pz and r~z=r^z,z=1,…,N-1, q~0=p0+q1 and r~0=r0-q1 and r~N=1 and p~N=0 are the new transition probabilities of X~n on {0,…,N} which can now be read from P~.
We now show thatP~ and P are intertwined through a stochastic link.Proposition 14.
(i) The matricesP~ and P are similar (with the same eigenvalues); that is,
(123)P~=ΛPΛ-1.
The linkΛ is given by Λ(x~,x)=(πx/πx~c)1(x≤x~), corresponding to the entries of a lower-triangular stochastic matrix. In other words, for all n≥0,
(124)Λ(x~,x)=P(Xn=x∣X~n=x~),
and πn′=π~n′Λ, where πn=Pπ0(Xn=·) and π~n=Pπ~0(X~n=·).
(ii) The linkΛ satisfies
(125)Λ(N,x)=πx,x=0,…,N.
(iii)π~0′=π0′=e0′:=(1,0,…,0) are admissible initial distributions of the chains X~n and Xn, satisfying
(126)π0′=π~0′Λ.
(iv)P~ is K-dual to P:(127)P~′=K-1PK,
where K(x,y)=(1/πyc)1 (x≤y).Proof.
(i) From (120),
(128)P~=Dπc-1P^Dπc,whereDπc=diag(π0c,…,πNc).
The random walk Xn being reversible, P′=DπPDπ-1. Thus,
(129)P~=Dπc-1H′P′H′-1Dπc=Dπc-1H′DπPDπ-1H′-1Dπc.
As a result, Λ=Dπc-1H′Dπ with the right entries Λ(x~,x)=(πx/πx~c)1  (x≤x~), satisfying ∑x≤x~(πx/πx~c)=1.
(ii) The last row ofΛ is given by Λ(N,x)=πx so that once X~n hits state N, the law of Xn is π.
(iii) The first row ofΛ is (1,0,…,0) so that e0′=e0′Λ.
(iv) Indeed,P~′=DπcP^′Dπc-1=DπcH-1PHDπc-1:=K-1PK with K=HDπc-1 with the right entries.SoP~ (as the composition of the Siegmund dual of P with a Doob-transform) can be obtained from P either from (i) through a stochastic link Λ or from (iv) through the duality kernel K. This is because BD chains such as the Moran model are reversible.Strong Stationary Time.The intertwining construction shows that the original positive recurrent BD chain Xn with transition matrix P may also be viewed as the output (through the link Λ) of an intertwined hidden Markov chain X~n with transition matrix P~. Once X~n hits its absorbing state {N}, the RW Xn is distributed like π, provided both Xn and X~n were started at 0. Furthermore, there exists a bivariate Markov chain (X~,X) with transition kernel
(130)P((x~,x),(y~,y))=P(x,y)·P~(x~,y~)·Λ(y~,y)(ΛP)(x~,y)1(ΛP)(x~,y)>0,
where x~∈{x,x±1},y~∈{y,y±1}. We have (ΛP)(x~,y)>0 if and only if y≤x~+1. We can check that ∑y‍∑y~‍P((x~,x),(y~,y))=1 as required for a stochastic matrix.Withx∈{0,…,N-1}, we will let
(131)τ~x~0,N=inf(n:X~n=N∣X~0=x~0)
be the first hitting time of {N} of X~n, starting from state x~0∈{0,…,N-1}. The random time τ~0,N gives some information on the speed of convergence of the law of the original process Xn to its invariant measure (is a strong stationary time in the sense of Diaconis and Fill [5]). The facts (123), (124), (125), and (126) indeed guarantee that τ~0,N is a strong stationary time of Xn in the sense that Xτ~0,N~dπ and is independent of τ~0,N (see [5, Theorems  2.4 and  2.17] or [6, Theorem  2.1]). Equivalently (see [4, Proposition3.2]), it holds that
(132)sep(πn,π)≤P(τ~0,N>n)≤E(τ~0,N)n,
where πn(·):=Pn(0,·) is the law of Xn started at 0, π its invariant measure. In (132), the separation discrepancy is defined by sep(πn,π)=supy[1-πn(y)/πy]. It satisfies sep(πn,π)≥∥πn-π∥TV, where ∥πn-π∥TV=(1/2)∑y‍|πn(y)-πy| is the total variation distance between πn and π.Furthermore, from (123), (126), there is a unique “witness” state say d=N such that either π~n(N)=0 or π~n(N)>0⇒πn(d)=π~n(N)πd>0 showing that this random time is stochastically the smallest since the first inequality in (132) turns out to be an equality (see [5, Remark  2.39] and Proposition 15 below).In the context of BD chains absorbed atN, the probability generating function (pgf) of τ~0,N≥N is [6, 28]:
(133)E(zτ~0,N)=∏k=1N(1-tk)z1-tkz,z∈[0,1],
where -1<tk<+1, k=1,…,N are the N distinct eigenvalues of both P~ and P, avoiding t0=1. The formula (133) also reads
(134)E(zτ~0,N)=∑l=1NAlz(1-tl)1-tlz,
where
(135)Al=∏k≠l1-tktl-tk,(136)P(τ~0,N>n)=∑l=1N‍∏k≠l‍1-tktl-tktln,n≥N-1.Thus,t1-nP(τ~0,N>n)→n↑∞∏k=2N‍((1-tk)/(t1-tk)), and τ~0,N has geometric tails with exponent t1. We have
(137)E(τ~0,N)=trace(I-P~)-1=∑k=1N(1-tk)-1,σ2(τ~0,N)=∑k=1Ntk(1-tk)-2=∑k=1N(1-tk)-2-∑k=1N(1-tk)-1.Note that, sincet1 is the dominant eigenvalue,
(138)σ2(τ~0,N)≤E(τ~0,N)1-t1.
When the eigenvalues tk are nonnegative, then τ~0,N=d∑k=1Nτk where the τks are independent with τk~d geom(1-tk) the geometric distribution with success parameter 1-tk on {1,2,…}. When the eigenvalues tk are not all positive, it is not obvious that the previous expression (133) of E(zτ~0,N) is indeed a pgf but it turns out that this is the case. Assuming tN<⋯<tl+1<0≤tl<⋯<t1<t0=1, (133) takes on the interpretation of
(139)τ~0,N-∑k=l+1Nbk=d∑k=1lτk,
where bk~d bernoulli (1/(1-tk)),τk~d geom(1-tk), and τ~0,N are all mutually independent. We can summarize these results that should mainly be attributed to [4–7] as follows.Proposition 15.
Suppose a Siegmund dual exists for a finite state space ergodic BD chainXn. Then, there exists a Markov chain X~n, intertwined with Xn, with {N} as an absorbing state and fully described in Proposition 14. The random time τ~0,N is a fastest strong stationary time for Xn whose law is characterized either by (133) or (136) involving the spectrum of either P or P~, the transition matrices governing the two processes.Remark 16.
When the RW is symmetric, except fortN=-1 and maybe {0}, the eigenvalues can be paired, because if tk∈(0,1) is an eigenvalue, then tN-k=-tk is one also. The pgf reads
(140)E(zτ~0,N)=2z21+z∏k=1N0-1(1-tk2)z21-tk2z2,ifN=2N0,E(zτ~0,N)=2z1+z∏k=1N0(1-tk2)z21-tk2z2,ifN=2N0+1,
depending on whether {0} is or is not an eigenvalue.
Computing the Mean and Variance of
τ
~
0
,
N
.Suppose that the tk are known explicitly. In this case, it is possible to compute E(τ~0,N) and σ2(τ~0,N) and find conditions under which
(141)E(τ~0,N)→∞,σ2(τ~0,NE(τ~0,N))→0,asN↑∞.
If this is the case, then τ~0,N/E(τ~0,N)→1 in probability and ⌊E(τ~0,N)/2⌋ is expected to be a cut-off time for Xn started at {0}.
We now apply this construction to the particular Moran model under study.
The Case of the Moran Model with Mutations. The eigenvalues tk are known, leading to: 1-tk=(k/N)(μ+μ¯((k-1)/N)). If we scale the characteristic time τ~0,N to N, the summations for large N can be replaced by integrals, leading to the estimate
(142)μN~N∫01dx(x+1/N)(μ+μ¯x)=N2Nμ-μ¯(∫01dxx+1/N-μ¯∫01dxμ+μ¯x)
and similarly for σ2(τ~0,N). From this, we easily get the following proposition.Proposition 17.
For the Moran model with mutations with0<μ<2, whenever the Siegmund dual exists,
(143)μN~Nμ(logN+logμ),σ2(τ~0,N)~(Nμ)2,
showing that σ2(τ~0,N/E(τ~0,N))~(logN)-2→0.We have(144)Al=∏k≠l‍1-tktl-tk=(∏k≠l‍(1-l(Nμ+(l-1)μ¯)k(Nμ+(k-1)μ¯)))-1,
and the tail distribution of τ~0,N is
(145)P(τ~0,N>n)=∑l=1N∏k≠l1-tktl-tktln,n≥N-1.Remark 18.
Whenμ=2, the Bernoulli-Laplace RW fails to be stochastically monotone, and so the above construction also fails as the chain has no Siegmund dual. However, as was shown in [29], similar conclusions on the fastest strong stationary time hold for the embedded continuous-time version of the Bernoulli-Laplace chain.Proposition 19.
We have the Gumbel limit law(146)τ~0,NScaled:=τ~0,N-(N/μ)logNN/μ→dX~de-(x+e-x),x∈ℝ.Proof.
Indeed, usingtk/(1-tk)~N/(kμ) to the dominant order in N,
(147)E(eitτ~0,NScaled)~N-it∏k=1N1-tk1-tk(1+itμ/N)~N-it∏k=1Nkk-it→Γ(1-it),asN↑∞,
which is the Fourier transform of e-(x+e-x), x∈ℝ.WithnN(θ)=⌊(N/2μ)(logN+θ)⌋, then
(148)∥PnN(θ)(0,·)-π∥TV→N↑∞c(θ),
where c(θ)→θ↑∞0 and c(θ)→θ↑-∞1. The expected mixing time is μN~(N/μ)logN, whereas the spectral gap is 1-t1=μ/N, the product of the 2 of which tends to ∞. Recalling σ2(τ~0,N)≤μN/(1-t1), σ2(τ~0,N/μN)=μN-2σ2(τ~0,N)≤1/((1-t1)μN), the condition (1-t1)μN→∞ is a sufficient condition for σ2(τ~0,N/μN)→0. If this holds, the contribution of ∑k=2N‍(1-tk)-1 to μN dominates the lead term (1-t1)-1 (see [30, 31] for recent developments and precisions).


## 5. Related Transition Times for the Moran Model with Mutations

Consider the ergodic Moran birth and death Markov chainsXn with mutations on the state space {0,…,N}. By (16), the mean return time to state {0} is E(τ0,0*)=1/π0. Thus, with ρ=(μ¯1μ¯1μμ/μ2μ2)1/μ¯<1, regardless of μ<1 or 2>μ>1(149)E(τ0,0*)~μμ¯1μ2(1ρ)N,asN→∞.Whenμ=1,E(τ0,0*)=μ¯1-N=μ2-N (=2N if μ1=μ2=1/2).Ifμ=2(μ1=μ2=1), E(τ0,0*)=1/π0=(2NN)~4N/πN.We wish now to estimate the expected time it takes forXn to move from one end of the state space to the other, that is, from state {0} to {N}. Let τ0,N be the random time of such a sweep (note that the expected fixation time of a wild-type allele is given by E(τ1,N)=E(τ0,N)-1/p0 as a result of τ0,N=dτ1,N+G, the sum of τ1,N and an independent geometrically distributed rv G≥1 with success probability p0). We will prove the following estimation of its mean value growing geometrically fast.Proposition 20.
Withρ¯=(μ¯2μ¯2μμ/μ1μ1)1/μ¯<1, regardless of μ<1 or 2>μ>1, it holds that
(150)E(τ0,N)~1μ2μμ¯2μ1(1ρ¯)N,asN→∞.
Whenμ=1, E(τ0,N)~μ1-N/μ2 (=2N+1 in the lazy Ehrenfest case: μ1=μ2=1/2). When μ=2, E(τ0,N)~4N/πN.Proof.
Let indeedτx,x+1 be the random time to first hit the state {x+1} starting from the state {x}. Depending on whether the move starting in x is up, down, or no move, τx,x+1 is either 1 or 1+τx,x+1′ (with τx,x+1′ a statistical copy of τx,x+1) or 1+τx-1,x+τx,x+1′. If we let μx be the mean value of τx,x+1, we thus get μx=px+rx(1+μx)+qx(1+μx-1+μx), x≥2. This leads to the recurrence (μ0=1/p0):
(151)μx=qxpxμx-1+1px,x≥1.
This recurrence can be solved explicitly to give(152)μx=1pxπx∑y=0xπy=:πxcpxπx.
Thus,(153)E(τ0,N)=∑x=0N-1‍μx=∑x=0N-1‍1pxπx∑y=0x‍πy=∑x=0N-1‍1px∑y=0x‍∏z=yx-1‍qz+1pz.
Looking at this sum formula, one expects that its leading term is(1/pN-1πN-1)∑y=0N-1‍πy because this is where ∑y=0x‍πy is the largest and πx the smallest. Observing πN-1pN-1=πNqN and qN=μ2,∑y=0N-1‍πy=1-πN→1, this estimation leads, recalling πN~μ1/μμ-2ρ-N, to
(154)E(τ0,N)~1pN-1πN-1~1μ2πN~1μ2μμ¯2μ1(1ρ-)N,
and this is indeed true because the remaining terms do not contribute significantly, see [32]. When μ=1,
(155)E(τ0,N)~1μ2πN=μ1-Nμ2.
Whenμ=2 (Laplace-Bernoulli), E(τ0,N)~1/μ2πN=(2NN).Remark 21.
Regardless of whetherμ≤1 or 2>μ>1, ρ>ρ- if and only if μ2>μ1. In this case, the limiting growth rate of the expected value of τ0,N dominates the one of τ0,0*. On the contrary, μ1>μ2⇒ρ<ρ-.Note that ifμ1=μ2=μ/2 (symmetric Moran model), ρ=ρ- and  (E(τ0,N)/E(τ0,0*))1/N→1  as N→∞, whereas E(τ0,N)/E(τ0,0*)→2/μ (>1 if μ<2, =1 if μ=2).The physical image of the situation is the following: it takes a small𝒪(N) amount of time to move from {0} to the stable equilibrium state {Nu*} because the drift is favorable to this displacement. Once at state {Nu*},Xn returns rapidly and very often to {Nu*} (it takes a 𝒪(N) time to do so), and coming back either to {0} or to {N} is hard because the drift now pushes Xn inside the domain {0,…,N}… As a result, both E(τ0,N) and E(τ0,0*) are, comparatively, geometrically very large, with rates computed just before. Recalling it takes about (NlogN)/μ steps for the chain to reach maximum entropy under the “sep-distance,” starting from {0}, although the return to {0} will eventually occur, the time it takes is so huge that such returns will not be observed. In other words, the chain will reach rapidly the invariant measure inside a single excursion. Let us finally check these last statements.First Return Time to Equilibrium State.Indeed, from (75), using (86),
(156)πNu*=(-Nν1Nu*)(-Nν2N(1-u*))(-NνN)~12πμN.Thus, by (16), E(τNu*,Nu**)~2πμN=𝒪(N), and the return time to the equilibrium state is relatively small.First Hitting Time of the Equilibrium State, Starting from
{
0
}. In this case, we have
(157)E(τ0,Nu*)=∑x=0Nu*-11px∑y=0x∏z=yx-1qz+1pz.If we scale this characteristic time toN, the summations for large N can be replaced by integrals, leading to the estimate
(158)E(τ0,Nu*)~N2∫0u*dv(1-v)(μ1+μ¯v)∫0vdueN(V(u)-V(v)).Recalling(159)∫0vdueNV(u)~eNV(u*)π2N|V′′(u*)|,ifv≤u*,
two successive appeals to the saddle point method show that E(τ0,Nu*)~(πN2/2N|V′′(u*)|)(μ2/μ1μ2)=(π/2μ)N, which is indeed 𝒪(N).Concerning the related first hitting time of stateNu*, starting from Nu*-1 (one step below), we have
(160)E(τNu*-1,Nu*)~Nμ2μ1μ2∫0u*dueN(V(u)-V(u*))~μ2μ1μ2πN2|V′′(u*)|=πμN2μ1μ2,
showing that E(τNu*-1,Nu*) is also 𝒪(N), in accordance with E(τNu*,Nu**)=𝒪(N).



---

*Source: 939308-2013-09-11.xml*