# A Wright-Fisher graph model and the impact of directional selection on
genetic variation

Ingemar Kaja,∗, Carina F. Mugalb,c, Rebekka Müllera

## Abstract

We introduce a multi-allele Wright-Fisher model with non-recurrent, reversible
mutation and directional selection. In this setting, the allele frequencies at
a single locus track the path of a hybrid jump-diffusion process with state
space given by the vertex and edge set of a graph. Vertices represent
monomorphic population states and edge-positions mark the biallelic
proportions of ancestral and derived alleles during polymorphic segments. We
derive the stationary distribution in mutation-selection-drift equilibrium and
obtain the expected allele frequency spectrum under large population size
scaling. For the extended model with multiple independent loci we derive
rigorous upper bounds for a wide class of associated measures of genetic
variation. Within this framework we present mathematically precise arguments
to conclude that the presence of directional selection reduces the magnitude
of genetic variation, as constrained by the bounds for neutral evolution.

### Keywords

Wright-Fisher jump-diffusion process, directional selection, mutation bias,
genetic diversity, effective mutation rate, theoretical population genetics

aDepartment of Mathematics, Uppsala University, Uppsala, Sweden
bDepartment of Ecology and Genetics, Uppsala University, Uppsala, Sweden
cLaboratory of Biometry and Evolutionary Biology, University of Lyon 1, UMR
CNRS 5558, Villeurbanne, France
∗Corresponding author<EMAIL_ADDRESS>
##  1 Introduction

The degree of genetic variation within a population is determined by the
interrelations of evolutionary processes such as mutation, genetic drift and
natural selection. Mutation, the fundamental source of genetic variation, is
frequently modeled as a non-recurrent mutation mechanism that initializes the
segregation of an allele in the population but otherwise does not influence
the population frequency of the allele (Kimura 1969; Sawyer and Hartl 1992;
McVean and Charlesworth 1999). Genetic drift and natural selection, on the
other hand, control the time span over which mutations segregate in a
population until eventually reaching fixation or extinction. While genetic
drift ultimately acts to eliminate genetic variation, different selection
mechanisms can either prolong or shorten the time to fixation or extinction.
Directional selection, where one of the alleles in a given pair of allelic
types has a selective advantage over the other, is commonly viewed as a force
to reduce the level of genetic variation. However, as pointed out in Novak and
Barton (2017), ”rigorous arguments for this idea are scarce”.

Population genetics modeling for the purpose of analyzing genetic variation
under the combined influence of different evolutionary processes typically
builds on some version of Wright-Fisher models (Fisher 1930; Wright 1931,
1938) or Moran type models (Moran 1958). As pioneered by Kimura (1964),
diffusion approximation techniques under scaling of evolutionary time and
large population size are instrumental and helped advance the understanding of
the distribution of inherited allele frequencies, both dynamically and under
steady-state, see e.g. Durrett (2008); Etheridge (2011). The original Wright-
Fisher model with mutation and selection (Wright 1931) concerns the allele
frequency distribution for two allelic types and distinguishes recurrent and
non-recurrent mutation. With regards to approximating allele frequencies with
trajectories of diffusion processes, non-recurrent mutation occur on the
boundary of the state space and provide the renewal of polymorphic segments.
Recurrent mutation is ongoing and appears through linear drift terms in the
diffusion generator (Etheridge 2011).

In this work we apply diffusion approximation methods to study a multi-allele
and multi-locus model with non-recurrent, reversible mutation and directional
selection in an isolated population assuming independence among loci. It is a
consequence of the assumption of the non-recurrent mutation mechanism that the
overall mutational input is small enough to prevent additional allelic types
at a locus that is already polymorphic. This is consistent with observations
in empirical data where multi-allelic single nucleotide variation is typically
rare (Cao et al. 2015; Phillips et al. 2015). In addition, mutation is
reversible since we work with a fixed, finite number of allelic types and all
mutation events involving a given pair of alleles may take place in both
directions. Within this framework our objective is essentially to show that
the presence of directional selection reduces the magnitude of genetic
variation, as constrained by the bounds for neutral evolution. To this aim we
derive stationary distributions over monomorphic and polymorphic states in
mutation-selection-drift equilibrium. Closed form expressions for the expected
allele frequency spectrum are obtained asymptotically under large population
size scaling and rigorous upper bounds are derived for a wide class of
associated measures of genetic variation under the influence of directional
selection.

To put our approach in context, the extension from studying the evolutionary
dynamics of genetic loci with two types to general multi-allele Wright-Fisher
models with a fixed number $K\geq 2$ possible allelic states for each genetic
locus, can be traced back to Wright (1949). For the case of recurrent mutation
mechanisms such $K$-allele models have been developed in much detail. The
state space for single locus frequencies is now (a subset of) the $K$-simplex,
which presents considerable challenges in extracting useful probabilistic
information. For a brief history of K-allele Wright-Fisher models with
recurrent mutation, we refer to Ferguson and Buzbas (2018), and for some of
the mathematical results to Etheridge (2011, Ch. 4 and 5). A recent
approximation approach to multi-allele models with recurrent mutation (Burden
and Tang 2016; Ferguson and Buzbas 2018) starts from the presumption that
mutation events are rare on the time scale of evolution relevant for the
diffusion approximation. Then, with sufficiently small mutation rates, the
allele frequencies will be mostly concentrated either on the vertices of the
$K$-simplex, or on the edges connecting a pair of mutating alleles and with
only a small fraction of probability mass remaining on simplex domains that
allow three or more alleles existing simultaneously.

In contrast, our approach towards modeling the multi-allelic case relies on
non-recurrent and reversible mutation in between a fixed number of genetic
types. A jump-diffusion process is introduced, biallelic by construction, with
state space consisting only of the vertices and edges of the graph subset of
the simplex. Some key features of the jump-diffusion process are already
implemented in Mugal et al. (2014) and Kaj and Mugal (2016) for the simpler
setting of arbitrary ancestral-and-derived alleles. In the graph model, the
vertices correspond to the presence of a specific fixed type (or allele) in a
genetic locus and the edges between two types represent continuous polymorphic
states. Similar boundary mutation multi-allele models have been discussed in
the context of synonymous codon usage (Zeng 2010) or so-called polymorphism-
aware phylogenetic models (De Maio et al. 2013; Borges et al. 2019), with a
focus on methodological development for statistical inference from genomic
data.

##  2 A Wright-Fisher graph model

A continuous time Markov process with state space given by a connected,
directed graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with vertices
$\mathcal{V}$ and unit length edges $\mathcal{E}$ captures the random change
of allelic types at a single locus. Such a locus can be of abstract nature
(mutant versus wild type) or specific, for instance consisting of $l$
nucleotides in the genome, where $l=1$ corresponds to a single site on the
genome or $l=3$ to nucleotide triplets, such as protein-coding codons. The
finite vertex set represents the available allelic types at the locus, while
the family of edges allows for keeping track of the possible mutations among
types and any polymorphic state. Each edge is a directed interval of length
one, starting in a vertex $u\in\mathcal{V}$ and leading to another vertex
$v\in\mathcal{V}$, such that the position on the edge records the relative
frequency of type $v$ as a mutant derived from ancestral type $u$. The
relevant graph-valued model is a hybrid jump and diffusion process with
compact state space, in which the open edges form a continuous interior and
the vertices are discrete boundary points. Mutation events occur only on the
boundary. Each mutation is succeeded by a polymorphic phase of two alleles co-
existing in the population, upon which the Wright-Fisher diffusion determines
the frequency and subsequent extinction and fixation probabilities of the
mutant. The graph in Fig. 1 illustrates an example state space on which the
graph-valued process moves.

Formally, we consider a Markov process $X=(X_{t})_{t\geq 0}$, with values in
the compact state space $D$ formed by the graph equipped with intervals
$[0,1)$ associated with each edge in $\mathcal{E}$, directed from $0$ to $1$,
and encoded by a triplet
$X_{t}=(U_{t},V_{t},Y_{t})\in\mathcal{V}\times\mathcal{V}\times[0,1)$. The
boundary of the state space consists of the point set $\partial
D=\\{(u,u,0)\colon u\in\mathcal{V}\\}$, where the boundary state $(u,u,0)$
represents a monomorphic locus at which the entire population has the same
allelic type $u\in\mathcal{V}$. A state $(u,v,y)$ within the interior of the
state space

$D^{\circ}=\\{(u,v,y)\colon u\in\mathcal{V},v\in\mathcal{V},u\not=v,0<y<1\\}$

lies on the directed edge leading from vertex state $(u,u,0)$ to vertex state
$(v,v,0)$. It arises when a mutation from $u$ to $v$ occurred and brought
mutants of type $v$ to be present in the population at relative frequency in
the infinitesimal interval $(y,y+dy)$. Consequently, the interior state
$(v,u,y)$, located on the complementary edge directed in the opposite
direction from $v$ to $u$, assigns relative frequency $y$ to mutant type $u$
derived from an ancestral $v$. Finally, the closure $D=D^{\circ}\cup\partial
D$ is reached along the limits

$(u,v,y)\to\left\\{\begin{array}[]{cc}(u,u,0)\in\partial D,&y\to 0\\\
(v,v,0)\in\partial D,&y\to 1\end{array}\right.,\quad(u,v,y)\in D^{\circ}.$

Figure 1: Panel A: A connected graph as state space for the graph-valued
process. Between each pair of connected vertices there are in fact two
directed edges. Panel B: Zoom-in of directed edges between types $u$ and $v$,
equipped with unit intervals. A mutation on the boundary from type $u$ to type
$v$ occurs with intensity $\lambda_{uv}$. The initial frequency of the mutant
$v$ is a fixed value $x\in(0,1)$ which is also the initial value for a Wright-
Fisher diffusion with selection coefficient $\gamma_{uv}$ that is started by
the mutation. The process either goes to fixation in type $v$ or to extinction
in $u$.

###  2.1 Reversible boundary mutation

The mutation mechanism of the process $X$ is specified by a fixed entry point
$x\in(0,1)$ and a family of nonnegative mutation rate parameters
$\\{\lambda_{uv}\colon v\not=u\\}$. Here, $x$ represents the fraction of a
population that is affected by a single mutation. We thus define
$\lambda_{uv}/x$ as the population mutation intensity from $u$ to $v$ per
evolutionary time unit, where the evolutionary time unit corresponds to
”$x^{-1}$ generations”. This means $\lambda_{uv}$ can be thought of as the
population mutation rate ”per generation”, which commonly represents the
macroscopic mutation rate. The motivation for this particular choice becomes
clear in the subsequent Section 3 when we introduce a parameter $N$ for
population size, take $x=1/N$, and run mutations at a rate of speed
$\lambda_{uv}/x=N\lambda_{uv}$.

The graph edge set consists of those edges $e_{uv}$ between types for which
the mutation intensity is positive, $\mathcal{E}=\\{e_{uv}\colon
u,v\in\mathcal{V},\,\lambda_{uv}>0\\}$. We postulate that every type is
essential, that all mutations are reversible, and that the graph $\mathcal{G}$
is irreducible, by assuming that the mutation rates satisfy the conditions

  * i)

$\lambda_{u}:=\sum_{v:v\not=u}\lambda_{uv}>0$ for each $u\in\mathcal{V}$,

  * ii)

$\lambda_{uv}>0\iff\lambda_{vu}>0$,

  * iii)

for each pair of allelic types, $u,v\in\mathcal{V}$, there is a sequence of
edge mutations $u\to w_{1}\to\dots\to w_{n}\to v$, such that
$\lambda_{u,w_{1}}\cdot\ldots\cdot\lambda_{w_{n},v}>0$.

The first assumption guarantees that no type is a mutational trap, as the
total intensity $\lambda_{u}$ of a mutation from $u$ to some other type is
strictly positive. The second condition ensures that any mutation from one
type to another may also occur in the reverse direction, i.e. mutation is
reversible. The third condition entails the assumption that every type may be
reached from every other type by a chain of non-recurrent mutation events all
occurring on the boundary.

The hybrid jump and diffusion mechanism of $X$ is such that, starting in a
boundary point $z=(u,u,0)\in\partial D$, the process holds during an
exponential time with intensity $\lambda_{u}/x$. It then jumps to an interior
point $z^{\prime}=(u,v,x)\in D^{\circ}$ governed by jump rates $\lambda_{uv}$,
which represents mutant type $v$ entering the population at (continuous)
fraction $x$. Assuming that a jump from $(u,u,0)$ to $(u,v,x)$ occurs at time
$r$, the interior trajectory of $X$ is a continuous path

$X_{t}=(u,v,Y_{t}^{r}),\quad r\leq t<\tau,$

where $Y^{r}_{t}$, $t\geq r$, is a diffusion in $D$ starting from
$Y^{r}_{r}=x$ such that the path $\xi_{s}^{r}=Y^{r}_{r+s}$, $s\geq 0$, is a
Wright-Fisher process with selection, which solves

$d\xi_{s}=\gamma_{uv}\,\xi_{s}(1-\xi_{s})\,ds+\sqrt{\xi_{s}(1-\xi_{s})}\,dB^{uv}_{s},\quad\xi_{0}=x.$
(1)

We denote by $\tau$ the first exit time of the interior state space
$D^{\circ}$, i.e. $X_{\tau}=(u,u,0)$ if $\xi_{s}$ is absorbed in 0 and
$X_{\tau}=(v,v,0)$ if absorbed in 1. The former case is extinction and the
latter case is fixation of the allelic type $v$. For edges
$e_{uv}\in\mathcal{E}$ the parameter $\gamma_{uv}$ denotes the selection
coefficient for mutations from $u$ to $v$. Moreover,
$\\{B^{uv},e_{uv}\in\mathcal{E}\\}$ are i.i.d. standard Brownian motions.
Figure 1B depicts the hybrid jump and diffusion setup.

###  2.2 Directional selection

To express genic selection in the model we let the family of selection
coefficients $\gamma_{uv}$ assigned to the edges $e_{uv}$ in the graph satisfy
the anti-symmetric condition

$\gamma_{uv}=-\gamma_{vu},\quad e_{uv}\in\mathcal{E}.$ (2)

The central instance is directional selection based on a static fitness
landscape, where each allelic type is assigned a (time-independent) fitness
level $F_{u}$, $u\in\mathcal{V}$, and each polymorphic pair of alleles $(u,v)$
has relative selection coefficient

$\gamma_{uv}=F_{v}-F_{u},\quad e_{uv}\in\mathcal{E},$ (3)

i.e. $\gamma_{uv}=-\gamma_{vu}$. Consequently, in each pair the type with the
highest fitness has the selective advantage at the expense and equal
disadvantage of the other type. Another relevant example of anti-symmetric
selection coefficients is the preferential fixation of strong (C and G) over
weak (A and T) nucleotides due to the process of GC-biased gene conversion
(gBGC) (Duret and Galtier 2009; Mugal et al. 2015), which analytically is
equivalent to directional selection (Nagylaki 1983).

