Invading and receding sharp-fronted travelling waves
IInvading and receding sharp–fronted travelling waves
Maud El-Hachem , Scott W. McCue , and Matthew J. Simpson ∗ School of Mathematical Sciences, Queensland University of Technology,Brisbane, Queensland 4001, AustraliaAugust 11, 2020
Abstract
Biological invasion, whereby populations of motile and proliferative individualslead to moving fronts that invade into vacant regions, are routinely studied usingpartial differential equation (PDE) models based upon the classical Fisher–KPPmodel. While the Fisher–KPP model and extensions have been successfully usedto model a range of invasive phenomena, including ecological and cellular invasion,an often–overlooked limitation of the Fisher–KPP model is that it cannot be usedto model biological recession where the spatial extent of the population decreaseswith time. In this work we study the
Fisher–Stefan model, which is a generalisationof the Fisher–KPP model obtained by reformulating the Fisher–KPP model as amoving boundary problem. The nondimensional Fisher–Stefan model involves justone single parameter, κ , which relates the shape of the density front at the movingboundary to the speed of the associated travelling wave, c . Using numerical simu-lation, phase plane and perturbation analysis, we construct approximate solutionsof the Fisher–Stefan model for both slowly invading and slowly receding travellingwaves, as well as for rapidly receding travelling waves. These approximations al-low us to determine the relationship between c and κ so that commonly–reported ∗ To whom correspondence should be addressed. E-mail: [email protected] a r X i v : . [ q - b i o . T O ] A ug xperimental estimates of c can be used to provide estimates of the unknown pa-rameter κ . Interestingly, when we reinterpret the Fisher–KPP model as a movingboundary problem, many disregarded features of the classical Fisher–KPP phaseplane take on a new interpretation since travelling waves solutions with c < Keywords:
Invasion; Reaction–diffusion; Partial differential equation; Stefan problem;Moving boundary problem. 2
Introduction
Biological invasion is normally associated with situations where individuals within a popu-lation undergo both movement and proliferation events [Murray 2002]. Such proliferationand movement, combined, can give rise to an invading front . An invading front involvesa population moving into a previously unoccupied space. Ecologists are particularly in-terested in biological invasion. For example, Skellam’s [Skellam 1951] work studies theinvasion of muskrats in Europe; similarly, Otto and coworkers [Otto et al. 2018] studythe spatial spreading of insects, whereas Bate and Hilker [Bate and Hilker 2019] studythe invasion of predators in a predator–prey system. These three studies all make use ofpartial differential equation (PDE) models of invasion.Another common application of biological invasion is the study of cell invasion, in-cluding wound healing and malignant spreading. Mathematical models of wound healingoften consider the closure of a wound space by populations of cells that are both mi-gratory and proliferative [Flegg et al. 2020, Jin et al. 2016, Maini et al. 2004, Sherrattand Murray 1990]. Malignant invasion involves combined migration and proliferation oftumour cells, which leads to tumour invasion into surrounding tissues [Byrne 2010,Curtinet al. 2020,Strobl et al. 2020,Swanson et al. 2003], as illustrated in Figure 1(a)–(b) whichshows the invasion of malignant melanoma cells into surrounding tissues. Regardless ofthe application, many continuum mathematical models of biological invasion involve thestudy of moving fronts, as shown schematically in Figure 1(c), using PDE models [Brown-ing et al. 2019, Sengers et al. 2007, Warne et al. 2019]. We interpret the schematic inFigure 1(c) by thinking of the population as being composed of individuals that undergodiffusive migration with diffusivity
D >
0, and carrying capacity–limited proliferation,such as the classical logistic model, with proliferation rate λ >
0. As indicated, these twoprocesses can lead to the spatial expansion as the population density profile moves in thepositive x –direction.The Fisher–KPP model [Canosa 1973, Fisher 1937, Kolmogorov et al. 1937, Murray2002] is probably the most commonly used reaction–diffusion equation to describe biolog-ical invasion in a single homogeneous population. The Fisher–KPP model assumes thatindividuals in the population proliferate logistically and move according to a simple linear3 a) (b) den s i t y x (c) invasionrecessionfront front f r on t Figure 1:
Biological motivation and schematic . (a) Malignant melanoma (dark)spreading superficially across the skin surface [NCI 1985] (reproduced with permission).(b) Vertical cross section through a human skin equivalent experimental model show-ing the inward invasion of a population of melanoma cells (dark) [Haridas et al. 2017](reproduced with permission). In (a)–(b) the region containing the leading edge of theinvading population is highlighted in a red rectangle and the location of the sharp front ishighlighted with blue arrows. (c) Schematic solution of a mathematical model showing asharp–fronted density profile that could either invade or recede, by moving in the positiveor negative x –direction, respectively. In the schematic the location of the sharp front isalso highlighted with a blue arrow.diffusion mechanism [Canosa 1973, Fisher 1937, Kolmogorov et al. 1937, Murray 2002].Travelling wave solutions of the Fisher–KPP model are often used to mimic biologicalinvasion [Maini et al. 2004, Maini et al. 2004b, Simpson et al. 2013]. Long time solutionsof the Fisher–KPP model that evolve from initial conditions with compact support leadto invasion waves that move with speed c = 2 √ λD [Murray 2002]. There are manyother popular choices of single–species mathematical models of biological invasion. Forexample, the Porous–Fisher model [McCue et al. 2019, Sanchez et al. 1995, Sherratt etal. 1996, Witelski 1995] is a generalisation of the Fisher–KPP model with a degeneratenonlinear diffusion term which results in sharp–fronted travelling wave solutions. Longtime solutions of the Porous–Fisher model that evolve from initial conditions with com-pact support lead to invasion waves that move with speed c = (cid:112) ( λD ) / Fisher–Stefan model [Du and Lin2010]. This approach involves reformulating the Fisher–KPP model as a moving bound-ary problem on 0 < x < L ( t ). Setting the density to zero at the moving front, x = L ( t ),means that the Fisher–Stefan model also gives rise to sharp–fronted solutions like thePorous–Fisher model [El–Hachem et al. 2019]. The motion of L ( t ) in the Fisher–Stefan4odel is given by a one–phase Stefan condition [Crank 1987, Dalwadi et al. 2020, Hill1987, Mitchell and O’Brien 2014].Populations of motile and proliferative individuals do not always invade new territory,in fact, sometimes motile and proliferative populations recede or retreat. The spatial re-cession of biological populations are often described in ecology. For example, populationsof desert locusts [Ibrahim et al. 2000], plants in grazed prairies [Sinkins and Otfinowski2012], Arctic foxes [Killengreen et al. 2007] and dung beetles [Horgan 2009] have allbeen observed to undergo both invasion and recession in different circumstances. Whilesome previous mathematical models of biological invasion and recession have been de-scribed [Chaplain et al. 2020, El–Hachem et al. 2020, Painter and Sherratt 2003, Simpsonet al. 2009], these previous models often focus on describing interactions between multiplesubpopulations in a heterogeneous community rather than classical single species models,such as the Fisher–KPP model. In fact, none of the three commonly–used single speciesmodels described here, the Fisher–KPP, Porous–Fisher or Fisher–Stefan models, havebeen used to study biological recession. This is probably because neither the classicalFisher–KPP or Porous–Fisher models ever give rise to receding populations. Given thatthe recession of population fronts is often observed, this limitation of the commonly–usedFisher–KPP and Porous–Fisher models is important and often overlooked.In this work we focus on the Fisher–Stefan model to study biological invasion andrecession. Unlike the classical Fisher–KPP and Porous–Fisher models, the Fisher–Stefanmodel can be used to simulate both biological invasion and recession. One way of inter-preting this difference is that the Fisher–Stefan model could be thought of as being morerealistic than the more commonly–used Fisher–KPP or Porous–Fisher models. As wewill show, travelling wave solutions of the Fisher–Stefan model can be used to representbiological invasion with a positive travelling wave speed, c >
0, as well as being able tomodel biological recession with a negative travelling wave speed, c <
0. We explore thesetravelling wave solutions using full time–dependent numerical solutions of the governingPDE, phase plane analysis, and perturbation approximations. A regular perturbationapproximation around c = 0 provides insight into both slowly invading and recedingtravelling waves, whereas a matched asymptotic expansion in the limit as c → −∞ pro-vides insight into rapidly receding waves. These perturbation solutions provide simple5elationships between κ and c . These relationships are useful because estimates of κ arenot available in the literature, whereas experimental measurements of c are relativelystraightforward to obtain [Maini et al. 2004, Maini et al. 2004b, Simpson et al. 2007]. In this work all dimensional variables and parameters are denoted with a circumflex, andnondimensional quantities are denoted using regular symbols. The Fisher–Stefan modelis a reformulation of the classical Fisher–KPP equation to include a moving boundary, ∂ ˆ u∂ ˆ t = ˆ D ∂ ˆ u∂ ˆ x + ˆ λ ˆ u (cid:18) − ˆ u ˆ K (cid:19) , < ˆ x < ˆ L (ˆ t ) , (1)where ˆ u (ˆ x, ˆ t ) ≥ x , and time,ˆ t >
0. Individuals in the population move according to a linear diffusion mechanism withdiffusivity ˆ
D >
0, the proliferation rate is ˆ λ >
K > < ˆ x < ˆ L (ˆ t ), with a zero flux condition atthe origin and the sharp front is modelled by setting the density to be zero at the leadingedge, giving ∂ ˆ u (0 , t ) ∂ ˆ x = 0 , ˆ u ( ˆ L (ˆ t ) , ˆ t ) = 0 . (2)The evolution of the domain is given by a classical one–phase Stefan condition that relatesthe speed of the moving front to the spatial gradient of the density profile at the movingboundary, d ˆ L (ˆ t )dˆ t = − ˆ κ ∂ ˆ u∂ ˆ x (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ L (ˆ t ) , (3)where ˆ κ is a constant to be specified [Crank 1987, Hill 1987].In the context of cell invasion, typical values of ˆ D are approximately 100–3000 µ m /h [John-ston et al. 2015, Johnston et al. 2016, Jin et al. 2016]; typical values ˆ λ are approximately0.04–0.06 /h [Johnston et al. 2015,Johnston et al. 2016,Jin et al. 2016]; and typical valuesof the carrying capacity density are 0.001–0.003 cells/ µ m [Johnston et al. 2015,Johnston6t al. 2016, Jin et al. 2016]. To simplify our analysis we will now nondimensionalise theFisher–Stefan model. Introducing dimensionless variables, x = ˆ x (cid:113) ˆ λ/ ˆ D , t = ˆ λ ˆ t , u = ˆ u/ ˆ K , L ( t ) = ˆ L (ˆ t ) (cid:113) ˆ λ/ ˆ D and κ = ˆ κ/ ˆ D , the Fisher–Stefan model can be simplified to give ∂u∂t = ∂ u∂x + u (1 − u ) , < x < L ( t ) , (4) ∂u (0 , t ) ∂x = 0 , u ( L ( t ) , t ) = 0 , (5)d L ( t )d t = − κ ∂u ( L ( t ) , t ) ∂x , (6)so that we only need to specify one parameter, κ . As mentioned previously, estimatesof diffusivity, proliferation rate and carrying capacity in the context of cell invasion areavailable in the literature [Jin et al. 2016, Maini et al. 2004]. In contrast, estimates of κ are not. Therefore, one of the aims of this work is to provide mathematical insight intohow estimates of κ can be obtained, and we will provide more discussion on this pointlater.In all cases where we consider time dependent solutions of Equations (4)–(6) we alwayschoose the initial condition to be u ( x,
0) = α (1 − H[ L (0)]) , (7)where α > · ] is the Heaviside function, so that u ( x,
0) = α for x < L (0) and u ( x, L (0)) = 0.To solve Equations (4)–(7) numerically, we transform the governing equations froman evolving domain, 0 < x < L ( t ), to a fixed domain, 0 < ξ <
1, by setting ξ = x/L ( t ). The transformed equations on the fixed domain are spatially discretised using auniform finite difference mesh and standard central finite difference approximations. Theresulting system of nonlinear ordinary differential equations (ODE) are integrated throughtime using an implicit Euler approximation. Newton–Raphson iteration and the Thomas7lgorithm are used to solve the resulting system of nonlinear algebraic equations [Simpsonet al. 2005]. Full details of the numerical method are given in the Appendix, and aMATLAB implementation of the algorithm is available on GitHub. We begin our analysis of the Fisher–Stefan model by presenting some time–dependentsolutions of Equations (4)–(7) before analysing these solutions using the phase plane andperturbation techniques.
Results in Figure 2 show a suite of numerical solutions of Equations (4)–(7) plotteda regular time intervals. These results suggest that the initial condition evolves intoconstant speed invading travelling waves in Figure 2(a)–(d) with κ >
0, and recedingtravelling waves in Figure 2(e)–(h) with κ <
0. To obtain these solutions we specifya value of κ , as indicated in each subfigure, and then measure the eventual speed ofthe travelling wave, c , by estimating d L ( t ) / d t using the numerical solution of the PDE.Therefore, in this approach to studying the travelling wave solutions we think of κ as aninput to the numerical algorithm, and c is an output from the numerical algorithm. Infact, in generating results in Figure 2 we took great care to choose κ carefully so thatour estimates of c are very simple, clean values such as c = 0 . , . , .
75 and 1.00. Atthis point it is not obvious how to choose κ to give such estimates and we will explainhow to make such choices later, in Section 3.2. Results in Figure 2 show that c is anincreasing function of κ , and the density profile at the leading edge is sharp in all cases.The shape of the density profile differs depending on whether we consider an invadingor receding travelling wave, since the receding travelling waves are much steeper thanthe invading travelling waves. These numerical results in Figure 2 are interesting sinceneither the Fisher–KPP nor the Porous–Fisher can be used to simulate this range ofbehaviours. The feature of the Fisher–Stefan model which enables us to simulate bothinvasion and retreat is the choice of κ . We will now explore the relationship between c κ by studying the travelling wave solutions in the phase plane.9 a) (e)(b) (f)(c) (g)(d) (h) Figure 2:
Invading and receding travelling wave solutions of the Fisher–Stefanmodel.
Numerical solutions of Equations (4)–(7) are given at t = 0 , , ,
30. Theinitial condition is given by Equation (7) with α = 0 . L (0) = 200. Results in(a)–(d) lead to invading travelling waves with c = 0 . , . , .
75 and 1 .
00, respectively.These travelling waves are obtained by choosing κ = 0 . , . , . . c = − . , − . , − .
00 and − .
99, respectively. These receding travelling waves are obtainedby choosing κ = − . , − . , − . − . c correspond are obtained at late time, t = 40. Note that estimates of κ are reportedhere in the caption to four decimal places, whereas the estimates given in the subfiguresare reported to two decimal places to keep the figure neat.10 .2 Phase plane analysis To analyse travelling wave solutions of the Fisher–Stefan model in the phase plane weconsider Equation (4) in terms of the travelling wave coordinate, z = x − ct and we seeksolutions of the form u ( x, t ) = U ( z ) which leads to the following ODE,d U d z + c d U d z + U (1 − U ) = 0 , −∞ < z < , (8)with boundary conditions U ( −∞ ) = 1 , U (0) = 0 , (9) c = − κ d U (0)d z , (10)where we choose z = 0 to correspond to the moving boundary.To study Equation (8) in the phase plane we rewrite this second order ODE as a firstorder dynamical system d U d z = V, (11)d V d z = − cV − U (1 − U ) , (12)where the equilibrium points are (0 ,
0) and (1 , ,
0) is a saddle point for all values of c , whereas (0 ,
0) is a stable node if c ≥ < c <
2; a centre if c = 0; an unstable spiral if − < c < c ≤ −
2. Typically, in the regular analysis of the Fisher–KPPmodel the possibility of travelling wave solutions with c < c < ,
0) is a stable spiral, implying that U ( z ) < κ as the input and c as the output of the numerical algorithm, here in the phaseplane we think of c as the input into the numerical algorithm to generate the phase planetrajectory and we use this numerically–generated trajectory to obtain the estimate of κ ,as we will now explain. Phase planes for c = 0 . , . , .
75 and 1 .
00 are given in Figure3 (a)–(d), respectively. Similarly, phase planes for c = − . , − . , − .
00 and − .
99 aregiven in Figure 3 (e)–(f), respectively. Each phase plane in Figure 3 corresponds to theparticular PDE solution in Figure 2. 12 a) (e)(b) (f)(c) (g)(d) (h)
Figure 3:
Phase plane for travelling wave solutions of the Fisher–Stefan model.
Equilibrium points are shown as black discs, and the point at which the trajectory in-tersects the V ( z ) axis are shown as pink discs. The numerical solution of the dynamicalsystem, Equations (11)–(12) is shown in dashed orange and the travelling wave solutionobtained from the numerical time–dependent PDE solutions, Equations (4)–(7) is super-imposed in solid purple for the invading travelling waves in (a)–(d) and in solid green forthe receding travelling waves in (e)–(h). The flow associated with the dynamical systemis shown with blue vectors obtained using Matlab’s quiver function.13he phase planes in Figure 3(a)–(d) correspond to invading fronts with various valuesof 0 < c <
2. As we previously describe, these trajectories in the phase plane are normallyneglected in the usual analysis of the Fisher–KPP model since these trajectories leave(1 ,
0) and eventually spiral into (0 ,
0) as z → ∞ , implying that U ( z ) < U ( z ) = 0, which means that we truncate the trajectoryat z = 0 and we only focus on that part of the trajectory in the fourth quadrant of thephase plane where U ( z ) >
0. Each trajectory in Figure 3(a)–(d) intersects the V ( z ) axisat a special point, (0 , V ∗ ), and this intersection point corresponds to the Stefan conditionwhere U = 0 and c = − κV ∗ . Estimating V ∗ from the numerical phase plane trajectorythen allows us to estimate κ . Following this approach we obtain estimates of κ for eachvalue of c , and these estimates compare very well with the estimates used to generate thetime–dependent PDE solutions in Figure 2. These phase planes explain why invadingtravelling waves for the Fisher–Stefan model are restricted to 0 < c < c > ,
0) and (0 ,
0) neverintersects the V ( z ) axis so that as c → − we have κ → ∞ [El–Hachem et al. 2019].For completeness we also show the remaining portion of the phase plane trajectoryin Figure 3(a)–(d) that eventually spirals into (0 ,
0) as z → ∞ . Further, for each phaseplane in Figure 3(a)–(d) we take the late time PDE solution from Figure 2(a)–(d) andtransform these PDE solutions into a ( U ( z ) , V ( z )) trajectory in the phase plane, and plotthese curves in the phase planes in Figure 3(a)–(d). In each case the trajectory obtainedby solving the dynamical system numerically is visually indistinguishable, at this scale,from the trajectory obtained by plotting the PDE solutions in the phase plane.The phase planes in Figure 3(e)–(h) correspond to receding fronts with various valuesof c <
0. As we previously described, these phase planes for c < , V ∗ ) and joins(1 ,
0) as z → ∞ . Again, we can use this trajectory to estimate κ and the estimatesfrom the phase plane compare well with the values used in the full time–dependent PDEsolutions in Figure 2(e)–(h). For completeness we take the late–time PDE solutions in14igure 2(e)–(h) and superimpose these trajectories in Figure 3(e)–(h) where we see thatthe numerical solution of the trajectory obtained from the dynamical system is againvisually indistinguishable from the trajectory obtained from the PDE solutions. Unlikethe invading travelling wave solutions where linear stability analysis in the phase planegives us the condition that 0 < c <
2, there is no restriction on c in the phase plane sothat the Fisher–Stefan model gives rise to receding travelling waves with −∞ < c < κ and c . This will be important because estimates of κ are not available in the literature,whereas estimates of c are easier to obtain experimentally [Maini et al. 2004, Maini et al.2004b, Simpson et al. 2007]. Here we aim to establish approximate solutions of the travelling wave solutions when | c | (cid:28)
1, and we begin by re–writing Equations (11)–(12) asd V d U = − cV − U (1 − U ) V , (13)where it is clear that an exact solution for V ( U ) can be obtained when c = 0. Thereforewe seek an perturbation solution [Murray 1984] of the form V ( U ) = V ( U ) + cV ( U ) + c V ( U ) + O ( c ) . (14)15ubstituting Equation (14) into Equation (13) gives,d V d U V + U (1 − U ) = 0 , V (1) = 0 , (15)d V d U V + d V d U V + V = 0 , V (1) = 0 , (16)d V d U V + d V d U V + V (cid:18) d V d U + 1 (cid:19) = 0 , V (1) = 0 . (17)The solutions of these differential equations are V ( U ) = (cid:112) U + 1)3 ( U − , (18) V ( U ) = − ( U − U ) / − √ U − √ U , (19) V ( U ) = − √ U + 1) / ( U − √ U + 3 − ( √ U + 3 + 3) × (cid:32) − U (6 U − U + 20) + 15 U ( U + 2) + 31+ √ U + 3 [(2 U + 1)(6 U + 3) − U − U − U + 30) ln (cid:20) ( √ U + 3 + 3)( U − √ U + 3 − (cid:21)(cid:19) , (20)and Maple code to generate these solutions are available on GitHub. These three solu-tions can be used to truncate Equation (14) at different orders, and in doing so we willmake use O (1), O ( c ) and O ( c ) perturbation solutions. Given our various approximateperturbation solutions for V ( U ), we can either directly plot these solution in the phaseplane and compare them with numerically–generated phase plane trajectories, or we canintegrate these perturbation solutions numerically to give an approximation for the shapeof the travelling wave, U ( z ). To estimate the shape of the travelling wave we integratethe perturbation solution for V ( U ) using Heun’s method with U (0) = 0, and we integratefrom z = 0 to z = − Z , where Z is taken to be sufficiently large.We now compare these various perturbation solutions with phase plane trajectoriesand time–dependent PDE solutions for both invading and receding travelling waves. Fig-ure 4 focuses on invading travelling wave with various c >
0. Results in Figure 4(a)–(c)show the phase plane for c = 0 . , .
50 and 0.75, respectively. The numerical solutionof the dynamical system is shown in green, and is superimposed on the O ( c ) and O ( c )16erturbation solutions in yellow and blue, respectively. In these results there is a visualdifference between the numerically–generated phase plane trajectories and the O ( c ) per-turbation solutions, however the O ( c ) perturbation solution compares very well with thenumerically–generated phase plane trajectories.17 a) (b) Figure 4:
Perturbation solutions for slow invading travelling waves. (a)–(c) show the phase plane for c = 0 . , .
50 and 0.75,respectively. Equilibrium points are shown with black discs. The numerical solution of Equations (11)–(12) are shown in green andthe point at which these trajectories intersect the V ( z ) axis are shown with a green disc. The O ( c ) and O ( c ) perturbation solutionsare shown in yellow and blue, respectively. The intersection of the V ( z ) for the O ( c ) and O ( c ) perturbation solutions are shown in ayellow and blue disc, respectively. Results in (d)–(f) compare the shape of the travelling wave profile, U ( z ), obtained using the numericalsolution of the phase plane trajectory (green) with the O ( c ) and O ( c ) perturbation solutions in yellow and blue, respectively. Results in(g)–(i) show magnified comparison of the three solutions in the regions highlighted by the dashed boxes in (d)–(f). esults in Figure 4(d)–(f) compare the shape of the travelling wave, U ( z ), using thenumerical solution of the dynamical system in the phase plane with the results obtainedfrom the O ( c ) and O ( c ) perturbation solutions. For the numerical solution of the dynam-ical system we deliberately show the invasion profile using the trajectory from z = − z = 5, which includes the unphysical part of the trajectory, z >
0, where U ( z ) isoscillatory. To make a clear distinction between the physical and unphysical parts ofthe invading profile we include a horizontal line at U ( z ) = 0 which makes it clear that U ( z ) > z <
0, and U ( z ) is oscillatory for z >
0. All three solutions are visuallyindistinguishable at the scale shown in Figure 4(d) where c = 0 .
25. For c = 0 .
50 and c = 0 .
75 we see a visually–distinct difference between the profiles from the phase planetrajectory and the O ( c ) perturbation solutions, whereas the O ( c ) perturbation solutiongives an excellent approximation for these larger speeds. Results in Figure 4(g)–(i) showmagnified comparisons of the shape of U ( z ) corresponding to the dashed inset regions inFigure 4(d)–(f) where it is easier to see the distinction between the three solutions.Results in Figure 5 for the receding travelling wave are presented in the exact sameformat as those in Figure 4. Here, in Figure 5 we consider c = − . , − .
75 and -1.00 andwe see that the O ( c ) perturbation solution provides a very accurate approximation ofboth the phase plane trajectory and the shape of the receding travelling wave.19 a) (f)(b) (e)(c) (d) (g)(h)(i) Figure 5:
Perturbation solutions for slow receding travelling waves. (a)–(c) show the phase plane for c = − . , − .
75 and-1.00, respectively. Equilibrium points are shown with black discs. The numerical solution of Equations (11)–(12) are shown in greenand the point at which these trajectories intersect the V ( z ) axis are shown with a green disc. The O ( c ) and O ( c ) perturbation solutionsare shown in yellow and blue, respectively. The intersection of the V ( z ) for the O ( c ) and O ( c ) perturbation solutions are shown in ayellow and blue disc, respectively. Results in (d)–(f) compare the shape of the travelling wave profile, U ( z ), obtained using the numericalsolution of the phase plane trajectory (green) with the O ( c ) and O ( c ) perturbation solutions in yellow and blue, respectively. Results in(g)–(i) show magnified comparison of the three solutions in the regions highlighted by the dashed boxes in (d)–(f). s we pointed out previously, one of the key conceptual limitations of using theFisher–Stefan model is that, unlike applications in physical and material sciences [Crank1987, Hill 1987], estimates of κ are not available and rarely discussed for biological ap-plications [Gaffney and Maini 1999]. One way to address this limitation is to use ouranalysis to provide simple relationships between κ and c , since the wave speed relativelystraightforward to measure [Maini et al. 2004, Maini et al. 2004b, Simpson et al. 2007].As we noted previously, all travelling wave solutions of the Fisher–Stefan model satisfy κ = − c/V (0), where V = V ( U ). When | c | (cid:28) V (0) using our pertur-bation solutions and this will provide various relationships between κ and c dependingon the order of the perturbation solution used for V (0), O (1) : κ = − cV (0) , (21) O ( c ) : κ = − cV (0) + cV (0) , (22) O ( c ) : κ = − cV (0) + cV (0) + c V (0) . (23)Now substituting expressions for V (0), V (0) and V (0) and expanding the resultingexpressions for | c | (cid:28) O (1) : κ ( c ) = √ c + O ( c ) , (24) O ( c ) : κ ( c ) = √ c −
35 (2 − √ c + O ( c ) , (25) O ( c ) : κ ( c ) = √ c −
35 (2 − √ c − √ (cid:20)
10 ln (cid:18)
62 + √ (cid:19) + 12 √ − (cid:21) c + O ( c ) , (26)which provides a simple way to relate c and κ for | c | (cid:28)
1. To explore the accuracy ofthese approximations we use numerical solutions in the phase plane to estimate κ in theinterval − < c < c and κ inFigure 6. We also superimpose the various approximations, given by Equations (24)–(26)in Figure 6, where we see that Equation (26) is particularly accurate for | c | (cid:28) . a) (f)(b) (e)(c) (d) (g)(h)(i) Figure 6:
Relationship between c and κ for | c | (cid:28) . The numerical estimate of κ asa function of c is given in green. Various perturbation approximations given by Equation(24)–(26) are given in red, yellow and blue, respectively. The various relationships be-tween c and κ are shown in two insets. The first inset, for − . < c < − .
1, is outlinedin black. The second inset, for 0 . < c < .
3, is outlined in pink.
As we noted in Section 3.3.1, Equation (15) can be solved exactly when c = 0, corre-sponding to a stationary wave. This solution can be written as V ( U ) = ± (cid:114) − U + 2 U + 13 , (27)where, we are primarily interested in the negative solution since V < U (0) = 0 can be integrated to give the shape of the stationary22ave, U ( z ) = 32 tanh (cid:32) z − arctanh √ (cid:33) − . (28)Results in Figure 7 compare these exact solutions for c = 0 with various numericalsolutions. Firstly, in Figure 7(a) we show a time–dependent solution of Equations (4)–(7)with κ = 0 which evolves into a stationary wave that is visually indistinguishable fromthe exact solution, Equation (28), at this scale. The phase plane in Figure 7(b) showsthe late–time PDE solution from Figure 7(a) plotted as a trajectory in the ( U ( z ) , V ( z ))phase plane. In this phase plane we superimpose the exact solution, Equation (27), whichforms a homoclinic orbit in the shape of a teardrop. The part of the homoclinic orbit inthe fourth quadrant of the phase plane corresponds to the stationary wave and we seethat the numerical trajectory and the exact solution are indistinguishable at this scale.23igure 7: Exact solution for the stationary travelling wave, c = 0 . (a) Com-parison of the exact solution, Equation (28), in dashed blue with the numerical solutionof Equations (4)–(7) with κ = 0 in solid green. The initial condition for the numericalsolution of the PDE is in orange. (b) Comparison of the exact solution of the phase planetrajectory, Equation (27), in dashed blue, with the trajectory obtained by plotting thePDE solution in the phase plane in solid green. Equilibrium points in the phase planeare shown with black discs. As we noted previously in Section 3.1, preliminary numerical simulations of travellingwaves with c < c d U d z + d U d z + 1 c U (1 − U ) = 0 , −∞ < z < , (29)24hich is singular as c → −∞ . Therefore, we will construct a matched asymptotic expan-sion [Murray 1984] by treating ε = 1 /c as a small parameter. The boundary conditionsfor this problem are U (0) = 0 and U ( z ) = 1 as z → −∞ . Setting ε = 0 and solving theresulting ODE gives the outer solution, U ( z ) = 1 , (30)which matches the boundary condition as z → −∞ . To construct the inner solution near z = 0 we rescale the independent variable ζ = z/(cid:15) , so that we have ζ = zc . Therefore, inthe boundary layer we haved U d ζ + d U d ζ + 1 c U (1 − U ) = 0 , −∞ < ζ < . (31)Now expanding U ( ζ ) in a series we obtain U ( ζ ) = U ( ζ ) + 1 c U ( ζ ) + 1 c U ( ζ ) + O (cid:18) c (cid:19) , (32)which we substitute into Equation (31) to give a family of boundary value problems,d U d ζ + d U d ζ = 0 , U (0) = 0 , U → ζ → −∞ , (33)d U d ζ + d U d ζ + U (1 − U ) = 0 , U (0) = 0 , U → ζ → −∞ , (34)d U d ζ + d U d ζ + U (1 − U ) = 0 , U (0) = 0 , U → ζ → −∞ . (35)The solution of these boundary value problems are U ( ζ ) = (1 − e − ζ ) , (36) U ( ζ ) = (cid:18) −
12 + ζ (cid:19) e − ζ + 12 e − ζ , (37) U ( ζ ) = e − ζ (cid:2) − e − ζ (cid:0) e − ζ (cid:1)(cid:3) − ζe − ζ (cid:18) e − ζ + 12 ζ + 12 (cid:19) , (38)and Maple code to generate these solutions are available on GitHub. Combining theinner and outer solution leads to U ( z ) = U ( z ) + c − U ( z ) + c − U ( z ) + O ( c − ), where U ( z ), U ( z ), U ( z ) correspond to Equations (36)–(38), respectively, written in terms of25he original variable z = ζ/c . By truncating this series at different orders we are able tocompare O (1), O ( c − ) and O ( c − ) perturbation solutions.Results in Figure 8 compare the numerical solutions of Equations (4)–(7) with variousperturbation solutions for fast receding travelling waves. Results in Figure 8(a)–(c) showlate–time numerical solutions of the PDE model in blue with c = − . , − .
49 and -2.99,respectively. In each subfigure, the O (1) and O ( c − ) perturbation solutions are plotted, inred and yellow, respectively. For these results we have not plotted the O ( c − ) perturbationsolution in order to keep Figure 8 easy to interpret. As expected we see that the matchbetween the numerical and perturbation solutions improves as c decreases, and we seethat the O ( c − ) perturbation solutions are more accurate than the O (1) perturbationsolutions. Results in Figure 8(d)–(f) show a magnified comparison of the three solutionsand the regions shown are highlighted in the dashed box in Figure 8(a)–(c).26 c)(b)(a) (d)(e)(f) Figure 8:
Perturbation solutions for slow receding travelling waves. (a)–(c) showplots of the shape of the travelling waves for c = − . , − .
49 and -2.99, respectively.Late time numerical solutions of Equations (4)–(7) are shown in blue, and the O (1) and O ( c − ) perturbation solutions are plotted in red and yellow, respectively. (d)–(f) showthe magnified regions highlighted by the dashed boxes in (a)–(c), respectively.For all travelling wave solutions we have κ = − c/V (0), and as c → −∞ we canestimate V (0) using our perturbation solutions to provide additional insight into the27elationship between κ and c by evaluating the following expressions, O (1) : κ = − c d U (0)d z , (39) O (cid:18) c (cid:19) : κ = − c d U (0)d z + 1 c d U (0)d z , (40) O (cid:18) c (cid:19) : κ = − c d U (0)d z + 1 c d U (0)d z + 1 c d U (0)d z , (41)where we must differentiate our expressions for where U ( z ), U ( z ) and U ( z ) with re-spect to z . Substituting our perturbation solutions into Equations (39)–(41) and thenexpanding the resulting terms as c → −∞ gives O (1) : κ ( c ) = − O (cid:18) c (cid:19) , (42) O (cid:18) c (cid:19) : κ ( c ) = − c + O (cid:18) c (cid:19) , (43) O (cid:18) c (cid:19) : κ ( c ) = − c − c + O (cid:18) c (cid:19) , (44)which provides us with a simple way to relate κ and c as c → −∞ . To explore the accuracyof these approximations we use numerical solutions in the phase plane to estimate κ inthe interval − < c < − c and κ in Figure 9. We also superimpose the various approximations, given by Equations(42)–(44) in Figure 9, where we see that κ → − + as c → −∞ , and that Equation (44)gives an excellent approximation of κ for c < − Relationship between c and κ near c → −∞ . The numerical estimateof κ as a function of c is given in green. Various perturbation approximations given byEquation (42)–(44) are given in red, yellow and blue, respectively. Various relationshipsbetween c and κ are shown in an inset, for − < c < − In this work we discuss approaches for modelling biological invasion and retreat. Themost commonly–used model to mimic biological invasion is the Fisher–KPP model [Mur-ray 2002], and generalisations of the Fisher–KPP model, such as the Porous–Fishermodel [Murray 2002, Witelski 1995]. While these single–species single PDE models havebeen used to model biological invasion in various contexts, they cannot be used to sim-ulate biological recession. As an alternative we highlight the use of the Fisher–Stefanmodel [El–Hachem et al. 2019], which is a different generalisation of the Fisher–KPPmodel that is obtained by reformulating the classical model as a moving boundary prob-lem. There are both advantages and disadvantages of reformulating the Fisher–KPP29odel as a moving boundary problem. One advantage of using the Fisher–Stefan modelis that it always predicts a well–defined sharp front and it has the ability to model bothbiological invasion and recession or retreat. These advantages are both very attractivebecause experimental observations of biological invasion typically observe well–definedsharp fronts [Maini et al. 2004, Maini et al. 2004b] and it is well–known that motile andproliferative populations can both invade and recede. The Fisher–KPP model cannotdescribe either of these features. A disadvantage of using the Fisher–Stefan model is theneed to specify the constant κ . While estimates of these kinds of parameters are well–known in the heat and mass transfer literature for modelling physical processes [Crank1987,Hill 1987], there are no such estimates for these parameters in a biological or ecolog-ical context that we are aware of. Part of the motivation for the analysis in this work is toprovide numerical and approximate analytical insight into the relationship between κ and c . We are motivated to do this because measurements of c are often reported [Maini etal. 2004, Maini et al. 2004b, Simpson et al. 2007] and so understanding how to interpretan estimate of c in terms of κ is of interest.In this work we compare the Fisher–KPP model and the Fisher–Stefan model andit is interesting to consider how these models can be used to interpret experimentalobservations. As we have already discussed, experimental estimates of c are the moststraightforward measurement to obtain. For example, Maini et al. [Maini et al. 2004]use a scratch assay to obtain an estimate of ˆ c , whereas Simpson et al. [Simpson et al.2007] report estimates of ˆ c using experimental observations of cell invasion within intactembryonic tissues. With these measurements of ˆ c it is possible to estimate the product ofthe diffusivity and the proliferation rate since ˆ c = 2 (cid:112) ˆ λ ˆ D for the Fisher–KPP model. Ifwe then assume a typical value for a cell proliferation rate by assuming a typical doublingtime is, say, 24 h, we have ˆ λ = ln (2) /
24 /h and this can then be used to estimate ˆ D byassuming that the Fisher–KPP model is relevant. This approach was followed by Mainiet al. [Maini et al. 2004, Maini et al. 2004b] and Simpson et al. [Simpson et al. 2007].Unfortunately this simple approach does not provide any certainty that the Fisher–KPPmodel is valid. Indeed, with more experimental effort it is possible to carefully analyse aproliferation assay to provide an estimate of ˆ λ [Browning et al. 2017], and to either trackindividual cells or to chemically–inhibit proliferation to obtain an independent estimateof ˆ D [Simpson et al. 2013]. If these more careful experiments are performed then it is30hen possible to examine if the relationship ˆ c = 2 (cid:112) ˆ λ ˆ D is indeed true. If it happens thatthis relationship does not hold then it is possible that the Fisher–Stefan model provides abetter explanation of the data since it is always possible to choose a value of ˆ κ to matchindependent estimates of ˆ D , ˆ λ and ˆ c .In conclusion we would like to mention that all of the models discussed in this workmake the very simple but extremely common assumption that the proliferation of individ-uals is given by a logistic source term. This assumption is widely invoked in many singlespecies models of invasion, including the Fisher–KPP model [Maini et al. 2004, Maini etal. 2004b, Simpson et al. 2007], the Porous–Fisher model [Buenzli et al. 2020, Sherrattand Murray 1990, Witelski 1995] and the Fisher–Stefan model [Du and Lin 2010, El–Hachem et al. 2019], as well as many more complicated multiple species analogues ofthese models [Chaplain et al. 2020, Painter and Sherratt 2003, Painter et al. 2015]. Weacknowledge that there is another class of models where different source terms are used,such as the bistable equation and various models that describe Allee effects [Courchampet al. 2008, Fadai and Simpson 2020, Fife 1979, Johnston et al. 2017, Lewis and Kareiva1993, Taylor and Hastings 2005]. These models are similar to the classical Fisher–KPPmodel except that the quadratic source term is generalised to a cubic polynomial, and itis well–known that such single species models can be used to simulate both biological andinvasion and retreat by changing the shape of the cubic source term. In this work we havedeliberately not focused on Allee–type models so that we do not conflate models of Alleeeffects with the Fisher–Stefan model. Of course, it would be very interesting to consideran extension of the Fisher–Stefan model with a more general source term [Browning etal. 2017, Tsoularis and Wallace 2002], such as an Allee effect. We anticipate many of thenumerical, phase plane and perturbation tools developed in this work would also play arole in the analysis of a Fisher–Stefan–type model with a generalised source term. Weleave this extension for future consideration.31 Appendix A: Numerical methods
To obtain numerical solutions of the Fisher–Stefan equation ∂u∂t = ∂ u∂x + u (1 − u ) , (45)for 0 < x < L ( t ) and t >
0, we first use a boundary fixing transformation ξ = x/L ( t ) sothat we have ∂u∂t = 1 L ( t ) ∂ u∂ξ + ξL ( t ) d L ( t )d t ∂u∂ξ + u (1 − u ) , (46)on the fixed domain, 0 < ξ <
1, for t >
0. Here L ( t ) is the length of the domain that wewill discuss later. To close the problem we also transform the boundary conditions giving ∂u∂ξ = 0 at ξ = 0 , (47) u = 0 at ξ = 1 . (48)We spatially discretise Equations (46)–(48) with a uniform finite difference mesh,with spacing ∆ ξ , approximating the spatial derivatives using a central finite differenceapproximation, giving u j +1 i − u ji ∆ t = 1( L j ) (cid:32) u j +1 i − − u j +1 i + u j +1 i +1 ∆ ξ (cid:33) + ξL j (cid:18) L j +1 − L 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 ) , (49)for i = 2 , . . . , m −
1, where m = 1 / ∆ ξ + 1 is the total number of spatial nodes on thefinite difference mesh, and the index j represents the time index so that u ji ≈ u ( ξ, t ),where ξ = ( i −
1) ∆ ξ and t = j ∆ t . 32iscretising Equations (47)–(48) leads to u j +12 − u j +11 = 0 , (50) u j +1 m = 0 . (51)To advance the discrete system from time t to t + ∆ t we solve the system of nonlin-ear algebraic equations, Equations (49)-(51), using Newton-Raphson iteration. Duringeach iteration of the Newton–Raphson algorithm we estimate the position of the movingboundary using the discretised Stefan condition, L j +1 = L j − ∆ tκL j (cid:32) u j +1 m − u j +1 m − ∆ ξ (cid:33) . (52)Within each time step the Newton–Raphson iterations continue until the maximumchange in the dependent variables is less than the tolerance (cid:15) . All results in this work areobtained by setting (cid:15) = 1 × − , ∆ ξ = 1 × − and ∆ t = 1 × − , and we find thatthese values are sufficient to produce grid–independent results. However, we recommendthat care be taken when using the algorithms on GitHub when considering larger valuesof κ , which can require a much denser mesh to give grid–independent results.We use the time–dependent solutions to provide an estimate of the travelling wavespeed c ∗ . The estimated wave speed is computed using the discretised position of themoving boundary such as c ∗ = ( L j +1 − L j ) / ∆ t . To construct the phase planes we solve Equations (11)–(12) numerically using Heun’smethod with a constant step size d z . In most cases we are interested in examining tra-jectories 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 solvingEquations (26)–(27) are on the appropriate stable or unstable manifold and sufficientlyclose to (1 , eig function [Mathworks 2020]to calculate the eigenvalues and eigenvectors for the particular choice of c of interest. The33ow of the dynamical system are plotted on the phase planes using the MATLAB quiver function [Mathworks 2020]. Acknowledgements:
We thank Stuart Johnston and Sean McElwain for helpful sug-gestions and feedback. This work is supported by the Australian Research Council(DP200100177, DP190103757). 34 eferences [Bate and Hilker 2019] Bate AM, Hilker FM (2019) Preytaxis and travelling waves in aneco–epidemiological model. Bulletin of Mathematical Biology. 81:995–1030.[Browning et al. 2017] Browning AP, McCue SW, Simpson MJ (2017) A Bayesian com-putational approach to explore the optimal the duration of a cell proliferation assay.Bulletin of Mathematical Biology. 79:188–1906.[Browning et al. 2019] Browning AP, Haridas P, Simpson MJ (2019) A Bayesian sequen-tial learning framework to parameterise continuum models of melanoma invasioninto human skin. Bulletin of Mathematical Biology. 81:676–698.[Buenzli et al. 2020] Buenzli PR, Lanaro M, Wong C, McLaughlin MP, Allenby MC,Woodruff MA, Simpson MJ (2020) Cell proliferation and migration explain porebridging dynamics in 3D printed scaffolds of different pore size. Acta Biomaterialia(available online) DOI.[Byrne 2010] Byrne HM (2010) Dissecting cancer through mathematics: from the cell tothe animal model. Nature Reviews Cancer. 10: 221–230.[Canosa 1973] Canosa J (1973) On a nonlinear diffusion equation describing populationgrowth. IBM Journal of Research and Development. 17:307–313.[Chaplain et al. 2020] Chaplain MAJ, Lorenzi T, Mcfarlane FR (2020) Bridging the gapbetween individual-based and continuum models of growing cell populations. Journalof Mathematical Biology. 80:343–371.[Courchamp et al. 2008] Courchamp F, Berec L, Gascoigne J (2008) Allee effects in ecol-ogy and conservation. Oxford University Press, Oxford.[Crank 1987] Crank J (1987) Free and moving boundary problems. Oxford UniversityPress, Oxford.[Curtin et al. 2020] Curtin L, Hawkins–Daarud A, van der Zee KG, Swanson KR, OwenMR (2020) Speed switch in glioblastoma growth rate due to enhanced hypoxia–induced migration. Bulletin of Mathematical Biology. 82:43.35Dalwadi et al. 2020] Dalwadi MP, Waters SL, Byrne HM, Hewitt IJ (2020). A Mathe-matical framework for developing freezing protocols in the cryopreservation of cells.SIAM Journal on Applied Mathematics. 80: 657–689.[Du and Lin 2010] Du Y, Lin Z (2010) Spreading–vanishing dichotomy in the diffusive lo-gistic model with a free boundary. SIAM Journal on Mathematical Analysis. 42:377–405.[El–Hachem et al. 2019] El–Hachem M, McCue SW, Jin W, Du Y, Simpson MJ (2019)Revisiting the Fisher–Kolmogorov–Petrovsky–Piskunov equation to interpret thespreading–extinction dichotomy. Proceedings of the Royal Society A – Mathematical,Physical and Engineering Sciences. 475:20190378.[El–Hachem et al. 2020] El–Hachem M, McCue SW, Simpson MJ (2020) A sharp–frontmoving boundary model for malignant invasion. Physica D: Nonlinear Phenomena.412:132639.[Fisher 1937] Fisher RA (1937) The wave of advance of advantageous genes. Annals ofEugenics. 7:355–369.[Fadai and Simpson 2020] Fadai NT, Simpson MJ (2020) Population dynamics withthreshold effects give rise to a diverse family of allee effects. Bulletin of MathematicalBiology. 82:74.[Fife 1979] Fife PC (1979) Long time behavior of solutions of bistable nonlinear diffusionequations. Archive for Rational Mechanics and Analysis. 70:31–36.[Flegg et al. 2020] Flegg JA, Menon SN, Byrne HM, McElwain DLS (2020) A currentperspective on wound healing and tumour–induced angiogenesis. Bulletin of Math-ematical Biology. 82:43.[Gaffney and Maini 1999] Gaffney EA, Maini PK (1999) Modelling corneal epithelialwound closure in the presence of physiological electric fields via a moving boundaryformalism. IMA Journal of Mathematics Applied in Medicine and Biology. 16:369–393.[Haridas et al. 2017] Haridas P, McGovern JA, McElwain DLS, Simpson MJ (2017)Quantitative comparison of the spreading and invasion of radial growth phase and36etastatic melanoma cells in a three–dimensional human skin equivalent model.PeerJ. 5:e3754.[Hill 1987] Hill JM (1987) One–dimensional Stefan problems: an introduction. LongmanScientific & Technical, Harlow.[Horgan 2009] Horgan FG (2009) Invasion and retreat: shifting assemblages of dungbeetles amidst changing agricultural landscapes in central Peru. Biodiversity andConservation. 18:3519.[Ibrahim et al. 2000] Ibrahim K, Sourrouille P, Hewitt GM (2000) Are recession popula-tions of the desert locust (Schistocerca gregaria) remnants of past swarms? MolecularEcology. 9:783–791.[Jin et al. 2016] 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.[Jin et al. 2017] Jin W, Shah ET, Penington CJ, McCue SW, Maini PK, Simpson MJ(2017) Logistic proliferation of cells in scratch assays is delayed. Bulletin of Mathe-matical Biology. 79:1028–1050.[Jin et al. 2020] Jin W, Lo K-Y, Sun Y-S, Ting Y-H, Simpson MJ (2020) Quantifying therole of different surface coatings in experimental models of wound healing. ChemicalEngineering Science. 220:115609.[Johnston et al. 2015] Johnston ST, Shah ET, Chopin LK, McElwain DLS, Simpson MJ(2015) Estimating cell diffusivity and cell proliferation rate by interpreting IncuCyteZOOM TM assay data using the Fisher–Kolmogorov model. BMC Systems Biology.9:38.[Johnston et al. 2016] 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.[Johnston et al. 2017] Johnston ST, Baker RE, McElwain DLS, Simpson MJ (2017) Co–operation, competition and crowding: a discrete framework linking allee kinetics,37onlinear diffusion, shocks and sharp–fronted travelling waves. Scientific Reports.7:42134.[Killengreen et al. 2007] Killengreen ST, Ims RA, Yoccoz NG, Br˚athen KA and HendenJ–A, Schott T (2007) Structural characteristics of a low Arctic tundra ecosystemand the retreat of the Arctic fox. Biological Conservation. 135: 459–472.[Kolmogorov et al. 1937] Kolmogorov AN, Petrovskii PG, Piskunov NS (1937) A studyof the diffusion equation with increase in the amount of substance, and its applicationto a biological problem. Moscow University Mathematics Bulletin. 1:1–26.[Lewis and Kareiva 1993] Lewis MA, Kareiva P (1993) Allee dynamics and the spread ofinvading organisms. Theoretical Population Biology. 43:141–158.[Mathworks 2020] MathWorks eig . Retrieved August 2020 fromhttps://au.mathworks.com/help/matlab/ref/eig.html.[Mathworks 2020] MathWorks quiverquiver