International Journal of Stochastic AnalysisVolume 2011 (2011), Article ID 605068, 37 pagesdoi:10.1155/2011/605068

Research Article

Nonconservative Diffusions on with Killing and Branching: Applications to Wright-Fisher Models with or without Selection

Thierry E. Huillet

Laboratoire de Physique Théorique et Modélisation, CNRS-UMR 8089 et Université de Cergy-Pontoise, 2 Avenue Adolphe Chauvin, 95302 Cergy-Pontoise, France

Received 22 November 2010; Accepted 9 May 2011

Academic Editor: Manuel O. Cáceres

Copyright © 2011 Thierry E. Huillet. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We consider nonconservative diffusion processes on the unit interval, so with absorbing barriers. Using Doob-transformation techniques involving superharmonic functions, we modify the original process to form a new diffusion process presenting an additional killing rate part . We limit ourselves to situations for which is itself nonconservative with upper bounded killing rate. For this transformed process, we study various conditionings on events pertaining to both the killing and the absorption times. We introduce the idea of a reciprocal Doob transform: we start from the process , apply the reciprocal Doob transform ending up in a new process which is but now with an additional branching rate , which is also upper bounded. For this supercritical binary branching diffusion, there is a tradeoff between branching events giving birth to new particles and absorption at the boundaries, killing the particles. Under our assumptions, the branching diffusion process gets eventually globally extinct in finite time. We apply these ideas to diffusion processes arising in population genetics. In this setup, the process is a Wright-Fisher diffusion with selection. Using an exponential Doob transform, we end up with a killed neutral Wright-Fisher diffusion . We give a detailed study of the binary branching diffusion process obtained by using the corresponding reciprocal Doob transform.

1. Introduction

We consider diffusion processes on the unit interval with a series of elementary stochastic models arising chiefly in population dynamics in mind. These connections found their way over the last sixty years, chiefly in mathematical population genetics. In this context, we refer to [1] and to its extensive and nonexhaustive list of references for historical issues in the development of modern mathematical population genetics (after Wright, Fisher, Crow, Kimura, Nagylaki, Maruyama, Ohta, Watterson, Ewens, Kingman, Griffiths, and Tavaré, to cite only a few). See also the general monographs [2–6].

Special emphasis is put on Doob-transformation techniques of the diffusion processes under concern. Most of the paper's content focuses on the specific Wright-Fisher (WF) diffusion model and some of its variations, describing the evolution of one two-locus colony undergoing random mating, possibly under the additional actions of mutation, selection, and so on. We now describe the content of this work in more detail.

Section 2 is devoted to generalities on one-dimensional diffusions on the unit interval [0, 1]. It is designed to fix the background and notations. Special emphasis is put on the Kolmogorov backward and forward equations, while stressing the crucial role played by the boundaries in such one-dimensional diffusion problems. Some questions such as the meaning of speed and scale functions, existence of an invariant measure, and validity of detailed balance are addressed in the light of the Feller classification of boundaries. When the boundaries are absorbing, the important problem of evaluating additive functionals along sample paths is then briefly discussed, emphasizing the prominent role played by the Green function of the model; several simple illustrative examples are supplied. So far, we have dealt with a given process, say , and recalled the various ingredients for computing the expectations of various quantities of interest, summing up over the history of paths. In this setup, there is no distinction among paths with different destinations, nor did we allow for annihilation or creation of paths inside the domain before the process reached one of the boundaries. The Doob transform of paths allows to do so. We, therefore, describe the transformation of sample paths techniques deriving from superharmonic additive functionals. Some Doob transformations of interest are then investigated, together with the problem of evaluating additive functionals of the transformed diffusion process itself. Roughly speaking, the transformation of paths procedure allows to select sample paths of the original process with, say, a fixed destination and/or, more generally, to kill certain sample paths that do not fit the integral criterion encoded by the additive functional. As a result, this selection of paths procedure leads to a new process described by an appropriate modification of the infinitesimal generator of the original process including a multiplicative killing part rate of the sample paths inside the interval. It turns out, therefore, that the same diffusion methods used in the previous discussions apply to the transformed processes, obtained after a change of measure.

Let us be more specific. In this work, we limit ourselves to nonconservative diffusion processes on the unit interval and so with absorbing barriers. Using Doob-transformation techniques involving superharmonic functions , we modify the original process to form a new diffusion process presenting an additional killing rate part . We further limit ourselves to situations for which is itself nonconservative with bounded above killing rate. For a large class of diffusion processes, the exponential function or some linear combinations of exponential functions are admissible superharmonic functions , leading to the required property on . The full transformed process has two stopping times: the time to absorption to the boundaries and the killing time inside the domain. We study various conditionings of the transformed process: conditioning on events leading to both random stopping times occurring after the current time or only in the remote future and conditioning on events leading to either killing or absorption time occurring first. We give the relevant quasistationary limit laws, in the spirit of Yaglom [7]. This is made possible thanks to the existence of an harmonic function for the full infinitesimal generator of the transformed process.

We next introduce the idea of a reciprocal Doob transform: we start now from the process , apply the reciprocal Doob transform ending up in a new process which is but now with an additional branching rate , which is bounded. Under this reciprocal technique, the particles are not killed, rather they are allowed either to survive or split. The transformed process is a binary branching diffusion. For this supercritical binary branching diffusion process, there is a tradeoff between branching events giving birth to new particles and absorption at the boundaries, killing the particles. Under our assumptions, the branching diffusion process gets eventually globally extinct in finite time.

We next apply these general ideas to diffusion processes arising in population genetics.

In Section 3 we start recalling that Wright-Fisher diffusion models with various drifts are continuous space-time models which can be obtained as scaling limits of a biased discrete Galton-Watson model with a conservative number of offspring over the generations. Sections 4 and 5 are devoted to a detailed study of both the neutral WF diffusion process and the WF diffusion with selection, respectively.

