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

Authors: Thierry E. Huillet
Journal: Journal of Probability and Statistics
 (2009)

---

## 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.


---

## Body

## 1. Introduction

Forward
evolution  of large populations in genetics has a long history, starting in the1920s; 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. Section2 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
generationsr∈ℤ. Assume the population size is constant, say n (n individuals (or genes)) over generations.
Assume the random reproduction law at generation 0 is νn:=(ν1,n,…,νn,n), satisfying(2.1)∑m=1nνm,n=n. Here, νm,n is the number of offspring of gene m. We avoid the trivial case: νm,n=1,m=1,…,n. One iterates the reproduction over
generations, while imposing the following additional assumptions:(i)
exchangeability:(ν1,n,…,νn,n)=d(νσ(1),n,…,νσ(n),n), for all permutations σ∈𝒮n;(ii)
time-homogeneity: reproduction laws are
independent and identically distributed (iid) at each generationr∈ℤ. This model, therefore, consists of a
conservative-conditioned branching Galton-Watson process in [n]ℤ,
where [n]:={0,1,…,n} (see Karlin-McGregor [18]).Famous reproduction laws are as follows.Example 2.1.
The multinomial Dirichlet model:νn~d Multin-Dirichlet (n;θ),
where θ>0 is a disorder parameter. With kn:=(k1,…,kn), νn admits the following joint exchangeable
distribution on the simplex |kn|:=∑m=1nkm=n: (2.2)ℙθ(νn=kn)=n![nθ]n∏m=1n[θ]kmkm!, where [θ]k=θ(θ+1)⋯(θ+k−1) is the rising factorial of θ.
This distribution can be obtained by conditioning n independent mean 1 Pòlya distributed random variables ξn=(ξ1,…,ξn) on summing to n,
that is to say, νn=d(ξn∣|ξn|=n), where (2.3)ℙθ(ξ1=k)=[θ]kk!(1+θ)−k(θ1+θ)θ,k∈ℕ. When θ↑∞,
this distribution reduces to the Wright-Fisher model for which νn~d Multin(n;1/n,…,1/n). Indeed, νn admits the following joint exchangeable
multinomial distribution on the simplex |kn|=n: (2.4)ℙ∞(νn=kn)=n!⋅n−n∏m=1nkm!. This distribution can be
obtained by conditioning n independent mean 1 Poisson distributed random variables ξn=(ξ1,…,ξn) on summing to n:νn=d(ξn∣|ξn|=n).
When n is large, using Stirling formula, n!~2πnn+1/2e−n; it follows that νn→n↑∞dξ∞ with joint finite-dimensional law: ℙ(ξn=kn)=∏m=1n(e−1/km!)=(e−n/∏m=1nkm!) on ℕn. Thanks to the product form of all
finite-dimensional laws of ξ∞,
we get an asymptotic independence property of νn.Example 2.2.
In the Moran model,νn~d 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 sizem from [n]:={0,1,…,n} at generation 0. Let (2.5)Nr(m)=#offspring  at  generationr∈ℕ+,forward  in  time. This sibship process is a
discrete-time homogeneous Markov chain, with transition
probability (2.6)ℙ(Nr+1(m)=k′∣Nr(m)=k)=ℙ(ν1,n+⋯+νk,n=k′). It is a martingale, with state space {0,…,n},
initial state m,
absorbing states {0,n}, and transient states {1,…,n−1}. The first hitting time of boundaries {0,n},
which is τ(m)=τ{0}(m)⋀τ{n}(m), is finite with probability 1 and has finite mean. Omitting reference to any
specific initial condition m,
the process (Nr;r∈ℕ) has the transition matrix Πn with entries Πn(k,k′)=ℙ(ν1,n+⋯+νk,n=k′) given by (2.6). We have Πn(0,k′)=δ0,k′ and Πn(n,k′)=δn,k′,
and Πn is not irreducible. However, Πn is aperiodic and (apart from absorbing states)
cannot be broken down into noncommunicating subsets; as a result it is
diagonalizable, with eigenvalues |λ0|≥|λ1|≥⋯≥|λn| and 1=λ0=λ1>|λ2|.Example 2.3 (Dirichlet binomial).
WithU a (0,1)-valued random variable with density beta(kθ,(n−k)θ)(2.7)ℙ(ν1,n+⋯+νk,n=k′)=(nk′)[kθ]k′[(n−k)θ]n−k′[nθ]n=𝔼[(nk′)Uk′(1−U)n−k′], which is a beta mixture of the
binomial distribution Bin(n,u).Example 2.4.
The Wright-Fisher model has a Bin(n,k/n) transition matrix (2.8)ℙ(Nr+1(m)=k′∣Nr(m)=k)=(nk′)(kn)k′(1−kn)n−k′.Remark 2.5 (statistical symmetry).
Due to exchangeability of the reproduction law, neutral
models are symmetric in the following sense. The transition probabilities ofN¯r(m):=n−Nr(m) are equal to the transition probabilities of Nr(m).


### 2.3. Backward in Time Process (neutral Coalescent)

Thecoalescent  backward  process can be defined as follows. Take a subsample of
size m from [n] at generation 0. Identify two individuals from [m] 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 [m].
Define the induced ancestral backward process (2.9)𝒜r(m)∈ℰm={equivalence  classes  (partitions)  of[m]},r∈ℕ,backwardintime. The ancestral process is a discrete-time Markov chain
with transition probability: (2.10)ℙ(𝒜r+1(m)=α∣𝒜r(m)=β)=Pβ;α;with(α,β)∈ℰm,α⊆β, where, with a=|α|= number of equivalence classes of α,b=|β|= number of equivalence classes of β,ba:=(b1,…,ba) clusters sizes of β, and (m)a:=m(m−1)⋯(m−a+1) a falling factorial (2.11)Pβ;α=Pb;a(n)(ba)=(n)a(n)b𝔼(∏l=1a(νl,n)bl) is the probability of a ba-merger. This is the probability that b randomly chosen individuals out of n have a≤b distinct parents, c merging classes and cluster sizes b1≥⋯≥bc≥2, bc+1=⋯=ba=1.Ifc=1 a unique multiple collision occurs of order b1≥2.Ifb1=2 a simple binary collision occurs involving
only two clusters.Ifc>1,
simultaneous multiple collisions of orders b1≥⋯≥bc≥2 occur. Thus, the jump's height of a transition b→a is b−a=∑i=1c(bi−1), corresponding to a partition of integer b−a into c summands, each ≥1.The chain's state-space is{equivalence relations on (partitions of) {1,…,m}};
it has dimension Bm:=∑l=0mSm,l (a Bell number), where Sm,l are the second-kind Stirling numbers.The chain has initial state𝒜0={(1),…,(m)},
and terminal absorbing state {(1,…,m)}.Examples
From the Dirichlet Example2.3, we get, Pb;a(n)(ba)=((n)a/[nθ]b)∏m=1a[θ]bm.
From the WF Example2.4: In this case, Pb;a(n)(ba)=(n)a/nb is the uniform distribution on {ba:b1+⋯+ba=b}.The Ancestral Count Process
Let(2.12)Ar(m)=#ancestors at generationr∈ℕ,backward-in-time.Then:Ar(m)=#blocks  of𝒜r(m). The backward ancestral count
process is a discrete-time Markov chain with transition probabilities [17, 19](2.13)ℙ(Ar+1(m)=a∣Ar(m)=b)=Pb,a(n):=b!a!∑b1,…,ba∈ℕ+b1+⋯+ba=bPb;a(n)(ba)b1!⋯ba!.=(na)(nb)∑b1,…,ba∈ℕ+b1+⋯+ba=b𝔼(∏l=1a(νl,nbl)). This Markov chain has
state-space {0,…,m},
initial state m,
absorbing states {0,1}. The process (Ar;r∈ℕ) has the transition matrix Pn with entries Pn(b,a)=Pb,a(n) given by (2.13). Note, by inclusion-exclusion principle, the
alternative alternating expression of Pb,a(n): (2.14)Pb,a(n):=(na)(nb)∑m=0a(−1)a−m(am)𝔼((ν1,n+⋯+νm,nb)).


