Discrete ancestral problems arising in population genetics are investigated. In the neutral case, the duality concept has been proved of particular interest in the understanding of backward in time ancestral process from the forward in time branching population dynamics. We show that duality formulae still are of great use when considering discrete nonneutral Wright-Fisher models. This concerns a large class of nonneutral models with completely monotone (CM) bias probabilities. We show that most classical bias probabilities used in the genetics literature fall within this CM class or are amenable to it through some “reciprocal mechanism" which we define. Next, using elementary algebra on CM functions, some suggested novel evolutionary mechanisms of potential interest are introduced and discussed.

1. Introduction

Forward evolution of large populations in genetics has a long history, starting in the 1920s; it is closely attached to the names of R. A. Fisher and S. Wright; see Nagylaki [1] for historical commentaries and on the role played by the French geneticist G. Malécot, starting shortly before the second world war. The book of Ewens [2] is an excellent modern presentation of the current mathematical theory. Coalescent theory is the corresponding backward problem obtained while running the forward evolution processes backward in time. It was discovered independently by several researchers in the 1980s, but definitive formalization is commonly attributed to Kingman [3]. Major contributions to the development of coalescent theory were made (among others) by P. Donnelly, R. Griffiths, R. Hudson, F. Tajima, and S. Tavaré (see the course of Tavaré in Saint-Flour [4] for a review). It included incorporating variations in population size, mutation, recombination, selection, and so forth. In (1999), Pitman [5] and Sagitov [6], independently, introduced coalescent processes with multiple collisions of ancestral lineages. Shortly later, the full class of exchangeable coalescent processes with simultaneous multiple mergers of ancestral lineages was discovered by Möhle and Sagitov [7], and Schweinsberg [8]. All these recent developments and improvements concern chiefly the discrete neutral case and its various scaling limits in continuous time and/or space. As was shown by Möhle [9, 10], neutral forward and backward theories learn much from one another by using a concept of duality introduced by Liggett [11]. Backward theory in the presence of mutations in the forward process is well understood, as it requires the study of a marked Kingman's tree (see Tavaré [4], for a review). In the works of Neuhauser and Krone [12], there is also some use of the duality concept in an attempt to understand the genealogies of a Wright-Fisher diffusion (as a limit of a discrete Wright-Fisher model) presenting a selection mechanism; this led these authors to the idea of the ancestral selection graph extending Kingman's coalescent tree of the neutral theory; see also Huillet [13] for related objectives in the context of Wright-Fisher diffusions with and without drifts. There is, therefore, some evidence that the concept of duality could help one understand the backward theory even in nonneutral situations when various evolutionary forces are the causes of deviation to neutrality (see Crow and Kimura [14], Maruyama [15], Gillespie [16], and Ewens [2], for a discussion on various models of utmost interest in population genetics).

In this note, we focus on discrete nonneutral Wright-Fisher (WF) models and on the conditions on the bias probabilities under which forward branching population dynamics is directly amenable to a dual discrete ancestral coalescent. We emphasize that duality formulae still are of great use when considering discrete nonneutral Wright-Fisher models, at least for specific deviation forces to neutrality. It is shown that it concerns a large class of nonneutral models involving completely monotone bias probabilities. Several classical examples are supplied in the light of complete monotonicity. In the process leading us to focus on these peculiar bias models, some unsuspected evolutionary mechanisms of potential interest are introduced and discussed, as suggested by elementary algebra on completely monotone functions. We emphasize that the relevance of these novel bias mechanisms in Biology seems to deserve additional work and confrontation with real-world problems is urged for to pinpoint their biological significance.

We will finally briefly outline the content of this manuscript. Section 2 is designed to fix the background and ideas. We introduce some basic facts about the discrete-time forward (Section 2.2) and backward processes (Section 2.3) arising from exchangeable reproduction laws (Section 2.1). In Section 2.4, we introduce a concept of duality and briefly recall its relevance to the study of the neutral case problem. The basic question we address in subsequent sections is whether this notion of duality still makes sense in nonneutral situations. We start supplying important nonneutral examples in Section 3. In Section 4, we show that duality does indeed make sense in the framework of discrete nonneutral Wright-Fisher models, but only for the class of completely-monotone state-dependent transition frequencies. In Section 5, we show that most nonneutrality mechanisms used in the literature fall within this class, or are amenable to it via some “reciprocal transformation,” starting with elementary mechanisms and ending up with more complex ones. In Section 6, we show that duality can be used in nonneutral situations to compute the extinction probabilities (invariant measure) of the dual backward ancestral processif one knows the invariant measure (resp., extinction probabilities) of the forward branching process.

2. Discrete-Time Neutral Coalescent

In this section, to fix the background and notations, we review some well-known facts from the cited literature.

2.1. Exchangeable Neutral Population Models: Reproduction Laws (the Cannings Model [17])

Consider a population with nonoverlapping generations 𝑟. Assume the population size is constant, say 𝑛 (𝑛 individuals (or genes)) over generations. Assume the random reproduction law at generation 0 is 𝝂𝑛=(𝜈1,𝑛,,𝜈𝑛,𝑛), satisfying𝑛𝑚=1𝜈𝑚,𝑛=𝑛.(2.1) Here, 𝜈𝑚,𝑛 is the number of offspring of gene 𝑚. We avoid the trivial case: 𝜈𝑚,𝑛=1,𝑚=1,,𝑛. One iterates the reproduction over generations, while imposing the following additional assumptions:

(i)exchangeability: (𝜈1,𝑛,,𝜈𝑛,𝑛)𝑑=(𝜈𝜎(1),𝑛,,𝜈𝜎(𝑛),𝑛), for all permutations 𝜎𝒮𝑛;(ii)time-homogeneity: reproduction laws are independent and identically distributed (iid) at each generation 𝑟. This model, therefore, consists of a conservative-conditioned branching Galton-Watson process in [𝑛], where [𝑛]={0,1,,𝑛} (see Karlin-McGregor [18]).

Famous reproduction laws are as follows.

Example 2.1. The multinomial Dirichlet model: 𝝂𝑛𝑑 Multin-Dirichlet (𝑛;𝜃), where 𝜃>0 is a disorder parameter. With 𝐤𝑛=(𝑘1,,𝑘𝑛), 𝝂𝑛 admits the following joint exchangeable distribution on the simplex |𝐤𝑛|=𝑛𝑚=1𝑘𝑚=𝑛: 𝜃𝝂𝑛=𝐤𝑛=𝑛![𝑛𝜃]𝑛𝑛𝑚=1[𝜃]𝑘𝑚𝑘𝑚!,(2.2) where [𝜃]𝑘=𝜃(𝜃+1)(𝜃+𝑘1) is the rising factorial of 𝜃. This distribution can be obtained by conditioning 𝑛 independent mean 1 Pòlya distributed random variables 𝝃𝑛=(𝜉1,,𝜉𝑛) on summing to 𝑛, that is to say, 𝝂𝑛𝑑=(𝝃𝑛|𝝃𝑛|=𝑛), where 𝜃𝜉1==𝑘[𝜃]𝑘𝑘!(1+𝜃)𝑘𝜃1+𝜃𝜃,𝑘.(2.3) When 𝜃, this distribution reduces to the Wright-Fisher model for which 𝝂𝑛𝑑 Multin(𝑛;1/𝑛,,1/𝑛). Indeed, 𝝂𝑛 admits the following joint exchangeable multinomial distribution on the simplex |𝐤𝑛|=𝑛: 𝝂𝑛=𝐤𝑛=𝑛!𝑛𝑛𝑛𝑚=1𝑘𝑚!.(2.4) This distribution can be obtained by conditioning 𝑛 independent mean 1 Poisson distributed random variables 𝝃𝑛=(𝜉1,,𝜉𝑛) on summing to 𝑛𝝂𝑛𝑑=(𝝃𝑛|𝝃𝑛|=𝑛). When 𝑛 is large, using Stirling formula, 𝑛!2𝜋𝑛𝑛+1/2𝑒𝑛; it follows that 𝝂𝑛𝑑𝑛𝝃 with joint finite-dimensional law: (𝝃𝑛=𝐤𝑛)=𝑛𝑚=1(𝑒1/𝑘𝑚!)=(𝑒𝑛/𝑛𝑚=1𝑘𝑚!) on 𝑛. Thanks to the product form of all finite-dimensional laws of 𝝃, we get an asymptotic independence property of 𝝂𝑛.