In Section 6, we apply the Doob-transformation techniques to these processes: The starting point process is a Wright-Fisher diffusion with selection differential . We use the exponential Doob kernel . The transformed process accounts for a neutral Wright-Fisher evolution for the allele 1 frequency, subject to the additional possibility of the extinction of the population itself due to killing at rate proportional to its heterozygosity. This model is of importance in population genetics as it first appeared in [8, Page 272] as a scaling limit of a discrete population genetics model of recombination. We particularize the relevant Yaglom limit laws obtained after conditionings on events pertaining to both the killing or the absorption times occurring first. The computations of the quasistationary distributions are explicit here. Our approach relies on the spectral expansion of the transition probability kernels of both and which are known (from the works of Kimura) to involve oblate spheroidal wave functions and Gegenbauer polynomials, respectively.

In Section 7, we follow the general reciprocal path indicated in Section 2 and apply it to the particular models under concern, thereby illustrating and developing the idea of a reciprocal Doob transform. We give a detailed study of the binary branching diffusion process obtained by using the corresponding reciprocal Doob transform when the starting point process is now a neutral Wright-Fisher diffusion process. We end up in a globally subcritical branching particle system, each diffusing according to the WF model with selection. This problem is amenable to the results obtained in [9, 10].

2. Diffusion Processes on The Unit Interval: A Reminder

We start with generalities on one-dimensional diffusions exemplifying our study to the Wright-Fisher model and its relatives. For more technical details, we refer to [8, 11–13].

2.1. Generalities on One-Dimensional Diffusions on the Interval

Let be a standard one-dimensional Brownian (Wiener) motion. Consider a one-dimensional Itô diffusion driven by on the interval say ; see [14]. We will let . Assume that it has locally Lipschitz continuous drift and local standard deviation (volatility) , namely, consider the stochastic differential equation (SDE) The condition on and guarantees in particular that there is no point in for which or would blow up and diverge as .

The Kolmogorov backward infinitesimal generator of (2.1) is . As a result, for all suitable in the domain of the operator , satisfies the Kolmogorov backward equation (KBE) In the definition of the mathematical expectation , we have , where indicates a random time at which the process should possibly be stopped (absorbed), given the process was started in . The description of this (adapted) absorption time is governed by the type of boundaries which are to .

2.2. Natural Coordinate, Scale, and Speed Measure

For such Markovian diffusions, it is interesting to consider the harmonic coordinate belonging to the kernel of that is, satisfying . For and its derivative , with , one finds One should choose a version of satisfying , . The function kills the drift of in the sense that considering the change of variable , The driftless diffusion is often termed the diffusion in natural coordinates with state-space . Its volatility is . The function is often called the scale function.

Whenever and , one can choose the integration constants defining so that with and . In this case, the state-space of is again , the same as for .

Finally, considering the random time change with inverse: defined by and the novel diffusion is easily checked to be identical in law to a standard Brownian motion. Let now weak- stand for the Dirac delta mass at . The random time can be expressed as where is the (positive) speed density at and the local time at of the Brownian motion before time . Both the scale function and the speed measure are, therefore, essential ingredients to reduce the original stochastic process to the standard Brownian motion . Indeed, it follows from the above arguments that if , then is a Brownian motion. The Kolmogorov backward infinitesimal generator may be written in Feller form

Examples 2.2 (from population genetics). (i) Assume that and . This is the neutral Wright-Fisher (WF) model discussed at length later. This diffusion is already in natural scale and , . The speed measure is not integrable.(ii) With , assume and . This is the Wright-Fisher model with mutation. The parameters can be interpreted as mutation rates. The drift vanishes when which is an attracting point for the dynamics. Here,,, with and if . The speed measure density is and so is always integrable.(iii) With , assume a model with quadratic logistic drift and local variance . This is the WF model with selection. For this diffusion (see [15]), and are not integrable. Here, is a selection or fitness parameter. We shall return at length to this model and its neutral version later.

2.3. The Transition Probability Density

Assume that and are now differentiable in . Let then stand for the transition probability density function of at given . Then, is the smallest solution to the Kolmogorov forward (Fokker-Planck) equation (KFE) where is the adjoint of ( acts on the terminal value , whereas acts on the initial value ). The way one can view this PDE depends on the type of boundaries that are.

We will next suppose that the boundaries or 1 are both exit (or absorbing) boundaries. From the Feller classification of boundaries, this will be the case if where a function if .

In this case, a sample path of can reach from the inside of in finite time but cannot reenter. The sample paths are absorbed at . There is an absorption at at time and . Whenever both boundaries are absorbing, the diffusion should be stopped at . Would none of the boundaries be absorbing, then , which we rule out.

Examples of diffusion with exit boundaries are WF model and WF model with selection. In the WF model including mutations, the boundaries are entrance boundaries and so are not absorbing.

When the boundaries are absorbing, is a subprobability. Letting , we clearly have . Such models are nonconservative.

For one-dimensional diffusions, the transition density is reversible with respect to the speed density ([8, Chapter 15, Section  13]) and so detailed balance holds The speed density satisfies . It may be written as a Gibbs measure with density: , where the potential function reads and with the measure standing for the reference measure.

Further, if is the transition probability density from to , , then , with terminal condition and so also satisfies the KBE when looking at it backward in time. The Feller evolution semigroup being time homogeneous, one may as well observe that with , operating the time substitution , itself solves the KBE In particular, integrating over , , with .

being a sub-probability, we may define the normalized conditional probability density , now with total mass 1. We get The term is the time-dependent birth rate at which mass should be created to compensate the loss of mass of the original process due to absorption of at the boundaries. In this creation of mass process, a diffusing particle started in dies at rate at point , where it is duplicated in two new independent particles both started at (resulting in a global birth) evolving in the same diffusive way (consider a diffusion process with forward infinitesimal generator governing the evolution of . Suppose that a sample path of this process has some probability that it will be killed or create a new copy of itself and that the killing and birth rates and depend on the current location of the path. Then, the process with the birth and death opportunities of a path has the infinitesimal generator , where . The rate can also depend on and ). The birth rate function depends here on and , not on .