### 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(Xt1,Xt2;t≥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  ∀x1∈ℰ1,∀x2∈ℰ2,∀t≥0: (2.15)𝔼x1Φ(Xt1,x2)=𝔼x2Φ(x1,Xt2).We  then recall
basic examples of dual processes from the neutral and exchangeable population
models (Möhle[10]). The neutral forward and backward processes (Nr,Ar;r∈ℕ) introduced in the two preceding subsections
are dual with respect to the hypergeometric sampling without replacement
kernels (2.16)(i)Φn1(m,k)=(mk)(nk),(ii)Φn2(m,k)=(n−mk)(nk)on{0,…,n}2. Namely, (i) reads(2.17)𝔼m[(Nrk)(nk)]=𝔼k[(mAr)(nAr)]=𝔼k[(n−Arn−m)(nn−m)]. Call type A individuals the descendants of the m first chosen individuals (allele A) in the study of the forward process; type a individuals are the remaining ones (allele a). The left-hand-side of the above equality is
an expression of the probability that k-samples (without
replacement) from population of size Nr at time r are all of type A,
given N0=m. If these k-samples are all
descendants of Ar ancestors at time −r,
this probability must be equal to the probability that a (n−m)-sample from population of size Ar at time −r are all of type a.
This is the meaning of the right-hand-side.Moreover, (ii) reads(2.18)𝔼m[(n−Nrk)(nk)]=𝔼k[(n−mAr)(nAr)]=𝔼k[(n−Arm)(nm)]. The left-hand-side is the
probability that a k-sample (without replacement) from population
of size Nr at time r are all of type a,
given N0=m. If these k-samples are all
descendants of Ar ancestors at time −r,
this probability must be equal to the probability that m-samples from
population of size Ar at time −r are themselves all of type a.WithPn′ the transpose of Pn,
a one-step (r=1) version of these formulae is (2.19)(i)ΠnΦn1=Φn1Pn′,(ii)ΠnΦn2=Φn2Pn′, where (Φn1,Φn2) are n×n matrices with entries Φn1(m,k) and Φn2(m,k), respectively, and (Πn,Pn) the transition matrices of forward and
backward processes. Note that the matrix Φn2 is symmetric and left-upper triangular. The
matrices Φn1 and Φn2 are both invertible, with respective
entries (2.20)[Φn1]−1(i,j)=(−1)i−j(ij)(ni),[Φn2]−1(i,j)=(−1)i+j−n(in−j)(ni)=(−1)i+j−n(jn−i)(nj). The matrix [Φn1]−1 is left-lower triangular, while [Φn2]−1 is symmetric right-lower triangular.
Thus, (2.21)(i)[Φn1]−1ΠnΦn1=Pn′,(ii)[Φn2]−1ΠnΦn2=Pn′. In any case, being similar
matrices, Πn and Pn′ (or Pn) both share the same eigenvalues. If Rn diagonalizing Πn is known so that Rn−1ΠnRn=Λn:=diag(λ0,…,λn), the diagonal matrix of the eigenvalues of Πn,
then, with Φn=Φn1 or Φn2, R˜n:=Φn−1Rn diagonalizes Pn′ and is obtained for free (and conversely). Rn is the matrix whose columns are the
right-eigenvectors of Πn and R˜n is the matrix whose columns (rows) are the
right eigenvectors (left-eigenvectors) of Pn′ (of Pn). Similarly, if Ln is the matrix whose rows are the left
eigenvectors of Πn, L˜n:=LnΦn is the matrix whose rows (columns) are the
left eigenvectors (right eigenvectors) of Pn′ (of Pn). With lk′ the kth row of Ln and rk the kth column of Rn, the spectral decomposition of Πn is (2.22)Πnr=∑k=0nλkrrklk′lk′rk,r∈ℕ, whereas, with l˜k the kth column of L˜n and r˜k′ the kth row of R˜n, the one of Pn reads (2.23)Pnr=∑k=0nλkrl˜kr˜k′r˜k′l˜k=∑k=0nλkrΦn′lk(Φn−1rk)′(Φn−1rk)′Φn′lk,r∈ℕ. 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.


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

Consider a population with nonoverlapping
generationsr∈ℤ. Assume the population size is constant, say n (n individuals (or genes)) over generations.
Assume the random reproduction law at generation 0 is νn:=(ν1,n,…,νn,n), satisfying(2.1)∑m=1nνm,n=n. Here, νm,n is the number of offspring of gene m. We avoid the trivial case: νm,n=1,m=1,…,n. One iterates the reproduction over
generations, while imposing the following additional assumptions:(i)
exchangeability:(ν1,n,…,νn,n)=d(νσ(1),n,…,νσ(n),n), for all permutations σ∈𝒮n;(ii)
time-homogeneity: reproduction laws are
independent and identically distributed (iid) at each generationr∈ℤ. This model, therefore, consists of a
conservative-conditioned branching Galton-Watson process in [n]ℤ,
where [n]:={0,1,…,n} (see Karlin-McGregor [18]).Famous reproduction laws are as follows.Example 2.1.
The multinomial Dirichlet model:νn~d Multin-Dirichlet (n;θ),
where θ>0 is a disorder parameter. With kn:=(k1,…,kn), νn admits the following joint exchangeable
distribution on the simplex |kn|:=∑m=1nkm=n: (2.2)ℙθ(νn=kn)=n![nθ]n∏m=1n[θ]kmkm!, where [θ]k=θ(θ+1)⋯(θ+k−1) is the rising factorial of θ.
This distribution can be obtained by conditioning n independent mean 1 Pòlya distributed random variables ξn=(ξ1,…,ξn) on summing to n,
that is to say, νn=d(ξn∣|ξn|=n), where (2.3)ℙθ(ξ1=k)=[θ]kk!(1+θ)−k(θ1+θ)θ,k∈ℕ. When θ↑∞,
this distribution reduces to the Wright-Fisher model for which νn~d Multin(n;1/n,…,1/n). Indeed, νn admits the following joint exchangeable
multinomial distribution on the simplex |kn|=n: (2.4)ℙ∞(νn=kn)=n!⋅n−n∏m=1nkm!. This distribution can be
obtained by conditioning n independent mean 1 Poisson distributed random variables ξn=(ξ1,…,ξn) on summing to n:νn=d(ξn∣|ξn|=n).
When n is large, using Stirling formula, n!~2πnn+1/2e−n; it follows that νn→n↑∞dξ∞ with joint finite-dimensional law: ℙ(ξn=kn)=∏m=1n(e−1/km!)=(e−n/∏m=1nkm!) on ℕn. Thanks to the product form of all
finite-dimensional laws of ξ∞,
we get an asymptotic independence property of νn.Example 2.2.
In the Moran model,νn~d 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 sizem from [n]:={0,1,…,n} at generation 0. Let (2.5)Nr(m)=#offspring  at  generationr∈ℕ+,forward  in  time. This sibship process is a
discrete-time homogeneous Markov chain, with transition
probability (2.6)ℙ(Nr+1(m)=k′∣Nr(m)=k)=ℙ(ν1,n+⋯+νk,n=k′). It is a martingale, with state space {0,…,n},
initial state m,
absorbing states {0,n}, and transient states {1,…,n−1}. The first hitting time of boundaries {0,n},
which is τ(m)=τ{0}(m)⋀τ{n}(m), is finite with probability 1 and has finite mean. Omitting reference to any
specific initial condition m,
the process (Nr;r∈ℕ) has the transition matrix Πn with entries Πn(k,k′)=ℙ(ν1,n+⋯+νk,n=k′) given by (2.6). We have Πn(0,k′)=δ0,k′ and Πn(n,k′)=δn,k′,
and Πn is not irreducible. However, Πn is aperiodic and (apart from absorbing states)
cannot be broken down into noncommunicating subsets; as a result it is
diagonalizable, with eigenvalues |λ0|≥|λ1|≥⋯≥|λn| and 1=λ0=λ1>|λ2|.Example 2.3 (Dirichlet binomial).
WithU a (0,1)-valued random variable with density beta(kθ,(n−k)θ)(2.7)ℙ(ν1,n+⋯+νk,n=k′)=(nk′)[kθ]k′[(n−k)θ]n−k′[nθ]n=𝔼[(nk′)Uk′(1−U)n−k′], which is a beta mixture of the
binomial distribution Bin(n,u).Example 2.4.
The Wright-Fisher model has a Bin(n,k/n) transition matrix (2.8)ℙ(Nr+1(m)=k′∣Nr(m)=k)=(nk′)(kn)k′(1−kn)n−k′.Remark 2.5 (statistical symmetry).
Due to exchangeability of the reproduction law, neutral
models are symmetric in the following sense. The transition probabilities ofN¯r(m):=n−Nr(m) are equal to the transition probabilities of Nr(m).