We note that under assumption (2) the distribution of the process $(X_{t})$
simplifies on each pair of edges $e_{uv}$ and $e_{vu}$ through the equality in
distribution

$(u,v,Y_{t}^{r})\stackrel{{\scriptstyle d}}{{=}}(v,u,1-Y_{t}^{r}),\quad r\leq
t<\tau.$

###  2.3 Properties of the Wright-Fisher graph process

#### Green’s function

The graph process $(X_{t})_{t\geq 0}$ restricted to a particular edge
$e_{uv}\in\mathcal{E}$, is a classical Wright-Fisher diffusion, described in
Eq. 1, with selection coefficient $\gamma_{uv}$. We write
$\mathbb{P}^{\gamma}_{x}$ for the probability measure and
$\mathbb{E}^{\gamma}_{x}$ for the expectation of the process starting at $x$,
and select the associated scale function $S_{\gamma}(x)$ and speed function
$m_{\gamma}(x)$ as

$S_{0}(x)=x,\quad S_{\gamma}(x)=\frac{1}{2\gamma}(1-e^{-2\gamma
x}),\quad\gamma\not=0,\qquad m_{\gamma}(x)=\frac{e^{2\gamma x}}{x(1-x)}.$

Since $m_{\gamma}$ is not integrable near $0$ or $1$, both points $\\{0,1\\}$
are exit boundary points, therefore accessible from the interior of the state
space. The diffusion can reach either of these boundaries but will stay at the
point hit first. The time $\tau_{0}$ required to reach $0$ is the extinction
time, the time $\tau_{1}$ to reach $1$ the fixation time, and
$\tau=\min(\tau_{0},\tau_{1})$ is the exit time of the interior interval
$(0,1)$. The corresponding fixation probability $q_{\gamma}(x)$ equals
$q_{\gamma}(x)=S_{\gamma}(x)/S_{\gamma}(1)$ (Kimura 1962), hence

$q_{\gamma}(x)=\mathbb{P}^{\gamma}_{x}(\tau_{1}<\tau_{0})=\frac{1-e^{-2\gamma
x}}{1-e^{-2\gamma}},\quad\gamma\not=0,\qquad q_{0}(x)=x.$

The occupation time functional

$\mathbb{E}^{\gamma}_{x}\Big{[}\int_{0}^{\tau}g(\xi_{s})\,ds\Big{]}=\int_{0}^{1}G_{\gamma}(x,y)g(y)\,dy,$
(4)

is obtained from the Green function $G_{\gamma}(x,y)$ defined by