When the boundaries of are absorbing, the spectra of both and are discrete (see [8, Page 330]): There exist positive eigenvalues ordered in ascending sizes and eigenvectors of both and satisfying and such that with and , the spectral representation holds.

Let be the smallest nonnull eigenvalue of the infinitesimal generator (and of ). Clearly, and by L' Hospital rule, therefore, . Putting in the latter evolution equation, independently of the initial condition where is the eigenvector of associated to , satisfying . The limiting probability norm (after a proper normalization) is called the Yaglom limit law of conditioned on being currently alive at all time (see [7]).

2.4. Additive Functionals Along Sample Paths

Let be the diffusion model defined by (2.1) on the interval , where both endpoints are assumed absorbing (exit). This process is, thus, transient and nonconservative. We wish to evaluate the nonnegative additive quantities where the functions and are both assumed nonnegative on and . The functional solves the Dirichlet problem and is a superharmonic function for , satisfying .

Some Examples.(1) Assume that and here, is the mean time of absorption (average time spent in before absorption), solution to (2) Whenever both are exit boundaries, it is of interest to evaluate the probability that first hits (say) at 1, given . This can be obtained by choosing and .Let then is a -harmonic function solution to , with boundary conditions and . Solving this problem, we get On the contrary, choosing to be a -harmonic function with boundary conditions and , .(3) Let and put and . As , converges weakly to and, is the Green function, solution to is, therefore, the mathematical expectation of the local time at , starting from (the sojourn time density at ). The solution is known to be (see [8, page 198] or [5, page 280]) The Green function is of particular interest to solve the general problem of evaluating additive functionals . Indeed, as is well known, see [8], for example, the integral operator with respect to the Green kernel inverts the second-order operator leading to Under this form, appears as a potential function and all potential function is superharmonic. Note that for all harmonic function satisfying , is again superharmonic because .(4) Also of interest are the additive functionals of the type where the functions and are again both assumed to be nonnegative. The functional solves the Dynkin problem, [8] involving the action of the resolvent operator on . Whenever , then is the potential function, solution to is, therefore, the mathematical expectation of the exponentially damped local time at , starting from (the temporal Laplace transform of the transition probability density from to at ), with . Then, it holds that The -potential function is also useful in the computation of the distribution of the first-passage time to starting from . From the convolution formula, and taking the Laplace transform of both sides with respect to time, we obtain the Laplace-Stieltjes transform (LST) of the law of as We have as a result of both terms in the ratio being finite and belonging to the same transience class of the process (under our assumptions that the boundaries are absorbing). Note that from the reversibility property

2.5. Transformation of Sample Paths (Doob-Transform) and Killing

In the preceding subsections, we have dealt with a given process and recalled the various ingredients for the expectations of various quantities of interest, summing over the history of paths. In this setup, there is no distinction among paths with different destinations nor did we allow for annihilation or creation of paths inside the domain before the process reached one of the boundaries. The Doob transform of paths allows to do so.

Consider a one-dimensional diffusion as in (2.1) with absorbing barriers. Let be its transition probability, and let be its absorption time at the boundaries.

Let be a nonnegative additive functional solving Recall the functions and are both chosen nonnegative so that so is .

Define a new transformed stochastic process by its transition probability In this construction of through a change of measure, sample paths of for which is large are favored. This is a selection of paths procedure due to Doob (see [11]).

Now, the KFE for clearly is , with and . The Kolmogorov backward operator of the transformed process is, therefore, by duality Developing, with and , we get and the new KB operator can be obtained from the latter by adding a drift term to the one in of the original process to form a new process with the KB operator and by killing its sample paths at death rate (provided ). Note that In others words, with , the novel time-homogeneous SDE to consider is possibly killed at rate as soon as . Whenever is killed, it enters conventionally into some coffin state added to the state-space. Let be the new absorption time at the boundaries of started at (with would the boundaries be inaccessible to the new process which we ruled out). Let be the killing time of started at (the hitting time of ), with if . Then, is the novel stopping time of . The SDE for , together with its global stopping time characterize the new process with generator to consider.

In the sequel, we shall limit ourselves to the cases for which the following additional conditions hold on the transformed process.

Nonconservativeness of .We will next suppose that the boundaries or 1 are both exit (or absorbing) boundaries for the new process in (2.38). From the Feller criterion for exit boundaries, this will be the case if where is the new speed measure density for and its scale function. Recalling and , we have So, we assume here that obeys itself a nonconservative diffusion.

Boundedness of the Killing Rate .In some examples, the killing rate is bounded above. For example, suppose that the drift of the diffusion process is bounded above by . (If the drift of is bounded below by , we are led to the same conclusions while considering the process instead of .) Then, choosing . Thus, is bounded above by . Because , all this makes sense if, , or (the opposite of the gradient of the potential function in (2.12) is bounded below). Let be a nonincreasing sequence of valued real numbers. Let be a sequence of nonnegative real numbers such that for all Whenever is bounded above and, , , we have Thus, is bounded above by . Therefore, for a large class of diffusion processes, the exponential function or some linear combinations of exponential functions are superharmonic functions , leading to a bounded above killing rate .

2.6. Normalizing and Conditioning

Because the transformed process is nonconservative, it is of interest to inspect various conditionings in the sense of Yaglom, [7].

Consider again the process with infinitesimal generator losing mass due to killing and/or absorption at the boundaries. Integrating over , with , we have with . This gives the tail distribution of the full stopping time .

