On Discrete-Time Multiallelic Evolutionary Dynamics Driven by Selection

Abstract

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

1. Introduction

Population genetics aims at elucidating the fate of genotype frequencies undergoing the basic evolutionary processes when various driving “forces” such as fitness, mutation, or recombination are at stake in the gene pool. This requires to clarify the updating mechanisms of the gene frequency-distributions over time. Another important additional driving source is the genetic drift whose nature is exclusively random. The corresponding field of interest is the statistical theory arising from this aspect of the gene replacement processes and it requires some use of the Markov chain theory (see [1]).

In this note, we revisit the basics of both the deterministic and stochastic dynamics arising in discrete-time asexual multiallele evolutionary genetics driven only by fitness. We do not touch at all upon other important mechanisms such as mutation. We start with the haploid case with 𝐾 alleles before switching to the more interesting diploid case.

Let us summarize and comment the material developed in Section 2. In the deterministic haploid case, a vector of fitness is attached to the alleles. The interest is on the evolution of the allelic frequency distribution over time. The updates of the allele frequency distributions are driven by the relative fitnesses of the alleles, ending up in a state where only the fittest monomorphic state will survive. This state is also an extremal point of the simplex over which the dynamics takes place. From this dynamics, it appears that the mean fitness increases as time passes by, the rate of increase being the variance in relative fitness (a well-known particular incarnation of the Fisher theorem of natural selection). We recall the background and the evolution equations.

The haploid replicator dynamics, as a nonlinear updating mapping from the simplex to the simplex, may be viewed as the discrete-time nonlinear Master equation of some random Markov process. We supply a construction for this process thereby giving a stochastic interpretation to the original deterministic formulation of the dynamics.

In the diploid case, there is a similar deterministic updating dynamics but now on the full array of the genotype frequencies. It involves the fitness matrix attached to the genotypes. When mating is random so that the Hardy-Weinberg law applies, we may look at the induced marginal allelic frequencies dynamics. The updating dynamics looks quite similar to the one occurring in the haploid case except that the mean fitness is now the mean fitness quadratic form in the current frequencies whereas marginal fitnesses are no longer constant but affine functions in these frequencies. As for the haploid case, it is possible to construct a Markov process whose nonlinear Master equation coincides with the diploid replicator dynamics. We supply this construction which we believe is new.

In the diploid context, the Fisher theorem still holds true but, as a result of the fitness landscape being more complex, there is a possibility for a polymorphic equilibrium state to emerge. Due to its major evolutionary interest, our subsequent concern is to identify examples of diploid dynamics leading to a unique polymorphic state on the simplex, either unstable or stable. We start with the unstable case and draw the attention on a class of fitness matrices leading to a unique unstable polymorphic equilibrium state: the class of strictly potential matrices. Strictly ultrametric matrices are particular instances of strictly potential matrices [2] which therefore will display unstable polymorphism as well. There is a useful canonical representation of strictly ultrametric matrices due to [3] which we recall which helps giving specific examples of strictly ultrametric matrices. When dealing with the class of strictly potential fitness matrices, the mean fitness quadratic form is definite-positive; we derive a related class of fitness matrices leading to a definite-negative mean fitness quadratic form. For this class of matrices, there will also be a unique polymorphic equilibrium state for the diploid dynamics and it will be stable. We also draw the attention on a subclass of the latter: the so-called “antistrictly ultrametric matrices.” For such matrices, among other things, the fitness of each homozygote should not exceed the ones of all the heterozygotes carrying the allele of the homozygote. To the best of our knowledge, the large class of potential fitness matrices as natural candidates for polymorphic states to emerge was not discussed in the literature.

Section 3 is devoted to the stochastic version of these considerations when the transitions in the constitutive allelic population sizes are given by a 𝐾-dimensional Wright-Fisher model with total constant-size 𝑁 (see [1, 4]). It takes into account an additional important driving source of evolution, namely, the genetic drift, whose nature is exclusively random. Under fitness only and in particular in the absence of mutations, the multiallelic Wright-Fisher model is a transient Markov chain on a (^{𝑁+𝐾−1}_𝐾−1)-dimensional state-space whose absorbing states are the monomorphic states. We give an expression for the fixation probabilities of this process. Then, we develop four conditioning problems: conditioning on fixating in a given monomorphic state, conditioning on avoiding the extremal states before the current instant, conditioning on nonfixation at each transition time, and conditioning on avoiding the extremal states in the remote future. For the last three conditioned processes, the equilibrium structure of the fitness mechanism shows up.

Finally, we run into similar considerations for the Moran model. When dealing with the fixation probabilities in this Moran context, we suggest a new mean-field approximation of these probabilities which is based on a well-known explicit formula for the 2-alleles case [1]. It concerns the case of multiplicative fitnesses only. Finally, we consider the Moran model conditioned on nonfixation at each transition time. We exploit the reversible character of this process to derive a new explicit product formula for its invariant probability measure.

2. Deterministic Evolutionary Dynamics

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

2.1. Single Locus: Haploid Population with 𝐾 Alleles

Consider 𝐾 alleles 𝐴_𝑘,𝑘=1,…,𝐾 attached to a single locus. Suppose that the current time-𝑡 allelic frequency distribution is given by the column vector 𝐱∶=𝑥_𝑘,𝑘=1,…,𝐾. (Throughout, a boldface variable, say 𝐱, will represent a column-vector so that its transpose, say 𝐱^∗, will be a line-vector.) We therefore have 𝐱∈𝑆_𝐾={𝐱≥𝟎∶|𝐱|∶=∑^𝐾_𝑘=1𝑥_𝑘=1}, the 𝐾-simplex as a convex subset of ℝ^𝐾 with dimension 𝐾−1. Let 𝐰∶=𝑤_𝑘>0,𝑘=1,…,𝐾 denote the absolute fitnesses of the alleles. Let 𝑤(𝐱)∶=_𝑙𝑤_𝑙𝑥_𝑙=𝐰^∗𝐱 (2.1) be the mean fitness of the population at time 𝑡. We will also need 𝜎²(𝐱)=^𝐾_𝑘=1𝑥_𝑘(𝑤_𝑘−𝑤(𝐱))², (2.2) the variance in absolute fitness, and𝜎²(𝐱)=^𝐾_𝑘=1𝑥_𝑘(𝑤_𝑘𝑤(𝐱)−1)²=𝜎²(𝐱)𝑤(𝐱)², (2.3) the variance in relative fitness 𝑤_𝑘/𝑤(𝐱).

2.1.1. Dynamics

From the deterministic evolutionary genetics point of view, the discrete-time update of the allele frequency distribution on the simplex 𝑆_𝐾 is given by 𝑥^′_𝑘=𝑝_𝑘𝑥(𝐱)∶=^𝑘𝑤_𝑘𝑤(𝐱),𝑘=1,…,𝐾. (2.4) (the symbol ′ is a common and useful notation to denote the updated frequency.) The quantity 𝑤_𝑘/𝑤(𝐱)−1 therefore interprets as the frequency-dependent Malthus growth rate parameter of 𝑥_𝑘. As required, the vector 𝐩(𝐱)∶=𝑝_𝑘(𝐱),𝑘=1,…,𝐾, maps 𝑆_𝐾 into 𝑆_𝐾. In vector form, with 𝐷_𝐱∶=diag(𝑥_𝑘,𝑘=1,…,𝐾), the nonlinear deterministic evolutionary dynamics reads 𝐱^′1=𝐩(𝐱)=𝐷𝑤(𝐱)^𝐰1𝐱=𝐷𝑤(𝐱)^𝐱𝐰, (2.5) or, with Δ𝐱∶=𝐱^′−𝐱, the increment of 𝐱1Δ𝐱=(𝐷𝑤(𝐱)^𝐰−𝐼)𝐱. (2.6) Avoiding the trivial case where fitnesses are all equal, without loss of generality, we can assume that either 𝑤₁≥⋯≥𝑤_𝐾=1 or 𝑤₁≤⋯≤𝑤_𝐾=1. Thus that allele 𝐴₁ or 𝐴_𝐾 has largest fitness.

2.1.2. Mean Fitness Increase

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

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

A vector 𝐱 of 𝑆_𝐾 can be thought of as a probability vector. The dynamical equation (2.4), as a nonlinear update mapping from 𝑆_𝐾 to 𝑆_𝐾, may be viewed as the discrete-time nonlinear Master equation of some Markov process whose construction we now give. Suppose that we have a population of 𝑁 haploid individuals each of which can be of one among 𝐾 types or colors (carrying one among the 𝐾 possible alleles). We will need to introduce an extra color-state, say 𝜕={0}, which will be absorbing for the process we will now construct. Let 𝐊(𝑡)∶=𝐾_𝑛(𝑡),𝑛=1,…,𝑁 be the random color distribution of the population at time 𝑡, therefore with enlarged state-space {0,1,…,𝐾}. Assume that the individuals are indistinguishable leading to the exchangeability property: 𝐾_𝑛(𝑡)=^𝑑𝐾₁(𝑡),𝑛=2,…,𝑁 (equality in distribution). Let 𝑈_𝑡,𝑚,𝑡=1,2,…;𝑚=1,…,𝑁 be 𝑁 i.i.d. driving sequences of uniformly distributed random variables on [0,1], independent of 𝐊(𝑡). Let 𝑤_𝑘∶=𝑤_𝑘/∑_𝑘𝑤_𝑘,𝑘=1,…,𝐾. To decide the allele 𝐾_𝑚(𝑡+1) carried by the individual number 𝑚∈{1,…,𝑁} at time 𝑡+1, with 𝟏(𝐴) being the indicator function of the event 𝐴, consider the random Markovian dynamical system:𝟏(𝐾_𝑚(𝑡+1)=𝑘,𝜏_𝑚>𝑡+1)=𝟏(𝑤_𝑘𝑁^𝑁_𝑛=1𝟏(𝐾_𝑛(𝑡)=𝑘,𝜏_𝑛>𝑡)>𝑈_𝑡+1,𝑚). (2.9) Here 𝑘∈{1,…,𝐾} and 𝜏_𝑛 is the first time that 𝐾_𝑛(𝑡) hits the absorbing state 𝜕. As a result of 𝐾_𝑛(𝑡)=^𝑑𝐾₁(𝑡), we naturally assume 𝜏=_𝑛^𝑑𝜏₁,𝑛=2,…,𝑁. For each 𝑛, our model therefore attributes a positive probability that 𝐾_𝑛(𝑡)=0 for all 𝑡≥𝜏_𝑛. Although, in principle, there is a possibility that the type of the 𝑚th particle is the one of the fictitious unobservable allele 𝐴₀, as a result of (2.9), the sample paths of 𝐊(𝑡) leading to this 𝐴₀ are ruled out because focus is on the observable states only.