$G_{\gamma}(x,y)=\left\\{\begin{array}[]{ll}2q_{\gamma}(x)(S_{\gamma}(1)-S_{\gamma}(y))m_{\gamma}(y),&0\leq
x\leq y\leq 1,\\\\[5.69054pt]
2(1-q_{\gamma}(x))(S_{\gamma}(y)-S_{\gamma}(0))m_{\gamma}(y),&0\leq y\leq
x\leq 1.\end{array}\right.$ (5)

Assuming we start from $z=(u,v,x)\in D^{\circ}$, the possible transitions from
$z$ to $z^{\prime}$ before exiting $D^{\circ}$ are those such that
$z^{\prime}=(u,v,y)$ for some $y$, $0<y<1$. For such a pair, Green’s function
is $G(z,z^{\prime})=G_{\gamma_{uv}}(x,y)$ and governs the transition of the
process from $z$ towards $z^{\prime}$. For background on mathematical
population genetics and more detailed properties of the Wright-Fisher
diffusion process with selection, we refer the reader to e.g. Maruyama (1977);
Karlin and Taylor (1981); Ewens (2004); Etheridge (2011).

#### Invariant measure on the boundary

By replacing the polymorphic excursions of $X$ with instantaneous jumps, we
obtain an embedded continuous time Markov chain. Indeed, starting in
$u\in\mathcal{V}$ the embedded chain holds during an exponential time with
rate $\lambda_{u}/x$, then with probability $\lambda_{uv}/\lambda_{u}$ picks
type $v$ and with probability $q_{\gamma_{uv}}(x)$ jumps to the new type $v$.
Taken together, the transition rate of the embedded chain from $u$ to $v$ is
$h_{uv}/x$, where $h_{uv}=\lambda_{uv}q_{\gamma_{uv}}(x)$. The invariant
measure on the boundary is a probability measure
$\eta^{x}=\\{\eta^{x}(z)=\eta_{u}^{x},z=(u,u,0)\in\partial D\\}$, which, since
the state space is finite, is a steady-state for the embedded chain. For a
pair of vertices $u$, $v$, the invariant distribution must satisfy

$\eta_{u}^{x}\lambda_{uv}q_{\gamma_{uv}}(x)=\eta_{v}^{x}\lambda_{vu}(1-q_{\gamma_{uv}}(1-x)),$

since fixation of a mutant $v$ starting from frequency $x$ is the same as
extinction of the ancestral type $u$ starting from frequency $1-x$. It is
straightforward to verify that assumption (2) for the selection coefficients
in the model now yields the relationship

$\eta_{u}^{x}h_{uv}=\eta_{u}^{x}\lambda_{uv}q_{\gamma_{uv}}(x)=\eta_{v}^{x}\lambda_{vu}q_{\gamma_{vu}}(x)=\eta_{v}^{x}h_{vu},\quad
e_{uv}\in\mathcal{E}.$ (6)

This shows that the embedded chain satisfies a detailed balance equation
across each pair of edges of the graph. The ”current of probability” from $u$
to $v$, given by $\eta_{u}^{x}h_{uv}$, equals the corresponding flow
$\eta_{v}^{x}h_{vu}$ from $v$ to $u$.

Hence there exists a unique invariant and reversible measure
$(\eta_{u}^{x})_{u\in\mathcal{V}}$ on the boundary, typically associated with
time-reversibility. The matrix $\mathcal{H}_{x}=(h_{uv}(x))/x$ with diagonal
elements $h_{uu}(x)=-\sum_{v}\lambda_{uv}q_{\gamma_{uv}}(x)$ is the
infinitesimal generator of the continuous time Markov chain
$\widehat{U}_{t}=U_{t}|Y_{t}=0$, for which the conditioning formalizes the
notion of instantaneous jumps and hence suppressed polymorphisms.

We are now in the position to connect the probability weights on the boundary
given by $\eta^{x}$ with the occupation measure in the interior of the state
space as provided by the Green function. If the process starts from the
boundary according to the reversible measure $\eta^{x}$, i.e. the initial
distribution of $X_{0}$ is $\eta^{x}$, then the dynamics of the first jump is
governed by

$\nu^{x}=\sum_{u,v\in\mathcal{V}}\eta_{u}^{x}\lambda_{uv}\delta_{(u,v,x)}.$
(7)

In words, the intensity of the first jump is determined by the accumulation of
jump intensities over all vertices and its outgoing edges, weighted by the
probability to be in a particular vertex. The subsequent relative position on
a particular edge is then given by the fixed entry point $x$. Under $\nu^{x}$,
the relevant Green function contribution to $z^{\prime}=(u,v,y)$ is

$G(\nu^{x},z^{\prime})=\eta_{u}^{x}\lambda_{uv}G_{\gamma_{uv}}(x,y),$ (8)

which also justifies writing

$G(\nu^{x},dz^{\prime})=\eta_{u}^{x}\lambda_{uv}G_{\gamma_{uv}}(x,y)\,dy.$ (9)

#### The generator

We consider real-valued functions $f$ defined on $D$, writing $f(z)=f_{uu}(0)$
for $z=(u,u,0)\in\partial D$ and $f(z)=f_{uv}(y)$ for $z=(u,v,y)\in
D^{\circ}$, and let $f_{uv}^{\prime}(y)$ and $f_{uv}^{\prime\prime}(y)$ denote
first and second order derivatives with respect to $y$ defined in the interior
$D^{\circ}$ of $D$. The infinitesimal generator of the Markov process
$(X_{t})_{t\geq 0}$ is the operator $\mathcal{L}$ which acts on a suitable
domain $\mathcal{D}$ of functions $f\colon D\to\mathbb{R}$, twice continuously
differentiable in the interior $D^{\circ}$, by

$\mathcal{L}f(z)=\sum_{v\in\mathcal{V}}\lambda_{uv}(f_{uv}(x)-f_{uu}(0)),\quad
z=(u,u,0)\in\partial D,$

and

$\mathcal{L}f(z)=\gamma_{uv}y(1-y)\,f_{uv}^{\prime}(y)+\frac{1}{2}y(1-y)\,f_{uv}^{\prime\prime}(y),\quad
z=(u,v,y)\in D^{\circ}.$

Let $\mathcal{F}$ denote the class of real-valued bounded functions on $D$
such that $f\in\mathcal{D}$, and such that $f_{uv}(y)\to f_{uu}(0)$ as $y\to
0$ and $f_{uv}(y)\to f_{vv}(0)$ as $y\to 1$, $u,v\in\mathcal{V}$.

###  2.4 Stationary distribution

In order to determine the equilibrium behavior of the Wright-Fisher graph
process we will construct a stationary distribution, i.e. a distribution which
is preserved under the time-dynamics of the model and hence represents the
typical probability weight assigned to the various parts of the graph in
steady-state. It can be shown in addition that the graph process satisfies
exponential ergodicity and that the stationary distribution is the unique
limit distribution. In this extended sense, the stationary distribution
measures how likely it is that the process visits a certain position on the
graph after allowing sufficient time to reach an equilibrium. However, the
proof of exponential ergodicity is extensive and outside the scope of the work
at hand.

For each edge $e_{uv}\in\mathcal{E}$ let $D_{e_{uv}}=\\{(u,v,y)\colon
u,v\in\mathcal{V},0\leq y<1\\}$ be the subset of $D$ which consists of the
monomorphic state $(u,u,0)$ and all polymorphic states $(u,v,y),0\leq y<1,$
between $u$ and $v$. Then $\cup_{e_{uv}\in\mathcal{E}}D_{e_{uv}}=D$ and the
intersection of two edge sets contains any shared vertex. A measure $\mu$ on
$D$ is stationary for $(X_{t})_{t\geq 0}$, by definition, if the balance
equations

$\int_{D}\mathcal{L}f(z)\,\mu(dz)=\sum_{e_{uv}\in\mathcal{E}}\int_{D_{e_{uv}}}\mathcal{L}f(z)\,\mu(dz)=0,\quad
f\in\mathcal{F},$ (10)

hold. We say that the measure $\mu$ is edge-reversible for the Wright-Fisher
graph process $(X_{t})_{t\geq 0}$, if the detailed balance edge equations

$\int_{D_{e_{uv}}}\mathcal{L}f(z)\,\mu(dz)+\int_{D_{e_{vu}}}\mathcal{L}f(z)\,\mu(dz)=0,\quad
f\in\mathcal{F},$

hold for every pair of edges $e_{uv},e_{vu}\in\mathcal{E}$. By summing this
relation over all pairs $u$ and $v$, linked by the two edges $e_{uv}$ and
$e_{vu}$, we recover Eq. 10. Thus, an edge-reversible measure $\mu$ yields a
stationary distribution of the Wright-Fisher graph process.

######  Theorem 1.

There exists an edge-reversible measure $\mu$ for $(X_{t})_{t\geq 0}$ on $D$,
which is given by

$\displaystyle\mu(z)$
$\displaystyle=\frac{\eta^{x}(z)}{1+\int_{D^{\circ}}G(\nu^{x},dz^{\prime})},\quad
z\in\partial D,$ $\displaystyle\mu(dz)$
$\displaystyle=\frac{G(\nu^{x},dz)}{1+\int_{D^{\circ}}G(\nu^{x},dz^{\prime})},\quad
z\in D^{\circ},$

where $\eta^{x}(z)$, $z\in\partial D$, is the unique boundary measure defined
by the detailed balance equations in Eq. 6, $\nu^{x}$ is be the averaged jump
measure in Eq. 7, both dependent on $x$, and $G(\nu^{x},dz^{\prime})$ is the
measure on $D^{\circ}$ introduced in Eq. 9.

###### Proof.

Put $\Omega=1+\int_{D^{0}}G(\nu^{x},dz^{\prime})$. To verify that $\mu$ is
edge-reversible we need to establish for each pair of edges
$e_{uv},e_{uv}\in\mathcal{E}$ the identity

$\begin{split}&\Omega\Big{(}\int_{D_{e_{uv}}}\mathcal{L}f(z)\,\mu(dz)+\int_{D_{e_{vu}}}\mathcal{L}f(z)\,\mu(dz)\Big{)}\\\
&=\eta_{u}^{x}\lambda_{uv}(f_{uv}(x)-f_{uu}(0))+\eta_{u}^{x}\lambda_{uv}\int_{0}^{1}\mathcal{L}f_{uv}(y)\,G_{\gamma_{uv}}(x,dy)\\\
&\quad+\eta_{v}^{x}\lambda_{vu}(f_{vu}(x)-f_{vv}(0))+\eta_{v}^{x}\lambda_{vu}\int_{0}^{1}\mathcal{L}f_{vu}(y)\,G_{\gamma_{vu}}(x,dy)=0,\end{split}$
(11)

$f\in\mathcal{F}$. Here, using Eq. 5,

$\int_{0}^{1}\mathcal{L}f_{uv}(y)\,G_{\gamma_{uv}}(x,dy)=q_{\gamma_{uv}}(x)A_{uv}(x)+(1-q_{\gamma_{uv}}(x))B_{uv}(x)$

with

$A_{uv}(x)=\int_{x}^{1}\Big{(}f_{uv}^{\prime}(y)+\frac{1}{2\gamma_{uv}}f_{uv}^{\prime\prime}(y)\Big{)}(1-e^{-2\gamma_{uv}(1-y)})\,dy$

and

$B_{uv}(x)=\int_{0}^{x}\Big{(}f_{uv}^{\prime}(y)+\frac{1}{2\gamma_{uv}}f_{uv}^{\prime\prime}(y)\Big{)}(e^{2\gamma_{uv}y}-1)\,dy.$

Partial integration twice in each of $A_{uv}(x)$ and $B_{uv}(x)$ yield,
noticing that $f_{uv}(y)\to f_{vv}(0)$ as $y\to 1$,

$A_{uv}(x)=f_{vv}(0)-f_{uv}(x)-\frac{1}{2\gamma_{uv}}f_{uv}^{\prime}(x)(1-e^{-2\gamma_{uv}(1-x)})$

and

$B_{uv}(x)=f_{uu}(0)-f_{uv}(x)+\frac{1}{2\gamma_{uv}}f_{uv}^{\prime}(x)(e^{2\gamma_{uv}x}-1).$

Since

$q_{\gamma_{uv}}(x)\frac{1-e^{-2\gamma_{uv}(1-x)}}{2\gamma_{uv}}-(1-q_{\gamma_{uv}}(x))\frac{e^{2\gamma_{uv}x}-1}{2\gamma_{uv}}=0,$

it follows that

$\displaystyle q_{\gamma_{uv}}(x)A_{uv}(x)+(1-q_{\gamma_{uv}}(x))B_{uv}(x)$
$\displaystyle\quad=q_{\gamma_{uv}}(x)(f_{vv}(0)-f_{uv}(x))+(1-q_{\gamma_{uv}}(x))(f_{uu}(0)-f_{uv}(x)),$

which is

$\displaystyle\int_{0}^{1}\mathcal{L}f_{uv}(y)\,G_{\gamma_{uv}}(x,dy)=-(f_{uv}(x)-f_{uu}(0))+q_{\gamma_{uv}}(x)(f_{vv}(0)-f_{uu}(0)).$

Similarly, by symmetry,

$\displaystyle\int_{0}^{1}\mathcal{L}f_{vu}(y)\,G_{\gamma_{vu}}(x,dy)=-(f_{vu}(x)-f_{vv}(0))+q_{\gamma_{vu}}(x)(f_{uu}(0)-f_{vv}(0)).$

Thus, by combining Eq. 11 with the detailed balance equation (6) for the
boundary measure $\eta^{x}$,

$\displaystyle\Omega\Big{(}\int_{D_{e_{uv}}}\mathcal{L}f(z)\,\mu(dz)+\int_{D_{e_{vu}}}\mathcal{L}f(z)\,\mu(dz)\Big{)}$
$\displaystyle=\eta_{u}^{x}\lambda_{uv}q_{\gamma_{uv}}(x)(f_{vv}(0)-f_{uu}(0))+\eta_{v}^{x}\lambda_{vu}q_{\gamma_{vu}}(x)(f_{uu}(0)-f_{vv}(0))=0,$

and therefore $\int_{D}\mathcal{L}f(z)\,\mu(dz)=0$ in view of Eq. 10. ∎

######  Remark 1.

Peng and Li (2013) study diffusion processes defined on an open, bounded
domain $D^{0}$ in $\mathbb{R}^{d}$ with holding and jumping from a regular
boundary $\partial D$, and provide existence and uniqueness of a stationary
distribution under suitable assumptions on the regularity of the coefficients
of the generator of the diffusion process. Our result for the case of the
graph-valued process is parallel to Peng and Li (2013, Theorem 4.2).

##  3 Large population size scaling

The polymorphic segments of the path of $(X_{t})_{t\geq 0}$ through the
interior $D^{\circ}$ run on the time scale of evolution, which is a
characteristic of the Wright-Fisher diffusion. The generic re-scaling approach
behind the Wright-Fisher diffusion approximation considers the change in
frequency in a population of size $N$ over the time span of $N$ generations.
Simultaneously, the relevant selection coefficient at the level of
generations, $s$, is of the order $s\sim\gamma/N\to 0$, where $\gamma$ is the
selection coefficient of the limiting diffusion process. To properly adapt the
holding time distribution in the present model to the evolutionary time scale
we therefore introduce a parameter $N$ as a proxy of population size and
prescribe that the jumps into the interior of the state space have size
$x=1/N$. The time scale of the system is then set by the speed of mutation
$\lambda_{uv}/x=\lambda_{uv}N$. Our goal in this section is to analyze the
stationary distribution $\mu$ in Theorem 1 with $x=1/N$ for large but fixed
$N$. Specifically we identify the dominant terms in the asymptotic expansion
of $\mu$ under scaling for large $N$ and drop remainder terms of order $\ln
N/N$ and smaller. During this procedure it is convenient to make a number of
simplifying approximations valid formally in the limit $N\to\infty$. It is
important to keep in mind however that the population size proxy $N$ is kept
as a finite model parameter.

###  3.1 Approximation of the stationary distribution

We recall that in our model two vertices $u$ and $v$ are always connected by
two directed edges $e_{uv}$ and $e_{vu}$ whenever the jump rates between $u$
and $v$ are positive. For each edge $e_{uv}$, as $N\to\infty$, we introduce
the scaled fixation probability $\omega_{\gamma}$, where $\gamma=\gamma_{uv}$,
by

$q_{\gamma}(1/N)=\frac{\omega_{\gamma}}{N}+O\Big{(}\frac{1}{N^{2}}\Big{)},\quad\omega_{\gamma}=\frac{2\gamma}{1-e^{-2\gamma}},\quad\gamma\not=0,\qquad\omega_{0}=1.$
(12)

Due to assumption (2) on directional selection, we obtain the symmetry
relation

$\omega_{\gamma_{vu}}=\omega_{-\gamma_{uv}}=e^{-2\gamma_{uv}}\,\omega_{\gamma_{uv}}.$
(13)

As before, the collection of jump rates $\\{\lambda_{uv}\\}$ and selection
coefficients $\\{\gamma_{uv}\\}$ associated with the vertices $\mathcal{V}$
and edges $\mathcal{E}$ again define an embedded scaled continuous time Markov
chain on $\mathcal{V}$. The scaled generator matrix arises as
$N\mathcal{H}_{1/N}\sim\mathcal{H}$, with

$\mathcal{H}=(h_{uv}),\quad h_{uv}=\lambda_{uv}\omega_{\gamma_{uv}},\quad-
h_{uu}=\sum_{v\in\mathcal{V}}\lambda_{uv}\omega_{\gamma_{uv}},$

and is irreducible. In analogy with the previous relation (6), the unique
solution $\eta=\\{\eta(z),z\in\partial D\\}$ of the detailed balance equations

$\eta_{v}\lambda_{vu}\omega_{\gamma_{vu}}=\eta_{u}\lambda_{uv}\omega_{\gamma_{uv}},\quad
e_{uv}\in\mathcal{E},$ (14)

is the scaled invariant boundary distribution of the embedded Markov chain
with generator $\mathcal{H}$. The solution $\eta$ of (14), that no longer
depends on $N$, is a convenient approximation of the solution $\eta^{x}$ of
(6) with $x=1/N$. The distribution of the first jump averaged over the scaled
invariant measure,

$\nu^{1/N}=\sum_{u,v\in\mathcal{V}}\eta_{u}\lambda_{uv}\delta_{(u,v,1/N)},$
(15)

still depends on the initial mutation frequency $1/N$. The next result records
the dominant terms in Theorem 1, where we have fixed all mutation and
selection parameters, and then choose $x=1/N$ and $N$ large enough so that
remainder terms of order $O(\ln N/N)$ and smaller are removed.

######  Proposition 1.

The stationary single site distribution $\mu$ in Theorem 1 satisfies for large
$N$ the approximation

$\mu(z)=\mu_{N}(z)+O(1/N),\quad z\in D,$

where the approximating distribution $\mu_{N}$ has monomorphic site
probabilities

$\mu_{N}(u,u,0)=\frac{\eta_{u}}{\Omega^{\prime}_{N}},\quad u\in\mathcal{V},$

for $z=(u,u,0)$, polymorphic density given by

$\begin{split}\mu_{N}(u,v,y)\,dy=\frac{2\eta_{u}\lambda_{uv}}{\Omega^{\prime}_{N}}\Big{\\{}&\frac{\omega_{\gamma_{uv}}(1-e^{-2\gamma_{uv}(1-y)})}{2\gamma_{uv}y(1-y)}\,1_{\\{1/N<y<1\\}}\\\
&\quad+(N-\omega_{\gamma_{uv}})\,1_{\\{0<y<1/N\\}}\Big{\\}}\,dy,\end{split}$

for $z=(u,v,y)$, and is normalized by $\Omega^{\prime}_{N}=\Omega_{N}+O(\ln
N/N)$, with

$\Omega_{N}=1+2\sum_{u,v\in\mathcal{V}}\eta_{u}\lambda_{uv}(1+\ln
N+K_{\gamma_{uv}})$ (16)

and

$K_{\gamma}=\omega_{\gamma}\int_{0}^{1}(-\ln y)(e^{-2\gamma
y}-e^{-2\gamma(1-y)})\,dy.$

Before proving Proposition 1, we comment on some properties of the function
$K_{\gamma}$ and state the approximate distribution $\mu_{N}$ for the case of
neutral evolution.

######  Remark 2.

First, we have $K_{\gamma}\geq 0$, $\gamma\geq 0$. Second, the function
$K_{\gamma}$ is odd, $K_{-\gamma}=-K_{\gamma}$. In particular,
$K_{\gamma_{vu}}=K_{-\gamma_{uv}}=-K_{\gamma_{uv}}$. Third, $K_{\gamma}$ grows
logarithmically for large $\gamma$: with $\gamma_{e}=0.5772\dots$ denoting
Euler’s constant,

$K_{\gamma}\leq\left\\{\begin{array}[]{ll}\gamma,&0\leq\gamma\leq 1\\\
\gamma_{e}+\ln 2\gamma,&\gamma\geq 1,\end{array}\right.\qquad
K_{\gamma}\sim\gamma_{e}+\ln 2\gamma\quad\mbox{for large $\gamma$}.$

######  Remark 3.

For the special case of neutral evolution, $\gamma_{uv}=0$ for all
$e_{uv}\in\mathcal{E}$, we have $\omega_{0}=1$ and $K_{0}=0$, so

$\mu_{N}(u,u,0)=\frac{\eta_{u}}{\Omega_{N}^{\prime}},\quad u\in\mathcal{V},$

and

$\mu_{N}(u,v,y)\,dy=\frac{2\eta_{u}\lambda_{uv}}{\Omega_{N}^{\prime}}\Big{\\{}y^{-1}\,1_{\\{1/N<y<1\\}}+(N-1)\,1_{\\{0<y<1/N\\}}\Big{\\}}\,dy,$

for $(u,v,y)\in D^{\circ}$, where $(\eta_{u})$ is the unique solution of the
balance equations

$\eta_{v}\lambda_{vu}=\eta_{u}\lambda_{uv},\quad e_{uv}\in\mathcal{E},$

and $\Omega_{N}^{\prime}=\Omega_{N}+O(1/N)$ the normalization factor under
neutrality with

$\Omega_{N}=1+2(1+\ln N)\sum_{u\in\mathcal{V}}\eta_{u}\lambda_{u}.$

###### Proof.

For fixed $z^{\prime}=(u,v,y)$,

$G(\nu^{1/N},z^{\prime})=\eta_{u}\lambda_{uv}NG_{\gamma_{uv}}(1/N,y).$

Here, using Eq. 12 with large $N$ and $\gamma=\gamma_{uv}$,

$\displaystyle NG_{\gamma}\Big{(}\frac{1}{N},y\Big{)}$
$\displaystyle=Nq_{\gamma}\Big{(}\frac{1}{N}\Big{)}\frac{1-e^{-2\gamma(1-y)}}{\gamma
y(1-y)}\,1_{\\{1/N<y<1\\}}$
$\displaystyle\qquad+N\Big{(}1-q_{\gamma}\Big{(}\frac{1}{N}\Big{)}\Big{)}\frac{e^{2\gamma
y}-1}{\gamma y(1-y)}\,1_{\\{0<y<1/N\\}}$
$\displaystyle=\omega_{\gamma}\frac{1-e^{-2\gamma(1-y)}}{\gamma
y(1-y)}\,1_{\\{1/N<y<1\\}}$
$\displaystyle\qquad+2(N-\omega_{\gamma})\,1_{\\{0<y<1/N\\}}+O\Big{(}\frac{1}{N}\Big{)},$

from which we obtain $\Omega^{\prime}_{N}\mu_{N}(z)$, $z\in D^{\circ}$.
Moreover,

$\begin{split}&N\int_{0}^{1}G_{\gamma}\Big{(}\frac{1}{N},y\Big{)}\,dy=\omega_{\gamma}\int_{1/N}^{1}\frac{1-e^{-2\gamma(1-y)}}{\gamma
y(1-y)}\,dy+2+O\Big{(}\frac{1}{N}\Big{)}.\end{split}$

By partial integration the remaining integral evaluates to

$\displaystyle\int_{1/N}^{1}\\{y+$
$\displaystyle(1-y)\\}\frac{1-e^{-2\gamma(1-y)}}{\gamma y(1-y)}\,dy$
$\displaystyle=\int_{0}^{1-1/N}\frac{1-e^{-2\gamma y}}{\gamma
y}\,dy+\int_{1/N}^{1}\frac{1-e^{-2\gamma(1-y)}}{\gamma y}\,dy$
$\displaystyle=\frac{2\ln N}{\omega_{\gamma}}+2\int_{0}^{1}(-\ln
y)(e^{-2\gamma y}-e^{-2\gamma(1-y)})\,dy+O\Big{(}\frac{\ln N}{N}\Big{)}.$

Hence,

$N\int_{0}^{1}G_{\gamma}\Big{(}\frac{1}{N},y\Big{)}\,dy=2(1+\ln
N+K_{\gamma})+O\Big{(}\frac{\ln N}{N}\Big{)}.$

Integration over $D^{\circ}$ yields

$\int_{D^{0}}G(\nu^{1/N},z^{\prime})\,dz^{\prime}=2\sum_{u,v\in\mathcal{V}}\eta_{u}\lambda_{uv}(1+\ln
N+K_{\gamma_{uv}})+O\Big{(}\frac{\ln N}{N}\Big{)}.$

The representation of an approximate distribution $\mu_{N}(z)$ as stated now
follows from Theorem 1.

Finally, to verify the claims in Remark 2, for $\gamma>0$,

$\begin{split}\frac{K_{\gamma}}{\omega_{\gamma}}&=\int_{0}^{1/2}(-\ln
y)(e^{-2\gamma y}-e^{-2\gamma(1-y)})\,dy\\\
&\quad-\int_{0}^{1/2}(-\ln(1-y))(e^{-2\gamma y}-e^{-2\gamma(1-y)})\,dy\geq
0.\end{split}$

The relation $\omega_{-\gamma}=e^{-\gamma}\omega_{\gamma}$ implies the
symmetry $K_{-\gamma}=-K_{\gamma}$. Furthermore, the change-of-variable
$x=2\gamma y$ yields

$K_{\gamma}=\frac{\omega_{\gamma}}{2\gamma}\int_{0}^{2\gamma}(-\ln
x+\ln(2\gamma))e^{-x}\,dx-\frac{\omega_{\gamma}}{2\gamma}\int_{0}^{2\gamma}-\ln(1-x/2\gamma)\,e^{-x}\,dx.$

The rightmost integral is positive and tends to zero as $\gamma\to\infty$, by
an application of the monotone convergence theorem. Also,
$\omega_{\gamma}/(2\gamma)\to 1$ as $\gamma\to\infty$. Thus,

$K_{\gamma}-\ln(2\gamma)\sim\frac{\omega_{\gamma}}{2\gamma}\int_{0}^{2\gamma}(-\ln
x)e^{-x}\,dx\to\gamma_{e},\quad\gamma\to\infty.$

Remark 3 follows directly from Proposition 1 for $\gamma=0$. ∎

###  3.2 Allele frequency spectra

The unfolded allele frequency spectrum (AFS) describes the allele frequency
distribution of the derived allele at a biallelic locus. The AFS is a summary
statistic of the stationary distribution and can be retrieved using the
modeling setup with two directed edges between each pair of types and the
large $N$ approximation in Proposition 1. For a given locus, the unfolded AFS
representing the density of the derived allele frequency of any type
corresponds to $\sum_{u,v\in\mathcal{V}}\mu_{N}(u,v,y),0<y<1$. We visualize
the polymorphic density on two directed edges $e_{uv}$ and $e_{vu}$ of a
multi-allele model for $\gamma_{uv}=1$ in Fig. 2A.

Figure 2: Polymorphic densities of two types $u$ and $v$ with selection
coefficient $\gamma_{uv}=1$. Panel A: Relative log-scaled densities for type
$u$ at frequency $1-y$ on the directed edge $e_{uv}$, $\mu_{N}(u,v,y)$ in red,
and on the directed edge $e_{vu}$, $\mu_{N}(v,u,1-y)$ in blue. Panel B: The
relative log-scaled density on the edges of the pair $\langle u,v\rangle$ for
type $u$ at frequency $1-y$, $\mu_{N}(\langle u,v\rangle,1-y,y)$ in purple,
which is the sum of the two curves in panel A. The relative log-scaled density
for type $v$ at frequency $1-y$, $\mu_{N}(\langle v,u\rangle,1-y,y)$ in pink,
is symmetric to the purple curve around $y=0.5$.

In practice, the unfolded AFS relies on a polarization of polymorphisms into
derived and ancestral types, knowledge that requires additional information
such as outgroup data, which is not always readily available. If this is the
case a folded AFS can be derived from data. The folded AFS takes biallelic
observations and typically measures the minor allele at some frequency
$y\in[0,0.5]$, and the other allelic type at complementary frequency
$1-y\in[0.5,1]$. To formalize representations of unfolded and folded allele
frequency spectra using the stationary distribution in the current model, we
introduce the set of unordered pairs of vertices
$\mathcal{P}\coloneqq\\{\langle u,v\rangle\colon u,v\in\mathcal{V}\\}$, with
$|\mathcal{P}|=\binom{|\mathcal{V}|}{2}$. Restricting to the edges of the pair
$\langle u,v\rangle$, the polymorphic density of the process when type $v$ has
frequency $y$ and type $u$ frequency $1-y$ equals

$\mu_{N}(\langle u,v\rangle,1-y,y)\coloneqq\mu_{N}(u,v,y)+\mu_{N}(v,u,1-y).$

The folded density of the minor allele on $\langle u,v\rangle$ is therefore

$\mu_{\mathrm{fold}}(\langle u,v\rangle,y)\coloneqq\mu_{N}(\langle
u,v\rangle,1-y,y)+\mu_{N}(\langle v,u\rangle,1-y,y),\quad 0<y\leq 1/2.$

Figure 2B depicts the polymorphic density $\mu_{N}(\langle u,v\rangle,1-y,y)$
on edges connecting a pair $\langle u,v\rangle$.

######  Corollary 1.

We have

$\begin{split}\mu_{N}(\langle
u,v\rangle,1-y,y)\,dy&=\frac{2\eta_{u}\lambda_{uv}}{\Omega_{N}^{\prime}}\Big{\\{}\frac{e^{2\gamma_{uv}y}}{y(1-y)}\,1_{\\{1/N<y<1-1/N\\}}\\\
&\qquad+N\,1_{\\{0<y<1/N\\}}+Ne^{2\gamma_{uv}}\,1_{\\{1-1/N<y<1\\}}\Big{\\}}\,dy,\end{split}$

and

$\begin{split}\mu_{\mathrm{fold}}(\langle
u,v\rangle,y)\,dy=\frac{2\eta_{u}\lambda_{uv}}{\Omega_{N}^{\prime}}\Big{\\{}&\frac{e^{2\gamma_{uv}y}+e^{2\gamma_{uv}(1-y)}}{y(1-y)}\,1_{\\{1/N<y\leq
1/2\\}}\\\
&+N(1+e^{2\gamma_{uv}})\,1_{\\{0<y<1/N\\}}\Big{\\}}\,dy,\end{split}$

with the normalization factor $\Omega_{N}^{\prime}$ in Proposition 1.

###### Proof.

Consider a fixed pair of types $\langle u,v\rangle\in\mathcal{P}$. The density
of the process when $v$ has frequency $y$ follows from adding up the two
densities on each directed edge derived in Proposition 1,

$\begin{split}&\Omega_{N}^{\prime}(\mu_{N}(u,v,y)+\mu_{N}(v,u,1-y))\\\
&=2\eta_{u}\lambda_{uv}\left\\{\frac{\omega_{\gamma_{uv}}(1-e^{-2\gamma_{uv}(1-y)})}{2\gamma_{uv}y(1-y)}1_{\\{1/N<y<1\\}}+(N-\omega_{\gamma_{uv}})1_{\\{0<y<1/N\\}}\right\\}\,dy\\\
&\,+2\eta_{v}\lambda_{vu}\left\\{\frac{\omega_{\gamma_{vu}}(1-e^{-2\gamma_{vu}y})}{2\gamma_{vu}y(1-y)}1_{\\{0<y<1-1/N\\}}+(N-\omega_{\gamma_{vu}})1_{\\{1-1/N<y<1\\}}\right\\}\,dy.\end{split}$

Using detailed balance, Eq. 14, and the relationship
$\gamma_{vu}=-\gamma_{uv}$, the contribution from the interior, $1/N<y<1-1/N$,
for large $N$ is

$\begin{split}&\eta_{u}\lambda_{uv}\omega_{\gamma_{uv}}\left\\{\frac{1-e^{-2\gamma_{uv}(1-y)}}{\gamma_{uv}y(1-y)}+\frac{1-e^{-2\gamma_{vu}y}}{\gamma_{vu}y(1-y)}\right\\}=2\eta_{u}\lambda_{uv}\frac{e^{2\gamma_{uv}y}}{y(1-y)}.\end{split}$

For large $N$, the contribution from close to the boundary at zero, $0<y<1/N$,
is

$2\eta_{u}\lambda_{uv}(N-\omega_{\gamma_{uv}})+\eta_{v}\lambda_{vu}\frac{\omega_{\gamma_{vu}}(1-e^{-2\gamma_{vu}y})}{\gamma_{vu}y(1-y)}=2\eta_{u}\lambda_{uv}N+O\Big{(}\frac{1}{N}\Big{)},$

and from close to the boundary at one, $1-1/N<y<1$,

$\begin{split}&\eta_{u}\lambda_{uv}\frac{\omega_{\gamma_{uv}}(1-e^{-2\gamma_{uv}(1-y)})}{\gamma_{uv}y(1-y)}+2\eta_{v}\lambda_{vu}(N-\omega_{\gamma_{vu}})\\\
&=2\eta_{u}\lambda_{uv}Ne^{2\gamma_{uv}}+O\Big{(}\frac{1}{N}\Big{)},\end{split}$

where we again used detailed balance and relation Eq. 13. As the alternative
view of merging the directed edges simply entails reshuffling contributions,
the normalization constant does not change. The expression for
$\mu_{\mathrm{fold}}(\langle u,v\rangle,y)\,dy$ follows from similar
calculations. ∎

An example of the unfolded AFS in a four type model with state space shown in
Fig. 3A is given as the orange curve in Fig. 3B. Similarly as for the unfolded
AFS, the sum $\sum_{\langle
u,v\rangle\in\mathcal{P}}\mu_{\textrm{fold}}(\langle u,v\rangle,y),0<y\leq
1/2,$ yields the folded, type-independent distribution of derived allele
frequencies at a locus (Fig. 3B, red curve).

Figure 3: Allele frequency spectra for a four type model. Panel A: State space
for a model with four types. Panel B: Log-scaled unfolded (blue) and folded
(red) AFS. Parameters: population size $N=10^{4}$, selection coefficients
defined as fitness differences with $F_{u}=F_{z}=0$ and $F_{v}=F_{w}=1$, and
equal mutation intensity $\lambda=N\times 10^{-8}$ among all types.

###  3.3 Extension to multiple loci

The model introduced here applies directly to a collection of $L$ loci which
evolve independently. For the special case that a locus represents a single
site, a collection of consecutive sites represents a DNA sequence. Thus, even
though the following considerations are in general about a collection of
independent loci, we may use the term “sequence” instead. The state of the
sequence is defined by a collection of independent holding and jumping
diffusion processes $X^{j}$, $j=1,\dots,L$, with values in the direct product
set $\prod_{j=1}^{L}\mathcal{G}^{j}$, where all graphs $\mathcal{G}^{j}$ have
the same vertex set $\mathcal{V}$ and edge set $\mathcal{E}$. We allow the set
of selection coefficients $\gamma^{j}_{uv}$, $1\leq j\leq L$, to differ from
one locus to another, but assume that the transition rates $\lambda_{uv}$ are
the same among loci or along the sequence. For this it is convenient to
introduce scaled intensities $\theta_{uv}=\lambda_{uv}L$, $u,v\in\mathcal{V}$.
The results in Propositions 1 and 1 then apply with
$\lambda_{uv}=\theta_{uv}/L$. Summing over $v$, the total intensity
$\theta_{u}=\sum_{v\in\mathcal{V}}\theta_{uv}$ of a mutation from type $u$
becomes $\theta_{u}=\lambda_{u}L$, $u\in\mathcal{V}$, and so $N\theta_{u}$ is
the total rate in the collection of loci per time unit of a mutation affecting
$u$. In contrast, the steady-state probabilities $(\mu^{j}(z))_{z\in D}$ and
the boundary probabilities $(\eta^{j}_{u})_{u\in\mathcal{V}}$ typically vary
between loci, $1\leq j\leq L$. The total mutation rate on the boundary,
$\widehat{\theta}$, averaged across loci, is

$\widehat{\theta}=\frac{1}{L}\sum_{j=1}^{L}\widehat{\theta}^{j},\quad\widehat{\theta}^{j}=\sum_{u\in\mathcal{V}}\eta_{u}^{j}\theta_{u}.$
(17)

The quantities $\widehat{\theta}^{j}$ are the jump rates of the distribution
$\nu^{1/N}$ in Eq. 15. Clearly,

$\theta_{\mathrm{min}}:=\min_{u\in\mathcal{V}}\theta_{u}\leq\widehat{\theta}\leq\max_{u\in\mathcal{V}}\theta_{u}=:\theta_{\mathrm{max}}.$
(18)

The closely related summation

$\widehat{\theta}_{\mathrm{eff}}=\sum_{u\in\mathcal{V}}\widehat{\mu}_{u}\theta_{u},\quad\widehat{\mu}_{u}=\frac{1}{L}\sum_{j=1}^{L}\frac{\eta^{j}_{u}}{\Omega_{N}^{j}}$
(19)

represents the _effective mutation rate_ of the sequence, weighted by the
average probability that loci are monomorphic. Of course,
$\widehat{\theta}_{\mathrm{eff}}\leq\widehat{\theta}$. We say that the
mutation mechanism on the graph is _homogeneous_ if the total rates in each
vertex coincide, i.e. $\theta_{u}=\theta$, $u\in\mathcal{V}$. Under the
stronger assumption of homogeneous mutation,

$\widehat{\theta}=\theta,\qquad\widehat{\theta}_{\mathrm{eff}}=\frac{1}{L}\sum_{j=1}^{L}\frac{1}{\Omega_{N}^{j}}\,\theta.$
(20)

We summarize the steady state of the collection of loci by considering the
measure-valued process
$\mathcal{X}_{\infty}=\sum_{j=1}^{L}\delta_{X_{\infty}^{j}}$, where
$X^{j}_{\infty}$ has the stationary single locus distribution $\mu^{j}$ on
$\mathcal{G}^{j}$. The corresponding unfolded AFS across multiple loci is
given by

$\frac{1}{L}\sum_{j=1}^{L}\sum_{u,v\in\mathcal{V}}\mu_{N}^{j}(u,v,y),\quad
0<y<1.$

For suitable functions $f$, $f(z)=f_{uv}(y)$, $z\in D$,

$\langle\mathcal{X}_{\infty},f\rangle=\sum_{j=1}^{L}f(X_{\infty}^{j}),\quad
f\in\mathcal{F},$

represents the sequence equilibrium distribution. The steady-state expectation
under the approximate large population size distribution $\mu_{N}$ in
Proposition 1 is

$\begin{split}\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle=\sum_{j=1}^{L}\Big{\\{}&\sum_{u\in\mathcal{V}}f_{uu}(0)\mu^{j}_{N}(u,u,0)\\\
&\qquad+\sum_{u,v\in\mathcal{V}}\int_{0}^{1}f_{uv}(y)\,\mu^{j}_{N}(u,v,y)\,dy\Big{\\}}.\end{split}$
(21)

##  4 Impact of directional selection on genetic variation

There exists a number of summary statistics to assess genetic variation in a
population or population sample. These arise as the result of evaluating
functionals $\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle$ for
specifically chosen functions $f$ and can be analyzed by using Eq. 21. The
first term in the sum over $L$ in (21) provides the weight of the boundary
probabilities $(\mu^{j}(u,u,0))_{1\leq j\leq L}$ over monomorphic loci, and
the second term adds the relevant contributions from the allele frequency
spectrum of the polymorphic loci. We begin with a list of the basic instances
of such statistics. As a reference for each case we specialize to neutral
evolution and derive the relevant neutral summary statistics. Under the
assumption $\gamma_{uv}=0$ for every $u,v\in\mathcal{V}$, Eq. 21 simplifies
into

$\begin{split}&\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle=\frac{L}{\Omega_{N}}\sum_{u\in\mathcal{V}}f_{uu}(0)\eta_{u}\\\
&\quad+\frac{2}{\Omega_{N}}\sum_{u,v\in\mathcal{V}}\eta_{u}\theta_{uv}\Big{(}f_{uv}(0+)+\int_{1/N}^{1}\frac{f_{uv}(y)}{y}\,dy\Big{)}+\mathcal{R}_{N}\end{split}$
(22)

with

$\Omega_{N}=1+\frac{2(1+\ln
N)\widehat{\theta}^{\,0}}{L},\quad\widehat{\theta}^{\,0}=\sum_{u\in\mathcal{V}}\eta_{u}\theta_{u}\leq\theta_{\mathrm{max}},$

and $\mathcal{R}_{N}=O(\ln N/N)$. Here, $(\eta_{u})_{u\in\mathcal{V}}$ is the
solution of the neutral detailed balance equation, that is,
$\eta_{u}\theta_{uv}=\eta_{v}\theta_{vu}$, for every $u,v\in\mathcal{V}$, c.f.
Remark 3. Under homogeneous mutation, $\widehat{\theta}^{\,0}=\theta$.

###  4.1 Summary statistics under neutral evolution

The following listing is derived from (22) with the remainder term
$\mathcal{R}_{N}$ suppressed. The additional approximation assuming $\ln N/L$
is not too large, indicated by writing $\sim$ instead of $=$, falls within the
approximation range of $\mathcal{R}_{N}$.

  * i)

Average effective mutation rate

Define $f$ by $f_{uu}(0)=\theta_{u}$ and $f_{uv}(y)=0$. The expected value
$\widehat{\theta}_{\mathrm{eff}}^{\,0}=\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle/L$
is the average effective mutation rate per sequence under neutral evolution,
specifically taking into effect that mutations only occur on the boundary of
the graph. The first term in (22) yields

$\widehat{\theta}_{\mathrm{eff}}^{\,0}=\frac{{\widehat{\theta}}^{\,0}}{1+2L^{-1}(1+\ln
N)\widehat{\theta}^{\,0}}\sim{\widehat{\theta}}^{\,0}\Big{(}1-\frac{2(1+\ln
N)}{L}\widehat{\theta}^{\,0}\Big{)}.$

  * ii)