Defining the conditional probability density , now with total mass 1, with , we get The term is the rate at which mass should be created to compensate the loss of mass of the process due to its possible absorption at the boundaries and/or killing. Again, we have , where is the smallest positive eigenvalue of , and therefore, putting in the latter evolution equation, we get that independently of the initial condition where is the solution to With the eigenvector of associated to , is of the product form where . This results directly from the fact that and that is the stated eigenvector of . A different way to see this is as follows. We have and the conditional density of given is, therefore, The rest follows from observing that, to the leading order in in (2.15), for large time where (.) is the eigenvector of (.) associated to and . From this, it is clear that and The limiting probability norm can, therefore, be interpreted as the Yaglom limit law of conditioned on the event .

Under our assumptions, in the transformation of paths process, the transformed process can both be absorbed at the boundaries and be killed. So, both and are finite with positive probability. We wish to understand the processes conditioned on the events or , (see [16]).

The probability mass cumulated at the boundaries by time clearly is [17]

As , this probability tends to . Note that

Now (assuming ),

Thus, is defined by [or ], with boundary conditions . It serves as a positive harmonic function for . This is a Sturm-Liouville problem to be solved for each case study.

The density of the process conditioned on the event is The density of the process conditioned on the event is Note that and, with , as required, because this is the probability that is neither in nor in state at time .

Note also that ] are the transition probability densities of conditioned on the event [of conditioned on the event ]. They are the Yaglom limits of both conditioned processes.

The backward infinitesimal generators of both processes with transition probability densities and are, respectively, given by We get, respectively, Thus, in , there is no multiplicative part (no killing) and a shift in the drift, showing that the associated conditioned process obeys the SDE with drift This process is ultimately absorbed at .

In , there is a killing multiplicative part which is enhanced and a shift in the drift, showing that the associated conditioned process exhibits a faster killing rate, but the drift shift guarantees that is not absorbed at the boundaries. We have

Additive Functionals of the Transformed Process.for the new process , it is also of interest to evaluate additive functionals along their own sample paths. Let then be such an additive functional where the functions and are themselves both nonnegative. It solves Then, recalling the expression of the Green function of in (2.22), we find explicitly

Specific transformations of interest. The case deserves a special treatment. Indeed, in this case, and so , the absorption time for the process governed by the new SDE (2.38). Here, . Assuming solves if with boundary conditions and ( and .), the new process is just conditioned on exiting at (at .). In the first case, the boundary 1 is exit, whereas 0 is entrance; reads with giving the new drift. In the second case, =/, and the boundary 0 is exit, whereas 1 is entrance. Thus, is just the exit time at (at .). Let . Then, solves , whose explicit solution is in terms of , the Green function of .Example 2.1. Consider the WF model on with selection for which, with and . Assume that solves if with and ; one gets, . The diffusion corresponding to (2.38) has the new drift: , independently of the sign of . It models the WF diffusion with selection conditioned on exit at . Assume that now solves if with boundary conditions . In this case study, one selects sample paths of with a large mean absorption time . Sample paths with large sojourn time in are favored. We have where is the Green function (2.22). The boundaries of are now both entrance boundaries and so . is not absorbed at the boundaries. The stopping time of is just its killing time . Let . Then, solves , , with explicit solution Assume that now solves if with boundary conditions . In this case study, one selects sample paths of with a large sojourn time density at recalling . The stopping time of occurs at rate . It is a killing time when the process is at for the last time after a geometrically distributed number of passages there with rate (or with success probability ). Let . Then, solves , with explicit solution when may be viewed as the age of a mutant currently observed to present frequency , see [18]. The Green function at of the transformed process is solution to . It takes the simple form Let be the smallest non-null eigenvalue of the infinitesimal generator . Let be the corresponding eigenvector, that is, satisfying, with boundary conditions . Then, . The new KB operator associated to the transformed process is obtained while killing the sample paths of the process governed by at constant death rate . The transition probability of the transformed stochastic process is Define . It is the transition probability of the process governed by ; it corresponds to the original process conditioned on never hitting the boundaries (the so-called process of , see [19]). It is simply obtained from by adding the additional drift term to , where is the eigenvector of associated to its smallest non-null eigenvalue. The determination of is a Sturm-Liouville problem. When is large, to the dominant order where is the Yaglom limit law of . Therefore Thus, the limit law of the process is the normalized Hadamard product of the eigenvectors and associated, respectively, to and . On the other hand, the limit law of is directly given by where is the appropriate normalizing constant. Comparing (2.77) and (2.78) The eigenvector associated to is, therefore, equal to the eigenvector associated to times the speed density of . When dealing for example with the neutral Wright-Fisher diffusion, it is known that with and (see Section 4.3, example ()). The process in this case obeys with the stabilizing drift toward : .The limit law of the process in this case is . The latter conditioning is more stringent than the Yaglom conditioning and so the limiting law has more mass away from the boundaries (compare with the uniform Yaglom limit). For additional similar examples in the context of WF diffusions and related ones, see [20].

2.7. Branching and the Reciprocal Doob Transform

Clearly, starting from the killed diffusion process with infinitesimal generator and applying the reciprocal Doob transform defined by leads to . Indeed, because Note that .

This suggests that starting from a diffusion process with infinitesimal generator (without its killing part) and applying the reciprocal Doob transform , one ends up with a modified process whose infinitesimal generator is where and is now a pure birth rate. Note that is now a subharmonic function for because .

Let . Because , we have and so is harmonic for Doob-transforming using , we get which is the infinitesimal generator of the original diffusion process.

Under our assumptions, both process and with infinitesimal generators and are nonconservative diffusion processes with absorbing barriers. Further, is bounded above. Therefore, may be written as where and .

The process with infinitesimal generator is now a pure binary branching diffusion process. For this class of models, an initial particle started at obeys a diffusion process with infinitesimal generator , absorbed when it hits the boundaries. At some random (mean ) exponential time, this particle dies, giving birth in the process to a random number (either 1 or 2) of daughter particles started where the mother particle died and diffusing independently as their mother did and so forth for the subsequent generation particles. We have .