In words, for the dynamics (2.9), the observable event 𝐾_𝑚(𝑡+1)=𝑘 is realized (together then with 𝜏_𝑚>𝑡+1) if the proportion at 𝑡 of type-𝑘 individuals, weighted by the corresponding scaled fitness 𝑤_𝑘, is large enough (compared to 𝑈_𝑡+1) and of course if the whole process was not absorbed at {0} in the previous step. Taking first the expectation with respect to the driving noise 𝑈_𝑡+1,𝑚 in (2.9), we get ℙ(𝐾_𝑚(𝑡+1)=𝑘,𝜏_𝑚>𝑡+1∣𝐊(𝑡)>𝟎)=𝑤_𝑘𝑁^𝑁_𝑛=1𝟏(𝐾_𝑛(𝑡)=𝑘,𝜏_𝑛>𝑡). (2.10) Putting 𝑧_𝑘(𝑡)∶=ℙ(𝐾₁(𝑡)=𝑘,𝜏₁>𝑡), recalling (𝐾_𝑛(𝑡);𝜏_𝑛)=^𝑑(𝐾₁(𝑡);𝜏₁),𝑛=2,…,𝑁, and taking the expectation with respect to 𝐊(𝑡), we get an unnormalized version of (2.4): 𝑧_𝑘(𝑡+1)=𝑤_𝑘𝑧_𝑘(𝑡),𝑘∈{1,…,𝐾}. (2.11) We have 1>ℙ(𝜏₁>𝑡)=∑^𝐾_𝑘=1𝑧_𝑘(𝑡)=∑^𝐾_𝑘=1𝑥_𝑘(0)𝑤^𝑡_𝑘→0, geometrically fast. Defining the normalized conditional probabilities 𝑥_𝑘𝑧(𝑡)=^𝑘(𝑡)∑^𝐾_𝑘=1𝑧_𝑘(𝑡)=ℙ(𝐾₁(𝑡)=𝑘∣𝜏₁>𝑡), (2.12) we obtain the normalized haploid dynamics (2.4): 𝑥^′_𝑘=𝑤_𝑘𝑥_𝑘∑^𝐾_𝑘=1𝑤_𝑘𝑥_𝑘,𝑘∈{1,…,𝐾}. (2.13) It now may be viewed as the nonlinear Master equation of some stochastic Markov process. Let us make some miscellaneous remarks.

(i) Clearly this construction makes also sense if 𝑁=1 (a single particle). (ii) When 𝑁 is finite, we should stress that the initial condition can be chosen to be deterministic, say with 𝑥_𝑘(0)=𝑧_𝑘(0)=𝑖_𝑘/𝑁, for some sequence of integers 𝑖_𝑘≥0 satisfying 𝑖₁+⋯+𝑖_𝐾=𝑁 (𝑖₀=0) and quantifying the initial population sizes. It could also be chosen to be random, with 𝑥_𝑘(0) defining the initial probability distribution of the alleles. This occurs in the large 𝑁 limit if 𝑖_𝑘=⌊𝑁𝑥_𝑘(0)⌋ so that 𝑖_𝑘/𝑁→𝑥_𝑘(0). The latter choice may therefore be interpreted as a large 𝑁 limit of the former one. (iii) In the stochastic interpretation (2.9) of the deterministic dynamics (2.4), 𝑥_𝑘(𝑡) can be interpreted either as the probability that the random allele carried by a typical individual is 𝐴_𝑘 or like the expected proportion of the individuals of type 𝑘 within the whole population (a frequentist point of view). (iv) The appeal to the coffin state 𝜕 was a necessary step to understand the normalization 𝑧_𝑘→𝑥_𝑘.(v) Even though 𝐊(𝑡) is exchangeable, it is not true that, with 𝑛₁≠𝑛₂ any two distinct individuals, their random labels 𝐾_𝑛1(𝑡) and 𝐾_𝑛2(𝑡) are independent. The random algorithm allowing to update the joint types of 𝐾_𝑛1(𝑡) and 𝐾_𝑛2(𝑡) could be written down but is much more involved.

2.2. Single Locus: Diploid Population with 𝐾 Alleles

We now run into similar considerations but with diploid populations.

2.2.1. Joint Evolutionary Dynamics

Let 𝑤_𝑘,𝑙≥0,𝑘,𝑙=1,…,𝐾 stand for the absolute fitness of the genotypes 𝐴_𝑘𝐴_𝑙 attached to a single locus. Assume 𝑤_𝑘,𝑙=𝑤_𝑙,𝑘 (𝑤_𝑘,𝑙 being proportional to the probability of an 𝐴_𝑘𝐴_𝑙 surviving to maturity, it is natural to take 𝑤_𝑘,𝑙=𝑤_𝑙,𝑘). Let then 𝑊 be the symmetric fitness matrix with 𝑘,𝑙-entry 𝑤_𝑘,𝑙.

Assume that the current frequency distribution at time 𝑡 of the genotypes 𝐴_𝑘𝐴_𝑙 is given by 𝑥_𝑘,𝑙. Let 𝑋 be the frequencies array with 𝑘,𝑙-entry 𝑥_𝑘,𝑙. The joint evolutionary dynamics in the diploid case is given by the updating:𝑥^′_𝑘,𝑙=𝑥_𝑘,𝑙𝑤_𝑘,𝑙,𝜔(𝑋) (2.14) where the mean fitness 𝜔 is now given by 𝜔(𝑋)=∑_𝑘,𝑙𝑥_𝑘,𝑙𝑤_𝑘,𝑙. The relative fitness of the genotype 𝐴_𝑘𝐴_𝑙 is 𝑤_𝑘,𝑙/𝜔(𝑋). The joint dynamics takes the matrix form: 𝑋^′=11𝜔(𝑋)𝑋∘𝑊=𝜔(𝑋)𝑊∘𝑋, (2.15) where ∘ stands for the (commutative) Hadamard product of matrices.

Let 𝐽 be the flat 𝐾×𝐾 matrix whose entries are all 1. Then Δ𝑋∶=𝑋^′1−𝑋=1𝜔(𝑋)(𝑋−𝐽)∘𝑊=𝜔(𝑋)𝑊∘(𝑋−𝐽). (2.16) We will also let𝜎²(𝑋)=^𝐾_𝑘,𝑙=1𝑥_𝑘,𝑙(𝑤_𝑘,𝑙−𝜔(𝑋))² (2.17) stand for the genotypic variance in absolute fitness and𝜎²(𝑋)=^𝐾_𝑘,𝑙=1𝑥_𝑘,𝑙(𝜔_𝑘,𝑙𝜔(𝑋)−1)²=𝜎²(𝑋)𝜔(𝑋)² (2.18) will stand for the diploid variance in relative fitness.

Consider the problem of evaluating the increase of the mean fitness. We haveΔ𝜔(𝑋)=_𝑘,𝑙Δ𝑥_𝑘,𝑙𝑤_𝑘,𝑙=_𝑘,𝑙𝑥_𝑘,𝑙(𝑤²_𝑘,𝑙𝜔(𝑋)−𝑤_𝑘,𝑙)=𝜔(𝑋)𝜎²(𝑋)>0 (2.19) with a relative rate of increase: Δ𝑤(𝑋)/𝑤(𝑋)=𝜎²(𝑋). This is the full diploid version of the Fisher theorem.

2.2.2. Marginal Allelic Dynamics

Assuming a Hardy-Weinberg equilibrium, the frequency distribution at time 𝑡, say 𝑥_𝑘,𝑙, of the genotypes 𝐴_𝑘𝐴_𝑙 is given by 𝑥_𝑘,𝑙=𝑥_𝑘𝑥_𝑙 where 𝑥_𝑘=∑_𝑙𝑥_𝑘,𝑙 is the marginal frequency of allele 𝐴_𝑘 in the whole genotypic population . The whole frequency information is now enclosed within 𝐱∶=𝑥_𝑘,𝑘=1,…,𝐾. For instance, the mean fitness is now given by the quadratic form: 𝜔(𝐱)∶=∑_𝑘,𝑙𝑥_𝑘𝑥_𝑙𝑤_𝑘,𝑙=𝐱^∗𝑊𝐱, with 𝐱^∗ being the transposed line vector of the column vector 𝐱=𝑋𝟏 (𝟏 the unit 𝐾-vector). We will also let𝜎²(𝐱)=^𝐾_𝑘,𝑙=1𝑥_𝑘𝑥_𝑙(𝑤_𝑘,𝑙−𝜔(𝐱))² (2.20) stand for the genotypic variance in absolute fitness and𝜎²(𝐱)=^𝐾_𝑘,𝑙=1𝑥_𝑘𝑥_𝑙(𝑤_𝑘,𝑙𝜔(𝐱)−1)²=𝜎²(𝐱)𝜔(𝐱)² (2.21) will stand for the diploid variance in relative fitness.

Consider now the update of the allelic marginal frequencies 𝐱 themselves. If we first define the frequency-dependent marginal fitness of 𝐴_𝑘 by 𝑤_𝑘(𝐱)=(𝑊𝐱)_𝑘∶=∑_𝑙𝑤_𝑘,𝑙𝑥_𝑙, the marginal dynamics is given as in (2.4) by:𝑥^′_𝑘=𝑥_𝑘𝑤_𝑘(𝐱)=1𝜔(𝐱)𝑥𝜔(𝐱)^𝑘(𝑊𝐱)_𝑘=∶𝑝_𝑘(𝐱),𝑘=1,…,𝐾. (2.22) This dynamics involves a multiplicative interaction between 𝑥_𝑘 and (𝑊𝐱)_𝑘, the 𝑘th entry of the image 𝑊𝐱 of 𝐱 by 𝑊. In (2.22) there is a normalization by the quadratic form 𝜔(𝐱)=𝐱^∗𝑊𝐱. In vector form (2.22) reads 𝐱^′=1𝐷𝜔(𝐱)^𝐱1𝑊𝐱=𝐷𝜔(𝐱)^𝑊𝐱𝐱=∶𝐩(𝐱), (2.23) where 𝐩 maps 𝑆_𝐾 into 𝑆_𝐾. Iterating, the time-𝑡 frequency distribution is𝐱(𝑡)=𝐩(𝐩(⋯𝑡times⋯(𝐩(𝐱₀)))). (2.24)

Except for the fact that the mean fitness 𝜔 in (2.22) is now a quadratic form in 𝐱 and that the marginal fitness of 𝐴_𝑘 is now frequency-dependent, depending linearly on 𝐱, as far as the marginal frequencies are concerned, the updating formalism (2.22) in the diploid case looks very similar to the one in (2.4) describing the haploid case.

In the diploid case, assuming fitnesses to be multiplicative, say with 𝑤_𝑘,𝑙=𝑤_𝑘𝑤_𝑙, then 𝑤_𝑘(𝐱)/𝐱^∗𝑊𝐱=𝑤_𝑘/∑_𝑙𝑤_𝑙𝑥_𝑙 and the dynamics (2.22) boils down to (2.4). However, the mean fitness in this case is 𝜔(𝐱)=(∑_𝑙𝑤_𝑙𝑥_𝑙)² and not 𝑤(𝐱)=∑_𝑙𝑤_𝑙𝑥_𝑙 as in the haploid case.

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