Polymorphic allele functionals

Using the effective mutation rate in i), we observe for functions $f$ that act
on edges only and are independent of the type, i.e. with $f_{uu}(0)=0$ and
$f_{uv}(y)=f(y)$, $u,v\in\mathcal{V}$, that (22) has the form

$\begin{split}\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle&=2\widehat{\theta}_{\mathrm{eff}}^{\,0}\Big{(}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,dy\Big{)}\\\
&\leq
2\widehat{\theta}^{\,0}\Big{(}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,dy\Big{)}.\end{split}$

  * iii)

Number of monomorphic sites

Let $f=f^{\partial}$ be the indicator function on the boundary $\partial D$.
The expected number of monomorphic sites out of $L$ is

$\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f^{\partial}\rangle=\frac{L}{\Omega_{N}}=\frac{L}{1+2L^{-1}(1+\ln
N)\widehat{\theta}^{\,0}}\sim L-2(1+\ln N)\,\widehat{\theta}^{\,0}.$

  * iv)

Number of polymorphic sites

Let $f^{\circ}=1-f^{\partial}$. For $L$ sufficiently large compared to $\ln N$
we obtain the familiar approximation of the expected number of polymorphic
sites as

$\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f^{\circ}\rangle=L-\frac{L}{\Omega_{N}}=2(1+\ln
N)\widehat{\theta}_{\mathrm{eff}}^{\,0}\leq 2(1+\ln N)\widehat{\theta}^{\,0}.$

  * v)