The process with infinitesimal generator is, thus, a branching diffusion with supercritical binary splitting mechanism (). There is, therefore, a competition between the branching phenomenon that leads to an exponential increase of the number of particles in the system and the absorption at the boundaries of the living particles.

Let be the global number of particles which are alive in the system at each time , descending from an Eve particle started at , and let be the global extinction time of the population. Under our assumptions, this branching model fits to the general formalism for branching diffusion developed in ([9, 10]) from which we conclude uniformly in . This means the global extinction of the particle system under concern: In the tradeoff between branching and absorption at the boundaries, the system gets eventually extinct with probability 1 in finite time. We shall develop a typical example arising in population genetics in the subsequent sections.

3. The Wright-Fisher Example

In this section, we briefly and informally recall that the celebrated WF diffusion process with or without a drift may be viewed as a scaling limit of a simple two alleles discrete space-time branching process preserving the total number of individuals in the subsequent generations (see [8, 12, 21] for example).

3.1. The Neutral Wright-Fisher Model

Consider a discrete-time Galton Watson branching process preserving the total number of individuals in each generation. We start with individuals. The initial reproduction law is defined as follows: let and be integers. Assume that the first-generation random offspring numbers admit the following joint exchangeable polynomial distribution on the discrete simplex : This distribution can be obtained by conditioning independent Poisson distributed random variables on summing to . Assume subsequent iterations of this reproduction law are independent so that the population is with constant size for all generations.

Let be the offspring number of the first individuals at the discrete generation corresponding to (say) allele (the remaining number counts the number of alleles at generation ). This sibship process is a discrete-time Markov chain with binomial transition probability given by Assume next that , where . Then, as well known, the dynamics of the continuous space-time rescaled process , can be approximated for large , to the leading term in , by a Wright-Fisher-Itô diffusion on (the purely random genetic drift case) Here, is a standard Wiener process. For this scaling limit process, a unit laps of time corresponds to a laps of time for the original discrete-time process, thus time is measured in units of . If the initial condition is is the diffusion approximation of the offspring frequency of a singleton at generation .

Equation (3.3) is a one-dimensional diffusion as in (2.1) on , with zero drift and volatility . This diffusion is already in natural coordinate, and so . The scale function is and the speed measure . One can check that both boundaries are exit in this case: the stopping time is where is the extinction time and the fixation time. The corresponding infinitesimal generators are and .

3.2. Nonneutral Cases

Two alleles Wright-Fisher models (with non-null drifts) can be obtained by considering the binomial transition probabilities bin where is now some state-dependent probability (which is different from the identity ) reflecting some deterministic evolutionary drift from the allele to the allele . For each , we have which is amenable to a diffusion approximation in terms of , under suitable conditions.

For instance, taking , where are small (-dependent) mutation probabilities from to ( to .). Assuming that , leads after scaling to the drift of WF model with positive mutations rates .

Taking where are small dependent selection parameter satisfying , leads, after scaling, to the WF model with selective drift , where . Essentially, the drift is a large approximation of the bias: . The WF diffusion with selection is thus where time is measured in units of . Letting define a new time scale with inverse , the time-changed process now obeys the SDE with a small diffusion term. Here, and time is the usual time clock.

The WF diffusion with selection (3.8) tends to drift to (.) if allele is selectively advantageous over (.) in the following sense: if ( .), the fixation probability at , which is [15] increases (decreases) with taking larger (smaller) values. Putting , the fixation probability at 1 of an allele mutant is of order: ; see [15].

4. The Neutral WF Model

In this section, we particularize the general ideas developed in the introductory Section 2 to the neutral WF diffusion (3.3) and draw some straightforward conclusions most of which are known which illustrate the use of Doob transforms.

4.1. Explicit Solutions of the Neutral KBE and KFE

As shown by Kimura in [22], it turns out that both Kolmogorov equations are exactly solvable, in this case, using spectral theory. Indeed, the solutions involve a series expansion in terms of eigenfunctions of the KB and KF infinitesimal generators with discrete eigenvalues spectrum. We now consider the specific neutral WF model.

With , let be the degree- Gegenbauer polynomials solving with ; we let . When , we have and so , where is a polynomial with degree satisfying and . With , let be defined by: . These polynomials clearly constitute a system of eigenfunctions for the KB operator with eigenvalues , thus with . In particular, , , ,  ,  . With , we have and and .

The eigenfunctions of the KF operator are given by , where the Radon measure of weights is the speed measure: , for the same eigenvalues. For instance, , ,  ,  .

Although really constitutes an eigenvalue, only is not a polynomial. When , from their definition, the polynomials satisfy in such a way that is a polynomial with degree .

Let . We note that if and the system is a complete orthogonal set of eigenvectors. Therefore, for any square-integrable function admitting a decomposition in the basis . where . This series expansion solves the KBE: ; where .

Moreover, the transition probability density of the neutral WF models admits the spectral expansion Starting from , the cumulated probability masses by time at the exit boundaries are, respectively, (see [17]) which tend as toward the extinction and fixation probabilities, namely, here and . Because and , we get the identities leading to the relationship .

The series expansion for solves the KFE of the WF model. The transition density is reversible with respect to the speed density since for The measures , are not probability measures because the are not necessarily positive over . This decomposition is not a mixture. We have the 2-norm for the weight function . We notice that so that although is indeed an eigenvalue, the above sums should be started at (expressing the lack of an invariant measure for the WF model as a result of its absorption at the boundaries).

We have and so is the exact tail distribution of the absorption time.

Since , to the leading order in , for large time which is independent of . Integrating over , so that the conditional probability is asymptotically uniform in the Yaglom limit. As time passes by, given absorption did not occur in the past, (as ) which is a uniformly distributed random variable on .

