Invasions Slow Down or Collapse in the Presence of Reactive Boundaries

Motivated by the propagation of thin bacterial films around planar obstacles, this paper considers the dynamics of travelling wave solutions to the Fisher–KPP equation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_t = u(1-u) + u_{xx} + u_{yy}$$\end{document}ut=u(1-u)+uxx+uyy in a planar strip \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\infty< x < \infty $$\end{document}-∞<x<∞, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 \le y \le L$$\end{document}0≤y≤L. We examine the propagation of fronts in the presence of a mixed boundary condition (also referred to as a ‘partially absorbing’ or ‘reactive’ boundary) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_y = \alpha u$$\end{document}uy=αu, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha >0$$\end{document}α>0, at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=0$$\end{document}y=0. The presence of boundary conditions of this kind leads to the development of front solutions that propagate in x but contain transverse structure in y. Motivated by the observation that the speed of propagation in the Fisher–KPP equation is determined (for exponentially decaying initial conditions) by the behaviour at the leading edge, we analyse the linearised Fisher–KPP equation in order to estimate the speed of the stable travelling front, a function of the width L and the imposed boundary conditions. For wide strips the speed estimate based on the linearised equation agrees well with the results of numerical simulations. For narrow channels numerical simulations indicate that the stable front propagates more slowly, and for sufficiently small L or sufficiently large \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α the front speed falls to zero and the front collapses. The reason for the collapse is the non-existence, far behind the front, of a stable positive equilibrium solution u(x, y). While existence of these equilibrium states can be demonstrated via phase plane arguments, the investigation of stability is similar to calculations of critical patch sizes carried out in similar ecological models.