As for the haploid case, there is a Markov chain governed by the Master equation (2.22). Consider a population of diploid individuals. The number of alleles 𝑁 in this population is therefore twice the number of genes. Each allele can be of one among 𝐾 types or colors (carrying one among the 𝐾 possible alleles). As before, we introduce an extra color-state, say 𝜕={0} which is absorbing for the process to be constructed. For 𝑐=1,2, let 𝐊^𝑐(𝑡)∶=𝐾^𝑐_𝑛(𝑡),𝑛=1,…,𝑁 be two independent copies of the random color distribution of the allelic population at time 𝑡. Let 𝐊(𝑡)=(𝐊¹(𝑡),𝐊²(𝑡)). Assume that the alleles are indistinguishable within each sample, leading to 𝐾^𝑐_𝑛(𝑡)=^𝑑𝐾^𝑐₁(𝑡),𝑛=2,…,𝑁,𝑐=1,2. For 𝑐=1,2, let 𝑈^𝑐_𝑡,𝑚,𝑡=1,2,…;𝑚=1,…,𝑁 be two mutually independent i.i.d. driving 𝑁-sequences of uniformly distributed random variables on [0,1] and independent of 𝐊(𝑡). To decide the type of the random allele 𝐾_𝑚(𝑡+1),𝑚=1,…,𝑁, at time 𝑡+1, consider now the Markovian dynamical system:𝟏(𝐾_𝑚(𝑡+1)=𝑘,𝜏_𝑚1>𝑡+1)=𝟏(𝑁^𝑁_𝑛=1𝟏(𝐾¹_𝑛(𝑡)=𝑘,𝜏¹_𝑛>𝑡)>𝑈¹_𝑡+1,𝑚1)𝟏(𝑁^𝑁_𝑛=1𝑊_𝑘,𝐾²_{𝑛(𝑡)𝟏(𝜏2𝑛>𝑡)}>𝑈²_𝑡+1,𝑚). (2.25) Here 𝑘∈{1,…,𝐾},𝜏^𝑐_𝑛 are the first hitting times of each 𝐾^𝑐_𝑛(𝑡) of the absorbing state 𝜕,𝑐=1,2 and 𝑊∶=𝑊/‖𝑊‖ for any matrix norm ‖𝑊‖, say, for example, ‖𝑊‖=∑_𝑘,𝑙𝑤_𝑘,𝑙. We assume 𝜏=^𝑐_𝑛^𝑑𝜏^𝑐₁,𝑛=2,…,𝑁,𝑐=1,2.

In words, for this new dynamics, the observable event 𝐾_𝑚(𝑡+1)=𝑘 is seen to be realized together with 𝜏_𝑚>𝑡+1 if two natural independent conditions are now satisfied which can be read from the two indicator functions in the right-hand side of (2.25):(i)first, the proportion at 𝑡 of type-𝑘 individuals of the first copy should be large enough (compared to 𝑈¹_𝑡+1);(ii)second, for the second sample copy 𝐊²(𝑡), the expected fitness of the genotypes 𝐴_𝑘𝐴_𝑙,𝑙=1,…,𝐾, containing allele 𝐴_𝑘 at 𝑡 should be large enough (compared to 𝑈²_𝑡+1).

As for the haploid case, it is necessary that both processes 𝐊(𝑡) should not be absorbed at {0} in the previous step. Taking first the expectation with respect to the independent driving noises 𝑈^𝑐_𝑡+1,𝑚 in (2.25), we get ℙ(𝐾_𝑚(𝑡+1)=𝑘,𝜏_𝑚=1>𝑡+1∣𝐊(𝑡)>𝟎)𝑁^𝑁_𝑛=1𝟏(𝐾¹_𝑛(𝑡)=𝑘,𝜏¹_𝑛1>𝑡)⋅𝑁^𝑁_𝑛=1𝑊_𝑘,𝐾²_{𝑛(𝑡)𝟏(𝜏2𝑛>𝑡)}. (2.26) Putting 𝑧_𝑘(𝑡)=ℙ(𝐾₁(𝑡)=𝑘,𝜏₁>𝑡), taking the expectation with respect to 𝐊(𝑡) and using our independence and exchangeability hypotheses, we get 𝑧_𝑘(𝑡+1)=𝑧_𝑘(𝑡)⋅𝔼[𝑊_𝑘,𝐾²_{1(𝑡)𝟏(𝜏}²_1>𝑡)]=𝑧_𝑘(𝑡)⋅^𝑘_𝑙=1𝑊_𝑘,𝑙𝑧_𝑙(𝑡), (2.27) corresponding to an unnormalized version of (2.22): 𝑧_𝑘(𝑡+1)=𝑧_𝑘(𝑡)(𝑊𝐳(𝑡))_𝑘,𝑘∈{1,…,𝐾}. (2.28) Defining the normalized conditional probabilities 𝑥_𝑘𝑧(𝑡)=^𝑘(𝑡)∑^𝐾_𝑘=1𝑧_𝑘(𝑡)=ℙ(𝐾₁(𝑡)=𝑘∣𝜏₁>𝑡), (2.29) we obtain the normalized nonlinear Markovian diploid dynamics (2.22): 𝑥^′_𝑘=𝑥_𝑘(𝑊𝐱)_𝑘∑^𝐾_𝑘=1𝑥_𝑘(𝑊𝐱)_𝑘,𝑘∈{1,…,𝐾}. (2.30) Note that this construction makes sense if 𝑁=2 (a single individual and its 2 alleles of 𝐾 possible types). The need for two copies of 𝐊(𝑡) was governed by the quadratic character of the interaction appearing in the numerator of (2.22).

2.2.4. Increase of Mean Fitness

Again, the mean fitness 𝜔(𝐱), as a Lyapunov function, increases as time passes by. We indeed have Δ𝜔(𝐱)=𝜔(𝐱^′1)−𝜔(𝐱)=𝜔(𝐱)²_𝑘,𝑙𝑥_𝑘𝑤_𝑘(𝐱)𝑤_𝑘,𝑙𝑥_𝑙𝑤_𝑙(𝐱)−_𝑘,𝑙𝑥_𝑘𝑤_𝑘,𝑙𝑥_𝑙>0, (2.31) because, defining 0<𝑋(𝐱)∶=∑_𝑘,𝑙𝑥_𝑘(1−𝑤_𝑘(𝐱)/𝜔(𝐱))𝑤_𝑘,𝑙(1−𝑤_𝑙(𝐱)/𝜔(𝐱))𝑥_𝑙, we have 2Δ𝜔(𝐱)=𝑋(𝐱)+𝜔(𝐱)(_𝑘𝑥_𝑘𝑤_𝑘(𝐱)²−𝜔(𝐱)²)>0. (2.32) Its partial rate of increase due to frequency shifts only is 𝛿𝜔(𝐱)∶=∑_𝑘Δ𝑥_𝑘𝑤_𝑘(𝐱). It satisfies𝛿𝜔(𝐱)𝜔(𝐱)=_𝑘𝑥_𝑘(𝑤_𝑘(𝐱)𝜔(𝐱)−1)²=_𝑘(Δ𝑥_𝑘)²𝑥_𝑘=𝜎²_𝐴(𝐱)2, (2.33) where 𝜎²_𝐴(𝐱) is the allelic variance in relative fitness:𝜎²_𝐴(𝐱)∶=2^𝐾_𝑘=1𝑥_𝑘(𝑤_𝑘(𝐱)𝜔(𝐱)−1)². (2.34)

2.2.5. An Alternative Representation of the Allelic Dynamics

There is an alternative vectorial representation of the dynamics (2.22). Define the symmetric positive-definite matrix 𝐺(𝐱)=𝐷_𝐱(𝐼−𝟏𝐱^∗) with quadratic entries in the frequencies: 𝐺(𝐱)_𝑘,𝑙=𝑥_𝑘(𝛿_𝑘,𝑙−𝑥_𝑙). (2.35) Introduce the quantity 𝑉_𝑊(𝐱)=(1/2)log𝜔(𝐱), which is half the logarithm of the mean fitness. Then, (2.22) may be recast as the gradient-like dynamics:1Δ𝐱=𝜔(𝐱)𝐺(𝐱)𝑊𝐱=𝐺(𝐱)∇𝑉_𝑊(𝐱), (2.36) with |Δ𝐱|=𝟏^∗Δ𝐱=0 as a result of 𝟏^∗𝐺(𝐱)=𝟎^∗. Note ∇𝑉_𝑊(𝐱)^∗Δ𝐱=∇𝑉_𝑊(𝐱)^∗𝐺(𝐱)∇𝑉_𝑊(𝐱)≥0. (2.37) The dynamics (2.36) is of pure gradient-type with respect to the Svirezhev-Shashahani distance metric 𝑑_𝐺(𝐱,𝐱^′); see [5, 6]. For this metric, the distance between 𝐱 and 𝐱^′=𝐱+Δ𝐱 of 𝑆_𝐾 is 𝑑_𝐺(𝐱,𝐱^′)=(Δ𝐱^∗𝐺⁻¹Δ𝐱)^1/2=(^𝐾_𝑘=1𝑥_𝑘⁻¹(Δ𝑥_𝑘)²)^1/2. (2.38) From (2.33) and (2.34), this quantity, which is the length of Δ𝐱 satisfying |Δ𝐱|=0, is also the square-root of half the allelic variance (the standard deviation) in relative fitness.

2.3. Equilibria (Diploid Case)

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

2.3.1. Preliminaries

In contrast with the haploid case, in the diploid situation, the dynamics (2.22) can have more complex equilibrium points, satisfying 𝑤_𝑘(𝐱_eq)=𝑤₁(𝐱_eq),𝑘=2,…,𝐾, and ∑_𝑙𝑥_eq,𝑙=1. To avoid linear manifolds of equilibria, we first assume that all principal minors of 𝑊 are nonsingular and also that the fitnesses of all homozygotes 𝑤_𝑘,𝑘 are positive. In this case, from the Bézout theorem, the number of equilibria is finite and less or equal than the number 2^𝐾−1 of faces of 𝑆_𝐾. Note that the 𝐾 extremal endpoints of 𝑆_𝐾 (0-faces) are always monomorphic fixed points of (2.22).

An instructive example fulfilling these preliminary conditions is 𝑊=𝐼. There are 2^𝐾−1 equilibrium points (the barycenters of the (^𝐾_𝑘+1)𝑘-dimensional faces, 𝑘=0,…,𝐾−1), but only one polymorphic equilibrium which is the barycenter 𝐱_𝐵 of 𝑆_𝐾. This point is the one with minimal fitness and it is unstable. The 0-faces are stable fixed points whereas the barycenters of the 𝑘 faces with 𝑘∈{1,…,𝐾−2} are saddle-points. The simplex 𝑆_𝐾 could be partitioned into pieces each of which being the attraction basins of the stable 0-face states: in contrast with the haploid case, the type of the survivor is not necessarily the one of the fittest; it will depend on the initial condition.

Similar conclusions can be drawn if instead of 𝑊=𝐼, we start with 𝑊=(𝐼−𝐷_𝝀)⁻¹ where 𝝀∶=(𝜆_𝑘,𝑘=1,…,𝐾) satisfies 𝟎≤𝝀<𝟏 (meaning 0≤𝜆_𝑘<1, for all 𝑘). In this case again, there is only one unstable polymorphic equilibrium which is easily seen to be 𝐱_eq=(𝟏−𝝀)/(𝐾−|𝝀|)∈𝑆_𝐾.