4.2. Additive Functionals for the Neutral WF

Let be the WF diffusion model defined by (3.3) on the interval , where both endpoints are absorbing (exit). We wish to evaluate the additive quantities where functions and are both nonnegative. With , solves Take and , when in this case, is the Green function. The solution takes the simple form The Green function solves the above general problem of evaluating additive functionals

As a Few Examples(1) Let and here, is the mean time of absorption (average time spent in before absorption). The solution is (the Crow and Kimura formula, see [2]) (2) Let and . Let . Then, is a harmonic function solution to , with boundary conditions and . The solution for WF model is: . Stated differently, is the probability that the exit time at is less than the one at , starting from . On the contrary, choosing to be a harmonic function with boundary conditions and , . Thus, .(3) Let measure the heterozygosity of the WF process at time and assume . A remarkable thing is that the average heterozygosity over the sample paths is which is the initial heterozygosity of the population.

4.3. Transformation of WF Sample Paths, [3]

With the transition probability density of WF model, define a new transformed stochastic process by its transition probability

Conditioning WF on exit at some boundary. Assume first solves with boundary conditions and ; hence, reads . In this case, (no killing), and so is the absorption time for a process governed by a new SDE with a drift term. The new process is just conditioned on exiting at . The boundary 1 is exit whereas 0 is entrance. Thus, the model for becomes , now with linear drift and . Its transition probability is where the subscript 1 indicates that this is the conditional transition probability of sample paths whose exit is necessarily at the boundary 1.

Assuming now solves if with boundary conditions and , the new process is just conditioned on exiting at . Boundary 0 is exit, whereas 1 is entrance; in this case, is . Thus, the model for becomes , with and . Its transition probability is where the subscript 0 indicates that this is the conditional transition probability of WF sample paths whose exit now is at . Recalling that, starting from , gets absorbed at (.) with probability (.), we recover that Using the solution to KFE for , we obtain an expression for both and , simply by premultiplying it by the corresponding right factor. Integrating the results over , we get the conditional tail distributions of the exit times at or 0, given the exit is at or 0.

Exploiting the large time behavior of , to the first order in , we get Integrating over , and are the large time behaviors of the absorption times at 1 and 0, respectively. Using this, we get the large time behaviors of the conditional probabilities where we recognize the densities of specific beta-distributed random variables. Specifically, we conclude that as time passes by, given absorption occurs at and given it has not occurred in the past, beta distribution on . Similarly, given absorption occurs at and given it has not occurred previously, beta distribution on .

In the previously displayed formula, is just the exit time at (at .) of the conditional transformed WF diffusions. Let . Then, with , solves , whose explicit solution is in terms of , the Green function of . For the WF model conditioned on exit at (.), we find, respectively, the Kimura and Ohta's formulae in [23] This result could have been guessed by observing that is the expected absorption time of the original WF model. When , (resp., ), it takes an average time 2 to reach 1 (.) for the WF model conditioned on exit at (.).

Selection of WF sample paths with large heterozygosity. Assume that now solves if with boundary conditions . Then, and this is the right eigenvector of associated to the smallest positive eigenvalue of the neutral WF model. In this case study, one selects sample paths of with large heterozygosity. The dynamics of in (2.38) is subject to a constant killing rate 1. The boundaries of are now both entrance boundaries and so . is not absorbed at the boundaries. The stopping time of is just its killing time which is mean 1 exponentially distributed, independently of the starting point . Indeed, recalling and observing if .

As time passes, killing of occurs, and given killing will never occur in the future, a random variable with density on which is a beta() density. In this selection of paths procedure, the conditional density of given is indeed , where . Therefore, Recalling and , we get , regardless of the initial condition . This is the beta limit law of the -process of the neutral WF diffusion.

Selection of WF sample paths with large sojourn time density at . Assume now that solves if and so . Using the Green function of the neutral WF model, the transition probability density of is Thus, given (), coincides with conditioned to exit in 1 (.) killed at rate when it passes through , necessarily at some time.

The stopping time of is just its killing time when the process is at for the last time with a geometrically number of passages at with rate 1 (or success probability ).

5. The WF Model with Selection

Now, we focus on the diffusion process (3.8). Let be the Gegenbauer eigen-polynomials of the KF operator corresponding to the neutral WF diffusion (3.3) and so with eigenvalues . Define the oblate spheroidal wave functions on as where obey the three-term recurrence defined in [24]. In the latter equality, the summation is over odd (even) values if is even (odd).

Define and where is the speed measure density of the WF model with selection (3.8).

The system constitute a system of eigenfunctions for the WF with selection generators and with eigenvalues implicitly defined in [24], thus with and . The eigenfunction expansion of the transition probability density of the WF model with selection is thus, [25], where . The WF model with selection can be viewed as a perturbation problem of the neutral WF model (see [3]). There exist perturbation developments of around with respect to , [25]. They are valid and useful for small .

The WF diffusion process with selection (3.8) is nonconservative, with finite hitting time of one of the boundaries. Following the general arguments developed in Section 2, the Yaglom limit of conditioned on is the normalized version of The limit law of conditioned on never hitting the boundaries in the remote future is the normalized version of Because the latter conditioning is more stringent than the former, the probability mass of (5.4) is more concentrated inside the interval than (5.3).

6. From the WF Model with Selection to the Neutral WF Model: Doob Transform and Killing

We shall consider the following transformation of paths for the WF model with selection. Consider the Wright-Fisher diffusion with selection (3.8): , . For this model, and both boundaries are exit.

Assume that so that the drift term is bounded above by , together with being bounded below (as a constant function here equal to ). We are then in the general framework of the problems under study in this paper. This suggests that for some admissible choice of a superharmonic exponential function , the -Doob transform of could lead to a transformed process with bounded killing rate . We shall choose for its interesting features.