Introduction
Many biological situations can be viewed as invasions of one population into another; either a new genetic trait or species out-competing an existing one for resources, or consuming a preexisting resource distributed over a spatial domain. For such situations a number of simple mathematical models have been proposed in order to quantitatively capture spatial aspects of the population dynamics and the propagation of the new species (Volpert et al. 1994;Shigesada and Kawasaki 1997;Murray 2002;Volpert and Petrovskii 2009;Lewis et al. 2016). The simplest situation considers the density of the invading population to be described by a single scalar variable u(x, t) that is a function of space and time. Further common simplifying assumptions are that the background state u = 0 is homogeneous, and locally unstable, i.e. that a small positive density of invaders u will grow over time. It is also plausible from biological considerations to argue that solutions should always remain bounded and hence the density u will saturate to a constant value at long times; without loss of generality we take this to be the state u = 1.
The simplest model of this kind describing propagation into a state where u is initially close to zero and then rises to saturate at u = 1 is the well-known Fisher-Kolmogorov-Petrovskii-Piskunov (FKPP) equation (Fisher 1937;Kolmogorov et al. 1937) which, in two spatial dimensions, takes the form where u = u(x, t) and x ≡ (x, y) ∈ L := (−∞, ∞) × [0, L] ⊆ R 2 . Our interest in this paper is in solutions in the form of 'fronts', i.e. solutions that for every fixed y and t are monotonically decreasing in x and which satisfy u(−∞, y, t) = 1 and u(∞, y, t) = 0. Instead of the domain L we occasionally consider the half-space := (−∞, ∞) × [0, ∞). Fronts therefore propagate from left to right along the x-axis. Travelling wave solutions are those for which u takes the form u (x, y, t) = u(x − ct, y) where c is the wave speed.
The corresponding problem in one spatial dimension, looking for y-independent solutions to (1) is discussed in detail by Volpert et al. (1994). For the case considered in this paper (the monostable case in which the form of the nonlinear term implies that the state u = 0 is unstable and u = 1 is stable), phase plane analysis shows that travelling waves cannot exist for wave speeds c < 2, but do exist for all speeds c ≥ 2 (Harris 1999).
Only at isolated specific values of the wave speed is it possible to find an exact analytic travelling wave solution to (1), one example is the exact solution u(x, t) = 1 2 tanh x − 5t/ √ 6 2 √ 6 + 1 2 2 , deduced by Ablowitz and Zeppetella (1979) which travels at a speed 5/ √ 6 ≈ 2.041 and which is an exact solution to (1) when u is taken to be independent of y. Murray (2002), citing earlier work by Canosa (1973), outlines an asymptotic approach to approximation of the form of travelling waves, taking the reciprocal of the wave speed 1/c to be a small parameter.
Stability issues have also been investigated in detail by many authors. In one spatial dimension it is well established that the solution to the initial value problem with exponentially decaying initial conditions (for example compactly supported initial conditions) will converge to a travelling wave solution that propagates with the minimum allowable speed c = 2. This was proved by Larson (1978) building on work of McKean (1975), and numerically verified by Manoranjan and Mitchell (1983). Initial conditions that decay more slowly converge to travelling waves that propagate faster, see for example Hamel and Roques (2010), Sherratt (1998) and the references therein.
Our interest in this paper is to consider extensions of the 1D Fisher-KPP problem to two-dimensional spatial domains. Clearly one-dimensional solutions are also exact solutions of the 2D Eq. (1) when it is posed either in an infinite spatial domain (x, y) ∈ R 2 , or in the domains or L defined above, with a Neumann boundary condition n · ∇u = 0 imposed everywhere.
In the case in which the domain is R 2 , families of exotic travelling wave solutions have been studied through analysis of the linearised problem by Brazhnik and Tyson (2000), see also Showalter (1995). In the case of the strip L , the evolution of fronts for the Fisher-KPP equation becomes connected to ecological questions of the critical patch size for the maintenance of the population since a necessary condition for the propagation of a travelling wave solution is that the population should be able to maintain itself at an equilibrium density in the region behind the travelling wave. This is an important question in the case we study here since reproduction and dispersal in a finite domain compete with extinction, or at least a reduced reproduction rate, outside the finite domain.
The effect of a finite domain size on population dynamical models in general has been studied by many authors, stretching back to Skellam (1951). It is well known that, for simple ecological models such as the 1D version of (1), say for u(y, t) independent of x, the imposition of Dirichlet boundary conditions u(0, t) = u(L , t) = 0 implies that a population cannot persist if the domain length L is less than a critical length L c := π (Kierstead and Slobodkin 1953;Ludwig et al. 1979).
In more general cases in which the environment for the population growth and dispersal changes unfavourably along the boundaries, rather than taking Dirichlet boundary conditions u = 0 at y = 0 and y = L, more realistic boundary conditions are the mixed-type (or 'reactive', or 'partially absorbing') boundary conditions u y = αu at y = 0 and u y = −βu at y = L, where α, β ≥ 0. This corresponds to situations where the Fisher-KPP equation is motivated as a continuum description arising from a microscopic model in which, when an individual crosses the boundary there is a positive probability of absorption onto the boundary rather than reflection back into the interior of the domain (Erban and Chapman 2007). Such partially absorbing behaviour has been studied experimentally observing the propagation of bacteriophages through populations of E. coli. Experimental observations indicated that the propagation fronts of the invading bacteriophage population developed a curved profile near the boundaries of the available bacterial population. This phenomenon could not be reproduced by assuming a no-flux boundary condition (Möbius et al. 2015), which provided the original motivation for the consideration of the partially absorbing boundary conditions that we use here.
The structure of this paper is as follows. In Sect. 2 we examine the linearised problem at the leading edge of a front in the Fisher-KPP equation. We solve the linear problem and hence deduce the speed c of the front in terms of the parameters that determine the reactive boundary condition. As a result we are able to estimate analytically how the presence of the reactive boundaries cause the travelling fronts to slow down in wide strips.
In Sect. 3 we discuss the form of the fully nonlinear travelling front far behind the leading edge. This also turns out to be analytically tractable. It also explains the shape of contours of constant population density when plotted in the (x, y) plane to describe or detect the travelling front, for example experimentally.
As the strip becomes narrower, the population density behind a stable propagating front is observed numerically to decrease. This motivates, in Sect. 3.2 an investigation of whether there is a critical width below which the strip is unable to sustain a positive population behind the front. This is indeed the case. In sufficiently narrow strips, when α and β are sufficiently large, the population density falls to zero behind the front and the solution collapses to zero everywhere: there appears to be no stable travelling wave in this regime.
Section 4 summarises these results, compares them with previous results for critical patch sizes, and outlines directions for future work.