Due to its evolutionary interest, we would like now to discuss the opportunity for a polymorphic state to be stable. Under the above assumptions on 𝑊, a unique stable internal (polymorphic) equilibrium state can exist, necessary and sufficient conditions being that (i) there is a unique 𝐳>𝟎 for which 𝑊𝐳=𝟏 and (ii)𝑊 has exactly one strictly positive dominant eigenvalue and at least one strictly negative eigenvalue or else the sequence of principal minors of 𝑊 alternates in sign (see Kingman [7]). If this is the case, the equilibrium polymorphic state is 𝐱_eq=𝐳/|𝐳|. It is stable in the sense that it is a local maximum of the mean fitness 𝜔(𝐱)=𝐱^∗𝑊𝐱, with 𝜔(𝐱_eq)=1/|𝐳|. Since |𝑊|≠0, the linearization of 𝐩(𝐱) at 𝐱_eq has no eigenvalue of modulus 1 and so 𝐱_eq is hyperbolic and/or isolated (see [8]). A stable isolated polymorphic state is asymptotically Lyapunov stable.

Under these additional assumptions therefore on 𝑊, starting from any initial condition in the interior of 𝑆_𝐾, all trajectories will be attracted by 𝐱_eq=𝐳/|𝐳| which is an isolated global maximum of 𝜔(𝐱).

When there is no such unique globally stable polymorphic equilibrium, all trajectories will still converge but perhaps to a local equilibrium state where some alleles get extinct. Which allele and how many alleles are concerned seems to be an unsolved problem in its full generality.

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

We now draw the attention on a particular class of fitness matrices that lead to a polymorphic state, either unstable or stable. We start with the unstable case, extending the above special diagonal case 𝑊=(𝐼−𝐷_𝝀)⁻¹ leading to a unique unstable polymorphic state.

(i) The Unstable Case
Let Λ≥0 be a symmetric irreducible strictly substochastic matrix satisfying Λ𝟏∶=𝐪<𝟏∶ the positive mass-defect vector of Λ is 𝟏−𝐪. Let 𝜆>0. Define the symmetric strictly potential matrix: 𝑊=𝜆⁻¹(𝐼−Λ)⁻¹≥0, with 𝑊⁻¹=𝜆(𝐼−Λ) defining a strictly row-diagonally dominant Stieltjes matrix with the properties [3]: (𝑊⁻¹)_𝑘,𝑘>0,(𝑊⁻¹)_𝑘,𝑙≤0 for 𝑘≠𝑙 and (𝑊⁻¹)_𝑘,𝑘+∑_𝑙≠𝑘(𝑊⁻¹)_𝑘,𝑙>0 for all 𝑘. Then 𝜆𝟏>𝐳=𝑊⁻¹𝟏=𝜆(𝐼−Λ)𝟏=𝜆(𝟏−𝐪)>𝟎. (2.39) The vector 𝐳=𝜆(𝟏−𝐪) is called the equilibrium potential of 𝑊. We have |𝐳|=𝜆(𝐾−|𝐪|). For this class of 𝑊 therefore, 𝑊𝐳=𝟏 admits a positive solution 𝐳.
Conversely, given a nonsingular matrix 𝑊≥0 satisfying 𝑊𝐳=𝟏 for some 𝐳≥𝟎, the matrix Λ=𝐼−𝜆⁻¹𝑊⁻¹ defines a substochastic matrix if and only 𝑊⁻¹ satisfies (𝑊⁻¹)_𝑘,𝑘>0,(𝑊⁻¹)_𝑘,𝑙≤0 for 𝑘≠𝑙 and 𝜆≥max_𝑘(𝑊⁻¹)_𝑘,𝑘. Then, 𝑊⁻¹ is row-diagonally dominant and 𝑊=𝜆⁻¹(𝐼−Λ)⁻¹ is a potential matrix.

Strictly ultrametric (sUm) matrices are special classes of positive-definite and symmetric strictly potential matrices [2]. An sUm matrix 𝑊 is symmetric with nonnegative entries, satisfying (i) 𝑤_𝑘,𝑙≥min{𝑤_𝑘,𝑗,𝑤_𝑗,𝑙}, for all 𝑗,𝑘,𝑙 and (ii) 𝑤_𝑘,𝑘>max_𝑙≠𝑘{𝑤_𝑘,𝑙}, for all 𝑘 (If in condition (ii),≥ is substituted for >, then 𝑊 is simply an ultrametric matrix and this new condition is implied by (i)). If 𝑊 is an sUm matrix, the fitness dynamics will admit an unstable polymorphic equilibrium state, as a result of 𝑊 being positive-definite.

Remark 2.1. Suppose that Λ is substochastic and primitive. Then 𝑊=𝜆⁻¹(𝐼−Λ)⁻¹>0 is an ultrametric matrix. If 𝑉 is the Hadamard reciprocal of 𝑊 with entries 𝑣_𝑘,𝑙=1/𝑤_𝑘,𝑙, it satisfies that 𝑣_𝑘,𝑙≤max{𝑣_𝑘,𝑗,𝑣_𝑗,𝑙}. Therefore 𝑉 is an ultrametric distance associated to the ultrametric potential 𝑊. Tree matrices are ultrametric matrices that are not sUm.

(ii) The Stable Case
To produce a stable equilibrium state from the sUm matrix construction, let 𝑊=𝜆⁻¹(𝐼−Λ)⁻¹ define a symmetric strictly potential matrix as before. Then, there is a 𝐳=𝜆(𝟏−𝐪)>𝟎 for which 𝑊𝐳=𝟏. With 𝛼>1, define 𝛼𝑊∶=|𝐳|𝐽−𝑊. (2.40) With 𝑊=[𝑤_𝑘,𝑙], we have min𝑤_𝑘,𝑙=𝛼/|𝐳|−max𝑤_𝑘,𝑙 and we can choose 𝛼>1 so that 𝑊≥0. We have (𝛼−1)𝑊𝐳=𝟏 and 𝜹^∗𝑊𝜹=−𝜹^∗𝑊𝜹<0 for all 𝜹≠𝟎 satisfying |𝜹|=0 showing that 𝐱_eq∶=𝐳/|𝐳|=(𝟏−𝐪)/(𝐾−|𝐪|) is now a stable polymorphic state for 𝑊. If 𝑊 is an sUm matrix, then clearly 𝑊 satisfies the “anti-sUm” property expressing a fitness domination of the heterozygotes 𝐴_𝑘𝐴_𝑙 over the homozygotes: 𝑤_𝑘,𝑙≤max{𝑤_𝑘,𝑗,𝑤_𝑗,𝑙},∀𝑗,𝑘,𝑙,𝑤_𝑘,𝑘<min_𝑙≠𝑘{𝑤_𝑘,𝑙},∀𝑘. (2.41) Nonnegative symmetric negative-definite anti-sUm fitness matrices 𝑊 will therefore display a stable polymorphic equilibrium state 𝐱_eq. Note that 𝐱^∗_eq𝑊𝐱_eq=(𝛼−1)/|𝐳| is now the maximal value of the mean fitness.

Example 2.2. When 𝐾=2, with 𝑠>−1, ℎ>0, and 𝑠ℎ>−1, let 𝑊=1+𝑠1+𝑠ℎ1+𝑠ℎ1>0 (2.42) define the fitness matrix with selection parameter 𝑠 and dominance ℎ. This 𝑊 is sUm if and only if 𝑠<0 and ℎ>1. The equilibrium state is 𝐱^∗_eq∶=(ℎ/(2ℎ−1);(ℎ−1)/(2ℎ−1)) and it is unstable. This 𝑊 is anti-sUm if and only if 𝑠>0 and ℎ>1. The equilibrium state is the same but it is now stable.

Note that a singular multiplicative fitness matrix of the form 𝑊=𝐰𝐰^∗ cannot be a strictly potential matrix because its determinant |𝑊| is zero.

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

Consider the problem consisting in splitting binarily and recursively the set {1,…,𝐾} till complete reduction to singletons (leaves) which are left behind in the process. For instance, consider the refinement sequence with 𝐾=6 of {1,…,6}≡(123456): (123456)→((23)(1456))→(((2)(3))((16)(45)))→(((1)(6))((4)(5))).

Starting from the left, there are 2𝐾−1=11 blocks of symbols (the total number of nodes in the splitting binary tree with 𝐾 leaves). To each encountered block, numbered from 𝑙=1 to 2𝐾−1, starting from the left, attach a vector 𝐮_𝑙 of size 𝐾 with 𝑖th entry 𝑢_𝑙(𝑖)=1 if symbol 𝑖 is in the block string,0 otherwise. For instance, from the above sequence, 𝐮^∗₁=(1,1,1,1,1,1),𝐮^∗₂=(0,1,1,0,0,0),𝐮^∗₃=(1,0,0,1,1,1),…, and 𝐮^∗₁₀=(0,0,0,1,0,0),𝐮^∗₁₁=(0,0,0,0,1,0). To each such 𝐮_𝑙, attach a number 𝑠_𝑙 which is >0 if |𝐮_𝑙|=1 (the leaves) and ≥0 if |𝐮_𝑙|>1 (the internal nodes, including the root). Then (see [3]), for any choice of 𝑠_𝑙 respecting these constraints 𝑊=^2𝐾−1_𝑙=1𝑠_𝑙𝐮_𝑙𝐮^∗_𝑙≥0 (2.43) is an sUm matrix and any sUm matrix can be represented in this way. Because for the 𝐮_𝑙 corresponding to the leaves 𝑠_𝑙>0, the diagonal terms of 𝑊 are necessarily >0.

Since for each set {1,…,𝐾}, there are 𝑏_𝐾 splitting tree sequences where 𝑏_𝐾 satisfies 𝑏_𝐾=∑^𝐾−1_𝑘=1(^𝐾_𝑘)𝑏_𝑘𝑏_𝐾−𝑘,𝑘≥2,𝑏₁=1, there are many ways to generate an sUm matrix.

Clearly, for each splitting procedure, with 𝜆⁻¹∶=∑^2𝐾−1_𝑙=1𝑠_𝑙>0,𝑊 may be written as 𝑊=𝜆⁻¹(𝐽−^2𝐾−1_𝑙=2𝜆𝑠_𝑙(𝐽−𝐮_𝑙𝐮^∗_𝑙))≥0, (2.44) where the matrices 𝐽−𝐮_𝑙𝐮^∗_𝑙 take values in {0,1}.

Now, with 𝛾⁻¹>0, any matrix of the form 𝑊=(𝜆⁻¹+𝛾⁻¹)𝐽−𝑊=𝛾⁻¹(𝐽+^2𝐾−1_𝑙=2𝛾𝑠_𝑙(𝐽−𝐮_𝑙𝐮^∗_𝑙))≥0 (2.45) is an anti-sUm matrix. Assuming 𝛾=1 and 𝑠_𝑙=𝑠ℎ_𝑙,𝑊 may be written under the form: 𝑊=𝐽+𝑠𝐴 where 𝐴∶=∑^2𝐾−1_𝑙=2ℎ_𝑙(𝐽−𝐮_𝑙𝐮^∗_𝑙). It has at most 2𝐾−1 independent parameters, namely, the ℎ_𝑙,𝑙=2,…,2𝐾−1, and 𝑠. If the ℎ_𝑙 are known, then 𝑊 is a one-parameter family of fitness matrices. (Fitness matrices of the form 𝐽+𝑠𝐴 were considered in [9] in the context of the estimation of 𝑠 problem).

