The speed of invasion in an advancing population

We derive rigorous estimates on the speed of invasion of an advantageous trait in a spatially advancing population in the context of a system of one-dimensional F-KPP equations. The model was introduced and studied heuristically and numerically in a paper by Venegas-Ortiz et al. (Genetics 196:497–507, 2014). In that paper, it was noted that the speed of invasion by the mutant trait is faster when the resident population is expanding in space compared to the speed when the resident population is already present everywhere. We use the Feynman–Kac representation to provide rigorous estimates that confirm these predictions.


INTRODUCTION
The present paper is motivated by an interesting paper by Venegas-Ortiz, Allen, and Evans [22] that investigates the invasion of a spatially expanding population by a new trait.The classical model for the invasion of a gene in a spatially extended population [9] or the expansion of a population in space [16] is the Fisher-Kolmogorov-Petrovsky-Piscounov (F-KPP) equation, that has been the subject of intense investigation for over 80 years 1 .The F-KPP equation is a non-linear reaction-diffusion equation that admits travelling wave solutions to which solutions starting with suitable initial conditions converge.This has been known since the early work of Kolmogorov et al. [16], but has been made both more precise and more general in the seminal book by Bramson [5].
The model discussed in [22] is a system of two coupled equations of the F-KPP type that describes the evolution of a population of two types (traits, alleles, ..) that diffuse, compete, and and switch between types.More specifically, they propose the system of equations N A , N B represent the masses of traits A and B, K is the carrying capacity, α, β, γ are parameters that satisfy α > γ > β/K ≥ 0. (1. 3) The different terms in these equations correspond to the following biological mechanisms: (i) The terms ∂ xx N model the spatial diffusion of the population.Note that the diffusion coefficients are the same for both types.This can be seen as biologically plausible, but this choice is mainly done to simplify the mathematical treatment.(ii) The terms proportional to α describe logistic growth with the quadratic terms corresponding to competitive pressure.Again it is assumed that the pressure exerted by both types and on each type are the same.This again simplifies the mathematics.(iii) The linear terms ±βN A can be interpreted as mutation rates from the A population to the B population.There effect is a net disadvantage of the A population.(iv) The non-linear terms ±γN A N B are interpreted as horizontal gene transfer from the B-types to the A-types.The idea is that when an A individual encounters a B individual, the genotype of the B individual can be switched to the A-type.
The choice of parameters in (1.3) ensures that the a priori disadvantaged A type can reemerge in a developed B-population and a stable equilibrium with co-existing types exists.The question addressed in [22] is to analyse how this effect leads to a hitchhiking of the A-type when the B-type is spreading in space.The authors of [22] make the following interesting and somewhat surprising observation.There are two easily derived travelling waves in the system.First, a population made purely of B individuals remains in that state and advances with a speed v B .Second, if B has invaded all space, and a A population is introduced, there is (with the choice of parameters that ensures the instability of the B population against the invasion of A individuals) a travelling wave of A particles that advances in the background of B particles with a speed v A < v B .If, however, one starts with initial conditions where A and B particles are present, say in the negative half-line, then the B population advances with speed v B again, but in some parameter range the A population advances with a speed v c that is strictly larger than the speed v A (and smaller than v B ). Somehow, the A individuals sense the empty space ahead of the B-wave and get attracted to it.Venegas-Ortiz et al [22] derive this result, and precise formulas for the speeds, using local linearisation and matching of solutions.These findings are supported by numerical simulations.In the present paper we derive rigorous estimates on the speeds using the Feynman-Kac representation, originally employed by Bramson [5] to control the precise speed of convergence to the travelling wave in the original F-KPP equation.It turns out that this point of view not only allows to give rigorous and precise bounds on the solutions of the system of equations, and hence the speeds, but also provides a clear and intuitive explanation for the fact that the empty space ahead of the B-wave allows for a faster advance of the A-wave.Namely, we will see that this is driven by large, unlikely excursion of the Brownian motion in the Feynman-Kac formula that reach ahead of the front of the B-wave.Mathematically, this involves some delicate estimates on probabilities of large excursions of Brownian bridges.
Systems of coupled F-KPP equations have been studied in different contexts in the literature, see e.g.[10,6,11,13,14,12,7,15,18].In particular, an analogous result to that in [22] and the present paper was derived rigorously in [14] using analytic methods.Rather recently, there has been interest in such systems in the context of dormancy, see e.g.[4].Applicable tools depend on the details of the equations.[10] use purely analytic methods involving sub-and super-solutions, while the equations appearing in [6] and [4] allow for a representation in terms of branching Brownian motion and the use of martingale methods.The equations in [22] (and [14,12,7]) are particularly nice, as they allow for the use of the Feynman-Kac representation.However, even the introduction of two different diffusion constants seems to spoil this feature, and it seems unclear (albeit interesting) to see how this method can be extended to more general settings.Outline.The remainder of this paper is organised as follows.In Section 2 we give a precise formulation of the model put forward in [22] and explain the special structure of the system that effectively reduces the problem to a time-dependent one-dimensional F-KPP equation.Afterwards we state our main result.Along the way we also recall some background on the standard F-KPP equation that will be needed.In Section 3 we present the Feynman-Kac representation, derive some first bounds, and give a heuristic explanation of the main result, based on the Feynman-Kac representation.Section 4 provides the necessary upper and lower bounds on the excursions of Brownian bridges.We compute fairly sharp bounds on the Laplace transforms of these excursions using the Laplace method.Armed with these estimates, we derive upper and lower bounds on solutions from which the wave speed v c is inferred in Section 5.At the end of the paper, in Section 6, we discuss our results and point to possible future extensions.