Example 2.2. In the Moran model, 𝝂𝑛𝑑 random permutation of (2,0,1,,1). In such a model, only one new gene per generation may come to life, at the expense of the simultaneous disappearance of some other genes.

2.2. Forward in Time Branching Process

Take a subsample of size 𝑚 from [𝑛]={0,1,,𝑛} at generation 0. Let 𝑁𝑟(𝑚)=#ospringatgeneration𝑟+,forwardintime.(2.5) This sibship process is a discrete-time homogeneous Markov chain, with transition probability 𝑁𝑟+1(𝑚)=𝑘𝑁𝑟𝜈(𝑚)=𝑘=1,𝑛++𝜈𝑘,𝑛=𝑘.(2.6) It is a martingale, with state space {0,,𝑛}, initial state 𝑚, absorbing states {0,𝑛}, and transient states {1,,𝑛1}. The first hitting time of boundaries {0,𝑛}, which is 𝜏(𝑚)=𝜏{0}𝜏(𝑚){𝑛}(𝑚), is finite with probability 1 and has finite mean. Omitting reference to any specific initial condition 𝑚, the process (𝑁𝑟;𝑟) has the transition matrix Π𝑛 with entries Π𝑛(𝑘,𝑘)=(𝜈1,𝑛++𝜈𝑘,𝑛=𝑘) given by (2.6). We have Π𝑛(0,𝑘)=𝛿0,𝑘 and Π𝑛(𝑛,𝑘)=𝛿𝑛,𝑘, and Π𝑛 is not irreducible. However, Π𝑛 is aperiodic and (apart from absorbing states) cannot be broken down into noncommunicating subsets; as a result it is diagonalizable, with eigenvalues |𝜆0||𝜆1||𝜆𝑛| and 1=𝜆0=𝜆1>|𝜆2|.

Example 2.3 (Dirichlet binomial). With 𝑈 a (0,1)-valued random variable with density beta(𝑘𝜃,(𝑛𝑘)𝜃)𝜈1,𝑛++𝜈𝑘,𝑛=𝑘=𝑛𝑘[𝑘𝜃]𝑘(𝑛𝑘)𝜃𝑛𝑘[𝑛𝜃]𝑛𝑛𝑘=𝔼𝑈𝑘(1𝑈)𝑛𝑘,(2.7) which is a beta mixture of the binomial distribution Bin(𝑛,𝑢).

Example 2.4. The Wright-Fisher model has a Bin(𝑛,𝑘/𝑛) transition matrix 𝑁𝑟+1(𝑚)=𝑘𝑁𝑟=𝑛𝑘(𝑚)=𝑘𝑘𝑛𝑘𝑘1𝑛𝑛𝑘.(2.8)

Remark 2.5 (statistical symmetry). Due to exchangeability of the reproduction law, neutral models are symmetric in the following sense. The transition probabilities of 𝑁𝑟(𝑚)=𝑛𝑁𝑟(𝑚) are equal to the transition probabilities of 𝑁𝑟(𝑚).

2.3. Backward in Time Process (neutral Coalescent)

Thecoalescentbackwardprocess can be defined as follows. Take a subsample of size 𝑚 from [𝑛] at generation 0. Identify two individuals from [𝑚] at each step if they share a common ancestor one generation backward in time. This defines an equivalence relation between 2 genes from the set [𝑚]. Define the induced ancestral backward process 𝒜𝑟(𝑚)𝑚={equivalenceclasses(partitions)of[𝑚]},𝑟,backwardintime.(2.9) The ancestral process is a discrete-time Markov chain with transition probability: 𝒜𝑟+1(𝑚)=𝛼𝒜𝑟(𝑚)=𝛽=𝑃𝛽;𝛼;with(𝛼,𝛽)𝑚,𝛼𝛽,(2.10) where, with 𝑎=|𝛼|= number of equivalence classes of 𝛼,𝑏=|𝛽|= number of equivalence classes of 𝛽,𝐛𝑎=(𝑏1,,𝑏𝑎) clusters sizes of 𝛽, and (𝑚)𝑎=𝑚(𝑚1)(𝑚𝑎+1) a falling factorial 𝑃𝛽;𝛼=𝑃(𝑛)𝑏;𝑎𝐛𝑎=(𝑛)𝑎(𝑛)𝑏𝔼𝑎𝑙=1𝜈𝑙,𝑛𝑏𝑙(2.11) is the probability of a 𝐛𝑎-merger. This is the probability that 𝑏 randomly chosen individuals out of 𝑛 have 𝑎𝑏 distinct parents, 𝑐 merging classes and cluster sizes 𝑏1𝑏𝑐2, 𝑏𝑐+1==𝑏𝑎=1.

If 𝑐=1 a unique multiple collision occurs of order 𝑏12.If 𝑏1=2 a simple binary collision occurs involving only two clusters.If 𝑐>1, simultaneous multiple collisions of orders 𝑏1𝑏𝑐2 occur. Thus, the jump's height of a transition 𝑏𝑎 is 𝑏𝑎=𝑐𝑖=1(𝑏𝑖1), corresponding to a partition of integer 𝑏𝑎 into 𝑐 summands, each 1.

The chain's state-space is {equivalence relations on (partitions of) {1,,𝑚}}; it has dimension 𝐵𝑚=𝑚𝑙=0𝑆𝑚,𝑙 (a Bell number), where 𝑆𝑚,𝑙 are the second-kind Stirling numbers.

The chain has initial state 𝒜0={(1),,(𝑚)}, and terminal absorbing state {(1,,𝑚)}.

From the Dirichlet Example 2.3, we get, 𝑃(𝑛)𝑏;𝑎(𝐛𝑎)=((𝑛)𝑎/[𝑛𝜃]𝑏)𝑎𝑚=1[𝜃]𝑏𝑚.
From the WF Example 2.4: In this case, 𝑃(𝑛)𝑏;𝑎(𝐛𝑎)=(𝑛)𝑎/𝑛𝑏 is the uniform distribution on {𝐛𝑎𝑏1++𝑏𝑎=𝑏}.

The Ancestral Count Process
Let 𝐴𝑟(𝑚)=#ancestorsatgeneration𝑟,backward-in-time.Then:𝐴𝑟(𝑚)=#blocksof𝒜𝑟(𝑚).(2.12) The backward ancestral count process is a discrete-time Markov chain with transition probabilities [17, 19]𝐴𝑟+1(𝑚)=𝑎𝐴𝑟(𝑚)=𝑏=𝑃(𝑛)𝑏,𝑎=𝑏!𝑎!𝑏1,,𝑏𝑎+𝑏1++𝑏𝑎=𝑏𝑃(𝑛)𝑏;𝑎𝐛𝑎𝑏1!𝑏𝑎!.=𝑛𝑎𝑛𝑏𝑏1,,𝑏𝑎+𝑏1++𝑏𝑎=𝑏𝔼𝑎𝑙=1𝜈𝑙,𝑛𝑏𝑙.(2.13) This Markov chain has state-space {0,,𝑚}, initial state 𝑚, absorbing states {0,1}. The process (𝐴𝑟;𝑟) has the transition matrix 𝑃𝑛 with entries 𝑃𝑛(𝑏,𝑎)=𝑃(𝑛)𝑏,𝑎 given by (2.13). Note, by inclusion-exclusion principle, the alternative alternating expression of 𝑃(𝑛)𝑏,𝑎: 𝑃(𝑛)𝑏,𝑎=𝑛𝑎𝑛𝑏𝑎𝑚=0(1)𝑎𝑚𝑎𝑚𝔼𝜈1,𝑛++𝜈𝑚,𝑛𝑏.(2.14)

2.4. Duality (neutral Case)

We start with a definition of the duality concept which is relevant to our purposes.