Examples 2.3. Assume also that 𝑠_𝑙=𝑠/(2(𝐾−1)), for all 𝑙, with 𝑠>0 being a selection parameter. Then, defining the (0,1]-valued matrix 1𝐴∶=2(𝐾−1)^2𝐾−1_𝑙=2(𝐽−𝐮_𝑙𝐮^∗_𝑙)>0, (2.46) an anti-sUm matrix of the form 𝑊=𝐽+𝑠𝐴 will admit a stable polymorphic equilibrium. Clearly, 𝐴 itself is anti-sUm. Because of this, there is a 𝐳_𝐴>𝟎 such that 𝐴𝐳_𝐴=𝟏. Thus (𝐽+𝑠𝐴)𝐳_𝐴=(|𝐳_𝐴|+𝑠)𝟏 showing that, with 𝐳=𝐳_𝐴/(|𝐳_𝐴|+𝑠),𝑊𝐳=(𝐽+𝑠𝐴)𝐳=𝟏. We thus have |𝐳|=|𝐳_𝐴|/(|𝐳_𝐴|+𝑠) and so 𝐱_eq=𝐳/|𝐳|=𝐳_𝐴/|𝐳_𝐴|. Furthermore, the equilibrium mean fitness for such models is 𝐱^∗_eq𝑊𝐱_eq=1/|𝐳|=(|𝐳_𝐴|+𝑠)/|𝐳_𝐴|>1.
For the following simple sequence example with 𝐾=4 (four alleles, say A,C,T,G), (1234)→((1)(234))→((24)(3))→((2)(4)) we find 1𝐴=65666635465456453 (2.47) which itself clearly is a symmetric anti-sUm matrix, together with 𝑊=𝐽+𝑠𝐴. For this example, 𝐱^∗_eq=(1/13)(8,1,3,1) and the equilibrium mean fitness is 1+𝑠/13.
Note that taking 𝑠_𝑙=𝑠/(2(𝐾−1)) just for the indices 𝑙 that were initially chosen to satisfy 𝑠_𝑙>0 would also lead to anti-sUm matrices 𝐴>0 and 𝑊=𝐽+𝑠𝐴>0. In this case, the sum (2.46) defining 𝐴 should then be restricted to the indices 𝑙 from 2 to 2𝐾−1 for which 𝑠_𝑙>0. Proceeding in this extreme way for the above simple example, we get the borderline anti-sUm shapes: 1𝐴=6,3444434444344443𝑊=𝑠+223+2𝑠33+2𝑠33+2𝑠33+2𝑠3𝑠+223+2𝑠33+2𝑠33+2𝑠33+2𝑠3𝑠+223+2𝑠33+2𝑠33+2𝑠33+2𝑠3𝑠+22. (2.48) Here, 𝐱^∗_eq is simply the barycenter of the 4-simplex.

Remark 2.4. The ultrametric conditions (2.41) should be compared with the so-called “triangle inequality” conditions (leading to stable polymorphism) pointed out in [10], which read 𝑤_𝑘,𝑙<𝑤_𝑘,𝑗+𝑤_𝑙,𝑗𝑤,∀𝑘≠𝑙andatleastone𝑗≠𝑘,𝑗≠𝑙,_𝑘,𝑙>𝑤_𝑘,𝑘+𝑤_𝑙,𝑙2,∀𝑘≠𝑙. (2.49) It is clear that the class of anti-sUm matrices is a particular subclass of the Lewontin one.

3. Stochastic Evolutionary Dynamics

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

3.1. The Wright-Fisher Model

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

3.1.1. The Model and Its First Properties

Consider an allelic population with constant size 𝑁. In the haploid (diploid) case, 𝑁 is (twice) the number of real individuals. Let 𝐢∶=𝑖_𝑘 and 𝐢^′∶=𝑖^′_𝑘,𝑘=1,…,𝐾 be two vectors of integers quantifying the size of the allelic populations at two consecutive generations 𝑡 and 𝑡+1. We will let 𝑆_𝐾,𝑁={𝐢integers∶|𝐢|=∑^𝐾_𝑘=1𝑖_𝑘=𝑁}. Suppose that the stochastic evolutionary dynamics is now given by a Markov chain whose one-step transition matrix 𝑃 from states 𝐈=𝐢 to 𝐈^′=𝐢^′ is given by the multinomial Wright-Fisher model:ℙ(𝐈^′_𝑡+1=𝐢^′∣𝐈_𝑡=𝐢)=∶𝑃(𝐢,𝐢^′𝑁𝑖)=_′1⋯𝑖^′_𝐾^𝐾𝑘=1𝑝_𝑘(𝐢𝑁)^𝑖′𝑘. (3.1) Therefore, given 𝐈_𝑡=𝐢, to form the next generation, each allele chooses its type independently of the others with probability 𝐩, where 𝐩=(𝑝_𝑘,𝑘=1,…,𝐾) is given either by (2.4) in the haploid case or by (2.22) in the diploid case. In the diploid case, the mechanism 𝐩 is assumed to present a unique polymorphic state, either stable or unstable.

The state-space dimension of the Markov chain governed by (3.1) is 𝑛=(^{𝑁+𝐾−1}_𝐾−1) (the number of compositions of integer 𝑁 into 𝐾 nonnegative parts which is also the number of ways to assign 𝑁 indistinguishable balls into 𝐾 distinguishable boxes). To view 𝑃(𝐢,𝐢^′) as a standard transition matrix of some Markov chain, we need first to order the states 𝐢 and 𝐢^′ in (3.1). Starting from the bottom right corner of 𝑃 states should be arranged in decreasing order when listing the entries of 𝑃 moving up and left along the lines and columns, respectively; or equivalently, starting from the top left corner of 𝑃 states should be arranged in increasing order when moving down and right.

For example, we can order the states 𝐢∈𝑆_𝐾,𝑁 in decreasing order from 𝑛 to 1, as follows. Let 𝐻(𝐢)=(1/𝑁)∑^𝐾_𝑘=1(𝑁+1)^{−(𝑘−1)}𝑖_𝑘 be a base-(𝑁+1) code of the state 𝐢 that will serve as a ranking function. The largest state 𝐢_𝑛 for which 𝐻(𝐢) is maximal (equal to 1) is 𝐢^∗_𝑛∶=(𝑁,0,…,0). Given a state 𝐢, define the subsequent state in decreasing order as 𝜎₋(𝐢)=argmin_{𝐣∶𝐻(𝐣)<𝐻(𝐢)}(𝐻(𝐢)−𝐻(𝐣)). (3.2) Then 𝐢^∗_𝑛=(𝑁,0,…,0),𝐢^∗_𝑛−1=(𝑁−1,1,0,…,0),…,𝐢^∗₁=(0,…,0,𝑁) and to pass from state 𝐢 to the next state 𝜎₋(𝐢) in this decreasing sequence, there is a unique 𝜹_𝐢 with entries in ℤ satisfying |𝜹_𝐢|=0 and such that 𝜎₋(𝐢)=𝐢−𝜹_𝐢. This way to order the 𝑛 states is consistent with the reverse lexicographic order. For instance, if 𝑁=3,𝐾=4, there are (⁶₃)=20 states ordered in decreasing order as follows:3000>2100>2010>2001>1200>1110>1101>1020>1011>1002>0300>0210>0201>0120>0111>0102>0030>0021>0012>0003.

The way the digits are propagated from left to right is clear and the consecutive 𝜹_𝐢 can easily be obtained. Proceeding in this way to order states, 𝑃(𝐢,𝐢^′) is a well-defined conventional object (matrix).

From (3.1), the marginal transition matrix from 𝐢 to 𝐼^′_𝑘=𝑖^′_𝑘 is binomial bin(𝑁,𝑝_𝑘(𝐢/𝑁)) with: 𝑃(𝐢,𝑖^′_𝑘𝑁𝑖)=_′𝑘𝑝_𝑘(𝐢𝑁)(1−𝑝_𝑘(𝐢𝑁))^{𝑁−𝑖′𝑘}. (3.3) Given 𝐈=𝐢, the 𝑘th component 𝐼^′_𝑘 of the updated state is random with 𝔼_𝐢(𝐼^′_𝑘𝑁)=𝑝_𝑘(𝐢𝑁),𝜎²_𝐢(𝐼^′_𝑘𝑁𝑝)=^𝑘(𝐢/𝑁)(1−𝑝_𝑘(𝐢/𝑁))𝑁, (3.4) suggesting that, in the large 𝑁 population limit, the deterministic evolutionary dynamics should be recovered. Indeed, from (3.1), the Laplace-Stieltjes transform of the joint law of 𝐈^′_𝑡+1∣𝐈_𝑡=𝐢 reads 𝔼_𝐢(𝑒^{−∑𝐾𝑘=1𝜆𝑘𝐼′𝑘/𝑁})=(^𝐾_𝑘=1𝑝_𝑘(𝐢𝑁)𝑒^{−𝜆𝑘/𝑁})^𝑁∼_𝑁↑∞𝑒^{−∑𝐾𝑘=1𝜆𝑘𝑝𝑘(𝐢/𝑁)}. (3.5) From the strong law of large numbers therefore, if 𝑖_𝑘∶=⌊𝑁𝑥_𝑘⌋,𝑘=1,…,𝐾, then, given 𝐱=𝑥_𝑘,𝑘=1,…,𝐾, 𝐼^′_𝑘𝑁^a.s._𝑁→_↑∞𝑥^′_𝑘=𝑝_𝑘(𝐱) (3.6) which is (2.22): when the population under study is very large, random fluctuations as modelled by (3.1) can be ignored so that the gene frequencies evolve deterministically.

3.1.2. Fixation Probabilities

Let 𝐞_𝑙 be the 𝐾-null vector except for its 𝑙th entry which is 1. The extremal pure states 𝑆^ex_𝐾,𝑁∶={𝐢^ex_𝑙∶=𝑁𝐞_𝑙,𝑙=1,…,𝐾} are all absorbing for this Markov chain because 𝑝_𝑘(𝐢^ex_𝑙/𝑁)=𝛿_𝑘,𝑙 and, from (3.1), any additional fixed point which 𝐩 could have on the boundary-faces of 𝑆_𝐾 which are not points does not give rise to an absorbing state for 𝑃. Under our assumptions on 𝐩, the chain is not recurrent, rather it is transient. Depending on the initial condition, say 𝐢₀, the chain will necessarily end up in one of the extremal states 𝐢^ex_𝑙, with some fixation probability, say 𝜋_𝑙(𝐢₀), which can be computed as follows. Let 𝝅_𝑙∶=𝜋_𝑙(𝐢),𝐢∈𝑆_𝐾,𝑁 be the harmonic function of the Wright-Fisher Markov chain which is the smallest solution to the boundary problem:(𝐼−𝑃)𝝅_𝑙=𝟎if𝐢∈𝑆_𝐾,𝑁⧵𝑆^ex_𝐾,𝑁,𝝅_𝑙=1(𝐢=𝐢^ex_𝑙)if𝐢∈𝑆^ex_𝐾,𝑁. (3.7) We also have ℙ(𝐈_𝜏=𝑁𝐞_𝑙∣𝐈₀=𝐢₀)=𝜋_𝑙(𝐢₀), (3.8) where 𝜏 (<∞ almost surely) is the random hitting time of 𝑆^ex_𝐾,𝑁 for 𝐈_𝑡 and the 𝜋_𝑙(𝐢₀)s are normalized so as ∑_𝑙𝜋_𝑙(𝐢₀)=1. Thus 𝜋_𝑙(𝐢₀) are the searched probabilities to end up in state 𝐢^ex_𝑙 starting from state 𝐢₀.