Number of segregating sites in a sample

To obtain the number of segregating sites in a sample of size $m$, we take
$f_{uu}(0)=0$ and $f_{uv}(y)=g_{m}(y)$, where

$g_{m}(y)=\sum_{k=1}^{m-1}\binom{m}{k}y^{k}(1-y)^{m-k}=1-y^{m}-(1-y)^{m},\quad
0\leq y\leq 1,$

is the probability that a sample of size $m\geq 2$ is polymorphic when drawn
from a population with derived frequency $y$. Then, it holds
$\int_{0}^{1/N}g_{m}(y)\,dy\sim O(1/N^{2})$ and

$\int_{1/N}^{1}y^{-1}g_{m}(y)\,dy\sim\int_{0}^{1}y^{-1}g_{m}(y)\,dy=\sum_{k=1}^{m-1}\frac{1}{k}.$

Hence

$\displaystyle S_{N,L}^{m}$
$\displaystyle=\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle\sim
2\widehat{\theta}_{\mathrm{eff}}^{\,0}\,\sum_{k=1}^{m-1}\frac{1}{k}\leq
2{\widehat{\theta}}^{\,0}\,\sum_{k=1}^{m-1}\frac{1}{k}.$

  * vi)

Pair-wise nucleotide differences

The standard measure of genetic diversity in the population, typically denoted
$\pi$, is the average number of pair-wise nucleotide differences normalized
per site (Nei and Li 1979). For sample size $m$ we take $f_{uu}(0)=0$ and
$f_{uv}(y)=h_{m}(y)$, where

$h_{m}(y)=\binom{m}{2}^{-1}\,\sum_{k=1}^{m-1}k(m-k)\binom{m}{k}y^{k}(1-y)^{m-k}=2y(1-y).$

Hence

$\pi=\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle/L=2\,\frac{\widehat{\theta}_{\mathrm{eff}}^{\,0}}{L}\leq
2\,\frac{\widehat{\theta}^{\,0}}{L},$

which turns out to be the average mutation load per site and is the same as
the expected proportion of segregating sites in a sample of size two,
$S_{N,L}^{2}/L$.

###  4.2 Allele frequency statistics and their upper bounds

We are interested in the overall effect of directional selection acting on the
functionals covered above and closely related quantities, as compared to their
counterparts under neutral evolution. Our main result shows that any presence
of directional selection among the alleles essentially benefits monomorphic
loci and constrains the number of polymorphic loci.

######  Theorem 2.

We consider the Wright-Fisher graph model extended to $L$ loci with fixed,
arbitrary parameters $(\theta_{uv})$ for mutation and
$\\{(\gamma_{uv}^{j})\\}_{1\leq j\leq L}$ for selection.

  * 1)

With $\widehat{\theta}$, $\widehat{\theta}_{\mathrm{eff}}$,
$\theta_{\mathrm{min}}$ and $\theta_{\mathrm{max}}$ defined in Eqs. 17, 19 and
18, we have

$\frac{\theta_{\mathrm{min}}}{1+2L^{-1}(1+\ln
N)\theta_{\mathrm{min}}}\leq\widehat{\theta}_{\mathrm{eff}}\leq\widehat{\theta}\leq\theta_{\mathrm{max}}.$

  * 2)

Let $f$ be a function on $D$ such that $f_{uv}(y)=f(y)\geq 0$ is a function
only of the frequency $y$ with $f(y)\to f(0+)\geq 0$, $y\to 0$, and $f(0)=0$.
Then

$\displaystyle 0\leq\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle$
$\displaystyle\leq
2\widehat{\theta}_{\mathrm{eff}}\,\Big{\\{}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,dy\Big{\\}}$
$\displaystyle\leq
2\widehat{\theta}\,\Big{\\{}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,dy\Big{\\}}.$
(23)

If, moreover, $\int_{0}^{1}y^{-1}f(y)\,dy<\infty$, then

$\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle\leq
2\widehat{\theta}_{\mathrm{eff}}\int_{0}^{1}\frac{f(y)}{y}\,dy\leq
2\widehat{\theta}\int_{0}^{1}\frac{f(y)}{y}\,dy.$ (24)

While the strength and direction of selection may vary arbitrarily within and
between sites, Theorem 2 illustrates that the effect of selective forces on
measures of genetic variation is only channeled through to the upper bounds
via the average mutation rates $\widehat{\theta}_{\mathrm{eff}}$ and
$\widehat{\theta}$, respectively. It is seen furthermore that the
proportionality constant $\widehat{\theta}_{\mathrm{eff}}$ is contained inside
an interval that does not depend on selection parameters, namely the interval
formed by the leftmost and the rightmost estimate in Theorem 2, 1). As a
corollary we observe that under the stronger assumption of homogeneous
mutation rates, introduced in Section 3.3, then $\widehat{\theta}$ will be
independent of any selective mechanisms, and the upper bounds in Items 2) and
24 will coincide with the corresponding expressions for neutral evolution in
Section 4.1 ii).

######  Corollary 2.

For the case when the mutation rates are homogeneous over all vertices, i.e.
$\theta_{u}=\theta$ for all $u\in\mathcal{V}$, then

$\widehat{\theta}=\widehat{\theta}^{\,0}=\theta$

and

$\frac{\theta}{1+2L^{-1}(1+\ln
N)\theta}\leq\widehat{\theta}_{\mathrm{eff}}\leq\theta.$

Hence,

$\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle\leq
2\theta\,\Big{\\{}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,dy\Big{\\}}$

and

$\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle\leq
2\theta\int_{0}^{1}\frac{f(y)}{y}\,dy,$

respectively.

###### Proof of Theorems 2 and 2.

1) For each single locus $j$,

$\displaystyle\Omega_{N}^{j}$
$\displaystyle=1+\frac{2}{L}\sum_{u,v\in\mathcal{V}}\eta_{u}^{j}\theta_{uv}(1+\ln
N+K_{\gamma^{j}_{uv}})$ $\displaystyle=1+\frac{2\widehat{\theta}^{j}}{L}(1+\ln
N)+\frac{2}{L}\sum_{u,v\in\mathcal{V}}\eta_{u}^{j}\theta_{uv}K_{\gamma^{j}_{uv}}.$

Here, by rewriting the double sum over all vertices in $\mathcal{V}$ as the
sum over all unordered pairs of vertices in $\mathcal{P}$ (see Section 3.2),

$\displaystyle\sum_{u,v\in\mathcal{V}}\eta_{u}^{j}\theta_{uv}K_{\gamma^{j}_{uv}}$
$\displaystyle=\sum_{\langle
u,v\rangle\in\mathcal{P}}\\{\eta_{u}^{j}\theta_{uv}K_{\gamma^{j}_{uv}}+\eta_{v}^{j}\theta_{vu}K_{\gamma^{j}_{vu}}\\}$
$\displaystyle=\sum_{\langle
u,v\rangle\in\mathcal{P}}\eta_{u}^{j}\theta_{uv}(1-e^{2\gamma_{uv}^{j}})K_{\gamma_{uv}^{j}}\leq
0,$

for every $\gamma_{uv}^{j}$. Hence $1\leq\Omega_{N}^{j}\leq
1+2\widehat{\theta}^{j}(1+\ln N)/L$ and therefore

$\begin{split}\widehat{\theta}\geq\widehat{\theta}_{\mathrm{eff}}&=\frac{1}{L}\sum_{j=1}^{L}\frac{\widehat{\theta}^{j}}{\Omega_{N}^{j}}\\\
&\geq\frac{1}{L}\sum_{j=1}^{L}\frac{\widehat{\theta}^{j}}{1+2\widehat{\theta}^{j}(1+\ln
N)/L}\geq\frac{\theta_{\mathrm{min}}}{1+2L^{-1}(1+\ln
N)\theta_{\mathrm{min}}}.\end{split}$

2) To prove the bound of $\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle$
in the second statement we take $f_{uu}(0)=f(0)=0$ in Eq. 21 and start from
the representation

$\displaystyle\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle$
$\displaystyle=\sum_{j=1}^{L}\sum_{u,v\in\mathcal{V}}\int_{0}^{1}f(y)\mu^{j}_{N}(u,v,y)\,dy.$

As we apply Proposition 1 it is convenient to have the auxiliary notation
$J_{uv}(y)$ (only used in this proof)

$J_{uv}(y)=\frac{1-e^{-2\gamma_{uv}(1-y)}}{2\gamma_{uv}(1-y)},\quad 0<y<1.$

Then

$\begin{split}&\int_{0}^{1}f(y)\mu^{j}_{N}(u,v,y)\,dy=\frac{2\eta_{u}^{j}\theta_{uv}}{L\Omega^{j}_{N}}\Big{\\{}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,\omega_{\gamma^{j}_{uv}}J_{uv}^{j}(y)\,dy\Big{\\}}.\end{split}$

We partition the right hand side as

$\displaystyle\int_{0}^{1}f(y)\mu^{j}_{N}(u,v,y)\,dy=\frac{2\eta_{u}^{j}\theta_{uv}}{L\Omega^{j}_{N}}\Big{\\{}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,dy\Big{\\}}-R_{N}^{j}(u,v),$

with

$R_{N}^{j}(u,v)=\frac{2\eta_{u}^{j}\theta_{uv}}{L\Omega^{j}_{N}}\int_{1/N}^{1}\frac{f(y)}{y}\Big{\\{}1-\omega_{\gamma_{uv}}^{j}J_{uv}^{j}(y)\Big{\\}}\,dy.$

Letting $R_{N}$ denote the sum

$R_{N}=\sum_{j=1}^{L}\sum_{u,v\in\mathcal{V}}R_{N}^{j}(u,v),$

these considerations imply

$\displaystyle\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle$
$\displaystyle=\sum_{u,v\in\mathcal{V}}\widehat{\mu}_{u}\,2\theta_{uv}\Big{\\{}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,dy\Big{\\}}-R_{N}.$

To complete the proof it remains to show that $R_{N}\geq 0$. Now,

$R_{N}=\sum_{j=1}^{L}\frac{2}{L\Omega_{N}^{j}}\int_{1/N}^{1}\frac{f(y)}{y}\sum_{u,v\in\mathcal{V}}\eta_{u}^{j}\theta_{uv}\Big{\\{}1-\omega_{\gamma^{j}_{uv}}J_{uv}^{j}(y)\Big{\\}}\,dy.$

Thus, it suffices to show, for each site $j$ and each frequency $y$,

$\sum_{u,v\in\mathcal{V}}\eta_{u}^{j}\theta_{uv}\big{\\{}1-\omega_{\gamma^{j}_{uv}}J_{uv}^{j}(y)\big{\\}}\geq
0.$

By rewriting the double summation as sum over all unordered pairs $\langle
u,v\rangle\in\mathcal{P}$, the previous inequality has the equivalent
representation

$\displaystyle\sum_{\langle
u,v\rangle\in\mathcal{P}}\Big{\\{}\eta_{u}^{j}\theta_{uv}\big{\\{}1-\omega_{\gamma^{j}_{uv}}J_{uv}^{j}(y)\big{\\}}+\eta_{v}^{j}\theta_{vu}\big{\\{}1-\omega_{\gamma^{j}_{vu}}J_{vu}^{j}(y)\big{\\}}\Big{\\}}\geq
0.$

By applying the detailed balance equation to each site and each edge, the task
is to show

$\displaystyle\sum_{\langle
u,v\rangle\in\mathcal{P}}\eta_{u}^{j}\theta_{uv}\Big{\\{}\big{\\{}1-\omega_{\gamma^{j}_{uv}}J_{uv}^{j}(y)\big{\\}}+\frac{\omega_{\gamma_{uv}^{j}}}{\omega_{\gamma_{vu}^{j}}}\big{\\{}1-\omega_{\gamma^{j}_{vu}}J_{vu}^{j}(y)\big{\\}}\Big{\\}}\geq
0.$