Definition 2.6 (see [11]). Two Markov processes (𝑋1𝑡,𝑋2𝑡;𝑡0), with state-spaces (1,2), are said to be dual with respect to some real-valued function Φ on the product space 1×2 if 𝑥11,𝑥22,𝑡0: 𝔼𝑥1Φ𝑋1𝑡,𝑥2=𝔼𝑥2Φ𝑥1,𝑋2𝑡.(2.15)

We then recall basic examples of dual processes from the neutral and exchangeable population models (Möhle[10]). The neutral forward and backward processes (𝑁𝑟,𝐴𝑟;𝑟) introduced in the two preceding subsections are dual with respect to the hypergeometric sampling without replacement kernels (i)Φ1𝑛(𝑚,𝑘)=𝑚𝑘𝑛𝑘,(ii)Φ2𝑛(𝑚,𝑘)=𝑘𝑛𝑚𝑛𝑘on{0,,𝑛}2.(2.16) Namely, (i) reads𝔼𝑚𝑁𝑟𝑘𝑛𝑘=𝔼𝑘𝑚𝐴𝑟𝑛𝐴𝑟=𝔼𝑘𝑛𝐴𝑟𝑛𝑚𝑛𝑛𝑚.(2.17) Call type 𝐴 individuals the descendants of the 𝑚 first chosen individuals (allele 𝐴) in the study of the forward process; type 𝑎 individuals are the remaining ones (allele 𝑎). The left-hand-side of the above equality is an expression of the probability that 𝑘-samples (without replacement) from population of size 𝑁𝑟 at time 𝑟 are all of type 𝐴, given 𝑁0=𝑚. If these 𝑘-samples are all descendants of 𝐴𝑟 ancestors at time 𝑟, this probability must be equal to the probability that a (𝑛𝑚)-sample from population of size 𝐴𝑟 at time 𝑟 are all of type 𝑎. This is the meaning of the right-hand-side.

Moreover, (ii) reads 𝔼𝑚𝑛𝑁𝑟𝑘𝑛𝑘=𝔼𝑘𝐴𝑛𝑚𝑟𝑛𝐴𝑟=𝔼𝑘𝑛𝐴𝑟𝑚𝑛𝑚.(2.18) The left-hand-side is the probability that a 𝑘-sample (without replacement) from population of size 𝑁𝑟 at time 𝑟 are all of type 𝑎, given 𝑁0=𝑚. If these 𝑘-samples are all descendants of 𝐴𝑟 ancestors at time 𝑟, this probability must be equal to the probability that 𝑚-samples from population of size 𝐴𝑟 at time 𝑟 are themselves all of type 𝑎.

With 𝑃𝑛 the transpose of 𝑃𝑛, a one-step (𝑟=1) version of these formulae is (i)Π𝑛Φ1𝑛=Φ1𝑛𝑃𝑛,(ii)Π𝑛Φ2𝑛=Φ2𝑛𝑃𝑛,(2.19) where (Φ1𝑛,Φ2𝑛) are 𝑛×𝑛 matrices with entries Φ1𝑛(𝑚,𝑘) and Φ2𝑛(𝑚,𝑘), respectively, and (Π𝑛,𝑃𝑛) the transition matrices of forward and backward processes. Note that the matrix Φ2𝑛 is symmetric and left-upper triangular. The matrices Φ1𝑛 and Φ2𝑛 are both invertible, with respective entries Φ1𝑛1(𝑖,𝑗)=(1)𝑖𝑗𝑖𝑗𝑛𝑖,Φ2𝑛1(𝑖,𝑗)=(1)𝑖+𝑗𝑛𝑖𝑛𝑖𝑛𝑗=(1)𝑖+𝑗𝑛𝑗𝑛𝑗.𝑛𝑖(2.20) The matrix [Φ1𝑛]1 is left-lower triangular, while [Φ2𝑛]1 is symmetric right-lower triangular. Thus, Φ(i)1𝑛1Π𝑛Φ1𝑛=𝑃𝑛,Φ(ii)2𝑛1Π𝑛Φ2𝑛=𝑃𝑛.(2.21) In any case, being similar matrices, Π𝑛 and 𝑃𝑛 (or 𝑃𝑛) both share the same eigenvalues. If 𝑅𝑛 diagonalizing Π𝑛 is known so that 𝑅𝑛1Π𝑛𝑅𝑛=Λ𝑛=diag(𝜆0,,𝜆𝑛), the diagonal matrix of the eigenvalues of Π𝑛, then, with Φ𝑛=Φ1𝑛 or Φ2𝑛, 𝑅𝑛=Φ𝑛1𝑅𝑛 diagonalizes 𝑃𝑛 and is obtained for free (and conversely). 𝑅𝑛 is the matrix whose columns are the right-eigenvectors of Π𝑛 and 𝑅𝑛 is the matrix whose columns (rows) are the right eigenvectors (left-eigenvectors) of 𝑃𝑛 (of 𝑃𝑛). Similarly, if 𝐿𝑛 is the matrix whose rows are the left eigenvectors of Π𝑛, 𝐿𝑛=𝐿𝑛Φ𝑛 is the matrix whose rows (columns) are the left eigenvectors (right eigenvectors) of 𝑃𝑛 (of 𝑃𝑛). With 𝑙𝑘 the 𝑘th row of 𝐿𝑛 and 𝑟𝑘 the 𝑘th column of 𝑅𝑛, the spectral decomposition of Π𝑛 is Π𝑟𝑛=𝑛𝑘=0𝜆𝑟𝑘𝑟𝑘𝑙𝑘𝑙𝑘𝑟𝑘,𝑟,(2.22) whereas, with 𝑙𝑘 the 𝑘th column of 𝐿𝑛 and 𝑟𝑘 the 𝑘th row of 𝑅𝑛, the one of 𝑃𝑛 reads 𝑃𝑟𝑛=𝑛𝑘=0𝜆𝑟𝑘𝑙𝑘𝑟𝑘𝑟𝑘𝑙𝑘=𝑛𝑘=0𝜆𝑟𝑘Φ𝑛𝑙𝑘Φ𝑛1𝑟𝑘Φ𝑛1𝑟𝑘Φ𝑛𝑙𝑘,𝑟.(2.23) In Möhle [10], a direct combinatorial proof of the duality result can be found (in the general exchangeable or neutral case); it was obtained by directly checking the consistency of (2.6), (2.13), and (2.16).

The duality formulae allow one to deduce the probabilistic structure of one process from the one of the other. The question we address now is does duality still make sense in nonneutral situations? We will see that it does in discrete nonneutral Wright-Fisher models, but only for some class of state-dependent transition frequencies.

3. beyond Neutrality (symmetry Breaking)

Discrete forward nonneutral models (with nonnull drifts) can be obtained by substituting 𝑘𝑘𝑛𝑝𝑛𝜈in1,𝑛++𝜈𝑘,𝑛=𝑘,(3.1) where 𝑝(𝑥)𝑥(0,1)(0,1)iscontinuous,increasing,with𝑝(0)=0,𝑝(1)=1.(3.2)𝑝(𝑥) is the state-dependent Bernoulli bias probability different from identity 𝑥 (as in neutral case).

When particularized to the WF model, this leads to the biased transition probabilities 𝑁𝑟+1(𝑚)=𝑘𝑁𝑟=𝑛𝑘(𝑚)=𝑘𝑝𝑘𝑛𝑘𝑘1𝑝𝑛𝑛𝑘.(3.3) In this binomial 𝑛-sampling with replacement model, a type 𝐴 individual is drawn with probability 𝑝(𝑘/𝑛) which is different from the uniform distribution 𝑘/𝑛 due to bias effects.

From this, we conclude (a symmetry breaking property) the following: the transition probabilities of 𝑁𝑟(𝑚)=𝑛𝑁𝑟(𝑚), 𝑟 are 𝑘Bin𝑛,1𝑝1𝑛𝑘Bin𝑛,𝑝𝑛,(3.4) and so, no longer coincide with the ones of (𝑁𝑟(𝑚);𝑟). The process 𝑁𝑟(𝑚), 𝑟 no longer is a martingale. Rather, if 𝑥𝑝(𝑥) is concave (convex), 𝑁𝑟(𝑚), 𝑟 is a submartingale (supermartingale) because 𝔼(𝑁𝑟+1(𝑚)𝑁𝑟(𝑚))=𝑛𝑝(𝑁𝑟(𝑚)/𝑛)𝑁𝑟(𝑚) (resp., 𝑁𝑟(𝑚)).