From (3.7), 𝝅_𝑙(𝐢₀) is known if 𝐢₀∈𝑆^ex_𝐾,𝑁. The remaining unknown restriction, say 𝝅^𝑄_𝑙, of 𝝅_𝑙=𝜋_𝑙(𝐢₀) to the nonextremal states is easily seen to be𝝅^𝑄_𝑙=(𝐼−𝑄)⁻¹𝐩_𝐢^ex_𝑙,𝐢₀∈𝑆_𝐾,𝑁⧵𝑆^ex_𝐾,𝑁. (3.9) In (3.9), 𝑄 is obtained from 𝑃 after erasing the lines and columns corresponding to all the 𝐾 extremal states and 𝐩_𝐢^ex_𝑙 is the 𝐢^ex_𝑙-column of 𝑃 where the entries corresponding to the extremal states have been deleted. When dealing with the Wright-Fisher model, 𝑄 is a positive matrix and also 𝐩_𝐢^ex_𝑙>𝟎, therefore 𝜋_𝑙(𝐢₀)>0 for all 𝐢₀∈𝑆_𝐾,𝑁⧵𝑆^ex_𝐾,𝑁 and this for each 𝑙: starting from any state 𝐢₀ which is not an extremal state, there is a positive probability to hit any of the extremal states. Fixation of the state 𝐢^ex_𝑙 means extinction of the remaining monomorphic states. It would therefore be of interest to understand the structure of the set 𝒜_𝑙={𝐢₀∶𝜋_𝑙(𝐢₀)>∑_𝑘≠𝑙𝜋_𝑘(𝐢₀)}={𝐢₀∶𝜋_𝑙(𝐢₀)>1/2}, for each 𝑙, which is the stochastic version of the attraction basin of 𝐢^ex_𝑙, especially when 𝑙 is the label of the extremal state with largest fitness 𝑤_𝑘,𝑘. If 𝐢₀∈𝒜_𝑙 indeed, the probability to end up in 𝐢^ex_𝑙 is larger than the probability to end up in any other extremal state.

Unfortunately, the development of the inverse of 𝐼−𝑄 in terms of its adjugate matrix in (3.9) shows that these fixation probabilities have a very complex determinantal alternating structure and the question of identifying 𝒜_𝑙 is very complex.

With 𝝅_𝑙 being the solution (3.9) to the Dirichlet problem with boundary conditions (3.7), the equilibrium measure of the chain therefore is 𝜋_eq(𝐢₀)∶=^𝐾_𝑙=1𝜋_𝑙(𝐢₀)𝛿_𝐢^ex_𝑙, (3.10) which depends on 𝐢₀. Necessarily, one allele will fixate and there is no polymorphic equilibrium state even when dealing with diploid populations. Which allele and with what probability will depend on the initial condition. Thanks to fluctuations, the picture therefore looks very different from the one pertaining to the deterministic theory.

Of importance also is the time it takes to get extinct. It relies on similar techniques. For instance, the expected overall fixation time 𝛼(𝐢₀)∶=𝔼_𝐢0(𝜏) solves the boundary problem: (𝐼−𝑃)𝜶=𝟏,𝐢₀∈𝑆_𝐾,𝑁⧵𝑆^ex_𝐾,𝑁,𝜶=0,𝐢₀∈𝑆^ex_𝐾,𝑁, (3.11) where 𝜶∶=𝛼(𝐢₀),𝐢₀∈𝑆_𝐾,𝑁. The restriction 𝜶^𝑄 of 𝜶 to the nonextremal states therefore is 𝜶^𝑄=(𝐼−𝑄)⁻¹𝟏,𝐢₀∈𝑆_𝐾,𝑁⧵𝑆^ex_𝐾,𝑁. (3.12)

3.1.3. Conditioning 𝐈_𝑡 on Nonfixation

There are four places where questions relative to conditioning on fixation are relevant in this context. (Similar conditioning problems were considered in [11] in the context of the 2-alleles Wright-Fisher diffusion).

(i) Consider the full fixation vector 𝝅_𝑙∶=𝜋_𝑙(𝐢₀). Remove from 𝝅_𝑙 the states 𝐢₀ for which 𝜋_𝑙(𝐢₀)=0 only, keeping the one, 𝐢^ex_𝑙, for which 𝜋_𝑙(𝐢^ex_𝑙)=1. The size of this vector, say 𝝅^𝑅_𝑙, is 𝑛_𝑅=(^{𝑁+𝐾−1}_𝐾−1)−𝐾+1. Let 𝑅 be the corresponding 𝑛_𝑅×𝑛_𝑅-reduced transition matrix obtained from 𝑃 after erasing the lines and columns corresponding to all the 𝐾 monomorphic states except 𝐢^ex_𝑙. The Markov chain 𝐈_𝑡 conditioned to exit in the extremal state 𝐢^ex_𝑙 only admits the stochastic transition matrix: 𝑅_𝑙∶=𝐷_𝝅⁻¹_𝑅𝑙𝑅𝐷_𝝅^𝑅_𝑙. (3.13) It is obtained from 𝑅 after a diagonal Doob transform based on 𝝅^𝑅_𝑙. The chain governed by 𝑅_𝑙 admits a unique absorbing state which is 𝐢^ex_𝑙. The entries of 𝑅_𝑙 are 𝑅_𝑙(𝐢,𝐢^′𝝅)=^𝑅𝑙(𝐢^′)𝝅^𝑅_𝑙(𝐢)𝑅(𝐢,𝐢^′), (3.14) and for this new conditioned Markov chain, transitions to states 𝐢^′ for which 𝝅^𝑅_𝑙(𝐢^′)>𝝅^𝑅_𝑙(𝐢) are favored.

(ii) Let us now consider again the fully reduced transition matrix 𝑄 obtained from 𝑃 after erasing the lines and columns corresponding to all the 𝐾 monomorphic states. The matrix 𝑄 is substochastic and irreducible, with 𝑄𝟏<𝟏. The law of the process corresponding to 𝐈_𝑡 conditioned on avoiding the monomorphic states before 𝑡 evolves as follows: with 𝜏∶=⋀^𝐾_𝑙=1𝜏_𝑙 being the time needed for first hitting one of the extremal states for 𝐈_𝑡, let 𝜋_𝑡(𝐢)=ℙ(𝐈_𝑡=𝐢∣𝜏>𝑡). Then, with 𝝅_𝑡=𝜋_𝑡(𝐢),𝐢∈𝑆_𝐾,𝑁⧵𝑆^ex_𝐾,𝑁, 𝝅^∗_𝑡+1=𝝅^∗_𝑡𝑄𝝅^∗_𝑡𝑄𝟏 (3.15) is the nonlinear Master Equation governing its evolution [12]. The reduced state-space dimension of this Markov chain is 𝑛_𝑄=(^{𝑁+𝐾−1}_𝐾−1)−𝐾 and 𝝅_𝑡→_𝑡↑∞𝝅_∞ where 𝝅_∞ is the left Perron probability eigenvector of 𝑄 associated to the dominant Perron eigenvalue 𝜌_𝑄<1, namely, 𝜌_𝑄𝝅^∗_∞=𝝅^∗_∞𝑄. (𝝅_∞ is called a Yaglom limit (see [13]) or a quasistationary distribution). If the process is started using this limiting quasistationary distribution, it remains in the same state over time and the fixation time 𝜏 is geometrically distributed with success probability 𝜌_𝑄.

(iii) One can define another stochastic process 𝐈_𝑡 which admits the stochastic transition matrix: 𝑄∶=𝐷⁻¹_𝑄𝟏𝑄 again defined on the reduced state-space. For each 𝑡, we have 𝑄(𝐢,𝐢^′)=ℙ(𝐈_𝑡+1=𝐢^′∣𝐈_𝑡=𝐢,𝜏>1), (3.16) and the conditioning on nonfixation occurs at each transition time. This process is an ergodic Markov chain with invariant probability measure solving 𝝅^∗_eq=𝝅^∗_eq𝑄. It has the following closed-form determinantal expression (see [14, Section  6] and [15, page 1559]):𝜋_eq(𝐢)=|(𝐼−𝑄)_[𝐢,𝐢]|∑_𝐢|(𝐼−𝑄)_[𝐢,𝐢]|,𝐢∈𝑆_𝐾,𝑁⧵𝑆^ex_𝐾,𝑁, (3.17) where (𝐼−𝑄)_[𝐢,𝐢] is the submatrix resulting from the deletion of row 𝐢 and column 𝐢 of 𝐼−𝑄. The question as to whether the process governed by 𝑄 is reversible or not arises. Defining the transition matrix of the time-reversed process ^⟵𝑄 by ^⟵𝑄^∗=𝐷_𝝅eq𝑄𝐷⁻¹_𝝅eq, (3.18) it does not hold that ^⟵𝑄=𝑄 and so detailed balance does not hold. Indeed, 𝑄 is similar to the transition matrix of an ergodic Wright-Fisher model and Wright-Fisher chains are not reversible.

(iv) If we condition on nonfixation in the remote future (see [16] for additional details), we get a Markov chain whose stochastic transition matrix is 𝑄=𝜌_𝑄⁻¹𝐷_𝝍⁻¹_∞𝑄𝐷_𝝍∞. (3.19) Here 𝝍_∞ is the positive right Perron eigenvector of 𝑄 associated to the Perron eigenvalue 𝜌_𝑄<1 satisfying 𝜌_𝑄𝝍_∞=𝑄𝝍_∞. This vector can be chosen so that ∑_𝑘𝜋_∞,𝑘𝜓_∞,𝑘=1, where 𝝅_∞ is again the left Perron probability eigenvector of 𝑄 associated to 𝜌_𝑄<1 (see [17]). With 𝝅_𝑡(⋅)=lim_𝑠↑∞ℙ(𝐈_𝑡=⋅∣𝜏>𝑡+𝑠), we have 𝝅^∗_𝑡+1=𝝅^∗_𝑡𝑄. The process governed by 𝑄 is an ergodic Markov chain whose invariant probability measure is 𝝅_∞=𝝅_∞∘𝝍_∞, the Schur product of 𝝅_∞ and 𝝍_∞ with 𝑘th entry 𝝅_∞,𝑘=𝜋_∞,𝑘𝜓_∞,𝑘.

For the last three conditionings, it is difficult to extract some information on the limiting distribution, either 𝝅_∞ or 𝝅_eq or 𝝅_∞, respectively. This is because it would suppose to solve the eigenvalue problems explicitly which is out of reach, at least theoretically. However, assuming a diploid population with a polymorphic equilibrium state 𝐱_eq for 𝐩, we expect that these distributions will present a global (local) maximum (minimum) near 𝐱_eq if 𝐱_eq is stable (unstable). These limiting distributions should be more sharply peaked around the extremum if we consider the conditioning (iv) compared to (ii) because, the latter conditioning being more stringent than the former, it should charge more heavily the sample paths that stay away from the monomorphic states.

Finally, we would like to stress that all these considerations are also relevant in the context of another fundamental stochastic model arising in the context of evolutionary genetics. We will give some elements of how to proceed with this model presenting very different properties.

3.2. The 𝐾-Alleles Moran Model

We now focus on the Moran model.