## 2.3. Backward in Time Process (neutral Coalescent)

Thecoalescent  backward  process can be defined as follows. Take a subsample of
size m from [n] at generation 0. Identify two individuals from [m] 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 [m].
Define the induced ancestral backward process (2.9)𝒜r(m)∈ℰm={equivalence  classes  (partitions)  of[m]},r∈ℕ,backwardintime. The ancestral process is a discrete-time Markov chain
with transition probability: (2.10)ℙ(𝒜r+1(m)=α∣𝒜r(m)=β)=Pβ;α;with(α,β)∈ℰm,α⊆β, where, with a=|α|= number of equivalence classes of α,b=|β|= number of equivalence classes of β,ba:=(b1,…,ba) clusters sizes of β, and (m)a:=m(m−1)⋯(m−a+1) a falling factorial (2.11)Pβ;α=Pb;a(n)(ba)=(n)a(n)b𝔼(∏l=1a(νl,n)bl) is the probability of a ba-merger. This is the probability that b randomly chosen individuals out of n have a≤b distinct parents, c merging classes and cluster sizes b1≥⋯≥bc≥2, bc+1=⋯=ba=1.Ifc=1 a unique multiple collision occurs of order b1≥2.Ifb1=2 a simple binary collision occurs involving
only two clusters.Ifc>1,
simultaneous multiple collisions of orders b1≥⋯≥bc≥2 occur. Thus, the jump's height of a transition b→a is b−a=∑i=1c(bi−1), corresponding to a partition of integer b−a into c summands, each ≥1.The chain's state-space is{equivalence relations on (partitions of) {1,…,m}};
it has dimension Bm:=∑l=0mSm,l (a Bell number), where Sm,l are the second-kind Stirling numbers.The chain has initial state𝒜0={(1),…,(m)},
and terminal absorbing state {(1,…,m)}.Examples
From the Dirichlet Example2.3, we get, Pb;a(n)(ba)=((n)a/[nθ]b)∏m=1a[θ]bm.
From the WF Example2.4: In this case, Pb;a(n)(ba)=(n)a/nb is the uniform distribution on {ba:b1+⋯+ba=b}.The Ancestral Count Process
Let(2.12)Ar(m)=#ancestors at generationr∈ℕ,backward-in-time.Then:Ar(m)=#blocks  of𝒜r(m). The backward ancestral count
process is a discrete-time Markov chain with transition probabilities [17, 19](2.13)ℙ(Ar+1(m)=a∣Ar(m)=b)=Pb,a(n):=b!a!∑b1,…,ba∈ℕ+b1+⋯+ba=bPb;a(n)(ba)b1!⋯ba!.=(na)(nb)∑b1,…,ba∈ℕ+b1+⋯+ba=b𝔼(∏l=1a(νl,nbl)). This Markov chain has
state-space {0,…,m},
initial state m,
absorbing states {0,1}. The process (Ar;r∈ℕ) has the transition matrix Pn with entries Pn(b,a)=Pb,a(n) given by (2.13). Note, by inclusion-exclusion principle, the
alternative alternating expression of Pb,a(n): (2.14)Pb,a(n):=(na)(nb)∑m=0a(−1)a−m(am)𝔼((ν1,n+⋯+νm,nb)).