In the binomial neutral Wright-Fisher transition probabilities, we replaced the success probability 𝑘/𝑛 by a more general function 𝑝(𝑘/𝑛). However, this replacement leaves open the question what model is in the background and what quantity the process (𝑁𝑟,𝑟) really counts? A concrete model in terms of offspring variables must be provided instead. To address this question, we emphasize that the reproduction law corresponding to the biased binomial model is multinomial and asymmetric, namely, 𝝂𝑛𝑑 Multin(𝑛;𝝅𝑛), where 𝝅𝑛=(𝜋1,𝑛,,𝜋𝑛,𝑛) and 𝜋𝑚,𝑛=𝑝(𝑚/𝑛)𝑝((𝑚1)/𝑛),𝑚=1,,𝑛. We note that under our hypothesis, 𝑛𝑚=1𝜋𝑚,𝑛=𝑝(1)𝑝(0)=1.(3.5) Due to its asymmetry, the law of the biased 𝝂𝑛 no longer is exchangeable.

We now recall some well-known bias examples arising in population genetics.

Example 3.1 (homographic model, selection). Assume 𝑝(𝑥)=(1+𝑠)𝑥1+𝑠𝑥,(3.6) where 𝑠>1 is a selection parameter. This model arises when gene 𝐴 (resp., 𝑎), with frequency 𝑥 (resp., 1𝑥), has fitness 1+𝑠 (resp., 1). The case 𝑠>0 arises when gene of type 𝐴 is selectively advantageous, whereas it is disadvantageous when 𝑠(1,0).

Example 3.2 (selection with dominance). Assume 𝑝(𝑥)=(1+𝑠)𝑥2+(1+𝑠)𝑥(1𝑥)1+𝑠𝑥2+2𝑠𝑥(1𝑥).(3.7) In this model, genotype 𝐴𝐴 (resp., 𝐴𝑎 and 𝑎𝑎), with frequency 𝑥2 (resp., 2𝑥(1𝑥) and (1𝑥)2) has fitness 1+𝑠 (resp., 1+𝑠 and 1). is a measure of the degree of dominance of heterozygote 𝐴𝑎. We impose 𝑠>1 and 𝑠>1. Note that the latter quantity can be put into the canonical form of deviation to neutrality: 𝑝(𝑥)=𝑥+𝑠𝑥(1𝑥)𝑥(21)1+𝑠𝑥2+2𝑠𝑥(1𝑥),(3.8) where the ratio appearing in the right-hand-side is the ratio of the difference of marginal fitnesses of 𝐴 and 𝑎 to their mean fitness. The case =1/2 corresponds to balancing selection with 𝑝(𝑥)=𝑥+(𝑠/2)(𝑥(1𝑥)/(1+𝑠𝑥)).

Example 3.3 (quadratic model). With 𝑐[1,1], a curvature parameter, one may choose 𝑝(𝑥)=𝑥(1+𝑐𝑐𝑥).(3.9) If 𝑐=1, 𝑝(𝑥)=𝑥(2𝑥)=1(1𝑥)2, this bias appears in a discrete 2-sex population model [9, 10]. We will give below an interpretation of this quadratic model when 𝑐(0,1] in terms of a joint one-way mutations and neutrality effects model.
We can relax the assumption 𝑝(0)=0,𝑝(1)=1 by assuming 0𝑝(0)𝑝(1)1, 𝑝(1)𝑝(0)[0,1).

Example 3.4 (affine model). Take, for example, 𝑝(𝑥)=(1𝜇2)𝑥+𝜇1(1𝑥),(3.10) where (𝜇1,𝜇2) are mutation probabilities, satisfying 𝜇11𝜇2. It corresponds to the mutation scheme 𝑎𝜇1𝜇2𝐴. To avoid discussions of intermediate cases, we will assume that 𝑝(0)=𝜇1>0 and 𝑝(1)<1(𝜇2>0). In this case, the matrix Π𝑛 is irreducible and even primitive and all states of this Markov chain are now recurrent. We have (𝑁𝑟+1>0𝑁𝑟=0)=1(1𝑝(0))𝑛>0 and (𝑁𝑟+1<𝑛𝑁𝑟=𝑛)=1𝑝(1)𝑛>0 and the boundaries {0} and {𝑛} no longer are strictly absorbing as there is a positive reflection probability inside the domain {0,1,,𝑛}.
For reasons to appear now, we will be only interested in functions 𝑞 such that 𝑞(𝑥)=1𝑝(𝑥) is a completely monotone function (CM) on (0,1), that is, satisfying (1)𝑙𝑞(𝑙)(𝑥)0,𝑥(0,1),(3.11) for all order-𝑙 derivatives 𝑞(𝑙) of 𝑞,𝑙0. If 𝑝(𝑥) is such that 𝑞 is CM, we will call it an admissible bias mechanism.

4. Nonneutral Wright-Fisher Models and Duality

Preliminaries. Let 𝐯𝑛=(𝑣(0),𝑣(1),,𝑣(𝑛)) be a (𝑛+1)-vector of [0,1]-valued numbers. Define the backward difference operator acting on 𝐯𝑛 by: 𝑣(𝑚)=𝑣(𝑚)𝑣(𝑚1),𝑚=1,,𝑛. We have 0𝑣(𝑚)=𝑣(𝑚),2𝑣(𝑚)=𝑣(𝑚)2𝑣(𝑚1)+𝑣(𝑚2), and so forth, and, starting from the endpoint 𝑣(𝑛)𝑗||𝑣(𝑚)𝑚=𝑛=𝑗𝑙=0(1)𝑗𝑙𝑗𝑙𝑣(𝑛𝑙),𝑗=0,,𝑛.(4.1) Let 𝑢 be some continuous function [0,1][0,1]. Consider the (𝑛+1)-vector 𝐮𝑛=(𝑢(0/𝑛),,𝑢(𝑚/𝑛),,𝑢(𝑛/𝑛)), sampling 𝑢 at points 𝑚/𝑛. The function 𝑢 is said to be -completely monotonic if (1)𝑗𝑗𝑢(𝑚/𝑛)|𝑚=𝑛0, for all 𝑗=0,,𝑛 and all 𝑛0. Let (𝑢1,𝑢2) be two continuous functions on [0,1]. Let 𝑢=𝑢1𝑢2. With 𝐮𝑛 the pointwise product of 𝐮1𝑛 and 𝐮2𝑛, assuming both functions (𝑢1,𝑢2) to be -completely monotonic, so will be 𝑢, by the Leibniz rule. In particular, if 𝑢 is -completely monotonic, so will be its integral powers 𝑢𝑖,𝑖. Our main result is as follows.

Theorem 4.2. Consider a nonneutral WF forward model (𝑁𝑟;𝑟) on {0,,𝑛}, with continuous, nondecreasing bias 𝑝(𝑥), satisfying 0𝑝(0)𝑝(1)1,𝑝(1)𝑝(0)[0,1].(4.2) This process has forward transition matrix Π𝑛𝑘,𝑘𝜈=1,𝑛++𝜈𝑘,𝑛=𝑘=𝑛𝑘𝑝𝑘𝑛𝑘𝑘1𝑝𝑛𝑛𝑘.(4.3) There exists a Markov chain (𝐴𝑟;𝑟) on {0,,𝑛} such that (𝑁𝑟,𝐴𝑟;𝑟) are dual with respect to Φ2𝑛(𝑚,𝑘)=(𝑘𝑛𝑚)/(𝑛𝑘) if and only if 𝑥𝑞(𝑥)=1𝑝(𝑥) is completely monotone on (0,1). In this case, the transition probability matrix of (𝐴𝑟;𝑟) is 𝑃𝑛𝑛𝑗(𝑖,𝑗)=𝑗𝑙=0(1)𝑗𝑙𝑗𝑙𝑞𝑙1𝑛𝑖0,(4.4)where𝑃𝑛 is a stochastic matrix if and only if 𝑝(0)=0; else, if 𝑝(0)>0, the matrix 𝑃𝑛 is substochastic.

