Table of Contents
ISRN Biomathematics
Volume 2013, Article ID 939308, 21 pages
http://dx.doi.org/10.1155/2013/939308
Research Article

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

Laboratoire de Physique Théorique et Modélisation, CNRS-UMR 8089 et Université de Cergy-Pontoise, 2 avenue Adolphe Chauvin, 95302 Cergy-Pontoise, France

Received 27 May 2013; Accepted 2 July 2013

Academic Editors: E. M. Cherry, B. Foy, and M. Glavinovic

Copyright © 2013 Thierry E. Huillet. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 size 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 , the first return time to the equilibrium state, and the first hitting time of starting from , together with the time needed for the walker to reach its invariant measure, again starting from . 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 . 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 (-theorem) and the recurrence of states.

1. Introduction

In Section 2, 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 -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 type alleles in an haploid population with two alleles whose total size 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 takes on the interpretation of the number of agents that belong to the first strategy (say right- or left-wing voters) at time . 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 , the first return time to the equilibrium state, and the first hitting time of starting from , together with the time needed for the Moran walker to reach its invariant measure, again starting from . 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 allele exceed the one favoring allele , the fixation probability of a type mutant tends to as grows large with a correction. On the contrary, this fixation probability tends to geometrically fast as grows large. In the critical case when the two mutation probabilities are egal, this fixation probability goes to from below.Under the additional presence of selection, the fixation condition can be achieved with probability even if the mutation probability favoring allele exceeds the one favoring allele . The condition is that the relative fitness of allele 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 Section 4, following the general stream of ideas developed in [47]. 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 limit (i.e., in the thermodynamic limit). From our analysis, it follows that the mutation Moran walker started at reaches stationarity in about time, far below the time needed to first return to the origin (or excursion’s length) which is geometrically large with . 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

Let be a tridiagonal (Jacobi) irreducible stochastic matrix with , , and , corresponding to the transition probability matrix of some ergodic discrete-time Markov chain on the state space . With a column vector of s, is stochastic if .

The left (row-) eigenvector associated to the eigenvalue satisfies , where   denotes transposition. Its components are given by and , , with chosen so as . It corresponds to the invariant probability measure of the chain.

The matrix is diagonally equivalent to the symmetric matrix so with real eigenvalues. Indeed, with , we have . As a symmetric matrix, is diagonalizable by an orthogonal transformation, and so, is diagonalizable.

Such nearest-neighbors random walks models are well known to be reversible (detailed balance holds). Indeed, , where is the transition matrix of the time-reversed process, given by (with denoting the transpose of ). The transition matrix of reversible Markov chains has real eigenvalues.

Under our assumptions, irreducible chains governed by are positive recurrent, meaning that every state recurs with probability , the mean time between recurrence to state being finite and equal to .

2.2. The Karlin-McGregor (KMG) Spectral Theory

Consider the polynomials in the variable , determined by and the 3-term recurrences: therefore, with ,, . The polynomials satisfy , and are of degree in . They are often called the random walk polynomials. Let , be the zeroes of the polynomial of degree , namely, with . The set constitutes the spectrum of (see [8, page ]). The quantity is the spectral gap.

Let be the spectral probability measure on with respect to which are orthogonal, namely, where , are the potential coefficients.

These polynomials are important in view of the KMG spectral representation theorem [8]. Indeed, with , we have and therefore, . Selecting the row , multiplying by , integrating with respect to , and applying the orthogonality relations (5), we get

Necessarily, because by the ergodic theorem, for all, as . From (7), because , we have , and so the generating function of is the Stieltjes transform of the spectral measure, satisfying as a result of . On the other side, let be the probability generating function (pgf) of the first return time to zero probability. As it can easily be checked by renewal arguments (see [9, page  330]), , showing that

We have and .

More generally, with , let be the Green potential function of the chain. We let , the expected sojourn time at , starting from . Note that and .(i)With the first hitting time of starting from , by renewal arguments, the pgf of reads with . We have With , we have (ii)With the first return time to , we have with