The transition density of admits the representation (5.2) in terms of their oblate spheroidal wave eigenfunctions. Let where is the skewed sample heterozygosity, damped by the factor . Then, solves , with solution . In this case study, one selects sample paths of with large . The dynamics of is the drift-less neutral WF dynamics , subject to quadratic killing at rate in , which is bounded above there. The boundaries of are still exit and the stopping time of is , where is its absorption time at the boundaries and its killing time. The density of the transformed process is . Its series expansion is exactly known using (5.2) for .

The transformed process accounts for a neutral evolution of the allele frequency subject to the additional extinction opportunity of the population itself due to killing at rate proportional to its heterozygosity. Leaving aside the fact that it can be obtained after a suitable Doob transformation, this model is of importance in population genetics: it first appeared in ([8, Page 272]) as a scaling limit of a population genetics model of recombination.

From the general study of Section 2, we obtain the following.

Conditioned on , the transformed process admits a Yaglom limit . With the first oblate spheroidal eigenvector of associated to the smallest positive eigenvalue , is of the product form This limiting probability is the Yaglom limit law of conditioned on the event that both the absorption and killing times exceed .

Let be the infinitesimal generator of with two stopping times. Now, there is a tradeoff between which of and occurs first. To solve it, we need to compute defined in (2.55) by , with boundary conditions . This is a Sturm-Liouville problem whose solution in our case is The function is minimal when , with value . This looks natural because when , the chance to hit before getting killed should be the lowest.

The density of the process conditioned on the event is and so is also explicitly known from the oblate spheroidal wave expansion (5.2) of . The tail distribution of given is obtained by integrating over .

Similarly, the density of the process conditioned on the event is The tail distribution of given is obtained by integrating over .

The associated conditioned on absorption first process obeys the SDE with drift and local variance unchanged . This process has no killing part and it gets eventually absorbed at .

In the generator of the conditioned on killing first process, there is a killing multiplicative part which is enhanced and a shift in the drift, showing that the associated conditioned process exhibits a faster killing rate, but the drift shift guarantees that is not absorbed at the boundaries. With and as in (6.3), the drift takes the peculiar explicit form

7. From the Neutral WF Model to the WF Model with Selection: Reciprocal Doob Transform and Branching

We now follow the general path indicated in Section 2.7 and apply it to the particular models under concern. We, therefore, illustrate and develop the idea of a reciprocal Doob transform on the specific example of interest.

The starting point is now the neutral Wright-Fisher diffusion: , . For this model, and both boundaries are exit. Its transition density admits the representation in terms of the Gegenbauer eigenpolynomials (see Section 4.1). We shall consider the following reciprocal transformation of paths for the neutral WF model: let and consider . We now have and .

In this case study, one selects sample paths of with large . The dynamics of is easily seen to be the WF with selection dynamics subject to quadratic branching at rate inside . We indeed have where is the KBE operator of the dynamics .

With , we clearly have and is an harmonic function for and as a result, Doob-transforming by , we get which is the infinitesimal generator of the original neutral WF martingale.

The birth (creating) rate in is bounded from above on . It may be put into the canonical form , where and whose range is the interval as .

The density of the transformed process is . It is exactly known because is known from (7.1).

The transformed process (with infinitesimal backward generator ) accounts for a branching diffusion (BD), where a diffusing mother particle (with generator and started at ) lives a random exponential time with constant rate . When the mother particle dies, it gives birth to a spatially dependent random number of particles (with mean ). If , independent daughter particles are started where their mother particle died; they move along a WF diffusion with selection and reproduce, independently and so on.

Because is bounded above by 2 and larger than 1 (indicating a supercritical branching process), we actually get a BD with binary scission whose random offspring number satisfies with (the event that 2 particles are generated in a splitting event is more probable than a single one).

For such a transformed process, the tradeoff is of a different nature: there is a competition between the boundaries which are still absorbing for the system of particles and the number of particles in the system at each time , which may grow due to branching events. The density of the transformed process has now the following interpretation: where is the density at of the th alive particle descending from the ancestral one (Eve), started at . In the latter formula, the sum vanishes if . A particle is alive at time if it came to birth before and has not been yet absorbed by the boundaries.

Let . Then, is the expected number of particle alive at time . We have But then, obeys the forward PDE as a result of . We have showing that is the average presence density at of the system of particles all descending from Eve started at .