Proof. Developing [Φ2𝑛]1Π𝑛Φ2𝑛=𝑃𝑛, we easily obtain 𝑃𝑛(𝑗,𝑖)=𝑃𝑛=𝑛𝑗(𝑖,𝑗)𝑗𝑙=0(1)𝑗𝑙𝑗𝑙1𝑝𝑛𝑙𝑛𝑖=𝑛𝑗(1)𝑗𝑗𝑞𝑚𝑛𝑖||||𝑚=𝑛.(4.5) This entry is nonnegative if and only if (1)𝑗𝑗(𝑞(𝑚/𝑛)𝑖)|𝑚=𝑛0, for all 𝑖,𝑗=0,,𝑛. But due to the above argument on -complete monotonicity of integral powers, this will be the case if and only if (1)𝑗𝑗(𝑞(𝑚/𝑛))|𝑚=𝑛0, for all 𝑗=0,,𝑛. As this must be true for arbitrary value of population size 𝑛, function 𝑞 has to be -completely monotonic. Adapting now the arguments of Theorem 2 developed in [20, page 223], for absolutely monotone functions on (0,1), this will be the case if and only if 𝑥𝑞(𝑥)=1𝑝(𝑥) is a completely monotone function on (0,1) in the sense that (1)𝑙𝑞(𝑙)(𝑥)0,𝑥(0,1),𝑙.(4.6) Next, since (𝐼)𝑢(𝑚)=𝑢(𝑚1) is a simple back-shift, 𝑛𝑗=0𝑃𝑛(𝑗,𝑖)=𝑛𝑗=0𝑃𝑛(𝑖,𝑗)=(𝐼)𝑛𝑞𝑚𝑛𝑖||||𝑚=𝑛=𝑞(0)𝑖,(4.7) and if 𝑞 is CM, 𝑃𝑛 is a stochastic matrix if and only if 𝑞(0)=1; else, if 𝑞(0)<1, the matrix 𝑃𝑛 is substochastic.
We note that the first column of the matrix 𝑃𝑛 is 𝑃𝑛(𝑖,0)=𝑞(1)𝑖 whereas its first line is 𝑃𝑛(0,𝑗)=𝛿0,𝑗, expressing, as required, that the state 0 of (𝐴𝑟;𝑟) is absorbing.

5. Examples

We show here that most of the simplest nonneutrality mechanisms used in the literature fall within the class which we would like to draw the attention on, or are amenable to it via some “reciprocal transformation” which we define. Elementary algebraic manipulations on CM functions allow to exhibit a vast class of unsuspected mechanisms. Note that in some cases, their biological relevance remains to be elucidated. The results presented in this section seem to be new. They serve as an illustration of our theorem.

5.1. Elementary Mechanisms

Assume first that 𝑝(𝑥)=𝑥 corresponding to the simple neutral case. Then 𝑞(𝑥)=1𝑥 is completely monotone on (0,1). With 𝑆𝑖,𝑗 the second kind Stirling numbers, we get a lower left triangular stochastic transition matrix 𝑃𝑛𝑛𝑗(𝑖,𝑗)=𝑗𝑙=0(1)𝑗𝑙𝑗𝑙𝑙𝑛𝑖=(𝑛)𝑗𝑛𝑖𝑆𝑖,𝑗,𝑗𝑖0,else.(5.1) The diagonal terms (eigenvalues) are all distinct with 𝑃𝑛(𝑖,𝑖)=(𝑛)𝑖𝑛𝑖. The matrix 𝑃𝑛 is stochastic. Due to triangularity, ancestral process is a pure death Markov process which may be viewed as a discrete coalescence tree.

From Example 3.4 (mutation). Assume that (3.10) holds 𝑝(𝑥)=(1𝜇2)𝑥+𝜇1(1𝑥) where (𝜇1,𝜇2) are mutation probabilities. Then, with 𝜅=1(𝜇1+𝜇2), 𝑞(𝑥)=1𝜇1𝜅𝑥 is completely monotone on (0,1) if and only 𝜇11𝜇2(𝜅0). In this case, 𝑃𝑛 is again lower left triangular (a pure death process). We have 𝑃𝑛𝑛𝑗(𝑖,𝑗)=𝑗𝑙=0(1)𝑗𝑙𝑗𝑙𝜇2𝑙+𝜅𝑛𝑖=(𝑛)𝑗𝑛𝑖𝑆𝜇2𝑖,𝑗𝜅𝑛,𝑗𝑖0,else.(5.2) in terms of generalized Stirling numbers 𝑆𝜇2𝑖,𝑗(𝜅/𝑛). We have 𝑃𝑛(𝑖,𝑖)=(𝑛)𝑖(𝜅/𝑛)𝑖 and the spectrum of 𝑃𝑛 is real. When 𝜇1>0, this matrix is substochastic with 𝑛𝑗=0𝑃𝑛(𝑖,𝑗)=(1𝜇1)𝑖.

A particular case deals with one-way mutations (𝜇1+𝜇2>0,𝜇1𝜇2=0).

If 𝜇2=0,𝑃𝑛(𝑖,𝑗)=(1𝜇1)𝑖(𝑛)𝑗𝑛𝑖𝑆𝑖,𝑗,𝑗𝑖,=0, else. Further, 𝑛𝑗=0𝑃𝑛(𝑖,𝑗)=(1𝜇1)𝑖<1.If 𝜇1=0,𝑃𝑛(𝑖,𝑗)=(𝑛)𝑗𝑛𝑖𝑆𝜇2𝑖,𝑗((1𝜇2)/𝑛),𝑗𝑖,=0, else. The corresponding matrix 𝑃𝑛 is stochastic.

From Example 3.3 (quadratic). Assume that𝑝(𝑥)=𝑥(1+𝑐𝑐𝑥), as in (3.9). Then 𝑞(𝑥)=(1𝑥)(1𝑐𝑥) which is completely monotone if and only if 𝑐[0,1]. The case 𝑐=0 is the neutral case, whereas 𝑐=1 appears in a 2-sex model of Möhle. In this quadratic case, since 𝑗(𝑞(𝑚/𝑛)𝑖)=0 if 𝑗>2𝑖, then 𝑃𝑛(𝑖,𝑗)=0 if 𝑗>2𝑖 and so 𝑃𝑛 is a Hessenberg-like matrix. Note that 𝑛𝑗=0𝑃𝑛(𝑖,𝑗)=𝑞(0)𝑖=1.

From the selection Example 3.1, when (3.6) holds 𝑝(𝑥)=(1+𝑠)𝑥1+𝑠𝑥,(5.3)where𝑞(𝑥)=1𝑝(𝑥)=(1𝑥)/(1+𝑠𝑥) is CM whenever selection parameter 𝑠>0. The induced matrix 𝑃𝑛 is stochastic. It is no longer lower left triangular so that the ancestral no longer is a pure death process, being rather a birth and death process. The induced coalescence pattern no longer is a discrete tree, but rather a graph (a discrete version of the ancestral selection graph of Neuhauser-Krone [12]).

From Example 3.2 (selection with dominance). The corresponding mechanism (3.7) with parameters (𝑠,) satisfying 𝑠>1 and 𝑠>1 is CM if and only if 𝑠𝑠=(12)/2>0 and (0,1/2). The case (0,1) corresponds to directional selection, where genotype 𝐴𝐴 has highest fitness compared to aa's and the heterozygote class 𝐴𝑎 has intermediate fitness compared to both homozygote classes. In this situation, marginal fitness of 𝐴 exceeds the one of 𝑎 and selective sweep is expected. When (0,1/2), allele 𝐴 is dominant to 𝑎, whereas when (1/2,1), allele 𝐴 is recessive to 𝑎 (a stabilizing effect slowing down the sweep). Critical value =1/2 is a case of pure genic balancing selection.