THE F-KPP EQUATIONS
It is convenient to introduce the total population mass N T ≡ N A + N B and to write the equations (1.1) and (1.2) in the the form We see that N T satisfies an autonomous F-KPP equation.Effectively, the second equation is a F-KPP equation with time dependent reaction rates.This structure is crucial for our analysis building on the Feynman-Kac formula.Equations of a similar structure have been also been studied in [12,7].It is furthermore convenient to eliminate the parameters K and α by rescaling.We define Then v and w solve where β = β/(αK) and γ = γ/α.Note that 1 > γ > β > 0. Note that the system of equations has four spatially constant fixpoints: The fixpoint (ii) is unphysical, since it corresponds to a negative mass for the population B. The fixpoints (i) and (iii) are unstable, and (iv) is the stable fixpoint.
The behaviour of v is well-known from Bramson's work [5], so solving for w amounts to solve the F-KPP equation with time dependent coefficients.A particularly simple situation arises if we choose initial conditions such that v(0, x) = 1, for all x ∈ R. In that case w solves the F-KPP equation In this case, with suitable initial conditions (e.g.Heaviside), w converges to a travelling wave solution that moves with speed 2(γ − β).A more interesting situation arises if the initial conditions are such that v(0, x) decays rapidly at +∞ and w(0, x) is non-zero.In that case, [22] observed that the w-wave follows behind the v-wave, but moves faster than it would in a fully established population.Recall that the standard F-KPP equation (2.5) admits travelling wave solutions where ω solves the ode 1 2 for all speeds ≥ √ 2. It was shown by Kolmogorov [16] that (2.9) has a unique solution up to translations such that lim x↓−∞ ω(x) = 1 and lim x↑∞ ω(x) = 0. We are only interested in the case λ = √ 2, since solutions with initial condition that converge rapidly to zero at infinity, and in particular with Heaviside initial conditions, converge to travelling waves with this speed (see [5] for more details).
Then for all δ > 0 sufficiently small there exist constants C Remark.Note that u c is strictly larger than 2 γ − β for β small enough.Notice that Venegas-Ortiz et al. derive in [22] a rather complicated looking equation, (Eq.8), and a simpler one (Eq.9), obtained by expanding in β.Our results show that the second version is exact, provided β is such that while the first seems incorrect.This is also in agreement with the finding in [10].An analogous result on an accelerated speed in a slightly different system of equations was derived by purely analytic methods by Holzer and Scheel [14], Lemma 11.
Remark.Note that in fact the result of Theorem 2.1 does not depend on the choice of a in the initial condition.This is not surprising as a finite shift of the initial condition does not affect the large time asymptotic of the solutions.
The remainder of this paper is devoted to proving Theorem 2.1.In the process, we will derive precise bounds on the behaviour of the solutions.

THE FEYNMAN-KAC REPRESENTATION
Bramson's analysis of the F-KPP equation [5] is based on the Feynman-Kac representation.We will do the same for the equation (2.16).