The Speed of the Front
The FKPP equation is a basic model for the propagation of a stable state u = 1 into an unstable state u = 0. The speed of propagation is expected to be controlled by the dynamics at the leading edge of the front since this is a pulled front (Saarloos 2003). We consider the FKPP equation posed in a strip (x, y) ∈ R × [0, L]. We solve the initial value problem always considering compactly supported initial conditions. We therefore expect (but have not proved) that solutions exist and are unique, and converge to the unique solution that travels with the slowest available speed.
Our interest is in the transverse structure (i.e. the variation of u(x, y, t) in the y direction) at the leading edge of the front. If the boundary conditions at y = 0 and y = L are zero flux (Neumann) conditions u y = 0 then we would expect that there is no transverse variation and the problem is one-dimensional: solutions converge to travelling wavefronts that are independent of y. It is not a priori clear that this should be the case with Neumann boundary conditions, but numerical solutions provide evidence that this is indeed the case.
Imposing reactive (or mixed, or Robin) boundary conditions of the form u y = αu at y = 0, and u y = −βu at y = L , does, however, generate transverse structure as we now show. In the region x 0 the amplitude u is small so that the FKPP equation can be well approximated by its linearisation u t = u x x + u yy + u. (3) We substitute the separable ansatz u = e −(x−ct) v(y) suitable for a constant coefficient linear PDE, and assuming that the leading edge of a stable front will decay exponentially in the x direction with exponent −1 as in the case in one dimension. This decay rate of the leading edge is confirmed by numerical investigations, as shown for example in Fig. 1b.
With this assumption, the speed c of propagation of the front is determined from the resulting eigenvalue equation for v(y): together with boundary conditions of the form (2). We now consider three cases in turn; in every case we impose the boundary condition v y = αv at y = 0, but we consider in turn imposing either a no-flux boundary condition u y = 0 at y = L, or a reactive boundary condition u y = −βu at y = L, or in the third case we consider (4) to be posed in the half-space y ≥ 0 with only a single boundary condition required (at y = 0).

