Critical length for the spreading-vanishing dichotomy in higher dimensions

Matthew J. Simpson¹¹1To whom correspondence should be addressed. E-mail: [email protected] School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland 4001, Australia

Abstract

We consider an extension of the classical Fisher-Kolmogorov equation, called the Fisher-Stefan model, which is a moving boundary problem on $0<x<L(t)$ . A key property of the Fisher-Stefan model is the spreading-vanishing dichotomy, where solutions with $L(t)>L_{\textrm{c}}$ will eventually spread as $t\to\infty$ , whereas solutions where $L(t)\ngtr L_{\textrm{c}}$ will vanish as $t\to\infty$ . In one dimension is it well-known that the critical length is $L_{\textrm{c}}=\pi/2$ . In this work we re-formulate the Fisher-Stefan model in higher dimensions and calculate $L_{\textrm{c}}$ as a function of spatial dimensions in a radially symmetric coordinate system. Our results show how $L_{\textrm{c}}$ depends upon the dimension of the problem and numerical solutions of the governing partial differential equation are consistent with our calculations.

Keywords:

Reaction-diffusion; Fisher-Kolmogorov; Fisher’s Equation; Stefan condition; Moving boundary problem; travelling waves; Invasion; Extinction.

1 Introduction

The well-known Fisher-Kolmogorov model [1, 2, 3, 4, 5] is a reaction-diffusion equation that is often used to describe the spatial and temporal spreading of a population density where individuals in that population undergo random diffusive migration and logistic proliferation. The Fisher-Kolmogorov model is often written as

\dfrac{\partial u(x,t)}{\partial t}=D\dfrac{\partial^{2}u(x,t)}{\partial x^{2}}+\lambda u(x,t)\left[1-\dfrac{u(x,t)}{K}\right],

(1)

