An introduction to the mathematical structure of the Wright–Fisher model of population genetics

In this paper, we develop the mathematical structure of the Wright–Fisher model for evolution of the relative frequencies of two alleles at a diploid locus under random genetic drift in a population of fixed size in its simplest form, that is, without mutation or selection. We establish a new concept of a global solution for the diffusion approximation (Fokker–Planck equation), prove its existence and uniqueness and then show how one can easily derive all the essential properties of this random genetic drift process from our solution. Thus, our solution turns out to be superior to the local solution constructed by Kimura.


Introduction
In population genetics, one considers the effects of recombination, selection, mutation, and perhaps others like migration on the distribution of alleles in a population, see e.g. (Ewens 2004;Bürger 2000;Rice 2004) as mathematical textbook references. The most basic and at the same time important model is the Wright-Fisher model for random genetic drift [developed implicitly by Fisher (1922) and explicitly by Wright (1931)]. In its simplest version-the one to be treated in the present paper-it is concerned with the evolution of the relative frequencies of two alleles at a single diploid locus in a finite population of fixed size with non-overlapping generations under the sole force of random genetic drift, without any other influences like mutations or selection. The model can be generalised-and so can our approach-to multiple alleles, several loci, with mutations, selections, spatial population structures, etc, see the above references. To find an exact solution (for the approximating diffusion process for the probability densities of the allele frequencies described by a Fokker-Planck equation) from which the properties of the resulting stochastic process can be deduced, however, is difficult. For the basic two-allele case, this was first achieved in the important work of , and he then went on to treat the case of several alleles (Kimura , 1956. His solution, however, is local in the sense that it does not naturally incorporate the transitions resulting from the irreversible loss of one or several of the alleles initially present in the population. Consequently, the resulting probability distribution does not integrate to 1, and it is difficult to read off the quantitative properties of the process from his solution. In the present paper, we introduce and describe a new global approach. This approach is mathematically more transparent than Kimura's scheme. We prove the existence of a unique such global solution (see Theorem 3.7), and we can deduce all desired quantities of the underlying stochastic process from our solution. The purpose of the present paper thus is to display the method in the simplest case, that of two alleles at a single locus, so that the structure becomes clear. The case of multiple alleles is presented in our companion paper (Tran et al. 2000) on the basis of the first author's thesis, and further generalisations will be systematically developed elsewhere within the mathematical framework of information geometry (Amari and Nagaoka 2000) and more specifically (Ay and Jost 2000;Jost 2000) on the basis of the second author's thesis.

The Wright-Fisher model
We consider a diploid population of size N. At a given locus, there could be either one of the two alleles A 1 ,A 2 . Thus, an individual can be a homozygote of type A 1 A 1 or A 2 A 2 or a heterozygote of type A 1 A 2 or A 2 A 1 -but we consider the latter two as the same-at the locus in question. The population reproduces in discrete time steps, and each individual in generation n ? 1 inherits one allele from each of its parents. When a parent is a heterozygote, each allele is chosen with probability 1/2. Here, for each individual in generation t ? 1, randomly two parents in generation n are chosen. Thus, the alleles in generation n ? 1 are chosen by random sampling with replacement from the ones in generation n. The quantity of interest is the number Y n of alleles A 1 in the population at time n. This number then varies between 0 and 2N. The transition probability then is whenever Y n takes the value 0 or 2N, that is, if either the allele A 1 or A 2 will disappear, it will stay there for all future times. Eventually, this will happen almost surely. This is the basic model. One can then derive expressions for the expected time for the allele A 1 to become either fixed, that is, Y n = 2N, or become extinct, Y n = 0, given its initial number Y 0 .
An important idea, first applied in Wright (1945), then is to rescale time and population size via and then consider the limit N ! 1: The rescaling of (2) yields a discrete Markov chain X t valued in f0; 1 2N ; . . .; 1g with t = 1 now corresponding to 2N generations. One readily verifies that the expectation values for the variation across generations satisfy A basic idea of our approach is to consider the kth moment m k (t) of the distribution about zero at the (2Nt)th generation, i.e.
We have Expanding the right hand side and noting (3) we obtain the following recursion formula when we assume that the population number N is so large that we can neglect all terms of order at least 1 N 2 : Under this assumption, the moments change very slowly per generation and we can replace the above system (6) by the system of differential equations _ m k ðtÞ ¼ À kðk À 1Þ 2 m k ðtÞ þ kðk À 1Þ 2 m kÀ1 ðtÞ; where the dot denotes a derivative w.r.t. the variable t.
These formulae now guide us in finding a continuous process that well approximates the above discrete process. We seek a continuous Markov process {X t } t C 0 valued in [0,1] with the same conditions as (3) and (7). The conditions (3) imply (see for example Ewens 2004, p. 137, for a derivation) that the probability density function u(x, t) of this continuous process is a solution of the Fokker-Planck (Kolmogorov forward) equation where we now use the notation u t :¼ o ot uðx; tÞ for the partial derivative w.r.t. the time variable t. The coefficient x(1x) in (8) comes from (3) and d p denotes the Dirac delta function at p. For the definition of this delta function, we use the product ðf ; gÞ :¼ for square integrable functions f ; g : ½0; 1 ! R on the unit interval (this will be described in more detail in ''Existence and uniqueness of solutions''), and we then put ðd p ; /Þ :¼ /ðpÞ whenever / : ½0; 1 ! R is a continuous function. 1 Let us also explain the interpretation of (8) for those not sufficiently versed in this mathematical formalism. The initial condition u(x,0) = d p (x) then simply says that at time 0, the relative frequency of allele A 1 is precisely p, without any uncertainty (this assumption is not essential, however, and the scheme works also for more general initial condition involving uncertainty about the initial distribution of the alleles). Subsequently, this allele frequence evolves stochastically, according to the equation and therefore, for t [ 0, we no longer know the precise value of this relative frequency, but only its probability density given by u(x, t). That is, for every x, the probability density that the allele frequency at time t has the value x is given by u(x, t).
In the continuum limit, the kth moment becomes R 1 0 uðx; tÞx k dx; and the condition (7) then implies Since the polynomials are dense in the space of (square integrable) functions, this yields for all square integrable functions / : ½0; 1 ! R that are twice differentiable in the open interval (0, 1), with the differential operator and its formal adjoint This solution concept will allow us to prove the existence of a unique solution from which we can then derive all features of interest of the Wright-Fisher process. We should point out that (12) is not just the integration by parts of (10), but also includes the boundary behaviour (of course, this may not be overt, but the mathematical trick here is to represent this boundary behaviour in an implicit form best suited for formal manipulation). It, thus, reflects transitions from the presence of both alleles to the irreversible loss of one of them. This is the crucial difference to  solution concept and the key for the properties of our solution.

Existence and uniqueness of solutions
We shall now apply a familiar mathematical scheme for the construction of a solution of a differential equation, an expansion in terms of eigenfunctions of the differential operator involved. For our problem, as formalised in Definition 2.1, these eigenfunctions can be constructed from a classical family of polynomials, the Gegenbauer polynomials, which we shall now introduce.

Preliminaries
For the sequel, we shall need some more notation. We need the function spaces To construct solutions in terms of expansions, we shall need a special case of the Gegenbauer polynomials [named after Leopold Gegenbauer (1849-1903]. 2 The polynomials Y m (z) we need are defined in terms of their generating function • The Gegenbauer polynomials satisfy the recurrence relation • The Gegenbauer polynomials solve the differential equation Lemma 3.2 [Abramowitz (1965), p. 774] The polynomials Y m are orthogonal polynomials on the interval [-1,1] with respect to the weight function (1z 2 ): Auxiliaries Lemma 3.3 For all m C 0 we have is a Gegenbauer polynomial and therefore solves (15), This is equivalent to This completes the proof. h In the sequel, we shall use the abbreviation wðxÞ :¼ xð1 À xÞ: Lemma 3.4 If X is an eigenvector of L corresponding to the eigenvalue k then wX is an eigenvector of L * corresponding to the eigenvalue k.
Proof Assume that X is an eigenvector of L for the eigenvalue k, i.e.
:¼ K; and the eigenvector of L corresponding to k m is the Gegenbauer polynomial X m (x) (up to a constant).
Proof From Lemma 3.3 we have L(X m ) = -k m X m in H 0 . So, K SpecðLÞ: Conversely, we shall prove that k 6 2 K is not an eigenvalue of L. In fact, assume that there is some For any n C 0, we can multiply this relation by wX n and then integrate on [0,1]. From the orthogonality (16) with respect to the weight function w, we obtain d n k n ðX n ; wX n Þ ¼ d n kðX n ; wX n Þ: Because (X n ,wX n )= 0 and k= k n , then d n = 0, V n C 0. Therefore, X = 0, i.e. k is not an eigenvalue of L. Thus

SpecðLÞ ¼ K:
Similarly, if X is an eigenvector of L for the eigenvalue k m , we will prove that X = cX m . In fact, representing X ¼ P 1 n¼0 d n X n ; it follows that X 1 n¼0 d n ðÀk n X n Þ ¼ X 1 n¼0 d n LðX n Þ ¼ L X 1 n¼0 d n X n ! ¼ Àk m X 1 n¼0 d n X n : For any k C 0, we multiply this relation by wX k and then integrate on [0,1] to obtain d k k k ðX k ; wX k Þ ¼ d k k m ðX k ; wX k Þ: Because (X k ,wX k ) = 0 and k m = k k for all k = m, then d k = 0, Vk = m. Hence X = d m X m . This completes the proof. h

Construction of the solution
In this subsection, we construct the solution and prove its uniqueness. We shall firstly find the general solution of the Fokker-Planck equation (10) by the separation of variables method. Then we shall construct a solution depending on parameters. We shall use (11, 12) to determine the parameters. Finally, we shall verify the solution.
Step 1 Assume that u 0 (x,t) = X(x)T(t) is a solution of the Fokker Planck equation (10). Then we have which implies that k is a constant which is independent of t, x. From Lemma 3.5, we obtain the general solution of the equation (10) as c m X m ðxÞe Àk m t : Remark 3.6 u 0 is the same as Kimura's solution (see for example Kimura 1955a,b).
Step 2 The general solution u 2 H of (10) then is where d 0 and d 1 are the Dirac delta functionals at 0 and 1.

Applications
Our global solution readily yields the quantities of interest of the evolution of the process (X t ) t C 0 such as the expectation and the second moment of the absorption time, mth moments, fixation probabilities, the probability of coexistence, or the probability of heterogeneity.

Absorption time
Let V 0 : = {0,1} be the domain representing a population of 1 allele. Here, 0 corresponds to the loss of A 1 , that is, the fixation of A 2 , and 1 corresponds to the opposite situation. Either of these irreverible events is called an absorption.
We denote by T 1 2 ðpÞ ¼ infft [ 0 : X t 2 V 0 jX 0 ¼ pg the first time when the population has only 1 allele left, that is, when absorption occurs. T 2 1 (p) is a continuous random variable valued in ½0; 1Þ with probability density function denoted by /(t, p). V 0 is invariant (absorption set) under the process X t , i.e. if X s 2 V 0 then X t 2 V 0 for all t C s. We have Therefore the expectation of the absorption time for having only one allele is Âð2m þ 2Þ 2 X 2m ðpÞ: and its second moment is t 2 e Àk m t dt dx Remark 4.1 EðT 1 2 ðpÞÞ ¼ À2fp lnðpÞ þ ð1 À pÞ lnð1 À pÞg is the unique solution of the one-dimensional boundary value problem Lv ¼ À1; in (0,1) vð0Þ ¼ vð1Þ ¼ 0: We easily check that this agrees with our formula above by using Mathematica (Fig. 4): nth moments By induction, it is easy to prove that Z ½0;1 x n X mÀ1 ðxÞ dx ¼ x n X m ðxÞ dx This nth moment coincides with  one.

Conclusion
We have constructed a unique global solution of the Fokker-Planck equation associated with the Wright-Fisher model. This solution leads to explicit formulae for the absorption time, fixation probabilities, the probability of coexistence, nth moments, heterogeneity, and other quantities.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.