# Nonconservative Diffusions on[0,1] with Killing and Branching: Applications to Wright-Fisher Models with or without Selection

Authors: Thierry E. Huillet
Journal: International Journal of Stochastic Analysis
 (2011)

---

## Abstract

We consider nonconservative diffusion processesxt 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 xt~ presenting an
additional killing rate part d>0. We limit ourselves to
situations for which xt~ 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 xt~, apply the reciprocal
Doob transform ending up in a new process which is xt but now with
an additional branching rate b>0, 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 xt is a Wright-Fisher diffusion with selection. Using an
exponential Doob transform, we end up with a killed neutral
Wright-Fisher diffusion xt~. We give a detailed study of the
binary branching diffusion process obtained by using the
corresponding reciprocal Doob transform.


---

## Body

## 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.Section2 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 xt, 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 processesxt 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 x̃t presenting an additional killing rate part d>0. We further limit ourselves to situations for which x̃t 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 d. 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 processx̃t, apply the reciprocal Doob transform ending up in a new process which is xt but now with an additional branching rate b>0, 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 Section3 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 Section6, we apply the Doob-transformation techniques to these processes: The starting point process xt is a Wright-Fisher diffusion with selection differential σ>0. We use the exponential Doob kernel α=e-σx. The transformed process accounts for a neutral Wright-Fisher evolution x̃t for the allele 1 frequency, subject to the additional possibility of the extinction of the population itself due to killing at rate d 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 xt and x̃t which are known (from the works of Kimura) to involve oblate spheroidal wave functions and Gegenbauer polynomials, respectively.In Section7, 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 eσx 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[0,1]

Let(wt;t≥0) be a standard one-dimensional Brownian (Wiener) motion. Consider a one-dimensional Itô diffusion driven by (wt;t≥0) on the interval say [0,1]; see [14]. We will let I=(0,1). Assume that it has locally Lipschitz continuous drift (x) and local standard deviation (volatility) g(x), namely, consider the stochastic differential equation (SDE)(2.1)dxt=f(xt)dt+g(xt)dwt,x0=x∈(0,1).
The condition on f(x) and g(x) guarantees in particular that there is no point x* in I for which |f(x)| or |g(x)| would blow up and diverge as |x-x*|→0.The Kolmogorov backward infinitesimal generator of (2.1) is G=f(x)∂x+(1/2)g2(x)∂x2. As a result, for all suitable ψ in the domain of the operator St:=etG, u:=u(x,t)=Eψ(xt∧τx) satisfies the Kolmogorov backward equation (KBE) (2.2)∂tu=G(u);u(x,0)=ψ(x).
In the definition of the mathematical expectation u, we have t∧τx:=inf(t,τx), where τx indicates a random time at which the process should possibly be stopped (absorbed), given the process was started in x. The description of this (adapted) absorption time is governed by the type of boundaries which {0,1} are to (xt;t≥0).


### 2.2. Natural Coordinate, Scale, and Speed Measure

For such Markovian diffusions, it is interesting to consider theG-harmonic coordinate φ∈C2 belonging to the kernel of G that is, satisfying G(φ)=0. For φ and its derivative φ′:=dφ/dy, with (x0,y0)∈(0,1), one finds (2.3)φ′(y)=φ′(y0)e-2∫y0y(f(z)/(g2(z)))dz,φ(x)=φ(x0)+φ′(y0)∫x0xe-2∫y0y(f(z)/(g2(z)))dzdy.
One should choose a version of φ satisfying φ′(y)>0, y∈I. The function φ kills the drift f of (xt;t≥0) in the sense that considering the change of variable yt=φ(xt),(2.4)dyt=(φ′g)(φ-1(yt))dwt,y0=φ(x).
The driftless diffusion (yt;t≥0) is often termed the diffusion in natural coordinates with state-space [φ(0),φ(1)]. Its volatility is g̃(y):=(φ′g)(φ-1(y)). The function φ is often called the scale function.Wheneverφ(0)>-∞ and φ(1)<+∞, one can choose the integration constants defining φ(x) so that (2.5)φ(x)=∫0xe-2∫0y(f(z)/(g2(z)))dzdy∫01e-2∫0y(f(z)/(g2(z)))dzdy,
with φ(0)=0 and φ(1)=1. In this case, the state-space of (yt;t≥0) is again [0,1], the same as for (xt;t≥0).Finally, considering the random time changet→θt with inverse: θ→tθ defined by θtθ=θ and (2.6)θ=∫0tθg̃2(ys)ds,
the novel diffusion (wθ:=ytθ;θ≥0) is easily checked to be identical in law to a standard Brownian motion. Let now δy(·)= weak-limε↓0(1/2ε)1(·∈(y-ε,y+ε)) stand for the Dirac delta mass at y. The random time θt can be expressed as (2.7)θt=∫01dx⋅m(x)∫0tδφ(x)(ws)ds=∫0tm(φ-1(ws))ds,
where m(x):=1/(g2φ′)(x) is the (positive) speed density at x=φ-1(y) and Lt(y):=limε↓0(1/2ε)∫0t1(ws∈(y-ε,y+ε))ds the local time at y of the Brownian motion before time t. Both the scale function φ and the speed measure dμ=m(x)·dx are, therefore, essential ingredients to reduce the original stochastic process (xt;t≥0) to the standard Brownian motion (wt;t≥0). Indeed, it follows from the above arguments that if θt=∫0tm(xs)ds, then (φ(xθt);t≥0) is a Brownian motion. The Kolmogorov backward infinitesimal generator G may be written in Feller form (2.8)G(⋅)=12ddμ(ddφ⋅).Examples 2.2 (from population genetics).
(i) Assume thatf(x)=0 and g2(x)=x(1-x). This is the neutral Wright-Fisher (WF) model discussed at length later. This diffusion is already in natural scale and φ(x)=x, m(x)=[x(1-x)]-1. The speed measure is not integrable.
(ii) Withu1,u2>0, assume f(x)=u1-(u1+u2)x and g2(x)=x(1-x). This is the Wright-Fisher model with mutation. The parameters u1,u2 can be interpreted as mutation rates. The drift vanishes when x=u1/(u1+u2) which is an attracting point for the dynamics. Here,φ′(y)=φ′(y0)y-2u1(1-y)-2u2,φ(x)=φ(x0)+φ′(y0)∫x0xy-2u1(1-y)-2u2dy, with φ(0)=-∞ and φ(1)=+∞ if u1,u2>1/2. The speed measure density is m(x)∝x2u1-1(1-x)2u2-1 and so is always integrable.
(iii) Withσ∈R, assume a model with quadratic logistic drift f(x)=σx(1-x) and local variance g2(x)=x(1-x). This is the WF model with selection. For this diffusion (see [15]), φ(x)=((1-e-2σx)/(1-e-2σ)) and m(x)∝[x(1-x)]-1e2σx 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 thatf(x) and g(x) are now differentiable in I. Let then p(x;t,y) stand for the transition probability density function of xt at y given x0=x. Then, p:=p(x;t,y) is the smallest solution to the Kolmogorov forward (Fokker-Planck) equation (KFE) (2.9)∂tp=G*(p),p(x;0,y)=δy(x),
where G*(·)=-∂y(f(y)·)+(1/2)∂y2(g2(y)·) is the adjoint of G (G* acts on the terminal value y, whereas G acts on the initial value x). The way one can view this PDE depends on the type of boundaries that {0,1} are.We will next suppose that the boundaries∘:=0 or 1 are both exit (or absorbing) boundaries. From the Feller classification of boundaries, this will be the case if forally0∈(0,1)(2.10)(i)m(y)∉L1(y0,∘),(ii)φ′(y)∫y0ym(z)dz∈L1(y0,∘),
where a function f(y)∈L1(y0,∘) if ∫y0∘|f(y)|dy<+∞.In this case, a sample path of(xt;t≥0) can reach ∘ from the inside of I in finite time but cannot reenter. The sample paths are absorbed at ∘. There is an absorption at ∘ at time τx,∘=inf(t>0:xt=∘|x0=x) and P(τx,∘<∞)=1. Whenever both boundaries {0,1} are absorbing, the diffusion xt should be stopped at τx:=τx,0∧τx,1. Would none of the boundaries {0,1} be absorbing, then τx=+∞, 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,p(x;t,y) is a subprobability. Letting ρt(x):=∫01p(x;t,y)dy, we clearly have ρt(x)=P(τx>t). Such models are nonconservative.For one-dimensional diffusions, the transition densityp(x;t,y) is reversible with respect to the speed density ([8, Chapter 15, Section  13]) and so detailed balance holds (2.11)m(x)p(x;t,y)=m(y)p(y;t,x),0<x,y<1.
The speed density m(y) satisfies G*(m)=0. It may be written as a Gibbs measure with density: m(y)∝(1/g2(y))e-U(y), where the potential function U(y) reads (2.12)U(y):=-2∫0yf(z)g2(z)dz,0<y<1,
and with the measure dy/g2(y) standing for the reference measure.Further, ifp(s,x;t,y) is the transition probability density from (s,x) to (t,y), s<t, then -∂sp=G(p), with terminal condition p(t,x;t,y)=δy(x) and so p(s,x;t,y) 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 p:=p(x;t,y), operating the time substitution t-s→t, p itself solves the KBE (2.13)∂tp=G(p),p(x;0,y)=δy(x).
In particular, integrating over y, ∂tρt(x)=G(ρt(x)), with ρ0(x)=1(x∈(0,1)).p(x;t,y) being a sub-probability, we may define the normalized conditional probability density q(x;t,y):=p(x;t,y)/ρt(x), now with total mass 1. We get (2.14)∂tq=-∂tρt(x)ρt(x)⋅q+G*(q),q(x;0,y)=δy(x).
The term bt(x):=-∂tρt(x)/ρt(x)>0 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 (xt;t≥0) at the boundaries. In this creation of mass process, a diffusing particle started in x dies at rate bt(x) at point (t,y), where it is duplicated in two new independent particles both started at y (resulting in a global birth) evolving in the same diffusive way  (consider a diffusion process with forward infinitesimal generator G* governing the evolution of p(x;t,y). 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 d and b depend on the current location y of the path. Then, the process with the birth and death opportunities of a path has the infinitesimal generator λ(y)·+G*(·), where λ(y)=b(y)-d(y). The rate can also depend on t and x).  The birth rate function bt(x) depends here on x and t, not on y.When the boundaries ofxt are absorbing, the spectra of both -G and -G* are discrete (see [8, Page 330]): There exist positive eigenvalues (λk)k≥1 ordered in ascending sizes and eigenvectors (vk,uk)k≥1 of both -G* and -G satisfying -G*(vk)=λkvk and -G(yk)=λkuk such that  with 〈uk,vk〉:=∫01uk(x)vk(x)dx and bk:=〈uk,vk〉-1, the spectral representation (2.15)p(x;t,y)=∑k≥1bke-λktuk(x)vk(y)
holds.Letλ1>λ0=0 be the smallest nonnull eigenvalue of the infinitesimal generator -G* (and of -G). Clearly, -(1/t)logρt(x)→t→∞λ1 and by L' Hospital rule, therefore, bt(x)→t→∞λ1. Putting ∂tq=0 in the latter evolution equation, independently of the initial condition x(2.16)q(x;t,y)⟶t→∞q∞(y)=v1(y),
where v1 is the eigenvector of -G* associated to λ1, satisfying -G*v1=λ1v1. The limiting probability v1/norm (after a proper normalization) is called the Yaglom limit law of (xt;t≥0) conditioned on being currently alive at all time t (see [7]).


### 2.4. Additive Functionals Along Sample Paths

Let(xt;t≥0) be the diffusion model defined by (2.1) on the interval I, where both endpoints are assumed absorbing (exit). This process is, thus, transient and nonconservative. We wish to evaluate the nonnegative additive quantities (2.17)α(x)=E(∫0τxc(xs)ds+d(xτx)),
where the functions c and d are both assumed nonnegative on I and ∂I={0,1}. The functional α(x)≥0 solves the Dirichlet problem (2.18)-G(α)=cifx∈I,α=difx∈∂I,
and α is a superharmonic function for G, satisfying -G(α)≥0.Some Examples.
(1) Assume thatc=1 and d=0: here, α=E(τx) is the mean time of absorption (average time spent in (0,1) before absorption), solution to
(2.19)-G(α)=1ifx∈I,α=0ifx∈∂I.
(2) Whenever both{0,1} are exit boundaries, it is of interest to evaluate the probability that xt first hits [0,1] (say) at 1, given x0=x. This can be obtained by choosing c=0 and d(∘)=1(∘=1).
Let thenα=:α1(x)=P(xtfirsthits[0,1]at1∣x0=x).α1(x) is a G-harmonic function solution to G(α1)=0, with boundary conditions α1(0)=0 and α1(1)=1. Solving this problem, we get(2.20)α1(x)=φ(x)-φ(0)φ(1)-φ(0)=∫0xdye-2∫0y(f(z)/(g2(z)))dz∫01dye-2∫0y(f(z)/(g2(z)))dz.
On the contrary, choosing α0(x) to be a G-harmonic function with boundary conditions α0(0)=1 and α0(1)=0, α0(x)=P(xtfirsthits[0,1]at0∣x0=x)=1-α1(x).
(3) Lety∈I and put c=(1/2ε)1(x∈(y-ε,y+ε)) and d=0. As ε→0, c converges weakly to δy(x) and, α=:𝔤(x,y)=E(lim(1/2ε)∫0τx1(xs∈(y-ε,y+ε))ds)=∫0∞p(x;s,y)ds is the Green function, solution to(2.21)-G(g)=δy(x)ifx∈I,g=0ifx∈∂I.𝔤 is, therefore, the mathematical expectation of the local time at y, starting from x (the sojourn time density at y). The solution is known to be (see [8, page 198] or [5, page 280])(2.22)g(x,y)=2(φ(x)-φ(0))(φ(1)-φ(y))(g2φ′)(y)(φ(1)-φ(0))ifx≤y,g(x,y)=2(φ(1)-φ(x))(φ(y)-φ(0))(g2φ′)(y)(φ(1)-φ(0))ifx≥y.
The Green function is of particular interest to solve the general problem of evaluating additive functionals α(x). Indeed, as is well known, see [8], for example, the integral operator with respect to the Green kernel inverts the second-order operator -G leading to(2.23)α(x)=∫Ig(x,y)c(y)dyifx∈I,α=difx∈∂I.
Under this form, α(x) appears as a potential function and all potential function is superharmonic. Note that for all harmonic function h≥0 satisfying -G(h)=0,
(2.24)αh(x):=∫Ig(x,y)c(y)dy+h(x)
is again superharmonic because -G(αh)=c≥0.
(4) Also of interest are the additive functionals of the type(2.25)αλ(x)=E(∫0τxe-λsc(xs)ds+d(xτx)),
where the functions c and d are again both assumed to be nonnegative. The functional αλ(x)≥0 solves the Dynkin problem, [8] (2.26)(λI-G)(αλ)=cifx∈I,αλ=difx∈∂I
involving the action of the resolvent operator (λI-G)-1 on c.
Wheneverc(x)=δy(x),d=0, then(2.27)αλ=:gλ(x,y)=E(∫0τxe-λsδy(xs)ds)=∫0∞e-λsp(x;s,y)ds
is the λ-potential function, solution to
(2.28)(λI-G)(gλ)=δy(x)ifx∈I,gλ=0ifx∈∂I.𝔤λ is, therefore, the mathematical expectation of the exponentially damped local time at y, starting from x (the temporal Laplace transform of the transition probability density from x to y at t), with g0=g. Then, it holds that
(2.29)αλ(x)=∫Igλ(x,y)c(y)dyifx∈I,αλ=difx∈∂I.
The λ-potential function is also useful in the computation of the distribution of the first-passage time τx,y to y starting from x. From the convolution formula,
(2.30)p(x;t,y)=∫0tP(τx,y∈ds)p(y;t-s,y),
and taking the Laplace transform of both sides with respect to time, we obtain the Laplace-Stieltjes transform (LST) of the law of τx,y as
(2.31)E(e-λτx,y)=gλ(x,y)gλ(y,y).
We have P(τx,y<∞)=𝔤(x,y)/(𝔤(y,y))∈(0,1) as a result of both terms in the ratio being finite and x,y belonging to the same transience class of the process (under our assumptions that the boundaries are absorbing). Note that from the reversibility property
(2.32)m(x)g(x,y)=m(y)g(y,x).