One Partially Absorbing Boundary and One No-Flux Boundary
For the case in which the boundary at y = 0 is partially absorbing: v y = αv, and for which v satisfies a no-flux boundary condition at y = L, i.e. v y = 0 at y = L, we would intuitively expect that the speed c is reduced below the value 2 which would correspond to the 1D FKPP front speed (i.e. the case α = 0. Given the form of (4) it is convenient to introduce the new variable p := √ 2 − c). Then the solution for v takes the form where A is an undetermined multiplicative constant, and p is the smallest positive solution of the equation Although (6) has a countable collection of solutions for p (eigenvalues), the physically relevant one is the leading eigenvalue (the smallest positive value of p) since this leads to the linear eigenmode with the largest positive growth rate and therefore this mode will dominate the solution for the leading edge of the front. We note that the corresponding eigenfunction v(y) is monotonically increasing, which agrees with our intuition about the biological interpretation of the reactive boundary condition being to reduce the population density in the vicinity of the boundary. The form (5) provides a very good approximation to the shape of the leading edge of the front, as shown in Fig. 1. Figure 1a shows the transverse (i.e. y-direction) structure  Fig. 1a is the theoretical form of the leading edge of the front; the dash-dotted line is very close to it since at this moment in time the leading edge of the front is at the position x = 90 at which the (rescaled) profile u(90, y, t)/u(90, L , t) is plotted. The dotted line is a slice through the domain at x = 80 which is through the middle of the front, departing from the profile computed from the leading edge. The dashed line, at x = 70 is well behind the front and shows the influence of the partially absorbing boundary condition at y = 0 where the population density is significantly reduced compared to its value at y = L where a no-flux boundary condition is imposed. Figure 1 and the similar subsequent figures were obtained from straightforward numerical simulations (spatial finite differences and explicit timestepping) and initial conditions that decayed exponentially in space. Implementation of the partially absorbing boundary conditions are straightforward in this case. For all results shown, the spatial and temporal resolution was sufficiently high that solutions were wellconverged. Part (b) of the figure plots the shape of the front u(x, y, t) at the two locations y = 0 and y = L, spaced by time intervals t = 2.0 showing the propagation into x > 0 at a constant velocity. At each timepoint and each point x, the value u(x, L , t) (shown by the red dash-dotted line) is larger than u(x, 0, t) (shown by the blue solid line). Figure 1c, d shows contours of constant u as time increases. In Fig. 1c, we observe that the front remains parallel to the boundary y = L where the no-flux condition is applied, but bends close to the partially absorbing boundary at y = 0. For the larger value for contours shown in Fig. 1d the contours do not extend all the way down to the boundary at y = 0 since the values of u there always remain below a value u ≈ 0.4423 in this case. This aspect is discussed further in Sect. 3.
Equation (6) can therefore be viewed as implicitly determining the front propagation speed c as a function of α. In particular it is interesting to note that (6) predicts that c = 0 when p = √ 2, i.e. on the line This result has two implications. Firstly, for any fixed positive α there is a minimum strip width L m (α) on which c = 0 and therefore below which the front cannot advance along the strip. Secondly, considering variations in L, there are two regimes: if then this calculation predicts that the front will always propagate forwards, for any value of α however large. But for L < L * m from this linearised theory at the leading edge we would predict that the front would propagate forwards for sufficiently small α and would propagate backwards for sufficiently large α, as long as a travelling wave solution exists for these parameter values. Note that lim α→∞ L m (α) = L * m . Figure 2 compares the results of the theoretical calculation at the leading edge (6) with speeds c estimated directly from numerical simulation. The front propagates significantly more slowly as α increases or when the strip becomes narrower, although for this range of parameter values the front always propagates forwards. It is interesting to note that, while the agreement is better for large strip widths L, the numerically computed value is always below the theoretical estimate. This is in contrast to the  (6) 1D FKPP problem for which numerical solutions (with compactly support initial conditions) produce stable fronts that travel at the lowest possible speed.

Two Partially Absorbing Boundaries
The solution of (4) in the case of a mixed boundary conditions v y = αv at y = 0 and v y = −βv at y = L (with α, β > 0) is similarly straightforward: the solution for v(y) still takes the form (5) but now the parameter p is the smallest positive solution of the equation For fixed β > 0 the line L 0 (α) in the (L , α) plane along which the speed c of the front is predicted to fall to zero is given by setting p = √ 2 in (9) which gives Thus the critical strip width L * m (β), below which the leading edge calculation indicates that the front speed could become zero for a sufficiently large value of α, is given by where tan −1 P denotes the principal value of the inverse tangent function. Figures 3 and 4 illustrate the propagation of fronts in narrow channels, with widths L = 3.0 and L = 2.0, respectively. The slowdown in the front propagation is clear from comparisons between parts (c) and (d) in each figure, while the linear approximation appears still to agree well with the numerically computed solution ahead of the front where u 1. It is also interesting to note that the shape of the contours of constant u have a very different appearance depending on the level at which they are plotted. Figures 3d and 4d show that if the contour value is high, the contour may not extend to anywhere near the partially absorbing boundary: the front more resembles a blob that propagates through the centre of the channel and does not appear to touch the boundaries at all. This is potentially important for observing experimental results where the imaging technique may similarly not allow direct observation of the shape of the front at low values of u and hence might not follow the front shape all the way down to the boundary itself.
For completeness, we note that if α = β then this calculation predicts that a front can propagate through the domain only if the domain size L satisfies

