1 Introduction
The WrightFisher (WF) model is one of the simplest model of genetic drift in population genetics. It has been developed independently by Sewall Wright [36] and R.A Fisher [27] in the 1930s. This model describes the dynamic of allele transmission over time for haploid or diploid individuals based on strict assumptions: nonoverlapping generations, constant population size over time, random (panmictic) mating.
One of the first applications to the WF model was proposed by Buri et al [4] 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 [10] 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 [32, 7, 29, 8, 31, 13], 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 [12, 21, 28], selection [24, 11, 2, 22, 28, 34, 20], the effect of demographic forces including the variation in population size through time and migration connecting populations [24, 12, 21, 34].
In population genetics, one focuses on different quantities of interest such as heterozygosity or fixation probabilities [26]. 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 [16]. The expression of the fixation time for WF model with no additional effect is wellknown [5, 33]. Solutions for fixation of a rare mutant gene (excluding the case of potential loss) by random drift and with selection [15] or with an additional mutation rate [18]
have also been derived. Moreover, the variance for fixation with no loss has been found in
[25] and [23] 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 [15] and for null initial probability (only) [19] 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 WrightFisher Model With Selection
2.1 Context and Notations
We consider an initial population size of haploid individuals and study a single locus with two possible alleles: and . This model is equivalent to the study of diploid individuals who can either be homozygote (  ) or heterozygote (  ), considering that and are neither dominant nor recessive. In that case, the dominance parameter of heterozygous individuals is [6, chapter 5.3].
Notice that, in the literature [3, 35], while studying genetic drift, an effective population size is usually determined. Denoted , 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 is not the subject of this article, we consider that our population size is equal to the effective population size ().
We introduce a selection effect that affects favorably viability of carriers of allele : their chance of survival is increased by a factor . There are many possibilities to introduce selection effects [26]. We chose to keep constant for all generations and we also define coefficient normalized by population size as usually done in the literature [15, 30, 1, 6]
(estimation of
from data in [1]).Let be the number of alleles A in generation , , and its proportion^{1}^{1}1We chose notation instead of a more natural to avoid confusion with initial allele A frequency, further denoted .
: they are realizations of random variables
and respectively. The frequency of alleles A obtained at time depends on frequency at time, but is moreover skewed by the selection effect, that is, we define a fictional proportion:
where the population of allele A is increased by an imaginary quantity
. The distribution of alleles in the new generation is then generated by the binomial distribution
. The generation represents the offspring of the generation obtained by random sampling with replacement. We also define the variation of the proportion over time by . The transition probabilities are defined asfor all integers , integers in and in .
Remark 1.
is an absorbing Markov chain
[14] where is the probability of transition from state to and and are the absorbing states.2.2 Fixation Time
The fixation time is defined as the time needed by the Markov chain to reach or with initial population in , that is:
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 . Our goal is to search for analytic expressions for its expectation when we consider the limit case for infinite population size , leading to a continuous variable in . We then redefine through the frequencies to remove sizedependency:
We decompose its expected value by the law of total expectation and introduce new notations as follows:
with expectation function , fixation probability function , expected fixation time to zero and one defined as and , respectively. Function is often called the fixation time of a mutant gene excluding the case of potential loss. We also define and . These quantities can be described in terms of infinite sums:
and we can now derive equations binding neighboring values of each function to itself using the transition probability operator , which is the density probability of the transition in one time step from allele probability to allele probability .
Proposition 1.
Functions , , and satisfy the following equations for initial probability (sometimes called master equations) :
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:
We multiply this equality by and sum for all integers to get in left hand side:
Results for and are obtained by considering the same proof replacing by or in the law of total probability. For equation in (the first equation) with sum all probabilities without the multiplicative factor . ∎
Remark 2.
With discrete proportions for the initial probability and a finite size population (, ), equations in Proposition 1 are valid with integration with respect to counting measures. They can be solved exactly by matrix inversion.
2.3 Second Order Differential Equations
Master equations can be approximated by Taylor expansions around small variations in proportions .
Proposition 2.
Equations in Proposition 1 can be approximated by the following second order differential equations:
(1) 
with natural boundary conditions: , , , and , which leads to and .
Proof.
A detailed proof is given in A. All the results are based on the following expansion given here for a generic function :
∎
Solution for unknown function in first equation is easy to found (see [16]): it is given by formula and we also have .
3 Analytic Solutions for Fixation Time
3.1 The Exponential Integral Function
Our solutions are based on the exponential integral function from to which is a special function defined by an integral (understood in terms of the Cauchy principal value):
We have the following relation:
(2) 
with the EulerMascheroni constant () and integral
(3) 
The shape of the EI function is presented in Figure 1.
3.2 Time to Fixation
We start by solving the differential equation in Proposition 2 for function and determine an analytic expression.
Theorem 1.
Time to fixation to both zero and one with initial allele frequency and a selection effect of order is given by formula
with
Proof.
Equation in in Proposition 1 can be written as
which can be integrated using the function:
We obtain by a second integration using 3 (see also [9]):
We determine the constants using boundary values and the limit and we get:
This system is easy to solve for constants and to find the formula presented in the theorem. The integral expression for comes from relation (2) as the terms involving cancel out. ∎
Remark 4.
Both functions and are solution to its associated equation in (1) as this equation is invariant by transformation .
3.3 Time to Zero and Time to One
Replacing by and by in equations (1) we get the following equations:
(4) 
showing the simple relation . Therefore, we only need to solve the equation for , as leading to:
(5) 
We now derive the complete solution for function by solving its differential equation in in Proposition 2. This result is already presented in [15] in integral form.
Theorem 2.
Time to fixation excluding allele loss with initial allele frequency and a selection effect of order is given by formula:
with
The proof is straightforward: it is similar to the proof of Theorem 1 by solving the equation for in Proposition 2 and then dividing by function 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
We plot the shape of functions and for various values of the selection coefficient in Figure 2. Population size is only present as a multiplicative factor for all fixation time formulas. It does not affect the general shape of these functions.
4 An Approximate Solution
In our search for simple solutions, we propose an approximate solution for the expression of involving only elementary functions.
Theorem 3.
An approximation solution satisfying the boundary conditions is
with
and the nth harmonic number (). An upper bound is given by
Proof.
We simplify the integral term by using the series expansion of the exponential function.
Using relations and we get:
The derivation of the upper bound is exposed in B. ∎
Remark 5.
By taking the limit when tends to zero or by considering we obtain the classic solution 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:
We denote the minimum number of elements in of our approximate solution such that . Figure 3 shows for a series of selection effect values (twenty values evenly spaced within the interval [1,30] and forty values evenly spaced within the interval ) and initial probabilities (all values of within the interval with a step of size 0.05). We observe that the increases with the selection effect value . At and marginally to , .
In the literature the selection coefficient is chosen less than in the simulation study of Kumira’ seminal work [15] (). It is also common in application studies to consider small selection effects [22]. This corresponds to an approximate solution with up to 40 terms valid at precision (see Figure 3 panel A). It should be noted that considering large eventually leads to large . For instance, in the experimental MS2 bacteriophage data [1] a reasonable set of is to . The corresponding alpha then ranges from to . is computationally intractable for the upper limit of this range of . However, [1] reports estimation of this coefficient close to . For greater than the approximate and the exact solutions are computationally unstable and some new ideas has to be developed. Nevertheless, the need for large selection coefficients is unclear in the literature as most papers rely on and with no understanding of the key role played by the product .
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 fixation times with population size for three chosen selection effect values () and four chosen initial probabilities () 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 ). The exact theoretical distribution is not known but some results are known with no selection effect [17]. Note that the distribution of fixation time to one seems Logsymmetric independently of the chosen initial probability and selection effect value (see C).
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 , selection effect values (twenty values evenly spaced within the interval ) and initial probabilities (all values of within the interval with a step of size ). For each combination of these parameters we simulated up to repetitions, interrupting the simulations once we get fixations at the desired state (one, zero or both). For some combinations of and for fixation to zero, the number of repetitions was not enough to reach 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 ).
Metric
The relative distance between the empirical mean of the fixation time to zero, denoted , and the exact solution is given by the following formula
The same way we define and , the relative distance involving the fixation time to one and both, respectively.
Empirical Results
In Figure 5 we observe that for a population size or above, marginally to the selection effects and initial probabilities, the median of is lower than 0.01, i.e 1% of errors. For and the same statistical control is reached at population size and , respectively. At population size the distributions are flat which suggests a strong dispersion of relative distances. The largest observed relative distance is equal to and is obtained by with values , and , leading to absolute values and . For this combination of parameters the exact solution we propose is quite far from the simulation results. Taking now , and , and are equal to and , respectively. The corresponding relative distance is equal to which is a much better estimate. These first results convinced us to explore the relative distance according to the range of tested parameters ( and ). In Figure 6 (panel A) we observe that the relative distance is driven by three main effects. (i) The population size . Consistent with Figure 5, we observe that when increases, the relative distance decreases. (ii) The selection effect value . We observe that when increases, the relative distance increases. (iii) The side effect. We observe that and are larger for small and large , respectively. is a combination of both and . 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 , and .
6 Conclusion
Exact and approximate formulas for time to fixation in WrightFisher 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
proportion in segment. 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 codeveloped 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] (2008) Estimation of 2 nes from temporal allele frequency data. Genetics 179 (1), pp. 497–502. Cited by: §2.1, §4.
 [2] (2008) Estimation of 2nes from temporal allele frequency data. Genetics 179 (1), pp. 497–502. Cited by: §1.
 [3] (2009) Understanding the multiple meanings of ‘inbreeding’and ‘effective size’for genetic management of african rhinoceros populations. African Journal of Ecology 47 (4), pp. 546–555. Cited by: §2.1.
 [4] (1956) Gene frequency in small populations of mutant drosophila. Evolution, pp. 367–402. Cited by: §1.
 [5] Diffusion process models in mathematical genetics. Oxford University http://www. stats. ox. ac. uk/~ etheridg/popgen/notes. pdf. Cited by: §1.
 [6] (2011) Some mathematical models from population genetics: École d’Été de probabilités de saintflour xxxix2009. Vol. 2012, Springer Science & Business Media. Cited by: §2.1, §2.1.
 [7] (2015) 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, pp. 101–113. Cited by: §1.