The representation and elementary bounds.
Lemma 3.1.The solution of (2.16) satisfies the equation ) where B is a Brownian motion starting in x.
Proof.The proof is identical to the one in [5].□ It is convenient to express the Brownian motion B in terms of its endpoint B t and a Brownian bridge from x to B t .Here z t 0,0 is a Brownian bridge from 0 to 0 in time t.Note that the bridge is independent of B t .This leads to the following reformulation of (3.1).
Lemma 3.2.The solution of (2.16) satisfies where E now refers to the expectation with respect to the Brownian bridges z t x,y resp.z t 0,0 .Proof.Elementary.□ The fact that 0 ≤ ω ≤ 1 and 0 ≤ w ≤ 1 − γ/ β yields the first bounds.
Lemma 3.3.The solution of (2.16) satisfies and For Heaviside initial conditions, this implies Proof.Eqs.(3.4) and (3.5) are immediate from the bounds on ω and w mentioned above.
Since w ≤ ω, we also have the lower bound To see how we can use these bounds, let us first ignore the possible excursions of the Brownian bridge and simply set z t 0,0 (s) = 0. We want to see where w(t, x) drops from 1 to zero.From (3.6) we already know that this must happen before x = 2(1 − β)t.Now assume that for some u ≤ 2(1 − β), w(t, ut + z) ≤ ϵ, for all z ≥ 0.Then, for z ≥ 0 independent of t, which tends to infinity if u < 2(γ − β).Hence, the hypothesis can only be true for u ≥ 2(γ − β).On the other hand, if , we get the corresponding upper bound which is decaying exponentially with t.This suggests a wave moving at speed u 0 = 2(γ − β), which is the speed we obtain if v(0, x) ≡ 1.This shows that the only way to move faster is to exploit the possibility of the Brownian bridge to make a forward excursion out of the region where ω = 1.

3.3.
Improved heuristics on the wave speed.First, note that in (3.3) y is negative, so that we cannot gain anything from it and pretend that it is equal to zero in this subsection.
To simplify the heuristics we also set a = 0.Moreover, as we are analysing the possible gain in ω by large Brownian bridge excursions to areas where ω is small, we will ignore w (which is always way smaller than ω) in (3.3).Hence, we are left with estimating For our heuristics we approximate ω by Hence, to further estimate the expectation in (3.11) we need an estimate on the time during which the indicator function takes the value 0. To this end, let with α = √ 2 − x/t, be the time the Brownian bridge spends above a line with slope α.Note that (3.11) Next, on the exponential scale where we used that heuristically the cheapest way to realise the event {T t > S } is to stay above this line up to roughly time S .This probability is roughly dominated by the event to be essentially on the line at time S .As we gain a factor (1 − γ) (on the exponential scale) as long as the Brownian bridge is above the line with slope ( √ 2 − x/t), to find the dominating event in the expectation in (3.11) we need to find the optimal S * , namely By differentiating the right-hand side of (3.16), we see that . (3.17) Now, we distinguish two cases.(Case 1) If S * is positive, we plug this back into (3.14).Then the exponent in (3.14) is to leading order equal to To see where w starts to decay to 0, we need to see for which x (3.18) is equal to zero (hence its exponential is of order 1).This leads to The exponent is zero if as a function of β, we observe that it is decreasing in β and there is exactly one critical value β * 1 such that Similarly, seeing x * 2 ( β) as a function if β we observe that it is decreasing in β and there is exactly one critical value β * 2 such that As the two critical values for β are the same, this suggests that for β > β * 1 the speed of the wave equals x * 2 /t and increases continuously to x * 1 for β < β * 1 .This will be made rigorous in the following sections.