Solution in the Half-Space y ≥ 0
Now we turn to the case of (4) subject to only one boundary condition, v y = αv at y = 0. Far away from y = 0 we would expect that the front solution resembles that for the 1D problem, and travels with speed c = 2. Substituting the solution ansatz u(x, y, t) = u 0 e −(x−2t) f (y) into (3) we require that f yy = 0, and f (0) = α f (0) which implies f (y) = f 0 (1 + αy) with f 0 an undetermined constant. Hence in the half-space y ≥ 0 the front always travels with speed c = 2 independent of the value of the parameter α.
The solution to (3) allows us to examine the shape of level sets u = const near y = 0, far ahead of the front. Setting u = u 0 and t = 0 without loss of generality, we find that the level set u = u 0 is described by the curve x = log(1 + αy) which has two features of note: even in the region y 1 away from the boundary, the level set curves do not become asymptotically parallel to the y-axis, and also, expanding the solution for small x and y near the y = 0 boundary we see that the level set takes the form y = 1 α x here, so that level set curves meet the x-axis at an angle θ determined by α through the relation tan θ = 1/α. This behaviour near the boundary, with the level sets making the angle tan −1 (1/α) with the x-axis, arises in the both the problems presented above for solutions in a strip 0 ≤ y ≤ L; this behaviour is determined by the partially absorbing boundary condition alone and does not depend, as one might intuitively hope, on the choice of 'far-field' boundary condition.

Behind the Front
Since there is no closed-form explicit solution even for the 1D FKPP equation, it is not possible to follow the evolution of the front solutions for the 2D problem into the fully nonlinear regime. However, it is possible to study the transverse (i.e. y-direction) structure of solutions far behind the front; we do this in Sect. 3.1 for each of the cases considered in the previous section. The simplest of these cases is the half-space y ≥ 0, which we consider first, and for which an explicit solution profile can be obtained. We then consider solutions confined to a strip 0 ≤ y ≤ L for which the existence of unique monotonic fronts is clear geometrically, and for which an implicit form of the solution is available. In Sect. 3.2 we consider the stability of solutions, looking in particular at the regime where the population density behind the front becomes very low.