[8]
(2016)
An approximate markov model for the wright–fisher diffusion and its application to time series data
. Genetics 203 (2), pp. 831–846. Cited by: §1.  [9] (1969) A table of integrals of the exponential integral. Journal of Research of the National Bureau of Standards 71, pp. 1–20. Cited by: §3.2.
 [10] (2013) Model and test in a fungus of the probability that beneficial mutations survive drift. Biology letters 9 (1), pp. 20120310. Cited by: §1.
 [11] (2015) Bayesian inference of selection in a heterogeneous environment from genetic timeseries data. Molecular Ecology 25 (1), pp. 121–134. Cited by: §1.
 [12] (2009) Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genetics 5 (10), pp. e1000695. Cited by: §1.
 [13] (2016) Estimating the effective population size from temporal allele frequency changes in experimental evolution. Genetics 204 (2), pp. 723–735. Cited by: §1.
 [14] (1976) Laurie finite markov chains springer verlag. New York. Cited by: Remark 1.
 [15] (1969) The average number of generations until fixation of a mutant gene in a finite population. Genetics 61 (3), pp. 763. Cited by: §1, §1, §2.1, §3.3, §4.
 [16] (1962) On the probability of fixation of mutant genes in a population. Genetics 47 (6), pp. 713. Cited by: §1, §2.3.
 [17] (1970) 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), pp. 131–133. Cited by: §5.1.
 [18] (1980) Average time until fixation of a mutant allele in a finite population under continued mutation pressure: studies by analytical, numerical, and pseudosampling methods. Proceedings of the National Academy of Sciences 77 (1), pp. 522–526. Cited by: §1.
 [19] (2017) Fisher–wright model with deterministic seed bank and selection. Theoretical Population Biology 114, pp. 29–39. Cited by: §1.
 [20] (2014) Population genetics inference for longitudinallysampled mutants under strong selection. Genetics 198 (3), pp. 1237–1250. Cited by: §1.
 [21] (2012) Demographic inference using spectral methods on SNP data, with an analysis of the human outofafrica expansion. Genetics 192 (2), pp. 619–639. Cited by: §1.
 [22] (2012) Estimating allele age and selection coefficient from timeserial data. Genetics 192 (2), pp. 599–607. Cited by: §1, §4.
 [23] (1972) 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), pp. 109–113. Cited by: §1, Remark 3.
 [24] (2013) Estimating selection coefficients in spatially structured populations from time series data of allele frequencies. Genetics 193 (3), pp. 973–984. Cited by: §1.
 [25] (1970) 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), pp. 251–255. Cited by: §1, Remark 3.
 [26] (2017) Exact results for the probability and stochastic dynamics of fixation in the wrightfisher model. Journal of theoretical biology 430, pp. 64–77. Cited by: §1, §2.1.
 [27] (1999) The genetical theory of natural selection: a complete variorum edition. Oxford University Press. Cited by: §1.
 [28] (2014) A novel spectral method for inferring general diploid selection from time series genetic data. The Annals of Applied Statistics 8 (4). Cited by: §1.
 [29] (2015) Inference under a wrightfisher model using an accurate beta approximation. Genetics 201 (3), pp. 1133–1141. Cited by: §1.
 [30] (2017) Statistical inference in the wright–fisher model using allele frequency data. Systematic biology 66 (1), pp. e30–e46. Cited by: §2.1.
 [31] (2017) Quantifying selection with poolseq time series data. Molecular Biology and Evolution 34 (11), pp. 3023–3034. Cited by: §1.
 [32] (2015) Multilocus analysis of genomic time series data from experimental evolution. PLOS Genetics 11 (4). Cited by: §1.
 [33] (2013) An introduction to the mathematical structure of the wright–fisher model of population genetics. Theory in Biosciences 132 (2), pp. 73–82. Cited by: §1.
 [34] (2014) Detecting and measuring selection from gene frequency data. Genetics 196 (3), pp. 799–817. Cited by: §1.
 [35] (2016) Prediction and estimation of effective population size. Heredity 117 (4), pp. 193–206. Cited by: §2.1.
 [36] (1931) Evolution in mendelian populations. Genetics 16 (2), pp. 97. Cited by: §1.
Appendix A Differential Equations for Fixation Time
We describe the proof for the equation in . All other results are proven the same way (we only have to change the constant term). Using a Taylor expansion, we get
(6)  
With , we have:
and
We derive approximate formulations for these two first moments by Taylor expansion in power of using relation . We get:
Similarly,
Removing the small order terms leads to the desired differential equation.
Appendix B Approximate Solution Upper Bound
We recall that is positive.
applying triangle inequalities and using relations , . With the loose upper bound and for positive , we get:
using the upper bound for both and . Eventually, we have:
and we get