where $u(x,t)$ is the population density, $x$ is position, $t$ is time, $D$ is the diffusivity, $\lambda$ is the proliferation rate and $K$ is the carrying capacity density. The Fisher-Kolmogorov equation is often studied on an infinite domain, $-\infty<x<\infty$ , where it is well known to give rise to travelling wave solutions that eventually propagate with speed $c=2\sqrt{D\lambda}$ for initial conditions with compact support [1, 2, 3, 4, 5]. The Fisher-Kolmogorov equation, and certain extensions, have been used to model invasion fronts in ecology [6, 7, 8, 9, 10, 11] and cell biology [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. Despite the widespread application of this fundamental model, there are several shortcomings. First, any localised initial condition with compact support will always lead to population growth and successful colonisation. Therefore, an implicit assumption in the Fisher-Kolmogorov model is that the population will always survive and never go extinct. Second, solutions of the Fisher-Kolmogorov equation are smooth and without compact support. This means that the Fisher-Kolmogorov model does not, strictly speaking, lead to a clear and unambiguous invasion front. This can be problematic because well defined invasion fronts are often observed in practice [16, 17].

There are several ways that these two shortcomings of the Fisher-Kolmogorov model have been addressed in the applied mathematics literature. One common approach is to consider an extension of Equation (1) that incorporates degenerate nonlinear diffusion since this leads to travelling wave solutions with a well-defined front [4, 16, 17, 32, 33, 34, 35, 36, 37, 38, 39]. Another extension that has received less attention in the applied mathematics literature, but far more attention in the analysis literature, is to re-cast Equation (1) as a moving boundary problem,

\dfrac{\partial u(x,t)}{\partial t}=D\dfrac{\partial^{2}u(x,t)}{\partial x^{2}}+\lambda u(x,t)\left[1-\dfrac{u(x,t)}{K}\right],

(2)

on $0<x<L(t)$ , with $\partial u(0,t)/\partial x=0$ and $u(L(t),t)=0$ , ensuring that there is a well-defined front at the moving boundary, $x=L(t)$ . To close the problem, the evolution of the moving boundary is taken to follow a classical Stefan condition, $\textrm{d}L(t)/\textrm{d}t=-\hat{\kappa}\partial u(L(t),t)/\partial x$ [40, 41, 42, 43, 44]. As we will show, the adoption of Equation (2) alleviates both the shortcomings of the classical Fisher-Kolmogorov model since the this moving boundary analogue leads to solutions with well-defined fronts, as well as permitting certain initial conditions to become extinct.

The moving boundary analogue of the Fisher-Kolmogorov model has been called the Fisher-Stefan model [45]. To simplify our analysis we rescale the variables: $x^{\prime}=x/\sqrt{D/\lambda}$ , $t^{\prime}=t\lambda$ and $u^{\prime}=u/K$ . Dropping the prime notation gives

\dfrac{\partial u(x,t)}{\partial t}=\dfrac{\partial^{2}u(x,t)}{\partial x^{2}}+u(x,t)\left[1-u(x,t)\right],

(3)

on $0<x<L(t)$ , with $\partial u(0,t)/\partial x=0$ , $u(L(t),t)=0$ with $\textrm{d}L(t)/\textrm{d}t=-\kappa\partial u(L(t),t)/\partial x$ . In this nondimensional framework there is just one parameter, $\kappa>0$ . Varying $\kappa$ allows us to control the finite speed of the moving boundary at $x=L(t)$ . This mathematical model has been studied extensively in the analysis literature [46, 47, 48, 49, 50, 51, 52]. Much of this analysis has centred upon proving properties of travelling wave solutions of Equation (3) as well as characterising the spreading-vanishing dichotomy. As we will show, it is very interesting that certain initial conditions and choices of $\kappa$ in Equation (3) lead to eventual spreading in the form of a travelling wave, whereas other choices of initial condition and $\kappa$ in Equation (3) lead to extinction. A key feature of the spreading-vanishing dichotomy is that there exists a critical length, $L_{\textrm{c}}$ , such that solutions of Equation (3) that evolve in such a way that $L(t)>L_{\textrm{c}}$ always lead to eventual spreading, whereas solutions that evolve in such a way so that $L(t)\ngtr L_{\textrm{c}}$ always leads to eventual extinction. A great deal of attention has been paid to the analysis of this spreading-vanishing dichotomy, and establishing that $L_{\textrm{c}}=\pi/2$ [46, 47, 48, 49, 50, 51, 52].

The aim of this work is to revisit the spreading-vanishing dichotomy in a more general setting by re-casting Equation (3) in a radially symmetric geometry. Using numerical simulations to guide our calculations, we provide a very straightforward interpretation of $L_{\textrm{c}}$ , and show how the critical length depends upon the dimension of the problem. All results are supported by numerical evidence. The details of the numerical scheme we use to solve Equation (3) is described in the Appendix, and MATLAB software to implement the numerical scheme and to explore and visualise various solutions is available on GitHub.

2 Results and Discussion

We will first revisit the spreading-vanishing dichotomy in a one-dimensional coordinate system before moving on to show that the calculations can be extended to higher dimensions.

2.1 Preliminary results in one dimension

We first begin by visualising numerical solutions of Equation (3). For simplicity, all numerical solutions in this work have the same initial condition: $u(x,0)=1$ for $0\leq x\leq 1$ with $L(0)=1$ . The solution of Equation (3) with $\kappa=1$ is given in Figure 1(a)-(d) showing that the initial density evolves into a constant speed, constant shape travelling wave. Insights into such travelling waves can be obtained by phase plane analysis [45]. The evolution of $L(t)$ is given in Figure 1(d) where we see that $L(t)$ eventually increases linearly as $t\to\infty$ .

Refer to caption — Figure 1: Numerical solutions of Equation (3). Results in (a)-(d) are for $\kappa=1$ and lead to a spreading travelling wave. Results in (e)-(h) are for $\kappa=0.2$ and lead to extinction. (d) and (h) show the evolution of $L(t)$ (blue) and $L_{\textrm{c}}=\pi/2$ (red).

With $\kappa=0.2$ the solution of Equation (3) is given in Figure 1(e)-(h). These solutions show that the initial density evolves in such a way that the population goes extinct as $t\to\infty$ . Solutions at $t=50$ and $t=100$ (Figure 1(f)-(g)) are plotted with an inset because these solutions are so close to zero that it is impossible to clearly visualise the solutions when they are plotted on the same scale as the other solutions in Figure 1. The evolution of $L(t)$ in Figure 1(h) shows that $L(t)$ initially increases, and then asymptotes to some constant value as $t\to\infty$ . In this case we always have $\partial u(x,t)/\partial x<0$ at $L(t)$ , and so the Stefan condition ensures that $\textrm{d}L(t)/\textrm{d}t\geq 0$ . This means that $L(t)$ never decreases.

This simple numerical demonstration in Figure 1 illustrates the key features of the spreading-vanishing dichotomy. Unlike the usual Fisher-Kolmogorov model where all positive initial conditions will eventually lead to spreading, here the moving boundary analogue of this model can either lead to an invading travelling wave or it can lead to extinction. Numerical solutions of Equation (2) in Figure 1 provide some intuition about how we may analyse this dichotomy. As illustrated in Figure 1, When extinction occurs we always have $u(x,t)\to 0$ as $t\to\infty$ for all $0\leq x\leq L(t)$ , suggesting that we can analyse the behaviour in this regime by linearising. If $v(x,t)\ll 1$ , an appropriate linear analogue of the Fisher-Stefan model is

\dfrac{\partial v(x,t)}{\partial t}=\dfrac{\partial^{2}v(x,t)}{\partial x^{2}}+v(x,t),

(4)

on $0<x<L(t)$ with $\partial v(0,t)/\partial x=0$ and $v(L(t),t)=0$ . For the purposes of this analysis we treat the domain length as fixed, $L(t)=L$ . The exact solution of Equation (4) is

v(x,t)=\sum_{n=1}^{\infty}A_{n}\cos(\lambda_{n}x)\textrm{e}^{-t(\lambda_{n}^{2}-1)},

(5)

where $\lambda_{n}=\pi(2n-1)/(2L)$ for $n=1,2,3,\ldots$ , and the coefficients $A_{n}$ can be chosen so that the solution matches an initial condition. As we will show, our analysis does not require that we specify these coefficients. We now approximate Equation (5) by retaining just the first term in the series. Such leading eigenvalue approximations can be particularly accurate as $t\to\infty$ [53, 54, 55],

V(x,t)\sim A_{1}\cos\left(\dfrac{\pi x}{2L}\right)\textrm{e}^{-t\left[\left(\dfrac{\pi}{2L}\right)^{2}-1\right]}.

(6)

With this solution we can specify a conservation statement for the time evolution of the total population in the domain,

\dfrac{\textrm{d}m}{\textrm{d}t}=\underbrace{\int_{0}^{L}V(x,t)\,\textrm{d}x}_{\textrm{accumulation due to source term}}-\underbrace{-\dfrac{\partial V(L,t)}{\partial x}}_{\textrm{loss due to the diffusive flux at $x=L$}}.

(7)

Setting $\textrm{d}m/\textrm{d}t=0$ we solve Equation (7) using $V(x,t)$ , giving $L_{\textrm{c}}=\pi/2\approx 1.570796327$ , which is independent of $A_{n}$ . This means that whenever a solution of Equation (3) evolves in such a way that $L(t)>\pi/2$ , the accumulation due to the source term is greater than the loss due to the diffusive flux at the moving boundary. In contrast, whenever $L(t)<\pi/2$ , the accumulation due to the source term is less than the loss at the moving boundary. This result has some interesting consequences and provides a straightforward explanation for our numerical results in Figure 1. For our initial condition we have $L(0)=1<\pi/2$ and it is unclear whether or not the transient solution will evolve such that $L(t)>\pi/2$ as $t\to\infty$ . When $\kappa$ is sufficiently large we see that the evolution of the solution is such that $L(t)$ eventually exceeds $\pi/2$ and at this point the population will continue to grow indefinitely, as shown in Figure 1(d). In contrast, for sufficiently small $\kappa$ , the evolution of the solution is such that $L(t)$ that never reaches $\pi/2$ and so the diffusive flux out of the domain is greater than the production due to the source term and we see eventual extinction where $\displaystyle{\lim_{t\to\infty}L(t)<\pi/2}$ , as shown in Figure 1(h).

Our numerical results and analysis of the linearised problem confirm that we have $L_{\textrm{c}}=\pi/2$ in one dimension and this explains our original choice of initial condition in Figure 1 where $L(0)=1<L_{\textrm{c}}$ . In this case $L(0)<L_{\textrm{c}}$ and without solving for the time dependent solution it is unclear whether $L(t)$ evolves in such a way that it ever exceeds $L_{\textrm{c}}$ . Therefore, with this initial condition we have the flexibility of choosing $\kappa$ to be sufficiently small such that $L(t)$ never reaches $L_{\textrm{c}}$ and the population goes extinct, or we can choose $\kappa$ to be sufficiently large that $L(t)$ will eventually exceed $L_{\textrm{c}}$ and lead to a travelling wave. Had we chosen a different initial condition with $L(0)>L_{\textrm{c}}$ we would always observe the formation of a travelling wave provided $\kappa>0$ . Since we are interested in studying the spreading-extinction dichotomy we restrict our numerical solutions to the case where $L(0)<L_{\textrm{c}}$ . Setting $L(0)=1$ is a convenient way to achieve this.

While it has long been known that $L_{\textrm{c}}=\pi/2$ in the more formal literature [46, 47, 48, 49, 50, 51, 52], here we offer a very simple and intuitive way to calculate and interpret the spreading-vanishing dichotomy in terms of this critical length. Our approach is mathematically and conceptually straightforward, and as we will now show, also applies in other geometries.

2.2 Critical length in higher dimensions

A generalisation of Equation (3) is to consider the Fisher-Stefan model in a radially symmetric coordinate system,

\dfrac{\partial u(x,t)}{\partial t}=\dfrac{\partial^{2}u(x,t)}{\partial x^{2}}+\dfrac{d-1}{x}\dfrac{\partial u(x,t)}{\partial x}+u(x,t)\left[1-u(x,t)\right],

(8)

for $t>0$ and $0<x<L(t)$ , with $\partial u(0,t)/\partial x=0$ , $u(L(t),t)=0$ , and $\textrm{d}L(t)/\textrm{d}t=-\kappa\partial u(L(t),t)/\partial x$ . Here $d=1,2,3$ is the dimension. We will now explore how $L_{\textrm{c}}$ depends upon $d$ .

Setting $d=2$ in Equation (8) allows us to consider the spreading-vanishing dichotomy on a disc. Preliminary numerical simulations indicate that for a particular choice of $u(x,0)$ the solution may either evolve into a moving front when $\kappa$ is sufficiently large, or the solution may go extinct when $\kappa$ is sufficiently small, just like we demonstrated in Figure 1 when $d=1$ . To explore this we linearise, with $v(x,t)\ll 1$ , to obtain

\dfrac{\partial v(x,t)}{\partial t}=\dfrac{\partial^{2}v(x,t)}{\partial x^{2}}+\dfrac{1}{x}\dfrac{\partial v(x,t)}{\partial x}+v(x,t),

(9)

on $0<x<L$ , with $\partial v(0,t)/\partial x=0$ and $v(L,t)=0$ . The exact solution of Equation (9) is

v(x,t)=\sum_{n=1}^{\infty}A_{n}J_{0}\left(x\sqrt{1+\lambda_{n}^{2}}\right)\textrm{e}^{-\lambda_{n}^{2}t},

(10)

where $J_{0}(x)$ is the zeroth-order Bessel function of the first kind. Here the eigenvalues, $\lambda_{n}$ , are obtained by setting $L\sqrt{1+\lambda_{n}^{2}}$ equal to the $n^{\textrm{th}}$ zeros of $J_{0}(x)$ . If we invoke a leading eigenvalue approximation we have $\sqrt{1+\lambda_{1}^{2}}=z/L$ , where $z$ is the first zero of $J_{0}(x)$ , or $z\approx 2.404825558$ . The leading eigenvalue approximation is

V(x,t)\sim A_{1}J_{0}\left(\dfrac{zx}{L}\right)\textrm{e}^{-t\left[\left(\dfrac{z}{L}\right)^{2}-1\right]}.

(11)

To proceed we formulate a condition where the accumulation due to the source term is precisely matched by the loss at the moving boundary. Taking care to account for the geometry we obtain,

\underbrace{\int_{0}^{L}V(x,t)2\pi x\,\textrm{d}x}_{\textrm{accumulation due to source term}}=\underbrace{-\dfrac{\partial V(L,t)}{\partial x}2\pi L}_{\textrm{loss due to the diffusive flux at $x=L$}},

(12)

in which we can substitute Equation (11) and solve for $L$ . This procedure shows that $L_{\textrm{c}}$ is the first zero of $J_{0}(x)$ , giving $L_{\textrm{c}}\approx 2.404825558$ . Numerical solutions in Figure 2 confirm this.

Setting $d=3$ in Equation (8) allows us to consider the spreading-vanishing dichotomy on a sphere. To explore the critical length in this case we linearise, with $v(x,t)\ll 1$ , to obtain

\dfrac{\partial v(x,t)}{\partial t}=\dfrac{\partial^{2}v(x,t)}{\partial x^{2}}+\dfrac{2}{x}\dfrac{\partial v(x,t)}{\partial x}+v(x,t),

(13)

on $0<x<L$ , with $\partial v(0,t)/\partial x=0$ and $v(L,t)=0$ . The exact solution of Equation (13) is

v(x,t)=\sum_{n=1}^{\infty}\dfrac{A_{n}}{x}\sin\left(x\sqrt{1+\lambda_{n}^{2}}\right)\textrm{e}^{-\lambda_{n}^{2}t},

(14)

where $\lambda_{n}=n\pi/L$ , for $n=1,2,3,\ldots$ . The associated leading eigenvalue approximation is

V(x,t)\sim\dfrac{A_{1}}{x}\sin\left[\dfrac{\pi x}{L}\right]\textrm{e}^{-t\left[\left(\dfrac{\pi}{L}\right)^{2}-1\right]}.

(15)

The associated conservation statement is,

\underbrace{\int_{0}^{L}V(x,t)4\pi x^{2}\,\textrm{d}x}_{\textrm{accumulation due to source term}}=\underbrace{-\dfrac{\partial V(L,t)}{\partial x}4\pi L^{2}}_{\textrm{loss due to the diffusive flux at $x=L$}},

(16)

in which we substitute Equation (15) and solve for $L$ , giving $L_{\textrm{c}}=\pi\approx 3.141592654$ Numerical solutions in Figure 3 confirm this.

3 Conclusions

In this work we explore an extension of the well-known Fisher-Kolmogorov model, called the Fisher-Stefan model. This extension involves reformulating the Fisher-Kolmogorov model as a moving boundary problem, and the evolution of the moving boundary is described by a classical Stefan condition. The Fisher-Stefan model has two key advantages over the usual Fisher-Kolmogorov model: (i) solutions of the Fisher-Stefan model are characterised by a sharp front, which is often relevant in practice; and, (ii) solutions of the Fisher-Stefan model do not always lead to successful colonisation since the population can go extinct as $t\to\infty$ . An important feature of the Fisher-Stefan model is the existence of a critical length, $L_{\textrm{c}}$ . As we show, solutions of the Fisher-Stefan model that evolve in such a way that $L(t)>L_{\textrm{c}}$ always lead to spreading, whereas solutions that evolve in such a way so that $L(t)\ngtr L_{\textrm{c}}$ always leads to extinction. Formal analysis of the spreading-vanishing dichotomy has established that $L_{\textrm{c}}=\pi/2$ in a standard one-dimensional Cartesian geometry [46, 47, 48, 49, 50, 51, 52]. This property is often called the spreading-vanishing dichotomy, and it has been extensively studied in one-dimension. Far less attention has been paid to similar behaviour in higher dimensions.

In this work we re-formulate the Fisher-Stefan model in a radially symmetric coordinate system and study the spreading-vanishing dichotomy in one, two and three dimensions. Using simple conservation arguments we confirm that $L_{\textrm{c}}=\pi/2$ in the usual one-dimensional geometry, and this result is confirmed by numerical solutions. By extending these ideas to consider the Fisher-Stefan model on a disc and sphere we find that the critical length depends upon the geometry of the problem. In particular the critical length for the Fisher-Stefan model on a disc is the first zero of the zeroth-order Bessel function of the first kind ( $L_{\textrm{c}}\approx 2.404825558$ ) and the critical length for the Fisher-Stefan model on a sphere is $\pi$ ( $L_{\textrm{c}}\approx 3.141592654$ ). Numerical solutions of the governing partial differential equation are consistent with these findings.

Our approach for calculating the critical length is conceptually straightforward, and shows how the critical length depends upon the dimension. A great deal of literature about the spreading-vanishing dichotomy involves rigorous analysis that can be both difficult to interpret and difficult to extend to more general problems. Our approach, in contrast, is quite flexible and there are many potential extensions. While we have not pursed further calculations for other choices of $d$ , it is straightforward to apply the methods described here for $d>3$ to explore the spreading-vanishing dichotomy on an hypersphere. For example, setting $d=4$ in Equation (8) and repeating the analysis shows that the critical length is the first zero of the first-order Bessel function of the first kind ( $L_{\textrm{c}}\approx 3.831705970$ ). Further calculations for different $d$ can also be performed.

We close by commenting on some extensions of Equation (8) that could be of interest. As we already pointed out in the Introduction, a very common extension of the classical Fisher-Kolmogorov model is to generalise the linear diffusion term to a degenerate nonlinear diffusion term [4, 16, 17, 32, 33, 34, 35, 36, 37, 38, 39]. This generalisation is often motivated by desire to seek solutions with a well-defined sharp front [4, 16, 17, 32, 33, 34, 35, 36, 37, 38, 39]. While it is certainly possible generalise Equation (8) to include degenerate nonlinear diffusion [56], an important consequence of this in the moving boundary context is that this extension does not lead to a spreading-vanishing dichotomy. This is because the flux at the moving boundary $x=L(t)$ , where $u(L(t),t)=0$ , is always zero thereby always preventing extinction. In contrast, the one-phase Stefan condition in Equation (8) leads to a loss of mass at the moving boundary and it is this loss of mass that enables the population to go extinct. Another interesting extension would be to consider the spreading-vanishing dichotomy in higher dimensions for non radial geometries and explore the relationship between the solution in radial and non-radial geometries [57, 58]. Such an extension to consider moving boundary problems in higher dimensions without radial symmetry would require a very different approach numerically and so we leave this for future consideration.

4 Appendix: Numerical Methods

We develop a simple, robust numerical solution of

\dfrac{\partial u(x,t)}{\partial t}=\dfrac{\partial^{2}u(x,t)}{\partial x^{2}}+\dfrac{d-1}{x}\dfrac{\partial u(x,t)}{\partial x}+u(x,t)\left[1-u(x,t)\right],

(17)

for $t>0$ and $0<x<L(t)$ , with $u(L(t),t)=0$ , $\partial u(0,t)/\partial x=0$ and $\textrm{d}L(t)/\textrm{d}t=-\kappa\partial u(L(t),t)/\partial x$ . To solve this partial differential equation we first re-cast the moving boundary problem on a fixed domain by setting $\xi=x/L(t)$ , giving

\dfrac{\partial u(\xi,t)}{\partial t}=\dfrac{1}{L^{2}(t)}\dfrac{\partial^{2}u(\xi,t)}{\partial\xi^{2}}+\left[\dfrac{\xi}{L(t)}\dfrac{\textrm{d}L(t)}{\textrm{d}t}+\dfrac{d-1}{\xi L^{2}(t)}\right]\dfrac{\partial u(\xi,t)}{\partial\xi}+u(\xi,t)\left[1-u(\xi,t)\right],

(18)

for $t>0$ and $0<\xi<1$ , with $u(1,t)=0$ , $\partial u(0,t)/\partial\xi=0$ and $\textrm{d}L(t)/\textrm{d}t=-(\kappa/L(t))\partial u(1,t)/\partial\xi$ .

We discretise Equation (18) on a uniform mesh with spacing $\Delta\xi$ , and the value of $u(\xi,t)$ at each node is indexed with a subscript $i=1,2,3,\ldots,N$ , where $N=1/\Delta\xi+1$ . Integrating through time using an implicit Euler approximation, using a subscript $m$ to denote the time step, leads to the following system of $N$ nonlinear algebraic equations,

$\displaystyle u_{1}^{m}$	$\displaystyle=u_{2}^{m},$
$\displaystyle\dfrac{u_{i}^{m}-u_{i}^{m-1}}{\Delta t}$	$\displaystyle=\dfrac{\left(u_{i-1}^{m}-2u_{i}^{m}+u_{i+1}^{m}\right)}{(\Delta\xi L^{m})^{2}}+\left[\dfrac{\xi_{i}}{L^{m}}\dfrac{\textrm{d}L^{m}}{\textrm{d}t}+\dfrac{d-1}{\xi_{i}(L^{m})^{2}}\right]\dfrac{\left(u_{i+1}^{m}-u_{i+1}^{m}\right)}{2\Delta\xi}$
	$\displaystyle+u_{i}^{m}\left[1-u_{i}^{m}\right],\quad i=2,3,4,\ldots N-1,$
$\displaystyle u_{N}^{m}$	$\displaystyle=0,$	(19)

where $\textrm{d}L^{m}/\textrm{d}t=\kappa u_{N-1}^{m}/(L^{m}\Delta\xi)$ . This system of nonlinear algebraic equations is solved using Newton-Raphson iteration until the change in estimates of $u_{i}^{m}$ per iteration fall below a tolerance of $\varepsilon$ . This algorithm provides estimates of $u(\xi,t)$ on $0<\xi<1$ which are then re-scaled to give $u(x,t)$ on $0<x<L(t)$ . A MATLAB implementation of this algorithm is available on GitHub. All results presented in this work use $\Delta\xi=1\times 10^{-3}$ , $\Delta t=1\times 10^{-4}$ and $\varepsilon=1\times 10^{-10}$ . This choice of discretisation appears to be suitable for our purposes, but numerical results should always be checked to test that they are independent of the mesh. Experience suggests that problems with larger $\kappa$ require particular care to ensure grid convergence.

Acknowledgements

This work was supported by the Australian Research Council (DP170100474). I appreciate many informal discussions on this topic with Scott McCue, Yihong Du, Wang Jin and Maud El-Hachem, together with helpful suggestions of three anonymous referees.

References

[1] Fisher RA (1937) The wave of advance of advantageous genes. Annals of Eugenics. 7, 355–369.
[2] Kolmogorov AN, Petrovskii IG, Piskunov NS (1937) A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem. Moscow University Mathematics Bulletin. 1, 1–26.
[3] Canosa J (1973) On a nonlinear diffusion equation describing population growth. IBM Journal of Research and Development. 17, 307–313.
[4] Murray JD (2002) Mathematical biology I. An introduction. New York: Springer.
[5] Grindrod P (2007) Patterns and waves. Oxford: Oxford University Press.
[6] Skellam JG (1951) Random dispersal in theoretical populations. Biometrika. 38, 196–218.
[7] Shigesada N, Kawasaki K, Takeda Y (1995) Modeling stratified diffusion in biological invasions. American Naturalist. 146, 229–251.
[8] Steele J, Adams J, Sluckin T (1998) Modelling paleoindian dispersals. World Archaeology. 30, 286–305.
[9] Forbes LK (1997) A two-dimensional model for large-scale bushfire spread. Journal of the Australian Mathematical Society Series B: Applied Mathematics. 39, 171-194.
[10] Broadbridge P, Bradshaw BH, Fulford GR, Aldis GK (2002) Huxley and Fisher equations for gene propagation: An exact solution. ANZIAM Journal. 44, 11-20.
[11] Bradshaw-Hajek BH, Broadbridge P (2004) A robust cubic reaction-diffusion system for gene propagation. Mathematical and Computer Modelling. 39, 1151-1163.
[12] Sherratt JA, Murray JD (1990) Models of epidermal wound healing. Proceedings of the Royal Socieyt of London: Series B. 241, 29–36.
[13] Gatenby RA, Gawlinski ET (1996) A reaction-diffusion model of cancer invasion. Cancer Research. 56, 5745–5753.
[14] Painter KJ, Sherratt JA (2003) Modelling the movement of interacting cell populations. Journal of Theoretical Biology. 225, 327–339.
[15] Swanson KR, Bridge C, Murray JD, Alvord EC Jr (2003) Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. Journal of Neurological Sciences. 216, 1–10.
[16] Maini PK, McElwain DLS, Leavesley DI (2004) Traveling wave model to interpret a wound-healing cell migration assay for human peritoneal mesothelial cells. Tissue Engineering. 10, 475–482.
[17] Maini PK, McElwain DLS, Leavesley D (2004) Travelling waves in a wound healing assay. Applied Mathematics Letters 17, 575–580.
[18] Simpson MJ, Landman KA, Hughes BD, Newgreen DF (2006) Looking inside an invasion wave of cells using continuum models: Proliferation is the key. Journal of Theoretical Biology. 243, 343–360.
[19] Sengers BG, Please CP, Oreffo ROC (2007) Experimental characterization and computational modelling of two-dimensional cell spreading for skeletal regeneration. Journal of the Royal Society Interface. 4, 1107–1117.
[20] Simpson MJ, Zhang DC, Mariani M, Landman KA, Newgreen DF (2007) Cell proliferation drives neural crest cell invasion of the intestine. Developmental Biology. 302, 553–568.
[21] Swanson KR, Rostomily RC, Alvord EC Jr (2008) A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle. British Journal of Cancer. 98, 113–119.
[22] Simpson MJ, Treloar KK, Binder BJ, Haridas P, Manton KJ, Leavesley DI, McElwain DLS, Baker RE (2013) Quantifying the roles of motility and proliferation in a circular barrier assay. Journal of the Royal Society Interface. 10, 20130007.
[23] Treloar KK, Simpson MJ, McElwain DLS, Baker RE (2014) Are in vitro estimates of cell diffusivity and cell proliferation rate sensitive to assay geometry? Journal of Theoretical Biology. 356, 71–84.
[24] Johnston ST, Shah ET, Chopin LK, McElwain DLS, Simpson MJ (2015) Estimating cell diffusivity and cell proliferation rate by interpreting IncuCyte ZOOM^TM assay data using the Fisher-Kolmogorov model. BMC Systems Biology. 9, 38.
[25] Johnston ST, Ross JV, Binder BJ, McElwain DLS, Haridas P, Simpson MJ (2016) Quantifying the effect of experimental design choices for in vitro scratch assays. Journal of Theoretical Biology. 400, 19–31.
[26] Jin W, Shah ET, Penington CJ, McCue SW, Chopin LK, Simpson MJ (2016) Reproducibility of scratch assays is affected by the initial degree of confluence: experiments, modelling and model selection. Journal of Theoretical Biology. 390, 136–145.
[27] Nardini JT, Chapnick DA, Liu X, Bortz DM (2016) Modeling keratinocyte wound healing dynamics: Cell–cell adhesion promotes sustained collective migration. Journal of Theoretical Biology. 400, 103–117.
[28] Haridas P, Penington CJ, McGovern JA, McElwain DLS, Simpson MJ (2017) Quantifying rates of cell migration and cell proliferation in co-culture barrier assays reveals how skin and melanoma cells interact during melanoma spreading and invasion. Journal of Theoretical Biology. 423, 13–25.
[29] Vittadello ST, McCue SW, Gunasingh G, Haass NK, Simpson MJ (2018) Mathematical models for cell migration with real-time cell cycle dynamics. Biophysical Journal. 114, 1241–1253.
[30] Browning AP, Haridas P, Simpson MJ (2019) A Bayesian sequential learning framework to parameterise continuum models of melanoma invasion into human skin. Bulletin of Mathematical Biology. 81, 676–698.
[31] Warne DJ, Baker RE, Simpson MJ (2019) Using experimental data and information criteria to guide model selection for reaction–diffusion problems in mathematical biology. Bulletin of Mathematical Biology. 81, 1760–1804.
[32] Harris S (2004) Fisher equation with density-dependent diffusion: special solutions. Journal of Physics A: Mathematical and Theoretical. 37, 6267.
[33] Sánchez Garduno F, Maini PK (1994) An approximation to a sharp type solution of a density-dependent reaction–diffusion equation. Applied Mathematics Letters. 7, 47-51.
[34] Sánchez Garduno F, Maini PK (1994) Traveling wave phenomena in some degenerate reaction–diffusion equations. Journal of Differential Equations. 117, 281-319.
[35] Witelski TP (1995) Merging traveling waves for the Porous–Fisher’s equation. Applied Mathematics Letters. 8, 57–62.
[36] Sherratt JA, Marchant BP (1996) Nonsharp travelling wave fronts in the Fisher equation with degenerate nonlinear diffusion. Applied Mathematics Letters. 9, 33–38.
[37] Simpson MJ, Baker RE, McCue SW (2011) Models of collective cell spreading with variable cell aspect ratio: A motivation for degenerate diffusion models. Physical Review E. 83, 021901.
[38] Baker RE, Simpson MJ (2012) Models of collective cell motion for cell populations with different aspect ratio: diffusion, proliferation and travelling waves. Physica A: Statistical Mechanics and its Applications. 391, 3729–3750.
[39] McCue SW, Jin W, Moroney TJ, Lo K-Y, Chou SE, Simpson MJ (2019) Hole-closing model reveals exponents for nonlinear degenerate diffusivity functions in cell biology. Physica D: Nonlinear Phenomena. 398, 130–140.
[40] Crank J (1987) Free and moving boundary problems. Oxford University Press, Oxford.
[41] Gupta SC (2017) The classical Stefan problem. Basic concepts, modelling and analysis with quasi-analytical solutions and methods. Second edition. Elsevier, Amsterdam.
[42] McCue SW, King JR, Riley DS (2003) Extinction behaviour for two-dimensional inward-solidification problems. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. 459, 977–999.
[43] McCue SW, King JR, Riley DS (2005) The extinction problem for three-dimensional inward solidification. Journal of Engineering Mathematics 52, 389–409.
[44] McCue SW, Wu B, Hill JM (2008) Classical two-phase Stefan problem for spheres. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. 464, 2055–2076.
[45] El-Hachem M, McCue SW, Jin W, Du Y, Simpson MJ (2019) Revisiting the Fisher-Kolmogorov-Petrovsky-Piskunov equation to interpret the spreading-extinction dichotomy. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. 475: 20190378.
[46] Du Y, Lin Z (2010) Spreading–vanishing dichotomy in the diffusive logistic model with a free boundary. SIAM Journal on Mathematical Analysis. 42, 377–405.
[47] Du Y, Guo Z (2011) Spreading–vanishing dichotomy in a diffusive logistic model with a free boundary, II. Journal of Differential Equations. 250, 4336–4366.
[48] G Bunting, Y Du, K Krakowski (2012) Spreading speed revisited: Analysis of a free boundary model. Networks & Heterogeneous Media. 7, 583–603.
[49] Du Y, Guo Z (2012) The Stefan problem for the Fisher-KPP equation. Journal of Differential Equations. 253, 996–1035.
[50] Du Y, Matano H, Wang K. (2014) Regularity and asymptotic behavior of nonlinear Stefan problems. Archive for Rational Mechanics and Analysis. 212, 957–1010.
[51] Du Y, Matsuzawa H, Zhou M (2014) Sharp estimate of the spreading speed determined by nonlinear free boundary problems. SIAM Journal on Mathematical Analysis. 46, 375–396.
[52] Du Y, Lou B (2015) Spreading and vanishing in nonlinear diffusion problems with free boundaries. Journal of the European Mathematical Society. 17, 2673–2724.
[53] Hickson RI, Barry SI, Mercer GN (2009a) Critical times in multilayer diffusion. Part 1: Exact solutions. International Journal of Heat and Mass Transfer. 52, 5776-5783.
[54] Hickson RI, Barry SI, Mercer GN (2009b) Critical times in multilayer diffusion. Part 2: Approximate solutions. International Journal of Heat and Mass Transfer. 52, 5784-5791.
[55] Simpson MJ, Baker RE (2015) Exact calculations of survival probability for diffusion on growing lines, disks and spheres: The role of dimension. Journal of Chemical Physics. 143, 094109.
[56] Fadai NT, Simpson MJ (2020) New travelling wave solutions for the Porous-Fisher model with a moving boundary. Journal of Physics A: Mathematical and Theoretical. 53, 095601.
[57] Vázquez JL (2007) The porous medium equation. Oxford University Press, Oxford.
[58] Barenblatt GI (1996) Scaling, self-similarity, and intermediate asymptotics. Cambridge University Press, Cambridge.