Example 5.1. Consider the mechanism 𝑝(𝑥)=𝑥𝛾 for some 𝛾>0. The function 𝑞(𝑥)=1𝑝(𝑥) is CM if and only if 𝛾(0,1). Although this model seems quite appealing, we could find no reference to it in the specialized mathematical genetics literature.

Example 5.2 (reciprocal mechanism). If 𝑝(𝑥) is not admissible in that 𝑞 is not CM, it can be that 𝑝(𝑥)=1𝑝(1𝑥) is itself admissible. As observed before, if 𝑁𝑟(𝑚) has transition probabilities given by Bin(𝑛,𝑝(𝑘/𝑛)),𝑝(𝑥) arises in the transition probabilities of 𝑁𝑟(𝑚)=𝑛𝑁𝑟(𝑚). Indeed, such transitions are Bin(𝑛,𝑝(𝑘/𝑛)) distributed.

If 𝑝(𝑥) is the selection mechanism of Example 3.1, (3.6), with 𝑠(1,0) (not admissible), 𝑝(𝑥)=(1+𝑠)𝑥/(1+𝑠𝑥) is itself an admissible selection mechanism because it has reciprocal selection parameter 𝑠=𝑠/(1+𝑠)>0. If 𝑝(𝑥) is the mechanism of Example 3.2, namely, (3.7), with parameters (𝑠,), then 𝑝(𝑥) is itself a selection with dominance mechanism with reciprocal parameters 𝑠=𝑠/(1+𝑠) and =1. Assuming that(𝑠𝑠<0,(1/2,1)), 𝑝(𝑥) is not admissible whereas 𝑝(𝑥) is because 𝑠𝑠>0 and (0,1/2). Similarly, when 𝛾(0,1), the mechanism 𝑝(𝑥)=1(1𝑥)𝛾 is not admissible but, from Example 5.1, 𝑝(𝑥)=1𝑝(1𝑥)=𝑥𝛾 is as follows.

5.2. Bias Mechanisms with Mutational Effects

Let 𝑝𝑀(𝑥)=(1𝜇2)𝑥+𝜇1(1𝑥) be the mutational bias mechanism (with 𝜅=1(𝜇1+𝜇2)0). Let 𝑝(𝑥) be a bias mechanism such that 𝑞(𝑥) is CM with 𝑝(1)𝑝(0)=1. Then, 𝑝𝑀(𝑥)=𝑝𝑀𝑝(𝑥)(5.4) is such that 𝑞𝑀(𝑥)=1𝑝𝑀(𝑥) is CM. It is, therefore, admissible and adds mutational effects to the primary mechanism 𝑝(𝑥). For example, 𝑝𝑀𝜇(𝑥)=1+𝑥(1+𝑠)1𝜇2𝜇11+𝑠𝑥(5.5) is a mechanism of selection combined with mutational effects. We have 𝑝𝑀(0)=𝜇1, 𝑝𝑀(1)=1𝜇2. The mechanisms 𝑝𝑀(𝑥) obtained in this way all share the specificity 𝑝𝑀(1)𝑝𝑀(0)=𝜅<1.

Note that, except for the mutational affine mechanism, it is not true in general that whenever 𝑝1(𝑥) and 𝑝2(𝑥) are two admissible bias mechanisms, then 𝑝1(𝑝2(𝑥)) is admissible.

5.3. Joint Bias Effects and Compound Bias

Let 𝑝1(𝑥) and 𝑝2(𝑥) be two admissible bias in that 𝑞1(𝑥)=1𝑝1(𝑥) and 𝑞2(𝑥)=1𝑝2(𝑥) are both completely monotone. Then 𝑞(𝑥)=𝑞1(𝑥)𝑞2(𝑥)isCM.(5.6) Thus, with 𝑥1𝑥2=𝑥1+𝑥2𝑥1𝑥2, the probabilistic product in [0,1]𝑝1(𝑥),𝑝2(𝑥)𝑝(𝑥)=𝑝1(𝑥)𝑝2(𝑥).(5.7) Whenever a WF model is considered with bias 𝑝(𝑥)=𝑝1(𝑥)𝑝2(𝑥) obtained from two distinct bias 𝑝1(𝑥) and 𝑝2(𝑥), we call it a WF model with joint bias effect.

Example 5.3 (joint selection and mutational effects). Let 𝑝1(𝑥)=𝑝𝑀(𝑥) and 𝑝2(𝑥)=(1+𝑠)𝑥/(1+𝑠𝑥). We get 𝑞(𝑥)=1𝜇1𝜅𝑥(1𝑥),𝜇1+𝑠𝑥𝑝(𝑥)=1+𝑥𝑠+1𝜇1+𝜅𝜅𝑥2,1+𝑠𝑥(5.8) with 𝑝(0)=𝜇1, 𝑝(1)=1. This mechanism differs from the traditional mechanism of selection combined with mutational effects.

Example 5.4 (joint mutation and neutral effects). Let 𝑝1(𝑥)=(1𝜇2)𝑥+𝜇1(1𝑥) and 𝑝2(𝑥)=𝑥. We get 𝑞(𝑥)=(1𝑥)(1𝜇1𝜅𝑥),𝑝(𝑥)=𝜇1+𝑥1𝜇1,+𝜅(1𝑥)(5.9) with 𝑝(0)=𝜇1, 𝑝(1)=1. When 𝜇1=0 (one-way mutations), we recover the quadratic mechanism with curvature parameter 𝑐=1𝜇2. This finding justifies some interest into the quadratic mechanisms with 𝑐1.
With 𝑗=1,2, the reproduction law of each elementary effect is 𝝂𝑗𝑛𝑑 Multin(𝑛;𝝅𝑗𝑛), where 𝜋𝑗𝑚,𝑛=𝑝𝑗(𝑚/𝑛)𝑝𝑗((𝑚1)/𝑛),𝑚=1,,𝑛. Then, 𝝂𝑛𝑑 Multin(𝑛;𝝅𝑛),𝜋𝑚,𝑛=𝑝(𝑚/𝑛)𝑝(𝑚1)/𝑛),𝑚=1,,𝑛, where 𝝅𝑛=𝝅1𝑛𝝅2𝑛 is easily obtained componentwise by 𝜋𝑚,𝑛=𝜋1𝑚𝑚,𝑛𝑙=1𝜋2𝑙,𝑛+𝜋2𝑚𝑚,𝑛𝑙=1𝜋1𝑙,𝑛,𝑚=1,,𝑛.(5.10) We let 𝝂𝑛=𝝂1𝑛𝝂2𝑛𝑑 Multin(𝑛;𝝅1𝑛𝝅2𝑛). It is the reproduction law of a WF model obtained jointly from the two bias 𝑝1(𝑥) and 𝑝2(𝑥).

Let 𝜙(𝑥)(0,1)(0,1) be an absolutely monotone function satisfying 𝜙(𝑙)(𝑥)0 for all 𝑙th derivatives 𝜙(𝑙) of 𝜙, for all 𝑥(0,1). Such functions are well known to be probability generating functions (pgfs) of -valued random variables, say 𝑁, that is to say, 𝜙(𝑥)=𝔼[𝑥𝑁]. Clearly, if 𝑞 is CM on (0,1), then so is 𝑞𝜙(𝑥)=𝜙(𝑞(𝑥)). Thus, 𝑝𝜙(𝑥)=1𝜙(1𝑝(𝑥)) is an admissible bias mechanism in that 𝑞𝜙(𝑥)=1𝑝𝜙(𝑥) is CM. We call it a compound bias.

Example 5.5. The general mechanism with mutational effects is in this class. Indeed, 𝑞𝑀(𝑥)=1𝑝𝑀(𝑥)=1𝑝𝑀(𝑝(𝑥))=1𝑝𝑀(1𝑞(𝑥)),(5.11) and so 𝜙(𝑥)=1𝑝𝑀(1𝑥)=1(1𝜇2)(1𝑥)𝜇1𝑥=𝜇2+𝜅𝑥 which is absolutely monotone as soon as 𝜅=1(𝜇1+𝜇2)0.