Existence of Solutions
Consider solutions to (1) with boundary condition (2); far behind the front the solution u(x, y, t) becomes independent of x and t since in the 1D problem solutions tend to the constant 1. Hence we write u(x, y, t) = w(y) where w(y) is required to solve the nonlinear ODE subject to the boundary conditions w y = αw at y = 0, and w → 1 as y → ∞.
Introducing the new variable z = dw/dy allows us to write (13) as the two-dimensional Hamiltonian system where the Hamiltonian H is given by and the Jacobian matrix at a point (w, z) in phase space is The 2D dynamical system (14) has two equilibrium points (0, 0) and (1, 0) which are a centre and a saddle point, respectively, and a sketch of the phase portrait is shown in Fig. 5. The required boundary conditions are satisfied for the unique trajectory of (14) that satisfies z = αw at y = 0 and (w, z) → (1, 0) as y → ∞ and for which both w and z are positive for all y ≥ 0. These conditions hold for the unique trajectory that forms a section of the stable manifold of the saddle point at (1, 0), starting at the point in the w > 0, z > 0 quadrant where the line z = αw intersects this stable manifold, illustrated in Fig. 5 by the intersection at the point (w 0 , z 0 ) of the (blue) dashed line and the homoclinic orbit (thick curve) which forms the stable manifold of (1, 0). The homoclinic orbit has the analytic expression for the transverse structure of the FKPP front, where y 0 is an as yet undetermined constant of integration. The form (17) clearly satisfies the condition w → 1 as y → ∞.
Fixing the condition w y = αw at y = 0 enables us to determine y 0 , implicitly, via  (17) where y 0 is determined as a function of α by solving (18) (solid black line) and numerical results (blue dots). a α = 0.2; b α = 0.9 Since the left-hand side of (18) increases monotonically from zero and asymptotically approaches 3/2 as y 0 increases to infinity, and the right-hand side increases monotonically from −α/2 at y 0 = 0 to very large positive values as y 0 increases, it is clear that for any fixed α > 0 Eq. (18) has a unique positive solution y 0 . Indeed in is clear by inspection of (18) that every such positive solution must satisfy y 0 ≥ 2cosh −1 √ 3/2 > 0. Figure 6 illustrates the form of the transverse structure of the front for two cases: α = 0.2 and α = 0.9. For the larger value of α the distortion is much more pronounced.
For the case of a strip 0 ≤ y ≤ L with either a no-flux or a partially absorbing boundary condition at y = L, the solution near y = 0 is close to the solution for a half-plane. The analytic solution in this case can be constructed implicitly using the phase portrait shown in Fig. 5 since the finite width L now gives an additional constraint which selects a unique (part of a) trajectory γ that must begin on the line z = αw (say at the point (w 2 , z 2 = αw 2 )) and end on the line z = −βw, (at the point, e.g. (w 3 , z 3 = −βw 3 )) as illustrated in Fig. 5.
The finite width constraint can be expressed through the requirement that the contour containing the trajectory must correspond to a value of H that satisfies which follows from the Hamiltonian (15), starting from the initial condition (w 2 , z 2 ) and setting H = H (w 2 , z 2 ) which is determined by the initial condition, and integrating along the trajectory γ starting at the point (w 2 , z 2 ) indicated by the blue solid circle and ending at the open circle at (w 3 , z 3 ) contained in the level set H = H (w 2 , z 2 ), as illustrated in Fig. 5. For a specified finite value of L, the constraint (19) enables us to select a trajectory in the phase plane that satisfies both this constraint and the required boundary conditions, showing that a solution indeed exists for any positive L.
The trajectory in Fig. 5 starting at the black square at (w 0 , z 0 ) and ending at the saddle point at (w 1 , z 1 ) = (1, 0) describes the solution for the half-space problem for which the solution at y = 0 must satisfy w y = αw, and for which w → 1 and w y → 0 as y → ∞.

Collapse Behind the Front
A necessary condition for the propagation of stable travelling waves is that the solution behind the front discussed in Sect. 3.1 is also stable. Numerical investigations reveal that this may not be the case in very narrow domains. A consequence of this instability is that travelling waves can no longer exist. In fact, in very narrow domains the speed of propagation of travelling waves departs significantly from the theoretical value computed at the leading edge, as illustrated in Fig. 7. The wavespeed falls to zero at a point (indicated in Fig. 7 by the black square) at which the solution behind the front also fails to remain positive as the front travels forwards.
For this value of α, the effect of the boundary condition becomes too significant for the strip to support a positive population behind the front. The critical point at which this occurs can be computed in a very similar way to the leading edge computation, since u is small so we can linearise, and also look for x-independent solutions, i.e. we investigate the linear stability of the solution u = 0 to the linearised equation subject to the boundary conditions (2). Setting u(y, t) = e σ t v(y) and defining q := √ 1 − σ leads to an eigenfunction of the form v(y) = A cos qy + B sin qy.  (21) After applying the boundary conditions at y = 0 and y = L we find that q must satisfy the condition The critical case is where the growth rate σ = 0, i.e. q = 1 which implies the relation for collapse to take place. For fixed β the line in the (L , α) plane along which collapse takes place can be written as It is instructive to compare the location of this line to the condition (10) deduced in Sect. 2 for the line on which the travelling wave speed is predicted to fall to zero, based on the calculation at the leading edge. A direct calculation shows that the inequality holds for all positive values of β and L, so that the collapse condition (21) always lies below the curve (10) in the (L , α) plane, as illustrated in Fig. 8. We conclude that the travelling wave ceases to exist due to the inability of the strip to sustain a positive population behind it, rather than because the invasion speed at the leading edge falls to zero.