The probability that the first return time to is finite is if and only if (a recurrence condition), and this happens if and only if (stochastic) is irreducible, in which case all states are in the same recurrence class. Else, if for some (a transience condition for ), this probability is smaller than , and one possible reason is that absorption occurs in at least one state ; all states then but the absorbing ones (typically here, either or , or both endpoints) are in the same transience class, and because the walker may or not return to before hitting one of the absorbing states. So 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 .

In view of and in case of recurrence,

Equations (6) and (7) suggest that, with the vectors and can be chosen to be the right (column) and left (row) eigenvectors of associated to the eigenvalue for which

In particular, and , as required. We have because , with the right normalization. Equation (19) yields another orthogonality relation which is dual to the one in (5), namely,

With , we now get the decomposition of into orthogonal projectors which is (7):

Note the biorthogonality relations: from which it follows that

There is a matrix version of all this. Calling (resp., ) the matrix whose column (row) is (resp., ), (18) reads and where   diag, . Since by (19),   diag, we have , which is (22). Since , . If we let be the matrix whose entry is , then and is the matrix formulation of (17). Let us finally check that and (i.e., (17)) are indeed admissible left and right eigen matrices. The relation is also which is which holds from (6). The relation is also or which is or again , because by the reversibility property, . If ( is involutive), then which means . If in addition , then : 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 of associated to its spectrum , in terms of the orthogonal polynomials evaluated at (the -matrix).

The interpretation of the right eigenvectors of is as follows. We have . Fixing the row , the left-hand-side of the last equality is . Therefore,

Suppose that we want to compute for some test function . Let then be the decomposition of in the basis . We have .

Thus, and

Whenever the spectral structure of the chain is known, the solution 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 that is even with . Assume that a random walk on is given so that is the transition matrix of some symmetric random walk reflected at the boundaries . The KMG theory is also relevant for such transition matrices . Let ,  be the zeroes of the polynomial of degree , namely, with . Then, is the probability measure on with respect to which the polynomial associated to is orthogonal. This measure is symmetric on . In particular, and , are eigenvalues of such s. When is odd, say with , the spectral measure is again symmetric, but no longer is an eigenvalue, and .

Example 1. The simplest such example is the simple gambler random walk for which , , () with a pure reflection at the endpoints (). The invariant measure is a truncated geometric distribution. It follows from [9, page ] that ,,,. The spectral gap tends to as if . The mass of the spectral measure can be seen to be equal to , , (with when ). It exhibits boundary effects causing deviation from the uniform measure on . The orthogonal polynomials involve two sine functions. A closed form expression of the generating functions of and follows.
Additional Comment. Let be the tridiagonal transition matrix associated to this RW . Let be the transition matrix associated to . Then, has the same structure as except that and are interchanged (a drift reversal property). is obtained from by applying a permutation reversing the order of the rows and then of the columns. So and share the same spectrum, and this is why is a symmetric function of and . However, the mass of the spectral measure of is different with: .

2.4. Special Cases: Spectral Positivity

Let be a subdiagonal matrix. It is the transpose of a superdiagonal matrix . Then, with diag

Thus, is the sum of a diagonal matrix and a symmetric positive definite matrix. From this, we conclude that if the holding probabilities , for all , then for all , and so and then are positive definite with real positive eigenvalues. The Jacobi matrix then has all its principal minors nonnegative and therefore is oscillatory (i.e., totally nonnegative and such that is totally positive) with all its minors nonnegative (see e.g., [10, page  99]). Oscillatory stochastic matrices have distinct positive eigenvalues, with: . As a sequence of numbers, the right row eigenvector associated to has exactly sign changes (see [10, page  101]). And so do the numbers: . The above condition is a simple sufficient condition allowing to exhibit examples of that are oscillatory but of course, it is not necessary; the full condition under our assumptions being for all or

Equivalently, if are the positive eigenvalues of , we need to check for all .