BROWNIAN BRIDGE ESTIMATES
In this section we provide the key input about Brownian bridges that is needed to make the heuristics above rigorous.
4.1.Probabilities of excursions.As ω is not exactly an indicator function, the key question is to know the distribution of the time a Brownian bridge z t 0,0 spends well above and well below a line ( Note that z t 0,0 (s) has the same law as z t 0,0 (t − s), and so we can replace T K t by for convenience.The following theorem provides precise tail asymptotic for T K t .Theorem 4.1.Let z t 0,0 be a Brownian bridge from zero to zero in time t.Let α > 0 and T t defined in (4.1).Then, for 0 < s ≤ 1, if K > 0.
Proof.To start, we define g t as the last time the Brownian bridge z t 0,0 is above the line αs + K, (see Figure 4.1) Then FIGURE 2. Schematic picture of the Brownian bridge in reversed time The conditional probability in (4.6) is known [21,20].A more convenient formula is given in [2], see Eq. ( 7) therein.For our setting this yields where Φ is the error function.Note that for K = 0, this simplifies to which recovers the result that the time spent by a Brownian bridge from 0 to 0 in time g t above 0 is uniformly distributed on [0, g t ].
Next we need to control the distribution of g t .Fortunately, this can be recovered from known results by Beghin and Orsingher [3].Lemma 4.2.With the notation above, (4.9) Proof.Looking back in time, we see that we can also interpret g t as By time reversal, this has the same law as t − h t where The latter probability can be computed using a result by Beghin and Orsingher [3] (Lemma 2.1).It yields that, for αt + K > 0, . (4.12) If α(t − r) + K ≤ 0, then this probability is equal to one.Note that, in particular, Note that the term in the second line in (4.12) is (asymptotically equal) and smaller than If αt(t − r) + K(t − 2r) > 0, the first term in (4.12) is asymptotically equal to and smaller than e − 2K(αt+K) Recalling that g t = t − h t , we get that P (g t ≥ q) = P (h t ≤ t − q) and hence the assertion of the lemma follows. □ We compute the probability density of the distribution of g t by differentiating (4.9).This gives the nice formula Thus, by (4.6), where we used (4.7) together with (4.16).We use the Laplace method to compute the integral in (4.17).The exponential term takes its maximum at v = s.Thus we need to compute the behaviour of the prefactor at s. Let us first consider the more complicated case K > 0. We get as x ↓ 0. Hence, Similarly, Inserting these asymptotics into (4.17),we find that, up to errors of order 1/t, so that finally In the remaining cases we get Therefore, using Lemma A.2, The control of the distribution of T K t given by Theorem 4.1 suffice to prove upper bounds on w and hence upper bounds on the wave speed.To prove lower bounds, it is also necessary to take possible fluctuations of the Brownian bridges in the negative direction into account.Therefore, we need on the distribution of T K t a lower bound where excursions of the Brownian bridge below zero are suppressed.We define, for b > 0, We want a lower bound on 27) The following lemma is not optimal but sufficient for our purposes.Lemma 4.3.For K > 0, b > 0, and L > 0,

.28)
Proof.Given g t , we use that The second event in turn contains the event that {S > g t − L}.
Hence, the main effort is to control the law of g t under the restriction that the bridge remains above −b.By the same reasoning as before, this amounts to proving a lower bound on To bound this, we have to revisit and alter the proof in [3].First, we note that The latter probability can be written up to normalisation as where B is a Brownian motion started in zero.Decomposing this over the values of B(r) gives Hence, Passing back to the Brownian bridge, this yields ) is monotone increasing in g t , it holds that Using (4.18), for L finite and S = st, for some C > 0. □ Remark.Note that, up to constants, the difference between the expression for P(T K t > st) is that a factor 1/s is missing; this is due to the lower bound in (4.39).To keep the difference in upper and lower bound of polynomial order in t one could choose L ∼ K and Lemma 4.4.Assume that 2λ > α 2 .Then, as t ↑ ∞, Proof.Note first that for any non-negative random variable T , From Theorem 4.1 we see that, for 0 < r < t, where P K (t, s) is polynomially bounded.Moreover, P (T ≥ r) = 1 for all r ≤ 0 and , s ∈ (0, 1).(4.46) f (s) takes its maximum in (0, 1) at provided that 2λ > α 2 .By an elementary computation, and the second derivative of f at s * is given by where and if K > 0. The claim of Lemma 4.4 follows.□ For the lower bound on w, we need to take negative excursions into account, as already mentioned in Section 4.1.The following Lemma provides a corresponding bound on the Laplace transform.Lemma 4.5.Assume that 2λ > α 2 and let K > 0.Then, as t ↑ ∞, Proof.The proof is a rerun of the proof of Lemma 4.4 using Lemma 4.3 instead of Theorem 4.1.□ Finally, we need bounds on the Laplace transform when 2λ ≤ α 2 .The following lemma confirms that in this case the Laplace transform is essentially of order one.Lemma 4.6.If 2λ ≤ α 2 and K < 0,  This finishes the proof of Lemma 4.6.□

CONTROLLING THE WAVE
We use the Brownian bridge estimates from the previous section to give a rigorous version of the heuristics outlined at the end of Section 3.
5.1.Bounds on the speed of the wave.We first control the behaviour of solutions on the exponential scale for large t.We begin with an upper bound.(i) Let u be such that Then, for all ϵ > 0 small enough, there exists a constant C > 0 such that In particular, w(t, ut) decays exponentially fast in t for (ii) Let u be such that Then, for all ϵ > 0 small enough, there exists a constant C > 0 such that (5.5) Remark.Lemma 5.1 implies that the solution is exponentially small if u > u c (given in (2.17)), hence the wave speed is not larger than u c .
Proof.We bound the integral in the Feynman-Kac representation (3.3) from above as follows.
Using the asymptotic of the lower tail (2.13), we see that on the second indicator function, for all y ≤ 0, which is larger than 1 − ϵ for −K large enough.On the first indicator function we use that ω ≥ 0. This leads to (5.8) , we obtain the upper bound (5.9) The exponential terms are (for x = ut), (5.12) This implies that w decays exponentially fast for u > 2(γ − β).□ Next, we need a corresponding lower bound.For this we use the lower bound from Lemma 4.5.
Lemma 5.2.Let b > 0. Assume that u is such that w(t, ut + z) ≤ ϵ, for all z ≥ −2b and t large enough.