Example 5.6. With 𝜃>0, taking 𝜙(𝑥)=𝑒𝜃(1𝑥) or (𝑒𝜃𝑥1)/(𝑒𝜃1), the pgf of a Poisson (or shifted-Poisson) random variable, 𝑝𝜙(𝑥)=1𝜙(1𝑝(𝑥)) is admissible if 𝑝(𝑥) is. Note that if 𝑞 is of the form 𝑞𝜙, where 𝜙 is the pgf of the Poisson random variable, then 𝑞𝜙(𝑥)𝛼 is admissible, for all 𝛼>0, a property reminiscent of infinite divisibility for pgfs. Taking 𝜙(𝑥)=(1𝜋)/(1𝜋𝑥) or 𝑥(1𝜋)/(1𝜋𝑥), 𝜋(0,1), the pgf of a geometric (or shifted-geometric) random variable, 𝑝𝜙(𝑥)=𝑠𝑝(𝑥)/(1+𝑠𝑝(𝑥)) or (𝑠+1)𝑝(𝑥)/(1+𝑠𝑝(𝑥)) is admissible if 𝑝(𝑥) is (with 𝑠=𝜋/(1𝜋)>0). In the external latter mechanism, one recognizes the one in (3.6) occurring in the model with selection of Example 3.1.

Example 5.7. Let 𝑝(𝑥)=𝑥𝛾 with 𝛾(0,1) as in Example 5.1. Then, 𝑝𝜙(𝑥)=1𝑞𝜙(𝑥), where 𝑞𝜙(𝑥)=𝑒𝜃(1𝑞(𝑥))=𝑒𝜃𝑥𝛾,𝜃>0, is admissible. Note that 𝑝𝜙(𝑥)𝑥0𝜃𝑥𝛾. The reciprocal function 𝑝𝜙(𝑥)=𝑞𝜙(1𝑥)=𝑒𝜃(1𝑥)𝛾 also interprets as an absolutely monotone discrete-stable pgf (see Steutel and van Harn [21]). It is not admissible.

Proceeding in this way, one can produce a wealth of admissible bias probabilities 𝑝𝜙, the signification of which in population genetics remaining though to be pinpointed, in each specific case study.

6. Limit Laws

Consider a WF model (𝑁𝑟;𝑟) on {0,,𝑛} with forward transition matrix Π𝑛𝑘,𝑘=𝑛𝑘𝑝𝑘𝑛𝑘𝑘1𝑝𝑛𝑛𝑘,(6.1) with admissible bias 𝑝(𝑥). Define (𝐴𝑟;𝑟) as the dual Markov chain on {0,,𝑛} with transition probability 𝑃𝑛𝑛𝑗(𝑖,𝑗)=𝑗𝑙=0(1)𝑗𝑙𝑗𝑙𝑞𝑙1𝑛𝑖.(6.2) Then, (𝑁𝑟,𝐴𝑟;𝑟) are dual with respect to Φ𝑛(𝑚,𝑘)=Φ2𝑛(𝑚,𝑘)=(𝑘𝑛𝑚)/(𝑛𝑘), to wit 𝔼𝑚𝑛𝑁𝑟𝑘𝑛𝑘=𝔼𝑘𝐴𝑛𝑚𝑟𝑛𝐴𝑟=𝔼𝑘𝑛𝐴𝑟𝑚𝑛𝑚.(6.3) We will distinguish two cases.

Case 1. Assume first that 𝑁𝑟𝑑𝑁as𝑟,independentlyof𝑁0=𝑚1.(6.4) Let 𝜋(𝑖)=(𝑁=𝑖) and 𝝅=(𝜋(0),,𝜋(𝑛)). The line vector 𝝅 is the left eigenvector of Π𝑛 associated to the eigenvalue 1𝝅=𝝅Π𝑛. It is the (unique) invariant probability measure (stationary distribution) of (𝑁𝑟;𝑟).
If this stationary distribution exists, then, using duality formula, necessarily, 𝐴𝑟0 as 𝑟 with probability 𝑘(𝐴=0)=𝜌(𝑘)<1. The numbers 𝜌(𝑘) are the extinction probabilities of the dual process started at 𝑘. As is well known, 𝜌=(𝜌(0),,𝜌(𝑛)) is the unique solution to (𝐼𝑃𝑛)𝝆=0 with 𝜌(0)=1.

Remark 6.1. Typical situations where (𝑁𝑟;𝑟) has an invariant measure is when mutational effects are present, and more generally when the bias mechanism satisfies 𝑝(0)>0 and 𝑝(1)<1. In this situation, the forward stochastic transition matrix Π𝑛 has an algebraically simple dominant eigenvalue 1. By Perron-Frobenius theorem lim𝑟Π𝑟𝑛=𝟏𝝅,(6.5) where 𝟏=(1,,1). The invariant probability measure can be approximated by subsequent iterates of Π𝑛, the convergence being exponentially fast, with rate governed by the second largest eigenvalue. Of course, detailed balance (stating that 𝜋𝑘Π𝑛(𝑘,𝑘)=𝜋𝑘Π𝑛(𝑘,𝑘)) does not hold here and the forward chain in equilibrium is not time-reversible.
In these recurrent cases, the dual ancestral process 𝐴𝑟 started at 𝑘 gets extinct with probability 𝜌(𝑘). The numbers 1𝜌(𝑘) are the probabilities that it gets killed before getting extinct; in other words, 1𝜌(𝑘) are the probabilities that 𝐴𝑟 first hits an extra coffin state, say {𝜕}, before hitting {0}.

In terms of moments, by the duality formula, we conclude that 𝑛𝑘1𝔼𝑛𝑁𝑘=𝜌(𝑘)=𝑘(𝐴=0),(6.6)relating 𝑘-factorial moments of 𝑛𝑁 to the extinction probabilities of 𝐴𝑟 given 𝐴0=𝑘. We also have 𝑛𝑘=0𝑣𝑘𝔼𝑛𝑁𝑘=𝔼(1+𝑣)𝑛𝑁=𝑛𝑘=0𝑛𝑘𝜌(𝑘)𝑣𝑘,(6.7) and so the probability generating function of 𝑁 can be expressed as (𝑢[0,1]): 𝔼𝑢𝑁=𝑛𝑘=0𝑛𝑘𝜌(𝑘)𝑢𝑛𝑘(1𝑢)𝑘,(6.8) in terms of the Bernstein-Bézier polynomial of (𝜌(𝑛𝑘);𝑘=0,,𝑛).

Let 𝝆=(𝜌(0),,𝜌(𝑛)). The vector 𝝆 is the right eigenvector of 𝑃𝑛 associated to the eigenvalue 1𝝆=𝑃𝑛𝝆. In this case, the matrix 𝑃𝑛 is sub-stochastic and the extinction probability of (𝐴𝑟;𝑟) given 𝐴0=𝑘 is less than one. Thanks to duality, we have Π𝑛Φ𝑛=Φ𝑛𝑃𝑛,(6.9) where the matrix Φ𝑛 is symmetric whereas the matrix Φ𝑛1 is symmetric right lower triangular, with Φ𝑛(𝑚,𝑘)=𝑘𝑛𝑚𝑛𝑘=𝑚𝑛𝑘𝑛𝑚,Φ𝑛1(𝑖,𝑗)=(1)𝑖+𝑗𝑛𝑖𝑛𝑖𝑛𝑗=(1)𝑖+𝑗𝑛𝑗𝑛𝑗.𝑛𝑖(6.10) Thus, 𝝅Π𝑛Φ𝑛=𝝅Φ𝑛=𝝅Φ𝑛𝑃𝑛,(6.11) showing that 𝝆 and 𝝅 are related through 𝝆=Φ𝑛𝝅or𝝅=Φ𝑛1𝝆.(6.12)The knowledge of the invariant measure 𝝅of the forward process allows one to compute the extinction probabilities 𝝆of the dual backward ancestral process (and conversely).