Equivalently,

$\displaystyle\sum_{\langle
u,v\rangle\in\mathcal{P}}\eta_{u}^{j}\theta_{uv}\omega_{\gamma^{j}_{uv}}R^{j}_{uv}\geq
0,\quad
R^{j}_{uv}=\omega_{\gamma^{j}_{uv}}^{-1}-J_{uv}^{j}(y)+\omega_{\gamma_{vu}^{j}}^{-1}-J_{vu}^{j}(y).$

Next we use the anti-symmetric relation (2) for the selection coefficients. If
we take a fixed site $j$ and an edge which connects two vertices, $u$ and $v$
say, and let $\gamma=\gamma^{j}_{uv}=-\gamma^{j}_{vu}$ be one of the relevant
selection coefficients, then it is straightforward to check that $R^{j}_{uv}$
is indeed nonnegative for any signed parameter $\gamma$ and $0<y<1$,

$R^{j}_{uv}=\frac{e^{2\gamma}-e^{-2\gamma}}{2\gamma}-\frac{e^{2\gamma(1-y)}-e^{-2\gamma(1-y)}}{2\gamma(1-y)}\geq
0.$

This verifies the claim $R_{N}\geq 0$ and yields

$\begin{split}\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle&\leq\sum_{j=1}^{L}\sum_{u,v\in\mathcal{V}}\frac{2\eta_{u}^{j}\theta_{uv}}{L\Omega^{j}_{N}}\Big{\\{}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,dy\Big{\\}}\\\
&=2\widehat{\theta}_{\mathrm{eff}}\Big{\\{}f(0+)+\int_{1/N}^{1}\frac{f(y)}{y}\,dy\Big{\\}}.\end{split}$

We note that this proof actually provides a more general result for functions
$f$ on $D$ that are not independent of $u$ and $v$ but fulfill
$f_{uv}(y)=f_{vu}(y)\geq 0$. Then

$0\leq\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle\leq
2\sum_{u,v\in\mathcal{V}}\widehat{\mu}_{u}\theta_{uv}\Big{\\{}f_{uv}(0+)+\int_{1/N}^{1}\frac{f_{uv}(y)}{y}\,dy\Big{\\}}.$

The other statement in part 2) of Theorem 2 is straightforward under the
additional assumption. For Corollary 2, a simple calculation verifies that
$\widehat{\theta}=\widehat{\theta}^{\,0}=\theta$ if $\theta_{u}=\theta$ for
every $u\in\mathcal{V}$. The rest of the corollary follows directly from
Theorem 2. ∎

###  4.3 The number of segregating sites and genetic diversity

We are now in position to consider concrete measures of genetic variation in a
population under the general model with selection and compare with the known
properties of these measures for neutral evolution as listed in Section 4.1.
Theorem 2 provides general estimates valid for arbitrary coefficients of
directional selection. First, Theorem 2 applied with the functions
$f^{\partial}$ and $f^{\circ}$ yield bounds which directly relate to the
listed items iii) and iv) of Section 4.1. In particular, the expected number
of segregating sites under selection satisfies

$\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f^{\circ}\rangle\leq
2\widehat{\theta}_{\mathrm{eff}}\,(1+\ln
N),\quad\widehat{\theta}_{\mathrm{eff}}=\sum_{u\in\mathcal{V}}\widehat{\mu}_{u}\theta_{u}.$
(25)

The parallel result for the number of segregating sites in a sample, i.e.
Theorem 2 applied with the function $f$ specified in item v), reads

$\displaystyle S_{N,L}^{m}$
$\displaystyle=\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle\leq
2\widehat{\theta}_{\mathrm{eff}}\,\sum_{k=1}^{m-1}\frac{1}{k}.$ (26)

Similarly, the expected genetic diversity in the population satisfies $\pi\leq
2\widehat{\theta}_{\mathrm{eff}}/L$, which extends vi) of Section 4.1.

For specific functions $f$ we may of course extract more detailed information
in addition to the upper bounds discussed here. It is again convenient to
carry out summation over unordered pairs of graph vertices. For this, let us
assume that $f_{uu}(0)=0$, $f_{uv}(y)=f(y)$, and
$\int_{0}^{1}y^{-1}f(y)\,dy<\infty$. Then

$\displaystyle\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle=$
$\displaystyle\quad\sum_{j=1}^{L}\sum_{\langle
u,v\rangle\in\mathcal{P}}\frac{2\eta_{u}^{j}\theta_{uv}\omega_{\gamma^{j}_{uv}}}{L\Omega_{N}^{j}}\int_{0}^{1}\frac{f(y)}{y(1-y)}\Big{(}\frac{e^{2\gamma_{uv}^{j}(1-y)}-e^{-2\gamma_{uv}^{j}(1-y)}}{2\gamma^{j}_{uv}}\Big{)}\,dy.$

In particular, for $f(y)=2y(1-y)$,

$\displaystyle\pi=\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f\rangle/L=\sum_{\langle
u,v\rangle\in\mathcal{P}}\frac{2\theta_{uv}}{L}\frac{1}{L}\sum_{j=1}^{L}\frac{2}{\Omega_{N}^{j}}\frac{\eta_{u}^{j}}{\omega_{\gamma_{vu}^{j}}}.$

The functional $f^{\circ}$ to obtain the expected number of segregating sites
does not fulfill the condition $\int_{0}^{1}y^{-1}f^{\circ}(y)\,dy<\infty$.
Nevertheless we obtain an explicit representation of the expected number of
segregating sites. Since the probability that a single site $j$ is polymorphic
is $(\Omega_{N}^{j}-1)/\Omega_{N}^{j}$, the expected number of segregating
sites in a sequence of length $L$ is

$\begin{split}&\mathbb{E}_{N}\langle\mathcal{X}_{\infty},f^{\circ}\rangle=\sum_{j=1}^{L}\frac{\Omega_{N}^{j}-1}{\Omega_{N}^{j}}\\\
&=\sum_{j=1}^{L}\frac{2\sum_{\langle
u,v\rangle\in\mathcal{P}}\eta_{u}^{j}\theta_{uv}\\{1+\ln
N+K_{\gamma_{uv}^{j}}+e^{2\gamma_{uv}^{j}}(1+\ln
N-K_{\gamma_{uv}^{j}})\\}}{L+2\sum_{\langle
u,v\rangle\in\mathcal{P}}\eta_{u}^{j}\theta_{uv}\\{1+\ln
N+K_{\gamma_{uv}^{j}}+e^{2\gamma_{uv}^{j}}(1+\ln
N-K_{\gamma_{uv}^{j}})\\}}.\end{split}$

##  5 Discussion

We have set up a multi-allele, multi-locus Wright-Fisher graph model to derive
rigorous upper bounds for a wide class of summary statistics of genetic
variation in Theorem 2. For any representative measure in this class the upper
bound is a multiple of the average effective mutation rate
$\widehat{\theta}_{\textrm{eff}}$. The multiplicative factor is independent of
directional selection and purely depends on the measure of genetic diversity.
Hence, mutation and directional selection only affect the upper bound through
$\widehat{\theta}_{\textrm{eff}}$ or $\widehat{\theta}$. To obtain selection-
independent upper bounds for arbitrary mutation rates, $\widehat{\theta}$ can
be replaced with e.g. $\theta_{\textrm{max}}$. The additional observation in
Corollary 2 that homogeneous mutation rates make $\widehat{\theta}$
independent of directional selection shows that the upper bounds are the same
as those for neutral evolution, and hence verifies the general presumption
that directional selection reduces genetic variation.

There exists a number of deterministic models to verify the reduction of
genetic variation due to directional selection (Feldman 1971; Novak and Barton
2017; Pontz and Feldman 2020), also referred to as ”constant frequency-
independent selection”. These models are based on replicator equations, that
were initially used by Feldman (1971) in this field. Within this deterministic
modeling approach analytical results on the interactions between loci due to
physical linkage and/or epistasis can be derived. On the other hand, mutation
and genetic drift as additional evolutionary forces are often not taken into
account. Here, we take a complementary approach and incorporate mutation and
genetic drift, while ignoring interactions among loci. In the following
section we will illustrate the relevance of this setting for studying the
interaction between mutation and fixation bias. Finally, we demonstrate that a
distribution of fitness effects of e.g. protein-coding sequences naturally
derives from our modeling setup.

###  5.1 Interaction between mutation and fixation bias

We say that there is a mutation bias between two allelic types
$u,v\in\mathcal{V}$, if mutation from one type to the other occurs more often
than in the reverse direction, i.e. if $\theta_{uv}\neq\theta_{vu}$, and there
is fixation bias between $u$ and $v$ whenever $\gamma_{uv}\not=0$. Thus,
fixation bias comprises non-zero directional selection. In addition to
selection, fixation bias can be caused by biased gene conversion, a mechanism
which is probabilistically equivalent to directional selection (Nagylaki
1983). The Wright-Fisher graph model allows to investigate the combined impact
of mutation bias and fixation bias on the stationary distribution and
population genetic measures derived as functionals of the stationary
distribution. Relevant combinations of mutation bias and fixation bias may act
in opposite directions and hence counterbalance their influence on the
stationary distribution or act in the same direction and thus reinforce each
other.

First of all we notice that given a collection of biased or unbiased mutation
rates $(\theta_{uv})$, any desired set of boundary distributions
$(\eta^{j}_{u})_{1\leq j\leq L}$ on $L$ loci is obtained from the detailed
balance equations (14) by putting

$\gamma^{j}_{uv}=\frac{1}{2}\ln\Big{(}\frac{\eta_{v}^{j}\theta_{vu}}{\eta_{u}^{j}\theta_{uv}}\Big{)}.$
(27)

As a reference case, the assignment
$2\gamma^{j}_{uv}=\ln(\theta_{vu}/\theta_{uv})$ yields uniform distributions
over all monomorphic states. These relations are well-known and appear in
studies of mutation-selection models such as Halpern and Bruno (1998); McVean
and Charlesworth (1999). Next, including polymorphic loci, for a scenario with
homogeneous mutation rates and prescribed selection parameters, Corollary 2
shows that genetic variation overall behaves much in the same way as for the
case with no fixation bias. Homogeneous mutation, however, is arguably not
necessarily realistic for genetic data, except perhaps on graphs of constant
degree, i.e., graphs with the same number of edges attached in each vertex.
The graph in Fig. 4 for the three stop codons provides a suitable example of a
graph of mixed degree. A particular mutation bias would be required to give
the relation $\theta_{TAG\to TAA}=\theta_{TAA\to TAG}+\theta_{TAA\to
TGA}=\theta_{TGA\to TAA}$ for homogeneous mutation.

Figure 4: The graph for stop codons. A specific mutation bias is required for
homogeneous mutation rates over all stop codons.

Hence, whenever the total mutation rates among types differ, it is worth
studying the combined impact of mutation bias and fixation bias on the upper
bounds in Theorem 2.

#### Four nucleotide model

The interaction of mutation and fixation bias is of interest in the four
nucleotide model (Fig. 3A, nodes representing nucleotides), as mutation rates
are frequently found to be different among the nucleotides (Stoltzfus and
Norris 2015; Long et al. 2018). Both biases impact for example the sequence
content, the abundance of each nucleotide in the sequence, which can be
obtained from Eq. 21.

The four nucleotide model can be reduced to a model with two types by grouping
the nucleotides into two classes. A natural classification arises when
studying nucleotide composition in the double-stranded DNA and considering A:T
base-pairs as one group and C:G base-pairs as another. This classification,
weak (A and T) and strong (C and G) bases, is commonly used to describe gBGC
(Duret and Galtier 2009; Mugal et al. 2015). The fixation bias towards GC over
AT nucleotides in the presence of gBGC interacts with the mutation bias
between the two classes, which acts in the opposite direction in several taxa,
i.e. mutations from GC to AT occur much more frequently than mutations from AT
to GC (Long et al. 2018). This illustrates that the two-type model can be
relevant to describe multiple alleles that can be classified into two types.
The traditionally studied two-type case is that of identically distributed
loci, $\gamma^{j}=\gamma$ for each $j=1,\ldots,L$ (Wright 1931; Li 1987;
Bulmer 1991; McVean and Charlesworth 1999). The resulting three-parameter
situation is frequently applied in models of gBGC (Muyle et al. 2011; De Maio
et al. 2013; Lachance and Tishkoff 2014) and analytical results of the
biallelic mutation-selection-drift model have been derived (Vogl and Clemente
2012; Vogl and Bergman 2015).

#### Interaction of mutation and fixation bias in a two-type model

We consider the graph with two types, $u$ and $v$, in which the dynamics at a
fixed locus $j$ are determined by three parameters $\theta_{u}$, $\theta_{v}$
and $\gamma^{j}=\gamma_{uv}^{j}$, $j=1,\ldots,L$. Clearly, by (14),

$(\eta^{j}_{u},\eta^{j}_{v})=\Big{(}\frac{\theta_{v}}{\theta_{v}+\theta_{u}e^{2\gamma^{j}}},\frac{\theta_{u}e^{2\gamma^{j}}}{\theta_{v}+\theta_{u}e^{2\gamma^{j}}}\Big{)}.$
(28)

Under the assumption $\gamma^{j}=\gamma$ for all sites, we may summarize the
effect of varying $\theta_{u}$, $\theta_{v}$, and $\gamma$ in Eq. 28, as

$\frac{\eta_{u}}{\eta_{v}}=\frac{\theta_{v}}{\theta_{u}}e^{-2\gamma}\begin{cases}>1&\mathrm{if}\,\gamma<\frac{1}{2}\ln(\theta_{v}/\theta_{u}),\\\
=1&\mathrm{if}\,\gamma=\frac{1}{2}\ln(\theta_{v}/\theta_{u}),\\\
<1&\mathrm{if}\,\gamma>\frac{1}{2}\ln(\theta_{v}/\theta_{u}).\end{cases}$

To illustrate the results obtained in Theorem 2, the upper bounds are
controlled by

$\displaystyle\widehat{\theta}_{\mathrm{eff}}=\frac{1}{L}\sum_{j=1}^{L}\frac{1}{\Omega_{N}^{j}}\frac{\theta_{u}\theta_{v}(1+e^{2\gamma^{j}})}{\theta_{v}+\theta_{u}e^{2\gamma^{j}}}\leq\widehat{\theta}=\frac{1}{L}\sum_{j=1}^{L}\frac{\theta_{u}\theta_{v}(1+e^{2\gamma^{j}})}{\theta_{v}+\theta_{u}e^{2\gamma^{j}}}.$

To some degree these averaged mutation rates still depend on the selection
coefficients. However,

$\widehat{\theta}_{\mathrm{eff}}\leq\theta_{\max}\,\frac{1}{L}\sum_{j=1}^{L}\frac{1}{\Omega_{N}^{j}}\leq\theta_{\max},$