Clearly, (and, therefore, also , because The expected number of particles in the system decays globally at rate .

The BD transformed process, therefore, admits an integrable Yaglom limit , solution to or . With , the first eigenvector of associated to the smallest positive eigenvalue , is of the product form This limiting probability is the Yaglom limiting average presence density at for the BD system of particles (it is also the ground state for ).

There is also a natural eigenvector of the backward operator , satisfying (the ground state for ). It is explicitly here that In the terminology of [26], both operators and its adjoint are critical ( () is said to be critical if there exists some function ( .), strictly positive in , such that: (.) and the operators do not possess a minimal positive Green function.). In this context, the constant is called the generalized principal eigenvalue. The eigenfunctions are their associated ground states.

We note that we have the -product property (see [26, Subsection 4.9]). With , let We have the condition We conclude (following [9, 10]) that, as a result of the condition (7.17) being trivially satisfied, global extinction holds in the following sense:

, uniformly in ,

there exists a constant , uniformly in ,

For all bounded measurable function on , From , it is clear that the process gets ultimately extinct with probability 1. In the tradeoff between branching and absorption at the boundaries, all particles get eventually absorbed and the global BD process turns out be subcritical (even though for all ): probability mass escapes out of although the BD survives with positive probability.

In the statement , the quantity is also where is the global extinction time of the particle system descending from an Eve particle started at . The number is the usual Malthus decay rate parameter. From has a natural interpretation in terms of the propensity of the particle system to survive to its extinction fate: the so-called reproductive value in demography.

with reads giving an interpretation of the constant (which may be hard to evaluate in practise).

The ground states of and its adjoint are, thus, and explicit here. It is useful to consider the process whose infinitesimal generator is given by the Doob transform because product-criticality is preserved under this transformation. The ground states associated to this new operator and its dual are . Developing, we obtain a process whose infinitesimal generator is with no multiplicative part. In our case study, we get adding a stabilizing drift towards to the original neutral WF model. The associated diffusion process is positive recurrent and so its invariant measure is integrable. It is the beta limit law of the -process (see (2.80) and (4.23)) for the neutral WF diffusion.

Remark 7.1. At time , let denote the positions of the BD particle system. Let stand for the functional generating function () of the measure-valued branching particle system. obeys the nonlinear (quadratic) Kolmogorov-Petrovsky-Piscounoff PDE, [27] where or is the shifted probability generating function of the branching law of . Thus, the nonlinear part reads , which is quadratic in .In particular, if , obeys the linear backward PDE involving . The latter evolution equation is the backward version of the forward PDE giving the evolution of as .

References

T. Nagylaki, “Gustave Malécot and the transition from classical to modern population genetics,” Genetics, vol. 122, no. 2, pp. 253–268, 1989.

J. F. Crow and M. Kimura, An Introduction to Population Genetics Theory, Harper & Row, New York, NY, USA, 1970.

T. Maruyama, Stochastic Problems in Population Genetics, vol. 17 of Lecture Notes in Biomathematics, Springer, Berlin, Germany, 1977.

W. J. Ewens, Mathematical Population Genetics: I. Theoretical Introduction, vol. 27 of Interdisciplinary Applied Mathematics, Springer, New York, NY, USA, 2nd edition, 2004.

R. Durrett, Probability Models for DNA Sequence Evolution, Probability and Its Applications, Springer, New York, NY, USA, 2nd edition, 2008.

J. H. Gillespie, The Causes of Molecular Evolution, Oxford University Press, New York, NY, USA, 1991.

A. M. Yaglom, “Certain limit theorems of the theory of branching random processes,” Doklady Akademii Nauk SSSR, vol. 56, pp. 795–798, 1947.

S. Karlin and H. M. Taylor, A Second Course in Stochastic Processes, Academic Press, New York, NY, USA, 1981.

S. Asmussen and H. Hering, “Strong limit theorems for general supercritical branching processes with applications to branching diffusions,” Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, vol. 36, no. 3, pp. 195–212, 1976.

S. Asmussen and H. Hering, “Some modified branching diffusion models,” Mathematical Biosciences, vol. 35, no. 3-4, pp. 281–299, 1977.

E. B. Dynkin, Markov Processes. Vols. I, II, vol. 122 of Die Grundlehren der Mathematischen Wissenschaften, Bände 121, Academic Press, New York, NY, USA; Springer, Berlin, Germany, 1965, translated with the Authorization and Assistance of the Author by J. Fabius, V. Greenberg, A. Maitra, G. Majone.

S. N. Ethier and T. G. Kurtz, Markov Processes: Characterization and Convergence, Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics, John Wiley & Sons, New York, NY, USA, 1986.

P. Mandl, Analytical Treatment of One-Dimensional Markov Processes, vol. 151 of Die Grundlehren der mathematischen Wissenschaften, Academia Publishing House of the Czechoslovak Academy of Sciences, Prague, Czech Republic; Springer, New York, NY, USA, 1968.

K. Itô, “On stochastic differential equations,” Memoirs of the American Mathematical Society, vol. 1951, no. 4, p. 51, 1951.

M. Kimura, “On the probability of fixation of mutant genes in a population,” Genetics, vol. 47, pp. 713–729, 1962.

D. Steinsaltz and S. N. Evans, “Quasistationary distributions for one-dimensional diffusions with killing,” Transactions of the American Mathematical Society, vol. 359, no. 3, pp. 1285–1324, 2007.

A. J. McKane and D. Waxman, “Singular solutions of the diffusion equation of population genetics,” Journal of Theoretical Biology, vol. 247, no. 4, pp. 849–858, 2007.

R. C. Griffiths, “The frequency spectrum of a mutation, and its age, in a general diffusion model,” Theoretical Population Biology, vol. 64, no. 2, pp. 241–251, 2003.

A. Lambert, “Population dynamics and random genealogies,” Stochastic Models, vol. 24, supplement 1, pp. 45–163, 2008.

T. Huillet, “On Wright-Fisher diffusion and its relatives,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2007, no. 11, Article ID P11006, 2007.

R. A. Blythe and A. J. McKane, “Stochastic models of evolution in genetics, ecology and linguistics,” Journal of Statistical Mechanics, vol. 7, Article ID P07018, 2007.

M. Kimura, “Stochastic processes and distribution of gene frequencies under natural selection,” Cold Spring Harbor Symposia on Quantitative Biology, vol. 20, pp. 33–53, 1955.

M. Kimura and T. Ohta, “The age of a neutral mutant persisting in a finite population,” Genetics, vol. 75, pp. 199–212, 1973.

S. Mano, “Duality, ancestral and diffusion processes in models with selection,” Theoretical Population Biology, vol. 75, no. 2-3, pp. 164–175, 2009.

M. Kimura, “Diffusion models in population genetics,” Journal of Applied Probability, vol. 1, pp. 177–232, 1964.

R. G. Pinsky, Positive Harmonic Functions and Diffusion, vol. 45 of Cambridge Studies in Advanced Mathematics, Cambridge University Press, Cambridge, UK, 1995.

A. Kolmogorov, I. Petrovsky, and N. Piscounov, “Étude de l'équation de la diffusion avec croissance de la quantité de matière et son application à un probléme biologique,” Moscow University Mathematics Bulletin, vol. 1, p. 125, 1937.