Exact sharp-fronted travelling wave solutions of the Fisher-KPP equation
EExact sharp-fronted travelling wave solutions of the Fisher-KPPequation
Scott W. McCue , Maud El-Hachem , and Matthew J. Simpson ∗ School of Mathematical Sciences, Queensland University of Technology, Brisbane,Queensland 4001, AustraliaSeptember 3, 2020
Abstract
A family of travelling wave solutions to the Fisher-KPP equation with speeds c = ± / √ c = 5 / √
6, which decays to zero in the far-field, is exceptional in the sense that it can bewritten simply in terms of an exponential function. This solution has the property that thephase-plane trajectory is a heteroclinic orbit beginning at a saddle point and ends at theorigin. For c = − / √
6, there is also a trajectory that begins at the saddle point, but thissolution is normally disregarded as being unphysical as it blows up for finite z . We reinter-pret 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 advanc-ing into, an empty space. By simulating the full moving boundary problem numerically,we demonstrate how time-dependent solutions evolve to this exact travelling solution forlarge time. The relevance of such receding travelling waves to mathematical models for cellmigration and cell proliferation is also discussed.
Keywords:
Fisher-Kolmogorov; Painlev´e property; Weierstraß elliptic functions; Moving bound-ary problem ∗ To whom correspondence should be addressed. E-mail: [email protected] a r X i v : . [ n li n . PS ] S e p Introduction
For various applications in ecology and cell biology, the Fisher-KPP equation [1–3] ∂u∂t = ∂ u∂x + u (1 − u ) , (1)provides a very well studied model for the growth and spread of a population of species orcell types [4–6]. One key mathematical result is that, with the associated initial and boundaryconditions u ( x,
0) = F ( x ) , < x < ∞ , (2) ∂u∂x = 0 on x = 0 , u → x → ∞ , (3)the time-dependent solution evolves towards a travelling wave profile U ( z ), where z = x − ct ,as t → ∞ [1–3]. A combination of phase-plane analysis and simple asymptotics demonstratesthat the travelling wave speed c satisfies c ≥ F ( x ) in (2) [4–6].In this work, we are motivated by our recent studies [7–9] where we have restricted theFisher-KPP equation (1) 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 [10] ∂u∂t = ∂ u∂x + u (1 − u ) , < x < s ( t ) , (4) ∂u∂x = 0 on x = 0 , (5) u = 0 , d s d t = − κ ∂u∂x on x = s ( t ) . (6)Here κ is a parameter that relates the leakage of the population at the boundary to the speedof the boundary. In this context, we and others have provided new interpretations for travellingwave solutions (1) for c <
2, including slowly moving fronts that advance with speed 0 1) e z/ √ (cid:17) − , (7)as shown in Figure 1(a)-(b). Other exact travelling wave solutions to (1) for c = ± / √ 6, whichcan be written in terms of Weierstraß elliptic functions [15], are normally disregarded as beingunphysical in the usual way [4]. However, in the context of the Fisher-Stefan model (4)-(6), weprovide a new physical interpretation of one of these solutions. In particular, we claim that oneof the profiles for c = − / √ κ = − . . . . . In this way, we illustrate a second physically realistic exact travellingwave solution to (1) for c = ± / √ c = ± / √ 6, taken from Ablowitz & Zeppetella [15]. This derivation involves Weierstraß ellipticfunctions [16]. By using phase-plane analysis, we demonstrate the qualitative behaviour of thisfamily of solutions and highlight the two special trajectories that evolve to ( U, V ) = (1 , 0) as z → −∞ . In section 3 we link this special trajectory to solutions of the Fisher-Stefan modeland demonstrate numerically that for κ = 0 . . . . , initial conditions evolve to this travellingwave solution for large time. Finally, we providing concluding remarks in section 4. c = ± / √ To study travelling wave solutions of (1) we write u = U ( z ), where z = x − ct , to gived U d z + c d U d z + U (1 − U ) = 0 . (8)We discuss the domain of interest and the boundary conditions below, but for the moment wehighlight the physically relevant boundary condition U → − , d U d z → − , as z → −∞ , (9)which applies for two special cases considered below. We rewrite (8) in the usual way asd U d z = V, (10)d V d z = − cV − U (1 − U ) . (11)3 a)(c) (b)(d) Figure 1: Exact solutions to the Fisher-KPP model . (a)-(b) shows the exact solutionsfor c = 5 / √ c = − / √ U ( z ) that satisfy (9) (thick black), horizontal arrowsare superimposed to emphasize the difference in direction. (b),(d) show various phase planetrajectories. Equilibrium points in the phase planes are shown with black discs, and trajectoriescorresponding to the special travelling wave solutions are plotted (thick black). Various othertrajectories (thin coloured lines) are superimposed, and the corresponding curves are given in(a),(c) using the same colours as in (b),(d). Some trajectories are numbered to emphasise thepoint that the numbered trajectories in (a)-(d) are identical.4ne 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 (10)-(11) for c < U -axis of the phase-plane for c > c = ± / √ 6. For these values, we may solve (8)exactly, as explained by Ablowitz & Zeppetella [15]. A summary of the working is as follows.We start by letting U = f ( z ) w ( z ) and substituting into (8). Then, by forcing f (cid:48)(cid:48) + cf (cid:48) + f = 0,we find f w (cid:48)(cid:48) + (2 f (cid:48) + cf ) w (cid:48) = f g . We choose the linearly independent solution f = e λz , where λ = ( − c + √ c − / 2, so that w (cid:48)(cid:48) + √ c − w (cid:48) = e λz g . 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 (cid:48)(cid:48) + √ c − h (cid:48) = 0, which suggests we choose h = e −√ c − z . Finally, with c = 25 / w to be d w/ d s = 6 w which, upon multiplyingboth sides by d w/ d s , integrates directly to (cid:18) d w d s (cid:19) = 4 w − ω, (12)where ω is a constant. The first-order ode (12) is separable and is solved exactly in terms ofthe Weierstraß p-function ℘ ( z ; 0; ω ) [16]. Rewriting the solution in terms of U and V gives U = e − z/ √ ℘ (cid:16) e − z/ √ − k ; 0; ω (cid:17) , (13) V = − √ 63 e − z/ √ (cid:16) ℘ (cid:16) e − z/ √ − k ; 0; ω (cid:17) + 2e − z/ √ ℘ (cid:48) (cid:16) e − z/ √ − k ; 0; ω (cid:17)(cid:17) , (14)where ℘ (cid:48) is the derivative of ℘ [16]. It is noteworthy that this exact solution is possible because(8) has the Painlev´e property for c = ± / √ 6. In other words, for these special values of c , themovable singularities of solutions to (8) are poles [17].It is worth plotting the exact solution(s) (13)-(14) both in the form U = U ( z ) and in thephase-plane. As the travelling wave solutions are invariant to translations in z , we fix eachsolution in the z -direction by setting U = U at z = 0. For a given point in the phase plane,( U, V ) = ( U , V ), we determine the two constants k and ω by solving the nonlinear algebraicsystem U = ℘ (1 − k ; 0; ω ), V = − ( √ / 3) ( ℘ (1 − k ; 0; ω ) + 2 ℘ (cid:48) (1 − k ; 0; ω )) numerically, forexample with Newton’s method. The system can be sensitive and may require a close initialguess.Fixing c = 5 / √ U ( z ) are shown in Figure 1(a),while corresponding trajectories in the phase-plane are shown in Figure 1(b). Each of these5urves could be drawn using the exact solution (13)-(14) or just as easily be generated usingnumerical solutions to (10)-(11). As is well known, the phase-plane in Figure 1(b) is charac-terised by two fixed points, namely the saddle point at (1 , 0) and the stable node at (0 , , , U ∼ const e − z/ √ as z → ∞ . As z decreases, we see in both Figures 1(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 isthe heteroclinic orbit (thick black curve) that joins the two fixed points; this trajectory cor-responds to the well-known exact solution (7), which notably satisfies the physically realisticboundary condition (9). The simplification from (13)-(14) to (7) in this case arises because thisspecial case corresponds to taking the limit k → 0, which is in effect pushing the singularity to z = −∞ . Numerical solutions to (1)-(3) with appropriate initial conditions evolve to (7), as wedemonstrate in section 3.Now turning to c = − / √ 6, we show results in Figure 1(c)-(d). Here it is convenient toreflect our phase-plane about the U -axis, so the heteroclinic orbit just mentioned is in theupper-half plane. The trajectories in the phase-plane are still given by (13)-(14), except thatwe must make the changes V → − V , z → − z . Five of the solutions shown in Figure 1(c)-(d)are the same as in Figure 1(a)-(b) (to enable a straightforward comparison across the figures,we have labelled these solutions ○ – ○ ); these are trajectories that start at the stable node(0 , U ∼ const e z/ √ as z → −∞ . As we follow these trajectories forincreasing z , we again note that they blow up at a finite value of z (with U → −∞ ). Othertrajectories also blow up for finite z (this time with U → ∞ ), except for the separatrix (solidblack curve) which eventually enters the saddle point (1 , k and ω such that ℘ ( ζ ; 0 , ω ) has one of itsdouble poles at ζ = − k . This is a numerical constraint which can be enforced by solving thesystem U = ℘ (1 − k ; 0; ω ) , ℘ (cid:48) ( k ; 0; ω ) /℘ ( k ; 0; ω ) = 0 , (15)where 0 < U < 1. In Figure 1(c) we have chosen U = 0 . 2, but of course this value is arbitraryand, in effect, corresponds to a translation in z . The key observation of this special solutionis that it satisfies the physically realistic boundary condition (9) as z → −∞ . As we see inFigure 1(c), as z increases, this solution for U is very flat until it decreases sharply to U = 06 a) (b) Figure 2: Time-dependent PDE solutions . (a) Numerical solution of (1)-(3) showing theinitial condition (green), intermediate-time solutions at t = 3, 6 and 9 (blue). The solution at t = 9 is superimposed with (7) (dashed orange). (b) Numerical solution of (4)-(6) showing aninitial condition (green), intermediate-time solutions at t = 3, 6 and 9 (blue). The solution at t = 9 is superimposed with (13)-(14) subject to (15) (dashed orange).at some finite value of z = z ∗ where V = V ∗ = − . . . . , and then continues to decrease as z increases further. While this solution satisfies (9), it is normally disregarded as it is negativefor all z > z c . In the following section we show how this profile is a receding travelling wavesolution for the Fisher-Stefan moving boundary problem (4)-(6). Figure 2(a) shows time-dependent solutions of (1)-(3). Here, we see that a carefully-choseninitial condition with the appropriate decay at infinity evolves to a travelling wave solution thatis visually indistinguishable from (7). To make this point we show the initial condition in green,together some intermediate-time solutions in blue. The latest solution is superimposed with (7)in dashed orange, confirming that the time-dependent solutions converge to the exact solutionreasonably rapidly.Similar results in Figure 2(b) show time-dependent solutions of (4)-(6). Here we have used asimple step function for an initial condition, which is shown in green, and have carefully chosenour parameter κ to be κ = c/V ∗ = − . . . . . Again we show three intermediate-time solutionsin blue. The superimposed exact solution (13)-(14) with (15) for c = − / √ t = 9.7 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 ap-plications including wound healing [18–21] and ecological invasion [22, 23]. For cell biologyapplications, the Fisher-KPP model is often used because cells are thought to move randomly,by diffusion, as well as proliferating logistically [19, 24]. Experimental observations of mov-ing cell fronts can be described by travelling wave solutions of the Fisher-KPP model [24], orgeneralisations of the Fisher-KPP model [25]. For the dimensional Fisher-KPP model withdiffusivity D , proliferation rate λ and carrying capacity density K , the speed of the travellingwave solution is c ≥ √ λD [4]. In reality, fronts of cells may move at a slower speed or evenretreat [9]. One way to deal with this is to write down a more complicated model with morethan one species [26] or with a different source term, like an Allee effect [27]. Even with asingle species model that retains a logistic growth term, we can introduce nonlinear degeneratediffusion, leading to the Porous-Fisher model [28], which gives rise to travelling wave solutionswith c ≥ (cid:112) λD/ λ ; allows for solutions of all speeds c < √ λD , including negativespeeds [7].Despite the apparent simplicity of the Fisher-KPP model, exact solutions are relatively elu-sive but of high interest since they provide important mathematical insight and can be usedas benchmarks for testing numerical methods [31]. It is well-known that the travelling wavesolution (7) can be written down exactly for a special wave speed c = 5 / √ c > 2, whereas solutions with c < ∞ < c < 2, which is more flexible that the usual Fisher-KPP and Porous-Fishermodels since it can be used to model both invasion and retreat; and (iii) the Fisher-Stefanmodel provides a simple physical interpretation for an exact solution with c = − / √ 6. This8verlooked solution can be expressed exactly using Weierstraß elliptic functions and, as we shownumerically, this solution is the long-time limit of our moving boundary problem (4)-(6).9 Appendix To construct the phase planes in the main document we solve the first order dynamical systemd U d z = V, (16)d V d z = − cV − U (1 − U ) , (17)numerically, using a two-stage Runge-Kutta method, also called Heun’s method [34], with aconstant step size, d z = 1 × − . This choice of discretisation leads to results that are visuallyindependent of the numerical discretisation. Our main interest is in examining phase planetrajectories that either enter or leave the saddle (1 , 0) along the stable or unstable manifold, re-spectively. Therefore, it is important that the initial condition we chose when solving Equations(16)–(17) are on the appropriate stable or unstable manifold, and sufficiently close to (1 , eig function [32] to calculate the eigenvalues and eigen-vectors when c = ± / √ 6. The vector field associated with Equations (16)–(17) are plotted onthe phase planes using the MATLAB quiver function [33]. MATLAB implementations of thesealgorithms are available on GitHub. Results in Figure 1 (main document) show regions of the phase plane that are deliberatelyscaled so that we can easily see the details of the trajectories associated with c = ± / √ 6. Sinceeach phase plane in Figure 1 (main document) is shown on a different scale, it is easy to forgetthat these phase planes are closely related, with the only difference being the sign of c . Tomake this point we show here, in Figure 3, the phase plane shown over a wider domain with c = ± / √ 6. We show, in thick red lines, the trajectories associated with the travelling waveswith c = ± / √ 6. For completeness, we also show many other phase plane trajectories, in thinblack lines. These additional trajectories are not associated with travelling waves. Note that thethick red trajectory in the third and fourth quadrant are associated with the invading travellingwave with c = 5 / √ 6, whereas the thick red trajectory in the first quadrant is associated withthe receding travelling wave with c = − / √ 6. To plot these trajectories on the same phaseplane, the trajectory associated with the receding travelling wave c = − / √ z = − × − .Figure 3: Phase plane visualisation . Phase plane showing the two equilibria (black discs).The two trajectories corresponding to the travelling wave solutions (thick red lines) and a rangeof other trajectories (thin black lines) are shown. We solve the Fisher–KPP equation ∂u∂t = ∂ u∂x + u (1 − u ) , (18)on 0 ≤ x ≤ L , with L chosen to be sufficiently large so that the time–dependent solutionsof the partial differential equation (PDE) model have sufficient time to approach a travellingwave. We discretise the domain with a uniform finite difference mesh, with spacing ∆ x , and weapproximate the spatial derivatives in Equation (18) using a central finite difference approxi-mation. The resulting system of coupled ordinary differential equations are integrated through11ime using an implicit Euler approximation, giving rise to u j +1 i − u ji ∆ t = (cid:32) u j +1 i − − u j +1 i + u j +1 i +1 ∆ x (cid:33) + u j +1 i (1 − u j +1 i ) , (19)for i = 2 , . . . , m − 1, where m = L/ ∆ x + 1 is the total number of spatial nodes, and j indexestime so that we approximate u ( x, t ) by u ji , where x = ( i − x and t = j ∆ t .For all numerical solutions of Equation (18) we enforce no–flux boundary conditions at x = 0and x = L , giving u j +12 − u j +11 = 0 , (20) u j +1 m − u j +1 m − = 0 . (21)In this work we wish to study the solution of Equation (18) with a travelling wave speedgreater than the minimum speed, c > 2. Therefore, care is taken to choose the initial conditionto achieve this. For the initial condition we set u ( x, 0) = 12 0 ≤ x < , e − √ ( x − ≤ x < L, (22)where the exponential decay rate of the initial condition is carefully chosen to be u ( x, ∼ e − x/ √ , because this choice leads to travelling waves with c = 5 / √ t to t + ∆ t we solve the system of nonlinear alge-braic system, Equations (19)–(21), 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 (18) are obtained by setting ∆ x = 1 × − , ∆ t = 1 × − , L = 60 and ε = 1 × − . We find that these choices lead to results that are independent of thenumerical discretisation. Using the numerically–generated time–dependent solutions we alsoestimate the travelling wave speed. To estimate c ∗ we specify a contour value, u ( x, t ) = u ∗ , andat the end of each time step, we use linear interpolation to find x ∗ such that u ( x ∗ , t ) = u ∗ , andwe then calculate c ∗ = [ x ∗ ( t + ∆ t ) − x ∗ ( t )] / ∆ t . MATLAB implementations of these algorithmsare available on GitHub. 12 .4 Numerical method for the Fisher–Stefan model We solve ∂u∂t = ∂ u∂x + u (1 − u ) , (23)for 0 < x < s ( t ) numerically by using a boundary fixing transformation ξ = x/s ( t ) so that wehave ∂u∂t = 1 s ( t ) ∂ u∂ξ + ξs ( t ) d s ( t )d t ∂u∂ξ + u (1 − u ) , (24)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 , (25) u = 0 at ξ = 1 . (26)Spatially discretising Equations (24)–(26) with a uniform finite difference mesh, with spac-ing ∆ ξ , approximating the spatial derivatives using central differences and an implicit Eulerapproximation for the temporal derivative gives, u j +1 i − u ji ∆ t = 1( s j ) (cid:32) u j +1 i − − u j +1 i + u j +1 i +1 ∆ ξ (cid:33) + ξs j (cid:18) s j +1 − s j ∆ t (cid:19) (cid:32) u j +1 i +1 − u j +1 i − ξ (cid:33) + u j +1 i (1 − u j +1 i ) , (27)for i = 2 , . . . , m − 1, where m = 1 / ∆ ξ + 1 is the total number of spatial nodes on the finitedifference mesh, and the index j represents the time index so that we approximate u ( ξ, t ) by u ji , where ξ = ( i − 1) ∆ ξ and t = j ∆ t .Discretising the boundary conditions, Equations (25)–(26), gives u j +12 − u j +11 = 0 , (28) u j +1 m = 0 . (29)The initial condition for the Fisher–Stefan problem is u ( x, 0) = 1 on 0 ≤ x ≤ s (0). Toadvance the discrete system from time t to t + ∆ t we solve the system of nonlinear algebraicequations, Equations (27)–(29), using Newton–Raphson iteration with convergence tolerance ε .13uring each iteration of the Newton–Raphson algorithm we update s ( t ) using the discretisedStefan condition s j +1 = s j − κ ∆ ts j ∆ ξ (cid:32) u j +1 m − − u j +1 m − + 3 u j +1 m (cid:33) , (30)where we set u j +1 m = 0 to be consistent with (29).All results in this work are obtained by setting κ = − . . . . , ∆ ξ = 1 × − , ∆ t = 1 × − , ε = 1 × − and s (0) = 100. Again, we find these choices lead to results that are independent ofthe numerical discretisation, and we also use the time–dependent PDE solution to estimate c ∗ bytracking the time evolution of the leading edge, c ∗ = ( s j +1 − s j ) / ∆ t . MATLAB implementationsof these algorithms are available on GitHub. Acknowledgements: This work is supported by the Australian Research Council (DP200100177).14 eferences [1] RA Fisher. The wave of advance of advantageous genes. Annals of Eugenics. 7 (1937)355-369.[2] AN Kolmogorov, PG Petrovskii, NS Piskunov. A study of the diffusion equation withincrease in the amount of substance, and its application to a biological problem. MoscowUniversity Mathematics Bulletin. 1 (1937) 1-26.[3] J Canosa. On a nonlinear diffusion equation describing population growth. IBM Journal ofResearch and Development. 17 (1973) 307-313.[4] JD Murray. Mathematical Biology I: An Introduction. Third edition, Springer, New York,(2002).[5] L Edelstein-Keshet. Mathematical Models in Biology. SIAM, Philadelphia, (2005).[6] M Kot. Elements of Mathematical Ecology. Camridge University Press, Cambridge, (2003).[7] M El–Hachem, SW McCue, W Jin, Y Du, MJ Simpson. Revisiting the Fisher-Kolmogorov-Petrovsky-Piskunov equation to interpret the spreading–extinction dichotomy. Proceed-ings of the Royal Society A: Mathematical, Physical and Engineering Sciences. 475 (2019)20190378.[8] M El–Hachem, SW McCue, MJ Simpson. A sharp-front moving boundary model for ma-lignant invasion. Physica D: Nonlinear Phenomena. 412, 132639.[9] M El–Hachem, SW McCue, MJ Simpson. Invading and receding sharp-fronted travellingwaves. arXiv:2008.00662[10] Y Du, Z Lin. Spreading–vanishing dichotomy in the diffusive logistic model with a freeboundary. SIAM Journal on Mathematical Analysis. 42 (2010) 377-405.[11] Y Du, H Matano, K Wang. Regularity and asymptotic behavior of nonlinear Stefan prob-lems. Archive for Rational Mechanics and Analysis 212 (2014) 957-1010.[12] Y Du, H Matsuzawa, M Zhou. Sharp estimate of the spreading speed determined by non-linear free boundary problems. SIAM Journal on Mathematical Analysis 46 (2014) 375-396.[13] Y Du, B Lou. Spreading and vanishing in nonlinear diffusion problems with free boundaries.Journal of the European Mathematical Society 17 (2015) 2673-2724.1514] P Kaliappan. An exact solution for travelling waves of u t = Du xx + u − u k . Physica D:Nonlinear Phenomena. 11 (1983) 368-374.[15] M Ablowitz, A Zeppetella. Explicit solutions of Fisher’s equation for a special wave speed.Bulletin of Mathematical Biology. 41 (1979) 835-840.[16] M Abramowitz, IA Stegun (Eds). Handbook of Mathematical Functions with Formulas,Graphs, and Mathematical Tables. Dover Publications, New York (1965).[17] PA Clarkson. Painlev´e Equations — Nonlinear Special Functions. In: Marcell´an F., VanAssche W. (eds) Orthogonal Polynomials and Special Functions. Lecture Notes in Mathe-matics, vol 1883. Springer, Berlin, Heidelberg (2006)[18] BG Sengers, CP Please, ROC Oreffo. Experimental characterization and computationalmodelling of two–dimensional cell spreading for skeletal regeneration. Journal of the RoyalSociety Interface. 4 (2007) 1107-1117.[19] JA Sherratt, JD Murray. Models of epidermal wound healing. Proceedings of the RoyalSociety of London: Series B. 241 (1990) 29-36.[20] ST Johnston, JV Ross, BJ Binder, DLS McElwain, P Haridas, MJ Simpson. Quantifyingthe effect of experimental design choices for in vitro scratch assays. Journal of TheoreticalBiology. 400 (2016) 19-31.[21] ST Vittadello, SW McCue, G Gunasingh, NK Haass, MJ Simpson. Mathematical modelsfor cell migration with real–time cell cycle dynamics. Biophysical Journal. 114 (2018) 1241-1253.[22] JG Skellam. Random dispersal in theoretical populations. Biometrika. 38 (1951) 196-218.[23] J Steele, J Adams, T Sluckin. Modelling paleoindian dispersals. World Archaeology. 30(1998) 286–305.[24] PK Maini, DLS McElwain, D Leavesley. Traveling waves in a wound healing assay. AppliedMathematics Letters. 17 (2004) 575-580.[25] DJ Warne, RE Baker, MJ Simpson. Using experimental data and information criteria toguide model selection for reaction–diffusion problems in mathematical biology. Bulletin ofMathematical Biology. 81 (2019) 1760-1804.[26] KJ Painter, JA Sherratt. Modelling the movement of interacting cell populations. Journalof Theoretical Biology. 225 (2003) 327-339.1627] NT Fadai, MJ Simpson. Population dynamics with threshold effects gives rise to a diversefamily of Allee effects. Bulletin of Mathematical Biology. 82 (2020) 74.[28] TP Witelski. Merging traveling waves for the porous–Fisher’s equation. Applied Mathe-matics Letters. 8 (1995) 57-62.[29] F S´anchez Garduno, PK Maini. An approximation to a sharp type solution of a density–dependent reaction–diffusion equation. Applied Mathematics Letters. 7 (1994) 47-51.[30] JA Sherratt, BP Marchant. Nonsharp travelling wave fronts in the Fisher equation withdegenerate nonlinear diffusion. Applied Mathematics Letters. 9 (1996) 33-38.[31] MJ Simpson, KA Landman. Characterizing and minimizing the operator spit error forFisher’s equation. Applied Mathematics Letters. 19 (2006) 604-612.[32] MathWorks eig . Retrieved August 2020 from https://au.mathworks.com.[33] MathWorks quiverquiver