Nonconservative Diffusions on with Killing and Branching: Applications to Wright-Fisher Models with or without Selection
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.
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  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 . 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
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 . 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 ), 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 ).
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 .
(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 , 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,  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 ).
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, .
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 ).
The probability mass cumulated at the boundaries by time clearly is 
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 .
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 ). 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