3.2.1. The Multiallelic Moran Model

Let 𝛼,𝛽∈{1,…,𝐾}. In the Moran version of the stochastic evolution, given 𝐈_𝑡=𝐈=𝐢, the only accessible values of 𝐈^′ are the neighboring states: 𝐢^′_𝛼,𝛽∶=𝐢+𝐝_𝛼,𝛽 where 𝐝^∗_𝛼,𝛽∶=(0,…,0,−1,0,…,1,0,…,0). Here −1 is in position 𝛼 and 1 in position 𝛽≠𝛼 corresponding to the transfer of an individual from the box 𝛼 (if nonempty) to the box 𝛽. With 𝑛(𝐢)=^#{𝑘∶𝑖_𝑘>0} being the number of nonempty entries of 𝐢, there are 𝑛(𝐢)(𝐾−1)≤𝐾(𝐾−1) accessible states from 𝐢. The Moran stochastic evolutionary dynamics is now given by a Markov chain whose one-step transition matrix 𝑃 from states 𝐈=𝐢 to 𝐈^′=𝐢^′ isℙ(𝐈_𝑡+1=𝐢^′∣𝐈_𝑡=𝐢)=0if𝐢^′≠𝐢^′_𝛼,𝛽,ℙ(𝐈_𝑡+1=𝐢^′_𝛼,𝛽∣𝐈_𝑡=𝐢)=∶𝑃(𝐢,𝐢^′_𝛼,𝛽𝑖)=^𝛼𝑁𝑝_𝛽(𝐢𝑁), (3.20) where 𝐩=(𝑝_𝛽(𝐢/𝑁),𝛽=1,…,𝐾) is given either by (2.4) in the haploid case or by (2.22) in the diploid case.

Summing 𝑃(𝐢,𝐢^′_𝛼,𝛽) over 𝛼,𝛽, 𝛽≠𝛼 in (3.20), we get the holding probability ℙ(𝐈_𝑡+1=𝐢∣𝐈_𝑡=𝐢)=1−_{𝛼,𝛽∶𝛽≠𝛼}𝑖_𝛼𝑁𝑝_𝛽(𝐢𝑁)=_𝛼𝑖_𝛼𝑁𝑝_𝛼(𝐢𝑁), (3.21) completing the characterization of the 𝐾-alleles Moran model. Under our assumptions on 𝐩, the holding probabilities are equal to 1 only for the extremal states 𝐢∈𝑆^ex_𝐾,𝑁 which are therefore the only absorbing states of the Moran chain, just like for the Wright-Fisher model. The drift at state 𝐢 is𝔼(𝐈_𝑡+1−𝐈_𝑡∣𝐈_𝑡=𝐢)=_𝛼𝑖_𝛼𝑁_𝛽≠𝛼𝑝_𝛽(𝐢𝑁)𝐝_𝛼,𝛽. (3.22)

Let us compute the Laplace-Stieltjes Transform of 𝐈^′ in the context of a Moran model. Omitting the argument 𝐢/𝑁 in 𝑝_𝛽(𝐢/𝑁), we get the factorized form: 𝔼_𝐢(𝑒^{−∑𝑘𝜆𝑘𝐼′𝑘})=_{𝛼,𝛽∶𝛼≠𝛽}𝑒^{−∑𝑘𝜆𝑘𝑖′𝛼,𝛽(𝑘)}𝑃(𝐢,𝐢^′_𝛼,𝛽)+𝑒^{−∑𝑘𝜆𝑘𝑖𝑘}_𝛽𝑖_𝛽𝑁𝑝_𝛽=𝑒^{−∑𝑘𝜆𝑘𝑖𝑘}(_{𝛼,𝛽∶𝛼≠𝛽}𝑒^{−∑𝑘𝜆𝑘𝐝𝛼,𝛽(𝑘)}𝑃(𝐢,𝐢^′_𝛼,𝛽)+_𝛽𝑖_𝛽𝑁𝑝_𝛽)=𝑒^{−∑𝑘𝜆𝑘𝑖𝑘}(_{𝛼,𝛽∶𝛼≠𝛽}𝑒^{−(𝜆𝛽−𝜆𝛼)}𝑖_𝛼𝑁𝑝_𝛽+_𝛽𝑖_𝛽𝑁𝑝_𝛽)=𝑒^{−∑𝑘𝜆𝑘𝑖𝑘}(_𝛽𝑒^−𝜆𝛽𝑝_{𝛽𝛼≠𝛽}𝑖_𝛼𝑁𝑒^𝜆𝛼+_𝛽𝑖_𝛽𝑁𝑝_𝛽)=𝑒^{−∑𝑘𝜆𝑘𝑖𝑘}(_𝛽𝑒^−𝜆𝛽𝑝_𝛽(_𝛼𝑖_𝛼𝑁𝑒^𝜆𝛼−𝑖_𝛽𝑁𝑒^𝜆𝛽)+_𝛽𝑖_𝛽𝑁𝑝_𝛽)=(𝑒^{−∑𝑘𝜆𝑘𝑖𝑘})(_𝛼𝑖_𝛼𝑁𝑒^𝜆𝛼)(_𝛽𝑒^−𝜆𝛽𝑝_𝛽). (3.23) Putting 𝜆_𝑙=0 if 𝑙≠𝑘, the 𝑘th marginal reads𝔼_𝐢(𝑒^{−𝜆𝑘𝐼′𝑘})=𝑒^{−𝜆𝑘𝑖𝑘}𝑖(1−^𝑘𝑁+𝑒^𝜆𝑘𝑖_𝑘𝑁)(1−𝑝_𝑘+𝑒^−𝜆𝑘𝑝_𝑘) (3.24) showing that 𝐼_𝑘(𝑡) is of the random walk type. Indeed, we get ℙ_𝐢(𝐼^′_𝑘=𝑖^′_𝑘)=0 if 𝑖^′_𝑘≠𝑖_𝑘±1 or 𝑖^′_𝑘≠𝑖_𝑘 and ℙ_𝐢(𝐼^′_𝑘=𝑖_𝑘𝑖)=(1−^𝑘𝑁)(1−𝑝_𝑘𝑖)+^𝑘𝑁𝑝_𝑘,ℙ_𝐢(𝐼^′_𝑘=𝑖_𝑘𝑖+1)=(1−^𝑘𝑁)𝑝_𝑘=_𝑙≠𝑘𝑃(𝐢,𝐢^′_𝑙,𝑘ℙ),_𝐢(𝐼^′_𝑘=𝑖_𝑘𝑖−1)=^𝑘𝑁(1−𝑝_𝑘)=_𝑙≠𝑘𝑃(𝐢,𝐢^′_𝑘,𝑙). (3.25) We have 𝔼_𝐢(𝐼^′_𝑘𝑁𝑖)=^𝑘𝑁+1𝑁(𝑝_𝑘−𝑖_𝑘𝑁);𝜎²_𝐢(𝐼^′_𝑘𝑁1)=𝑁²(𝑖_𝑘𝑁𝑖(1−^𝑘𝑁)+𝑝_𝑘(1−𝑝_𝑘)). (3.26) There is only a small correction (of order 𝑁⁻¹) of the updated mean to its current value and fluctuations around the mean are small too (of order 𝑁⁻¹). The evolution process is very slow.

3.2.2. Fixation Probabilities

As a random walk model, the Moran model has a much simpler transition matrix 𝑃 of the Jacobi type. The equilibrium measure of the chain again is𝜋_eq∶=^𝐾_𝑙=1𝜋_𝑙(𝐢₀)𝛿_𝐢^ex_𝑙, (3.27) where 𝝅_𝑙 again solves the Dirichlet problem (3.7) but with this new simpler Jacobi 𝑃. For the Moran model, the explicit expression (3.9) of the fixation probability simplifies a little bit because 𝑝_𝐢^ex_𝑙(𝐢)≠0 only for the 𝐾−1 neighboring states of 𝐢^ex_𝑙 that is {𝐢∶𝐢+𝐝_𝛼,𝛽=𝐢^ex_𝑙} for some 𝐝_𝛼,𝛽.

When 𝐾=2 (2 alleles), the random walk transition probabilities (𝑝₁+𝑝₂=1) ℙ(𝐈_𝑡+1=𝐢^′_1,2∣𝐈_𝑡1𝑖=𝐢)=∶𝑃(𝐢,𝐢+−1)=¹𝑁𝑝₂(𝐢𝑁),ℙ(𝐈_𝑡+1=𝐢^′_2,1∣𝐈_𝑡1𝑖=𝐢)=∶𝑃(𝐢,𝐢+−1)=²𝑁𝑝₁(𝐢𝑁) (3.28) are the probabilities that the first component 𝐼_𝑡,1 of 𝐈_𝑡 moves down and up by one unit, respectively. In this case, the Dirichlet problem giving the fixation probabilities solves explicitly. With 𝜙(𝑖₀)=1+∑^𝑖0−1_𝑖=1∏^𝑖_𝑖1=1(𝑖₁𝑝₂(𝐢/𝑁)/𝑖₂𝑝₁(𝐢/𝑁)) (𝜙(0)=0) being the harmonic function of the 2-alleles chain, we easily get that 𝜋₁(𝑖₀,𝑁−𝑖₀)=𝜙(𝑖₀)𝜙(𝑁) (3.29) is the probability that the extremal state 𝐢^ex₁=(𝑁,0) is reached given 𝐢₀=(𝑖₀,𝑁−𝑖₀). Assuming a model with multiplicative fitnesses: 𝑝_𝛼(𝐢/𝑁)=(𝑖_𝛼/𝑁)(𝑤_𝛼/𝜔(𝐢/𝑁)), then 𝜙 takes the simple form (𝑖₁+𝑖₂=𝑁) 𝜙(𝑖₀)=1+^{𝑖0−1𝑖}_{𝑖=1𝑖1=1}𝑖₁𝑝₂(𝐢/𝑁)𝑖₂𝑝₁(𝐢/𝑁)=1+^𝑖0−1_𝑖=1(𝑤₂𝑤₁)^𝑖 (3.30) showing (see [1, page 109]) that𝜋₁(𝑖₀,𝑁−𝑖₀)=1−(𝑤₂/𝑤₁)^𝑖01−(𝑤₂/𝑤₁)^𝑁. (3.31) Assuming 𝑤₁=1+𝑠/𝑁 and 𝑤₂=1, putting 𝑖₀=[𝑁𝑥₀], for large 𝑁, we get [18] 𝜋₁(𝑁𝑥₀,𝑁−𝑖₀)∼1−𝑒^−𝑠𝑥01−𝑒^−𝑠. (3.32)

In the general fitness case 𝜙(𝑖₀)=1+^{𝑖0−1𝑖}_{𝑖=1𝑖1=1}𝑖₁𝑝₂(𝐢/𝑁)𝑖₂𝑝₁(𝐢/𝑁)=1+^{𝑖0−1𝑖}_{𝑖=1𝑖1=1}(𝑊𝐢/𝑁)₂(𝑊𝐢/𝑁)₁, (3.33) where, as usual, (𝑊𝐢/𝑁)_𝑘=∑²_𝑙=1𝑤_𝑘,𝑙𝑖_𝑙/𝑁,𝑘,𝑙=1,2, leading to (𝑊𝐢𝑁)₁=𝑤_1,2+(𝑤_1,1−𝑤_1,2)𝑖₁𝑁,(𝑊𝐢𝑁)₂=𝑤_2,2+(𝑤_2,1−𝑤_2,2)𝑖₁𝑁. (3.34)

