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

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.

READ FULL TEXT VIEW PDF

page 11

page 14

page 21

11/23/2019

Hypergeometric functions and moments of operators

We consider some hypergeometric functions and prove that they are elemen...
09/12/2020

Inflation of poorly conditioned zeros of systems of analytic functions

Given a system of analytic functions and an approximate zero, we transfo...
06/11/2018

Solutions of New Potential Integral Equations Using MLFMA Based on the Approximate Stable Diagonalization

We present efficient solutions of recently developed potential integral ...
10/20/2020

Model-specific Data Subsampling with Influence Functions

Model selection requires repeatedly evaluating models on a given dataset...
07/31/2018

A simple analysis of flying capacitor converter

The paper aims to emphasise how switched systems can be analysed with el...

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 [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: 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 [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 well-known [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 Wright-Fisher 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 proportion111We 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 as

for 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 re-define through the frequencies to remove size-dependency:

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.

Remark 3.

Similar results are available for the second moment:

, leading to equation

and we would get the variance as . 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 [25, 23].

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 Euler-Mascheroni constant () and integral

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

Figure 2: (A) Time-to-fixation function and (B) fixation-to-one function for various selection effect values and arbitrary population size .

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 n-th 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 .

Figure 3: Heatmap of calculated on a series of selection effects (A) and (B) and initial probabilities (see Quality of the Approximate Solution). Note that is not observed for these series of .

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 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 simulations with population size of . 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 , 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 .

Figure 5: Histogram of relative distances , and 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 , and 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 observed fixations for the same series of parameters. The coefficient of variation is independent of the population size (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

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 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] J. P. Bollback, T. L. York, and R. Nielsen (2008) Estimation of 2 nes from temporal allele frequency data. Genetics 179 (1), pp. 497–502. Cited by: §2.1, §4.
  • [2] J. P. Bollback, T. L. York, and R. Nielsen (2008) Estimation of 2nes from temporal allele frequency data. Genetics 179 (1), pp. 497–502. Cited by: §1.
  • [3] S. Braude and A. R. Templeton (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] P. Buri (1956) Gene frequency in small populations of mutant drosophila. Evolution, pp. 367–402. Cited by: §1.
  • [5] A. Etheridge Diffusion process models in mathematical genetics. Oxford University http://www. stats. ox. ac. uk/~ etheridg/popgen/notes. pdf. Cited by: §1.
  • [6] A. Etheridge (2011) Some mathematical models from population genetics: École d’Été de probabilités de saint-flour xxxix-2009. Vol. 2012, Springer Science & Business Media. Cited by: §2.1, §2.1.
  • [7] M. A. Fares (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] A. Ferrer-Admetlla, C. Leuenberger, J. D. Jensen, and D. Wegmann (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] M. Geller and E. W. Ng (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] D. R. Gifford, J. A. G. de Visser, and L. M. Wahl (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] Z. Gompert (2015) Bayesian inference of selection in a heterogeneous environment from genetic time-series data. Molecular Ecology 25 (1), pp. 121–134. Cited by: §1.
  • [12] R. N. Gutenkunst, R. D. Hernandez, S. H. Williamson, and C. D. Bustamante (2009) Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genetics 5 (10), pp. e1000695. Cited by: §1.
  • [13] A. Jonas, T. Taus, C. Kosiol, C. Schlötterer, and A. Futschik (2016) Estimating the effective population size from temporal allele frequency changes in experimental evolution. Genetics 204 (2), pp. 723–735. Cited by: §1.
  • [14] J. G. Kemeny and J. Snell (1976) Laurie finite markov chains springer verlag. New York. Cited by: Remark 1.
  • [15] M. Kimura and T. Ohta (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] M. Kimura (1962) On the probability of fixation of mutant genes in a population. Genetics 47 (6), pp. 713. Cited by: §1, §2.3.
  • [17] M. Kimura (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] M. Kimura (1980) 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), pp. 522–526. Cited by: §1.
  • [19] B. Koopmann, J. Müller, A. Tellier, and D. Živković (2017) Fisher–wright model with deterministic seed bank and selection. Theoretical Population Biology 114, pp. 29–39. Cited by: §1.
  • [20] M. Lacerda and C. Seoighe (2014) Population genetics inference for longitudinally-sampled mutants under strong selection. Genetics 198 (3), pp. 1237–1250. Cited by: §1.
  • [21] S. Lukić and J. Hey (2012) Demographic inference using spectral methods on SNP data, with an analysis of the human out-of-africa expansion. Genetics 192 (2), pp. 619–639. Cited by: §1.
  • [22] A. Malaspinas, O. Malaspinas, S. N. Evans, and M. Slatkin (2012) Estimating allele age and selection coefficient from time-serial data. Genetics 192 (2), pp. 599–607. Cited by: §1, §4.
  • [23] T. Maruyama (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] I. Mathieson and G. McVean (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] P. Narain (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] H. Shafiey and D. Waxman (2017) Exact results for the probability and stochastic dynamics of fixation in the wright-fisher model. Journal of theoretical biology 430, pp. 64–77. Cited by: §1, §2.1.
  • [27] R. A. F. Sir and R. Fisher (1999) The genetical theory of natural selection: a complete variorum edition. Oxford University Press. Cited by: §1.
  • [28] M. Steinrücken, A. Bhaskar, and Y. S. Song (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] P. Tataru, T. Bataillon, and A. Hobolth (2015) Inference under a wright-fisher model using an accurate beta approximation. Genetics 201 (3), pp. 1133–1141. Cited by: §1.
  • [30] P. Tataru, M. Simonsen, T. Bataillon, and A. Hobolth (2017) Statistical inference in the wright–fisher model using allele frequency data. Systematic biology 66 (1), pp. e30–e46. Cited by: §2.1.
  • [31] T. Taus, A. Futschik, and C. Schlötterer (2017) Quantifying selection with pool-seq time series data. Molecular Biology and Evolution 34 (11), pp. 3023–3034. Cited by: §1.
  • [32] J. Terhorst, C. Schlötterer, and Y. S. Song (2015) Multi-locus analysis of genomic time series data from experimental evolution. PLOS Genetics 11 (4). Cited by: §1.
  • [33] T. D. Tran, J. Hofrichter, and J. Jost (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] R. Vitalis, M. Gautier, K. J. Dawson, and M. A. Beaumont (2014) Detecting and measuring selection from gene frequency data. Genetics 196 (3), pp. 799–817. Cited by: §1.
  • [35] J. Wang, E. Santiago, and A. Caballero (2016) Prediction and estimation of effective population size. Heredity 117 (4), pp. 193–206. Cited by: §2.1.
  • [36] S. Wright (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

Appendix C Fixation time distribution on the 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 simulations with population size of .

Appendix D Coefficient of variation is independent of the population size

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