# Exponential Integral Solutions for Fixation Time
in Wright-Fisher Model With Selection

Vincent Runge Email<EMAIL_ADDRESS>Université Paris-Saclay, CNRS,
Univ Evry, Laboratoire de Mathématiques et Modélisation d’Evry, 91037, Evry-
Courcouronnes, France. Arnaud Liehrmann  Université Paris-Saclay, CNRS, Univ
Evry, Laboratoire de Mathématiques et Modélisation d’Evry, 91037, Evry-
Courcouronnes, France. Université Paris-Saclay, CNRS, INRAE, Université Evry,
Institute of Plant Sciences Paris-Saclay (IPS2), 91190, Gif sur Yvette,
France.  Pauline Spinga  Université Paris-Saclay, CNRS, Univ Evry, Laboratoire
de Mathématiques et Modélisation d’Evry, 91037, Evry-Courcouronnes, France.

###### Abstract

In this work we derive new analytic expressions for fixation time in Wright-
Fisher model with selection. The three standard cases for fixation are
considered: fixation to zero, to one or both. Second order differential
equations for fixation time are obtained by a simplified approach using only
the law of total probability and Taylor expansions. The obtained solutions are
given by a combination of exponential integral functions with elementary
functions. We then state approximate formulas involving only elementary
functions valid for small selection effects. The quality of our results are
explored throughout an extensive simulation study. We show that our results
approximate the discrete problem very accurately even for small population
size (a few hundreds) and large selection coefficients.

Keywords: Wright-Fisher model, fixation time, selection effect, exponential
integral functions

##  1 Introduction

The Wright-Fisher (WF) model is one of the simplest model of genetic drift in
population genetics. It has been developed independently by Sewall Wright [1]
and R.A Fisher [2] in the 1930s. This model describes the dynamic of allele
transmission over time for haploid or diploid individuals based on strict
assumptions: non-overlapping generations, constant population size over time,
random (panmictic) mating.

One of the first applications to the WF model was proposed by Buri et al. [3]
who studied the distribution of gene frequency under selection due to a
mutation in Drosophilia Melanogaster populations. This idea was undertaken
later by Grifford et al. [4] who studied the probability of fixation in
mutated strains of fungus Aspergillus nidulans impacting their viability. More
recent examples among many others include the analysis of genomics data [5, 6,
7, 8, 9, 10], reflecting the relevance of such a simple model to this day.

In addition to genetic drift, the WF model can account for the effects of
other evolutionary forces such as mutation [11, 12, 13], selection [14, 15,
16, 17, 13, 18, 19], the effect of demographic forces including the variation
in population size through time and migration connecting populations [14, 11,
12, 18].

In population genetics, one focuses on different quantities of interest such
as heterozygosity or fixation probabilities [20]. Here we focus on the average
time until fixation (or simply fixation time) in WF model. Many solutions for
fixation time have been found using the diffusion approximation method [21].
The expression of the fixation time for WF model with no additional effect is
well-known [22, 23]. Solutions for fixation of a rare mutant gene (excluding
the case of potential loss) by random drift and with selection [24] or with an
additional mutation rate [25] have also been derived. Moreover, the variance
for fixation with no loss has been found in [26] and [27] with a selection
coefficient.

In this work we intended to contribute to the study of the fixation time for
WF model with selection by stating new analytic expressions. We analyzed the
three possible fixation times: fixation to zero, to one or both. To the best
of our knowledge, only some partial results have been described in literature:
fixation to one for rare mutant for all initial probabilities in complex
integral form [24] and for null initial probability (only) [28] in exponential
integral form. No author has yet highlighted the use of the exponential
integral function for describing its complete solution.

This paper is organized as follows. In Section 2 we present the WF model with
selection effect and the derivation of its associated ordinary differential
equations. In Section 3 we find expressions for fixation time using
exponential integral (EI) functions. In Section 4 we study approximation
results to bypass the use of EI functions for small selection effect. We
eventually provide an extensive simulation study in Section 5 to explore the
validity range of our theoretical results with respect to population size and
selection effect strength.

##  2 The Wright-Fisher Model With Selection

###  2.1 Context and Notations

We consider an initial population size of $N$ haploid individuals and study a
single locus with two possible alleles: $A$ and $a$. This model is equivalent
to the study of diploid individuals who can either be homozygote ($AA$ \-
$aa$) or heterozygote ($Aa$ \- $aA$), considering that $A$ and $a$ are neither
dominant nor recessive. In that case, the dominance parameter of heterozygous
individuals is $h=1/2$ [29, chapter 5.3].

Notice that, in the literature [30, 31], while studying genetic drift, an
effective population size is usually determined. Denoted $N_{e}$, it is
defined as the size of an ideal population that would have the same rate of
loss of genetic diversity - or that would lose heterozygosity at an equal rate
- than the real studied population. The determination of $N_{e}$ is not the
subject of this article, we consider that our population size is equal to the
effective population size ($N~{}=~{}N_{e}$).

We introduce a selection effect $s$ that affects favorably viability of
carriers of allele $A$: their chance of survival is increased by a factor
$1+s$. There are many possibilities to introduce selection effects [20]. We
chose to keep $s$ constant for all generations and we also define coefficient
$\alpha=2sN$ normalized by population size as usually done in the literature
[24, 32, 33, 29] (estimation of $s$ from data in [33]).

Let $x_{t}$ be the number of alleles A in generation $t$,
$x_{t}\in\\{0,\dots,N\\}$, and $y_{t}=x_{t}/N$ its proportion111We chose
notation $y_{t}$ instead of a more natural $p_{t}$ to avoid confusion with
initial allele A frequency, further denoted $p$.: they are realizations of
random variables $X_{t}$ and $Y_{t}$ respectively. The frequency $y_{t+1}$ of
alleles A obtained at time $t+1$ depends on frequency $y_{t}$ at time $t$, but
is moreover skewed by the selection effect, that is, we define a fictional
proportion:

$\pi_{t}=\frac{x_{t}+sx_{t}}{N+sx_{t}}=\frac{(1+s)y_{t}}{1+sy_{t}}\,,$

where the population of allele A is increased by an imaginary quantity
$sx_{t}$. The distribution of alleles in the new generation is then generated
by the binomial distribution $X_{t+1}\backsim\mathcal{B}(N,\pi_{t})$. The
generation $t+1$ represents the offspring of the generation $t$ obtained by
random sampling with replacement. We also define the variation of the
proportion $Y_{t}$ over time by $\Delta P_{t}=Y_{t+1}-Y_{t}$. The transition
probabilities are defined as

$p_{ij}=\mathds{P}(X_{t+1}=j\mid
X_{t}=i)=\binom{N}{j}\Big{(}\frac{(1+s)i}{N+si}\Big{)}^{j}\Big{(}1-\frac{(1+s)i}{N+si}\Big{)}^{N-j}\,,$

for all integers $t$, integers $i$ in $\left\\{1,...,N-1\right\\}$ and $j$ in
$\left\\{0,...,N\right\\}$.

######  Remark 1.

$\left(X_{t}\right)_{t\geq 0}$ is an absorbing Markov chain [34] where
$p_{ij}$ is the probability of transition from state $i$ to $j$ and $X_{t}=0$
and $X_{t}=N$ are the absorbing states.

###  2.2 Fixation Time

The fixation time is defined as the time $T=T_{p}$ needed by the Markov chain
$(X_{t})_{t\geq 0}$ to reach $0$ or $N$ with initial population $X_{0}=pN$ in
$\\{0,\dots,N\\}$, that is:

$T_{p}=\min_{t\geq 0}\\{t\,\,\hbox{such
that}\,\,X_{t}=0\,\,\hbox{or}\,\,X_{t}=N\,\,\hbox{with}\,\,X_{0}=pN\\}\,.$

This is a random variable returning the number of iterations needed to get an
homogeneous population. It is parametrized by its initial allele A frequency
$p$. Our goal is to search for analytic expressions for its expectation when
we consider the limit case for infinite population size $N$, leading to a
continuous variable $p$ in $[0,1]$. We then re-define $T_{p}$ through the
frequencies $Y_{t}$ to remove size-dependency:

$T_{p}=\min_{t\geq 0}\\{t\,\,\hbox{such
that}\,\,Y_{t}=0\,\,\hbox{or}\,\,Y_{t}=1\,\,\hbox{with}\,\,Y_{0}=p\\}\,.$

We decompose its expected value by the law of total expectation and introduce
new notations as follows:

$\displaystyle\mathds{E}[T_{p}]$
$\displaystyle=\mathds{P}(Y_{T_{p}}=0)\mathds{E}[T_{p}\mid
Y_{T_{p}}=0]+\mathds{P}(Y_{T_{p}}=1)\mathds{E}[T_{p}\mid Y_{T_{p}}=1]\,,$
$\displaystyle m(p)$ $\displaystyle=(1-u(p))m_{0}(p)+u(p)m_{1}(p)\,,$

with expectation function $p\mapsto m(p)=\mathds{E}[T_{p}]$, fixation
probability function $p\mapsto u(p)=\mathds{P}(Y_{T_{p}}=1)$, expected
fixation time to zero and one defined as $p\mapsto
m_{1}(p)=\mathds{E}[T_{p}\mid Y_{T_{p}}=1]$ and $p\mapsto
m_{0}(p)=\mathds{E}[T_{p}\mid Y_{T_{p}}=0]$, respectively. Function
$m_{1}(\cdot)$ is often called the fixation time of a mutant gene excluding
the case of potential loss. We also define $M_{1}(p)=u(p)m_{1}(p)$ and
$M_{0}(p)=(1-u(p))m_{0}(p)$. These quantities can be described in terms of
infinite sums:

$\left\\{\begin{aligned}
u(p)&=\mathds{P}(Y_{T_{p}}=1)=\sum_{t=0}^{\infty}\mathds{P}(T_{p}=t\,,\,Y_{T_{p}}=1)\,,\\\
M_{0}(p)&=\sum_{t=0}^{\infty}t\mathds{P}(T_{p}=t\,,\,Y_{T_{p}}=0)\,,\\\
M_{1}(p)&=\sum_{t=0}^{\infty}t\mathds{P}(T_{p}=t\,,\,Y_{T_{p}}=1)\,,\\\
m(p)&=\mathds{E}[T_{p}]=\sum_{t=0}^{\infty}t\mathds{P}(T_{p}=t)=M_{0}(p)+M_{1}(p)\,,\\\
\end{aligned}\right.$

and we can now derive equations binding neighboring values of each function to
itself using the transition probability operator $\mathds{P}_{p\rightarrow
p+\Delta p}=\mathds{P}(Y_{1}=p+\Delta p\mid Y_{0}=p)$, which is the density
probability of the transition in one time step from allele probability $p$ to
allele probability $p+\Delta p$.

######  Proposition 1.

Functions $u(\cdot)$, $M_{0}(\cdot)$, $M_{1}(\cdot)$ and $m(\cdot)$ satisfy
the following equations for initial probability $p$ (sometimes called master
equations) :

$\left\\{\begin{aligned} u(p)&=\int_{\Delta p}\mathds{P}_{p\rightarrow
p+\Delta p}u(p+\Delta p)d(\Delta p)\,,\\\ M_{0}(p)&=1-u(p)+\int_{\Delta
p}\mathds{P}_{p\rightarrow p+\Delta p}M_{0}(p+\Delta p)d(\Delta p)\,,\\\
M_{1}(p)&=u(p)+\int_{\Delta p}\mathds{P}_{p\rightarrow p+\Delta
p}M_{1}(p+\Delta p)d(\Delta p)\,,\\\ m(p)&=1+\int_{\Delta
p}\mathds{P}_{p\rightarrow p+\Delta p}m(p+\Delta p)d(\Delta p)\,,\\\
\end{aligned}\right.$

with integration performed with respect to the Lebesgue measure.

###### Proof.

We first provide details for the last equation. The proof is based on the law
of total probability:

$\mathds{P}(T_{p}=t)=\int_{\Delta p}\mathds{P}_{p\rightarrow p+\Delta
p}\mathds{P}(T_{p}=t\mid Y_{1}=p+\Delta p)d(\Delta p)\,.$

We multiply this equality by $t$ and sum for all integers to get $m(p)$ in
left hand side:

$\displaystyle m(p)$ $\displaystyle=\int_{\Delta p}\mathds{P}_{p\rightarrow
p+\Delta p}\sum_{t=0}^{\infty}t\mathds{P}(T_{p}=t\mid Y_{1}=p+\Delta
p)d(\Delta p)\,,$ $\displaystyle=\int_{\Delta p}\mathds{P}_{p\rightarrow
p+\Delta p}\sum_{t=1}^{\infty}(t-1+1)\mathds{P}(T_{p}=t\mid Y_{1}=p+\Delta
p)d(\Delta p)\,,$ $\displaystyle=\int_{\Delta p}\mathds{P}_{p\rightarrow
p+\Delta p}\sum_{t=1}^{\infty}(t-1)\mathds{P}(T_{p}=t\mid Y_{1}=p+\Delta
p)d(\Delta p)+1\,,$ $\displaystyle=\int_{\Delta p}\mathds{P}_{p\rightarrow
p+\Delta p}\sum_{t=1}^{\infty}(t-1)\mathds{P}(T_{p+\Delta p}=t-1)d(\Delta
p)+1\,,$ $\displaystyle=\int_{\Delta p}\mathds{P}_{p\rightarrow p+\Delta
p}\sum_{t=0}^{\infty}t\mathds{P}(T_{p+\Delta p}=t)d(\Delta p)+1\,,$
$\displaystyle=\int_{\Delta p}\mathds{P}_{p\rightarrow p+\Delta p}m(p+\Delta
p)d(\Delta p)+1\,.$

Results for $M_{0}(\cdot)$ and $M_{1}(\cdot)$ are obtained by considering the
same proof replacing $\mathds{P}(T_{p}=t)$ by
$\mathds{P}(T_{p}=t,Y_{T_{p}}=0)$ or $\mathds{P}(T_{p}=t,Y_{T_{p}}=1)$ in the
law of total probability. For equation in $u(\cdot)$ (the first equation) with
sum all probabilities $\mathds{P}(T_{p}=t,Y_{T_{p}}=1)$ without the
multiplicative factor $t$. ∎

######  Remark 2.

With discrete proportions for the initial probability and a finite size
population ($p=i/N$, $i=0,\dots,N$), equations in Proposition 1 are valid with
integration with respect to counting measures. They can be solved exactly by
matrix inversion.

######  Remark 3.

Similar results are available for the second moment:
$w(p)=\mathds{E}[T^{2}_{p}]=\sum_{t=0}^{\infty}t^{2}\mathds{P}(T_{p}=t)$,
leading to equation

$w(p)=1+2m(p)+\int_{\Delta p}\mathds{P}_{p\rightarrow p+\Delta p}w(p+\Delta
p)d(\Delta p)\,,$

and we would get the variance $v(p)$ as $w(p)-(m(p))^{2}$. We did not succeed
in finding explicit solutions for second moments. Explicit solutions are
known, to the best of our knowledge, only in non-selection mode [26, 27].

###  2.3 Second Order Differential Equations

Master equations can be approximated by Taylor expansions around small
variations in proportions $p$.

######  Proposition 2.

Equations in Proposition 1 can be approximated by the following second order
differential equations:

$\left\\{\begin{aligned} &u^{\prime\prime}(p)+\alpha\,u^{\prime}(p)=0\,,\\\
&M_{0}^{\prime\prime}(p)+\alpha\,M_{0}^{\prime}(p)+\frac{2N}{p(1-p)}(1-u(p))=0\,,\\\
&M_{1}^{\prime\prime}(p)+\alpha\,M_{1}^{\prime}(p)+\frac{2N}{p(1-p)}u(p)=0\,,\\\
&m^{\prime\prime}(p)+\alpha\,m^{\prime}(p)+\frac{2N}{p(1-p)}=0\,,\\\
\end{aligned}\right.$ (1)

with natural boundary conditions: $u(0)=0$, $u(1)=1$, $M_{0}(0)=0$,
$M_{1}(1)=0$ and $m(0)=m(1)=0$, which leads to $M_{0}(1)=(1-u(1))m_{0}(1)=0$
and $M_{1}(0)=u(0)m_{1}(0)=0$.

###### Proof.

A detailed proof is given in A. All the results are based on the following
expansion given here for a generic function $f$:

$\displaystyle\int\mathds{P}_{p\rightarrow p+q}$ $\displaystyle f(p+\Delta
p)d(\Delta p)$ $\displaystyle=\int\mathds{P}_{p\rightarrow p+\Delta
p}\left(f(p)+\frac{f^{\prime}(p)}{1!}\Delta
p+\frac{f^{\prime\prime}(p)}{2!}\Delta p^{2}+O(\Delta p^{3})\right)d(\Delta
p)$ $\displaystyle\approx f(p)+\mathds{E}[\Delta
P]f^{\prime}(p)+\frac{1}{2}\mathds{E}[(\Delta P)^{2}]f^{\prime\prime}(p)\,.$

∎

Solution for unknown function $u(\cdot)$ in first equation is easy to found
(see [21]): it is given by formula $u(p)=\frac{1-e^{-\alpha
p}}{1-e^{-\alpha}}$ and we also have
$1-u(p)=\frac{1-e^{\alpha(1-p)}}{1-e^{\alpha}}$.

##  3 Analytic Solutions for Fixation Time

###  3.1 The Exponential Integral Function

Our solutions are based on the exponential integral function
$x\mapsto\mathrm{Ei}(x)$ from $\mathbb{R}^{*}$ to $\mathbb{R}$ which is a
special function defined by an integral (understood in terms of the Cauchy
principal value):

$\mathrm{Ei}(x)=-\int_{-x}^{+\infty}\frac{e^{-t}}{t}dt=\int_{-\infty}^{x}\frac{e^{t}}{t}dt\,.$

We have the following relation:

$\mathrm{Ei}(x)=\ln|x|e^{x}-\int_{0}^{x}\ln|z|e^{z}dz+\gamma\,,$ (2)

with $\gamma$ the Euler-Mascheroni constant ($\gamma\approx 0.5772156649$) and
integral

$\int e^{-\alpha x}\mathrm{Ei}(\alpha x)dx=\frac{\ln\alpha
x}{\alpha}+\frac{e^{-\alpha x}\mathrm{Ei}(\alpha x)}{\alpha}\,.$ (3)

The shape of the EI function is presented in Figure 1.

Figure 1: Shape of the exponential integral function near zero.

###  3.2 Time to Fixation

We start by solving the differential equation in Proposition 2 for function
$m(\cdot)$ and determine an analytic expression.

######  Theorem 1.

Time to fixation to both zero and one with initial allele frequency $p$ and a
selection effect of order $\alpha=2sN$ is given by formula

$m(p)=2N\Big{(}\omega(\alpha,p)+\omega(-\alpha,1-p)\Big{)}\,,$

with

$\displaystyle\alpha\,\omega(\alpha,p)$
$\displaystyle=u(p)\Big{(}\ln|\alpha|+\gamma-e^{-\alpha}\mathrm{Ei}(\alpha)\Big{)}-\Big{(}\ln|\alpha
p|+\gamma-e^{-\alpha p}\mathrm{Ei}(\alpha p)\Big{)}$
$\displaystyle=\Big{[}\frac{1-e^{-\alpha p}}{1-e^{-\alpha
x}}\Big{(}-e^{-\alpha x}\mathrm{Ei}(\alpha x)+\ln|\alpha
x|+\gamma\Big{)}\Big{]}_{x=p}^{x=1}$ $\displaystyle=\Big{[}\frac{1-e^{-\alpha
p}}{1-e^{-\alpha x}}\Big{(}e^{-\alpha x}\int_{0}^{\alpha
x}\ln|z|e^{z}dz\Big{)}\Big{]}_{x=p}^{x=1}\,,$

###### Proof.

Equation in $m(\cdot)$ in Proposition 1 can be written as

$\left[e^{\alpha p}m^{\prime}(p)\right]^{\prime}=-2N\frac{e^{\alpha
p}}{p(1-p)}=-2N\Big{(}\frac{e^{\alpha
p}}{p}+e^{\alpha}\frac{e^{-\alpha(1-p)}}{1-p}\Big{)}\,,$

which can be integrated using the $\mathrm{Ei}$ function:

$e^{\alpha
p}m^{\prime}(p)=2N\Big{(}-\mathrm{Ei}({\alpha}p)+e^{{\alpha}}\mathrm{Ei}({\alpha}(p-1))\Big{)}+C_{1}\,.$

We obtain by a second integration using 3 (see also [35]):

$m(p)=\frac{2N}{\alpha}\left[e^{-\alpha p}\mathrm{Ei}(\alpha
p)-e^{-\alpha(p-1)}\mathrm{Ei}(\alpha(p-1))+\ln(1-p)-\ln(p)\right]$
$+\,C_{2}-C_{1}\frac{e^{-\alpha p}}{\alpha}\,.$

We determine the constants using boundary values $m(0)=m(1)=0$ and the limit
$\lim\limits_{x\rightarrow 0}\Big{(}\mathrm{Ei}(\alpha
x)-\ln|x|\Big{)}=\gamma+\ln|\alpha|$ and we get:

$\left\\{\begin{aligned}
0&=C_{2}-\frac{C_{1}}{\alpha}+\frac{2N}{\alpha}\Big{(}\gamma+\ln|\alpha|-e^{\alpha}\mathrm{Ei}(-\alpha)\Big{)}\,,\\\
0&=C_{2}-C_{1}\frac{e^{-\alpha}}{\alpha}-\frac{2N}{\alpha}\Big{(}\gamma+\ln|\alpha|-e^{-\alpha}\mathrm{Ei}(\alpha)\Big{)}\,.\\\
\end{aligned}\right.$

This system is easy to solve for constants $C_{1}$ and $C_{2}$ to find the
formula presented in the theorem. The integral expression for
$\omega(\alpha,p)$ comes from relation (2) as the terms involving $\gamma$
cancel out. ∎

######  Remark 4.

Both functions $p\mapsto 2N\omega(\alpha,p)$ and $p\mapsto
2N\omega(-\alpha,1-p)$ are solution to its associated equation in (1) as this
equation is invariant by transformation $(\alpha,p)\mapsto(-\alpha,1-p)$.

###  3.3 Time to Zero and Time to One

Replacing $M_{0}(\cdot)$ by $(1-u(\cdot))m_{0}(\cdot)$ and $M_{1}(\cdot)$ by
$u(\cdot)m_{1}(\cdot)$ in equations (1) we get the following equations:

$\left\\{\begin{aligned}
&m_{0}^{\prime\prime}(p)+\alpha\frac{1+e^{\alpha(1-p)}}{1-e^{\alpha(1-p)}}\,m_{0}^{\prime}(p)+\frac{2N}{p(1-p)}=0\,,\\\
&m_{1}^{\prime\prime}(p)+\alpha\frac{1+e^{-\alpha p}}{1-e^{-\alpha
p}}\,\,m_{1}^{\prime}(p)+\frac{2N}{p(1-p)}=0\,,\\\ \end{aligned}\right.$ (4)

showing the simple relation $m_{0}(p)=m_{1}(1-p)$. Therefore, we only need to
solve the equation for $M_{1}(\cdot)$, as
$M_{0}(p)=(1-u(p))m_{1}(1-p)=\frac{1-u(p)}{u(1-p)}M_{1}(1-p)$ leading to:

$M_{0}(p)=\frac{1-e^{\alpha(1-p)}}{1-e^{-\alpha(1-p)}}\frac{1-e^{-\alpha}}{1-e^{\alpha}}M_{1}(1-p)=e^{-\alpha
p}M_{1}(1-p)\,.$ (5)

We now derive the complete solution for function $m_{1}(\cdot)$ by solving its
differential equation in $M_{1}(\cdot)$ in Proposition 2. This result is
already presented in [24] in integral form.

######  Theorem 2.

Time to fixation excluding allele $A$ loss with initial allele frequency $p$
and a selection effect of order $\alpha=2sN$ is given by formula:

$m_{1}(p)=\frac{2N}{\alpha}\Big{(}e_{1}(p)+e_{2}(p)+e_{3}(p)+\frac{C}{u(p)}+\gamma+\ln|\alpha|-e^{-\alpha}\mathrm{Ei}(\alpha))\Big{)}\,,$

with

$\left\\{\begin{aligned} &e_{1}(p)=\frac{e^{-\alpha p}\mathrm{Ei}(\alpha
p)-e^{\alpha(1-p)}\mathrm{Ei}(\alpha(p-1))}{1-e^{-\alpha p}}\,,\\\
&e_{2}(p)=\frac{\mathrm{Ei}(-\alpha
p)-e^{-\alpha}\mathrm{Ei}(\alpha(1-p))}{1-e^{-\alpha p}}\,,\\\
&e_{3}(p)=\frac{1+e^{-\alpha p}}{1-e^{-\alpha
p}}\ln\Big{(}\frac{1-p}{p}\Big{)}\,,\\\
&C=\frac{e^{-\alpha}}{1-e^{-\alpha}}\Big{(}2(\gamma+\ln|\alpha|)-e^{\alpha}\mathrm{Ei}(-\alpha)-e^{-\alpha}\mathrm{Ei}(\alpha)\Big{)}\,.\end{aligned}\right.$

The proof is straightforward: it is similar to the proof of Theorem 1 by
solving the equation for $M_{1}(\cdot)$ in Proposition 2 and then dividing by
function $u$ to get the result. We can use relation (5) to notice that it is
quite easy to come back to solution of Theorem 1 using relations

$\left\\{\begin{aligned} &u(p)e_{1}(p)+e^{-\alpha p}u(1-p)e_{1}(1-p)=0\,,\\\
&u(p)e_{2}(p)+e^{-\alpha p}u(1-p)e_{2}(1-p)=e^{-\alpha p}\mathrm{Ei}(\alpha
p)-e^{\alpha(1-p)}\mathrm{Ei}(\alpha(p-1))\,,\\\ &u(p)e_{3}(p)+e^{-\alpha
p}u(1-p)e_{3}(1-p)=\ln\Big{(}\frac{1-p}{p}\Big{)}\,.\\\ \end{aligned}\right.$

We plot the shape of functions $m(\cdot)$ and $m_{1}(\cdot)$ for various
values of the selection coefficient in Figure 2. Population size $N$ is only
present as a multiplicative factor for all fixation time formulas. It does not
affect the general shape of these functions.

Figure 2: (A) Time-to-fixation function $p\mapsto m(p)$ and (B) fixation-to-
one function $p\mapsto m_{1}(p)$ for various selection effect values and
arbitrary population size $N=100$.

##  4 An Approximate Solution

In our search for simple solutions, we propose an approximate solution for the
expression of $m(\cdot)$ involving only elementary functions.

######  Theorem 3.

An approximation solution satisfying the boundary conditions is

$m_{(k)}(p)=2N\Big{(}\omega_{k}(\alpha,p)+\omega_{k}(-\alpha,1-p)\Big{)}\,,$

with

$\displaystyle\alpha\omega_{k}(\alpha,p)$
$\displaystyle=\Big{[}\frac{1-e^{-\alpha p}}{1-e^{-\alpha
x}}\Big{(}-\sum_{n=1}^{k}\frac{(-\alpha x)^{n}}{n!}\Big{(}\ln|\alpha
x|-H_{n}\Big{)}\Big{)}\Big{]}_{x=p}^{x=1}\,,$

and $H_{n}$ the n-th harmonic number ($H_{n}=\sum_{i=1}^{n}\frac{1}{i}$). An
upper bound is given by

$|m(p)-m_{(k)}(p)|\leq
2N\frac{e^{\alpha}\alpha^{k}}{k!}\Big{(}\frac{2|\ln\alpha|+e^{-1}}{k+1}+3\Big{)}\,.$

###### Proof.

We simplify the integral term by using the series expansion of the exponential
function.

$e^{-y}\int_{0}^{y}\ln|u|e^{u}du=\int_{0}^{y}\ln|u|e^{u-y}du=\sum_{n=0}^{+\infty}\frac{1}{n!}\int_{0}^{y}\ln|u|(u-y)^{n}du$
$=\sum_{n=0}^{+\infty}\frac{1}{n!}\left(\left[\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{n-k}}{k+1}\right]y^{n+1}\ln|y|-\left[\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{n-k}}{(k+1)^{2}}\right]y^{n+1}\right)$

Using relations
$\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{n-k}}{k+1}=\frac{(-1)^{n}}{n+1}$ and
$\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{n-k}}{(k+1)^{2}}=\frac{(-1)^{n}H_{n+1}}{n+1}$
we get:

$e^{-y}\int_{0}^{y}\ln|u|e^{u}du=-\sum_{n=1}^{+\infty}\frac{(-y)^{n}}{n!}\Big{(}\ln|y|-H_{n}\Big{)}\,.$

The derivation of the upper bound is exposed in B. ∎

######  Remark 5.

By taking the limit when $\alpha$ tends to zero or by considering $m_{(1)}(p)$
we obtain the classic solution $-2N(p\ln(p)+(1-p)\ln(1-p))$ corresponding to
the time to fixation with no selection effect.

#### Quality of the Approximate Solution

We compare the relative distance between the exact solution of the fixation
time to both zero and one (see Theorem 1) with our approximate solution. The
relative distance is given by the following formula:

$D_{k}(p)=\frac{|m(p)-m_{(k)}(p)|}{m(p)}.$

We denote $k_{0.01}$ the minimum number $k$ of elements in $m_{(k)}(p)$ of our
approximate solution such that $D_{k}(p)<0.01$. Figure 3 shows $k_{0.01}$ for
a series of selection effect values (twenty $\alpha$ values evenly spaced
within the interval [1,30] and forty $\alpha$ values evenly spaced within the
interval $[0.1,3]$) and initial probabilities (all values of $p$ within the
interval $[0.05,0.95]$ with a step of size 0.05). We observe that the
$k_{0.01}$ increases with the selection effect value $\alpha$. At $\alpha=0.1$
and marginally to $p$, $k_{0.01}=1$.

In the literature the selection coefficient $\alpha$ is chosen less than $16$
in the simulation study of Kumira’ seminal work [24] ($sN\in[0,8]$). It is
also common in application studies to consider small selection effects
$\alpha<10$ [17]. This corresponds to an approximate solution with up to 40
terms valid at precision $0.01$ (see Figure 3 panel A). It should be noted
that considering large $N$ eventually leads to large $\alpha$. For instance,
in the experimental MS2 bacteriophage data [33] a reasonable set of $N$ is
$1\times 10^{7}$ to $2\times 10^{8}$. The corresponding alpha then ranges from
$0$ to $1\times 10^{9}$. $k_{0.01}$ is computationally intractable for the
upper limit of this range of $\alpha$. However, [33] reports estimation of
this coefficient close to $1$. For $\alpha$ greater than $35-40$ the
approximate and the exact solutions are computationally unstable and some new
ideas has to be developed. Nevertheless, the need for large selection
coefficients $\alpha$ is unclear in the literature as most papers rely on $N$
and $s$ with no understanding of the key role played by the product $sN$.

Figure 3: Heatmap of $k_{0.01}$ calculated on a series of selection effects
$-$ (A) $\alpha\in[1,30]$ and (B) $\alpha\in[0.1,3]$ $-$ and initial
probabilities (see Quality of the Approximate Solution). Note that
$k_{0.01}=2$ is not observed for these series of $\alpha$.

##  5 Simulation Study

###  5.1 Fixation Time Distribution

To give a better understanding of our problem, we have decided to present the
empirical distribution of the fixation time using histograms, separating the
two types of fixation (zero and one). We simulated $10^{5}$ fixation times
with population size $N=10^{4}$ for three chosen selection effect values
($\alpha=1,5,10$) and four chosen initial probabilities ($p=0.1,0.3,0.5,0.7$)
in order to appreciate the shape of the fixation time distribution and its
variation between fixations to zero and one. These results presented in Figure
4 show the impact of the selection effect reducing the probability fixation of
allele loss (see function $1-u(\cdot)$). The exact theoretical distribution is
not known but some results are known with no selection effect [36]. Note that
the distribution of fixation time to one seems Log-symmetric independently of
the chosen initial probability and selection effect value (see C).

Figure 4: Histograms of relative fixation time to zero (in red) and fixation
time to one (in blue). Each histogram represents the result of $10^{4}$
simulations with population size of $N=10^{4}$. The histograms are truncated
after a fixation time of 40000.

###  5.2 Simulations With Respect to The Selection Coefficient, Population
Size and Initial Probability

#### Design of Simulations

We compare the empirical mean of the fixation time to zero, one and both,
obtained by simulations against our exact solution in Theorem 1 (both) and
Theorem 2 (zero or one). We test different population sizes
$(N=100,200,300,1000)$, selection effect values (twenty $\alpha$ values evenly
spaced within the interval $[1,30]$) and initial probabilities (all values of
$p$ within the interval $[0.05,0.95]$ with a step of size $0.05$). For each
combination of these parameters we simulated up to $10^{5}$ repetitions,
interrupting the simulations once we get $10^{4}$ fixations at the desired
state (one, zero or both). For some combinations of $\alpha$ and $p$ for
fixation to zero, the number of repetitions was not enough to reach $10^{4}$
fixations. In that case, we decide to not calculate the empirical mean and set
the corresponding values to NA in Figure 6. This limitation is well understood
by the theoretical probability of fixation to zero (function $1-u(\cdot)$).

#### Metric

The relative distance between the empirical mean of the fixation time to zero,
denoted $\hat{m}_{0}(p)$, and the exact solution $m_{0}(p)$ is given by the
following formula

$d_{0}(p)=\frac{|\hat{m}_{0}(p)-m_{0}(p)|}{\hat{m}_{0}(p)}.$

The same way we define $d_{1}(p)$ and $d(p)$, the relative distance involving
the fixation time to one and both, respectively.

#### Empirical Results

In Figure 5 we observe that for a population size $N=200$ or above, marginally
to the selection effects and initial probabilities, the median of $d_{1}(p)$
is lower than 0.01, i.e 1% of errors. For $d_{0}(p)$ and $d(p)$ the same
statistical control is reached at population size $N=1000$ and $N=300$,
respectively. At population size $N=100$ the distributions are flat which
suggests a strong dispersion of relative distances. The largest observed
relative distance is equal to $0.17$ and is obtained by $d_{0}(p)$ with values
$\alpha=28.47$, $p=0.05$ and $N=100$, leading to absolute values
$\hat{m}_{0}(p)=8.70$ and $m_{0}(p)=10.18$. For this combination of parameters
the exact solution we propose is quite far from the simulation results. Taking
now $\alpha=2.53$, $p=0.05$ and $N=100$, $\hat{m}_{0}(p)$ and $m_{0}(p)$ are
equal to $138.60$ and $138.48$, respectively. The corresponding relative
distance $d_{0}(p)$ is equal to $10^{-3}$ which is a much better estimate.
These first results convinced us to explore the relative distance according to
the range of tested parameters ($\alpha$ and $p$). In Figure 6 (panel A) we
observe that the relative distance is driven by three main effects. (i) The
population size $N$. Consistent with Figure 5, we observe that when $N$
increases, the relative distance decreases. (ii) The selection effect value
$\alpha$. We observe that when $\alpha$ increases, the relative distance
increases. (iii) The side effect. We observe that $d_{0}(p)$ and $d_{1}(p)$
are larger for small and large $p$, respectively. $d(p)$ is a combination of
both $d_{0}(p)$ and $d_{1}(p)$. The side effect correlates well with the
coefficient of variation (see Figure 6 panel B). Hence, the increases of
relative distances in these regions could be explained by an increased
uncertainty on the estimate of $\hat{m}(p)$, $\hat{m}_{0}(p)$ and
$\hat{m}_{1}(p)$.

Figure 5: Histogram of relative distances $d_{0}(p)$, $d_{1}(p)$ and $d(p)$
calculated on a series of population sizes, selection effects and initial
probabilities (see Design of Simulations). When the median (blue vertical
line) of the relative distance distribution is lower than 0.01 (red vertical
line), i.e 1% of errors, the track background color is set to yellow. Figure
6: (A) Heatmap of relative distances $d_{0}(p)$, $d_{1}(p)$ and $d(p)$
calculated on a series of population sizes, selection effects and initial
probabilities (see Design of Simulations). (B) Heatmap of the coefficient of
variation estimated on $10^{4}$ observed fixations for the same series of
parameters. The coefficient of variation is independent of the population size
$N$ (see D).

##  6 Conclusion

Exact and approximate formulas for time to fixation in Wright-Fisher problem
with selection have been found. All exact formulas are derived by only using
the law of total probability and Taylor expansions. The derivation of these
results have then been simplified in comparison to the standard diffusion
approximation method based on the Kolmogorov backward partial differential
equation. Time to fixation to zero, one and both are only expressed in terms
of elementary and exponential integral functions with respect to the initial
allele $A$ proportion $p$ in segment $[0,1]$. For small selection effects an
approximate solution is possible and highlight its link with the standard
solution with no selection. A simulation study have been conducted to
highlight the large validity range of our formulas, even for small population
size. For future developments building estimators for the selection
coefficient based on our exact results would be a possible next step. A better
understanding of the higher moments for time to fixation to zero and one as
the variance would lead to explicit confidence intervals for this model with
selection.

## Contributions and Acknowledgment

Pauline Spinga has been working three months on this project during her
undergraduate internship. She developed an R Shiny app comparing simulation
results to the exact fixation time formula
https://github.com/PaulineSpinga/wrightfisher. She also helped finding the
approximate solution of Theorem 3 and made most of the literature review.
Arnaud Liehrmann was responsible for the simulation study and co-developed its
related R package https://github.com/vrunge/WrightFisherSelection. Vincent
Runge supervised the project, found the theoretical results and wrote the
manuscript. We also would like to thank our colleague Guillem Rigaill (INRAE
researcher and Arnaud’s PhD supervisor) for helpful comments concerning the
simulation study.

## References

  * [1] S. Wright, Evolution in mendelian populations, Genetics 16 (2) (1931) 97.
  * [2] R. A. F. Sir, R. Fisher, The genetical theory of natural selection: a complete variorum edition, Oxford University Press, 1999.
  * [3] P. Buri, Gene frequency in small populations of mutant drosophila, Evolution (1956) 367–402.
  * [4] D. R. Gifford, J. A. G. de Visser, L. M. Wahl, Model and test in a fungus of the probability that beneficial mutations survive drift, Biology letters 9 (1) (2013) 20120310.
  * [5] J. Terhorst, C. Schlötterer, Y. S. Song, Multi-locus analysis of genomic time series data from experimental evolution, PLOS Genetics 11 (4) (2015).
  * [6] M. A. Fares, Experimental evolution and next generation sequencing illuminate the evolutionary trajectories of microbes, in: Advances in the Understanding of Biological Sciences Using Next Generation Sequencing (NGS) Approaches, Springer, 2015, pp. 101–113.
  * [7] P. Tataru, T. Bataillon, A. Hobolth, Inference under a wright-fisher model using an accurate beta approximation, Genetics 201 (3) (2015) 1133–1141.
  * [8] A. Ferrer-Admetlla, C. Leuenberger, J. D. Jensen, D. Wegmann, An approximate markov model for the wright–fisher diffusion and its application to time series data, Genetics 203 (2) (2016) 831–846.
  * [9] T. Taus, A. Futschik, C. Schlötterer, Quantifying selection with pool-seq time series data, Molecular Biology and Evolution 34 (11) (2017) 3023–3034.
  * [10] A. Jonas, T. Taus, C. Kosiol, C. Schlötterer, A. Futschik, Estimating the effective population size from temporal allele frequency changes in experimental evolution, Genetics 204 (2) (2016) 723–735.
  * [11] R. N. Gutenkunst, R. D. Hernandez, S. H. Williamson, C. D. Bustamante, Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data, PLoS Genetics 5 (10) (2009) e1000695.
  * [12] S. Lukić, J. Hey, Demographic inference using spectral methods on SNP data, with an analysis of the human out-of-africa expansion, Genetics 192 (2) (2012) 619–639.
  * [13] M. Steinrücken, A. Bhaskar, Y. S. Song, A novel spectral method for inferring general diploid selection from time series genetic data, The Annals of Applied Statistics 8 (4) (2014).
  * [14] I. Mathieson, G. McVean, Estimating selection coefficients in spatially structured populations from time series data of allele frequencies, Genetics 193 (3) (2013) 973–984.
  * [15] Z. Gompert, Bayesian inference of selection in a heterogeneous environment from genetic time-series data, Molecular Ecology 25 (1) (2015) 121–134.
  * [16] J. P. Bollback, T. L. York, R. Nielsen, Estimation of 2nes from temporal allele frequency data, Genetics 179 (1) (2008) 497–502.
  * [17] A.-S. Malaspinas, O. Malaspinas, S. N. Evans, M. Slatkin, Estimating allele age and selection coefficient from time-serial data, Genetics 192 (2) (2012) 599–607.
  * [18] R. Vitalis, M. Gautier, K. J. Dawson, M. A. Beaumont, Detecting and measuring selection from gene frequency data, Genetics 196 (3) (2014) 799–817.
  * [19] M. Lacerda, C. Seoighe, Population genetics inference for longitudinally-sampled mutants under strong selection, Genetics 198 (3) (2014) 1237–1250.
  * [20] H. Shafiey, D. Waxman, Exact results for the probability and stochastic dynamics of fixation in the wright-fisher model, Journal of theoretical biology 430 (2017) 64–77.
  * [21] M. Kimura, On the probability of fixation of mutant genes in a population, Genetics 47 (6) (1962) 713.
  * [22] A. Etheridge, Diffusion process models in mathematical genetics.
  * [23] T. D. Tran, J. Hofrichter, J. Jost, An introduction to the mathematical structure of the wright–fisher model of population genetics, Theory in Biosciences 132 (2) (2013) 73–82.
  * [24] M. Kimura, T. Ohta, The average number of generations until fixation of a mutant gene in a finite population, Genetics 61 (3) (1969) 763.
  * [25] M. Kimura, Average time until fixation of a mutant allele in a finite population under continued mutation pressure: Studies by analytical, numerical, and pseudo-sampling methods, Proceedings of the National Academy of Sciences 77 (1) (1980) 522–526.
  * [26] P. Narain, A note on the diffusion approximation for the variance of the number of generations until fixation of a neutral mutant gene, Genetics Research 15 (2) (1970) 251–255.
  * [27] T. Maruyama, The average number and the variance of generations at particular gene frequency in the course of fixation of a mutant gene in a finite population, Genetics Research 19 (2) (1972) 109–113.
  * [28] B. Koopmann, J. Müller, A. Tellier, D. Živković, Fisher–wright model with deterministic seed bank and selection, Theoretical Population Biology 114 (2017) 29–39.
  * [29] A. Etheridge, Some Mathematical Models from Population Genetics: École D’Été de Probabilités de Saint-Flour XXXIX-2009, Vol. 2012, Springer Science & Business Media, 2011.
  * [30] S. Braude, A. R. Templeton, Understanding the multiple meanings of ‘inbreeding’and ‘effective size’for genetic management of african rhinoceros populations, African Journal of Ecology 47 (4) (2009) 546–555.
  * [31] J. Wang, E. Santiago, A. Caballero, Prediction and estimation of effective population size, Heredity 117 (4) (2016) 193–206.
  * [32] P. Tataru, M. Simonsen, T. Bataillon, A. Hobolth, Statistical inference in the wright–fisher model using allele frequency data, Systematic biology 66 (1) (2017) e30–e46.
  * [33] J. P. Bollback, T. L. York, R. Nielsen, Estimation of 2 nes from temporal allele frequency data, Genetics 179 (1) (2008) 497–502.
  * [34] J. G. Kemeny, J. Snell, Laurie finite markov chains springer verlag, New York (1976).
  * [35] M. Geller, E. W. Ng, A table of integrals of the exponential integral, Journal of Research of the National Bureau of Standards 71 (1969) 1–20.
  * [36] M. Kimura, The length of time required for a selectively neutral mutant to reach fixation through random frequency drift in a finite population, Genetics Research 15 (1) (1970) 131–133.

##  Appendix A Differential Equations for Fixation Time

We describe the proof for the equation in $m(\cdot)$. All other results are
proven the same way (we only have to change the constant term). Using a Taylor
expansion, we get

$\displaystyle m(p)$ $\displaystyle=1+\int_{\Delta p}\mathds{P}_{p\rightarrow
p+q}m(p+\Delta p)d(\Delta p)$ (6)
$\displaystyle=1+\int\mathds{P}_{p\rightarrow p+\Delta
p}\left(m(p)+\frac{m^{\prime}(p)}{1!}\Delta
p+\frac{m^{\prime\prime}(p)}{2!}\Delta p^{2}+O(\Delta p^{3})\right)d(\Delta
p)$ $\displaystyle\approx 1+m(p)+\mathds{E}[\Delta
P]m^{\prime}(p)+\frac{1}{2}\mathds{E}[(\Delta P)^{2}]m^{\prime\prime}(p)\,.$

With $\Delta P=(Y_{t+1}-Y_{t})\left.\right|_{Y_{t}=p}$, we have:

$\mathds{E}[Y_{t}]=p\,,\quad\mathds{E}[Y_{t+1}]=\frac{(1+s)p}{1+sp}\,,\quad\mathds{E}[\Delta
P]=\frac{sp(1-p)}{1+ps}\,,$

and

$\operatorname*{\rm{Var}}(\Delta
P)=\operatorname*{\rm{Var}}(Y_{t+1})=\frac{1}{N}\frac{(1+s)p}{1+sp}(1-\frac{(1+s)p}{1+sp})=\frac{1}{N}\frac{(1+s)p(1-p)}{(1+sp)^{2}}\,,$
$\mathds{E}[(\Delta P)^{2}]=\operatorname*{\rm{Var}}(\Delta
P)+\mathds{E}[\Delta
P]^{2}=\frac{p(1-p)}{(1+ps)^{2}}\left[\frac{1+s}{N}+s^{2}p(1-p)\right]\,.$

We derive approximate formulations for these two first moments by Taylor
expansion in power of $\alpha/2N$ using relation $s=\frac{\alpha}{2N}$. We
get:

$\mathds{E}[\Delta
P]=p(1-p)\frac{\frac{\alpha}{2N}}{1+p\frac{\alpha}{2N}}\simeq
p(1-p)\frac{\alpha}{2N}(1-p\frac{\alpha}{2N})=p(1-p)\frac{\alpha}{2N}+o\Big{(}\frac{\alpha}{2N}\Big{)}\,.$

Similarly,

$\displaystyle\mathds{E}[(\Delta P)^{2}]$
$\displaystyle=\frac{p(1-p)}{N}\frac{(1+\frac{\alpha}{2N})}{(1+p\frac{\alpha}{2N})^{2}}+o\Big{(}\frac{\alpha}{2N}\Big{)}\,$
$\displaystyle=\frac{p(1-p)}{N}(1+\frac{\alpha}{2N})(1-2p\frac{\alpha}{2N}))+o\Big{(}\frac{\alpha}{2N}\Big{)}=\frac{p(1-p)}{N}+o\Big{(}\frac{\alpha}{2N}\Big{)}\,.$

Removing the small order terms leads to the desired differential equation.

##  Appendix B Approximate Solution Upper Bound

We recall that $\alpha$ is positive.

$\frac{\alpha}{2N}|m(p)-m_{(k)}(p)|=\Big{|}\Big{[}\frac{1-e^{-\alpha
p}}{1-e^{-\alpha x}}\Big{(}-\sum_{n=k+1}^{+\infty}\frac{(-\alpha
x)^{n}}{n!}\Big{(}\ln|\alpha x|-H_{n}\Big{)}\Big{)}\Big{]}_{x=p}^{x=1}$
$-\Big{[}\frac{1-e^{\alpha(1-p)}}{1-e^{\alpha
x}}\Big{(}-\sum_{n=k+1}^{+\infty}\frac{(\alpha x)^{n}}{n!}\Big{(}\ln|\alpha
x|-H_{n}\Big{)}\Big{)}\Big{]}_{x=1-p}^{x=1}\Big{|}$
$\leq\sum_{n=k+1}^{+\infty}\frac{\alpha^{n}}{n!}\Big{(}|\ln\alpha|+H_{n}\Big{)}+\sum_{n=k+1}^{+\infty}\frac{(\alpha
p)^{n}}{n!}\Big{(}|\ln\alpha p|+H_{n}\Big{)}$
$+(1-p)\sum_{n=k+1}^{+\infty}\frac{\alpha^{n}}{n!}\Big{(}|\ln\alpha|+H_{n}\Big{)}+\sum_{n=k+1}^{+\infty}\frac{(\alpha(1-p))^{n}}{n!}\Big{(}|\ln\alpha(1-p)|+H_{n}\Big{)}\,,$

applying triangle inequalities and using relations $\frac{1-e^{-\alpha
p}}{1-e^{-\alpha}}\leq 1$, $\frac{1-e^{\alpha(1-p)}}{1-e^{\alpha}}\leq 1-p$.
With the loose upper bound $H_{n}\leq n$ and
$\sum_{n=k+1}^{+\infty}\frac{z^{n}}{n!}\leq e^{z}\frac{z^{k+1}}{(k+1)!}$ for
positive $z$, we get:

$\frac{1}{2N}|m(p)-m_{(k)}(p)|\leq\alpha^{k}\Big{[}\frac{e^{\alpha}}{k!}\Big{(}\frac{|\ln\alpha|}{k+1}+1\Big{)}+$
$\frac{e^{\alpha p}p^{k+1}}{k!}\Big{(}\frac{|\ln\alpha
p|}{k+1}+1\Big{)}+\frac{e^{\alpha(1-p)}(1-p)^{k+1}}{k!}\Big{(}\frac{|\ln\alpha(1-p)|}{k+1}+1\Big{)}\Big{]}$
$\leq\frac{e^{\alpha}\alpha^{k}}{k!}\Big{(}\frac{|\ln\alpha|+p^{k+1}|\ln\alpha
p|+(1-p)^{k+1}|\ln\alpha(1-p)|}{k+1}+3\Big{)}\,,$

using the upper bound $e^{\alpha}$ for both $e^{\alpha p}$ and
$e^{\alpha(1-p)}$. Eventually, we have:

$p^{k+1}|\ln p|\leq pp^{k}|\ln p|\leq\frac{p}{Me}\leq pe^{-1}\,,$

and we get

$|m(p)-m_{(k)}(p)|\leq
2N\frac{e^{\alpha}\alpha^{k}}{k!}\Big{(}\frac{2|\ln\alpha|+e^{-1}}{k+1}+3\Big{)}\,.$

##  Appendix C Fixation time distribution on the $\log$ scale.

Figure 1: The distribution of fixation time to one seems Log-symmetric.
Histograms of relative fixation time to zero (in red) and fixation time to one
(in blue) on the log scale. Each histogram represents the result of $10^{4}$
simulations with population size of $N=10^{4}$.

##  Appendix D Coefficient of variation is independent of the population size
$N$

Figure 2: The coefficient of variation is independent of the population size
N. Heatmap of the coefficient of variation estimated on $10^{4}$ observed
fixations for a series of population sizes, selection effects and initial
probabilities.


[affiliation=1]Hye-jinShim [affiliation=2]MdSahidullah [affiliation=1]Jee-
weonJung [affiliation=1]ShinjiWatanabe [affiliation=3]TomiKinnunen

# Beyond Silence:
Bias Analysis through Loss and Asymmetric Approach in Audio Anti-Spoofing

###### Abstract

Current trends in audio anti-spoofing detection research strive to improve
models' ability to generalize across unseen attacks by learning to identify a
variety of spoofing artifacts. This emphasis has primarily focused on the
spoof class. Recently, several studies have noted that the distribution of
silence differs between the two classes, which can serve as a shortcut. In
this paper, we extend class-wise interpretations beyond silence. We employ
loss analysis and asymmetric methodologies to move away from traditional
attack-focused and result-oriented evaluations towards a deeper examination of
model behaviors. Our investigations highlight the significant differences in
training dynamics between the two classes, emphasizing the need for future
research to focus on robust modeling of the bonafide class.

###### keywords:

anti-spoofing, deepfake detection, spoofing detection, shortcut learning,
ASVspoof

##  1 Introduction

Recent progress in voice conversion (VC) and text-to-speech (TTS) technologies
have intensified concerns regarding their potential for malicious use,
emphasizing the critical role of audio anti-spoofing systems. Audio spoofing
detection systems, as binary classifiers, distinguish genuine human speech
(bonafide) from artificially generated (spoofed) speech. The main direction in
this field is toward advancing the systems' ability to generalize across
unseen spoofing attacks by focusing on learning diverse spoofing artifacts. To
this end, a pivotal shift from traditional hand-crafted features [1, 2] to
data-driven approaches [3, 4, 5] and data augmentation [6, 7, 8] have been
pursued. Research on diversifying training data [9, 10, 11], incorporating
domain adaptation strategies [12, 13, 14], and leveraging large-scale pre-
trained models [15, 16, 17] also follow this trend.

Despite these efforts, the evolution of TTS and VC systems, especially those
utilizing state-of-the-art methods like diffusion, are likely to involve fewer
spoofing artifacts and create subtle variations in them without hard effort.
This suggests that learning diverse spoofing artifacts does not guarantee the
detection of emerging unknown attacks. To tackle this, several studies have
explored methods for modeling robust bonafide features and distinguishing them
from spoof features within the latent space [18, 19, 20, 21]. Class-wise
analyses have sought to determine the differences between bonafide and spoof
classes as detailed in [22, 23, 20] and their findings primarily converge on
the impact of `silence', which has a spurious correlation with spoofing
detection.

Figure 1: Comparison of training loss. The left and right figures illustrate
bonafide and spoof classes. x-axis and y-axis indicate training epochs and
loss magnitude. Regardless of the implementation of data augmentation, the two
class losses differ on a large scale.

These kinds of external factors that could unintentionally lead the biased
model predictions are known as shortcuts [24], recently getting a lot of
attention in the broader deep learning literature. These shortcuts challenge
the ability to determine if a model genuinely distinguishes between classes or
relies on irrelevant cues associated with class labels. In the context of
audio anti-spoofing, `silence' is a well-known data bias; for instance, the
model trained with ASVspoof2017 dataset [25] is vulnerable to silence [26] and
the ASVspoof2019 dataset [27] exhibits an unequal distribution of silence
between spoofed and bonafide samples [28, 29]. The impact of silence varies
not only with the data source but also with its distinct types/patterns
introduced during data collection. Efforts to understand and mitigate the
influence of silence, including the analysis of silence trimming and its
effects, have been undertaken [28, 20]. Other than silence, unveiled factors
also have been studied in [30], recently.

While in-depth analyses in audio anti-spoofing have significantly advanced the
field and provided invaluable insights, several areas remain for further
explorations: (i) the majority of studies have focused on spoof class,
specifically attack-wise interpretations; (ii) class-wise analyses exist, yet
the emphasis has primarily been on silence, except for [30] exploring other
data bias factors; (iii) existing analyses predominantly rely on evaluating
outcomes, such as performance metrics and score distributions, without a
thorough examination of the internal workings of models during training. Our
research diverges from existing studies by delving into the training process
through loss analysis and adopting a novel asymmetric approach for
class/phase-wise interpretation.

Our analysis starts with examining the loss curves for bonafide and spoof
classes separately, both with and without the application of RawBoost [7] data
augmentation, as illustrated in Figure 1. Note that illustrated loss values
are raw values. For the model update, a loss weight of 0.9 (for bonafide) and
0.1 (for spoof) was employed, considering the imbalanced number of samples in
each class. Interestingly, the results indicate that the loss associated with
the bonafide class is significantly lower than that for the spoof class. Even
if one considers the loss weights in the training phase, a substantial gap
still remains between the two classes. These results could signify that the
bonafide class is inherently easier to train than the spoof class.
Alternatively, this might suggest that modeling the bona fide class is not
trivial; however, there exists a shortcut that significantly reduces the
training loss. We deploy various loss functions to reveal the meaning behind
the low magnitude of the bonafide loss, between the two scenarios mentioned.
In particular, the objective functions that can consider the difficulty of
samples enable us to understand how the model deals with each class based on
our findings.

Furthermore, our approach involves training and testing the model with an
asymmetric intervention to assess the bonafide and spoof classes separately.
Unlike the methodology in [30], which applies interventions across both
classes and phases for a comprehensive model-level interpretation, our
strategy focuses interventions on one side only. This allows a more precise
understanding of how different phases and class-specific traits affect model
performance. Additionally, to mitigate potential concerns regarding the
influence of silence, we also conduct our analysis combined with silence
trimming as well, demonstrating that our findings are not biased by the
presence of silence. Our findings pave the way for new directions of future
research, shifting focus from the currently predominant studies centered
around the spoof class.

##  2 Method

In this section, we introduce two primary methodologies that facilitate our
in-depth investigations of anti-spoofing systems: (i) a loss-based analysis
that contrasts objective functions, distinguishing between those that
prioritize either easier or more challenging samples, and (ii) an asymmetric
intervention analysis that examines the impact of intervention on phases
(train or test) and classes (bonafide or spoof).

###  2.1 Loss based analysis

Previous analyses of silence within audio anti-spoofing contexts have paved
the way for separate examinations of bonafide and spoof classes. This study
extends these insights by exploring the various objective functions that can
help us understand the model's behavior in each class; loss analysis is widely
adopted in machine learning for understanding the model behavior during the
training [31]. Furthermore, the loss has been directly employed in the image
domain to infer data bias in [32, 33, 34] by analyzing its correlation with
sample difficulty and data bias.

In this work, we conduct two types of analysis using loss functions.
Initially, we calculate the loss for bonafide and spoof classes separately
during training to observe their convergence patterns and magnitude
differences. Significant differences in these areas may indicate model bias,
as shown in Figure 1. Subsequently, we differentiate between two categories of
loss functions for assessing model behavior: prioritizing hard or easy
samples. Among diverse loss functions utilized in hard negative mining, we
deliberately choose a few that dynamically modulate each sample's impact on
the overall loss, rather than explicitly selecting samples. This strategy
enables an in-depth exploration of the model's inherent reactions to varying
sample types.

FocalLoss [35] prioritizes hard samples by modifying the categorical cross-
entropy loss. It amplifies the loss for samples with inaccurate model
predictions or where the model exhibits uncertainty (i.e., low predicted
probability for the correct class), thus prioritizing hard or incorrectly
classified samples without resorting to specific sample selection strategies.

SuperLoss [36] assigns a weight to each sample based on a moving average of
its past losses, targeting potentially noisy or outlier samples. This method
progressively focuses on challenging or ambiguous samples.

CurricularFace [37] emphasizes more challenging samples by increasing their
relative loss compared to easier ones using class-specific margins. It adjusts
the target margins for classes, making it easier or harder for the model to
classify them correctly as training progresses.

Generalized cross entropy (GCE) [38], in contrast, focuses on easier samples.
It enhances the penalty for misclassifying classes with lower probabilities,
thereby prioritizing minority classes. This method can be particularly useful
in audio anti-spoofing to direct the model's attention towards bonafide
samples. In addition, GCE is employed in model debiasing research to
concentrate on biased samples, aligning with empirical observations that such
samples exhibit lower loss values during training.

###  2.2 Asymmetric intervention analysis

Our analysis follows an interventional approach proposed recently in [30] to
reveal shortcuts in speech anti-spoofing beyond silence. The essence of this
approach is to intentionally modify (intervene) an existing dataset to provoke
the classifier to rely on shortcuts. We start by reviewing the approach and
then explain how we modify it in this study.

To explore potential shortcuts in audio anti-spoofing, [30] employed standard
audio manipulations (e.g., MP3 compression or additive noise) with randomized
parameters applied to datasets to create artificial statistical associations
between the audio and the class label. For example, applying MP3 compression
to only bonafide data in both training and test data while leaving spoof data
unaltered resulted in equal error rates (EERs) of 0%. Conversely, applying MP3
compression exclusively to spoofed test data while keeping all other
conditions the same led to an opposite outcome (EER $>99\%$), illustrating a
complete label flip. Such extreme interventions reveal the vulnerability of
the spoofing detection system against potential bias factors unrelated to
spoofing artifacts residing within the data.

Table 1: Five experimental configurations using eight subsets of a dataset. O
is the original (unintervened) configuration. All other four sets have one
particular subset intervened.

Intervened phase | Configuration | Train set | Test set
---|---|---|---
- | O | $\mathcal{D}_{\text{trn, bona}}\cup\mathcal{D}_{\text{trn,spf}}$ | $\mathcal{D}_{\text{test, bona}}\cup\mathcal{D}_{\text{test, spf}}$
Train |  Tr_B | $\mathcal{D}^{\textbf{intervened}}_{\text{trn, bona}}\cup\mathcal{D}_{\text{trn,spf}}$ | $\mathcal{D}_{\text{test, bona}}\cup\mathcal{D}_{\text{test, spf}}$
Tr_S | $\mathcal{D}_{\text{trn, bona}}\cup\mathcal{D}^{\textbf{intervened}}_{\text{trn,spf}}$ | $\mathcal{D}_{\text{test, bona}}\cup\mathcal{D}_{\text{test, spf}}$
Test |  Te_B | $\mathcal{D}_{\text{trn, bona}}\cup\mathcal{D}_{\text{trn,spf}}$ | $\mathcal{D}^{\textbf{intervened}}_{\text{test, bona}}\cup\mathcal{D}_{\text{test, spf}}$
Te_S | $\mathcal{D}_{\text{trn, bona}}\cup\mathcal{D}_{\text{trn,spf}}$ | $\mathcal{D}_{\text{test, bona}}\cup\mathcal{D}^{\textbf{intervened}}_{\text{test, spf}}$

The previous approach concurrently applied interventions during both training
and test phases to understand the overall model's behaviors under different
interventions, producing the boundary cases of `near-perfect' outcomes and
`worse than coin flip' (label flip) results. In contrast, our investigation
takes a class-wise interpretive approach based on the understanding that each
class responds differently. Consequently, we strategically apply interventions
at only one phase or one class at a time, leveraging two binary dimensions:
phase (training or testing) and class (bonafide or spoof). Our approach
enables (i) to reveal the interventions associated with each class separately,
(ii) to evaluate the robustness of the model's representation for each class
separately, and (iii) to compare the robustness of class modeling by observing
the outcomes of intervention.

In particular, our methodology of employing class- and phase-wise intervention
in either class in either phase results in four distinct intervention
configurations. Formally, let
$\mathcal{D}=\\{(x_{i},y_{i})|x_{i}\in\mathbb{X},y_{i}\in\mathbb{Y},\text{for}\
i=1,2,...,N$, where $\mathcal{D}$ represents the dataset, ($x_{i},y_{i}$)
denotes an individual sample, and $N$ is the number of samples. Typically, the
corpus is divided into two subsets: the training ($\mathcal{D}_{\text{trn}}$)
and the test ($\mathcal{D}_{\text{test}}$) sets. We introduce another
categorization based on class, dividing the dataset into four subsets:
$\mathcal{D}_{\text{trn, bona}}$, $\mathcal{D}_{\text{trn, spf}}$,
$\mathcal{D}_{\text{test, bona}}$, $\mathcal{D}_{\text{test, spf}}$.
Furthermore, we generate four additional subsets with applied interventions,
where the labels for each sample remain:
$\mathcal{D}^{\textbf{intervened}}_{\text{trn, bona}}$,
$\mathcal{D}^{\textbf{intervened}}_{\text{trn, spf}}$,
$\mathcal{D}^{\textbf{intervened}}_{\text{test, bona}}$,
$\mathcal{D}^{\textbf{intervened}}_{\text{test, spf}}$. Note that the
intervention does not affect the ground-truth label $\mathbb{Y}$. Composing
the eight total subsets, we present five experimental configurations in Table
1: O, $\textbf{Tr}\\_\textbf{B}$, $\textbf{Tr}\\_\textbf{S}$,
$\textbf{Te}\\_\textbf{B}$, and $\textbf{Te}\\_\textbf{S}$. Therefore, the
comparison of B and S, modification on the class side, enables us to analyze
the class-wise effect of interventions. Comparison of Tr and Te, on the other
hand, helps to understand the class-wise difference and how the model operates
differently depending on whether the model learns such intervened condition.

##  3 Experimental setup

###  3.1 Datasets

ASVspoof2019 [27] is a dataset dedicated to developing audio anti-spoofing and
spoofing-robust automatic speaker verification systems. It includes 6 and 13
types of spoofing attacks in train/dev and test. The number of utterances for
train/dev and test is $50,224$ and $71,237$, respectively.

Table 2: Performance comparison on different loss functions. Apart from
original categorical cross entropy (CCE), FocalLoss, SuperLoss, and
CurricularFace are designed to concentrate on hard samples based on empirical
loss while generalized cross entropy (GCE) focuses on easier samples. All
evaluations are conducted on the ASVspoof2019 LA dataset.

Loss function |  | EER (%)
---|---|---
Original (CCE) |  | 1.39
FocalLoss[35] |  | 1.67
SuperLoss[36] |  | 1.45
CurricularFace[37] |  | 1.53
GCE[38] |  | 1.35

ASVspoof2021 [39] includes the latest spoofing attacks, which consists of
different test sets, each for LA and DF scenarios. The attacks in the LA
subset are processed with real telephony systems introducing artifacts related
to encoding and transmission. The attacks in the DF subset includes various
audio compressions.

###  3.2 Implementation details

Model architectures used in this paper are AASIST, AASIST-L [40] models with
data augmentation and multi-dataset co-training. Those two models are state-
of-the-art models that directly operate on raw waveform and they only differ
in the number of parameters.

Table 3: Performance comparison when the model is tested on asymmetric ways.
All models are trained on the original training dataset but evaluated in
different configurations. Here, the ratio indicates a relative change of Te_B
compared to that of Te_S.

System | O | Te_S | Te_B | Ratio
---|---|---|---|---
MP3 |  |  |  |
AASIST | 0.83 | 0.77 | 9.07 | 137.33
AASIST-L | 1.39 | 1.21 | 8.78 | 41.06
AASIST-L w/ RawBoost | 1.59 | 1.62 | 3.96 | 79.00
AASIST-L w/ MDL | 0.99 | 1.29 | 3.63 | 8.80
Noise |  |  |  |
AASIST | 0.83 | 0.39 | 31.08 | 68.75
AASIST-L | 1.39 | 0.65 | 35.50 | 46.09
AASIST-L w/ RawBoost | 1.59 | 1.06 | 13.27 | 22.04
AASIST-L w/ MDL | 0.99 | 0.69 | 8.24 | 24.17
Loudness |  |  |  |
AASIST | 0.83 | 1.18 | 5.33 | 12.86
AASIST-L | 1.39 | 2.23 | 6.73 | 6.36
AASIST-L w/ RawBoost | 1.59 | 1.52 | 2.38 | 11.29
AASIST-L w/ MDL | 0.99 | 1.44 | 4.09 | 6.89

Intervention types are selected among five different interventions in [30]. We
employ three interventions: _MP3 compression_ , _additive white noise_ , and
_loudness normalization_. Those interventions are considered since MP3
compression and white noise are discovered as the most influential ones, while
loudness normalization is the least effective one.

Data augmentation is implemented to check the difference when we employ the
model considered more robust. We utilize RawBoost [7] which includes three
different augmentation techniques: linear and non-linear convolutive noise,
multi-band filters, and Hammerstein systems [41]. We deploy three of them
simultaneously as it showed the best result in [7].

Table 4: Performance comparison in EER when the model is tested on asymmetric
ways with silence trimming. Silence trimming is applied in both the training
and the test phases; we additionally report the results when the silence is
only removed from either phase.

System | Silence trimming | O | MP3 | Noise | Loudness
---|---|---|---|---|---
Te_S | Te_B | Ratio | Te_S | Te_B | Ratio | Te_S | Te_B | Ratio
AASIST-L | - | 1.39 | 1.21 | 8.78 | 41.06 | 0.65 | 35.50 | 46.09 | 2.23 | 6.73 | 6.36
AASIST-L w/ RawBoost | - | 1.59 | 1.62 | 3.96 | 79.00 | 1.06 | 13.27 | 22.04 | 1.52 | 2.38 | 11.29
AASIST-L | train | 35.05 | 28.36 | 39.32 | 0.64 | 47.09 | 28.14 | 0.57 | 27.44 | 47.28 | 1.61
| test | 25.14 | 20.07 | 37.7 | 2.48 | 8.27 | 59.04 | 2.01 | 13.92 | 40.83 | 1.40
| train&test | 18.65 | 17.52 | 30.65 | 10.62 | 17.24 | 34.49 | 11.23 | 14.79 | 26.99 | 2.16
AASIST-L w/ RawBoost | train | 45.52 | 40.65 | 48.47 | 0.61 | 66.93 | 31.14 | 0.47 | 36.10 | 55.35 | 1.04
| test | 27.79 | 29.5 | 30 | 1.29 | 30.74 | 26.85 | 0.32 | 21.82 | 37.05 | 1.55
| train&test | 19.73 | 17.56 | 31.69 | 5.51 | 21.07 | 32.35 | 9.42 | 12.75 | 29.73 | 1.43

Multi-dataset trained model with sharpness optimization [11] is utilized to
investigate the robust model similar to data augmentation. To further the
enhance generalization capability, the model is trained using multiple
datasets at the same time and optimized by sharpness-aware optimization [42].
We select the model trained by both ASVspoof2015 and ASVspoof2019 LA with
adaptive sharpness-aware minimization (ASAM) [43] as it showed the best
performance in ASVspoof2019 LA evaluation by 0.99% of EER.

Silence trimming is implemented to mitigate the influence of silence that
might distort the imbalance results. Our silence trimming algorithm works as
follows. First, we detect the speech frames using a simple energy-based
algorithm as described in [44]. We have used a frame size of $25$ ms with an
$8$ms shift. Then we remove the silences where the silence length is more than
$50$ ms.

Table 5: Performance comparison in EER when the model trained on asymmetric
ways. All evaluations are conducted using original test set without
interventions. The bold type indicates the case that spoof side intervention
lower EER than the bonafide one.

System |  | Intervention
---
type
Config. |  | 2019
---
LA
| 2021
---
LA
| 2021
---
DF
Original | - | - | 1.39 | 12.18 | 21.8
RawBoost | - | - | 1.59 | 6.42 | 17.48
RawBoost | MP3 | Tr_B | 16.21 | 15.45 | 30.86
Tr_S | 3.68 | 17.81 | 23.1
Noise | Tr_B | 71.02 | 62.45 | 62.66
Tr_S | 15.48 | 24.12 | 30.06
Loudness | Tr_B | 13.05 | 10.66 | 18.7
Tr_S | 1.82 | 9.49 | 17.97

##  4 Results and Analysis

###  4.1 Comparison of class-wise loss and different objective functions

As previously introduced in Section 1 with Figure 1, the observed loss curve
reveals an unexpected pattern, especially considering the predominance of
spoof data in the training set; in general, more frequently represented spoof
classes are expected to have lower losses, consistent with the general bias of
neural networks to focus on more frequently represented classes. The results
indicate that effective training of the bonafide class is hindered not only by
its fewer amount of samples in the training dataset but also by its
comparatively lower loss. Consequently, the neural network would prioritize
minimizing the spoof class's loss in which the loss scale is much higher,
leading to an inherent bias, which is unintended. This research marks the
first to identify this particular bias within the context of audio spoofing
detection, shedding light on the challenges posed by class imbalance and its
impact on model training dynamics.

Building upon these insights, we further examine the potential class bias in
model training towards the spoof class using diverse loss functions outlined
in Section 2.1. Our comparison, shown in Table 2, involves four different loss
functions – three (FocalLoss [35], SuperLoss [36], and CurricularFace [37])
designed for challenging (spoof) samples and one (GCE[38] for easier
(bonafide) samples. The results reveal a decline in performance with the three
loss functions aimed at the spoof class, while the GCE, which prefers the
bonafide class, slightly improves model performance. Despite adjustments using
tunable parameters in FocalLoss, SuperLoss, and CurricularFace, all outcomes
showed no significant enhancement111Presented performances are the best
results among our experiments.. Results highlight the training bias towards
the spoof class. This pattern of bias and its impact on performance encourages
us to call for a significant shift in focus for future audio anti-spoofing
efforts, emphasizing the importance of robust modeling the bonafide class over
solely concentrating on detecting spoofing artifacts. In addition, the field
of anomaly sound detection [45] further demonstrates the advantages of
focusing on bonafide modeling. State-of-the-art systems prioritize modeling
the normal sound class and effectively identify deviations significantly
distant from these norms as anomalies.

###  4.2 Asymmetric test results

Table 3 presents the results of applying asymmetric interventions to a
specific class during the test phase. These results are analyzed from two
angles: (i) evaluating the intervention's impact by comparing the original
results (O) with those after interventions on the spoof (Te_S) or bonafide
(Te_B) classes (columns 2 and 3), and (ii) examining how interventions
differently affect the bonafide class compared to the spoof class, as
indicated in column 4. From the first angle, enhancements in performance post-
intervention hint at possible model biases due to shortcuts. Conversely,
declines indicate the intervention's irrelevance to the target class or might
reflect domain differences or unknown factors. The ``Ratio'' in column 4,
calculated as the relative performance impact
$(|\textbf{O}-\textbf{Te}\\_\textbf{B}|/\textbf{O})/(|\textbf{O}-\textbf{Te}\\_\textbf{S}|/\textbf{O})$,
offers deeper insight into which class the intervention impacts more. Here,
$\text{Ratio}>1$ refers to a greater influence on the bonafide class, while
$\text{Ratio}<1$ implies spoof class is more affected.

Our findings highlight two main observations that align with our analysis of
loss functions. First, interventions are prone to improve Te_S (e.g., AASIST
with MP3 or AASIST-L with Noise on Te_S), while consistently lowering Te_B
performance. This implies that factors leading to quality degradation in
utterances tend to be linked with the spoof class, helping the model to
classify such inputs as spoofed. The drop in Te_B performance may result from
interventions causing bonafide utterances to resemble spoofed ones more
closely. Second, the ``Ratio'' always surpasses 1, as a universal trend
regardless of which intervention applied. This emphasizes the bonafide class's
modeling fragility, suggesting that future research should more focus on the
robust modeling of the bonafide class.

###  4.3 Further test interventions with silence trimming

As introduced in Section 1, the silence could significantly influence the
class-wise imbalance results. This prompts an important question: ``Do our
observations remain valid when silence is excluded from the training or
testing datasets or both?" To address this, we conducted experiments with
silence removed – referred to as silence trimming – to assess its effect on
our findings and show the results in Table 4. The results show a remarkable
consistency across most of our results. When silence was eliminated from both
phases, effectively eliminating its influence, the ratio consistently exceeded
1, reaching a peak of 46.09. In additional experiments where silence was
removed from only one phase, there were a few instances where the ratio fell
below 1; however, the majority still exhibited ratios significantly greater
than 1. These consistent outcomes robustly support the conclusion that the
presence of silence does not skew our findings.

###  4.4 Asymmetric training results

Lastly, we focus on asymmetric interventions during the training phase to
explore two key questions: (i) ``What happens when interventions are applied
solely during training?" and (ii) ``Does the class-specific effect of
interventions remain consistent if the test set remains untouched?"222The
second question is related to the analysis in Section 4.2, where we want to
clarify whether performance degradation stems from suddenly appearing unseen
domains in the test phase. According to the findings presented in Table 5, an
improvement was observed in only one out of nine instances for the Tr_S
condition. In contrast, interventions generally led to diminished performance
in the bonafide class, highlighting the necessity for more advanced methods to
accurately model genuine speech. Our findings diverge from those in [30], with
EER deteriorating under all asymmetric intervention conditions except for one.
This discrepancy is likely due to our one-side distinct strategy of applying
interventions exclusively during either the training or testing phase,
creating a mismatch. Nevertheless, similar to findings in [30], we also
observed that interventions involving loudness were generally less impactful
compared to other types of interventions across all evaluations.

##  5 Conclusion

This paper has conducted an in-depth investigation into the behavior of audio
anti-spoofing models through various experiments focused on loss analysis and
asymmetric interventions. Our analyses expand the perspective beyond attack-
centric or silence-focused interpretations. The findings suggest that current
training practices, which primarily aim to detect spoofing artifacts in known
attacks, may neglect the robust modeling of bona fide speech, potentially
introducing bias in model learning. By advocating for a more balanced focus on
understanding both bona fide and spoofed classes, our research paves the way
for future studies to enhance the efficacy of audio anti-spoofing systems.

## References

  * [1] M. Sahidullah, T. Kinnunen, and C. Hanilçi, ``A comparison of features for synthetic speech detection,'' in _Proc. Interspeech_ , 2015.
  * [2] M. Todisco, H. Delgado, and N. Evans, ``Constant q cepstral coefficients: A spoofing countermeasure for automatic speaker verification,'' _Computer Speech & Language_, vol. 45, pp. 516–535, 2017.
  * [3] H. Tak, J. Patino, M. Todisco, A. Nautsch, N. Evans, and A. Larcher, ``End-to-end anti-spoofing with Rawnet2,'' in _Proc. ICASSP_ , 2021.
  * [4] J.-w. Jung, H.-j. Shim, H.-S. Heo, and H.-J. Yu, ``Replay attack detection with complementary high-resolution information using end-to-end dnn for the asvspoof 2019 challenge,'' _Proc. Interspeech_ , 2019.
  * [5] H. Tak, J. Patino, M. Todisco, A. Nautsch, N. Evans, and A. Larcher, ``End-to-end anti-spoofing with RawNet2,'' in _Proc. ICASSP_ , 2021.
  * [6] R. K. Das, J. Yang, and H. Li, ``Data augmentation with signal companding for detection of logical access attacks,'' in _Proc. ICASSP_ , 2021.
  * [7] H. Tak, M. Kamble, J. Patino, M. Todisco, and N. Evans, ``RawBoost: A raw data boosting and augmentation method applied to automatic speaker verification anti-spoofing,'' in _Proc. ICASSP_ , 2022.
  * [8] A. Cohen, I. Rimon, E. Aflalo, and H. H. Permuter, ``A study on data augmentation in voice anti-spoofing,'' _Speech Communication_ , vol. 141, pp. 56–67, 2022.
  * [9] R. K. Das, J. Yang, and H. Li, ``Assessing the scope of generalized countermeasures for anti-spoofing,'' in _Proc. ICASSP_ , 2020.
  * [10] D. Paul, M. Sahidullah, and G. Saha, ``Generalization of spoofing countermeasures: A case study with asvspoof 2015 and btas 2016 corpora,'' in _Proc. ICASSP_ , 2017.
  * [11] H.-j. Shim, J.-w. Jung, and T. Kinnunen, ``Multi-dataset co-training with sharpness-aware optimization for audio anti-spoofing,'' in _Proc. Interspeech_ , 2023.
  * [12] I. Himawan, F. Villavicencio, S. Sridharan, and C. Fookes, ``Deep domain adaptation for anti-spoofing in speaker verification systems,'' _Computer Speech & Language_, vol. 58, pp. 377–402, 2019.
  * [13] H. Ma, J. Yi, J. Tao, Y. Bai, Z. Tian, and C. Wang, ``Continual learning for fake audio detection,'' in _Proc. Interspeech_ , 2021.
  * [14] Y. Xie, H. Cheng, Y. Wang, and L. Ye, ``Domain generalization via aggregation and separation for audio deepfake detection,'' _IEEE Transactions on Information Forensics and Security_ , 2023.
  * [15] X. Wang and J. Yamagishi, ``Investigating active-learning-based training data selection for speech spoofing countermeasure,'' in _Proc. SLT_ , 2022.
  * [16] H. Tak, M. Todisco, X. Wang, J.-w. Jung, J. Yamagishi, and N. Evans, ``Automatic speaker verification spoofing and deepfake detection using wav2vec 2.0 and data augmentation,'' in _Proc. Interspeech_ , 2022.
  * [17] X. Wang and J. Yamagishi, ``Can large-scale vocoded spoofed data improve speech spoofing countermeasure with a self-supervised front end?'' in _Proc. ICASSP_ , 2024.
  * [18] T. Chen, A. Kumar, P. Nagarsheth, G. Sivaraman, and E. Khoury, ``Generalization of audio deepfake detection,'' _Proc. Speaker Odyssey Workshop_ , 2020.
  * [19] Y. Zhang, F. Jiang, and Z. Duan, ``One-class learning towards synthetic voice spoofing detection,'' _IEEE Signal Processing Letters_ , vol. 28, pp. 937–941, 2021.
  * [20] Y. Zhang, Z. Li, J. Lu, H. Hua, W. Wang, and P. Zhang, ``The impact of silence on speech anti-spoofing,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 2023.
  * [21] Y. Ren, H. Peng, L. Li, and Y. Yang, ``Lightweight voice spoofing detection using improved one-class learning and knowledge distillation,'' _IEEE Transactions on Multimedia_ , 2023.
  * [22] B. Chettri, R. G. Hautamäki, M. Sahidullah, and T. Kinnunen, ``Data quality as predictor of voice anti-spoofing generalization,'' in _Proc. Interspeech_ , 2021.
  * [23] X. Liu, X. Wang, M. Sahidullah, J. Patino, H. Delgado, T. Kinnunen, M. Todisco, J. Yamagishi, N. Evans, A. Nautsch _et al._ , ``Asvspoof 2021: Towards spoofed and deepfake speech detection in the wild,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 2023.
  * [24] R. Geirhos, J.-H. Jacobsen, C. Michaelis, R. Zemel, W. Brendel, M. Bethge, and F. A. Wichmann, ``Shortcut learning in deep neural networks,'' _Nature Machine Intelligence_ , vol. 2, no. 11, pp. 665–673, 2020.
  * [25] T. Kinnunen, M. Sahidullah, H. Delgado, M. Todisco, N. Evans, J. Yamagishi, and K. A. Lee, ``The asvspoof 2017 challenge: Assessing the limits of replay spoofing attack detection,'' in _Proc. Interspeech_ , 2017.
  * [26] B. Chettri, E. Benetos, and B. L. Sturm, ``Dataset artefacts in anti-spoofing systems: a case study on the asvspoof 2017 benchmark,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 28, pp. 3018–3028, 2020.
  * [27] M. Todisco, X. Wang, V. Vestman, M. Sahidullah, H. Delgado, A. Nautsch, J. Yamagishi, N. Evans, T. Kinnunen, and K. Lee, ``ASVspoof 2019: Future horizons in spoofed and fake audio detection,'' in _Proc. Interspeech_ , 2019.
  * [28] N. M. Müller _et al._ , ``Speech is silver, silence is golden: What do asvspoof-trained models really learn?'' _Proc. ASVspoof 2021 Workshop (Interspeech satellite)_ , pp. 55–60, 2021.
  * [29] Y. Zhang, W. Wang, and P. Zhang, ``The effect of silence and dual-band fusion in anti-spoofing system,'' in _Proc. Interspeech_ , 2021.
  * [30] H.-j. Shim, R. G. Hautamäki, M. Sahidullah, and T. Kinnunen, ``How to construct perfect and worse-than-coin-flip spoofing countermeasures: A word of warning on shortcut learning,'' in _Proc. Interspeech_ , 2023.
  * [31] L. Wu, Z. Zhu _et al._ , ``Towards understanding generalization of deep learning: Perspective of loss landscapes,'' _ICML Workshop_ , 2017.
  * [32] J. Nam, H. Cha, S. Ahn, J. Lee, and J. Shin, ``Learning from failure: De-biasing classifier from biased classifier,'' _Advances in Neural Information Processing Systems_ , vol. 33, pp. 20 673–20 684, 2020.
  * [33] J. Lee, E. Kim, J. Lee, J. Lee, and J. Choo, ``Learning debiased representation via disentangled feature augmentation,'' _Advances in Neural Information Processing Systems_ , vol. 34, pp. 25 123–25 133, 2021.
  * [34] I. Hwang, S. Lee, Y. Kwak, S. J. Oh, D. Teney, J.-H. Kim, and B.-T. Zhang, ``Selecmix: Debiased learning by contradicting-pair sampling,'' _Advances in Neural Information Processing Systems_ , vol. 35, pp. 14 345–14 357, 2022.
  * [35] T.-Y. Lin, P. Goyal, R. Girshick, K. He, and P. Dollár, ``Focal loss for dense object detection,'' in _Proc. ICCV_ , 2017.
  * [36] T. Castells, P. Weinzaepfel, and J. Revaud, ``Superloss: A generic loss for robust curriculum learning,'' _Advances in Neural Information Processing Systems_ , vol. 33, pp. 4308–4319, 2020.
  * [37] Y. Huang, Y. Wang, Y. Tai, X. Liu, P. Shen, S. Li, J. Li, and F. Huang, ``Curricularface: adaptive curriculum learning loss for deep face recognition,'' in _Proc. CVPR_ , 2020.
  * [38] Z. Zhang and M. Sabuncu, ``Generalized cross entropy loss for training deep neural networks with noisy labels,'' _Advances in neural information processing systems_ , vol. 31, 2018.
  * [39] J. Yamagishi, X. Wang _et al._ , ``ASVspoof2021: Accelerating progress in spoofed and deep fake speech detection,'' in _Proc. ASVspoof 2021 Workshop (Interspeech satellite)_ , 2021.
  * [40] J.-w. Jung, H.-S. Heo, H. Tak, H.-j. Shim, J. S. Chung, B.-J. Lee, H.-J. Yu, and N. Evans, ``Aasist: Audio anti-spoofing using integrated spectro-temporal graph attention networks,'' in _Proc. ICASSP_ , 2022.
  * [41] A. Y. Kibangou and G. Favier, ``Wiener-Hammerstein systems modeling using diagonal volterra kernels coefficients,'' _IEEE signal processing letters_ , vol. 13, no. 6, pp. 381–384, 2006.
  * [42] P. Foret, A. Kleiner, H. Mobahi, and B. Neyshabur, ``Sharpness-aware minimization for efficiently improving generalization,'' in _Proc. ICLR_ , 2020.
  * [43] J. Kwon, J. Kim, H. Park, and I. K. Choi, ``Asam: Adaptive sharpness-aware minimization for scale-invariant learning of deep neural networks,'' in _Proc. ICML_ , 2021.
  * [44] T. Kinnunen and H. Li, ``An overview of text-independent speaker recognition: From features to supervectors,'' _Speech Communication_ , vol. 52, no. 1, pp. 12–40, 2010.
  * [45] V. Chandola, A. Banerjee, and V. Kumar, ``Anomaly detection: A survey,'' _ACM computing surveys (CSUR)_ , vol. 41, no. 3, pp. 1–58, 2009.