Consider the discrete WF model with mutations of Example 3.4. In this case, 𝑁𝑟𝑑𝑁 as 𝑟, regardless of 𝑁0=𝑚 and (𝑁𝑟;𝑟) has an invariant measure which is difficult to compute. Looking at the backward process, the matrix 𝑃𝑛 is substochastic (if 𝜇1>0) and lower left triangular. Due to triangularity, the right eigenvector 𝝆 of 𝑃𝑛 can easily be computed explicitly in terms of (𝑃𝑛(𝑖,𝑗);𝑗𝑖),𝑖=0,,𝑛. We, therefore, get the following alternating expression for the invariant measure: 𝜋𝑛𝑖(𝑖)=𝑖𝑗=0(1)𝑖𝑗𝑖𝑗𝜌(𝑛𝑗).(6.13) Concerning moments, for instance, we have 𝜌(1)=𝜇2/(𝜇1+𝜇2) so that 𝔼[𝑁]=𝑛𝜇1/(𝜇1+𝜇2); from (5.2) we also have 𝜌𝜇(2)=2𝑛𝜇2(1+𝜅)+𝜅2(1𝜅)𝑛(𝑛1)𝜅2=1𝑛(𝑛1)𝑛(𝑛1)𝑛(2𝑛1)𝜇1𝜇1+𝜇2𝑁+𝔼2,(6.14) allowing to compute 𝔼[𝑁2] and then the variance of 𝑁. We get 𝜎2(𝑁𝑛)=2𝜇1𝜇2𝜇1+𝜇22𝜇2𝑛1+𝜇2+1+𝑜(𝑛)𝑛𝜇1𝜇22𝜇1+𝜇23𝑛,(6.15) suggesting (when 𝜇1𝜇2>0) a central limit theorem for 𝑁 as 𝑛 grows large: 1𝑛𝑁𝜇𝑛1𝜇1+𝜇2𝑑𝑛𝒩𝜇0,1𝜇22𝜇1+𝜇23.(6.16)

Case 2. Conversely, assume now that given 𝑁0=𝑚𝑁𝑟𝑑0as𝑟,withprobability𝑚(𝑁=0)=𝜌(𝑚),(6.17) so that boundaries {0,𝑛} are absorbing. Then, the ancestral process (𝐴𝑟;𝑟) possesses an invariant distribution, in that 𝐴𝑟𝑑𝐴as𝑟,independentlyof𝐴0=𝑘[𝑛].(6.18) In terms of moments, the duality formula means that 𝑛𝑚1𝔼𝑛𝐴𝑚=𝜌(𝑚)=𝑚𝑁=0,(6.19)relating 𝑚-factorial moments of 𝑛𝐴 to the extinction probabilities of 𝑁𝑟 given 𝑁0=𝑚. Stated differently, the probability generating function of 𝐴 is (𝑢[0,1]) 𝔼𝑢𝐴=𝑛𝑚=0𝑛𝑚𝜌(𝑚)𝑢𝑛𝑚(1𝑢)𝑚.(6.20) Let 𝜋(𝑖)=(𝐴=𝑖), with 𝝅=𝝅𝑃𝑛. Then, using duality, 𝝆 is the right eigenvector of Π𝑛 associated to the eigenvalue 1𝝆=Π𝑛𝝆. Thus, 𝝆 and 𝝅 are related through 𝝆=Φ𝑛𝝅or𝝅=Φ𝑛1𝝆.(6.21)The knowledge of the extinction probabilities 𝝆of the forward process allows one to compute the invariant measure 𝝅of the dual backward ancestral process (and conversely).

Typical situations where boundaries {0,𝑛} are absorbing to (𝑁𝑟;𝑟) occur when 𝑝(0)=0 and 𝑝(1)=1. The simplest case is the neutral case, but the nonneutral selection and selection with dominance mechanisms or the quadratic mechanism (Examples 3.1, 3.2 and 3.3) are also in this class. For instance,
(i)in the neutral case, 𝜌(𝑚)=1𝑚/𝑛. Thus, 𝜋(𝑖)=(𝑛𝑖)𝑖𝑗=0(1)𝑖𝑗(𝑖𝑗)(𝑗/𝑛)=𝛿𝑖,1 and 𝐴𝑟𝑑1 as 𝑟, the degenerate state reached when the most recent common ancestor (MRCA) is attained,(ii)nondegenerate solutions of 𝐴 are obtained when considering bias mechanisms with 𝑝(0)=0 and 𝑝(1)=1,(iii)consider any biased WF model with 𝑝(0)=0 and 𝑝(1)=1 for which 𝑝(𝑥)𝑥0𝜆𝑥, 𝜆>1. Then, due to large sample asymptotic independence,𝝂𝑛𝑑𝝃,(6.22)where 𝝃 is an iid sequence with 𝜉1𝑑 Poisson (𝜆) (as it can easily be checked by the Poisson limit to the binomial distribution). In this case, the limiting extinction probability of (𝑁𝑟;𝑟) given 𝑁0=𝑚 is lim𝑛𝜌(𝑚)=𝜌𝑚,𝑚=1,2,, where 0<𝜌<1 is the smallest solution to the fixed point equation𝑥=𝑒𝜆(1𝑥),(6.23)where𝜌 is the singleton extinction probability of a supercritical Galton-Watson process with offspring distribution Poisson(𝜆). More precisely, proceeding as in Möhle [9, Theorem 4.5], we have𝑛𝜌𝑚𝜌(𝑚)=𝜌𝑚1𝜌𝑚1+𝜆𝜌2+𝜆(1𝜌)𝜌1(𝜆𝜌)2𝑚,(6.24) showing that the convergence of 𝜌(𝑚) to 𝜌𝑚 is of order 𝑛1. As a result, we get the asymptotic normality 1𝑛𝐴𝑛(1𝜌)𝑑𝑛𝒩0,𝜌(1𝜌)1+𝜆𝜌.(6.25) Intuitively, (1/𝑛)𝔼[𝑛𝐴]=𝜌(1)=1(𝑁=0)𝜌, showing that 𝔼[𝐴]𝑛𝑛(1𝜌) and 1𝔼𝑛(𝑛1)𝑛𝐴𝑛1𝐴=𝜌(2)=2𝑁=0)𝜌2,(6.26) showing that 𝜎2(𝐴)𝑛𝑛𝜌(1𝜌)/(1+𝜆𝜌).
For quadratic Example 3.3, 𝑝(𝑥)=𝑥(1+𝑐𝑐𝑥), with 𝑐[0,1], 𝜆=1+𝑐>1 as soon as 𝑐>0. When 𝑐(0,1], we thus always have asymptotic normality. For Example 3.1, with selection, 𝑝(𝑥)=(1+𝑠)𝑥/(1+𝑠𝑥), with 𝑠>1𝑝(𝑥)𝑥0(1+𝑠)𝑥 and so 𝜆=1+𝑠. We have asymptotic normality only when 𝑠>0, that is, when the fitness is advantageous (corresponding as required to complete monotonicity of corresponding 𝑞=1𝑝).

Note that this asymptotic behavior does not hold for the Lipshitz continuous admissible mechanism 𝑝(𝑥)=𝑥𝛾 of Example 5.1 (or more generally for mechanisms satisfying 𝑝(𝑥)𝑥0𝜃𝑥𝛾, 𝜃>0 as in compound bias Example 5.7) with 𝛾(0,1) because its behavior near 0 is not linear. This puzzling class of models seems to deserve a special study as deviation to normality is expected. We postpone it to a future work.

7. Concluding Remarks

In this note, we focused on discrete nonneutral Wright-Fisher models and on the conditions on the bias probabilities under which the forward branching dynamics is amenable to a dual discrete ancestral coalescent. It was shown that it concerns a large class of nonneutral models involving completely monotone bias probabilities. Several examples were supplied, some standard, some less classical. The Wright-Fisher model with forward binomial transition matrix is a particular instance of the Dirichlet model with Dirichlet-binomial transition matrix. Following the same lines, using the representation of the Dirichlet binomial distribution as a beta mixture of the binomial distribution, it would be interesting to exhibit that the corresponding conditions on the bias mechanism were the starting point to be a forward Dirichlet branching process. Also of particular interest in this respect would be the discrete nonneutral Moran models whose forward transition matrices are simpler because of their tridiagonal Jacobi structure. We hope to be able to consider shortly these cases (and maybe others) in a future work.