where the averaging sum is an estimate of the fraction of monomorphic loci in
the sequence. It is of interest nonetheless to compare with the harmonic mean
of the mutation rates appearing under neutral evolution, namely

$\displaystyle\widehat{\theta}^{\,0}_{\mathrm{eff}}=\frac{1}{\Omega_{N}}\frac{2\theta_{u}\theta_{v}}{\theta_{v}+\theta_{u}}\leq\widehat{\theta}^{\,0}=\frac{2\theta_{u}\theta_{v}}{\theta_{v}+\theta_{u}}.$

Considering the ratio $\widehat{\theta}/\widehat{\theta}^{\,0}$ for
$\gamma^{j}=\gamma$, $j=1,\ldots,L$, we obtain the relations

$\frac{\widehat{\theta}}{\widehat{\theta}^{\,0}}=\frac{(\theta_{u}+\theta_{v})(1+e^{-2\gamma})}{2(\theta_{u}+\theta_{v}e^{-2\gamma})}\begin{cases}>1&\mathrm{if}\,(\theta_{u}>\theta_{v}\land\gamma<0)\lor(\theta_{u}<\theta_{v}\land\gamma>0),\\\
=1&\mathrm{if}\,\gamma=0\lor\theta_{u}=\theta_{v},\\\
<1&\mathrm{if}\,(\theta_{u}<\theta_{v}\land\gamma<0)\lor(\theta_{u}>\theta_{v}\land\gamma>0).\end{cases}$

This implies that $\widehat{\theta}>\widehat{\theta}^{\,0}$ if mutation bias
is opposing fixation bias and $\widehat{\theta}<\widehat{\theta}^{\,0}$ if
mutation and fixation biases enhance each other.

Figure 5: Genetic diversity $\pi$ (solid lines) and its upper bound
$2\widehat{\theta}_{\mathrm{eff}}/L$ (dashed lines) in panel A and the
expected number of segregating sites and its upper bound
$2\widehat{\theta}_{\mathrm{eff}}(1+\ln N)$ (dashed lines) in panel B for
different combinations of mutation parameters in a two-type model with equal
selection coefficient among loci. Gray curves: $\theta_{u}=\theta_{v}=1.5\cdot
10^{-3}$, green curves: $\theta_{u}=3\cdot 10^{-3}>\theta_{v}=1\cdot 10^{-3}$,
pink curves: $\theta_{u}=1\cdot 10^{-3}<\theta_{v}=3\cdot 10^{-3}$. Other
parameters: population size $N=500$ and number of loci $L=1500$.

Such insights can be used together with the results of Theorem 2, for
instance, considering the case of genetic diversity (Fig. 5),

$\begin{split}\pi&=\frac{1}{L\Omega_{N}}\frac{4\theta_{u}\theta_{v}}{\theta_{v}+\theta_{u}e^{2\gamma}}\frac{e^{2\gamma}-1}{2\gamma}\leq
2\frac{\widehat{\theta}_{\mathrm{eff}}}{L}=\frac{1}{L\Omega_{N}}\frac{2\theta_{u}\theta_{v}(1+e^{2\gamma})}{\theta_{v}+\theta_{u}e^{2\gamma}}.\end{split}$

Without mutation bias or with a mutation bias that enhances the fixation bias,
genetic diversity decreases monotonically as selection becomes stronger (gray
and green solid curves in Fig. 5A). If mutation bias counteracts fixation
bias, genetic diversity first increases in the weak selection regime compared
to neutral evolution until a maximum is reached for an intermediate selection
coefficient, and decreases thereafter for stronger selection (pink solid curve
in Fig. 5A). A similar behavior is observed and discussed in McVean and
Charlesworth (1999). The upper bound (dashed lines in Fig. 5A) is constant for
equal mutation rates, decreases monotonically if mutation and fixation bias
reinforce each other, and increases monotonically for counterbalancing biases.
The behavior of the expected number of segregating sites and its upper bounds
under the different combinations of mutation rates is very akin to the curves
for genetic diversity (Fig. 5B).

The scenario depicted here where all loci have equal selective pressure that
can become arbitrarily large is rather artificial. In many taxa the genome-
wide average of fixation bias in gBGC takes a value in the weak selection
regime (De Maio et al. 2013; Glémin et al. 2015; Galtier et al. 2018; Boman et
al. 2021). Likewise, according to the nearly neutral theory (Ohta 1973, 1976,
1992) polymorphisms segregate in a population if selection is neutral or
nearly neutral. In this selection regime the upper bounds capture the behavior
of the measure of genetic variation well. Only when selection becomes strong,
the upper bounds become more generous. However, strong selection immediately
removes genetic variation and consequently, the interaction of mutation and
fixation bias in the strong selection regime is less relevant when considering
a large collection of loci.

###  5.2 Distribution of fitness effects

The present model equipped in each locus with a static fitness landscape as in
Eq. 3 can be applied for example to protein-coding sequences. In this case
every locus represents a nucleotide triplet where the fitness of each possible
codon type is particular to the specific locus. However, the number of fitness
parameters required to capture the selection effects over many loci and all
codons may quickly grow out of hand (for a recent review on this topic see
Youssef et al. 2021). An alternative view of modeling natural selection is
that of choosing fitness parameters from a representative distribution. In
this approach, the distribution of fitness effects (DFE) should reduce the
parameter space but preserve some of the characteristic features. Generally,
the DFE is composed of a distribution of negative selection coefficients, a
distribution of positive selection coefficients, and a proportion of selection
coefficients at zero representing neutral evolution. The steady-state of the
process $(X_{t})_{t\geq 0}$ in our modeling framework allows deriving a DFE
which is informative about the fraction of mutations in equilibrium that are
beneficial, deleterious and neutral, respectively.

#### The distribution of fitness effects of novel mutations

The common understanding of a DFE is the distribution of fitness effects of
all novel mutations that occur in a population (Eyre-Walker and Keightley
2007). As we apply a boundary mutation model, the weights of the selection
coefficients are given by the steady states of the Wright-Fisher graph model
on the boundary, $\eta_{u}^{j}/\Omega_{N}^{j},u\in\mathcal{V}$, and the
corresponding mutation intensities, $\theta_{uv},u,v\in\mathcal{V}$. We may
assume that the mutation intensities have been reduced in advance to discount
for the presence of any strongly deleterious mutations. Hence, the DFE on the
graph is a discrete probability distribution function on the real line,

$H_{\textrm{dfe}}(\gamma)=\widehat{\theta}_{\mathrm{eff}}^{-1}\,\frac{1}{L}\sum_{j=1}^{L}\sum_{u,v\in\mathcal{V}}\frac{\eta_{u}^{j}}{\Omega_{N}^{j}}\theta_{uv}1_{\\{\gamma_{uv}^{j}\leq\gamma\\}},\quad-\infty<\gamma<\infty,$
(29)

with jumps at each of the finite number of values
$\gamma_{uv}^{j},u,v\in\mathcal{V},j=1,\ldots,L$. In Proposition 2 we obtain
some key properties of the DFE.

######  Proposition 2.

For any choice of selection coefficients on the Wright-Fisher graph model,

  * i)

the probability weight for each negative selection coefficient,
$\gamma_{uv}^{j}<0$, is always larger or equal the weight of the corresponding
positive value, that is,

$H_{\mathrm{dfe}}(\gamma_{uv}^{j})-H_{\mathrm{dfe}}(\gamma_{uv}^{j}-)\geq
H_{\mathrm{dfe}}(-\gamma_{uv}^{j})-H_{\mathrm{dfe}}(-\gamma_{uv}^{j}-),$

  * ii)

the average selection load is deleterious, in the sense

$\widehat{\gamma}=\widehat{\theta}_{\mathrm{eff}}^{-1}\,\frac{1}{L}\sum_{j=1}^{L}\sum_{u,v\in\mathcal{V}}\frac{\eta_{u}^{j}}{\Omega_{N}^{j}}\theta_{uv}\gamma_{uv}^{j}\leq
0,$

  * iii)

the total contribution to positive selection is

$1-H_{\mathrm{dfe}}(0)<1/2.$

###### Proof.

To prove i) we use the alternative representation of $H_{\mathrm{dfe}}$
summing over pairs of vertices and apply detailed balance on the edges between
them,

$\displaystyle H_{\mathrm{dfe}}(\gamma)$
$\displaystyle=\widehat{\theta}_{\mathrm{eff}}^{-1}\frac{1}{L}\sum_{j=1}^{L}\sum_{\langle
u,v\rangle\in\mathcal{P}}\frac{1}{\Omega_{N}^{j}}\Big{\\{}\eta_{u}^{j}\theta_{uv}1_{\\{\gamma_{uv}^{j}\leq\gamma\\}}+\eta_{v}^{j}\theta_{vu}1_{\\{\gamma_{vu}^{j}\leq\gamma\\}}\Big{\\}}$
$\displaystyle=\widehat{\theta}_{\mathrm{eff}}^{-1}\frac{1}{L}\sum_{j=1}^{L}\sum_{\langle
u,v\rangle\in\mathcal{P}}\frac{\eta_{u}^{j}}{\Omega_{N}^{j}}\theta_{uv}\Big{\\{}1_{\\{\gamma_{uv}^{j}\leq\gamma\\}}+e^{2\gamma_{uv}^{j}}1_{\\{-\gamma_{uv}^{j}\leq\gamma\\}}\Big{\\}}.$

For each fixed negative selection coefficient $\gamma=\gamma_{uv}^{j}<0$, we
conclude that the jump size
$H_{\mathrm{dfe}}(\gamma)-H_{\mathrm{dfe}}(\gamma-)$ is proportional to $1$,
while the accompanying jump on the side of positive selection coefficients,
$H_{\mathrm{dfe}}(-\gamma)-H_{\mathrm{dfe}}(-\gamma-)$, equals the same
proportionality constant times $e^{2\gamma}$. Hence property i) follows since
$1\geq e^{2\gamma}$. For claim ii) we have, similarly,

$\displaystyle\widehat{\gamma}$
$\displaystyle=\widehat{\theta}_{\mathrm{eff}}^{-1}\frac{1}{L}\sum_{j=1}^{L}\sum_{\langle
u,v\rangle\in\mathcal{P}}\frac{\eta_{u}^{j}}{\Omega_{N}^{j}}\theta_{uv}\gamma_{uv}^{j}\\{1-e^{2\gamma^{j}_{uv}}\\}.$

This quantity is nonpositive, since
$\gamma_{uv}^{j}(1-e^{2\gamma_{uv}^{j}})\leq 0$ for all
$\gamma_{uv}^{j},j=1,\ldots,L$. It follows that $\widehat{\gamma}\leq 0$.
Claim iii) is merely a rephrased conclusion of i). ∎

To help interpret the shape of $H_{\mathrm{dfe}}(\gamma)$, let us suppose that
the $\gamma$-values for negative selection can be well approximated by a
continuous distribution with density function $g(x)$, $x<0$, on the negative
half line. It then follows from the proof of Proposition 2 that the associated
positive $\gamma$-values have a density proportional to $g(-x)e^{-2x}$, $x>0$.
The full approximation on the real line is then obtained by normalizing the
contributions for $x<0$, $x>0$ and possibly an atom for neutral mutations at
$x=0$. Figure 6 displays the resulting densities for the case of $g$ being an
exponential, $g(x)=e^{-|x|}$ (red), or a gamma distribution,
$g(x)=|x|^{a-1}e^{-|x|}/\Gamma(a)$ (blue for $a=2$ and light blue for
$a=0.15$).

Figure 6: Approximative density of the DFE based on an exponential
distribution with parameter 1 (red), a gamma distribution with shape parameter
2 and scale parameter 1 (blue), and a gamma distribution with shape parameter
0.15 and scale parameter 1 (light blue) for the negative selection
coefficients. The average selection load is $\widehat{\gamma}=-2/3$ in case of
the exponential distribution and $\widehat{\gamma}=-26/15$ and
$\widehat{\gamma}\approx-0.058$, respectively, for the gamma distributions.

Moreover, Proposition 2 provides an approximation of the present model, where
ancestral and derived alleles are kept distinct but not the genetic types. We
consider the general Wright-Fisher graph model with average mutation load
$\widehat{\theta}_{\mathrm{eff}}$ and DFE represented by
$H_{\mathrm{dfe}}(\gamma)$. In each locus, lump together all boundary states
into one generic vertex state $0$, which represents ”ancestral”. Attach to $0$
a single, directed edge $e_{0}$ of length $1$ which returns to $0$ at the
endpoint. At exponential mutation rate $\widehat{\theta}_{\mathrm{eff}}$ a
derived allele appears at frequency $1/N$ on the edge, following the path of a
Wright-Fisher diffusion with selection coefficient drawn randomly from the
distribution $H_{\mathrm{dfe}}$. At the time of extinction or fixation the
derived returns to the state of ancestral. The same dynamics applies
independently over the $L$ loci. The unfolded AFS of derived alleles in
steady-state should then be well approximated by

$2\widehat{\theta}_{\mathrm{eff}}\int_{\mathbb{R}}\omega_{\gamma}\frac{1-e^{-2\gamma(1-y)}}{2\gamma
y(1-y)}\,H_{\mathrm{dfe}}(d\gamma),\quad 0<y<1.$

The category of probabilistic models briefly described here is known as the
Poisson random field approach in population genetics, see e.g. Sawyer and
Hartl (1992); Mugal et al. (2014); Kaj and Mugal (2016). Proposition 2
provides general support for the Poisson random field, and Proposition 2 ii)
even justifies the further simplified approximation where $H_{\mathrm{dfe}}$
is a one-point distribution with unit mass on a fixed
$\gamma=\widehat{\gamma}<0$ for each mutation.

#### The distribution of fitness effects of polymorphisms

While in our setting mutation is linked to the boundary states, selection
naturally acts on the interior of the graph. This suggests introducing the
relevant distribution function for selection coefficients of segregating
polymorphisms

$\tilde{H}_{\mathrm{pdfe}}(\gamma)=\frac{1}{L}\sum_{j=1}^{L}\sum_{u,v\in\mathcal{V}}1_{\\{\gamma_{uv}^{j}\leq\gamma\\}}\int_{0}^{1}\mu^{j}_{N}(u,v,y)\,dy,\quad-\infty<\gamma<\infty.$
(30)

The function $\tilde{H}_{\mathrm{pdfe}}$ is a (improper) discrete distribution
on the real line with a finite number of jumps at each of the values
$\pm\gamma_{uv}^{j}$ of the graph. Since

$\tilde{H}_{\mathrm{pdfe}}(\gamma)\to\tilde{H}_{\mathrm{pdfe}}(\infty)=\frac{1}{L}\sum_{j=1}^{L}\frac{\Omega^{j}_{N}-1}{\Omega^{j}_{N}}<1,\quad$

