A Duality Approach to the Genealogies of Discrete Nonneutral Wright-Fisher Models

Abstract

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 𝐤_𝑛∶=(𝑘₁,…,𝑘_𝑛), 𝝂_𝑛 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 𝝃_𝑛=(𝜉₁,…,𝜉_𝑛) on summing to 𝑛, that is to say, 𝝂=_𝑛𝑑(𝝃_𝑛∣|𝝃_𝑛|=𝑛), where ℙ_𝜃𝜉₁=𝑘=[𝜃]^𝑘𝑘!(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 𝝃_𝑛=(𝜉₁,…,𝜉_𝑛) on summing to 𝑛∶𝝂=_𝑛𝑑(𝝃_𝑛∣|𝝃_𝑛|=𝑛). When 𝑛 is large, using Stirling formula, 𝑛!∼√2𝜋𝑛^𝑛+1/2𝑒^−𝑛; it follows that 𝝂→_{𝑛𝑑𝑛↑∞}𝝃_∞ with joint finite-dimensional law: ℙ(𝝃_𝑛=𝐤_𝑛)=∏^𝑛𝑚=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 𝑁_𝑟(𝑚)=#offspringatgeneration𝑟∈ℕ₊,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 |𝜆₀|≥|𝜆₁|≥⋯≥|𝜆_𝑛| and 1=𝜆₀=𝜆₁>|𝜆₂|.

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 𝛽,𝐛_𝑎∶=(𝑏₁,…,𝑏_𝑎) 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 𝑏₁≥⋯≥𝑏_𝑐≥2, 𝑏_𝑐+1=⋯=𝑏_𝑎=1.

If 𝑐=1 a unique multiple collision occurs of order 𝑏₁≥2.If 𝑏₁=2 a simple binary collision occurs involving only two clusters.If 𝑐>1, simultaneous multiple collisions of orders 𝑏₁≥⋯≥𝑏_𝑐≥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 𝒜₀={(1),…,(𝑚)}, and terminal absorbing state {(1,…,𝑚)}.

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

The Ancestral Count Process
Let 𝐴_𝑟(𝑚)=#anTchesetno:r𝐴sa_𝑟tgeneration𝑟∈ℕ,backward-in-time.(𝑚)=#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𝜈𝑏^𝑙,𝑛_𝑙.(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 𝑃^{(𝑛)𝑏,𝑎}: 𝑃^{(𝑛)𝑏,𝑎}∶=^𝑛𝑎^𝑛𝑏^𝑎_𝑚₌₀(−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 (ℰ₁,ℰ₂), are said to be dual with respect to some real-valued function Φ on the product space ℰ₁×ℰ₂ if ∀𝑥₁∈ℰ₁,∀𝑥₂∈ℰ₂,∀𝑡≥0: 𝔼_𝑥1Φ𝑋^1𝑡,𝑥₂=𝔼_𝑥2Φ𝑥₁,𝑋^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.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 𝑁₀=𝑚. 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 𝑁₀=𝑚. 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 ((ii)i)ΠΠ_𝑛𝑛Φ^1𝑛=Φ^1𝑛𝑃^𝑛,Φ^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)^𝑖−𝑗𝑗𝑖_2𝑛𝑛𝑖,⁻¹(𝑖,𝑗)==((−−11))^{𝑖+𝑗−𝑛}_{𝑖+𝑗−𝑛}𝑛−𝑖𝑗𝑛𝑖𝑛−𝑗𝑖𝑛𝑗.(2.20) The matrix [Φ^1𝑛]⁻¹ is left-lower triangular, while [Φ^2𝑛]⁻¹ is symmetric right-lower triangular. Thus, ((ii)i)ΦΦ^1𝑛_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(𝜆₀,…,𝜆_𝑛), 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 𝑘⟶𝑛𝑝𝑘𝑛inℙ𝜈^1,𝑛+⋯+𝜈_𝑘,𝑛=𝑘^,(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+𝑠)𝑥²+(1+𝑠ℎ)𝑥(1−𝑥)1+𝑠𝑥²+2𝑠ℎ𝑥(1−𝑥).(3.7) In this model, genotype 𝐴𝐴 (resp., 𝐴𝑎 and 𝑎𝑎), with frequency 𝑥² (resp., 2𝑥(1−𝑥) and (1−𝑥)²) 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−𝑥)ℎ−𝑥(2ℎ−1)1+𝑠𝑥²+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−𝑥)², 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−𝜇₂)𝑥+𝜇₁(1−𝑥),(3.10) where (𝜇₁,𝜇₂) are mutation probabilities, satisfying 𝜇₁≤1−𝜇₂. It corresponds to the mutation scheme 𝑎⇄^𝜇1_𝜇2𝐴. To avoid discussions of intermediate cases, we will assume that 𝑝(0)=𝜇₁>0 and 𝑝(1)<1(𝜇₂>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 ∇⁰𝑣(𝑚)=𝑣(𝑚),∇²𝑣(𝑚)=𝑣(𝑚)−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 (𝑢¹,𝑢²) be two continuous functions on [0,1]. Let 𝑢=𝑢¹⋅𝑢². With 𝐮_𝑛 the pointwise product of 𝐮^1𝑛 and 𝐮^2𝑛, assuming both functions (𝑢¹,𝑢²) 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𝑛]⁻¹Π_𝑛Φ^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,=(𝑛)_𝑗⋅𝑛^−𝑖⋅𝑆_𝑖,𝑗,𝑗el≤se.𝑖(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−𝜇₂)𝑥+𝜇₁(1−𝑥) where (𝜇₁,𝜇₂) are mutation probabilities. Then, with 𝜅∶=1−(𝜇₁+𝜇₂), 𝑞(𝑥)=1−𝜇₁−𝜅𝑥 is completely monotone on (0,1) if and only 𝜇₁≤1−𝜇₂(𝜅≥0). In this case, 𝑃_𝑛 is again lower left triangular (a pure death process). We have 𝑃_𝑛(𝑖,𝑗)=⎧⎪⎨⎪⎩𝑛𝑗^𝑗_𝑙=0(−1)^𝑗−𝑙𝑗𝑙𝜇²+𝜅𝑙𝑛^𝑖=∶(𝑛)_𝑗⋅𝑛^−𝑖⋅𝑆^𝜇2_𝑖,𝑗𝜅0,𝑛,𝑗el≤se.𝑖(5.2) in terms of generalized Stirling numbers 𝑆^𝜇2_𝑖,𝑗(𝜅/𝑛). We have 𝑃_𝑛(𝑖,𝑖)=(𝑛)_𝑖(𝜅/𝑛)^𝑖 and the spectrum of 𝑃_𝑛 is real. When 𝜇₁>0, this matrix is substochastic with ∑^𝑛𝑗=0𝑃_𝑛(𝑖,𝑗)=(1−𝜇₁)^𝑖.

A particular case deals with one-way mutations (𝜇₁+𝜇₂>0,𝜇₁⋅𝜇₂=0).

If 𝜇₂=0,𝑃_𝑛(𝑖,𝑗)=(1−𝜇₁)^𝑖⋅(𝑛)_𝑗⋅𝑛^−𝑖⋅𝑆_𝑖,𝑗,𝑗≤𝑖,=0, else. Further, ∑^𝑛𝑗=0𝑃_𝑛(𝑖,𝑗)=(1−𝜇₁)^𝑖<1.If 𝜇₁=0,𝑃_𝑛(𝑖,𝑗)=(𝑛)_𝑗⋅𝑛^−𝑖⋅𝑆^𝜇2_𝑖,𝑗((1−𝜇₂)/𝑛),𝑗≤𝑖,=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 𝑠≥𝑠_ℎ∶=(1−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−𝜇₂)𝑥+𝜇₁(1−𝑥) be the mutational bias mechanism (with 𝜅=1−(𝜇₁+𝜇₂)≥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+𝑠𝑥(5.5) is a mechanism of selection combined with mutational effects. We have 𝑝_𝑀(0)=𝜇₁, 𝑝_𝑀(1)=1−𝜇₂. 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 𝑝¹(𝑥) and 𝑝²(𝑥) are two admissible bias mechanisms, then 𝑝¹(𝑝²(𝑥)) is admissible.

5.3. Joint Bias Effects and Compound Bias

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

Example 5.3 (joint selection and mutational effects). Let 𝑝¹(𝑥)=𝑝_𝑀(𝑥) and 𝑝²(𝑥)=(1+𝑠)𝑥/(1+𝑠𝑥). We get 𝑞(𝑥)=1−𝜇¹−𝜅𝑥(1−𝑥)𝑝(𝑥)=𝜇₁1+𝑠𝑥,+𝑥𝑠+1−𝜇₁+𝜅−𝜅𝑥²1+𝑠𝑥,(5.8) with 𝑝(0)=𝜇₁, 𝑝(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−𝑥) and 𝑝²(𝑥)=𝑥. We get 𝑝𝑞((𝑥𝑥))==𝜇(1₁−𝑥)(1−𝜇₁−𝜅𝑥),+𝑥1−𝜇₁+𝜅(1−𝑥),(5.9) with 𝑝(0)=𝜇₁, 𝑝(1)=1. When 𝜇₁=0 (one-way mutations), we recover the quadratic mechanism with curvature parameter 𝑐=1−𝜇₂. 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 𝑝¹(𝑥) and 𝑝²(𝑥).

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, 𝑞_𝑀(𝑥)==11−−𝑝𝑝_𝑀𝑀(𝑥)=1−𝑝_𝑀(𝑝(𝑥))(1−𝑞(𝑥)),(5.11) and so 𝜙(𝑥)=1−𝑝_𝑀(1−𝑥)=1−(1−𝜇₂)(1−𝑥)−𝜇₁𝑥=𝜇₂+𝜅𝑥 which is absolutely monotone as soon as 𝜅=1−(𝜇₁+𝜇₂)≥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𝑁₀=𝑚≥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 𝑘𝑛⁻¹𝔼𝑛−𝑘𝑁^∞=𝜌^∞(𝑘)=ℙ_𝑘(𝐴_∞=0),(6.6)relating 𝑘-factorial moments of 𝑛−𝑁_∞ to the extinction probabilities of 𝐴_𝑟 given 𝐴₀=𝑘. 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 𝐴₀=𝑘 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^𝑛𝑚,(𝑖,𝑗)==((−−11))^{𝑖+𝑗−𝑛}_{𝑖+𝑗−𝑛}𝑛−𝑖𝑗𝑛𝑖𝑛−𝑗𝑖𝑛𝑗.(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).

Example
Consider the discrete WF model with mutations of Example 3.4. In this case, 𝑁−→_𝑟𝑑𝑁_∞ as 𝑟↑∞, regardless of 𝑁₀=𝑚 and (𝑁_𝑟;𝑟∈ℕ) has an invariant measure which is difficult to compute. Looking at the backward process, the matrix 𝑃_𝑛 is substochastic (if 𝜇₁>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)=𝜇₂/(𝜇₁+𝜇₂) so that 𝔼[𝑁_∞]=𝑛𝜇₁/(𝜇₁+𝜇₂); from (5.2) we also have 𝜌_∞(2)=𝜇²𝑛𝜇₂(1+𝜅)+𝜅²=(11−𝜅)𝑛−(𝑛−1)𝜅²𝑛(𝑛−1)𝑛(𝑛−1)−𝑛(2𝜇𝑛₁−+1𝜇)𝜇¹₂+𝔼𝑁^2∞,(6.14) allowing to compute 𝔼[𝑁^2∞] and then the variance of 𝑁_∞. We get 𝜎²(𝑁_∞)=𝑛²𝜇₁𝜇₂𝜇₁+𝜇₂²2𝑛𝜇₁+𝜇₂+1+𝑜(𝑛)∼^𝑛↑∞𝜇₁𝜇₂2𝜇₁+𝜇₂³𝑛,(6.15) suggesting (when 𝜇₁𝜇₂>0) a central limit theorem for 𝑁_∞ as 𝑛 grows large: 1√𝑛𝑁^∞−𝑛𝜇¹𝜇₁+𝜇₂^𝑑⟶_𝑛↑∞𝒩0,𝜇¹𝜇₂2𝜇₁+𝜇₂³.(6.16)

Case 2. Conversely, assume now that given 𝑁₀=𝑚𝑁⟶_𝑟𝑑0as𝑟↑∞,withprobabilityℙ_𝑚(𝑁_∞=0)=∶𝜌_∞(𝑚),(6.17) so that boundaries {0,𝑛} are absorbing. Then, the ancestral process (𝐴_𝑟;𝑟∈ℕ) possesses an invariant distribution, in that 𝐴⟶_𝑟𝑑𝐴_∞as𝑟↑∞,independentlyof𝐴₀=𝑘∈[𝑛].(6.18) In terms of moments, the duality formula means that 𝑚𝑛⁻¹𝔼𝑛−𝑚𝐴^∞=𝜌^∞(𝑚)=ℙ_𝑚𝑁_∞=0,(6.19)relating 𝑚-factorial moments of 𝑛−𝐴_∞ to the extinction probabilities of 𝑁_𝑟 given 𝑁₀=𝑚. Stated differently, the probability generating function of 𝐴_∞ is (𝑢∈[0,1]) 𝔼𝑢^𝐴∞=^𝑛_𝑚₌₀𝑚𝑛𝜌^∞(𝑚)𝑢^𝑛−𝑚(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).

Examples
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 𝑁₀=𝑚 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𝑛𝜌^𝑚−𝜌_∞(𝑚)=𝜌^𝑚11+−𝜆𝜌𝜌𝑚²+𝜆(1−𝜌)𝜌1−(𝜆𝜌)²𝑚,(6.24) showing that the convergence of 𝜌_∞(𝑚) to 𝜌^𝑚 is of order 𝑛⁻¹. As a result, we get the asymptotic normality 1√𝑛𝐴^∞−𝑛(1−𝜌)⟶^𝑑_𝑛↑∞𝒩0,𝜌(1−𝜌)1+𝜆𝜌.(6.25) Intuitively, (1/𝑛)𝔼[𝑛−𝐴_∞]=𝜌_∞(1)=ℙ₁(𝑁_∞=0)→𝜌, showing that 𝔼[𝐴_∞]∼_𝑛↑∞𝑛(1−𝜌) and 1𝑛(𝑛−1)𝔼𝑛−𝐴^∞𝑛−1−𝐴_∞=𝜌_∞(2)=ℙ₂𝑁_∞=0)⟶𝜌²,(6.26) showing that 𝜎²(𝐴_∞)∼_𝑛↑∞𝑛𝜌(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.