Discussion
In this paper we have analysed the widely used FKPP model of front propagation to look at 2D front propagation near a boundary that models a constant fractional loss of the local population density. This extension of the 1D FKPP model has received very little attention in the literature previously. We are motivated in particular by the experiments reported by Möbius et al. (2015) who observed a change in the shape of the propagating front near boundaries that could be modelled as 'partially absorbing' in the way we set out above.
Our results show that in wide strips, the speed of propagation depends heavily on the nature of the boundary conditions, but is still well approximated by the linearised regime ahead of the front. The wavespeed that is computed directly from the linearised problem in a bounded strip tends to overestimate that obtained from numerical simulations. In a half-space the wavespeed remains constant but the shape of the wave near the boundary depends on the boundary condition parameter. Explicitly we derive the expression tan θ = 1/α for the angle θ that contours of constant u make with the reactive boundary at which we impose the boundary condition u y = αu.
If the strip is sufficiently narrow, then when the influence of the boundary conditions is sufficiently strong, the travelling wave ceases to exist: it collapses at the point at which the regime behind the front can no longer sustain a positive population. Moreover, in narrow strips the wavespeed predicted by the linearised calculation at the leading edge ceases to be a good approximation to the speed of travelling waves com-puted numerically: the wavespeed is clearly determined as part of the fully nonlinear nature of the problem.
It is worth commenting briefly on the collapse phenomenon since the linear stability of the trivial state behind the front appears to be closely related to discussions in the literature of 'critical patch sizes' for populations in finite domains, as discussed by various authors including Skellam (1951), Kierstead and Slobodkin (1953), and more recently Ryabov and Blasius (2008). In particular the KiSS model (using the notation of these authors for a population density P(x, t)) P t = μ 0 P + D P x x , on 0 ≤ x ≤ L P(0, t) = P(L , t) = 0, which these authors discuss is essentially the linearised FKPP equation with Dirichlet boundary conditions, i.e. the limiting version of (20) with reactive boundaries as α, β → ∞. That the KiSS model is able to sustain positive population density P(x, t) only if L ≥ π √ D/μ 0 is for exactly the same reasons as our model in the regime far behind the front. The nonlinear extension of the KiSS model, in the form with a logistic nonlinearity, has a stationary solution that can be expressed in terms of elliptic functions, and the critical patch size (i.e. domain length L) for persistence remains the same: we require L ≥ π √ D/μ 0 for persistence. Ludwig et al. (1979) considered an extension of the KiSS model in which the domain is the real line, but the mortality μ is negative outside a finite interval: For this model, Ludwig et al. (1979) find that a population can persist if the domain size L satisfies which for positive m and μ 0 implies that the population cannot exist in domains strictly smaller than √ D/μ 0 . Even for small m, a critical domain size exists and can be approximated by L ≈ 2 √ Dm/μ 0 when 0 < m 1. It is interesting to compare these results with conditions (7) and (11) for the case of reactive boundaries considered here. Condition (23) implies that for every choice of the mortality m and the other problem parameters there exists a critical patch size, tending to zero as m → 0, and L ≈ 2 √ m for small m, and tending to π as m → ∞. The same is true in the problem considered here: for any positive value of α, there is a minimum domain size given by solving (7) below which the front will collapse and no travelling wave solution appears to be possible. Hence (21) implicitly yields a critical patch size in the sense of the KiSS problem, and the critical patch size depends on the value of the boundary parameters α and β: the front collapses if these boundary parameters are sufficiently large.
In terms of future directions there are three avenues that we intend to investigate in further work. The first is to compare these solutions for the Fisher-KPP PDE with the results of stochastic simulations for front propagation, incorporating the effects of the partially absorbing boundary at the microscopic level (Erban and Chapman 2007). The second direction is to consider further extensions of these kinds of boundary condition, for example introducing spatially inhomogeneous boundary conditions. The third direction is to vary the functional form of the PDE model (1), as has been done by many authors; either replacing the nonlinear term with a continuous but piecewise-affine construction, or considering more general nonlinearities. A further avenue is of course to make rigorous the statements that we conjecture here, and to provide the same level of proof for the 2D FKPP equation as is available for its 1D counterpart. We intend to report on these directions in future papers.