### 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(xt;t≥0) as in (2.1) with absorbing barriers. Let p(x;t,y) be its transition probability, and let τx be its absorption time at the boundaries.Letα(x):=E(∫0τxc(xs)ds+d(xτx)) be a nonnegative additive functional solving (2.33)-G(α)=cifx∈I,α=difx∈∂I.
Recall the functions c and d are both chosen nonnegative so that so is α.Define a new transformed stochastic process(x¯t;t≥0) by its transition probability (2.34)p¯(x;t,y)=α(y)α(x)p(x;t,y).
In this construction of (x¯t;t≥0) through a change of measure, sample paths of (xt;t≥0) for which α(y) is large are favored. This is a selection of paths procedure due to Doob (see [11]).Now, the KFE forp¯ clearly is ∂tp¯=G¯*(p¯), with p(x;0,y)=δy(x) and G¯*(p¯)=α(y)G*(p¯/α(y)). The Kolmogorov backward operator of the transformed process is, therefore, by duality (2.35)G¯(⋅)=1α(x)G(α(x)⋅).
Developing, with α′(x):=dα(x)/dx and G̃(·):=(α′/α)g2∂x(·)+G(·), we get (2.36)G¯(⋅)=1αG(α)⋅+G̃(⋅)=-cα⋅+G̃(⋅),
and the new KB operator can be obtained from the latter by adding a drift term (α′/α)g2∂x to the one in G of the original process to form a new process (x̃t;t≥0) with the KB operator G̃ and by killing its sample paths at death rate d(x):=(c/α)(x) (provided c≠0). Note that (2.37)d(x)=1(1/(c(x)))∫01g(x,y)c(y)dy.
In others words, with f̃(x):=f(x)+(α′/α)g2(x), the novel time-homogeneous SDE to consider is (2.38)dx̃t=f̃(x̃t)dt+g(x̃t)dwt,x̃0=x∈(0,1),
possibly killed at rate d=c/α as soon as c≠0. Whenever (x̃t;t≥0) is killed, it enters conventionally into some coffin state {∂} added to the state-space. Let τ̃x be the new absorption time at the boundaries of (x̃t;t≥0) started at x (with τ̃x=∞ would the boundaries be inaccessible to the new process x̃t which we ruled out). Let τ̃x,∂ be the killing time of (x̃t;t≥0) started at x (the hitting time of ∂), with τ̃x,∂=∞ if c=0. Then, τ¯x:=τ̃x∧τ̃x,∂ is the novel stopping time of (x̃t;t≥0). The SDE for (x̃t;t≥0), together with its global stopping time τ¯x characterize the new process (x¯t;t≥0) with generator G¯ to consider.In the sequel, we shall limit ourselves to the cases for which the following additional conditions hold on the transformed process.(i) Nonconservativeness of x̃t.
We will next suppose that the boundaries∘:=0 or 1 are both exit (or absorbing) boundaries for the new process x̃t in (2.38). From the Feller criterion for exit boundaries, this will be the case if forally0∈(0,1)(2.39)m̃(y)∉L1(y0,∘),φ̃′(y)∫y0ym̃(z)dz∈L1(y0,∘),
where m̃(y)=1/((g2φ̃′)(y)) is the new speed measure density for x̃t and φ̃ its scale function. Recalling f̃=f+(α′/α)g2 and g̃2=g2, we have
(2.40)φ̃′(y)=φ̃′(y0)e-2∫y0y(f̃(z)/(g2(z)))dz=φ̃′(y0)α2(y)φ′(y),φ̃(x)=φ̃(x0)+φ̃′(y0)∫x0xα-2(y)e-2∫y0y(f(z)/(g2(z)))dzdy.
So, we assume here that x̃t obeys itself a nonconservative diffusion.(ii) Boundedness of the Killing Rate d.
In some examples, the killing rated=-G(α)/α is bounded above. For example, suppose that the drift of the diffusion process (xt;t≥0) is bounded above by f*=maxx(f(x))>0. (If the drift of (xt;t≥0) is bounded below by f*<0, we are led to the same conclusions while considering the process 1-xt instead of xt.) Then, choosing α(x)=e-ax,a>0,-G(α)=(af-(a2/2)g2)α<af*α.  Thus, d=-G(α)/α is bounded above by af*. Because -G(α)=c≥0, all this makes sense if, forallx, af(x)-(a2/2)g2(x)≥0 or -∂U:=2f/g2≥a (the opposite of the gradient of the potential function U in (2.12) is bounded below).
Let(ak;k≥1) be a nonincreasing sequence of [0,1]-valued real numbers. Let (αk;k≥1) be a sequence of nonnegative real numbers such that for all x∈(0,1)(2.41)α(x)=∑k≥1αke-akx<∞.
Whenever f is bounded above and, forallx, 2f/g2≥a1, we have
(2.42)-G(α)=∑k≥1αk(akf-ak22g2)e-akx≤f*∑k≥1αkake-akx<a1f*∑k≥1αke-akx=a1f*α.
Thus, d=-G(α)/α is bounded above by a1f*.
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 d=-G(α)/α.


### 2.6. Normalizing and Conditioning