For a simple random walk whose -transition matrix is spectrally nonnegative (resp., spectrally positive), there exists a symmetric random walk on (resp., on ), reflected at the endpoints, started at an even integer (resp., odd integer ), such that , (resp., ). 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 is symmetric on , and by passing to , the spectrum is being folded: If (resp., ) is the symmetric spectral measure of with (resp., ), then is the spectral measure of . Let and be the up and down probabilities that in one step given is in state different from the endpoints, , then

This, together with and (resp., ), allows to determine recursively the transition matrix of from the one of . From these facts, we get the following proposition.

Proposition 2. If a BD chain is spectrally nonnegative, then it is stochastically monotone.

Proof. Under our hypothesis indeed, because this is also . This condition means that 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 chain with transition matrix is not spectrally nonnegative, so has strictly negative eigenvalues (). With , consider the lazy modified stochastic matrix

The eigenvalues of are , with the same left and right eigenvectors as the ones of . One can always choose small enough so that . For this , is now spectrally nonnegative, and so governs a stochastically monotone chain . The RW , with shifted spectrum, clearly has the same invariant measure as the one of .

For example, solving (34) with the boundary conditions and gives the forward and backward recurrences of the unknown odd and even probability terms in the bulk:

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 height for some excursion (sample paths of between two consecutive visits to state ). This event will be realized if and only if (i) downward paths started from hit state before hitting state and (ii) upward paths started at first reach (with probability ) and then, paths started at hit without returning to again in the intervening time. These two events are independent. Therefore, with , the first hitting time of starting from

Assume . Let denote the random walk stopped when it first hits . Its transition matrix is except for its first row which is , translating that state is now absorbing. Define the scale (or harmonic) function of this random walk as the function which makes a martingale. The function is important because, as is well-known, for all , with , the first hitting time of starting from

Using this remark, we get

Note that if . The event may also occur, with probability .

The interpretation the event takes on is twofold: either it is the maximal height of an excursion if the state remains reflecting, either, if state is itself assumed absorbing, constitutes the probability of a fixation event at , starting from state . We can thus compute the probability of a fixation before extinction at occurs. It remains to compute (defined up to an arbitrary multiplicative constant). We wish to have , leading to and where , . Imposing (fixing the arbitrary multiplicative constant), the searched “harmonic” function is where satisfies . Thus, and

Equations (39) and (41) characterize the law of the excursion height of the random walk . We note that increases with . From (38) and (41), we get an explicit expression of the probability of the event : where is the first hitting time of starting from , with . In particular, if with also absorbing, is the fixation probability starting from , whereas is the extinction probability.

3. A BD Example: The Moran Model

The example we have in mind is the -allele Moran model with bias mechanism . Let be continuous and .

The Moran model is characterized by the following transition probabilities, given the walker is in state :

