title: Exact sharp-fronted travelling wave solutions of the Fisher-KPP equation
authors: Scott W McCue Maud El-Hachem Matthew J Simpson
text: 
  
  Cosmic expansion with matter creation and bulk viscosity
      Samuel Lepe^3
      December 30, 2023
  ========================================================
  
  
  A family of travelling wave solutions to the Fisher-KPP equation with speeds c=± 5/√(6) can be expressed exactly using Weierstraß elliptic functions.  The well-known solution for c=5/√(6), which decays to zero in the far-field, is exceptional in the sense that it can be written simply in terms of an exponential function.  This solution has the property that the phase-plane trajectory is a heteroclinic orbit beginning at a saddle point and ends at the origin.  For c=-5/√(6), there is also a trajectory that begins at the saddle point, but this solution is normally disregarded as being unphysical as it blows up for finite z.  We reinterpret this special trajectory as an exact sharp-fronted travelling solution to a Fisher-Stefan type moving boundary problem, where the population is receding from, instead of advancing into, an empty space.  By simulating the full moving boundary problem numerically, we demonstrate how time-dependent solutions evolve to this exact travelling solution for large time.  The relevance of such receding travelling waves to mathematical models for cell migration and cell proliferation is also discussed.
  
  
  
  
    
  Keywords: Fisher-Kolmogorov; Painlevé property;  Weierstraß elliptic functions; Moving boundary problem
  
  
  
  § INTRODUCTION
   
  
  For various applications in ecology and cell biology, the Fisher-KPP equation <cit.>
  
      ∂ u/∂ t=∂^2 u/∂ x^2+u(1-u),
  
  provides a very well studied model for the growth and spread of a population of species or cell types <cit.>.  One key mathematical result is that, with the associated initial and boundary conditions
  
      u(x,0)=F(x),    0<x<∞,
  
  
      ∂ u/∂ x=0      x=0,
         u→ 0     x→∞,
  
  the time-dependent solution evolves towards a travelling wave profile U(z), where z=x-ct, as t→∞ <cit.>.  A combination of phase-plane analysis and simple asymptotics demonstrates that the travelling wave speed c satisfies c≥ 2 and is selected by the far-field behaviour of the initial condition F(x) in (<ref>) <cit.>.
  
  In this work, we are motivated by our recent studies <cit.> where we have restricted the Fisher-KPP equation (<ref>) to hold on the moving domain 0<x<s(t), together with a Stefan-type moving boundary condition, to give the so-called Fisher-Stefan model <cit.>
  
      ∂ u/∂ t=∂^2 u/∂ x^2+u(1-u),
         0<x<s(t),
  
  
      ∂ u/∂ x=0      x=0,
  
  
      u=0,   ds/dt=-κ∂ u/∂ x     x=s(t).
  
  Here κ is a parameter that relates the leakage of the population at the boundary to the speed of the boundary.  In this context, we and others have provided new interpretations for travelling wave solutions (<ref>) for c<2, including slowly moving fronts that advance with speed 0<c<2 <cit.>, stationary profiles for c=0 and receding fronts with speed c<0 <cit.>.  Travelling wave solutions for c<2 are interesting because they are normally disregarded as being unphysical (because they do not satisfy the boundary conditions and/or are not restricted to 0<U<1 for all z∈ℝ) <cit.>.
  
  In this letter we focus on travelling wave solutions to (<ref>) for the special values c=± 5/√(6).  For c=5/√(6) there is a well known exact solution <cit.>
  
      U=(1+(√(5)-1) e^z/√(6))^-2,
  
  as shown in Figure <ref>(a)-(b).  Other exact travelling wave solutions to (<ref>) for c=± 5/√(6), which can be written in terms of Weierstraß elliptic functions <cit.>, are normally disregarded as being unphysical in the usual way <cit.>.  However, in the context of the Fisher-Stefan model (<ref>)-(<ref>), we provide a new physical interpretation of one of these solutions. In particular, we claim that one of the profiles for c=-5/√(6) corresponds to a receding travelling wave to (<ref>)-(<ref>) with a special value of κ=-0.906….  In this way, we illustrate a second physically realistic exact travelling wave solution to (<ref>) for c=± 5/√(6).
  
  
  
  In section <ref> we review the exact travelling solutions to the Fisher-KPP equation for c=± 5/√(6), taken from Ablowitz & Zeppetella <cit.>.  This derivation involves Weierstraß elliptic functions <cit.>.  By using phase-plane analysis, we demonstrate the qualitative behaviour of this family of solutions and highlight the two special trajectories that evolve to (U,V)=(1,0) as z→ -∞.  In section <ref> we link this special trajectory to solutions of the Fisher-Stefan model and demonstrate numerically that for κ=0.906…, initial conditions evolve to this travelling wave solution for large time. Finally, we providing concluding remarks in section <ref>.
  
  
  
  § TRAVELLING WAVE SOLUTIONS FOR C=± 5/√(6)
  
  
  To study travelling wave solutions of (<ref>) we write u=U(z), where z=x-ct, to give
  
      d^2U/dz^2+cdU/dz+U(1-U)=0.
  
  We discuss the domain of interest and the boundary conditions below, but for the moment we highlight the physically relevant boundary condition
  
      U→ 1^-,   dU/dz→ 0^-,
           z→ -∞,
  
  which applies for two special cases considered below.  We rewrite (<ref>) in the usual way as
  
      dU/d z   = V, 
      dV/d z   = -cV - U(1-U).
  
  One point to note here is that this system is reversible under the substitution z→ -z, V→ -V, c→ -c, which means that the phase-plane for (<ref>)-(<ref>) for c<0 is simply a reflection about the U-axis of the phase-plane for c>0.
  
  The concern here is with special values c=± 5/√(6).  For these values, we may solve (<ref>) exactly, as explained by Ablowitz & Zeppetella <cit.>.  A summary of the working is as follows.  We start by letting U=f(z)w(z) and substituting into (<ref>).  Then, by forcing f”+cf'+f=0, we find fw”+(2f'+cf)w'=f^2g^2.  We choose the linearly independent solution f=e^λ z, where λ=(-c+√(c^2-4))/2, so that w”+√(c^2-4)w'=e^λ zg^2.  The equations simplify by setting w=w(s), s=h(z).  The left-hand side of this differential equation reduces to a single term if h”+√(c^2-4)h'=0, which suggests we choose h=e^-√(c^2-4)z.  Finally, with c^2=25/6 we end up with the second-order differential equation for w to be d^2w/ds^2=6w which, upon multiplying both sides by dw/ds, integrates directly to
  
      (dw/ds)^2=4w^3-ω,
  
  where ω is a constant.  The first-order ode (<ref>) is separable and is solved exactly in terms of the Weierstraß p-function ℘(z;0;ω) <cit.>.  Rewriting the solution in terms of U and V gives
  
      U    = e^-2z/√(6)℘(e^-z/√(6)-k;0;ω), 
      
      V    =  -√(6)/3e^-2z/√(6)(
      ℘(e^-z/√(6)-k;0;ω)
      +2e^-z/√(6)℘'(e^-z/√(6)-k;0;ω)
      ),
  
  where ℘' is the derivative of ℘ <cit.>.  It is noteworthy that this exact solution is possible because (<ref>) has the Painlevé property for c=± 5/√(6).  In other words, for these special values of c, the movable singularities of solutions to (<ref>) are poles <cit.>.
  
  It is worth plotting the exact solution(s) (<ref>)-(<ref>) both in the form U=U(z) and in the phase-plane.  As the travelling wave solutions are invariant to translations in z, we fix each solution in the z-direction by setting U=U_0 at z=0.  For a given point in the phase plane, (U,V)=(U_0,V_0), we determine the two constants k and ω by solving the nonlinear algebraic system U_0=℘(1-k;0;ω), V_0=-(√(6)/3)(℘(1-k;0;ω)+2℘'(1-k;0;ω)
  ) numerically, for example with Newton's method.  The system can be sensitive and may require a close initial guess.
  
  Fixing c= 5/√(6) for the moment, travelling wave profiles U(z) are shown in Figure <ref>(a), while corresponding trajectories in the phase-plane are shown in Figure <ref>(b).  Each of these curves could be drawn using the exact solution (<ref>)-(<ref>) or just as easily be generated using numerical solutions to (<ref>)-(<ref>).  As is well known, the phase-plane in Figure <ref>(b) is characterised by two fixed points, namely the saddle point at (1,0) and the stable node at (0,0).  All of the trajectories in Figure <ref>(b) enter the stable node (0,0); a linearisation about (0,0) demonstrates that U∼ e^-2z/√(6) as z→∞.  As z decreases, we see in both Figures <ref>(a) and (b) that, with one exception, each solution (thin coloured curves) blows up (with U→±∞) at a finite value of z.  This is because ℘(ζ;0,ω) has a double pole at ζ=0.  The exception is the heteroclinic orbit (thick black curve) that joins the two fixed points; this trajectory corresponds to the well-known exact solution (<ref>), which notably satisfies the physically realistic boundary condition (<ref>).  The simplification from (<ref>)-(<ref>) to (<ref>) in this case arises because this special case corresponds to taking the limit k→ 0, which is in effect pushing the singularity to z=-∞.  Numerical solutions to (<ref>)-(<ref>) with appropriate initial conditions evolve to (<ref>), as we demonstrate in section <ref>.
  
  Now turning to c= -5/√(6), we show results in Figure <ref>(c)-(d).  Here it is convenient to reflect our phase-plane about the U-axis, so the heteroclinic orbit just mentioned is in the upper-half plane.  The trajectories in the phase-plane are still given by (<ref>)-(<ref>), except that we must make the changes V→ -V, z→ -z.  Five of the solutions shown in Figure <ref>(c)-(d) are the same as in Figure <ref>(a)-(b) (to enable a straightforward comparison across the figures, we have labelled these solutions 1–5); these are trajectories that start at the stable node (0,0),  except now we see that U∼ e^2z/√(6) as z→ -∞.  As we follow these trajectories for increasing z, we again note that they blow up at a finite value of z (with U→ -∞).  Other trajectories also blow up for finite z (this time with U→∞), except for the separatrix (solid black curve) which eventually enters the saddle point (1,0).
  
  We now interpret the separatrix in Figure <ref>(d) in terms of an exact sharp-fronted travelling wave solution of (<ref>).  As explained in Ablowitz & Zeppetella <cit.>, in order to extract this special case from (<ref>)-(<ref>), we must choose the constants k and ω such that ℘(ζ; 0,ω) has one of its double poles at ζ=-k. This is a numerical constraint which can be enforced by solving the system
  
      U_0=℘(1-k;0;ω),
        ℘'(k;0;ω)/℘(k;0;ω)^2=0,
  
  where 0<U_0<1.  In Figure <ref>(c) we have chosen U_0=0.2, but of course this value is arbitrary and, in effect, corresponds to a translation in z. The key observation of this special solution is that it satisfies the physically realistic boundary condition (<ref>) as z→ -∞.  As we see in Figure <ref>(c), as z increases, this solution for U is very flat until it decreases sharply to U=0 at some finite value of z=z^* where V=V^*=-2.251…, and then continues to decrease as z increases further.  While this solution satisfies (<ref>), it is normally disregarded as it is negative for all z>z_c.  In the following section we show how this profile is a receding travelling wave solution for the Fisher-Stefan moving boundary problem (<ref>)-(<ref>).
  
  
  
  
  
  
  
  § TIME-DEPENDENT SOLUTIONS TO THE FISHER-STEFAN MODEL
  
  
  Figure <ref>(a) shows time-dependent solutions of (<ref>)-(<ref>). Here, we see that a carefully-chosen initial condition with the appropriate decay at infinity evolves to a travelling wave solution that is visually indistinguishable from  (<ref>).  To make this point we show the initial condition in green, together some intermediate-time solutions in blue.  The latest solution is superimposed with (<ref>) in dashed orange, confirming that the time-dependent solutions converge to the exact solution reasonably rapidly.
  
  
  
  
  Similar results in Figure <ref>(b) show time-dependent solutions of (<ref>)-(<ref>).  Here we have used a simple step function for an initial condition, which is shown in green, and have carefully chosen our parameter κ to be κ =  c/V^*=-0.906….  Again we show three intermediate-time solutions in blue.  The superimposed exact solution (<ref>)-(<ref>) with (<ref>) for c=-5/√(6) in orange compares extremely well with the numerical solution at t=9.
  
  
  
  
  
  
  
  § CONCLUSION
   
  
  We present a new interpretation of an exact travelling wave solution of the Fisher-KPP model.  The Fisher-KPP model is one of the most well-studied reaction-diffusion equations with applications including wound healing <cit.> and ecological invasion <cit.>.  For cell biology applications, the Fisher-KPP model is often used because cells are thought to move randomly, by diffusion, as well as proliferating logistically <cit.>.  Experimental observations of moving cell fronts can be described by travelling wave solutions of the Fisher-KPP model <cit.>, or generalisations of the Fisher-KPP model <cit.>.  For the dimensional Fisher-KPP model with diffusivity D, proliferation rate λ and carrying capacity density K,  the speed of the travelling wave solution is c≥ 2√(λ D) <cit.>.  In reality, fronts of cells may move at a slower speed or even retreat <cit.>.  One way to deal with this is to write down a more complicated model with more than one species <cit.> or with a different source term, like an Allee effect <cit.>.  Even with a single species model that retains a logistic growth term, we can introduce nonlinear degenerate diffusion, leading to the Porous-Fisher model <cit.>, which gives rise to travelling wave solutions with c≥√(λ D/2) <cit.>.   An alternative to all of these modifications is to simply retain Fisher (<ref>) but include a moving boundary as in (<ref>)-(<ref>) <cit.>.  This has appealing features, namely: a sharp moving front, which seems biologically reasonable; retains the classical logistic growth term with an easy to measure λ; allows for solutions of all speeds c<2√(λ D), including negative speeds <cit.>.
  
  Despite the apparent simplicity of the Fisher-KPP model, exact solutions are relatively elusive but of high interest since they provide important mathematical insight and can be used as benchmarks for testing numerical methods <cit.>.  It is well-known that the travelling wave solution (<ref>) can be written down exactly for a special wave speed c=5/√(6) <cit.>.  This special travelling wave is consistent with the usual view that travelling wave solutions of the nondimensional Fisher-KPP have positive speed  c > 2, whereas solutions with c< 2 are normally disregarded on the grounds of being unphysical <cit.>.  In our work we take a different point of view and re-formulate the Fisher-KPP model with a moving boundary, often called the Fisher-Stefan model <cit.>.  The Fisher-Stefan model has several attractive features: (i) travelling wave solutions of the Fisher-Stefan model have a well-defined front without needing to introduce the complication degenerate nonlinear diffusion; (ii) the Fisher-Stefan model gives rise to travelling wave solutions with ∞ < c< 2, which is more flexible that the usual Fisher-KPP and Porous-Fisher models since it can be used to model both invasion and retreat; and (iii)  the Fisher-Stefan model provides a simple physical interpretation for an exact solution with c=-5/√(6).  This overlooked  solution can be expressed exactly using Weierstraß elliptic functions and, as we show numerically, this solution is the long-time limit of our moving boundary problem (<ref>)-(<ref>). 
  
  
  
  
  
  
  
  
  
  § APPENDIX 
   
  
  
   §.§ Numerical method for the phase plane
   
  To construct the phase planes in the main document we solve the first order dynamical system
  
      dU/d z    = V, 
      dV/d z   = -cV - U(1-U),
  
  numerically, using a two-stage Runge-Kutta method, also called Heun's method <cit.>, with a constant step size, dz = 1 × 10^-3.  This choice of discretisation leads to results that are visually independent of the numerical discretisation.  Our main interest is in examining phase plane trajectories that either enter or leave the saddle (1,0) along the stable or unstable manifold, respectively.  Therefore, it is important that the initial condition we chose when solving Equations (<ref>)–(<ref>) are on the appropriate stable or unstable manifold, and sufficiently close to (1,0).  To choose this point we use the MATLAB eig function <cit.> to calculate the eigenvalues and eigenvectors when c=± 5/√(6). The vector field associated with  Equations (<ref>)–(<ref>) are plotted on the phase planes using the MATLAB quiver function <cit.>.  MATLAB implementations of these algorithms are available on https://github.com/maudelhachem/McCue2020GitHub.
  
  
  
  
   §.§ Phase plane visualisation
   
  Results in Figure 1 (main document) show regions of the phase plane that are deliberately scaled so that we can easily see the details of the trajectories associated with c = ± 5 / √(6).  Since each phase plane in Figure 1 (main document) is shown on a different scale, it is easy to forget that these phase planes are closely related, with the only difference being the sign of c.  To make this point we show here, in Figure <ref>, the phase plane shown over a wider domain with c = ± 5 / √(6).  We show, in thick red lines, the trajectories associated with the travelling waves with c = ± 5 / √(6).  For completeness, we also show many other phase plane trajectories, in thin black lines.  These additional trajectories are not associated with travelling waves.  Note that the thick red trajectory in the third and fourth quadrant are associated with the invading travelling wave with c=5/√(6), whereas the thick red trajectory in the first quadrant is associated with the receding travelling wave with c=-5/√(6).  To plot these trajectories on the same phase plane, the trajectory associated with the receding travelling wave c=-5/√(6) is obtained by integrating the dynamical system in the other direction, effectively by setting dz = -1 × 10^-3.
  
  
  
  
  
  
   §.§ Numerical method for the Fisher–KPP model
  
  We solve the Fisher–KPP equation
  
      ∂ u/∂ t =∂^2 u/∂ x^2 +  u(1-u),
  
  on 0 ≤ x ≤ L, with L chosen to be sufficiently large so that the time–dependent solutions of the partial differential equation (PDE) model have sufficient time to approach a travelling wave.  We discretise the domain with a uniform finite difference mesh, with spacing Δ x, and we approximate the spatial derivatives in Equation (<ref>) using a central finite difference approximation.  The resulting system of coupled ordinary differential equations are integrated through time using an implicit Euler approximation, giving rise to
  
      u_i^j+1 - u_i^j/Δ t = ( u_i-1^j+1- 2 u_i^j+1 + u_i+1^j+1/Δ x^2) +  u_i^j+1 (1-u_i^j+1),
  
  for i = 2, …, m-1, where m = L/Δ x + 1 is the total number of spatial nodes, and j indexes time so that we approximate u(x,t) by u_i^j, where x=(i-1)Δ x and t=jΔ t.
  
  For all numerical solutions of Equation (<ref>) we enforce no–flux boundary conditions at x=0 and x=L, giving
  
      u_2^j+1-u_1^j+1 = 0,  
         u_m^j+1- u_m-1^j+1 = 0.
  
  
  
  In this work we wish to study the solution of Equation (<ref>) with a travelling wave speed greater than the minimum speed, c > 2.  Therefore, care is taken to choose the initial condition to achieve this.  For the initial condition we set
  
      u(x,0) =
      12      0 ≤ x < 10, 
      
      e^-2√(6)(x-10)2      10 ≤ x < L,
  
  where the exponential decay rate of the initial condition is carefully chosen to be u(x,0) ∼e^-2x/√(6), because this choice leads to travelling waves with c=5/√(6) <cit.>.
  
  To advance the discrete system from time t to t + Δ t we solve the system of nonlinear algebraic system, Equations (<ref>)–(<ref>), using Newton–Raphson iteration with convergence tolerance ε. The resulting systems of linear equations are solved efficiently using the Thomas algorithm. All numerical solutions of  Equation (<ref>) are obtained by setting Δ x = 1 × 10^-4, Δ t = 1 × 10^-2, L=60 and ε= 1 × 10^-8.  We find that these choices lead to results that are independent of the numerical discretisation.   Using the numerically–generated time–dependent solutions we  also estimate the travelling wave speed.  To estimate c^* we specify a contour value, u(x,t) = u^*, and at the end of each time step, we use linear interpolation to find x^* such that u(x^*,t) = u^*, and we then calculatec^* = [x^*(t + Δ t)-x^*(t)]/Δ t. MATLAB implementations of these algorithms are available on https://github.com/maudelhachem/McCue2020GitHub.
  
  
  
   §.§ Numerical method for the Fisher–Stefan model
  
  We solve
  
      ∂ u/∂ t =∂^2 u/∂ x^2 +  u(1-u),
  
  for 0 < x < s(t) numerically by using a boundary fixing transformation ξ = x / s(t) so that we have
  
      ∂ u/∂ t = 1/s^2(t)∂^2 u/∂ξ^2+ξ/s(t)ds(t)/d t∂ u/∂ξ + u(1-u),
  
  on the fixed domain, 0 < ξ < 1.  Here s(t) is the length of the domain that we will discuss later.  To close the problem we also transform the appropriate boundary conditions
  
      ∂ u∂ξ = 0   at  ξ=0,  
         u = 0   at  ξ=1.
  
  
  Spatially discretising Equations (<ref>)–(<ref>) with a uniform finite difference mesh, with spacing Δξ, approximating the spatial derivatives using central differences and an implicit Euler approximation for the temporal derivative gives,
  
      u_i^j+1 - u_i^jΔ t   =  1(s^j)^2( u_i-1^j+1 - 2 u_i^j+1 + u_i+1^j+1/Δξ^2)+ ξs^j(s^j+1 - s^jΔ t) ( u_i+1^j+1 - u_i-1^j+1/2 Δξ) +  u_i^j+1(1 - u_i^j+1) ,
  
  for i = 2, …, m-1, where m = 1/Δξ + 1 is the total number of spatial nodes on the finite difference mesh, and the index j represents the time index so that we approximate u(ξ,t) by u_i^j, where ξ = (i - 1) Δ ξ and t=jΔ t.
  
  Discretising the boundary conditions, Equations (<ref>)–(<ref>), gives
  
      u_2^j+1-u_1^j+1 = 0,  
         u_m^j+1 = 0.
  
  
  The initial condition for the Fisher–Stefan problem is u(x,0)=1 on 0 ≤ x ≤ s(0). To advance the discrete system from time t to t + Δ t we solve the system of nonlinear algebraic equations, Equations (<ref>)–(<ref>), using Newton–Raphson iteration with convergence tolerance ε. During each iteration of the Newton–Raphson algorithm we update s(t) using the discretised Stefan condition
  
      s^j+1 = s^j - κΔ ts^jΔξ(u_m-2^j+12-2u_m-1^j+1+ 3u_m^j+12),
  
  where we set  u_m^j+1 = 0 to be consistent with (<ref>).
  
  All results in this work are obtained by setting κ=-0.906…, Δξ = 1 × 10^-5, Δ t = 1 × 10^-2, ε = 1 × 10^-8 and s(0)=100. Again, we find these choices lead to results that are independent of the numerical discretisation, and we also use the time–dependent PDE solution to estimate c^* by tracking the time evolution of the leading edge, c^* = (s^j+1 - s^j)/Δ t. MATLAB implementations of these algorithms are available on https://github.com/maudelhachem/McCue2020GitHub.
  
  
  
  
  
    
  Acknowledgements:This work is supported by the Australian Research Council (DP200100177).
  
  99
  
  Fisher1937 RA Fisher. The wave of advance of advantageous genes. Annals of Eugenics. 7 (1937) 355-369.
  
  Kolmogorov1937 AN Kolmogorov, PG Petrovskii, NS Piskunov. 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 (1937) 1-26.
  
  Canosa1973 J Canosa. On a nonlinear diffusion equation describing population growth. IBM Journal of Research and Development. 17 (1973) 307-313.
  
  Murray02 JD Murray. Mathematical Biology I: An Introduction. Third edition, Springer, New York, (2002).
  
  Edelstein2005 L Edelstein-Keshet. Mathematical Models in Biology. SIAM, Philadelphia,  (2005).
  
  Kot03 M Kot. Elements of Mathematical Ecology. Camridge University Press, Cambridge, (2003).
  
  Elhachem2019 M El–Hachem, SW McCue, W Jin, Y Du, MJ Simpson. 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 (2019) 20190378.
  
  Elhachem2020a M El–Hachem, SW McCue, MJ Simpson. A sharp-front moving boundary model for malignant invasion. Physica D: Nonlinear Phenomena. 412, 132639.
  
  Elhachem2020b M El–Hachem, SW McCue, MJ Simpson. Invading and receding sharp-fronted travelling waves.  https://arxiv.org/abs/2008.00662arXiv:2008.00662
  
  Du2010 Y Du, Z Lin. Spreading–vanishing dichotomy in the diffusive logistic model with a free boundary. SIAM Journal on Mathematical Analysis. 42 (2010) 377-405.
  
  Du2014a Y Du, H Matano, K Wang. Regularity and asymptotic behavior of nonlinear Stefan problems. Archive for Rational Mechanics and Analysis 212 (2014) 957-1010.
  
  Du2014b Y Du, H Matsuzawa, M Zhou. Sharp estimate of the spreading speed determined by nonlinear free boundary problems. SIAM Journal on Mathematical Analysis 46 (2014) 375-396.
  
  Du2015 Y Du,  B Lou. Spreading and vanishing in nonlinear diffusion problems with free boundaries. Journal of the European Mathematical Society 17 (2015) 2673-2724.
  
  Kaliappan1984 P Kaliappan. An exact solution for travelling waves of u_t=Du_xx+u - u^k.  Physica D: Nonlinear Phenomena. 11 (1983) 368-374.
  
  Ablowitz1979 M Ablowitz, A Zeppetella. Explicit solutions of Fisher's equation for a special wave speed. Bulletin of Mathematical Biology. 41 (1979) 835-840.
  
  Abramowitz1965 M Abramowitz, IA Stegun (Eds). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover Publications, New York (1965).
  
  Clarkson2006 PA Clarkson.  Painlevé Equations — Nonlinear Special Functions.  In: Marcellán F., Van Assche W. (eds) Orthogonal Polynomials and Special Functions. Lecture Notes in Mathematics, vol 1883. Springer, Berlin, Heidelberg (2006)
  
  Sengers2007 BG Sengers, CP Please, ROC Oreffo. Experimental characterization and computational modelling of two–dimensional cell spreading for skeletal regeneration. Journal of the Royal Society Interface.  4 (2007) 1107-1117.
  
  Sherratt90 JA Sherratt, JD Murray. Models of epidermal wound healing. Proceedings of the Royal Society of London: Series B. 241 (1990) 29-36.
  
  Johston2015 ST Johnston, JV Ross, BJ Binder, DLS McElwain, P Haridas, MJ Simpson. Quantifying the effect of experimental design choices for in vitro scratch assays.  Journal of Theoretical Biology. 400 (2016) 19-31.
  
  Vittadello2018 ST Vittadello, SW McCue, G Gunasingh, NK Haass, MJ Simpson. Mathematical models for cell migration with real–time cell cycle dynamics.  Biophysical Journal.  114 (2018) 1241-1253.
  
  Skellam1951 JG Skellam. Random dispersal in theoretical populations. Biometrika. 38 (1951) 196-218.
  
  Steele1998  J Steele, J Adams, T Sluckin. Modelling paleoindian dispersals. World Archaeology. 30 (1998) 286–305.
  
  Maini2004 PK Maini, DLS McElwain, D Leavesley. Traveling waves in a wound healing assay. Applied Mathematics Letters. 17 (2004) 575-580.
  
  Warne2019 DJ Warne, RE Baker, MJ Simpson. Using experimental data and information criteria to guide model selection for reaction–diffusion problems in mathematical biology. Bulletin of Mathematical Biology. 81 (2019) 1760-1804.
  
  Painter2003 KJ Painter, JA Sherratt. Modelling the movement of interacting cell populations. Journal of Theoretical Biology. 225 (2003) 327-339.
  
  Fadai2020 NT Fadai, MJ Simpson.   Population dynamics with threshold effects gives rise to a diverse family of Allee effects. Bulletin of Mathematical Biology. 82 (2020) 74.
  
  Witelski1995 TP Witelski. Merging traveling waves for the porous–Fisher's equation. Applied Mathematics Letters. 8 (1995) 57-62.
  
  Sanchez1994 F Sánchez Garduno, PK Maini. An approximation to a sharp type solution of a density–dependent reaction–diffusion equation. Applied Mathematics Letters. 7 (1994) 47-51.
  
  Sherratt1996 JA Sherratt, BP Marchant. Nonsharp travelling wave fronts in the Fisher equation with degenerate nonlinear diffusion.  Applied Mathematics Letters. 9 (1996) 33-38.
  
  
  Simpson2006 MJ Simpson, KA Landman. Characterizing and minimizing the operator spit error for Fisher’s equation.  Applied Mathematics Letters. 19 (2006) 604-612.
  
  
  eig MathWorks eig. Retrieved August 2020 from  https://au.mathworks.com/help/matlab/ref/eig.html https://au.mathworks.com.
  
  quiver MathWorks quiver. Retrieved August 2020 from  https://www.mathworks.com/help/matlab/ref/quiver.html https://www.mathworks.com.
  
  
  Simpson2005 M.J. Simpson, K.A. Landman, T.P. Clement. Assessment of a non-traditional operator split algorithm for simulation of reactive transport. Mathematics and Computers in Simulation. 70 (2005) 44-60.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  