This 2-alleles exact solution can be used in the full diploid 𝐾-alleles Moran case with multiplicative fitnesses. Indeed, from this, with 𝐢₀=(𝑖₁,…,𝑖_𝐾), the fixation probability of 𝐴_𝑙 can be conjectured to be approximated qualitatively by𝜋_𝑙(𝐢^ex_𝑙)=1,𝜋_𝑙(𝐢₀)=1−((∑_𝑘≠𝑙𝑖_𝑘𝑤_𝑘/(𝑁−𝑖_𝑙))/𝑤_𝑙)^𝑖𝑙1−((∑_𝑘≠𝑙𝑖_𝑘𝑤_𝑘/(𝑁−𝑖_𝑙))/𝑤_𝑙)^𝑁if𝐢₀≠𝐢^ex_𝑙. (3.35) This can be justified as follows: mark one particular box with size 𝑖_𝑙, corresponding to the allele 𝐴_𝑙 with fitness 𝑤_𝑙. Then clump the 𝐾−1 remaining boxes into a single box with size 𝑁−𝑖_𝑙, corresponding to a fictitious allele with average fitness ∑_𝑘≠𝑙𝑖_𝑘𝑤_𝑘/(𝑁−𝑖_𝑙). We are left with a 2-alleles Moran multiplicative fitness model for which (from the 2-alleles exact solution) the fixation probability of 𝐴_𝑙 is given by (3.35). This formula constitutes sort of a mean field approximation to the full Dirichlet problem associated to the Moran model.

Assuming 𝑤_𝑘∼1+𝑠_𝑘/𝑁, a Kimura-like approximation of (3.35) would lead for large 𝑁 to 𝜋_𝑙(𝑁𝐱₀)∼1−𝑒^{(𝑥𝑙/(1−𝑥𝑙))∑𝑘≠𝑙𝑥𝑘𝑠𝑘}1−𝑒^{(1/(1−𝑥𝑙))∑𝑘≠𝑙𝑥𝑘𝑠𝑘}, (3.36) where 𝐱₀=(𝑥₁,…,𝑥_𝐾) is now a point of the continuous 𝐾-simplex different from the 𝑙-unit vector 𝐞_𝑙∶=(0,…,0,1,0,…,0).

From (3.35), when the fitness of allele 𝐴_𝑙 is large (small), compared to the average fitness of the remaining alleles, then the fixation probability of 𝐴_𝑙 gets close to 1 (resp., to 0). Note also that the larger 𝑖_𝑙 is, the larger the fixation probability is.

As required also, for all 𝑘≠𝑙,𝜋_𝑙(𝐢^ex_𝑘)=(1−(𝑤_𝑘/𝑤_𝑙)⁰)/(1−(𝑤_𝑘/𝑤_𝑙)^𝑁)=0. As another particular initial configuration case, suppose that we start from the 2-alleles type state: 𝐢₀=∶𝐢₀(𝑘,𝑙)=(0,…,0,𝑁−1,0,…,0,1,0,…,0) where the 1 is in position 𝑙 (i.e., 𝑖_𝑙=1) and the entry 𝑁−1 in position 𝑘≠𝑙 (i.e., 𝑖_𝑘=𝑁−1). Although for this choice of the initial state, the fixation of 𝐴_𝑘 is very likely, there still is a positive probability that allele 𝐴_𝑙 gets fixed which is seen to be from (3.35): 𝜋_𝑙(𝐢₀(𝑘,𝑙))=1−𝑤_𝑘/𝑤_𝑙1−(𝑤_𝑘/𝑤_𝑙)^𝑁, (3.37) depending only on the relative fitnesses of 𝐴_𝑙 and 𝐴_𝑘 (see [19, 20] for a similar expression). As 𝑁 gets large, this probability gets close to 1−𝑤_𝑘/𝑤_𝑙 if 𝑤_𝑘<𝑤_𝑙 and close to (𝑤_𝑙/𝑤_𝑘)^𝑁−1∼0 if 𝑤_𝑘>𝑤_𝑙.

Conversely, assuming the 2-alleles type state 𝐢₀ to be defined by 𝑖_𝑙=𝑁−1 and 𝑖_𝑘=1, 𝜋_𝑙(𝐢₀)=1−(𝑤_𝑘/𝑤_𝑙)^𝑁−11−(𝑤_𝑘/𝑤_𝑙)^𝑁, (3.38) which, as required, gets close to 1 as 𝑁 gets large if 𝑤_𝑘<𝑤_𝑙 and close to 𝑤_𝑙/𝑤_𝑘 if 𝑤_𝑘>𝑤_𝑙.

3.2.3. Conditioning 𝐈_𝑡 on Nonfixation

The conditioning developments discussed for the Wright-Fisher model are also relevant in the Moran model context substituting the 𝑃 of Moran for the 𝑃 of Wright-Fisher. Let us revisit condition (iii).

(iii) With 𝑄 being now the reduced substochastic matrix of the full Moran transition matrix defined in (3.20), consider the stochastic process 𝐈_𝑡 with stochastic transition matrix: 𝑄∶=𝐷⁻¹_𝑄𝟏𝑄 defined on the reduced state-space with dimension 𝑛_𝑄. The process 𝐈_𝑡 is again an ergodic Markov chain with invariant probability measure solving 𝝅^∗_eq=𝝅^∗_eq𝑄. With ^⟵𝑄 being the transition matrix of the time-reversed process it now holds that ^⟵𝑄=𝑄 and so detailed balance holds when dealing with the Moran case (see also [21]). This will be proved if we can exhibit an equilibrium probability measure 𝝅_eq such that𝜋_eq(𝐢)𝑄(𝐢,𝐢^′)=𝑄(𝐢^′,𝐢)𝜋_eq(𝐢^′), (3.39) for all neighboring state 𝐢^′=𝐢+𝐝_𝐢 of 𝐢.

With 𝐣=(𝑗₁,…,𝑗_𝐾) any terminal state, suppose that we want to use (3.39) to compute 𝜋_eq(𝐣) starting from the smallest available state in the system which is 𝐣₀=(0,…,0,1,𝑁−1). This is possible because 𝜋_eq(𝐣) may be represented as𝜋_eq(𝐣)=𝜋_eq(𝐣₀)^{𝐣−𝐝𝐾,1}_𝐢=𝐣0𝑄(𝐢,𝐢^′)𝑄(𝐢^′,,𝐢) (3.40) where 𝜋_eq(𝐣₀) can be chosen so that ∑_𝐣𝜋_eq(𝐣)=1. Let us give some details.

(i) Note first, by reversing path, that for two consecutive states (𝐢,𝐢^′), the ratio 𝑄(𝐢,𝐢^′)/𝑄(𝐢^′,𝐢) can be computed. We have𝑄(𝐢,𝐢^′)𝑄(𝐢^′=,𝐢)(𝑄𝟏)_𝐢′(𝑄𝟏)_𝐢𝑄(𝐢,𝐢^′)𝑄(𝐢^′,,𝐢) (3.41) where 𝐢^′=𝐢+𝐝_𝐢 for some 𝐝_𝐢 of the form 𝐝_𝑘,𝑙 and therefore𝑄(𝐢,𝐢^′𝑖)=^𝑘𝑁𝑝_𝑙(𝐢𝑁),𝑄(𝐢^′,𝐢)=𝑄(𝐢^′,𝐢^′+𝐝_𝑙,𝑘𝑖)=^𝑙+1𝑁𝑝_𝑘(𝐢+𝐝_𝑘,𝑙𝑁). (3.42) Clearly, from such a structure of the entries of 𝑄, for each possible transition 𝐢→𝐢^′, the ratio 𝑄(𝐢,𝐢^′)/𝑄(𝐢^′,𝐢) will only depend separately on a ratio involving the terminal and the initial states 𝐢^′ and 𝐢 (the detailed balance condition holds).

(ii) Second, the sequence of 𝐢-s in (3.40) is governed by the following path starting from 𝐣₀ and ending up in the target state 𝐣: we can use the following sequence of 𝐝_𝑘,𝑙s∶(𝐝⋯^𝑗_{𝐾,𝐾−1}^𝐾−1−1𝐝_{𝐾,𝐾−1})(𝐝⋯^𝑗_{𝐾,𝐾−2}^𝐾−2𝐝_{𝐾,𝐾−1})⋯(𝐝⋯^𝑗_𝐾,1¹𝐝_𝐾,1), (3.43) filling up successively the entries of 𝐣 to the left of the last entry of 𝐣₀ by using the 𝑁−1 individuals of the reservoir state 𝐣₀=(0,…,0,1,𝑁−1). By doing so, each intermediate state 𝐢 is separated from the next 𝐢^′ by some clearly identified 𝐝_𝐢, and after evaluating the probability ratio 𝑄(𝐢,𝐢^′)/𝑄(𝐢^′,𝐢) for each consecutive states of this sequence, we are done. Because there exists a probability distribution 𝜋_eq(𝐣) such that the reversibility identity (3.39) holds, then this Moran process is reversible with (3.40) as its stationary distribution. Using (3.42), this constitutes an exact explicit product-form formula.

4. Concluding Remarks

This paper studies a classical population genetic model describing a one-locus multiallelic population subject to natural selection, random mating and then random genetic drifts as from the Wright-Fisher and Moran models. Although the two-alleles version of this model is fairly well studied and understood, this is not so much the case of the multiallelic one, especially in the discrete-time context which we adopt here. Let us summarize our results emphasizing that the ones which we believe are new.

Considering first the deterministic updating mechanisms driven by selection, we underline that it has the form of a nonlinear Master equation suggesting that it is possible to construct an underlying Markov process governed by this Master equation. We briefly and intuitively supply such a construction.

In the diploid context, we pay attention on a class of fitness matrices that leads to polymorphism. Would the equilibrium polymorphic state be unstable, we suggest that the class of potential matrices constitute a large such admissible class of fitness matrices. It contains the class of strictly ultrametric matrices which therefore deserves some interest. To the best of our knowledge, there is no discussion of such fitness models in the population genetics context. Would the polymorphic state be stable, we derive a related class of fitness matrices leading to a definite-negative mean fitness quadratic form. Some simple examples are supplied and detailed.

The last section is devoted to the stochastic version of these considerations taking into account an additional important driving source of evolution, namely, the random genetic drift. When driven by selection only and in particular in the absence of mutations, the multiallelic Wright-Fisher model is a transient Markov chain whose absorbing states are the monomorphic states. We give an expression for the fixation probabilities for this process. Then, we develop four conditioning problems: conditioning on fixating in a given monomorphic state, conditioning on avoiding the extremal states before the current instant, conditioning on nonfixation at each transition time, and conditioning on avoiding the extremal states in the remote future. Finally we run into similar considerations but for the Moran model. When dealing with the fixation probabilities in this Moran context, we suggest a mean-field approximation of these probabilities which is based on a well-known explicit formula for the 2-alleles case. It concerns the case of multiplicative fitnesses only. Finally, we consider the Moran model conditioned on nonfixation at each transition time. We exploit the reversible character of this process to derive an explicit product formula for its invariant probability measure.