## 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(Xt1,Xt2;t≥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  ∀x1∈ℰ1,∀x2∈ℰ2,∀t≥0: (2.15)𝔼x1Φ(Xt1,x2)=𝔼x2Φ(x1,Xt2).We  then recall
basic examples of dual processes from the neutral and exchangeable population
models (Möhle[10]). The neutral forward and backward processes (Nr,Ar;r∈ℕ) introduced in the two preceding subsections
are dual with respect to the hypergeometric sampling without replacement
kernels (2.16)(i)Φn1(m,k)=(mk)(nk),(ii)Φn2(m,k)=(n−mk)(nk)on{0,…,n}2. Namely, (i) reads(2.17)𝔼m[(Nrk)(nk)]=𝔼k[(mAr)(nAr)]=𝔼k[(n−Arn−m)(nn−m)]. Call type A individuals the descendants of the m first chosen individuals (allele A) in the study of the forward process; type a individuals are the remaining ones (allele a). The left-hand-side of the above equality is
an expression of the probability that k-samples (without
replacement) from population of size Nr at time r are all of type A,
given N0=m. If these k-samples are all
descendants of Ar ancestors at time −r,
this probability must be equal to the probability that a (n−m)-sample from population of size Ar at time −r are all of type a.
This is the meaning of the right-hand-side.Moreover, (ii) reads(2.18)𝔼m[(n−Nrk)(nk)]=𝔼k[(n−mAr)(nAr)]=𝔼k[(n−Arm)(nm)]. The left-hand-side is the
probability that a k-sample (without replacement) from population
of size Nr at time r are all of type a,
given N0=m. If these k-samples are all
descendants of Ar ancestors at time −r,
this probability must be equal to the probability that m-samples from
population of size Ar at time −r are themselves all of type a.WithPn′ the transpose of Pn,
a one-step (r=1) version of these formulae is (2.19)(i)ΠnΦn1=Φn1Pn′,(ii)ΠnΦn2=Φn2Pn′, where (Φn1,Φn2) are n×n matrices with entries Φn1(m,k) and Φn2(m,k), respectively, and (Πn,Pn) the transition matrices of forward and
backward processes. Note that the matrix Φn2 is symmetric and left-upper triangular. The
matrices Φn1 and Φn2 are both invertible, with respective
entries (2.20)[Φn1]−1(i,j)=(−1)i−j(ij)(ni),[Φn2]−1(i,j)=(−1)i+j−n(in−j)(ni)=(−1)i+j−n(jn−i)(nj). The matrix [Φn1]−1 is left-lower triangular, while [Φn2]−1 is symmetric right-lower triangular.
Thus, (2.21)(i)[Φn1]−1ΠnΦn1=Pn′,(ii)[Φn2]−1ΠnΦn2=Pn′. In any case, being similar
matrices, Πn and Pn′ (or Pn) both share the same eigenvalues. If Rn diagonalizing Πn is known so that Rn−1ΠnRn=Λn:=diag(λ0,…,λn), the diagonal matrix of the eigenvalues of Πn,
then, with Φn=Φn1 or Φn2, R˜n:=Φn−1Rn diagonalizes Pn′ and is obtained for free (and conversely). Rn is the matrix whose columns are the
right-eigenvectors of Πn and R˜n is the matrix whose columns (rows) are the
right eigenvectors (left-eigenvectors) of Pn′ (of Pn). Similarly, if Ln is the matrix whose rows are the left
eigenvectors of Πn, L˜n:=LnΦn is the matrix whose rows (columns) are the
left eigenvectors (right eigenvectors) of Pn′ (of Pn). With lk′ the kth row of Ln and rk the kth column of Rn, the spectral decomposition of Πn is (2.22)Πnr=∑k=0nλkrrklk′lk′rk,r∈ℕ, whereas, with l˜k the kth column of L˜n and r˜k′ the kth row of R˜n, the one of Pn reads (2.23)Pnr=∑k=0nλkrl˜kr˜k′r˜k′l˜k=∑k=0nλkrΦn′lk(Φn−1rk)′(Φn−1rk)′Φn′lk,r∈ℕ. 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(3.1)k→np(kn)inℙ(ν1,n+⋯+νk,n=k′), where (3.2)p(x):x∈(0,1)→(0,1)is  continuous,  increasing,withp(0)=0,p(1)=1.p(x) is the state-dependent Bernoulli bias probability
different from identity x (as in neutral case).When particularized to the WF model, this leads to the
biased transition probabilities(3.3)ℙ(Nr+1(m)=k′∣Nr(m)=k)=(nk′)p(kn)k′(1−p(kn))n−k′. In this binomial n-sampling with replacement model, a type A individual is drawn with probability p(k/n) which is different from the uniform
distribution k/n due to bias effects.From this, we conclude (a symmetry breaking property) the following:
the transition probabilities ofN¯r(m):=n−Nr(m), r∈ℕ are (3.4)Bin(n,1−p(1−kn))≠Bin(n,p(kn)), and so, no longer coincide with
the ones of (Nr(m);r∈ℕ). The process Nr(m), r∈ℕ no longer is a martingale. Rather, if x→p(x) is concave (convex), Nr(m), r∈ℕ is a submartingale (supermartingale) because 𝔼(Nr+1(m)∣Nr(m))=np(Nr(m)/n)≥Nr(m) (resp., ≤Nr(m)).In the binomial neutral Wright-Fisher transition
probabilities, we replaced the success probabilityk/n by a more general function p(k/n).
However, this replacement leaves open the question what model is in the
background and what quantity the process (Nr,r∈ℕ) 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, νn~d Multin(n;πn),
where πn:=(π1,n,…,πn,n) and πm,n=p(m/n)−p((m−1)/n),m=1,…,n. We note that under our hypothesis, (3.5)∑m=1nπm,n=p(1)−p(0)=1. Due to its asymmetry, the law of
the biased νn no longer is exchangeable.We now recall some well-known bias examples arising in
population genetics.Example 3.1 (homographic model, selection).
Assume(3.6)p(x)=(1+s)x1+sx, where s>−1 is a selection parameter. This model arises
when gene A (resp., a), with frequency x (resp., 1−x), has fitness 1+s (resp., 1). The case s>0 arises when gene of type A is selectively advantageous, whereas it is
disadvantageous when s∈(−1,0).Example 3.2 (selection with dominance).
Assume(3.7)p(x)=(1+s)x2+(1+sh)x(1−x)1+sx2+2shx(1−x). In this model, genotype AA (resp., Aa and aa), with frequency x2 (resp., 2x(1−x) and (1−x)2) has fitness 1+s (resp., 1+sh and 1). h is a measure of the degree of dominance of
heterozygote Aa.
We impose s>−1 and sh>−1. Note that the latter quantity can be put into
the canonical form of deviation to neutrality: (3.8)p(x)=x+sx(1−x)h−x(2h−1)1+sx2+2shx(1−x), where the ratio appearing in the
right-hand-side is the ratio of the difference of marginal fitnesses of A and a to their mean fitness. The case h=1/2 corresponds to balancing selection with p(x)=x+(s/2)(x(1−x)/(1+sx)).Example 3.3 (quadratic model).
Withc∈[−1,1], a curvature parameter, one may
choose (3.9)p(x)=x(1+c−cx). If c=1, p(x)=x(2−x)=1−(1−x)2,
this bias appears in a discrete 2-sex population model [9, 10]. We will give below an interpretation of
this quadratic model when c∈(0,1] in terms of a joint one-way mutations and
neutrality effects model.
We can relax the assumptionp(0)=0,p(1)=1 by assuming 0≤p(0)≤p(1)≤1, p(1)−p(0)∈[0,1).Example 3.4 (affine model).
Take, for
example,(3.10)p(x)=(1−μ2)x+μ1(1−x), where (μ1,μ2) are mutation probabilities, satisfying μ1≤1−μ2. It corresponds to the mutation scheme a⇄μ1μ2A.
To avoid discussions of intermediate cases, we will assume that p(0)=μ1>0 and p(1)<1(μ2>0). In this case, the matrix Πn is irreducible and even primitive and all
states of this Markov chain are now recurrent. We have ℙ(Nr+1>0∣Nr=0)=1−(1−p(0))n>0 and ℙ(Nr+1<n∣Nr=n)=1−p(1)n>0 and the boundaries {0} and {n} no longer are strictly absorbing as there is a
positive reflection probability inside the domain {0,1,…,n}.
For reasons to appear now, we will be only interested
in functionsq such that q(x):=1−p(x) is a completely monotone function (CM) on (0,1), that is, satisfying (3.11)(−1)lq(l)(x)≥0,∀x∈(0,1), for all order-l derivatives q(l) of q,l≥0.
If p(x) is such that q is CM, we will call it an admissible bias
mechanism.


## 4. Nonneutral Wright-Fisher Models and Duality