Because the transformed processx¯t is nonconservative, it is of interest to inspect various conditionings in the sense of Yaglom, [7].(i) Consider again the process with infinitesimal generator G¯ losing mass due to killing and/or absorption at the boundaries. Integrating over y, with ρ¯t(x):=∫Ip¯(x;t,y)dy=P(τ¯x>t), we have (2.43)∂tρ¯t(x)=G¯(ρ¯t(x))=-d(x)ρ¯t(x)+G̃(ρ¯t(x)),
with ρ¯0(x)=1(x∈(0,1)). This gives the tail distribution of the full stopping time τ¯x.Defining the conditional probability densityq¯(x;t,y):=p¯(x;t,y)/ρ¯t(x), now with total mass 1, with q¯(x;0,y)=δy(x), we get (2.44)∂tq¯=-∂tρ¯t(x)ρ¯t(x)⋅q¯+G¯*(q¯)=(b¯t(x)-d(y))⋅q¯+G̃*(q¯).
The term b¯t(x)=-∂tρ¯t(x)/ρ¯t(x)>0 is the rate at which mass should be created to compensate the loss of mass of the process (x̃t;t≥0) due to its possible absorption at the boundaries and/or killing. Again, we have b¯t(x)→λ1, where λ1 is the smallest positive eigenvalue of -G, and therefore, putting ∂tq¯=0 in the latter evolution equation, we get that independently of the initial condition x(2.45)q¯(x;t,y)⟶t→∞q¯∞(y),
where q¯∞(y) is the solution to (2.46)-G̃*(q¯∞)=(λ1-d(y))⋅q¯∞,or-G¯*(q¯∞)=λ1⋅q¯∞.
With v1 the eigenvector of -G* associated to λ1, q¯∞(y) is of the product form (2.47)q¯∞(y)=α(y)v1(y)〈α,v1〉,
where 〈α,v1〉=∫01α(y)v1(y)dy. This results directly from the fact that G¯*(·¯)=α(y)G*(·/α(y)) and that v1 is the stated eigenvector of -G*. A different way to see this is as follows. We have (2.48)ρ¯t(x)=1α(x)∫01α(y)p(x;t,y)dy,
and the conditional density of x¯t given τ¯x>t is, therefore,(2.49)q¯(x;t,y)=α(y)p(x;t,y)∫01α(y)p(x;t,y)dy.
The rest follows from observing that, to the leading order in t in (2.15), for large time (2.50)p(x;t,y)~b1e-λ1t⋅u1(x)v1(y),
where u1 (v1,resp.) is the eigenvector of -G (-G*,resp.) associated to λ1 and b1=〈u1,v1〉-1. From this, it is clear that -(1/t)logρ¯t(x)→t→∞λ1 and (2.51)q¯(x;t,y)~e-λ1t⋅α(y)v1(y)e-λ1t⋅〈α,v1〉=q¯∞(y).
The limiting probability q¯∞=αv1/norm can, therefore, be interpreted as the Yaglom limit law of (x¯t;t≥0) conditioned on the event τ¯x>t.(ii) Under our assumptions, in the transformation of paths process, the transformed process (x¯t;t≥0) can both be absorbed at the boundaries and be killed. So, both τ̃x and τ̃x,∂ are finite with positive probability. We wish to understand the processes (x¯t;t≥0) conditioned on the events {τ̃x<τ̃x,∂} or {τ̃x,∂<τ̃x}, (see [16]).The probability mass cumulated at the boundaries{0,1} by time t clearly is [17] (2.52)P(τ̃x≤t)=P(τ̃x,∂>t)=12(∫0t[p¯(x;s,0)+p¯(x;s,1)]ds).Ast→∞, this probability tends to P(τ̃x,∂=∞)=P(τ̃x<τ̃x,∂)=:β(x). Note that (2.53)β(x)=P(x̃τ¯x∈{0,1})=E(1(x̃τ¯x∈{0,1})).Now (assumingx≠{0,1}), (2.54)G¯(P(τ̃x,∂>t))=12(∫0t[G¯(p¯(x;s,0))+G¯(p¯(x;s,1))]ds)=12(∫0t[∂s(p¯(x;s,0))+∂s(p¯(x;s,1))]ds)=12(p¯(x;t,0)+p¯(x;t,1))⟶t→∞0.Thus,β is defined by (2.55)G¯(β(x))=0,
[or G̃(β(x))=d(x)β(x)], with boundary conditions β(0)=β(1)=1. It serves as a positive harmonic function for G¯. This is a Sturm-Liouville problem to be solved for each case study.The density of the process(x¯t;t≥0) conditioned on the event {τ̃x<τ̃x,∂} is (2.56)p¯a(x;t,y)=β(y)β(x)p¯(x;t,y)=β(y)β(x)α(y)α(x)p(x;t,y).
The density of the process (x¯t;t≥0) conditioned on the event {τ̃x,∂<τ̃x} is (2.57)p¯∂(x;t,y)=1-β(y)1-β(x)p¯(x;t,y)=1-β(y)1-β(x)α(y)α(x)p(x;t,y).
Note that (2.58)∫01p¯a(x;t,y)dy=P(t̃x>t∣τ̃x<τ̃x,∂),∫01p¯∂(x;t,y)dy=P(τ̃x,∂>t∣τ̃x,∂<τ̃x),
and, with τ̅x=τ̃x∧τ̃x,∂,(2.59)P(τ¯x>t)=β(x)P(τ̃x>t∣τ̃x<τ̃x,∂)+(1-β(x))P(τ̃x,∂>t∣τ̃x,∂<τ̃x)=β(x)∫01β(y)β(x)p¯(x;t,y)dy+(1-β(x))∫011-β(y)1-β(x)p¯(x;t,y)dy=∫01p¯(x;t,y)dy,
as required, because this is the probability that x¯· is neither in {0,1} nor in state ∂ at time t.Note also thatq¯a(x;t,y)=p¯a(x;t,y)/∫01p¯a(x;t,y)dy[resp.,q¯∂(x;t,y)=p¯∂(x;t,y)/∫01p¯∂(x;t,y)dy] are the transition probability densities of (x¯t;t≥0) conditioned on the event {τ̃x<τ̃x,∂andτ̃x>t} [of (x¯t;t≥0) conditioned on the event {τ̃x,∂<τ̃xandτ̃x,∂>t,resp.}]. They are the Yaglom limits of both conditioned processes.The backward infinitesimal generators of both processes with transition probability densitiesp¯a and p¯∂ are, respectively, given by (2.60)G¯a(⋅)=1β(x)G¯(β(x)⋅),G¯∂(⋅)=11-β(x)G¯((1-β(x))⋅).
We get, respectively, (2.61)G¯a(⋅)=G̃(⋅)+β′g2β(x)∂x(⋅),G¯∂(⋅)=-d1-β⋅+G̃(⋅)-β′g21-β(x)∂x(⋅).
Thus, in G¯a(·), there is no multiplicative part (no killing) and a shift in the drift, showing that the associated conditioned process (x̃a,t;t≥0) obeys the SDE (2.62)dx̃a=f̃a(x̃a)dt+g(x̃a)dwt,
with drift (2.63)f̃a(x)=f̃(x)+β′g2β(x)=f(x)+α′g2α(x)+β′g2β(x)=f(x)+g2(x)[α′α(x)+β′β(x)].
This process is ultimately absorbed at {0,1}.InG¯∂(·), there is a killing multiplicative part which is enhanced d/(1-β)>d and a shift in the drift, showing that the associated conditioned process (x̃∂,t;t≥0) exhibits a faster killing rate, but the drift shift guarantees that (x̃∂,t;t≥0) is not absorbed at the boundaries. We have (2.64)f̃∂(x)=f(x)+g2(x)[α′α(x)-β′1-β(x)].Additive Functionals of the Transformed Process.
for the new process(x̃t;t≥0), it is also of interest to evaluate additive functionals along their own sample paths. Let then α̃(x):=E(∫0τ¯xc̃(x̃s)ds+d̃(x̃τ¯x)) be such an additive functional where the functions c̃ and d̃ are themselves both nonnegative. It solves
(2.65)-G¯(α̃)=c̃ifx∈I,α̃=d̃ifx∈∂I.
Then, recalling the expression of the Green function 𝔤(x,y) of (xt;t≥0) in (2.22), we find explicitly
(2.66)α̃(x)=1α(x)∫Ig(x,y)α(y)c̃(y)dy.Specific transformations of interest.
(i) The case c=0 deserves a special treatment. Indeed, in this case, τ̃x,∂=∞ and so τ¯x:=τ̃x, the absorption time for the process (x̃t;t≥0) governed by the new SDE (2.38). Here, G¯=G̃. Assuming α solves -G(α)=0 if x∈I with boundary conditions α(0)=0 and α(1)=1 (α(0)=1 and α(1)=0,resp.), the new process (x̃t;t≥0) is just (xt;t≥0) conditioned on exiting at x=1 (at x=0,resp.). In the first case, the boundary 1 is exit, whereas 0 is entrance; α reads
(2.67)α(x)=∫0xe-2∫0y(f(z)/(g2(z)))dzdy∫01e-2∫0y(f(z)/(g2(z)))dzdy,
with
(2.68)f̃(x)=f(x)+g2(x)e-2∫0x(f(z)/(g2(z)))dz∫0xe-2∫0y(f(z)/(g2(z)))dzdy,
giving the new drift. In the second case, α(x)=(∫x1e-2∫0y(f(z)/g2(z))dzdy/(∫01e-2∫0y(f(z)/g2(z))dzdy)), and the boundary 0 is exit, whereas 1 is entrance. Thus, τ̃x is just the exit time at x=1 (at x=0,resp.). Let α̃(x):=E(τ̃x). Then, α̃(x) solves -G̃(α̃)=1, whose explicit solution is
(2.69)α̃(x)=1α(x)∫Ig(x,y)α(y)dy,
in terms of 𝔤(x,y), the Green function of (xt;t≥0).
Example 2.1.
Consider the WF model on[0,1] with selection for which, with σ∈R,f(x)=σx(1-x) and g2(x)=x(1-x). Assume that  α solves -G(α)=0 if x∈(0,1) with α(0)=0 and α(1)=1; one gets, α(x)=(1-e-2σx)/(1-e-2σ). The diffusion corresponding to (2.38) has the new drift: f̃(x)=σx(1-x)coth(σx), independently of the sign of σ. It models the WF diffusion with selection conditioned on exit at ∘=1.
(ii) Assume that α now solves -G(α)=1 if x∈I with boundary conditions α(0)=α(1)=0. In this case study, one selects sample paths of (xt;t≥0) with a large mean absorption time α(x)=E(τx). Sample paths with large sojourn time in I are favored. We have
(2.70)α(x)=∫Ig(x,y)dy,
where 𝔤(x,y) is the Green function (2.22). The boundaries of (x̃t;t≥0) are now both entrance boundaries and so τ̃x=∞. (x̃t;t≥0) is not absorbed at the boundaries. The stopping time τ¯x of (x̃t;t≥0) is just its killing time τ̃x,∂. Let α̃(x):=E(τ̃x,∂). Then, α̃(x) solves -G¯(α̃)=1, α̃(0)=α̃(1)=0, with explicit solution
(2.71)α̃(x)=1α(x)∫Ig(x,y)α(y)dy.
(iii) Assume that α now solves -G(α)=δy(x) if x∈I with boundary conditions α(0)=α(1)=0. In this case study, one selects sample paths of (xt;t≥0) with a large sojourn time density at y recalling α(x)=:𝔤(x,y)=E(∫0τxδy(xs)ds). The stopping time τ̃y(x) of (x̃t;t≥0) occurs at rate δy(x)/𝔤(x,y). It is a killing time when the process is at y for the last time after a geometrically distributed number of passages there with rate 1/𝔤(x,y) (or with success probability 1/(1+𝔤(x,y))). Let α̃y(x):=E(τ̃y(x)). Then, α̃y(x) solves -G¯(α̃)=1, with explicit solution
(2.72)α̃y(x)=1g(x,y)∫Ig(x,z)g(z,y)dz.
when x=1/N,α̃y(1/N) may be viewed as the age of a mutant currently observed to present frequency y, see [18].
The Green function aty0∈(0,1) of the transformed process (x̃t;t≥0) is 𝔤̃y(x,y0) solution to -G¯(𝔤̃y)=δy0(x). It takes the simple form
(2.73)g̃y(x,y0)=1g(x,y)∫Ig(x,z)g(z,y)δy0(z)dz=g(y0,y)g(x,y)g(x,y0).
(iv) Let λ1 be the smallest non-null eigenvalue of the infinitesimal generator G. Let α=u1 be the corresponding eigenvector, that is, satisfying, -Gu1=λ1u1 with boundary conditions u1(0)=u1(1)=0. Then, c=λu1. The new KB operator associated to the transformed process (x¯t;t≥0) is
(2.74)G¯(⋅)=1αG(α)⋅+G̃(⋅)=-λ1⋅+G̃(⋅),
obtained while killing the sample paths of the process (x̃t;t≥0) governed by G̃ at constant death rate d=λ1. The transition probability of the transformed stochastic process (x¯t;t≥0) is
(2.75)p¯(x;t,y)=u1(y)u1(x)p(x;t,y).
Define p̃(x;t,y)=eλ1tp¯(x;t,y). It is the transition probability of the process (x̃t;t≥0) governed by G̃; it corresponds to the original process (xt;t≥0) conditioned on never hitting the boundaries {0,1} (the so-called Q-process of (xt;t≥0), see [19]). It is simply obtained from (xt;t≥0) by adding the additional drift term (u1′/u1)g2 to f, where u1 is the eigenvector of G associated to its smallest non-null eigenvalue. The determination of α=u1 is a Sturm-Liouville problem. When t is large, to the dominant order
(2.76)p(x;t,y)~e-λ1tu1(x)v1(y)〈u1,v1〉,
where v1 is the Yaglom limit law of (xt;t≥0). Therefore
(2.77)p̃(x;t,y)~eλ1tu1(y)u1(x)e-λ1tu1(x)v1(y)〈u1,v1〉=u1(y)v1(y)〈u1,v1〉.
Thus, the limit law of the Q-process (x̃t;t≥0) is the normalized Hadamard product of the eigenvectors u1 and v1 associated, respectively, to G and G*. On the other hand, the limit law of (x̃t;t≥0) is directly given by
(2.78)p̃(x;t,y)⟶t→∞p̃(y)=1Zg2(y)e2∫0y((f(z)+((u1′/u1)g2)(z))/(g2(z)))dz=u12(y)Zg2(y)e2∫0y(f(z)/(g2(z)))dz,,
where Z is the appropriate normalizing constant. Comparing (2.77) and (2.78)
(2.79)v1(y)=u1(y)g2(y)e2∫0y(f(z)/(g2(z)))dz=u1(y)m(y).
The eigenvector v1 associated to G* is, therefore, equal to the eigenvector u1 associated to G times the speed density of (xt;t≥0).
When dealing for example with the neutral Wright-Fisher diffusion, it is known thatλ1=1 with u1=x(1-x) and v1≡1 (see Section 4.3, example (ii)). The Q-process (x̃t;t≥0) in this case obeys
(2.80)dx̃t=(1-2x̃t)dt̅+x̃t(1-x̃t)dwt,
with the stabilizing drift toward 1/2: f̃(x)=(u1′/u1)g2(x)=1-2x.The limit law of the Q-process (x̃t;t≥0) in this case is 6y(1-y). 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 generatorG¯(·)=-d(x)·+G̃(·) and applying the reciprocal Doob transform defined by (2.81)α̃(x)=1α(x)
leads to α̃(x)-1G¯(α̃(x)·)=G(·). Indeed, (2.82)α̃(x)-1G¯(α̃(x)⋅)=-d(x)⋅+α̃(x)-1G̃(α̃(x)⋅),α̃(x)-1G̃(α̃(x)⋅)=α̃(x)-1G̃(α̃(x))⋅+G̃(⋅)+1α̃g2α̃′∂x=d(x)⋅+G(⋅),
because (2.83)α̃(x)-1G̃(α̃(x))=+d(x),1α̃g2α̃′∂x=αg2(1α)′∂x=-α′αg2∂x.
Note that G¯(α̃(x))=G̃(α̃(x))-d(x)α̃(x)=0.This suggests that starting from a diffusion process with infinitesimal generatorG̃ (without its killing part) and applying the reciprocal Doob transform α̃(x)=1/(α(x)), one ends up with a modified process whose infinitesimal generator is (2.84)G̃¯(⋅)=G(⋅)+b(x),
where G(·)=G̃(·)-(α′/α)g2(x)∂x and (2.85)b(x)=α̃(x)-1G̃(α̃(x))=d(x)>0
is now a pure birth rate. Note that α̃ is now a subharmonic function for G̃ because -G̃(α̃)=-bα̃≤0.Letβ̃(x):=α̃(x)-1=α(x). Because G(α(x))=-α(x)b(x), we have (2.86)G̃¯(β̃(x))=G(β̃(x))+b(x)β̃(x)=0,
and so β̃(x)>0 is harmonic for G̃¯: Doob-transforming G̃¯ using β̃, we get (2.87)β̃-1G̃¯(β̃⋅)=b(x)⋅+β̃-1G(β̃⋅)=G̃(⋅),
which is the infinitesimal generator of the original diffusion process.Under our assumptions, both processxt and x̃t with infinitesimal generators G and G̃ are nonconservative diffusion processes with absorbing barriers. Further, b(x)>0 is bounded above. Therefore, b(x) may be written as (2.88)b(x)=b*(μ(x)-1),
where b*=maxxb(x)>0 and μ(x)∈[1,2].The process with infinitesimal generatorG̃¯ is now a pure binary branching diffusion process. For this class of models, an initial particle started at x obeys a diffusion process with infinitesimal generator G, absorbed when it hits the boundaries. At some random (mean b*) exponential time, this particle dies, giving birth in the process to a random number M(x) (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 EM(x)=μ(x).The process with infinitesimal generatorG̃¯ is, thus, a branching diffusion with supercritical binary splitting mechanism (μ(x)>1). 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.LetNt(x) be the global number of particles which are alive in the system at each time t, descending from an Eve particle started at x, and let (2.89)T(x)=inf⁡(t>0:Nt(x)=0)
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 (2.90)P(T(x)<∞,Nt(x)=0∀t≥T(x))=1,
uniformly in x. 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.


## 2.1. Generalities on One-Dimensional Diffusions on the Interval[0,1]

Let(wt;t≥0) be a standard one-dimensional Brownian (Wiener) motion. Consider a one-dimensional Itô diffusion driven by (wt;t≥0) on the interval say [0,1]; see [14]. We will let I=(0,1). Assume that it has locally Lipschitz continuous drift (x) and local standard deviation (volatility) g(x), namely, consider the stochastic differential equation (SDE)(2.1)dxt=f(xt)dt+g(xt)dwt,x0=x∈(0,1).
The condition on f(x) and g(x) guarantees in particular that there is no point x* in I for which |f(x)| or |g(x)| would blow up and diverge as |x-x*|→0.The Kolmogorov backward infinitesimal generator of (2.1) is G=f(x)∂x+(1/2)g2(x)∂x2. As a result, for all suitable ψ in the domain of the operator St:=etG, u:=u(x,t)=Eψ(xt∧τx) satisfies the Kolmogorov backward equation (KBE) (2.2)∂tu=G(u);u(x,0)=ψ(x).
In the definition of the mathematical expectation u, we have t∧τx:=inf(t,τx), where τx indicates a random time at which the process should possibly be stopped (absorbed), given the process was started in x. The description of this (adapted) absorption time is governed by the type of boundaries which {0,1} are to (xt;t≥0).


## 2.2. Natural Coordinate, Scale, and Speed Measure

For such Markovian diffusions, it is interesting to consider theG-harmonic coordinate φ∈C2 belonging to the kernel of G that is, satisfying G(φ)=0. For φ and its derivative φ′:=dφ/dy, with (x0,y0)∈(0,1), one finds (2.3)φ′(y)=φ′(y0)e-2∫y0y(f(z)/(g2(z)))dz,φ(x)=φ(x0)+φ′(y0)∫x0xe-2∫y0y(f(z)/(g2(z)))dzdy.
One should choose a version of φ satisfying φ′(y)>0, y∈I. The function φ kills the drift f of (xt;t≥0) in the sense that considering the change of variable yt=φ(xt),(2.4)dyt=(φ′g)(φ-1(yt))dwt,y0=φ(x).
The driftless diffusion (yt;t≥0) is often termed the diffusion in natural coordinates with state-space [φ(0),φ(1)]. Its volatility is g̃(y):=(φ′g)(φ-1(y)). The function φ is often called the scale function.Wheneverφ(0)>-∞ and φ(1)<+∞, one can choose the integration constants defining φ(x) so that (2.5)φ(x)=∫0xe-2∫0y(f(z)/(g2(z)))dzdy∫01e-2∫0y(f(z)/(g2(z)))dzdy,
with φ(0)=0 and φ(1)=1. In this case, the state-space of (yt;t≥0) is again [0,1], the same as for (xt;t≥0).Finally, considering the random time changet→θt with inverse: θ→tθ defined by θtθ=θ and (2.6)θ=∫0tθg̃2(ys)ds,
the novel diffusion (wθ:=ytθ;θ≥0) is easily checked to be identical in law to a standard Brownian motion. Let now δy(·)= weak-limε↓0(1/2ε)1(·∈(y-ε,y+ε)) stand for the Dirac delta mass at y. The random time θt can be expressed as (2.7)θt=∫01dx⋅m(x)∫0tδφ(x)(ws)ds=∫0tm(φ-1(ws))ds,
where m(x):=1/(g2φ′)(x) is the (positive) speed density at x=φ-1(y) and Lt(y):=limε↓0(1/2ε)∫0t1(ws∈(y-ε,y+ε))ds the local time at y of the Brownian motion before time t. Both the scale function φ and the speed measure dμ=m(x)·dx are, therefore, essential ingredients to reduce the original stochastic process (xt;t≥0) to the standard Brownian motion (wt;t≥0). Indeed, it follows from the above arguments that if θt=∫0tm(xs)ds, then (φ(xθt);t≥0) is a Brownian motion. The Kolmogorov backward infinitesimal generator G may be written in Feller form (2.8)G(⋅)=12ddμ(ddφ⋅).Examples 2.2 (from population genetics).
(i) Assume thatf(x)=0 and g2(x)=x(1-x). This is the neutral Wright-Fisher (WF) model discussed at length later. This diffusion is already in natural scale and φ(x)=x, m(x)=[x(1-x)]-1. The speed measure is not integrable.
(ii) Withu1,u2>0, assume f(x)=u1-(u1+u2)x and g2(x)=x(1-x). This is the Wright-Fisher model with mutation. The parameters u1,u2 can be interpreted as mutation rates. The drift vanishes when x=u1/(u1+u2) which is an attracting point for the dynamics. Here,φ′(y)=φ′(y0)y-2u1(1-y)-2u2,φ(x)=φ(x0)+φ′(y0)∫x0xy-2u1(1-y)-2u2dy, with φ(0)=-∞ and φ(1)=+∞ if u1,u2>1/2. The speed measure density is m(x)∝x2u1-1(1-x)2u2-1 and so is always integrable.
(iii) Withσ∈R, assume a model with quadratic logistic drift f(x)=σx(1-x) and local variance g2(x)=x(1-x). This is the WF model with selection. For this diffusion (see [15]), φ(x)=((1-e-2σx)/(1-e-2σ)) and m(x)∝[x(1-x)]-1e2σx 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 thatf(x) and g(x) are now differentiable in I. Let then p(x;t,y) stand for the transition probability density function of xt at y given x0=x. Then, p:=p(x;t,y) is the smallest solution to the Kolmogorov forward (Fokker-Planck) equation (KFE) (2.9)∂tp=G*(p),p(x;0,y)=δy(x),
where G*(·)=-∂y(f(y)·)+(1/2)∂y2(g2(y)·) is the adjoint of G (G* acts on the terminal value y, whereas G acts on the initial value x). The way one can view this PDE depends on the type of boundaries that {0,1} are.We will next suppose that the boundaries∘:=0 or 1 are both exit (or absorbing) boundaries. From the Feller classification of boundaries, this will be the case if forally0∈(0,1)(2.10)(i)m(y)∉L1(y0,∘),(ii)φ′(y)∫y0ym(z)dz∈L1(y0,∘),
where a function f(y)∈L1(y0,∘) if ∫y0∘|f(y)|dy<+∞.In this case, a sample path of(xt;t≥0) can reach ∘ from the inside of I in finite time but cannot reenter. The sample paths are absorbed at ∘. There is an absorption at ∘ at time τx,∘=inf(t>0:xt=∘|x0=x) and P(τx,∘<∞)=1. Whenever both boundaries {0,1} are absorbing, the diffusion xt should be stopped at τx:=τx,0∧τx,1. Would none of the boundaries {0,1} be absorbing, then τx=+∞, 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,p(x;t,y) is a subprobability. Letting ρt(x):=∫01p(x;t,y)dy, we clearly have ρt(x)=P(τx>t). Such models are nonconservative.For one-dimensional diffusions, the transition densityp(x;t,y) is reversible with respect to the speed density ([8, Chapter 15, Section  13]) and so detailed balance holds (2.11)m(x)p(x;t,y)=m(y)p(y;t,x),0<x,y<1.
The speed density m(y) satisfies G*(m)=0. It may be written as a Gibbs measure with density: m(y)∝(1/g2(y))e-U(y), where the potential function U(y) reads (2.12)U(y):=-2∫0yf(z)g2(z)dz,0<y<1,
and with the measure dy/g2(y) standing for the reference measure.Further, ifp(s,x;t,y) is the transition probability density from (s,x) to (t,y), s<t, then -∂sp=G(p), with terminal condition p(t,x;t,y)=δy(x) and so p(s,x;t,y) 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 p:=p(x;t,y), operating the time substitution t-s→t, p itself solves the KBE (2.13)∂tp=G(p),p(x;0,y)=δy(x).
In particular, integrating over y, ∂tρt(x)=G(ρt(x)), with ρ0(x)=1(x∈(0,1)).p(x;t,y) being a sub-probability, we may define the normalized conditional probability density q(x;t,y):=p(x;t,y)/ρt(x), now with total mass 1. We get (2.14)∂tq=-∂tρt(x)ρt(x)⋅q+G*(q),q(x;0,y)=δy(x).
The term bt(x):=-∂tρt(x)/ρt(x)>0 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 (xt;t≥0) at the boundaries. In this creation of mass process, a diffusing particle started in x dies at rate bt(x) at point (t,y), where it is duplicated in two new independent particles both started at y (resulting in a global birth) evolving in the same diffusive way  (consider a diffusion process with forward infinitesimal generator G* governing the evolution of p(x;t,y). 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 d and b depend on the current location y of the path. Then, the process with the birth and death opportunities of a path has the infinitesimal generator λ(y)·+G*(·), where λ(y)=b(y)-d(y). The rate can also depend on t and x).  The birth rate function bt(x) depends here on x and t, not on y.When the boundaries ofxt are absorbing, the spectra of both -G and -G* are discrete (see [8, Page 330]): There exist positive eigenvalues (λk)k≥1 ordered in ascending sizes and eigenvectors (vk,uk)k≥1 of both -G* and -G satisfying -G*(vk)=λkvk and -G(yk)=λkuk such that  with 〈uk,vk〉:=∫01uk(x)vk(x)dx and bk:=〈uk,vk〉-1, the spectral representation (2.15)p(x;t,y)=∑k≥1bke-λktuk(x)vk(y)
holds.Letλ1>λ0=0 be the smallest nonnull eigenvalue of the infinitesimal generator -G* (and of -G). Clearly, -(1/t)logρt(x)→t→∞λ1 and by L' Hospital rule, therefore, bt(x)→t→∞λ1. Putting ∂tq=0 in the latter evolution equation, independently of the initial condition x(2.16)q(x;t,y)⟶t→∞q∞(y)=v1(y),
where v1 is the eigenvector of -G* associated to λ1, satisfying -G*v1=λ1v1. The limiting probability v1/norm (after a proper normalization) is called the Yaglom limit law of (xt;t≥0) conditioned on being currently alive at all time t (see [7]).


## 2.4. Additive Functionals Along Sample Paths

Let(xt;t≥0) be the diffusion model defined by (2.1) on the interval I, where both endpoints are assumed absorbing (exit). This process is, thus, transient and nonconservative. We wish to evaluate the nonnegative additive quantities (2.17)α(x)=E(∫0τxc(xs)ds+d(xτx)),
where the functions c and d are both assumed nonnegative on I and ∂I={0,1}. The functional α(x)≥0 solves the Dirichlet problem (2.18)-G(α)=cifx∈I,α=difx∈∂I,
and α is a superharmonic function for G, satisfying -G(α)≥0.Some Examples.
(1) Assume thatc=1 and d=0: here, α=E(τx) is the mean time of absorption (average time spent in (0,1) before absorption), solution to
(2.19)-G(α)=1ifx∈I,α=0ifx∈∂I.
(2) Whenever both{0,1} are exit boundaries, it is of interest to evaluate the probability that xt first hits [0,1] (say) at 1, given x0=x. This can be obtained by choosing c=0 and d(∘)=1(∘=1).
Let thenα=:α1(x)=P(xtfirsthits[0,1]at1∣x0=x).α1(x) is a G-harmonic function solution to G(α1)=0, with boundary conditions α1(0)=0 and α1(1)=1. Solving this problem, we get(2.20)α1(x)=φ(x)-φ(0)φ(1)-φ(0)=∫0xdye-2∫0y(f(z)/(g2(z)))dz∫01dye-2∫0y(f(z)/(g2(z)))dz.
On the contrary, choosing α0(x) to be a G-harmonic function with boundary conditions α0(0)=1 and α0(1)=0, α0(x)=P(xtfirsthits[0,1]at0∣x0=x)=1-α1(x).
(3) Lety∈I and put c=(1/2ε)1(x∈(y-ε,y+ε)) and d=0. As ε→0, c converges weakly to δy(x) and, α=:𝔤(x,y)=E(lim(1/2ε)∫0τx1(xs∈(y-ε,y+ε))ds)=∫0∞p(x;s,y)ds is the Green function, solution to(2.21)-G(g)=δy(x)ifx∈I,g=0ifx∈∂I.𝔤 is, therefore, the mathematical expectation of the local time at y, starting from x (the sojourn time density at y). The solution is known to be (see [8, page 198] or [5, page 280])(2.22)g(x,y)=2(φ(x)-φ(0))(φ(1)-φ(y))(g2φ′)(y)(φ(1)-φ(0))ifx≤y,g(x,y)=2(φ(1)-φ(x))(φ(y)-φ(0))(g2φ′)(y)(φ(1)-φ(0))ifx≥y.
The Green function is of particular interest to solve the general problem of evaluating additive functionals α(x). Indeed, as is well known, see [8], for example, the integral operator with respect to the Green kernel inverts the second-order operator -G leading to(2.23)α(x)=∫Ig(x,y)c(y)dyifx∈I,α=difx∈∂I.
Under this form, α(x) appears as a potential function and all potential function is superharmonic. Note that for all harmonic function h≥0 satisfying -G(h)=0,
(2.24)αh(x):=∫Ig(x,y)c(y)dy+h(x)
is again superharmonic because -G(αh)=c≥0.
(4) Also of interest are the additive functionals of the type(2.25)αλ(x)=E(∫0τxe-λsc(xs)ds+d(xτx)),
where the functions c and d are again both assumed to be nonnegative. The functional αλ(x)≥0 solves the Dynkin problem, [8] (2.26)(λI-G)(αλ)=cifx∈I,αλ=difx∈∂I
involving the action of the resolvent operator (λI-G)-1 on c.
Wheneverc(x)=δy(x),d=0, then(2.27)αλ=:gλ(x,y)=E(∫0τxe-λsδy(xs)ds)=∫0∞e-λsp(x;s,y)ds
is the λ-potential function, solution to
(2.28)(λI-G)(gλ)=δy(x)ifx∈I,gλ=0ifx∈∂I.𝔤λ is, therefore, the mathematical expectation of the exponentially damped local time at y, starting from x (the temporal Laplace transform of the transition probability density from x to y at t), with g0=g. Then, it holds that
(2.29)αλ(x)=∫Igλ(x,y)c(y)dyifx∈I,αλ=difx∈∂I.
The λ-potential function is also useful in the computation of the distribution of the first-passage time τx,y to y starting from x. From the convolution formula,
(2.30)p(x;t,y)=∫0tP(τx,y∈ds)p(y;t-s,y),
and taking the Laplace transform of both sides with respect to time, we obtain the Laplace-Stieltjes transform (LST) of the law of τx,y as
(2.31)E(e-λτx,y)=gλ(x,y)gλ(y,y).
We have P(τx,y<∞)=𝔤(x,y)/(𝔤(y,y))∈(0,1) as a result of both terms in the ratio being finite and x,y belonging to the same transience class of the process (under our assumptions that the boundaries are absorbing). Note that from the reversibility property
(2.32)m(x)g(x,y)=m(y)g(y,x).


## 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(xt;t≥0) as in (2.1) with absorbing barriers. Let p(x;t,y) be its transition probability, and let τx be its absorption time at the boundaries.Letα(x):=E(∫0τxc(xs)ds+d(xτx)) be a nonnegative additive functional solving (2.33)-G(α)=cifx∈I,α=difx∈∂I.
Recall the functions c and d are both chosen nonnegative so that so is α.Define a new transformed stochastic process(x¯t;t≥0) by its transition probability (2.34)p¯(x;t,y)=α(y)α(x)p(x;t,y).
In this construction of (x¯t;t≥0) through a change of measure, sample paths of (xt;t≥0) for which α(y) is large are favored. This is a selection of paths procedure due to Doob (see [11]).Now, the KFE forp¯ clearly is ∂tp¯=G¯*(p¯), with p(x;0,y)=δy(x) and G¯*(p¯)=α(y)G*(p¯/α(y)). The Kolmogorov backward operator of the transformed process is, therefore, by duality (2.35)G¯(⋅)=1α(x)G(α(x)⋅).
Developing, with α′(x):=dα(x)/dx and G̃(·):=(α′/α)g2∂x(·)+G(·), we get (2.36)G¯(⋅)=1αG(α)⋅+G̃(⋅)=-cα⋅+G̃(⋅),
and the new KB operator can be obtained from the latter by adding a drift term (α′/α)g2∂x to the one in G of the original process to form a new process (x̃t;t≥0) with the KB operator G̃ and by killing its sample paths at death rate d(x):=(c/α)(x) (provided c≠0). Note that (2.37)d(x)=1(1/(c(x)))∫01g(x,y)c(y)dy.
In others words, with f̃(x):=f(x)+(α′/α)g2(x), the novel time-homogeneous SDE to consider is (2.38)dx̃t=f̃(x̃t)dt+g(x̃t)dwt,x̃0=x∈(0,1),
possibly killed at rate d=c/α as soon as c≠0. Whenever (x̃t;t≥0) is killed, it enters conventionally into some coffin state {∂} added to the state-space. Let τ̃x be the new absorption time at the boundaries of (x̃t;t≥0) started at x (with τ̃x=∞ would the boundaries be inaccessible to the new process x̃t which we ruled out). Let τ̃x,∂ be the killing time of (x̃t;t≥0) started at x (the hitting time of ∂), with τ̃x,∂=∞ if c=0. Then, τ¯x:=τ̃x∧τ̃x,∂ is the novel stopping time of (x̃t;t≥0). The SDE for (x̃t;t≥0), together with its global stopping time τ¯x characterize the new process (x¯t;t≥0) with generator G¯ to consider.In the sequel, we shall limit ourselves to the cases for which the following additional conditions hold on the transformed process.(i) Nonconservativeness of x̃t.
We will next suppose that the boundaries∘:=0 or 1 are both exit (or absorbing) boundaries for the new process x̃t in (2.38). From the Feller criterion for exit boundaries, this will be the case if forally0∈(0,1)(2.39)m̃(y)∉L1(y0,∘),φ̃′(y)∫y0ym̃(z)dz∈L1(y0,∘),
where m̃(y)=1/((g2φ̃′)(y)) is the new speed measure density for x̃t and φ̃ its scale function. Recalling f̃=f+(α′/α)g2 and g̃2=g2, we have
(2.40)φ̃′(y)=φ̃′(y0)e-2∫y0y(f̃(z)/(g2(z)))dz=φ̃′(y0)α2(y)φ′(y),φ̃(x)=φ̃(x0)+φ̃′(y0)∫x0xα-2(y)e-2∫y0y(f(z)/(g2(z)))dzdy.
So, we assume here that x̃t obeys itself a nonconservative diffusion.(ii) Boundedness of the Killing Rate d.
In some examples, the killing rated=-G(α)/α is bounded above. For example, suppose that the drift of the diffusion process (xt;t≥0) is bounded above by f*=maxx(f(x))>0. (If the drift of (xt;t≥0) is bounded below by f*<0, we are led to the same conclusions while considering the process 1-xt instead of xt.) Then, choosing α(x)=e-ax,a>0,-G(α)=(af-(a2/2)g2)α<af*α.  Thus, d=-G(α)/α is bounded above by af*. Because -G(α)=c≥0, all this makes sense if, forallx, af(x)-(a2/2)g2(x)≥0 or -∂U:=2f/g2≥a (the opposite of the gradient of the potential function U in (2.12) is bounded below).
Let(ak;k≥1) be a nonincreasing sequence of [0,1]-valued real numbers. Let (αk;k≥1) be a sequence of nonnegative real numbers such that for all x∈(0,1)(2.41)α(x)=∑k≥1αke-akx<∞.
Whenever f is bounded above and, forallx, 2f/g2≥a1, we have
(2.42)-G(α)=∑k≥1αk(akf-ak22g2)e-akx≤f*∑k≥1αkake-akx<a1f*∑k≥1αke-akx=a1f*α.
Thus, d=-G(α)/α is bounded above by a1f*.
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 d=-G(α)/α.


## 2.6. Normalizing and Conditioning

Because the transformed processx¯t is nonconservative, it is of interest to inspect various conditionings in the sense of Yaglom, [7].(i) Consider again the process with infinitesimal generator G¯ losing mass due to killing and/or absorption at the boundaries. Integrating over y, with ρ¯t(x):=∫Ip¯(x;t,y)dy=P(τ¯x>t), we have (2.43)∂tρ¯t(x)=G¯(ρ¯t(x))=-d(x)ρ¯t(x)+G̃(ρ¯t(x)),
with ρ¯0(x)=1(x∈(0,1)). This gives the tail distribution of the full stopping time τ¯x.Defining the conditional probability densityq¯(x;t,y):=p¯(x;t,y)/ρ¯t(x), now with total mass 1, with q¯(x;0,y)=δy(x), we get (2.44)∂tq¯=-∂tρ¯t(x)ρ¯t(x)⋅q¯+G¯*(q¯)=(b¯t(x)-d(y))⋅q¯+G̃*(q¯).
The term b¯t(x)=-∂tρ¯t(x)/ρ¯t(x)>0 is the rate at which mass should be created to compensate the loss of mass of the process (x̃t;t≥0) due to its possible absorption at the boundaries and/or killing. Again, we have b¯t(x)→λ1, where λ1 is the smallest positive eigenvalue of -G, and therefore, putting ∂tq¯=0 in the latter evolution equation, we get that independently of the initial condition x(2.45)q¯(x;t,y)⟶t→∞q¯∞(y),
where q¯∞(y) is the solution to (2.46)-G̃*(q¯∞)=(λ1-d(y))⋅q¯∞,or-G¯*(q¯∞)=λ1⋅q¯∞.
With v1 the eigenvector of -G* associated to λ1, q¯∞(y) is of the product form (2.47)q¯∞(y)=α(y)v1(y)〈α,v1〉,
where 〈α,v1〉=∫01α(y)v1(y)dy. This results directly from the fact that G¯*(·¯)=α(y)G*(·/α(y)) and that v1 is the stated eigenvector of -G*. A different way to see this is as follows. We have (2.48)ρ¯t(x)=1α(x)∫01α(y)p(x;t,y)dy,
and the conditional density of x¯t given τ¯x>t is, therefore,(2.49)q¯(x;t,y)=α(y)p(x;t,y)∫01α(y)p(x;t,y)dy.
The rest follows from observing that, to the leading order in t in (2.15), for large time (2.50)p(x;t,y)~b1e-λ1t⋅u1(x)v1(y),
where u1 (v1,resp.) is the eigenvector of -G (-G*,resp.) associated to λ1 and b1=〈u1,v1〉-1. From this, it is clear that -(1/t)logρ¯t(x)→t→∞λ1 and (2.51)q¯(x;t,y)~e-λ1t⋅α(y)v1(y)e-λ1t⋅〈α,v1〉=q¯∞(y).
The limiting probability q¯∞=αv1/norm can, therefore, be interpreted as the Yaglom limit law of (x¯t;t≥0) conditioned on the event τ¯x>t.(ii) Under our assumptions, in the transformation of paths process, the transformed process (x¯t;t≥0) can both be absorbed at the boundaries and be killed. So, both τ̃x and τ̃x,∂ are finite with positive probability. We wish to understand the processes (x¯t;t≥0) conditioned on the events {τ̃x<τ̃x,∂} or {τ̃x,∂<τ̃x}, (see [16]).The probability mass cumulated at the boundaries{0,1} by time t clearly is [17] (2.52)P(τ̃x≤t)=P(τ̃x,∂>t)=12(∫0t[p¯(x;s,0)+p¯(x;s,1)]ds).Ast→∞, this probability tends to P(τ̃x,∂=∞)=P(τ̃x<τ̃x,∂)=:β(x). Note that (2.53)β(x)=P(x̃τ¯x∈{0,1})=E(1(x̃τ¯x∈{0,1})).Now (assumingx≠{0,1}), (2.54)G¯(P(τ̃x,∂>t))=12(∫0t[G¯(p¯(x;s,0))+G¯(p¯(x;s,1))]ds)=12(∫0t[∂s(p¯(x;s,0))+∂s(p¯(x;s,1))]ds)=12(p¯(x;t,0)+p¯(x;t,1))⟶t→∞0.Thus,β is defined by (2.55)G¯(β(x))=0,
[or G̃(β(x))=d(x)β(x)], with boundary conditions β(0)=β(1)=1. It serves as a positive harmonic function for G¯. This is a Sturm-Liouville problem to be solved for each case study.The density of the process(x¯t;t≥0) conditioned on the event {τ̃x<τ̃x,∂} is (2.56)p¯a(x;t,y)=β(y)β(x)p¯(x;t,y)=β(y)β(x)α(y)α(x)p(x;t,y).
The density of the process (x¯t;t≥0) conditioned on the event {τ̃x,∂<τ̃x} is (2.57)p¯∂(x;t,y)=1-β(y)1-β(x)p¯(x;t,y)=1-β(y)1-β(x)α(y)α(x)p(x;t,y).
Note that (2.58)∫01p¯a(x;t,y)dy=P(t̃x>t∣τ̃x<τ̃x,∂),∫01p¯∂(x;t,y)dy=P(τ̃x,∂>t∣τ̃x,∂<τ̃x),
and, with τ̅x=τ̃x∧τ̃x,∂,(2.59)P(τ¯x>t)=β(x)P(τ̃x>t∣τ̃x<τ̃x,∂)+(1-β(x))P(τ̃x,∂>t∣τ̃x,∂<τ̃x)=β(x)∫01β(y)β(x)p¯(x;t,y)dy+(1-β(x))∫011-β(y)1-β(x)p¯(x;t,y)dy=∫01p¯(x;t,y)dy,
as required, because this is the probability that x¯· is neither in {0,1} nor in state ∂ at time t.Note also thatq¯a(x;t,y)=p¯a(x;t,y)/∫01p¯a(x;t,y)dy[resp.,q¯∂(x;t,y)=p¯∂(x;t,y)/∫01p¯∂(x;t,y)dy] are the transition probability densities of (x¯t;t≥0) conditioned on the event {τ̃x<τ̃x,∂andτ̃x>t} [of (x¯t;t≥0) conditioned on the event {τ̃x,∂<τ̃xandτ̃x,∂>t,resp.}]. They are the Yaglom limits of both conditioned processes.The backward infinitesimal generators of both processes with transition probability densitiesp¯a and p¯∂ are, respectively, given by (2.60)G¯a(⋅)=1β(x)G¯(β(x)⋅),G¯∂(⋅)=11-β(x)G¯((1-β(x))⋅).
We get, respectively, (2.61)G¯a(⋅)=G̃(⋅)+β′g2β(x)∂x(⋅),G¯∂(⋅)=-d1-β⋅+G̃(⋅)-β′g21-β(x)∂x(⋅).
Thus, in G¯a(·), there is no multiplicative part (no killing) and a shift in the drift, showing that the associated conditioned process (x̃a,t;t≥0) obeys the SDE (2.62)dx̃a=f̃a(x̃a)dt+g(x̃a)dwt,
with drift (2.63)f̃a(x)=f̃(x)+β′g2β(x)=f(x)+α′g2α(x)+β′g2β(x)=f(x)+g2(x)[α′α(x)+β′β(x)].
This process is ultimately absorbed at {0,1}.InG¯∂(·), there is a killing multiplicative part which is enhanced d/(1-β)>d and a shift in the drift, showing that the associated conditioned process (x̃∂,t;t≥0) exhibits a faster killing rate, but the drift shift guarantees that (x̃∂,t;t≥0) is not absorbed at the boundaries. We have (2.64)f̃∂(x)=f(x)+g2(x)[α′α(x)-β′1-β(x)].Additive Functionals of the Transformed Process.
for the new process(x̃t;t≥0), it is also of interest to evaluate additive functionals along their own sample paths. Let then α̃(x):=E(∫0τ¯xc̃(x̃s)ds+d̃(x̃τ¯x)) be such an additive functional where the functions c̃ and d̃ are themselves both nonnegative. It solves
(2.65)-G¯(α̃)=c̃ifx∈I,α̃=d̃ifx∈∂I.
Then, recalling the expression of the Green function 𝔤(x,y) of (xt;t≥0) in (2.22), we find explicitly
(2.66)α̃(x)=1α(x)∫Ig(x,y)α(y)c̃(y)dy.Specific transformations of interest.
(i) The case c=0 deserves a special treatment. Indeed, in this case, τ̃x,∂=∞ and so τ¯x:=τ̃x, the absorption time for the process (x̃t;t≥0) governed by the new SDE (2.38). Here, G¯=G̃. Assuming α solves -G(α)=0 if x∈I with boundary conditions α(0)=0 and α(1)=1 (α(0)=1 and α(1)=0,resp.), the new process (x̃t;t≥0) is just (xt;t≥0) conditioned on exiting at x=1 (at x=0,resp.). In the first case, the boundary 1 is exit, whereas 0 is entrance; α reads
(2.67)α(x)=∫0xe-2∫0y(f(z)/(g2(z)))dzdy∫01e-2∫0y(f(z)/(g2(z)))dzdy,
with
(2.68)f̃(x)=f(x)+g2(x)e-2∫0x(f(z)/(g2(z)))dz∫0xe-2∫0y(f(z)/(g2(z)))dzdy,
giving the new drift. In the second case, α(x)=(∫x1e-2∫0y(f(z)/g2(z))dzdy/(∫01e-2∫0y(f(z)/g2(z))dzdy)), and the boundary 0 is exit, whereas 1 is entrance. Thus, τ̃x is just the exit time at x=1 (at x=0,resp.). Let α̃(x):=E(τ̃x). Then, α̃(x) solves -G̃(α̃)=1, whose explicit solution is
(2.69)α̃(x)=1α(x)∫Ig(x,y)α(y)dy,
in terms of 𝔤(x,y), the Green function of (xt;t≥0).
Example 2.1.
Consider the WF model on[0,1] with selection for which, with σ∈R,f(x)=σx(1-x) and g2(x)=x(1-x). Assume that  α solves -G(α)=0 if x∈(0,1) with α(0)=0 and α(1)=1; one gets, α(x)=(1-e-2σx)/(1-e-2σ). The diffusion corresponding to (2.38) has the new drift: f̃(x)=σx(1-x)coth(σx), independently of the sign of σ. It models the WF diffusion with selection conditioned on exit at ∘=1.
(ii) Assume that α now solves -G(α)=1 if x∈I with boundary conditions α(0)=α(1)=0. In this case study, one selects sample paths of (xt;t≥0) with a large mean absorption time α(x)=E(τx). Sample paths with large sojourn time in I are favored. We have
(2.70)α(x)=∫Ig(x,y)dy,
where 𝔤(x,y) is the Green function (2.22). The boundaries of (x̃t;t≥0) are now both entrance boundaries and so τ̃x=∞. (x̃t;t≥0) is not absorbed at the boundaries. The stopping time τ¯x of (x̃t;t≥0) is just its killing time τ̃x,∂. Let α̃(x):=E(τ̃x,∂). Then, α̃(x) solves -G¯(α̃)=1, α̃(0)=α̃(1)=0, with explicit solution
(2.71)α̃(x)=1α(x)∫Ig(x,y)α(y)dy.
(iii) Assume that α now solves -G(α)=δy(x) if x∈I with boundary conditions α(0)=α(1)=0. In this case study, one selects sample paths of (xt;t≥0) with a large sojourn time density at y recalling α(x)=:𝔤(x,y)=E(∫0τxδy(xs)ds). The stopping time τ̃y(x) of (x̃t;t≥0) occurs at rate δy(x)/𝔤(x,y). It is a killing time when the process is at y for the last time after a geometrically distributed number of passages there with rate 1/𝔤(x,y) (or with success probability 1/(1+𝔤(x,y))). Let α̃y(x):=E(τ̃y(x)). Then, α̃y(x) solves -G¯(α̃)=1, with explicit solution
(2.72)α̃y(x)=1g(x,y)∫Ig(x,z)g(z,y)dz.
when x=1/N,α̃y(1/N) may be viewed as the age of a mutant currently observed to present frequency y, see [18].
The Green function aty0∈(0,1) of the transformed process (x̃t;t≥0) is 𝔤̃y(x,y0) solution to -G¯(𝔤̃y)=δy0(x). It takes the simple form
(2.73)g̃y(x,y0)=1g(x,y)∫Ig(x,z)g(z,y)δy0(z)dz=g(y0,y)g(x,y)g(x,y0).
(iv) Let λ1 be the smallest non-null eigenvalue of the infinitesimal generator G. Let α=u1 be the corresponding eigenvector, that is, satisfying, -Gu1=λ1u1 with boundary conditions u1(0)=u1(1)=0. Then, c=λu1. The new KB operator associated to the transformed process (x¯t;t≥0) is
(2.74)G¯(⋅)=1αG(α)⋅+G̃(⋅)=-λ1⋅+G̃(⋅),
obtained while killing the sample paths of the process (x̃t;t≥0) governed by G̃ at constant death rate d=λ1. The transition probability of the transformed stochastic process (x¯t;t≥0) is
(2.75)p¯(x;t,y)=u1(y)u1(x)p(x;t,y).
Define p̃(x;t,y)=eλ1tp¯(x;t,y). It is the transition probability of the process (x̃t;t≥0) governed by G̃; it corresponds to the original process (xt;t≥0) conditioned on never hitting the boundaries {0,1} (the so-called Q-process of (xt;t≥0), see [19]). It is simply obtained from (xt;t≥0) by adding the additional drift term (u1′/u1)g2 to f, where u1 is the eigenvector of G associated to its smallest non-null eigenvalue. The determination of α=u1 is a Sturm-Liouville problem. When t is large, to the dominant order
(2.76)p(x;t,y)~e-λ1tu1(x)v1(y)〈u1,v1〉,
where v1 is the Yaglom limit law of (xt;t≥0). Therefore
(2.77)p̃(x;t,y)~eλ1tu1(y)u1(x)e-λ1tu1(x)v1(y)〈u1,v1〉=u1(y)v1(y)〈u1,v1〉.
Thus, the limit law of the Q-process (x̃t;t≥0) is the normalized Hadamard product of the eigenvectors u1 and v1 associated, respectively, to G and G*. On the other hand, the limit law of (x̃t;t≥0) is directly given by
(2.78)p̃(x;t,y)⟶t→∞p̃(y)=1Zg2(y)e2∫0y((f(z)+((u1′/u1)g2)(z))/(g2(z)))dz=u12(y)Zg2(y)e2∫0y(f(z)/(g2(z)))dz,,
where Z is the appropriate normalizing constant. Comparing (2.77) and (2.78)
(2.79)v1(y)=u1(y)g2(y)e2∫0y(f(z)/(g2(z)))dz=u1(y)m(y).
The eigenvector v1 associated to G* is, therefore, equal to the eigenvector u1 associated to G times the speed density of (xt;t≥0).
When dealing for example with the neutral Wright-Fisher diffusion, it is known thatλ1=1 with u1=x(1-x) and v1≡1 (see Section 4.3, example (ii)). The Q-process (x̃t;t≥0) in this case obeys
(2.80)dx̃t=(1-2x̃t)dt̅+x̃t(1-x̃t)dwt,
with the stabilizing drift toward 1/2: f̃(x)=(u1′/u1)g2(x)=1-2x.The limit law of the Q-process (x̃t;t≥0) in this case is 6y(1-y). 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 generatorG¯(·)=-d(x)·+G̃(·) and applying the reciprocal Doob transform defined by (2.81)α̃(x)=1α(x)
leads to α̃(x)-1G¯(α̃(x)·)=G(·). Indeed, (2.82)α̃(x)-1G¯(α̃(x)⋅)=-d(x)⋅+α̃(x)-1G̃(α̃(x)⋅),α̃(x)-1G̃(α̃(x)⋅)=α̃(x)-1G̃(α̃(x))⋅+G̃(⋅)+1α̃g2α̃′∂x=d(x)⋅+G(⋅),
because (2.83)α̃(x)-1G̃(α̃(x))=+d(x),1α̃g2α̃′∂x=αg2(1α)′∂x=-α′αg2∂x.
Note that G¯(α̃(x))=G̃(α̃(x))-d(x)α̃(x)=0.This suggests that starting from a diffusion process with infinitesimal generatorG̃ (without its killing part) and applying the reciprocal Doob transform α̃(x)=1/(α(x)), one ends up with a modified process whose infinitesimal generator is (2.84)G̃¯(⋅)=G(⋅)+b(x),
where G(·)=G̃(·)-(α′/α)g2(x)∂x and (2.85)b(x)=α̃(x)-1G̃(α̃(x))=d(x)>0
is now a pure birth rate. Note that α̃ is now a subharmonic function for G̃ because -G̃(α̃)=-bα̃≤0.Letβ̃(x):=α̃(x)-1=α(x). Because G(α(x))=-α(x)b(x), we have (2.86)G̃¯(β̃(x))=G(β̃(x))+b(x)β̃(x)=0,
and so β̃(x)>0 is harmonic for G̃¯: Doob-transforming G̃¯ using β̃, we get (2.87)β̃-1G̃¯(β̃⋅)=b(x)⋅+β̃-1G(β̃⋅)=G̃(⋅),
which is the infinitesimal generator of the original diffusion process.Under our assumptions, both processxt and x̃t with infinitesimal generators G and G̃ are nonconservative diffusion processes with absorbing barriers. Further, b(x)>0 is bounded above. Therefore, b(x) may be written as (2.88)b(x)=b*(μ(x)-1),
where b*=maxxb(x)>0 and μ(x)∈[1,2].The process with infinitesimal generatorG̃¯ is now a pure binary branching diffusion process. For this class of models, an initial particle started at x obeys a diffusion process with infinitesimal generator G, absorbed when it hits the boundaries. At some random (mean b*) exponential time, this particle dies, giving birth in the process to a random number M(x) (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 EM(x)=μ(x).The process with infinitesimal generatorG̃¯ is, thus, a branching diffusion with supercritical binary splitting mechanism (μ(x)>1). 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.LetNt(x) be the global number of particles which are alive in the system at each time t, descending from an Eve particle started at x, and let (2.89)T(x)=inf⁡(t>0:Nt(x)=0)
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 (2.90)P(T(x)<∞,Nt(x)=0∀t≥T(x))=1,
uniformly in x. 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 numberN 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 withN individuals. The initial reproduction law is defined as follows: let |kN|:=∑m=1Nkm=N and kN:=(k1,…,kN) be integers. Assume that the first-generation random offspring numbers νN:=(νN(1),…,νN(N)) admit the following joint exchangeable polynomial distribution on the discrete simplex |kN|=N: (3.1)P(νN=kN)=N!⋅N-N∏n=1Nkn!.
This distribution can be obtained by conditioning N independent Poisson distributed random variables on summing to N. Assume subsequent iterations of this reproduction law are independent so that the population is with constant size for all generations.LetNr(n) be the offspring number of the n first individuals at the discrete generation r∈N0 corresponding to (say) allele A1 (the remaining number N-Nr(n) counts the number of alleles A2 at generation r). This sibship process is a discrete-time Markov chain with binomial transition probability given by (3.2)P(Nr+1(n)=k′∣Nr(n)=k)=(Nk′)(kN)k′(1-kN)N-k′.
Assume next that n=[Nx], where x∈(0,1). Then, as well known, the dynamics of the continuous space-time rescaled process xt:=N[Nt](n)/N, t∈R+ can be approximated for large N, to the leading term in N-1, by a Wright-Fisher-Itô diffusion on [0,1] (the purely random genetic drift case)(3.3)dxt=xt(1-xt)dwt,x0=x.
Here, (wt;t≥0) is a standard Wiener process. For this scaling limit process, a unit laps of time t=1 corresponds to a laps of time N for the original discrete-time process, thus time is measured in units of N. If the initial condition is x=N-1,xt is the diffusion approximation of the offspring frequency of a singleton at generation [Nt].Equation (3.3) is a one-dimensional diffusion as in (2.1) on [0,1], with zero drift f(x)=0 and volatility g(x)=x(1-x). This diffusion is already in natural coordinate, and so φ(x)=x. The scale function is x and the speed measure [x(1-x)]-1dx. One can check that both boundaries are exit in this case: the stopping time is τx=τx,0∧τx,1 where τx,0 is the extinction time and τx,1 the fixation time. The corresponding infinitesimal generators are G(·)=(1/2)x(1-x)∂x2(·) and G*(·)=(1/2)∂y2(y(1-y)·).


### 3.2. Nonneutral Cases

Two alleles Wright-Fisher models (with non-null drifts) can be obtained by considering the binomial transition probabilities bin(N,pN)(3.4)P(Nr+1(n)=k′∣Nr(n)=k)=(Nk′)(pN(kN))k′(1-pN(kN))N-k′,
where (3.5)pN(x):x∈(0,1)⟶(0,1)
is now some state-dependent probability (which is different from the identity x) reflecting some deterministic evolutionary drift from the allele A1 to the allele A2. For each r, we have (3.6)E(Nr+1(n)∣Nr(n)=k)=NpN(kN),σ2(Nr+1(n)∣Nr(n)=k)=NpN(kN)(1-pN(kN)),
which is amenable to a diffusion approximation in terms of xt:=N[Nt](n)/N, t∈R+ under suitable conditions.For instance, takingpN(x)=(1-π2,N)x+π1,N(1-x), where (π1,N,π2,N) are small (N-dependent) mutation probabilities from A1 to A2 (A2 to A1,resp.). Assuming that (N·π1,N,N·π2,N)→N→∞(u1,u2), leads after scaling to the drift of WF model with positive mutations rates (u1,u2).Taking(3.7)pN(x)=(1+s1,N)x(1+s1,Nx+(1-x)(1+s2,N)),
where si,N>0 are small N-dependent selection parameter satisfying N·si,N→N→∞σi>0,i=1,2, leads, after scaling, to the WF model with selective drift σx(1-x), where σ:=σ1-σ2. Essentially, the drift f(x) is a large N approximation of the bias: N(pN(x)-x). The WF diffusion with selection is thus (3.8)dxt=σxt(1-xt)dt+xt(1-xt)dwt,
where time is measured in units of N. Letting θt=Nt define a new time scale with inverse tθ=θ/N, the time-changed process yθ=xθ/N now obeys the SDE (3.9)dyθ=syθ(1-yθ)dθ+1Nyθ(1-yθ)dwθ,
with a small diffusion term. Here, s=s1-s2 and time θ is the usual time clock.The WF diffusion with selection (3.8) tends to drift to ∘=1 (∘=0,resp.) if allele A1 is selectively advantageous over A2:σ1>σ2 (σ1<σ2,resp.) in the following sense: if σ>0 (<0,resp.), the fixation probability at ∘=1, which is [15] (3.10)P(τx,1<τx,0)=1-e-2σx1-e-2σ
increases (decreases) with σ taking larger (smaller) values. Putting x=1/N, the fixation probability at 1 of an allele A1 mutant is of order: 2σ/N; see [15].


## 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 withN individuals. The initial reproduction law is defined as follows: let |kN|:=∑m=1Nkm=N and kN:=(k1,…,kN) be integers. Assume that the first-generation random offspring numbers νN:=(νN(1),…,νN(N)) admit the following joint exchangeable polynomial distribution on the discrete simplex |kN|=N: (3.1)P(νN=kN)=N!⋅N-N∏n=1Nkn!.
This distribution can be obtained by conditioning N independent Poisson distributed random variables on summing to N. Assume subsequent iterations of this reproduction law are independent so that the population is with constant size for all generations.LetNr(n) be the offspring number of the n first individuals at the discrete generation r∈N0 corresponding to (say) allele A1 (the remaining number N-Nr(n) counts the number of alleles A2 at generation r). This sibship process is a discrete-time Markov chain with binomial transition probability given by (3.2)P(Nr+1(n)=k′∣Nr(n)=k)=(Nk′)(kN)k′(1-kN)N-k′.
Assume next that n=[Nx], where x∈(0,1). Then, as well known, the dynamics of the continuous space-time rescaled process xt:=N[Nt](n)/N, t∈R+ can be approximated for large N, to the leading term in N-1, by a Wright-Fisher-Itô diffusion on [0,1] (the purely random genetic drift case)(3.3)dxt=xt(1-xt)dwt,x0=x.
Here, (wt;t≥0) is a standard Wiener process. For this scaling limit process, a unit laps of time t=1 corresponds to a laps of time N for the original discrete-time process, thus time is measured in units of N. If the initial condition is x=N-1,xt is the diffusion approximation of the offspring frequency of a singleton at generation [Nt].Equation (3.3) is a one-dimensional diffusion as in (2.1) on [0,1], with zero drift f(x)=0 and volatility g(x)=x(1-x). This diffusion is already in natural coordinate, and so φ(x)=x. The scale function is x and the speed measure [x(1-x)]-1dx. One can check that both boundaries are exit in this case: the stopping time is τx=τx,0∧τx,1 where τx,0 is the extinction time and τx,1 the fixation time. The corresponding infinitesimal generators are G(·)=(1/2)x(1-x)∂x2(·) and G*(·)=(1/2)∂y2(y(1-y)·).


## 3.2. Nonneutral Cases

Two alleles Wright-Fisher models (with non-null drifts) can be obtained by considering the binomial transition probabilities bin(N,pN)(3.4)P(Nr+1(n)=k′∣Nr(n)=k)=(Nk′)(pN(kN))k′(1-pN(kN))N-k′,
where (3.5)pN(x):x∈(0,1)⟶(0,1)
is now some state-dependent probability (which is different from the identity x) reflecting some deterministic evolutionary drift from the allele A1 to the allele A2. For each r, we have (3.6)E(Nr+1(n)∣Nr(n)=k)=NpN(kN),σ2(Nr+1(n)∣Nr(n)=k)=NpN(kN)(1-pN(kN)),
which is amenable to a diffusion approximation in terms of xt:=N[Nt](n)/N, t∈R+ under suitable conditions.For instance, takingpN(x)=(1-π2,N)x+π1,N(1-x), where (π1,N,π2,N) are small (N-dependent) mutation probabilities from A1 to A2 (A2 to A1,resp.). Assuming that (N·π1,N,N·π2,N)→N→∞(u1,u2), leads after scaling to the drift of WF model with positive mutations rates (u1,u2).Taking(3.7)pN(x)=(1+s1,N)x(1+s1,Nx+(1-x)(1+s2,N)),
where si,N>0 are small N-dependent selection parameter satisfying N·si,N→N→∞σi>0,i=1,2, leads, after scaling, to the WF model with selective drift σx(1-x), where σ:=σ1-σ2. Essentially, the drift f(x) is a large N approximation of the bias: N(pN(x)-x). The WF diffusion with selection is thus (3.8)dxt=σxt(1-xt)dt+xt(1-xt)dwt,
where time is measured in units of N. Letting θt=Nt define a new time scale with inverse tθ=θ/N, the time-changed process yθ=xθ/N now obeys the SDE (3.9)dyθ=syθ(1-yθ)dθ+1Nyθ(1-yθ)dwθ,
with a small diffusion term. Here, s=s1-s2 and time θ is the usual time clock.The WF diffusion with selection (3.8) tends to drift to ∘=1 (∘=0,resp.) if allele A1 is selectively advantageous over A2:σ1>σ2 (σ1<σ2,resp.) in the following sense: if σ>0 (<0,resp.), the fixation probability at ∘=1, which is [15] (3.10)P(τx,1<τx,0)=1-e-2σx1-e-2σ
increases (decreases) with σ taking larger (smaller) values. Putting x=1/N, the fixation probability at 1 of an allele A1 mutant is of order: 2σ/N; see [15].


## 4. The Neutral WF Model

In this section, we particularize the general ideas developed in the introductory Section2 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.Withz∈(-1,1), let (Pk(z);k≥0) be the degree-(k+1) Gegenbauer polynomials solving (1-z2)Pk′′(z)+k(k+1)Pk(z)=0 with Pk′(±1)=∓(1/2),k≥1; we let P0(z):=(1-z)/2. When k≥1, we have Pk(±1)=0 and so Pk(z)=(1-z2)Qk(z), where Qk(z) is a polynomial with degree k-1 satisfying Qk(-1)=(-1)k-1 and Qk(1)=1. With x∈(0,1), let (uk(x);k≥0) be defined by: uk(x)=Pk(1-2x). These polynomials clearly constitute a system of eigenfunctions for the KB operator -G=-(1/2)x(1-x)∂x2 with eigenvalues λk=(k(k+1))/2,k≥0, thus with -G(uk(x))=λkuk(x). In particular, u0(x)=x, u1(x)=x-x2, u2(x)=x-3x2+2x3,  u3(x)=x-6x2+10x3-5x4,  u4(x)=x-10x2+30x3-35x4+14x5,…. With k≥1, we have uk(0)=uk(1)=0 and uk′(0)=1 and uk′(1)=-1.The eigenfunctions of the KF operatorG*(·)=(1/2)∂x2[y(1-y)·] are given by vk(y)=m(y)·uk(y),k≥0, where the Radon measure of weights m(y)dy is the speed measure: m(y)dy=dy/(y(1-y)), for the same eigenvalues. For instance, v0(y)=1/(1-y), v1(y)=1,  v2(y)=1-2y,  v3(y)=1-5y+5y2,v4(y)=1-9y+21y2-14y3,….Althoughλ0=0 really constitutes an eigenvalue, only v0(y) is not a polynomial. When k≥1, from their definition, the uk(x) polynomials satisfy uk(0)=uk(1)=0 in such a way that vk(y)=m(y)·uk(y),k≥1 is a polynomial with degree k-1.Let〈f,g〉m=∫01f(x)g(x)m(x)dx. We note that 〈vj,uk〉=〈uj,uk〉m=0 if j≠k and the system uk(x);k≥1 is a complete orthogonal set of eigenvectors. Therefore, for any square-integrable function ψ(x)∈L2([0,1],m(y)dy) admitting a decomposition in the basis uk(x),k≥1.(4.1)Eψ(xt)=∑k≥1cke-λktuk(x),whereck=〈ψ,uk〉m〈vk,uk〉,
where ψ(x)=∑k≥1ckuk(x). This series expansion solves the KBE: ∂tu=G(u); u(x,0)=ψ(x) where u=u(x,t):=Eψ(xt).Moreover, the transition probability densityp(x;t,y) of the neutral WF models admits the spectral expansion (4.2)p(x;t,y)=∑k≥1bke-λktuk(x)vk(y),wherebk=1〈vk,uk〉.
Starting from x, the cumulated probability masses by time t at the exit boundaries {0,1} are, respectively, (see [17]) (4.3)12∫0tp(x;s,0)ds=∑k≥1bk2λk(1-e-λkt)uk(x)vk(0),12∫0tp(x;s,1)ds=∑k≥1bk2λk(1-e-λkt)uk(x)vk(1),
which tend as t→∞ toward the extinction and fixation probabilities, namely, here P(τx,0<∞)=P(τx,0<τx,1)=1-x and P(τx,1<∞)=x. Because vk(0)=1 and vk(1)=(-1)k-1, we get the identities (4.4)1-x=∑k≥1bk2λkuk(x)or1x=∑k≥1bk2λkvk(x),x=∑k≥1(-1)k-1bk2λkuk(x)or11-x=∑k≥1(-1)k-1bk2λkvk(x),
leading to the relationship  ∑k≥1((-1)k-1bk/2λk)=1.The series expansion forp(x;t,y) solves the KFE of the WF model. The transition density p(x;t,y) is reversible with respect to the speed density since for 0<x,y<1(4.5)m(x)p(x;t,y)=m(y)p(y;t,x)=∑k≥1bke-λktvk(x)vk(y).
The measures vk(y)dy, k≥1 are not probability measures because the vk(y) are not necessarily positive over [0,1]. This decomposition is not a mixture. We have 〈vk,uk〉=∥uk∥2,m2 the 2-norm for the weight function m. We notice that 〈v0,u0〉=∫01(y/(1-y))dy=∞ so that c0=b0=0 although λ0=0 is indeed an eigenvalue, the above sums should be started at k=1 (expressing the lack of an invariant measure for the WF model as a result of its absorption at the boundaries).We haveP(τx>t)=∫01P(xt∈dy) and so (4.6)ρt(x):=P(τx>t)=∑k≥1∫01vk(y)dy〈vk,uk〉e-λktuk(x)
is the exact tail distribution of the absorption time.Sincev1(y)=1, to the leading order in t, for large time (4.7)P(xt∈dy)=6e-t⋅x(1-x)dy+O(e-3t),
which is independent of y. Integrating over y, ρt(x):=P(τx>t)~6e-t·x(1-x) so that the conditional probability (4.8)P(xt∈dy∣τx>t)~t→∞dy
is asymptotically uniform in the Yaglom limit. As time passes by, given absorption did not occur in the past, xt→dx∞ (as t→∞) which is a uniformly distributed random variable on [0,1].


### 4.2. Additive Functionals for the Neutral WF

Let(xt;t≥0) be the WF diffusion model defined by (3.3) on the interval I=[0,1], where both endpoints are absorbing (exit). We wish to evaluate the additive quantities (4.9)α(x)=E(∫0τxc(xs)ds+d(xτx)),
where functions c and d are both nonnegative. With G=(1/2)x(1-x)∂x2, α(x) solves (4.10)-G(α)=cifx∈I,α=difx∈∂I.
Take c=limε↓0(1/2ε)1(x∈(y-ε,y+ε))=:δy(x) and d=0, when y∈I: in this case, α:=𝔤(x,y) is the Green function. The solution takes the simple form (4.11)g(x,y)=2xyifx<y,g(x,y)=21-x1-yifx>y.
The Green function solves the above general problem of evaluating additive functionals α(x)(4.12)α(x)=∫Ig(x,y)c(y)dyifx∈I,α=difx∈∂I.As a Few Examples
(1) Letc=1 and d=0: here, α(x)=E(τx) is the mean time of absorption (average time spent in I before absorption). The solution is (the Crow and Kimura formula, see [2])
(4.13)α(x)=2x∫x1dyy+2(1-x)∫0xdy1-y=-2(xlog⁡x+(1-x)log⁡(1-x)).
(2) Letc=0 and d(∘)=1(∘=1). Let α(x)=P(xtfirsthits[0,1]at1|x0=x). Then, α(x) is a G-harmonic function solution to G(α)=0, with boundary conditions α(0)=0 and α(1)=1. The solution for WF model is: α(x)=x. Stated differently, x=P(τx,1<τx,0) is the probability that the exit time at ∘=1 is less than the one at ∘=0, starting from x.
On the contrary, choosingα(x) to be a G-harmonic function with boundary conditions α(0)=1 and α(1)=0, α(x)=P(xtfirsthits[0,1]at0|x0=x)=1-x. Thus, 1-x=P(τx,0<τx,1).
(3) Letc(xs)=2xs(1-xs) measure the heterozygosity of the WF process at time s and assume d(0)=d(1)=1. A remarkable thing is that the average heterozygosity over the sample paths is
(4.14)α(x)=E(∫0τxc(xs)ds)=4x∫x1(1-y)dy+4(1-x)∫0xydy=2x(1-x),
which is the initial heterozygosity of the population.


### 4.3. Transformation of WF Sample Paths, [3]

Withp(x;t,y) the transition probability density of WF model, define a new α-transformed stochastic process (x¯t;t≥0) by its transition probability (4.15)p¯(x;t,y)=α(y)α(x)p(x;t,y).(i) Conditioning WF on exit at some boundary. Assume first α solves -G(α)=0 with boundary conditions α(0)=0 and α(1)=1; hence, α reads α(x)=x. In this case, τ̃x,∂=∞ (no killing), and so τ¯x:=τ̃x is the absorption time for a process (x̃t;t≥0) governed by a new SDE with a drift term. The new process (x̃t;t≥0) is just (xt;t≥0) conditioned on exiting at ∘=1. The boundary 1 is exit whereas 0 is entrance. Thus, the model for (x̃t;t≥0) becomes dx̃t=(1-x̃t)dt+x̃t(1-x̃t)dwt, x̃0=x∈(0,1) now with linear drift f̃(x)=1-x and g(x)=x(1-x). Its transition probability is (4.16)p¯1(x;t,y)=yxp(x;t,y),
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 -G(α)=0 if x∈I with boundary conditions α(0)=1 and α(1)=0, the new process (x̃t;t≥0) is just (xt;t≥0) conditioned on exiting at x=0. Boundary 0 is exit, whereas 1 is entrance; in this case, α is α(x)=1-x. Thus, the model for (x̃t;t≥0) becomes dx̃t=-x̃tdt+x̃t(1-x̃t)dwt, x̃0=x∈(0,1) with f̃(x)=-x and g(x)=x(1-x). Its transition probability is (4.17)p¯0(x;t,y)=1-y1-xp(x;t,y),
where the subscript 0 indicates that this is the conditional transition probability of WF sample paths whose exit now is at ∘=0. Recalling that, starting from x, (xt;t≥0) gets absorbed at ∘=1 (0,resp.) with probability x (1-x,resp.), we recover that (4.18)p(x;t,y)=x⋅p¯1(x;t,y)+(1-x)⋅p¯0(x;t,y).
Using the solution to KFE for p, we obtain an expression for both p¯1(x;t,y) and p¯0(x;t,y), simply by premultiplying it by the corresponding right factor. Integrating the results over y, we get the conditional tail distributions of the exit times at ∘=1 or 0, given the exit is at ∘=1 or 0.Exploiting the large time behavior ofp(x;t,y), to the first order in t, we get (4.19)p¯1(x;t,y)~6e-t⋅(1-x)y,p¯0(x;t,y)=6e-t⋅x(1-y).
Integrating over y, ρ¯t,1(x):=P1(τ̃x>t)~3e-t·(1-x) and ρ¯t,0(x):=P0(τ̃x>t)~3e-t·x 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 (4.20)P1(x̃t∈dy∣τ̃x>t)~2ydy,P0(x̃t∈dy∣τ̃x>t)~2(1-y)dy,
where we recognize the densities of specific beta-distributed random variables. Specifically, we conclude that as time passes by, given absorption occurs at ∘=1 and given it has not occurred in the past, x̃t→d beta(2,1) distribution on [0,1]. Similarly, given absorption occurs at ∘=0 and given it has not occurred previously, x̃t→d beta(1,2) distribution on [0,1].In the previously displayed formula,τ̃x is just the exit time at ∘=1 (at ∘=0,resp.) of the conditional transformed WF diffusions. Let α̃(x):=E(τ̃x). Then, with G¯(·)=(1/α(x))G(α(x)·), α̃(x) solves -G¯(α̃)=1, whose explicit solution is (4.21)α̃(x)=1α(x)∫01g(x,y)α(y)dy,
in terms of 𝔤(x,y), the Green function of (xt;t≥0). For the WF model conditioned on exit at ∘=1 (0,resp.), we find, respectively, the Kimura and Ohta's formulae in [23] (4.22)α̃1(x)=-2x(1-x)log⁡(1-x),α̃0(x)=-21-xxlog⁡x.
This result could have been guessed by observing that xα̃1(x)+(1-x)α̃0(x) is the expected absorption time of the original WF model. When x→0+, (resp.,  x→1-), it takes an average time 2 to reach 1 (0,resp.) for the WF model conditioned on exit at ∘=1 (0,resp.).(ii) Selection of WF sample paths with large heterozygosity. Assume that α now solves -G(α)=2x(1-x) if x∈I with boundary conditions α(0)=α(1)=0. Then, α=2x(1-x) and this α is the right eigenvector of -G associated to the smallest positive eigenvalue λ1=1 of the neutral WF model. In this case study, one selects sample paths of (xt;t≥0) with large heterozygosity. The dynamics of (x̃t;t≥0) in (2.38) is (4.23)dx̃t=(1-2x̃t)dt+x̃t(1-x̃t)dwt,
subject to a constant killing rate 1. The boundaries of (x̃t;t≥0) are now both entrance boundaries and so τ̃x=∞. (x̃t;t≥0) is not absorbed at the boundaries. The stopping time τ¯x of (x̃t;t≥0) is just its killing time τ̃x,∂ which is mean 1 exponentially distributed, independently of the starting point x. Indeed, (4.24)P(τ̃x,∂>t)=∫01p¯(x;t,y)dy=1α(x)∫01α(y)p(x;t,y)dy=1x(1-x)∑k≥1bke-λktuk(x)∫01y(1-y)vk(y)dy=∑k≥1bke-λktvk(x)∫01uk(y)dy=6e-t16=e-t,
recalling x(1-x)vk(x)=uk(x) and observing ∫01uk(y)dy=0 if k≥2.As time passes, killing ofx̃t occurs, and given killing will never occur in the future, x̃t→dx̃∞ a random variable with density 6y(1-y) on [0,1] which is a beta(2,2) density. In this selection of paths procedure, the conditional density of (x̃t;t≥0) given τ̃x,∂=∞ is indeed p̃(x;t,y)=etp¯(x;t,y), where p¯(x;t,y)=(y(1-y)/(x(1-x)))p(x;t,y). Therefore, (4.25)p̃(x;t,y)=∑k≥1bke-(λk-1)tuk(x)x(1-x)y(1-y)vk(y).
Recalling u1(x)=x(1-x),v1(y)=1 and b1=6, we get p̃(x;t,y)→t→∞6y(1-y), regardless of the initial condition x. This is the beta(2,2) limit law of the Q-process of the neutral WF diffusion.(iii) Selection of WF sample paths with large sojourn time density at y. Assume now that α solves -G(α)=δy(x) if x∈I and so α(x)=:𝔤(x,y). Using the Green function of the neutral WF model, the transition probability density of (x¯t;t≥0) is (4.26)p¯(x;t,y)=yxp(x;,y)ifx<y,p¯(x;t,y)=1-y1-xp(x;t,y)ifx>y.
Thus, given x<y (x>y), (x¯t;t≥0) coincides with (xt;t≥0) conditioned to exit in 1 (0,resp.) killed at rate δy(x) when it passes through y, necessarily at some time.The stopping timeτ̃y(x) of (x̃t;t≥0) is just its killing time when the process is at y for the last time with a geometrically number of passages at y with rate 1 (or success probability 1/2).


## 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.Withz∈(-1,1), let (Pk(z);k≥0) be the degree-(k+1) Gegenbauer polynomials solving (1-z2)Pk′′(z)+k(k+1)Pk(z)=0 with Pk′(±1)=∓(1/2),k≥1; we let P0(z):=(1-z)/2. When k≥1, we have Pk(±1)=0 and so Pk(z)=(1-z2)Qk(z), where Qk(z) is a polynomial with degree k-1 satisfying Qk(-1)=(-1)k-1 and Qk(1)=1. With x∈(0,1), let (uk(x);k≥0) be defined by: uk(x)=Pk(1-2x). These polynomials clearly constitute a system of eigenfunctions for the KB operator -G=-(1/2)x(1-x)∂x2 with eigenvalues λk=(k(k+1))/2,k≥0, thus with -G(uk(x))=λkuk(x). In particular, u0(x)=x, u1(x)=x-x2, u2(x)=x-3x2+2x3,  u3(x)=x-6x2+10x3-5x4,  u4(x)=x-10x2+30x3-35x4+14x5,…. With k≥1, we have uk(0)=uk(1)=0 and uk′(0)=1 and uk′(1)=-1.The eigenfunctions of the KF operatorG*(·)=(1/2)∂x2[y(1-y)·] are given by vk(y)=m(y)·uk(y),k≥0, where the Radon measure of weights m(y)dy is the speed measure: m(y)dy=dy/(y(1-y)), for the same eigenvalues. For instance, v0(y)=1/(1-y), v1(y)=1,  v2(y)=1-2y,  v3(y)=1-5y+5y2,v4(y)=1-9y+21y2-14y3,….Althoughλ0=0 really constitutes an eigenvalue, only v0(y) is not a polynomial. When k≥1, from their definition, the uk(x) polynomials satisfy uk(0)=uk(1)=0 in such a way that vk(y)=m(y)·uk(y),k≥1 is a polynomial with degree k-1.Let〈f,g〉m=∫01f(x)g(x)m(x)dx. We note that 〈vj,uk〉=〈uj,uk〉m=0 if j≠k and the system uk(x);k≥1 is a complete orthogonal set of eigenvectors. Therefore, for any square-integrable function ψ(x)∈L2([0,1],m(y)dy) admitting a decomposition in the basis uk(x),k≥1.(4.1)Eψ(xt)=∑k≥1cke-λktuk(x),whereck=〈ψ,uk〉m〈vk,uk〉,
where ψ(x)=∑k≥1ckuk(x). This series expansion solves the KBE: ∂tu=G(u); u(x,0)=ψ(x) where u=u(x,t):=Eψ(xt).Moreover, the transition probability densityp(x;t,y) of the neutral WF models admits the spectral expansion (4.2)p(x;t,y)=∑k≥1bke-λktuk(x)vk(y),wherebk=1〈vk,uk〉.
Starting from x, the cumulated probability masses by time t at the exit boundaries {0,1} are, respectively, (see [17]) (4.3)12∫0tp(x;s,0)ds=∑k≥1bk2λk(1-e-λkt)uk(x)vk(0),12∫0tp(x;s,1)ds=∑k≥1bk2λk(1-e-λkt)uk(x)vk(1),
which tend as t→∞ toward the extinction and fixation probabilities, namely, here P(τx,0<∞)=P(τx,0<τx,1)=1-x and P(τx,1<∞)=x. Because vk(0)=1 and vk(1)=(-1)k-1, we get the identities (4.4)1-x=∑k≥1bk2λkuk(x)or1x=∑k≥1bk2λkvk(x),x=∑k≥1(-1)k-1bk2λkuk(x)or11-x=∑k≥1(-1)k-1bk2λkvk(x),
leading to the relationship  ∑k≥1((-1)k-1bk/2λk)=1.The series expansion forp(x;t,y) solves the KFE of the WF model. The transition density p(x;t,y) is reversible with respect to the speed density since for 0<x,y<1(4.5)m(x)p(x;t,y)=m(y)p(y;t,x)=∑k≥1bke-λktvk(x)vk(y).
The measures vk(y)dy, k≥1 are not probability measures because the vk(y) are not necessarily positive over [0,1]. This decomposition is not a mixture. We have 〈vk,uk〉=∥uk∥2,m2 the 2-norm for the weight function m. We notice that 〈v0,u0〉=∫01(y/(1-y))dy=∞ so that c0=b0=0 although λ0=0 is indeed an eigenvalue, the above sums should be started at k=1 (expressing the lack of an invariant measure for the WF model as a result of its absorption at the boundaries).We haveP(τx>t)=∫01P(xt∈dy) and so (4.6)ρt(x):=P(τx>t)=∑k≥1∫01vk(y)dy〈vk,uk〉e-λktuk(x)
is the exact tail distribution of the absorption time.Sincev1(y)=1, to the leading order in t, for large time (4.7)P(xt∈dy)=6e-t⋅x(1-x)dy+O(e-3t),
which is independent of y. Integrating over y, ρt(x):=P(τx>t)~6e-t·x(1-x) so that the conditional probability (4.8)P(xt∈dy∣τx>t)~t→∞dy
is asymptotically uniform in the Yaglom limit. As time passes by, given absorption did not occur in the past, xt→dx∞ (as t→∞) which is a uniformly distributed random variable on [0,1].


## 4.2. Additive Functionals for the Neutral WF

Let(xt;t≥0) be the WF diffusion model defined by (3.3) on the interval I=[0,1], where both endpoints are absorbing (exit). We wish to evaluate the additive quantities (4.9)α(x)=E(∫0τxc(xs)ds+d(xτx)),
where functions c and d are both nonnegative. With G=(1/2)x(1-x)∂x2, α(x) solves (4.10)-G(α)=cifx∈I,α=difx∈∂I.
Take c=limε↓0(1/2ε)1(x∈(y-ε,y+ε))=:δy(x) and d=0, when y∈I: in this case, α:=𝔤(x,y) is the Green function. The solution takes the simple form (4.11)g(x,y)=2xyifx<y,g(x,y)=21-x1-yifx>y.
The Green function solves the above general problem of evaluating additive functionals α(x)(4.12)α(x)=∫Ig(x,y)c(y)dyifx∈I,α=difx∈∂I.As a Few Examples
(1) Letc=1 and d=0: here, α(x)=E(τx) is the mean time of absorption (average time spent in I before absorption). The solution is (the Crow and Kimura formula, see [2])
(4.13)α(x)=2x∫x1dyy+2(1-x)∫0xdy1-y=-2(xlog⁡x+(1-x)log⁡(1-x)).
(2) Letc=0 and d(∘)=1(∘=1). Let α(x)=P(xtfirsthits[0,1]at1|x0=x). Then, α(x) is a G-harmonic function solution to G(α)=0, with boundary conditions α(0)=0 and α(1)=1. The solution for WF model is: α(x)=x. Stated differently, x=P(τx,1<τx,0) is the probability that the exit time at ∘=1 is less than the one at ∘=0, starting from x.
On the contrary, choosingα(x) to be a G-harmonic function with boundary conditions α(0)=1 and α(1)=0, α(x)=P(xtfirsthits[0,1]at0|x0=x)=1-x. Thus, 1-x=P(τx,0<τx,1).
(3) Letc(xs)=2xs(1-xs) measure the heterozygosity of the WF process at time s and assume d(0)=d(1)=1. A remarkable thing is that the average heterozygosity over the sample paths is
(4.14)α(x)=E(∫0τxc(xs)ds)=4x∫x1(1-y)dy+4(1-x)∫0xydy=2x(1-x),
which is the initial heterozygosity of the population.


## 4.3. Transformation of WF Sample Paths, [3]

Withp(x;t,y) the transition probability density of WF model, define a new α-transformed stochastic process (x¯t;t≥0) by its transition probability (4.15)p¯(x;t,y)=α(y)α(x)p(x;t,y).(i) Conditioning WF on exit at some boundary. Assume first α solves -G(α)=0 with boundary conditions α(0)=0 and α(1)=1; hence, α reads α(x)=x. In this case, τ̃x,∂=∞ (no killing), and so τ¯x:=τ̃x is the absorption time for a process (x̃t;t≥0) governed by a new SDE with a drift term. The new process (x̃t;t≥0) is just (xt;t≥0) conditioned on exiting at ∘=1. The boundary 1 is exit whereas 0 is entrance. Thus, the model for (x̃t;t≥0) becomes dx̃t=(1-x̃t)dt+x̃t(1-x̃t)dwt, x̃0=x∈(0,1) now with linear drift f̃(x)=1-x and g(x)=x(1-x). Its transition probability is (4.16)p¯1(x;t,y)=yxp(x;t,y),
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 -G(α)=0 if x∈I with boundary conditions α(0)=1 and α(1)=0, the new process (x̃t;t≥0) is just (xt;t≥0) conditioned on exiting at x=0. Boundary 0 is exit, whereas 1 is entrance; in this case, α is α(x)=1-x. Thus, the model for (x̃t;t≥0) becomes dx̃t=-x̃tdt+x̃t(1-x̃t)dwt, x̃0=x∈(0,1) with f̃(x)=-x and g(x)=x(1-x). Its transition probability is (4.17)p¯0(x;t,y)=1-y1-xp(x;t,y),
where the subscript 0 indicates that this is the conditional transition probability of WF sample paths whose exit now is at ∘=0. Recalling that, starting from x, (xt;t≥0) gets absorbed at ∘=1 (0,resp.) with probability x (1-x,resp.), we recover that (4.18)p(x;t,y)=x⋅p¯1(x;t,y)+(1-x)⋅p¯0(x;t,y).
Using the solution to KFE for p, we obtain an expression for both p¯1(x;t,y) and p¯0(x;t,y), simply by premultiplying it by the corresponding right factor. Integrating the results over y, we get the conditional tail distributions of the exit times at ∘=1 or 0, given the exit is at ∘=1 or 0.Exploiting the large time behavior ofp(x;t,y), to the first order in t, we get (4.19)p¯1(x;t,y)~6e-t⋅(1-x)y,p¯0(x;t,y)=6e-t⋅x(1-y).
Integrating over y, ρ¯t,1(x):=P1(τ̃x>t)~3e-t·(1-x) and ρ¯t,0(x):=P0(τ̃x>t)~3e-t·x 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 (4.20)P1(x̃t∈dy∣τ̃x>t)~2ydy,P0(x̃t∈dy∣τ̃x>t)~2(1-y)dy,
where we recognize the densities of specific beta-distributed random variables. Specifically, we conclude that as time passes by, given absorption occurs at ∘=1 and given it has not occurred in the past, x̃t→d beta(2,1) distribution on [0,1]. Similarly, given absorption occurs at ∘=0 and given it has not occurred previously, x̃t→d beta(1,2) distribution on [0,1].In the previously displayed formula,τ̃x is just the exit time at ∘=1 (at ∘=0,resp.) of the conditional transformed WF diffusions. Let α̃(x):=E(τ̃x). Then, with G¯(·)=(1/α(x))G(α(x)·), α̃(x) solves -G¯(α̃)=1, whose explicit solution is (4.21)α̃(x)=1α(x)∫01g(x,y)α(y)dy,
in terms of 𝔤(x,y), the Green function of (xt;t≥0). For the WF model conditioned on exit at ∘=1 (0,resp.), we find, respectively, the Kimura and Ohta's formulae in [23] (4.22)α̃1(x)=-2x(1-x)log⁡(1-x),α̃0(x)=-21-xxlog⁡x.
This result could have been guessed by observing that xα̃1(x)+(1-x)α̃0(x) is the expected absorption time of the original WF model. When x→0+, (resp.,  x→1-), it takes an average time 2 to reach 1 (0,resp.) for the WF model conditioned on exit at ∘=1 (0,resp.).(ii) Selection of WF sample paths with large heterozygosity. Assume that α now solves -G(α)=2x(1-x) if x∈I with boundary conditions α(0)=α(1)=0. Then, α=2x(1-x) and this α is the right eigenvector of -G associated to the smallest positive eigenvalue λ1=1 of the neutral WF model. In this case study, one selects sample paths of (xt;t≥0) with large heterozygosity. The dynamics of (x̃t;t≥0) in (2.38) is (4.23)dx̃t=(1-2x̃t)dt+x̃t(1-x̃t)dwt,
subject to a constant killing rate 1. The boundaries of (x̃t;t≥0) are now both entrance boundaries and so τ̃x=∞. (x̃t;t≥0) is not absorbed at the boundaries. The stopping time τ¯x of (x̃t;t≥0) is just its killing time τ̃x,∂ which is mean 1 exponentially distributed, independently of the starting point x. Indeed, (4.24)P(τ̃x,∂>t)=∫01p¯(x;t,y)dy=1α(x)∫01α(y)p(x;t,y)dy=1x(1-x)∑k≥1bke-λktuk(x)∫01y(1-y)vk(y)dy=∑k≥1bke-λktvk(x)∫01uk(y)dy=6e-t16=e-t,
recalling x(1-x)vk(x)=uk(x) and observing ∫01uk(y)dy=0 if k≥2.As time passes, killing ofx̃t occurs, and given killing will never occur in the future, x̃t→dx̃∞ a random variable with density 6y(1-y) on [0,1] which is a beta(2,2) density. In this selection of paths procedure, the conditional density of (x̃t;t≥0) given τ̃x,∂=∞ is indeed p̃(x;t,y)=etp¯(x;t,y), where p¯(x;t,y)=(y(1-y)/(x(1-x)))p(x;t,y). Therefore, (4.25)p̃(x;t,y)=∑k≥1bke-(λk-1)tuk(x)x(1-x)y(1-y)vk(y).
Recalling u1(x)=x(1-x),v1(y)=1 and b1=6, we get p̃(x;t,y)→t→∞6y(1-y), regardless of the initial condition x. This is the beta(2,2) limit law of the Q-process of the neutral WF diffusion.(iii) Selection of WF sample paths with large sojourn time density at y. Assume now that α solves -G(α)=δy(x) if x∈I and so α(x)=:𝔤(x,y). Using the Green function of the neutral WF model, the transition probability density of (x¯t;t≥0) is (4.26)p¯(x;t,y)=yxp(x;,y)ifx<y,p¯(x;t,y)=1-y1-xp(x;t,y)ifx>y.
Thus, given x<y (x>y), (x¯t;t≥0) coincides with (xt;t≥0) conditioned to exit in 1 (0,resp.) killed at rate δy(x) when it passes through y, necessarily at some time.The stopping timeτ̃y(x) of (x̃t;t≥0) is just its killing time when the process is at y for the last time with a geometrically number of passages at y with rate 1 (or success probability 1/2).


## 5. The WF Model with Selection

Now, we focus on the diffusion process (3.8). Let (vk(y))k≥1 be the Gegenbauer eigen-polynomials of the KF operator corresponding to the neutral WF diffusion (3.3) and so with eigenvalues λk=k(k+1)/2,k≥1. Define the oblate spheroidal wave functions on [0,1] as (5.1)wkσ(y)=∑l≥1′fklvl(y),
where fkl obey the three-term recurrence defined in [24]. In the latter equality, the l summation is over odd (even) values if k is even (odd).Definevkσ(y)=eσywkσ(y) and ukσ(x)=(1/m(x))vkσ(x) where m(x)=e2σx/(x(1-x)) is the speed measure density of the WF model with selection (3.8).The system(ukσ(x),vkσ(x))k≥1 constitute a system of eigenfunctions for the WF with selection generators -G and -G* with eigenvalues λkσ implicitly defined in [24], thus with -G(ukσ(x))=λkσukσ(x) and -G*(vkσ(y))=λkσvkσ(y). The eigenfunction expansion of the transition probability density of the WF model with selection is thus, [25],(5.2)p(x;t,y)=∑k≥1bkσe-λkσtukσ(x)vkσ(y),
where bkσ=〈vkσ,ukσ〉-1. The WF model with selection can be viewed as a perturbation problem of the neutral WF model (see [3]). There exist perturbation developments of λkσ around λk with respect to σ2, [25]. They are valid and useful for small σ.The WF diffusion processxt with selection (3.8) is nonconservative, with finite hitting time τx of one of the boundaries. Following the general arguments developed in Section 2, the Yaglom limit of xt conditioned on τx>t is the normalized version of (5.3)v1σ(y)=eσyw1σ(y).
The limit law of xt conditioned on never hitting the boundaries in the remote future is the normalized version of (5.4)u1σ(y)v1σ(y)=1m(y)v1σ(x)2=y(1-y)w1σ(y)2.
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): dxt=σxt(1-xt)dt+xt(1-xt)dwt, x0=x∈(0,1). For this model, G=σx(1-x)∂x+(1/2)x(1-x)∂x2 and both boundaries are exit.Assume thatσ>0 so that the drift term is bounded above by f*=σ/4, together with 2f/g2 being bounded below (as a constant function here equal to 2σ). 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 α=e-ax, the α-Doob transform of xt could lead to a transformed process with bounded killing rate d=-G(α)/α. We shall choose a=σ for its interesting features.The transition densityp(x;t,y) of xt admits the representation (5.2) in terms of their oblate spheroidal wave eigenfunctions. Let (6.1)α(x)=E(∫0τxc(xs)ds),
where c(xs)=2xs(1-xs)e-σxs/4 is the skewed sample heterozygosity, damped by the factor e-σxs/4. Then, α solves -G(α)=(1/2)σ2x(1-x)e-σx, with solution α(x)=e-σx. In this case study, one selects sample paths of (xt;t≥0) with large α(y). The dynamics of (x̃t;t≥0) is the drift-less neutral WF dynamics dx̃t=x̃t(1-x̃t)dwt, subject to quadratic killing at rate d(x)=(1/2)σ2x(1-x) in I, which is bounded above there. The boundaries of (x̃t;t≥0) are still exit and the stopping time τ¯x of (x̃t;t≥0) is τ¯x=τ̃x∧τ̃x,∂, where τ̃x is its absorption time at the boundaries and τ̃x,∂ its killing time. The density of the transformed process is p¯(x;t,y)=(α(y)/(α(x)))p(x;t,y). Its series expansion is exactly known using (5.2) for p.The transformed process(x¯t;t≥0) accounts for a neutral evolution of the allele A1 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 Section2, we obtain the following.(i) Conditioned on τ¯x>t, the transformed process (x¯t;t≥0) admits a Yaglom limit q¯∞. With v1σ the first oblate spheroidal eigenvector of -G* associated to the smallest positive eigenvalue λ1σ, q¯∞ is of the product form (6.2)q¯∞(y)=e-σyv1σ(y)∫01e-σyv1σ(y)dy=w1σ(y)∫01w1σ(y)dy.
This limiting probability q¯∞ is the Yaglom limit law of (x¯t;t≥0) conditioned on the event τ¯x>t that both the absorption and killing times exceed t.(ii) Let G¯(·)=eσxG(e-σx·) be the infinitesimal generator of (x¯t;t≥0) with two stopping times. Now, there is a tradeoff between which of τ̃x and τ̃x,∂ occurs first. To solve it, we need to compute β defined in (2.55) by G¯(β(x))=0, with boundary conditions β(0)=β(1)=1. This is a Sturm-Liouville problem whose solution in our case is (6.3)β(x)=e-σx+e-σ(1-x)1+e-σ.
The function β(x)=P(τ̃x<τ̃x,∂) is minimal when x=1/2, with value β(1/2)=1/(cosh(σ/2)). This looks natural because when x̃0=x=1/2, the chance to hit {0,1} before getting killed should be the lowest.The density of the process(x¯t;t≥0) conditioned on the event {τ̃x<τ̃x,∂} is (6.4)p¯a(x;t,y)=β(y)β(x)p¯(x;t,y)=β(y)β(x)α(y)α(x)p(x;t,y),
and so is also explicitly known from the oblate spheroidal wave expansion (5.2) of p(x;t,y). The tail distribution of τ̃x given {τ̃x<τ̃x,∂} is obtained by integrating p¯a over y.Similarly, the density of the process(x¯t;t≥0) conditioned on the event {τ̃x,∂<τ̃x} is (6.5)p¯∂(x;t,y)=1-β(y)1-β(x)p¯(x;t,y)=1-β(y)1-β(x)α(y)α(x)p(x;t,y).
The tail distribution of τ̃x,∂ given {τ̃x,∂<τ̃x} is obtained by integrating p¯∂ over y.The associated conditioned on absorption first process(x̃a,t;t≥0) obeys the SDE (6.6)dx̃a=f̃a(x̃a)dt+g(x̃a)dwt,
with drift (6.7)f̃a(x)=f(x)+g2(x)[α′α(x)+β′β(x)]=σx(1-x)1-eσ(1-2x)1+eσ(1-2x),
and local variance unchanged g2(x)=x(1-x). This process has no killing part and it gets eventually absorbed at {0,1}.In the generatorG¯∂(·) of the conditioned on killing first process, there is a killing multiplicative part which is enhanced d/(1-β)>d and a shift in the drift, showing that the associated conditioned process (x̃∂,t;t≥0) exhibits a faster killing rate, but the drift shift guarantees that (x̃∂,t;t≥0) is not absorbed at the boundaries. With g2(x)=x(1-x) and β as in (6.3), the drift takes the peculiar explicit form (6.8)f̃∂(x)=-g2β′1-β(x).


## 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 Section2.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:dxt=xt(1-xt)dwt, x0=x∈(0,1). For this model, G=(1/2)x(1-x)∂x2 and both boundaries are exit. Its transition density p(x;t,y) admits the representation (7.1)p(x;t,y)=∑k≥1bke-λktuk(x)vk(y),
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 α(x)=eσx and consider G¯(·)=α-1G(α·). We now have G(α)=(1/2)σ2x(1-x)eσx and b(x)=G(α)/α>0.In this case study, one selects sample paths of(xt;t≥0) with large α(y). The dynamics of (x̃t;t≥0) is easily seen to be the WF with selection dynamics (7.2)dx̃t=σx̃t(1-x̃)dt+x̃t(1-x̃t)dwt,
subject to quadratic branching at rate b(x)=(1/2)σ2x(1-x) inside I. We indeed have(7.3)G¯(⋅)=e-σxG(eσx⋅)=b(x)⋅+G̃(⋅),
where G̃ is the KBE operator of the dynamics (x̃t;t≥0).Withβ(x)=α(x)-1=e-σx, we clearly have (7.4)G¯(β(x))=0,
and β is an harmonic function for G¯ and as a result, Doob-transforming G¯ by β, we get (7.5)β-1G¯(β⋅)=G(⋅),
which is the infinitesimal generator of the original neutral WF martingale.The birth (creating) rateb in G¯ is bounded from above on (0,1). It may be put into the canonical form b(x)=b*(μ(x)-1), where b*=maxx∈[0,1](b(x))=σ2/8>0 and (7.6)μ(x)=1+4x(1-x),
whose range is the interval [1,2] as x∈[0,1].The density of the transformed process isp¯(x;t,y)=(α(y)/(α(x)))p(x;t,y). It is exactly known because p is known from (7.1).The transformed process (with infinitesimal backward generatorG¯) accounts for a branching diffusion (BD), where a diffusing mother particle (with generator G̃ and started at x) lives a random exponential time with constant rate b*. When the mother particle dies, it gives birth to a spatially dependent random number M(x) of particles (with mean μ(x)). If M(x)≠0, M(x) 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μ(x) 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 (7.7)M(x)=0w.p.p0(x)=0,M(x)=1w.p.p1(x)=2-μ(x),M(x)=2w.p.p2(x)=μ(x)-1,
with p2(x)≥p1(x) (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{0,1} which are still absorbing for the system of particles and the number of particles Nt(x) in the system at each time t, which may grow due to branching events. The density p¯ of the transformed process has now the following interpretation: (7.8)p¯(x;t,y)=E[∑n=1Nt(x)p(n)(x;t,y)],
where p(n)(x;t,y) is the density at (t,y) of the nth alive particle descending from the ancestral one (Eve), started at x. In the latter formula, the sum vanishes if Nt(x)=0. A particle is alive at time t if it came to birth before t and has not been yet absorbed by the boundaries.Letρ¯t(x)=∫(0,1)p¯(x;t,y)dy. Then, ρ¯t(x) is the expected number of particle alive at time t. We have (7.9)∂tρ¯t(x)=G¯(ρ¯t(x)),ρ¯0(x)=1(x∈(0,1)).
But then, q¯(x;t,y):=p¯(x;t,y)/ρ¯t(x) obeys the forward PDE (7.10)∂tq¯(x;t,y)=(-∂tρ¯t(x)ρ¯t(x)+b(y))q¯(x;t,y)+G̃*(q¯(x;t,y)),
as a result of ∂tp¯(x;t,y)=G¯*(p¯(x;t,y)). We have (7.11)q¯(x;t,y)=E[∑n=1Nt(x)p(n)(x;t,y)]E[Nt(x)]
showing that q¯(x;t,y) is the average presence density at (t,y) of the system of particles all descending from Eve started at x.Clearly,-(logρ¯t(x)/t)→t→∞λ1=1 (and, therefore, also -(∂tρ¯t(x)/(ρ¯t(x))), because (7.12)ρ¯t(x)=1α(x)∑k≥1bke-λktuk(x)∫01α(y)vk(y)dy.
The expected number of particles in the system decays globally at rate λ1.The BD transformed process, therefore, admits an integrable Yaglom limitq¯∞, solution to -G̃*(q¯∞)=(λ1+b(y))q¯∞ or -G¯*(q¯∞)=λ1q¯∞. With v1(y)=1, the first eigenvector of -G* associated to the smallest positive eigenvalue λ1=1, q¯∞ is of the product form (7.13)q¯∞(y)=eσyv1(y)∫01eσyv1(y)dy=σeσyeσ-1.
This limiting probability q¯∞ is the Yaglom limiting average presence density at (t,y) for the BD system of particles (it is also the ground state for G¯*).There is also a natural eigenvectorϕ¯∞ of the backward operator -G¯, satisfying -G¯(ϕ¯∞)=λ1ϕ¯∞ (the ground state for G¯). It is explicitly here that (7.14)ϕ¯∞(x)=1α(x)u1(x)=e-σxx(1-x).
In the terminology of [26], both operators G¯(·)+λ1· and its adjoint are critical  (G¯(·)+λ1· (G¯*(·)+λ1·) is said to be critical if there exists some function ϕ¯∞∈C2 (q¯∞∈C2,resp.), strictly positive in (0,1), such that: G¯(ϕ¯∞)+λ1ϕ¯∞=0 (G¯*(q¯∞)+λ1q¯∞=0,resp.) and the operators do not possess a minimal positive Green function.). In this context, the constant λ1 is called the generalized principal eigenvalue. The eigenfunctions (ϕ¯∞,q¯∞) are their associated ground states.We note that we have theL1-product property (see [26, Subsection 4.9]). (7.15)∫01u1(x)v1(x)dx=∫01ϕ¯∞(x)q¯∞(x)dx<∞.
With pn(x)=P(M(x)=n), let (7.16)l(x)=∑n≥1pn(x)nlog⁡n=2log⁡2p2(x).
We have the xlogx condition (7.17)∫01l(x)ϕ¯∞(x)q¯∞(x)dx=8log⁡2∫01x(1-x)u1(x)v1(x)dx<∞.
We conclude (following [9, 10]) that, as a result of the condition (7.17) being trivially satisfied, global extinction holds in the following sense:(i)P(Nt(x)=0)→t→∞1, uniformly in x,(ii) there exists a constant γ>0:eλ1t[1-P(Nt(x)=0)]→t→∞γϕ¯∞(x), uniformly in x,(iii) For all bounded measurable function ψ on I,(7.18)E[∑n=1Nt(x)ψ(x̃t(n))∣Nt(x)>0]⟶t→∞γ-1∫(0,1)ψ(y)q¯∞(y)dy.
From (i), 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 μ(x)=EM(x)>1 for all x∈(0,1)): probability mass escapes out of I although the BD survives with positive probability.In the statement(ii), the quantity 1-P(Nt(x)=0)=P(Nt(x)>0) is also P(T(x)>t) where T(x) is the global extinction time of the particle system descending from an Eve particle started at x. The number -λ1 is the usual Malthus decay rate parameter. From  (ii),ϕ¯∞(x) 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.(iii) with ψ=1 reads E[Nt(x)|Nt(x)>0]→t→∞γ-1 giving an interpretation of the constant γ (which may be hard to evaluate in practise).The ground states ofG¯+λ1 and its adjoint are, thus, (ϕ¯∞,q¯∞) and explicit here. It is useful to consider the process whose infinitesimal generator is given by the Doob transform (7.19)ϕ¯∞-1(G¯+λ1)(ϕ¯∞⋅)=ϕ¯∞-1(G̃+b+λ1)(ϕ¯∞⋅),
because product-criticality is preserved under this transformation. The ground states associated to this new operator and its dual are (1,ϕ¯∞q¯∞). Developing, we obtain a process whose infinitesimal generator is (7.20)G̃+ϕ¯∞′ϕ¯∞g2∂x=G+u1′u1g2∂x,
with no multiplicative part. In our case study, we get (1/2)x(1-x)∂x2+(1-2x)∂x adding a stabilizing drift towards 1/2 to the original neutral WF model. The associated diffusion process is positive recurrent and so its invariant measure ϕ¯∞q¯∞=u1v1∝y(1-y) is integrable. It is the beta(2,2) limit law of the Q-process (see (2.80) and (4.23)) for the neutral WF diffusion.Remark 7.1.
At timet, let (x̃t(n))n=1Nt(x) denote the positions of the BD particle system. Let u(x,t;z)=E[∏n=1Nt(x)zψ(x̃t(n))] stand for the functional generating function (|z|≤1) of the measure-valued branching particle system. u(x,t;z) obeys the nonlinear (quadratic) Kolmogorov-Petrovsky-Piscounoff PDE, [27] (7.21)∂tu(x,t;z)=b*θ(x,u(x,t;z))+G̃(u(x,t;z));u(x,0;z)=zψ(x),
where θ(x,z)=E[zM(x)]-z=(p2(x)z2+p1(x)z)-z or(7.22)θ(x,z)=4x(1-x)z(z-1)
is the shifted probability generating function of the branching law of M(x). Thus, the nonlinear part reads b*θ(x,u(x,t;z))=b(x)u(x,t;z)(u(x,t;z)-1), which is quadratic in u.
In particular, ifu(x,t):=∂zu(x,t;z)z=1=E[∑n=1Nt(x)ψ(x̃t(n))], u(x,t) obeys the linear backward PDE(7.23)∂tu(x,t)=b(x)u(x,t)+G̃(u(x,t));u(x,0)=ψ(x),
involving G¯(·)=G̃(·)+b(x). The latter evolution equation is the backward version of the forward PDE giving the evolution of p¯(x;t,y) as ∂tp¯(x;t,y)=G¯*(p¯(x;t,y)),p¯(x;0,y)=δx(y).



---

*Source: 605068-2011-07-18.xml*