as $\gamma\to\infty$, the distribution is defect. Conditioning on polymorphic
states yields a proper probability distribution
$H_{\textrm{pdfe}}(\gamma)=\tilde{H}_{\mathrm{pdfe}}(\gamma)/\tilde{H}_{\mathrm{pdfe}}(\infty)$.
It can be shown that Proposition 2 holds analogously for
$H_{\textrm{pdfe}}(\gamma)$.

###### Proof of Proposition 2 for $H_{\textrm{pdfe}}$.

Summing over pairs of vertices, $H_{\mathrm{pdfe}}$ can be represented as

$\displaystyle H_{\mathrm{pdfe}}(\gamma)$
$\displaystyle=\tilde{H}_{\mathrm{pdfe}}(\infty)^{-1}\frac{1}{L}\sum_{j=1}^{L}\sum_{\langle
u,v\rangle\in\mathcal{P}}\frac{2}{L\Omega^{j}_{N}}\eta_{u}^{j}\theta_{uv}\Big{\\{}(1+\ln
N+K_{\gamma_{uv}^{j}})1_{\\{\gamma_{uv}^{j}\leq\gamma\\}}$
$\displaystyle\hskip 85.35826pt+(1+\ln
N-K_{\gamma_{uv}^{j}})e^{2\gamma_{uv}^{j}}1_{\\{-\gamma_{uv}^{j}\leq\gamma\\}}\Big{\\}}.$

For each fixed negative selection coefficient $\gamma=\gamma_{uv}^{j}<0$, the
jump size $H_{\mathrm{pdfe}}(\gamma)-H_{\mathrm{pdfe}}(\gamma-)$ is
proportional to $1+\ln N+K_{\gamma}$, whereas the jump on the side of positive
selection coefficients,
$H_{\mathrm{pdfe}}(-\gamma)-H_{\mathrm{pdfe}}(-\gamma-)$, equals the same
proportionality constant times $(1+\ln N-K_{\gamma})e^{2\gamma}$. Hence
property i) in Proposition 2 for $H_{\mathrm{pdfe}}$ follows if we can show

$1+\ln N+K_{\gamma}\geq(1+\ln N-K_{\gamma})e^{2\gamma},$

that is

$-K_{\gamma}\leq\frac{1-e^{2\gamma}}{1+e^{2\gamma}}\,(1+\ln N),\quad\gamma<0,$
(31)

or, equivalently,

$K_{\gamma}\leq\frac{1-e^{-2\gamma}}{1+e^{-2\gamma}}\,(1+\ln
N),\quad\gamma>0.$ (32)

For $0<\gamma<1$, using Remark 2, we have e.g.

$K_{\gamma}\leq\gamma\leq\frac{1-e^{-2\gamma}}{1+e^{-2\gamma}}\,(1+\ln 2)$

and hence (32) holds for $N\geq 2$. Let $C$ be a constant such that
$|\gamma_{uv}^{j}|\leq C$ for all $1\leq j\leq L$ and $u,v\in\mathcal{V}$. For
$1\leq\gamma\leq C$, by Remark 2,

$K_{\gamma}\leq\gamma_{e}+\ln
2\gamma\leq\frac{1-e^{-2\gamma}}{1+e^{-2\gamma}}+\ln
2\gamma\leq\frac{1-e^{-2\gamma}}{1+e^{-2\gamma}}\big{(}1+\ln(4C^{2})\big{)},$

and hence (32) holds for $N\geq 4C^{2}$.

For claim ii) we have,

$\displaystyle\widehat{\gamma}$
$\displaystyle=\tilde{H}_{\mathrm{pdfe}}(\infty)^{-1}\frac{1}{L}\sum_{j=1}^{L}\sum_{\langle
u,v\rangle\in\mathcal{P}}\frac{2}{L\Omega^{j}_{N}}\eta_{u}^{j}\theta_{uv}\Big{\\{}(1+\ln
N+K_{\gamma^{j}_{uv}})\gamma_{uv}^{j}$ $\displaystyle\hskip 85.35826pt-(1+\ln
N-K_{\gamma^{j}_{uv}})\gamma_{uv}^{j}e^{2\gamma^{j}_{uv}}\Big{\\}}.$

The quantity in curly brackets is nonpositive if

$\gamma K_{\gamma}\leq\gamma\,\frac{1-e^{-2\gamma}}{1+e^{-2\gamma}}\,(1+\ln
N),$

which is true for all $\gamma=\gamma_{uv}^{j},j=1,\ldots,L$, and $N$
sufficiently large by the inequalities (31) and (32). Hence, it follows that
$\widehat{\gamma}\leq 0$. ∎

In conclusion, Proposition 2 demonstrates that both DFEs, that of novel
mutations as well as that of segregating polymorphisms, are at equilibrium
skewed towards the negative selection regime.

### Acknowledgements

The authors thank Nicolas Lartillot and Thibault Latrille for valuable
discussions about the use of mutation-selection models for protein-coding
sequence evolution. CFM has received financial support from the Knut and Alice
Wallenberg Foundation (2014/0044 to Hans Ellegren) and the Swedish Research
Council (2013-8271 to Hans Ellegren).

## References

  * Boman et al. (2021) Boman, J., Mugal, C. F., and Backström, N. 2021.  The effects of GC-biased gene conversion on patterns of genetic diversity among and across butterfly genomes.  Genome Biology and Evolution, 13(5): evab064.
  * Borges et al. (2019) Borges, R., Szöllősi, G. J., and Kosiol, C. 2019.  Quantifying GC-biased gene conversion in great ape genomes using polymorphism-aware models.  Genetics, 212(4): 1321–1336.
  * Bulmer (1991) Bulmer, M. 1991.  The selection-mutation-drift theory of synonymous codon usage.  Genetics, 129(3): 897–907.
  * Burden and Tang (2016) Burden, C. J. and Tang, Y. 2016.  An approximate stationary solution for multi-allele neutral diffusion with low mutation rates.  Theoretical Population Biology, 112: 22–32.
  * Cao et al. (2015) Cao, M., et al. 2015\.  Analysis of human triallelic SNPs by next-generation sequencing.  Annals of Human Genetics, 79(4): 275–281.
  * De Maio et al. (2013) De Maio, N., Schlötterer, C., and Kosiol, C. 2013.  Linking great apes genome evolution across time scales using polymorphism-aware phylogenetic models.  Molecular Biology and Evolution, 30(10): 2249–2262.
  * Duret and Galtier (2009) Duret, L. and Galtier, N. 2009.  Biased gene conversion and the evolution of mammalian genomic landscapes.  Annual Review of Genomics and Human Genetics, 10(1): 285–311.
  * Durrett (2008) Durrett, R. 2008.  Probability Models for DNA Sequence Evolution.  Springer, New York.
  * Etheridge (2011) Etheridge, A. 2011.  Some mathematical models from population genetics: École d’Été de Probabilités de Saint-Flour XXXIX-2009, volume 2012.  Springer Science & Business Media.
  * Ewens (2004) Ewens, W. J. 2004.  Mathematical Population Genetics 1: Theoretical Introduction.  Springer Verlag, Berlin.
  * Eyre-Walker and Keightley (2007) Eyre-Walker, A. and Keightley, P. D. 2007.  The distribution of fitness effects of new mutations.  Nature Reviews Genetics, 8(8): 610–618.
  * Feldman (1971) Feldman, M. W. 1971.  Equilibrium studies of two locus haploid populations with recombination.  Theoretical Population Biology, 2(3): 299–318.
  * Ferguson and Buzbas (2018) Ferguson, J. M. and Buzbas, E. O. 2018.  Inference from the stationary distribution of allele frequencies in a family of Wright–Fisher models with two levels of genetic variability.  Theoretical Population Biology, 122: 78–87.
  * Fisher (1930) Fisher, R. A. 1930.  The genetical theory of natural selection.  Clarendon Press, Oxford.
  * Galtier et al. (2018) Galtier, N., et al. 2018\.  Codon usage bias in animals: Disentangling the effects of natural selection, effective population size, and GC-biased gene conversion.  Molecular Biology and Evolution, 35(5): 1092–1103.
  * Glémin et al. (2015) Glémin, S., et al. 2015\.  Quantification of GC-biased gene conversion in the human genome.  Genome Research, 25(8): 1215–1228.
  * Halpern and Bruno (1998) Halpern, A. L. and Bruno, W. J. 1998.  Evolutionary distances for protein-coding sequences: modeling site-specific residue frequencies.  Molecular Biology and Evolution, 15(7): 910–917.
  * Kaj and Mugal (2016) Kaj, I. and Mugal, C. F. 2016.  The non-equilibrium allele frequency spectrum in a Poisson random field framework.  Theoretical Population Biology, 111: 51–64.
  * Karlin and Taylor (1981) Karlin, S. and Taylor, H. E. 1981.  A second course in stochastic processes.  Academic Press, New York.
  * Kimura (1962) Kimura, M. 1962.  On the probability of fixation of mutant genes in a population.  Genetics, 47(6): 713–719.
  * Kimura (1964) Kimura, M. 1964.  Diffusion models in population genetics.  Journal of Applied Probability, 1(2): 177–232.
  * Kimura (1969) Kimura, M. 1969.  The number of heterozygous nucleotide sites maintained in a finite population due to steady flux of mutations.  Genetics, 61(4): 893–903.
  * Lachance and Tishkoff (2014) Lachance, J. and Tishkoff, S. A. 2014.  Biased gene conversion skews allele frequencies in human populations, increasing the disease burden of recessive alleles.  The American Journal of Human Genetics, 95(4): 408–420.
  * Li (1987) Li, W.-H. 1987.  Models of nearly neutral mutations with particular implications for nonrandom usage of synonymous codons.  Journal of Molecular Evolution, 24(4): 337–345.
  * Long et al. (2018) Long, H., et al. 2018\.  Evolutionary determinants of genome-wide nucleotide composition.  Nature Ecology & Evolution, 2(2): 237–240.
  * Maruyama (1977) Maruyama, T. 1977.  Stochastic Problems in Population Genetics.  Springer, Berlin Heidelberg.
  * McVean and Charlesworth (1999) McVean, G. A. T. and Charlesworth, B. 1999.  A population genetic model for the evolution of synonymous codon usage: patterns and predictions.  Genetical Research, 74(2): 145–158.
  * Moran (1958) Moran, P. A. P. 1958.  Random processes in genetics.  Mathematical Proceedings of the Cambridge Philosophical Society, 54(1): 60–71.
  * Mugal et al. (2014) Mugal, C. F., Wolf, J. B., and Kaj, I. 2014.  Why time matters: codon evolution and the temporal dynamics of dN/dS.  Molecular Biology and Evolution, 31(1): 212–31.
  * Mugal et al. (2015) Mugal, C. F., Weber, C. C., and Ellegren, H. 2015.  GC-biased gene conversion links the recombination landscape and demography to genomic base composition.  BioEssays, 37(12): 1317–1326.
  * Muyle et al. (2011) Muyle, A., Serres-Giardi, L., Ressayre, A., Escobar, J., and Glémin, S. 2011\.  GC-biased gene conversion and selection affect GC content in the _Oryza_ genus (rice).  Molecular Biology and Evolution, 28(9): 2695–2706.
  * Nagylaki (1983) Nagylaki, T. 1983.  Evolution of a finite population under gene conversion.  Proceedings of the National Academy of Sciences, 80(20): 6278–6281.
  * Nei and Li (1979) Nei, M. and Li, W. H. 1979.  Mathematical model for studying genetic variation in terms of restriction endonucleases.  Proceedings of the National Academy of Sciences, 76(10): 5269–5273.
  * Novak and Barton (2017) Novak, S. and Barton, N. H. 2017.  When does frequency-independent selection maintain genetic variation?  Genetics, 207(2): 653–668.
  * Ohta (1973) Ohta, T. 1973.  Slightly deleterious mutant substitutions in evolution.  Nature, 246(5428): 96–98.
  * Ohta (1976) Ohta, T. 1976.  Role of very slightly deleterious mutations in molecular evolution and polymorphism.  Theoretical Population Biology, 10(3): 254–275.
  * Ohta (1992) Ohta, T. 1992.  The nearly neutral theory of molecular evolution.  Annual Review of Ecology and Systematics, 23(1): 263–286.
  * Peng and Li (2013) Peng, J. and Li, W. V. 2013.  Diffusions with holding and jumping boundary.  Science China Mathematics, 56(1): 161–176.
  * Phillips et al. (2015) Phillips, C., Amigo, J., Carracedo, Á., and Lareu, M. 2015.  Tetra-allelic SNPs: Informative forensic markers compiled from public whole-genome sequence data.  Forensic Science International: Genetics, 19: 100–106.
  * Pontz and Feldman (2020) Pontz, M. and Feldman, M. W. 2020.  Loss of genetic variation in the two-locus multiallelic haploid model.  Theoretical Population Biology, 136: 12–21.
  * Sawyer and Hartl (1992) Sawyer, S. A. and Hartl, D. L. 1992.  Population genetics of polymorphism and divergence.  Genetics, 132(4): 1161–1176.
  * Stoltzfus and Norris (2015) Stoltzfus, A. and Norris, R. W. 2015.  On the causes of evolutionary transition:transversion bias.  Molecular Biology and Evolution, 33(3): 595–602.
  * Vogl and Bergman (2015) Vogl, C. and Bergman, J. 2015.  Inference of directional selection and mutation parameters assuming equilibrium.  Theoretical Population Biology, 106: 71–82.
  * Vogl and Clemente (2012) Vogl, C. and Clemente, F. 2012.  The allele-frequency spectrum in a decoupled moran model with mutation, drift, and directional selection, assuming small mutation rates.  Theoretical Population Biology, 81(3): 197–209.
  * Wright (1931) Wright, S. 1931.  Evolution in mendelian populations.  Genetics, 16(2): 97–159.
  * Wright (1938) Wright, S. 1938.  The distribution of gene frequencies under irreversible mutation.  Proceedings of the National Academy of Sciences, 24(7): 253–259.
  * Wright (1949) Wright, S. 1949.  Adaptation and Selection.  In G. L. Jepsen, E. Mayr, and G. G. Simpson, editors, Genetics, Palaeontology and Evolution, pages 365–389. Princeton University Press, Princeton (NJ).
  * Youssef et al. (2021) Youssef, N., Susko, E., Roger, A. J., and Bielawski, J. P. 2021.  Shifts in amino acid preferences as proteins evolve: A synthesis of experimental and theoretical work.  Protein Science, 30(10): 2009–2028.
  * Zeng (2010) Zeng, K. 2010.  A simple multiallele model and its application to identifying preferred-unpreferred codons using polymorphism data.  Molecular Biology and Evolution, 27(6): 1327–1337.