Preliminaries.
Letvn:=(v(0),v(1),…,v(n)) be a (n+1)-vector of [0,1]-valued numbers. Define the backward difference
operator ∇ acting on vn by: ∇v(m)=v(m)−v(m−1),m=1,…,n. We have ∇0v(m)=v(m),∇2v(m)=v(m)−2v(m−1)+v(m−2), and so forth, and, starting from the endpoint v(n)(4.1)∇jv(m)|m=n=∑l=0j(−1)j−l(jl)v(n−l),j=0,…,n. Let u be some continuous function [0,1]→[0,1]. Consider the (n+1)-vector un:=(u(0/n),…,u(m/n),…,u(n/n)),
sampling u at points m/n. The function u is said to be ∇-completely monotonic if (−1)j∇ju(m/n)|m=n≥0,
for all j=0,…,n and all n≥0. Let (u1,u2) be two continuous functions on [0,1].
Let u=u1⋅u2.
With un the pointwise product of un1 and un2,
assuming both functions (u1,u2) to be ∇-completely monotonic, so will be u, by the Leibniz rule. In particular, if u is ∇-completely monotonic, so will be its integral
powers ui,i∈ℕ.
Our main result is as follows.Theorem 4.2.
Consider
a nonneutral WF forward model(Nr;r∈ℕ)  on  {0,…,n},
with continuous, nondecreasing bias  p(x),  satisfying (4.2)0≤p(0)≤p(1)≤1,p(1)−p(0)∈[0,1].
This process has forward
transition matrix (4.3)Πn(k,k′)=ℙ(ν1,n+⋯+νk,n=k′)=(nk′)p(kn)k′(1−p(kn))n−k′. There exists a Markov chain (Ar;r∈ℕ) on {0,…,n}  such that  (Nr,Ar;r∈ℕ)  are dual with respect to  Φn2(m,k)=(n−mk)/(nk)  if and only if  x→q(x)=1−p(x)  is completely monotone on  (0,1).
In this case, the transition probability matrix of  (Ar;r∈ℕ) is (4.4)Pn(i,j)=(nj)∑l=0j(−1)j−l(jl)q(1−ln)i≥0,wherePn  is a stochastic matrix if and only if  p(0)=0;  else, if  p(0)>0,
the matrix  Pn  is substochastic.Proof.
Developing[Φn2]−1ΠnΦn2=Pn′, we easily obtain (4.5)Pn′(j,i)=Pn(i,j)=(nj)∑l=0j(−1)j−l(jl)[1−p(n−ln)]i=(nj)(−1)j∇j(q(mn)i)|m=n. This entry is nonnegative if and
only if (−1)j∇j(q(m/n)i)|m=n≥0,
for all i,j=0,…,n. But due to the above argument on ∇-complete monotonicity of integral powers, this
will be the case if and only if (−1)j∇j(q(m/n))|m=n≥0,
for all j=0,…,n.
As this must be true for arbitrary value of population size n,
function q 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 x→q(x):=1−p(x) is a completely monotone function on (0,1) in the sense that (4.6)(−1)lq(l)(x)≥0,∀x∈(0,1),l∈ℕ. Next, since (I−∇)u(m)=u(m−1) is a simple back-shift, (4.7)∑j=0nPn′(j,i)=∑j=0nPn(i,j)=(I−∇)n(q(mn)i)|m=n=q(0)i, and if q is CM, Pn is a stochastic matrix if and only if q(0)=1; else, if q(0)<1,
the matrix Pn is substochastic.
We note that the first column of the matrixPn is Pn(i,0)=q(1)i whereas its first line is Pn(0,j)=δ0,j,
expressing, as required, that the state 0 of (Ar;r∈ℕ) 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 thatp(x)=x corresponding to the simple neutral case. Then q(x)=1−x is completely monotone on (0,1).
With Si,j the second kind Stirling numbers, we get a
lower left triangular stochastic transition matrix (5.1)Pn(i,j)={(nj)∑l=0j(−1)j−l(jl)(ln)i=(n)j⋅n−i⋅Si,j,j≤i,0,else. The diagonal terms (eigenvalues)
are all distinct with Pn(i,i)=(n)i⋅n−i. The matrix Pn is stochastic. Due to triangularity, ancestral
process is a pure death Markov process which may be viewed as a discrete
coalescence tree.From Example3.4 (mutation). Assume that (3.10) holds p(x)=(1−μ2)x+μ1(1−x) where (μ1,μ2) are mutation probabilities. Then, with κ:=1−(μ1+μ2), q(x)=1−μ1−κx is completely monotone on (0,1) if and only μ1≤1−μ2(κ≥0). In this case, Pn is again lower left triangular (a pure death
process). We have (5.2)Pn(i,j)={(nj)∑l=0j(−1)j−l(jl)(μ2+κln)i=:(n)j⋅n−i⋅Si,jμ2(κn),j≤i,0,else, in terms of generalized Stirling
numbers Si,jμ2(κ/n).
We have Pn(i,i)=(n)i(κ/n)i and the spectrum of Pn is real. When μ1>0,
this matrix is substochastic with ∑j=0nPn(i,j)=(1−μ1)i.A particular case deals with one-way mutations(μ1+μ2>0,μ1⋅μ2=0).Ifμ2=0,Pn(i,j)=(1−μ1)i⋅(n)j⋅n−i⋅Si,j,j≤i,=0,
else. Further, ∑j=0nPn(i,j)=(1−μ1)i<1.Ifμ1=0,Pn(i,j)=(n)j⋅n−i⋅Si,jμ2((1−μ2)/n),j≤i,=0,
else. The corresponding matrix Pn is stochastic.From Example3.3 (quadratic). Assume thatp(x)=x(1+c−cx),
as in (3.9). Then q(x)=(1−x)(1−cx) which is completely monotone if and only if c∈[0,1].
The case c=0 is the neutral case, whereas c=1 appears in a 2-sex model of Möhle. In this quadratic case,
since ∇j(q(m/n)i)=0 if j>2i, then Pn(i,j)=0 if j>2i and so Pn is a Hessenberg-like matrix. Note that ∑j=0nPn(i,j)=q(0)i=1.From the selection Example3.1, when (3.6)
holds (5.3)p(x)=(1+s)x1+sx,whereq(x)=1−p(x)=(1−x)/(1+sx) is CM whenever selection parameter s>0.
The induced matrix Pn 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 Example3.2 (selection with dominance). The
corresponding mechanism (3.7) with parameters (s,h) satisfying s>−1 and sh>−1 is CM if and only if s≥sh:=(1−2h)/h2>0 and h∈(0,1/2).
The case h∈(0,1) corresponds to directional selection, where
genotype AA has highest fitness compared to aa's and the heterozygote class Aa has intermediate fitness compared to both
homozygote classes. In this situation, marginal fitness of A exceeds the one of a and selective sweep is expected. When h∈(0,1/2),
allele A is dominant to a,
whereas when h∈(1/2,1),
allele A is recessive to a (a stabilizing effect slowing down the sweep).
Critical value h=1/2 is a case of pure genic balancing selection.Example 5.1.
Consider the mechanismp(x)=xγ for some γ>0.
The function q(x)=1−p(x) 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).
Ifp(x) is not admissible in that q is not CM, it can be that p¯(x):=1−p(1−x) is itself admissible. As observed before, if Nr(m) has transition probabilities given by Bin(n,p(k/n)),p¯(x) arises in the transition probabilities of N¯r(m):=n−Nr(m). Indeed, such transitions are Bin(n,p¯(k/n)) distributed.Ifp(x) is the selection mechanism of Example 3.1,
(3.6), with s∈(−1,0) (not admissible), p¯(x)=(1+s¯)x/(1+s¯x) is itself an admissible selection mechanism
because it has reciprocal selection parameter s¯=−s/(1+s)>0. If p(x) is the mechanism of Example 3.2, namely, (3.7),
with parameters (s,h),
then p¯(x) is itself a selection with dominance mechanism
with reciprocal parameters s¯=−s/(1+s) and h¯=1−h. Assuming that(s≤sh<0,h∈(1/2,1)), p(x) is not admissible whereas p¯(x) is because s¯≥sh¯>0 and h¯∈(0,1/2).
Similarly, when γ∈(0,1),
the mechanism p(x)=1−(1−x)γ is not admissible but, from Example 5.1, p¯(x):=1−p(1−x)=xγ is as follows.


### 5.2. Bias Mechanisms with Mutational Effects

LetpM(x)=(1−μ2)x+μ1(1−x) be the mutational bias mechanism (with κ=1−(μ1+μ2)≥0). Let p(x) be a bias mechanism such that q(x) is CM with p(1)−p(0)=1.
Then, (5.4)p˜M(x)=pM(p(x)) is such that q˜M(x):=1−p˜M(x) is CM. It is, therefore, admissible and adds
mutational effects to the primary mechanism p(x). For example, (5.5)p˜M(x)=μ1+x((1+s)(1−μ2)−μ1)1+sx is a mechanism of selection
combined with mutational effects. We have p˜M(0)=μ1, p˜M(1)=1−μ2. The mechanisms p˜M(x) obtained in this way all share the specificity p˜M(1)−p˜M(0)=:κ<1.Note that, except for the mutational affine mechanism,
it is not true in general that wheneverp1(x) and p2(x) are two admissible bias mechanisms, then p1(p2(x)) is admissible.


### 5.3. Joint Bias Effects and Compound Bias

Letp1(x) and p2(x) be two admissible bias in that q1(x):=1−p1(x) and q2(x):=1−p2(x) are both completely monotone.
Then (5.6)q(x)=q1(x)q2(x)is  CM. Thus, with x1∘x2:=x1+x2−x1x2, the probabilistic product in [0,1](5.7)(p1(x),p2(x))→p(x)=p1(x)∘p2(x). Whenever a WF model is
considered with bias p(x)=p1(x)∘p2(x) obtained from two distinct bias p1(x) and p2(x),
we call it a WF model with joint bias effect.Example 5.3 (joint selection and mutational effects).
Letp1(x)=pM(x) and p2(x)=(1+s)x/(1+sx).
We get (5.8)q(x)=(1−μ1−κx)(1−x)1+sx,p(x)=μ1+x(s+1−μ1+κ)−κx21+sx, with p(0)=μ1, p(1)=1. This mechanism differs from the traditional
mechanism of selection combined with mutational effects.Example 5.4 (joint mutation and neutral effects).
Letp1(x)=(1−μ2)x+μ1(1−x) and p2(x)=x.
We get (5.9)q(x)=(1−x)(1−μ1−κx),p(x)=μ1+x(1−μ1+κ(1−x)), with p(0)=μ1, p(1)=1. When μ1=0 (one-way mutations), we recover the quadratic
mechanism with curvature parameter c=1−μ2. This finding justifies some interest into the
quadratic mechanisms with c≠1.
Withj=1,2,
the reproduction law of each elementary effect is νnj~d Multin(n;πnj), where πm,nj=pj(m/n)−pj((m−1)/n),m=1,…,n. Then, νn~d Multin(n;πn),πm,n=p(m/n)−p(m−1)/n),m=1,…,n,
where πn:=πn1⊙πn2 is easily obtained componentwise
by (5.10)πm,n=πm,n1∑l=1mπl,n2+πm,n2∑l=1mπl,n1,m=1,…,n. We let νn:=νn1⊙νn2~d Multin(n;πn1⊙πn2).
It is the reproduction law of a WF model obtained jointly from the two bias p1(x) and p2(x).Letϕ(x):(0,1)→(0,1) be an absolutely monotone function satisfying ϕ(l)(x)≥0 for all lth derivatives ϕ(l) of ϕ,
for all x∈(0,1).
Such functions are  well known to be probability generating functions (pgfs) of ℕ-valued random variables, say N,
that is to say, ϕ(x)=𝔼[xN]. Clearly, if q is CM on (0,1),
then so is qϕ(x):=ϕ(q(x)). Thus, pϕ(x):=1−ϕ(1−p(x)) is an admissible bias mechanism in that qϕ(x):=1−pϕ(x) is CM. We call it a compound bias.Example 5.5.
The general mechanism with mutational effects is in this class.
Indeed,(5.11)q˜M(x)=1−p˜M(x)=1−pM(p(x))=1−pM(1−q(x)), and so ϕ(x)=1−pM(1−x)=1−(1−μ2)(1−x)−μ1x=μ2+κx which is absolutely monotone as soon as κ=1−(μ1+μ2)≥0.Example 5.6.
Withθ>0,
taking ϕ(x)=e−θ(1−x) or (eθx−1)/(eθ−1),
the pgf of a Poisson (or shifted-Poisson) random variable, pϕ(x)=1−ϕ(1−p(x)) is admissible if p(x) is. Note that if q is of the form qϕ,
where ϕ is the pgf of the Poisson random variable,
then qϕ(x)α is admissible, for all α>0,
a property reminiscent of infinite divisibility for pgfs. Taking ϕ(x)=(1−π)/(1−πx) or x(1−π)/(1−πx), π∈(0,1),
the pgf of a geometric (or shifted-geometric) random variable, pϕ(x)=sp(x)/(1+sp(x)) or (s+1)p(x)/(1+sp(x)) is admissible if p(x) is (with s=π/(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.
Letp(x)=xγ with γ∈(0,1) as in Example 5.1. Then, pϕ(x)=1−qϕ(x),
where qϕ(x)=e−θ(1−q(x))=e−θxγ,θ>0,
is admissible. Note that pϕ(x)~x↓0θxγ.
The reciprocal function p¯ϕ(x)=qϕ(1−x)=e−θ(1−x)γ 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 probabilitiespϕ,
the signification of which in population genetics remaining though to be
pinpointed, in each specific case study.


## 5.1. Elementary Mechanisms

Assume first thatp(x)=x corresponding to the simple neutral case. Then q(x)=1−x is completely monotone on (0,1).
With Si,j the second kind Stirling numbers, we get a
lower left triangular stochastic transition matrix (5.1)Pn(i,j)={(nj)∑l=0j(−1)j−l(jl)(ln)i=(n)j⋅n−i⋅Si,j,j≤i,0,else. The diagonal terms (eigenvalues)
are all distinct with Pn(i,i)=(n)i⋅n−i. The matrix Pn is stochastic. Due to triangularity, ancestral
process is a pure death Markov process which may be viewed as a discrete
coalescence tree.From Example3.4 (mutation). Assume that (3.10) holds p(x)=(1−μ2)x+μ1(1−x) where (μ1,μ2) are mutation probabilities. Then, with κ:=1−(μ1+μ2), q(x)=1−μ1−κx is completely monotone on (0,1) if and only μ1≤1−μ2(κ≥0). In this case, Pn is again lower left triangular (a pure death
process). We have (5.2)Pn(i,j)={(nj)∑l=0j(−1)j−l(jl)(μ2+κln)i=:(n)j⋅n−i⋅Si,jμ2(κn),j≤i,0,else, in terms of generalized Stirling
numbers Si,jμ2(κ/n).
We have Pn(i,i)=(n)i(κ/n)i and the spectrum of Pn is real. When μ1>0,
this matrix is substochastic with ∑j=0nPn(i,j)=(1−μ1)i.A particular case deals with one-way mutations(μ1+μ2>0,μ1⋅μ2=0).Ifμ2=0,Pn(i,j)=(1−μ1)i⋅(n)j⋅n−i⋅Si,j,j≤i,=0,
else. Further, ∑j=0nPn(i,j)=(1−μ1)i<1.Ifμ1=0,Pn(i,j)=(n)j⋅n−i⋅Si,jμ2((1−μ2)/n),j≤i,=0,
else. The corresponding matrix Pn is stochastic.From Example3.3 (quadratic). Assume thatp(x)=x(1+c−cx),
as in (3.9). Then q(x)=(1−x)(1−cx) which is completely monotone if and only if c∈[0,1].
The case c=0 is the neutral case, whereas c=1 appears in a 2-sex model of Möhle. In this quadratic case,
since ∇j(q(m/n)i)=0 if j>2i, then Pn(i,j)=0 if j>2i and so Pn is a Hessenberg-like matrix. Note that ∑j=0nPn(i,j)=q(0)i=1.From the selection Example3.1, when (3.6)
holds (5.3)p(x)=(1+s)x1+sx,whereq(x)=1−p(x)=(1−x)/(1+sx) is CM whenever selection parameter s>0.
The induced matrix Pn 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 Example3.2 (selection with dominance). The
corresponding mechanism (3.7) with parameters (s,h) satisfying s>−1 and sh>−1 is CM if and only if s≥sh:=(1−2h)/h2>0 and h∈(0,1/2).
The case h∈(0,1) corresponds to directional selection, where
genotype AA has highest fitness compared to aa's and the heterozygote class Aa has intermediate fitness compared to both
homozygote classes. In this situation, marginal fitness of A exceeds the one of a and selective sweep is expected. When h∈(0,1/2),
allele A is dominant to a,
whereas when h∈(1/2,1),
allele A is recessive to a (a stabilizing effect slowing down the sweep).
Critical value h=1/2 is a case of pure genic balancing selection.Example 5.1.
Consider the mechanismp(x)=xγ for some γ>0.
The function q(x)=1−p(x) 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).
Ifp(x) is not admissible in that q is not CM, it can be that p¯(x):=1−p(1−x) is itself admissible. As observed before, if Nr(m) has transition probabilities given by Bin(n,p(k/n)),p¯(x) arises in the transition probabilities of N¯r(m):=n−Nr(m). Indeed, such transitions are Bin(n,p¯(k/n)) distributed.Ifp(x) is the selection mechanism of Example 3.1,
(3.6), with s∈(−1,0) (not admissible), p¯(x)=(1+s¯)x/(1+s¯x) is itself an admissible selection mechanism
because it has reciprocal selection parameter s¯=−s/(1+s)>0. If p(x) is the mechanism of Example 3.2, namely, (3.7),
with parameters (s,h),
then p¯(x) is itself a selection with dominance mechanism
with reciprocal parameters s¯=−s/(1+s) and h¯=1−h. Assuming that(s≤sh<0,h∈(1/2,1)), p(x) is not admissible whereas p¯(x) is because s¯≥sh¯>0 and h¯∈(0,1/2).
Similarly, when γ∈(0,1),
the mechanism p(x)=1−(1−x)γ is not admissible but, from Example 5.1, p¯(x):=1−p(1−x)=xγ is as follows.


## 5.2. Bias Mechanisms with Mutational Effects

LetpM(x)=(1−μ2)x+μ1(1−x) be the mutational bias mechanism (with κ=1−(μ1+μ2)≥0). Let p(x) be a bias mechanism such that q(x) is CM with p(1)−p(0)=1.
Then, (5.4)p˜M(x)=pM(p(x)) is such that q˜M(x):=1−p˜M(x) is CM. It is, therefore, admissible and adds
mutational effects to the primary mechanism p(x). For example, (5.5)p˜M(x)=μ1+x((1+s)(1−μ2)−μ1)1+sx is a mechanism of selection
combined with mutational effects. We have p˜M(0)=μ1, p˜M(1)=1−μ2. The mechanisms p˜M(x) obtained in this way all share the specificity p˜M(1)−p˜M(0)=:κ<1.Note that, except for the mutational affine mechanism,
it is not true in general that wheneverp1(x) and p2(x) are two admissible bias mechanisms, then p1(p2(x)) is admissible.


## 5.3. Joint Bias Effects and Compound Bias

Letp1(x) and p2(x) be two admissible bias in that q1(x):=1−p1(x) and q2(x):=1−p2(x) are both completely monotone.
Then (5.6)q(x)=q1(x)q2(x)is  CM. Thus, with x1∘x2:=x1+x2−x1x2, the probabilistic product in [0,1](5.7)(p1(x),p2(x))→p(x)=p1(x)∘p2(x). Whenever a WF model is
considered with bias p(x)=p1(x)∘p2(x) obtained from two distinct bias p1(x) and p2(x),
we call it a WF model with joint bias effect.Example 5.3 (joint selection and mutational effects).
Letp1(x)=pM(x) and p2(x)=(1+s)x/(1+sx).
We get (5.8)q(x)=(1−μ1−κx)(1−x)1+sx,p(x)=μ1+x(s+1−μ1+κ)−κx21+sx, with p(0)=μ1, p(1)=1. This mechanism differs from the traditional
mechanism of selection combined with mutational effects.Example 5.4 (joint mutation and neutral effects).
Letp1(x)=(1−μ2)x+μ1(1−x) and p2(x)=x.
We get (5.9)q(x)=(1−x)(1−μ1−κx),p(x)=μ1+x(1−μ1+κ(1−x)), with p(0)=μ1, p(1)=1. When μ1=0 (one-way mutations), we recover the quadratic
mechanism with curvature parameter c=1−μ2. This finding justifies some interest into the
quadratic mechanisms with c≠1.
Withj=1,2,
the reproduction law of each elementary effect is νnj~d Multin(n;πnj), where πm,nj=pj(m/n)−pj((m−1)/n),m=1,…,n. Then, νn~d Multin(n;πn),πm,n=p(m/n)−p(m−1)/n),m=1,…,n,
where πn:=πn1⊙πn2 is easily obtained componentwise
by (5.10)πm,n=πm,n1∑l=1mπl,n2+πm,n2∑l=1mπl,n1,m=1,…,n. We let νn:=νn1⊙νn2~d Multin(n;πn1⊙πn2).
It is the reproduction law of a WF model obtained jointly from the two bias p1(x) and p2(x).Letϕ(x):(0,1)→(0,1) be an absolutely monotone function satisfying ϕ(l)(x)≥0 for all lth derivatives ϕ(l) of ϕ,
for all x∈(0,1).
Such functions are  well known to be probability generating functions (pgfs) of ℕ-valued random variables, say N,
that is to say, ϕ(x)=𝔼[xN]. Clearly, if q is CM on (0,1),
then so is qϕ(x):=ϕ(q(x)). Thus, pϕ(x):=1−ϕ(1−p(x)) is an admissible bias mechanism in that qϕ(x):=1−pϕ(x) is CM. We call it a compound bias.Example 5.5.
The general mechanism with mutational effects is in this class.
Indeed,(5.11)q˜M(x)=1−p˜M(x)=1−pM(p(x))=1−pM(1−q(x)), and so ϕ(x)=1−pM(1−x)=1−(1−μ2)(1−x)−μ1x=μ2+κx which is absolutely monotone as soon as κ=1−(μ1+μ2)≥0.Example 5.6.
Withθ>0,
taking ϕ(x)=e−θ(1−x) or (eθx−1)/(eθ−1),
the pgf of a Poisson (or shifted-Poisson) random variable, pϕ(x)=1−ϕ(1−p(x)) is admissible if p(x) is. Note that if q is of the form qϕ,
where ϕ is the pgf of the Poisson random variable,
then qϕ(x)α is admissible, for all α>0,
a property reminiscent of infinite divisibility for pgfs. Taking ϕ(x)=(1−π)/(1−πx) or x(1−π)/(1−πx), π∈(0,1),
the pgf of a geometric (or shifted-geometric) random variable, pϕ(x)=sp(x)/(1+sp(x)) or (s+1)p(x)/(1+sp(x)) is admissible if p(x) is (with s=π/(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.
Letp(x)=xγ with γ∈(0,1) as in Example 5.1. Then, pϕ(x)=1−qϕ(x),
where qϕ(x)=e−θ(1−q(x))=e−θxγ,θ>0,
is admissible. Note that pϕ(x)~x↓0θxγ.
The reciprocal function p¯ϕ(x)=qϕ(1−x)=e−θ(1−x)γ 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 probabilitiespϕ,
the signification of which in population genetics remaining though to be
pinpointed, in each specific case study.


## 6. Limit Laws

Consider a WF model(Nr;r∈ℕ) on {0,…,n} with forward transition matrix (6.1)Πn(k,k′)=(nk′)p(kn)k′(1−p(kn))n−k′, with admissible bias p(x). Define (Ar;r∈ℕ) as the dual Markov chain on {0,…,n} with transition probability (6.2)Pn(i,j)=(nj)∑l=0j(−1)j−l(jl)q(1−ln)i. Then, (Nr,Ar;r∈ℕ) are dual with respect to Φn(m,k):=Φn2(m,k)=(n−mk)/(nk),
to wit (6.3)𝔼m[(n−Nrk)(nk)]=𝔼k[(n−mAr)(nAr)]=𝔼k[(n−Arm)(nm)]. We will distinguish two cases.Case 1.
Assume first that(6.4)Nr→dN∞asr↑∞,independently  ofN0=m≥1. Let π∞(i)=ℙ(N∞=i) and π∞=(π∞(0),…,π∞(n))′. The line vector π∞′ is the left eigenvector of Πn associated to the eigenvalue 1:π∞′=π∞′Πn. It is the (unique) invariant probability
measure (stationary distribution) of (Nr;r∈ℕ).
If this stationary distribution exists, then, using
duality formula, necessarily,Ar→0 as r↑∞ with probability ℙk(A∞=0)=:ρ∞(k)<1. The numbers ρ∞(k) are the extinction probabilities of the dual
process started at k. As is well known, ρ∞=(ρ∞(0),…,ρ∞(n))′ is the unique solution to (I−Pn)ρ∞=0 with ρ∞(0)=1.Remark 6.1.
Typical situations where(Nr;r∈ℕ) has an invariant measure is when mutational
effects are present, and more generally when the bias mechanism satisfies p(0)>0 and p(1)<1.
In this situation, the forward stochastic transition matrix Πn has an algebraically simple dominant
eigenvalue 1.
By Perron-Frobenius theorem (6.5)limr↑∞Πnr=1π∞′, where 1′=(1,…,1). The invariant probability measure can be
approximated by subsequent iterates of Πn,
the convergence being exponentially fast, with rate governed by the second
largest eigenvalue. Of course, detailed balance (stating that πkΠn(k,k′)=πk′Πn(k′,k)) does not hold here and the forward chain in
equilibrium is not time-reversible.
In these recurrent cases, the dual ancestral processAr started at k gets extinct with probability ρ∞(k).
The numbers 1−ρ∞(k) are the probabilities that it gets killed
before getting extinct; in other words, 1−ρ∞(k) are the probabilities that Ar first hits an extra coffin state, say {∂}, before hitting {0}.In terms of moments, by the duality formula, we
conclude that(6.6)(nk)−1𝔼[(n−N∞k)]=ρ∞(k)=ℙk(A∞=0),relating k-factorial moments of n−N∞ to the extinction probabilities of Ar given A0=k. We also have (6.7)∑k=0nvk𝔼[(n−N∞k)]=𝔼[(1+v)n−N∞]=∑k=0n(nk)ρ∞(k)vk, and so the probability
generating function of N∞ can be expressed as (u∈[0,1]): (6.8)𝔼[uN∞]=∑k=0n(nk)ρ∞(k)un−k(1−u)k, in terms of the Bernstein-Bézier
polynomial of (ρ∞(n−k);k=0,…,n).Letρ∞=(ρ∞(0),…,ρ∞(n))′. The vector ρ∞ is the right eigenvector of Pn associated to the eigenvalue 1:ρ∞=Pnρ∞. In this case, the matrix Pn is sub-stochastic and the extinction
probability of (Ar;r∈ℕ) given A0=k is less than one. Thanks to duality, we
have (6.9)ΠnΦn=ΦnPn′, where the matrix Φn is symmetric whereas the matrix Φn−1 is symmetric right lower triangular,
with (6.10)Φn(m,k)=(n−mk)(nk)=(n−km)(nm),Φn−1(i,j)=(−1)i+j−n(in−j)(ni)=(−1)i+j−n(jn−i)(nj). Thus, (6.11)π∞′ΠnΦn=π∞′Φn=π∞′ΦnPn′, showing that ρ∞ and π∞ are related through (6.12)ρ∞=Φnπ∞orπ∞=Φn−1ρ∞.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 Example3.4. In this case, Nr→dN∞ as r↑∞, regardless of N0=m and (Nr;r∈ℕ) has an invariant measure which is difficult to
compute. Looking at the backward process, the matrix Pn is substochastic (if μ1>0) and lower left triangular. Due to
triangularity, the right eigenvector ρ∞ of Pn can easily be computed explicitly in terms of (Pn(i,j);j≤i),i=0,…,n. We, therefore, get the following alternating
expression for the invariant measure: (6.13)π∞(i)=(ni)∑j=0i(−1)i−j(ij)ρ∞(n−j). Concerning moments, for
instance, we have ρ∞(1)=μ2/(μ1+μ2) so that 𝔼[N∞]=nμ1/(μ1+μ2);
from (5.2) we also have  (6.14)ρ∞(2)=μ2[nμ2(1+κ)+κ2][(1−κ)(n−(n−1)κ2)]=1n(n−1)(n(n−1)−n(2n−1)μ1μ1+μ2+𝔼[N∞2]), allowing to compute 𝔼[N∞2] and then the variance of N∞. We get (6.15)σ2(N∞)=n2μ1μ2(μ1+μ2)2(2n(μ1+μ2)+1)+o(n)~n↑∞μ1μ22(μ1+μ2)3n, suggesting (when μ1μ2>0) a central limit theorem for N∞ as n grows large: (6.16)1n(N∞−nμ1μ1+μ2)→n↑∞d𝒩(0,μ1μ22(μ1+μ2)3).Case 2.
Conversely, assume now that givenN0=m(6.17)Nr→d0asr↑∞,with  probabilityℙm(N∞=0)=:ρ∞(m), so that boundaries {0,n} are absorbing. Then, the ancestral process (Ar;r∈ℕ) possesses an invariant distribution, in
that (6.18)Ar→dA∞asr↑∞,independently  ofA0=k∈[n]. In terms of moments, the duality
formula means that (6.19)(nm)−1𝔼[(n−A∞m)]=ρ∞(m)=ℙm(N∞=0),relating m-factorial moments of n−A∞ to the extinction probabilities of Nr given N0=m. Stated differently, the probability generating
function of A∞ is (u∈[0,1]) (6.20)𝔼[uA∞]=∑m=0n(nm)ρ∞(m)un−m(1−u)m. Let π∞(i)=ℙ(A∞=i),
with π∞′=π∞′Pn. Then, using duality, ρ∞ is the right eigenvector of Πn associated to the eigenvalue 1:ρ∞=Πnρ∞. Thus, ρ∞ and π∞ are related through (6.21)ρ∞=Φnπ∞orπ∞=Φn−1ρ∞.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,n} are absorbing to (Nr;r∈ℕ) occur when p(0)=0 and p(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,ρ∞(m)=1−m/n. Thus, π∞(i)=(ni)∑j=0i(−1)i−j(ij)(j/n)=δi,1 and Ar→d1 as r↑∞, the degenerate state reached when the most
recent common ancestor (MRCA) is attained,(ii)
nondegenerate solutions ofA∞ are obtained when considering bias mechanisms
with p(0)=0 and p(1)=1,(iii)
consider any biased WF model withp(0)=0 and p(1)=1 for which p(x)~x↑0λx, λ>1. Then, due to large sample asymptotic
independence,(6.22)νn→dξ∞,where ξ∞ is an iid sequence with ξ1~d Poisson (λ) (as it can easily be checked by the Poisson
limit to the binomial distribution). In this case, the limiting extinction
probability of (Nr;r∈ℕ) given N0=m is limn↑∞ρ∞(m)=ρm,m=1,2,…,
where 0<ρ<1 is the smallest solution to the fixed point
equation(6.23)x=e−λ(1−x),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(6.24)n(ρm−ρ∞(m))=ρm(1−ρ1+λρm2+λ(1−ρ)ρ1−(λρ)2m), showing that the convergence of ρ∞(m) to ρm is of order n−1.
As a result, we get the asymptotic normality (6.25)1n(A∞−n(1−ρ))→n↑∞d𝒩(0,ρ(1−ρ)1+λρ). Intuitively, (1/n)𝔼[n−A∞]=ρ∞(1)=ℙ1(N∞=0)→ρ,
showing that 𝔼[A∞]~n↑∞n(1−ρ) and (6.26)1n(n−1)𝔼[(n−A∞)(n−1−A∞)]=ρ∞(2)=ℙ2(N∞=0)→ρ2, showing that σ2(A∞)~n↑∞nρ(1−ρ)/(1+λρ).
For quadratic Example3.3, p(x)=x(1+c−cx),
with c∈[0,1], λ=1+c>1 as soon as c>0.
When c∈(0,1],
we thus always have asymptotic normality. For Example 3.1, with selection, p(x)=(1+s)x/(1+sx), with s>−1:p(x)~x↑0(1+s)x and so λ=1+s.
We have asymptotic normality only when s>0, that
is, when the fitness is advantageous (corresponding as required to complete
monotonicity of corresponding q=1−p).Note that this asymptotic behavior does not hold for
the Lipshitz continuous admissible mechanismp(x)=xγ of Example 5.1 (or more generally for
mechanisms satisfying p(x)~x↓0θxγ, θ>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.



---

*Source: 714701-2009-02-15.xml*