.13)
This contradicts the hypothesis, unless , then there exists a constant C > 0 depending on K, b and u w(t, ut) ≥ C t 3/2 e (γ(1−ϵ)− β)t− u 2 t 2 (5.15)This contradicts the hypothesis unless u > 2 γ(1 − ϵ) − β . (5.16) Remark.From (5.14) it follows that the speed is not smaller than (5.16) it follows that it is not smaller than 2 γ − β .Altogether, this implies that the speed is not smaller than the maximum of the two, i.e. it is at least Proof.In the representation (3.3) of w(t, x) we would like to use the assumption of the lemma to argue that the term involving w in the exponent is negligible but this could be spoiled by large negative excursions of the Brownian bridge.To avoid this problem we restrict the expectation in (3.3) on the Brownian bridge to a subset {U b t ≤ L}, For any L > 0 In particular, w(t, u * t + z) ≤ C/t if, for some δ > 0, (5.30) Proof.The proof is straightforward from (5.9), choosing K = −c + ln t. □ The next lemma shows that this upper bound is not too bad.
Lemma 5.4.Assume that β and γ are such that u * > 2 γ − β .Let a be a constant independent of t.Then, with c − = √ 2, there exists a constant 0 < C < ∞, independent of t, such that, for x = u * t + z, w(t, u * t + z) ≤ C/t, then In particular, this is in contradiction with the assumption if Proof.This is straightforward from (5.21), choosing K = c − ln t. □ Lemma 5.4 tells us that w(t, u * t + z − ) is greater than O(1/t).The next lemma tells us that at time of order ln t later, this will have grown to O(1).Lemma 5.5.Let ϵ(t) > 0. Assume that w(t, ut + z) ≥ ϵ(t) ∀z ≤ 0. (5.33) Then, for all δ > 0 sufficiently small, there exists a constant c such that (5.34) Proof.Starting from the Feynman-Kac formula, we have In this paper we have used the Feynman-Kac representation to derive the speed of advance of a hitch-hiking subpopulation within an advancing population.Apart from the fact that this allowed fairly sharp control of the precise behaviour of the wave fronts, the method provides a very clear intuitive understanding of the reason for the acceleration in an advancing population compared to a fully established one.Namely, the acceleration is driven by rare excursion of a Brownian bridge reaching ahead of the B-population.Translating this back into an underlying individual based model, heuristically this may be interpreted as having excursions of A particles into the empty space ahead of the bulk wave taking advantage of higher growth rate in the absence of competition.
Technically, we took advantage of the special features of the model that allowed to reduce the analysis to that of a scalar F-KPP equation with time-dependent parameters.This is a delicate property that gets spoiled already if the diffusion coefficients of the two types are different.An explicit useable Feynman-Kac representation for systems of pdes does not exist.Still, we are optimistic that the Feynman-Kac representation (used for each one-dimensional component of the system) can be used in such situations.This is subject of ongoing research.

APPENDIX A. THE LAPLACE METHOD WITH PREFACTOR
To evaluate the asymptotic behaviour of integrals, we use Laplace's method.Below we state for convenience the results we use.Proofs can be found in many places in the literature, e.g.[8].
Lemma A.1.Let f : (0, 1) → R be a twice differentiable function with a unique maximiser at s * and f ′′ (s * ) < 0.Moreover, let P(s, t) be a rational function in s and t.Then We also need a similar statement for integrals that get their contribution from the boundary of the domain of integration.Lemma A.2. Let f : (y, 1) → R be a twice differentiable function, which is monotone decreasing.Moreover, let P(s, t) be a polynomially bounded function in s and t in [y, 1].If, for some δ ≥ 0, lim x↓y (y − x) −δ P(y, t) = c(t, y), then lim t↑∞ (t f ′ (y)) 1+δ e t f (y) 1 y P(s, t) c(t, y) e −t f (s) ds = Γ(1 + δ).(A.2)

FIGURE 1 .
FIGURE 1. Schematic picture of the Brownian bridge spending time T K t

2 .
The Laplace transforms.As seen in(3.14),we need to control the Laplace transform of T K t .The behaviour of the Laplace transform is very different weither 2λ > α 2 or 2λ ≤ α 2 .