For a Moran model is the amount of say allele or type individuals in a population of size , and the reproduction law is obtained while choosing two individuals at random among , 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 over the generations. Given then , the transition () occurs with probability () translating the fact that the individual bound to die is chosen among type (type ) individuals with probability    , and the one which splits into two new born is chosen among the other type, with probability (resp., ), after deforming the neutral frequency by the bias mechanism . These bias take into account various evolutionary forces which are the causes of deviation to neutrality (see [1214] 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 (no external deformation). Extended Moran models are studied in [15].

Assuming and , with , the corresponding Markov chain is ergodic with invariant distribution: where

We have the following obvious facts.

Proposition 4. Let be a Moran process with transition probabilities (44) for some bias . Then, is again a Moran process with transition probabilities of the type (44) but with the new bias , namely, , . The spectrum of the transition matrix of is the same as the one of . If the bias is such that , then the transition matrix of coincides with the one of .

Proof. We have , and is obtained from applying the permutation reversing the order of the rows and then of the columns.

3.1. Fixation Probability

Forcing the states to be absorbing, the scale function reads with being the fixation probability, whereas is the extinction probability, starting from (for the neutral model with , the scale function is , and the fixation probability is trivially ). The probability of a fixation before extinction occurs in an excursion reads

We have where is the fixation probability of a mutant-type allele , for example, the probability that, starting from , the RW hits state (fixation of the mutant-type allele ) before hitting state (extinction of allele ).

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 where are mutation probabilities in . We  let  ,   , , and  .

Schematically, type alleles are being produced whenever type alleles do not mutate to type (with probability ) or type alleles mutate to type (with probability )

The transition matrix of this model has therefore the particular BD structure with , . When , we introduce the auxiliary parameters with   and.

The drift of this RW is

The deterministic part of the dynamics of has a stable point at , cancelling the drift. When is nondecreasing, we have . We note that in the expression of , the roles of and are interchanged.   So, if and only if (symmetric mutations).

If , we recover if is even the heat-exchange Bernoulli-Laplace model [9] as a borderline example, but is strictly decreasing in this case. If , then and we obtain an aperiodic model amenable (through a suitable time substitution) to the Ehrenfest urn model provided is even.

For one-way mutations, or leading to the simpler forms and of, but with and , respectively, corresponding to (resp., ) being absorbing. We avoid this situation because it deserves another treatment due to transience.

Moran Model with Mutations and Selection. Another example could be which is the composition of the selection bias mechanism with the mutation mechanism (50). When considering the fixation probability 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 leading to in which case the RW is spectrally positive. Assume is continuous and nondecreasing, with .

Then, as can easily be checked: for all if and only if with . Indeed, imposing for all leads to if and if . So, if is nondecreasing with , then .(i)Take, for example, the mutation mechanism now satisfying . The condition leads to (fair mutation). In this case, with . If , then and , , and .  For example, if , has three distinct real nonnegative eigenvalues ordered in decreasing size. We note however that it is not necessary that for to be spectrally positive; indeed, for the full Moran model with mutations , under the sole condition ( is nondecreasing), the eigenvalues can be seen to be nonnegative. Still, we were not able to find a sharp general condition on that would lead to a spectrally nonnegative Moran chain.(ii)For the mechanism of (55) combining mutations with selection, assuming to guarantee that is nondecreasing, the condition relates to give , where both . 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 for large and then will follow.(i) Suppose first that and let .

Proposition 5. The fixation probability satisfies the sharp threshold property: (i)if (), then ; (ii)if (), then ; (iii)if (), then .

Proof. With , we first have the identity for sums of inverse binomial coefficients (see e.g., [16, page  197])
There is therefore a closed-form expression of the fixation probability . Observing now
we get(i)  (): . So: ;(ii) (): . So: ;(iii) (): . So: .
If , because if : and is very close to . The fixation probability tends to from below.
If , + + , and is very close to if and very close to up to a factor if .
When , the probability, starting from to hit the state (fixation) starting from before returning to , tends to geometrically fast. On the contrary, would , this probability would converge to ; mutations producing type individuals are sufficiently large to enhance fixation of a mutant allele with probability . If , the fixation probability approaches .

(ii)Suppose now that . With , , Proceeding as before, we have showing that . The dominating term in is thus either the one corresponding to or the one corresponding to .

With , by Stirling formula, we get

Now, whenever , regardless of whether or ,

Indeed, let . The condition is also

If , then if and only if (or ) and then .

If , then if and only if (or ) and then .

We obtained the following proposition.

Proposition 6. When , the fixation probability also satisfies the following sharp threshold property.(i)If , the fixation probability tends to 0 geometrically fast because , leading to . The fixation probability of a mutant tends to 0, algebraically fast at rate given by (63). (ii)If , regardless of whether or , . Thus, leading to from below.(iii)If (balanced mutations), , and therefore, : the fixation probability tends to from below.

(i) For the Moran model with selection mechanism as in (56), gene (resp., ) has fitness (resp., ) so that, with , the mutant relative fitness showing that:. Whenever (), the mutant gene is selectively advantageous (deleterious) compared to . When (), tends to , whereas when (), this probability tends to exponentially fast like . When (neutral case), tends to 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 the Moran model with bias mechanism (55) involving both mutation and selection, where, with ,