A sharp-front moving boundary model for malignant invasion
AA sharp-front moving boundary model for malignantinvasion
Maud El-Hachem , Scott W. McCue , and Matthew J. Simpson ∗ School of Mathematical Sciences, Queensland University of Technology,Brisbane, Queensland 4001, AustraliaMay 11, 2020
Abstract
We analyse a novel mathematical model of malignant invasion which takes theform of a two-phase moving boundary problem describing the invasion of a popula-tion of malignant cells into a population of background tissue, such as skin. Cells inboth populations undergo diffusive migration and logistic proliferation. The inter-face between the two populations moves according to a two-phase Stefan condition.Unlike many reaction-diffusion models of malignant invasion, the moving boundarymodel explicitly describes the motion of the sharp front between the cancer andsurrounding tissues without needing to introduce degenerate nonlinear diffusion.Numerical simulations suggest the model gives rise to very interesting travellingwave solutions that move with speed c , and the model supports both malignantinvasion and malignant retreat, where the travelling wave can move in either thepositive or negative x -directions. Unlike the well-studied Fisher-Kolmogorov andPorous-Fisher models where travelling waves move with a minimum wave speed c ≥ c ∗ >
0, the moving boundary model leads to travelling wave solutions with | c | < c ∗∗ . We interpret these travelling wave solutions in the phase plane and show ∗ To whom correspondence should be addressed. E-mail: [email protected] a r X i v : . [ n li n . PS ] M a y hat they are associated with several features of the classical Fisher-Kolmogorovphase plane that are often disregarded as being nonphysical. We show, numerically,that the phase plane analysis compares well with long time solutions from the fullpartial differential equation model as well as providing accurate perturbation ap-proximations for the shape of the travelling waves. Keywords:
Travelling wave; Reaction-diffusion; Stefan problem; Phase Plane; Cancer;Cell invasion. 2
Introduction
Populations of motile and proliferative cells can give rise to moving fronts that are asso-ciated with cancer progression and malignant invasion [1–4]. Similar invasive phenomenaare associated with wound healing [5–7], development [8, 9] and ecology [10–13]. Mathe-matically, these fronts are often studied using reaction-diffusion equations that are basedupon the well-known Fisher-Kolmogorov model or extensions [14–18]. While such modelsare able to capture certain important features, such as the formation of constant speedtravelling wave solutions, there are other features of the standard Fisher-Kolmogorovmodel that are inconsistent with biological observations. For example, classical travel-ling wave solutions of the Fisher-Kolmogorov model on −∞ < x < ∞ do not involvea well-defined front because the travelling wave solutions do not have compact supportand the cell density is always positive, with u ( x, t ) → x → ∞ . Solutions of theFisher-Kolmogorov model on −∞ < x < ∞ always lead to travelling waves for initialconditions with compact support, and these travelling waves lead to the colonisation ofinitially-vacant regions without ever retreating. These two features are inconsistent withmany experimental observations. Experimental images in Figure 1(a)–(b) show key fea-tures of malignant invasion. Here a population of motile and proliferative melanoma cellsis placed onto the surface of human skin tissues in Figure 1(a). A vertical section throughthe skin tissues show the melanoma invading vertically downward into the surroundingskin cells and we see a clear sharp front between the two subpopulations [19, 20]. Inreality, such fronts can either invade into, or retreat from, the surrounding tissues [21].Neither of these biological features are consistent with travelling wave solutions of theclassical Fisher-Kolmogorov model.One way to extend the Fisher-Kolmogorov model to produce a well-defined front isto introduce nonlinear degenerate diffusion [9, 17, 18, 22–30]. Such models, often calledthe Porous-Fisher model, give rise to travelling wave solutions with a well-defined sharpfront that always lead to advancing travelling waves that never retreat. One potentialweakness of this approach is that the introduction of nonlinear degenerate diffusion leadsto additional model parameters that can be difficult to estimate and interpret [17, 29–31].Another way to introduce a sharp front into the Fisher–Kolmogorov model is to recastthe problem as a moving boundary problem [32–35]. This approach involves studying3 b) (c) (cid:2) (cid:2) (cid:2) (cid:3) D (cid:2) D (cid:3) c e ll den s i t y Cancer Skin(a) position
Figure 1:
Experimental motivation and model schematic . (a) Experimental pro-tocol where a population of motile and proliferative melanoma cells are placed onto thesurface of human skin tissues kept at an air-liquid interface to simulate the in vivo en-vironment. Scale bar is 6 mm. (b) Vertical cross section through the tissues in (a)highlighting the vertical downward invasion of melanoma cells (dark) into surroundingskin tissue (light). The sharp front separating the invading malignant population fromthe surrounding tissues is visually distinct and highlighted in the red rectangle. Imagesin (a)-(b) are reproduced from Haridas [19] with permission. (b) Schematic solution of aone-dimensional partial differential equation solution showing the spatial distribution ofa population of cancer cells and skin cells separated by a sharp front. The cancer cellshave density u ( x, t ), diffusivity D u and proliferation rate λ u . The skin cells have density v ( x, t ), diffusivity D v and proliferation rate λ v .the Fisher-Kolmogorov model on 0 < x < L ( t ), and specifying that u ( L ( t ) , t ) = 0 togive a well-defined front. In this approach a Stefan-condition is applied to determinethe speed of the moving front [32–35]. Such models, sometimes called the Fisher-Stefanmodel [36–38], have been extensively studied using rigorous analysis [39–45] but havereceived far less attention in terms of how the solutions of such free boundary problemsrelate to biological observations. Interestingly, while free boundary problems are routinelyused to study many problems in industrial and applied mathematics [46–48], they areless frequently encountered in the mathematical biology literature.Of course, a key difference between the classical Fisher-Kolmogorov model and thekinds of applications in Figure 1(a)-(b) is that the usual Fisher-Kolmogorov model dealswith just one population of cells, whereas malignant invasion involves one population ofcells invading into another population of cells. To model such applications, the Fisher-Kolmogorov model can be extended to a system of partial differential equations to repre-sent the different cell types present [2, 49–53]. While the Fisher-Kolmogorov and Porous-Fisher models have been extended to deal with multiple interacting populations, theunderlying issues associated with the single population models, described above, also4pply to the multiple population analogue [52].In this work we study a mathematical model of cell invasion that involves describingtwo populations of cells as a moving boundary problem. A schematic of this model inFigure 1(c) shows that we consider two cell populations, such as a population of cancercells invading into a population of skin cells, which is consistent with the experimen-tal images in Figure 1(a)-(b). Cells in both populations undergo linear diffusion andproliferate logistically. The motion of the sharp front is governed by a two-phase Stefancondition [32–35,54–56]. As we will show, various properties of the solutions of this modelare consistent with experimental observations. Namely, this model leads to a well-definedfront and travelling wave solutions that represent either malignant advance or retreat. Itis interesting that the travelling wave analysis of this model is intimately related with theclassical phase plane associated with travelling wave solutions of the Fisher-Kolmogorovmodel. However, for our model we make use of certain trajectories in the classical phaseplane that are normally discarded on the grounds of being nonphysical. Here, in thecontext of a moving boundary problem, these normally-discarded features play key rolesin determining the travelling wave solutions. From this point forward all dimensional variables and parameters are denoted with acircumflex, whereas nondimensional quantities are denoted using regular symbols.
We consider a reaction-diffusion model of a population of cancer cells with density ˆ u (ˆ x, ˆ t ),and a population of skin cells with density ˆ v (ˆ x, ˆ t ). The system of equations can be written5s ∂ ˆ u∂ ˆ t = ˆ D u ∂ ˆ u∂ ˆ x + ˆ λ u ˆ u (cid:18) − ˆ u ˆ K u (cid:19) , − ˆ L u < ˆ x < ˆ s (ˆ t ) , (1) ∂ ˆ v∂ ˆ t = ˆ D v ∂ ˆ v∂ ˆ x + ˆ λ v ˆ v (cid:18) − ˆ v ˆ K v (cid:19) , ˆ s (ˆ t ) < ˆ x < ˆ L v , (2)where the densities are functions of position, ˆ x , and time, ˆ t . Cancer cells undergo diffusivemigration with diffusivity ˆ D u >
0, and proliferate logistically with rate ˆ λ u > K u >
0. Similarly, skin cells undergo diffusive migration withdiffusivity ˆ D v > λ v > K v >
0. The model is defined on the ˆ L u < ˆ x < ˆ L v , with a moving boundaryˆ x = ˆ s (ˆ t ) separating the population of cancer cells, ˆ x < ˆ s (ˆ t ), from the population of skincells, ˆ x > ˆ s (ˆ t ).Since we are interested in cell invasion we focus on travelling wave solutions of Equa-tions (1)-(2) by setting ˆ L u and ˆ L v to be sufficiently large to model an infinite domainproblem. The boundary conditions we consider are ∂ ˆ u∂ ˆ x (cid:12)(cid:12)(cid:12)(cid:12) ˆ x = − ˆ L u = 0 , ∂ ˆ v∂ ˆ x (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ L v = 0 , (3)ˆ u (ˆ s (ˆ t ) , ˆ t ) = 0 , ˆ v (ˆ s (ˆ t ) , ˆ t ) = 0 . (4)This means that we have no flux of cancer cells at the left-most boundary and no flux ofskin cells at the right-most boundary, and the density of both populations is zero at themoving boundary, as in Figure 1(a).We describe the motion of the moving boundary by a two-phase Stefan condition,dˆ s (ˆ t )dˆ t = − ˆ κ u ∂ ˆ u∂ ˆ x (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ s (ˆ t ) − ˆ κ v ∂ ˆ v∂ ˆ x (cid:12)(cid:12)(cid:12)(cid:12) ˆ x =ˆ s (ˆ t ) . (5)Here the speed of the moving boundary is the sum of two terms: the first term on theright of Equation (5) is proportional to the spatial gradient of the cancer cell densityat the moving boundary, ˆ x = ˆ s (ˆ t ), and the second term on the right of Equation (5)is proportional to the spatial gradient of the skin cell density at the moving boundary,ˆ x = ˆ s (ˆ t ). The constants of proportionality, ˆ κ u and ˆ κ u , play an important role in relating6he shape of the density profiles to the speed of the interface. We will consider therelationship between these constants and the speed of the interface later.In this work we consider initial conditions given byˆ u (ˆ x,
0) = ˆ φ (ˆ x ) on − ˆ L u < ˆ x < ˆ s (ˆ t ) , (6)ˆ v (ˆ x,
0) = ˆ ψ (ˆ x ) on ˆ s (ˆ t ) < ˆ x < ˆ L v , (7)such that ˆ φ (ˆ s (0)) = ˆ ψ (ˆ s (0)) = 0. We nondimensionalise the dependent variables by writing u = ˆ u/ ˆ K u and v = ˆ v/ ˆ K v , andwe nondimensionalise the independent variables by writing x = ˆ x (cid:113) ˆ λ u / ˆ D u and t = ˆ λ u ˆ t .In this nondimensional framework our model can be written as ∂u∂t = ∂ u∂x + u (1 − u ) , − L u < x < s ( t ) , (8) ∂v∂t = D ∂ v∂x + λv (1 − v ) , s ( t ) < x < L v , (9)where the boundary conditions are given by ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) x = − L u = 0 , ∂v∂x (cid:12)(cid:12)(cid:12)(cid:12) x = L v = 0 , (10) u ( s ( t ) , t ) = 0 , v ( s ( t ) , t ) = 0 , (11)d s ( t )d t = − κ u ∂v∂x (cid:12)(cid:12)(cid:12)(cid:12) x = s ( t ) − κ v ∂u∂x (cid:12)(cid:12)(cid:12)(cid:12) x = s ( t ) . (12)The nondimensional model has four parameters, D = ˆ D v ˆ D u , λ = ˆ λ v ˆ λ u , κ u = ˆ κ u ˆ K u ˆ D u , κ v = ˆ κ v ˆ K v ˆ D u . (13)In this framework, D is a relative diffusivity, and setting D = 1 means that the cancercells and skin cells are equally motile. In contrast, setting D >
D < λ .We consider numerical solutions of Equations (8)-(9) on a domain with L u = 0 and L v = L , where L is chosen to be sufficiently large to facilitate the numerical simulationof travelling wave solutions. We chose piecewise initial conditions given by u ( x,
0) = φ ( x ) = α, < x < β,α (cid:18) − x − βs (0) − β (cid:19) , β < x < s (0) , (14) v ( x,
0) = ψ ( x ) = α (cid:18) x − s (0) L − β − s (0) (cid:19) , s (0) < x < L − β,α, L − β < x < L, (15)where the parameters α > β > α when we are wellaway from the interface, x = s (0). Near the interface we set the density to be a linearfunction of position. Typical initial conditions in Figure 2 show how varying α , β and s (0) affects the shape of the initial condition. (a) (b) (c) Figure 2:
Initial condition . Three initial conditions on 0 < x <
60 are shown for: (a) α = 0 . β = 20 and s (0) = 30; (b) α = 0 . β = 10 and s (0) = 30; and (c) α = 1, β = 10 and s (0) = 20. To solve Equations (8)-(9) we use boundary fixing transformations to recast the movingboundary problem on two fixed domains. These transformations, ξ = x/s ( t ) and η =8 x − s ( t )) / ( L − s ( t )) + 1, allow us to re-write Equations (8)-(9) as, ∂u∂t = 1 s ( t ) ∂ u∂ξ + ξs ( t ) d s ( t )d t ∂u∂ξ + u (1 − u ) , < ξ < , (16) ∂v∂t = D ( L − s ( t )) ∂ v∂η + (cid:18) − ηL − s ( t ) (cid:19) d s ( t )d t ∂v∂η + λv (1 − v ) , < η < , (17)so that we now have u ( ξ, t ) on 0 < ξ < v ( η, t ) on 1 < η <
2. The transformedboundary conditions are ∂u∂ξ (cid:12)(cid:12)(cid:12)(cid:12) ξ =0 = 0 , ∂v∂η (cid:12)(cid:12)(cid:12)(cid:12) η =2 = 0 , (18) u (1 , t ) = 0 , v (1 , t ) = 0 , (19)d s ( t )d t = − κ u s ( t ) ∂u∂ξ (cid:12)(cid:12)(cid:12)(cid:12) ξ =1 − κ v L − s ( t ) ∂v∂η (cid:12)(cid:12)(cid:12)(cid:12) η =1 . (20)Equations (16)-(17) and the associated boundary conditions can now be solved numeri-cally using a standard central difference approximation for the transformed spatial deriva-tives and a backward Euler approximation for the temporal derivatives. These detailsare given in Appendix A. Typically, we find that numerical solutions of Equations (8)-(12) evolve into constantspeed, constant shape travelling waves, such as those shown in Figure 3(a). In thiscase we have D = λ = 1 so that the cancer cells and skin cells are equally motileand proliferative. The travelling wave profiles in Figure 3(a) are generated by choosingparticular values of κ u and κ v that leads to an invading malignant population movingwith positive speed, c = 0 .
2. In contrast, choosing different values of κ u and κ v canlead to a retreating malignant front, as in Figure 3(e), where we have a travelling wavewith c = − .
2. These two numerical travelling wave solutions in Figure 3(a) and (e) areinteresting, especially when we compare the properties of these travelling waves with themore familiar properties of the travelling wave solutions of the Fisher-Kolmogorov modelwhere there are three important differences:9. The moving boundary model (8)-(12) supports travelling wave solutions with well-defined sharp front whereas the Fisher-Kolmogorov model does not;2. Travelling wave solutions of the moving boundary model (8)-(12) can either advanceor retreat, whereas analogous travelling wave solutions of the Fisher-Kolmogorovmodel only ever advance;3. Travelling wave solutions of the nondimensional moving boundary model (8)-(12)move with speed | c | < c ≥ f)(e)(a) (f)(b) (c)(g) (d)(h) Figure 3:
Travelling wave solutions for D = λ = 1. All partial differential equation solutions are obtained with L = 60, β = 1, α = 0 . s (0) = 30. Results in (a) correspond to κ u = 1 . κ v = 0 .
5. Results in (e) correspond to κ u = 0 . κ v = 1 . c = 0 . c = − .
2. Solutions of Equations (8)-(9) in (a) and (e) show u ( x, t ) in solid yellow and v ( x, t ) in solid green, at t = 20 ,
30 and 40. Phase planes in (b) and (f), and (c) and (g) show the trajectoriescorresponding to the U ( z ) and V ( z ) travelling waves, respectively. Relevant trajectories in (b)-(c) and (f)-(g) are shown in dashed linesupon which we superimpose the solid lines from the numerical solution of Equations (8)-(9) transformed into travelling wave coordinates.Results in (d) and (h) show U ( z ) and V ( z ) as a function of z where results from the phase plane are given in dashed lines superimposedupon the solutions of Equations (8)-(9) shifted so that the moving boundary is at z = 0. In the phase planes the equilibrium points areshown with a black disc and the intersection of the phase plane trajectory with the vertical axis is shown with a pink disc. .5 Phase plane analysis To study travelling wave solutions of the moving boundary model we re-write the gov-erning equations in the travelling wave coordinate, z = x − ct , and seek solutions of theform U ( z ) = u ( x − ct ) and V ( z ) = v ( x − ct ). Writing Equations (8)-(9) in the travellingwave coordinates leads tod U d z + c d U d z + U (1 − U ) = 0 , −∞ < z < , (21) D d V d z + c d V d z + λV (1 − V ) = 0 , < z < ∞ , (22)where the relevant boundary conditions are U ( −∞ ) = 1 , U (0) = 0 , (23) V (0) = 0 , V ( ∞ ) = 1 , (24) c = − κ v d V (0)d z − κ u d U (0)d z . (25)The travelling wave solution for the cancer population, U ( z ), is described by Equations(21) and (23), while the travelling wave solution for the skin population, V ( z ), is describedby Equations (22) and (24). These two travelling waves are coupled through Equation(25), which is associated with the Stefan condition at the moving interface. This meansthat the travelling wave solutions for U ( z ) and V ( z ) can be studied in two separate phaseplanes, and these two phase planes are coupled by Equation (25).To simplify our study of these two phase planes we note that Equation (21) for U ( z )is identical to Equation (22) for V ( z ) when D = λ = 1. Therefore, it is sufficient for usto study Equation (22) for V ( z ) and to recall that setting D = λ = 1 means that ouranalysis of this phase plane corresponds to U ( z ). To make progress we re-write Equation(22) as a first-order system d V d z = X, (26)d X d z = − cD X − λD V (1 − V ) . (27)At this point we remark that Equation (26)-(27) defines a two-dimensional phase plane12or ( V ( z ) , X ( z )) that is identical to the phase plane associated with the well-studiedtravelling wave solutions of the Fisher-Kolmogorov model [16, 18]. Therefore, all thewell-known properties of that phase plane will play a role here in our study of Equations(8)-(9). In particular, the equilibrium points are (0 ,
0) and (1 , ,
0) is a saddle for all c whereas (0 ,
0) is a stable node if c ≥ √ λD or a stablespiral for c < √ λD . Normally, when considering travelling wave solutions of the Fisher-Kolmogorov model, we are interested in the heteroclinic trajectory between (1 ,
0) and(0 , ,
0) when c < √ λD is ruled out on the basis of requiring V ( z ) >
0. This classical argumentgives rise to the well-known condition that c ≥ √ λD for travelling wave solutions of theFisher-Kolmogorov model [16, 18]. In contrast, for our moving boundary model we havea very different situation where, for example, the travelling wave in Figure 3(a) leads to c = 0 . < √ λD .To explore these solutions we show the phase plane corresponding to the travellingwave in Figure 3(a) in Figure 3(b)-(c) for U ( z ) and V ( z ) population, respectively. Inall phase planes, we generate the trajectories numerically using techniques described inAppendix A. Figure 3(b) shows the ( U ( z ) , W ( z )) phase plane, where W ( z ) = d U ( z ) / d z and c = 0 . ,
0) and (0 ,
0) are highlighted, and the heteroclinic trajectory that leaves (1 ,
0) andspirals into (0 ,
0) is shown with a dotted line. Normally, when considering travelling wavesolutions of the Fisher-Kolmogorov model, this heteroclinic trajectory would be regardedas nonphysical since it implies that U ( z ) < z along that trajectory.However, instead of rejecting this trajectory, the travelling wave solution for U ( z ) in ourmoving boundary model simply corresponds to the portion of that heteroclinic trajectoryin the fourth quadrant where U ( z ) ≥
0. The point where the trajectory intersects the U ( z ) = 0 axis corresponds to the slope of the travelling wave at the moving boundary,(0 , W ∗ ( z )). This point of intersection is important because it plays a role in satisfyingEquation (25). To provide an additional check on our phase plane in Figure 3(b) we takethe u ( x, t ) travelling wave profile in Figure 3(a) and superimpose the ( U ( z ) , W ( z )) profilecalculated from that travelling wave as a solid line in the phase plane. This exercise showsthat this solid curve is visually indistinguishable from the first part of the heteroclinictrajectory where U ( z ) ≥
0. 13igure 3(c) shows the ( V ( z ) , X ( z )) phase plane associated with the v ( x, t ) travellingwave profile in Figure 3(a). Again, we highlight the equilibrium points at (1 ,
0) and(0 ,
0) and we show the trajectory moving towards the saddle point at (1 ,
0) along thestable manifold. In the usual study of the Fisher-Kolmogorov model this trajectory is notnormally considered because it is not associated with a heteroclinic trajectory, and indeedthe phase plane in Figure 3(c) indicates that this trajectory originates far away from therelevant region of the phase plane. However, we find that part of the trajectory in the firstquadrant, just before (1 ,
0) where V ( z ) ≥
0, corresponds to the travelling wave solutionfor the v ( x, t ) population. The point at which this trajectory intersects the V ( z ) = 0axis, (0 , X ∗ ( z )), corresponds to the slope of the travelling wave at the moving boundary.Taking the two phase planes in Figure 3(b)-(c) together, the two intersection points W ∗ ( z ) and X ∗ ( z ) are such that they satisfy Equation (25), c = − κ v X ∗ ( z ) − κ u W ∗ ( z ).Therefore, these two intersection points play a critical role in relating the speed of thetravelling wave solution with the constants κ u and κ v .To summarise the results in Figure 3(a)-(c), and to make an explicit connection be-tween the physical solutions of the partial differential equation model and the nonphysicalfeatures of the phase plane trajectories, we superimpose various solutions in Figure 3(d).The solid green and solid yellow lines in Figure 3(d) show long time solutions of Equa-tions (8)-(12) that are shifted so that the moving boundary is at z = 0. The dashedlines in Figure 3(d) shows the U ( z ) and V ( z ) associated with the relevant phase planetrajectories from Figure 3(b)-(c), respectively. In the case of the U ( z ) trajectory we seethat the shape of the trajectory matches the solution from Equations (8)-(9) where z ≤ U ( z ) ≥
0. The phase plane trajectory of U ( z ) for z > U ( z )oscillates about U ( z ) = 0 and this does not correspond to any part of the solution ofEquations (8)-(9). In the case of the V ( z ) profile we see that the shape of the phaseplane trajectory matches the solution from Equations (8)-(9) where z ≥ V ( z ) ≥ V ( z ) for z < V ( z ) < κ u and κ v that lead to c =0 .
2. Results in Figure 3(e)-(h) correspond to different choices of κ u and κ v such thatthe travelling wave leads to a receding front with c = − .
2. Numerical solutions of14quations (8)-(9) in Figure 3(e) show the travelling wave solutions and the phase planesin Figure 3(f)-(g) show the phase plane trajectories associated with the U ( z ) and V ( z )travelling waves. Again, a summary comparing the physical travelling wave solutionsfrom Equations (8)-(9) with the phase plane trajectories is given in Figure 3(h). Thiscomparison shows that the travelling wave solutions of Equations (8)-(9) compare verywell with the physical portion of the phase plane trajectories in Figure 3(f)-(g).The first set of travelling wave solutions we report in Figure 3 correspond to thesimplest possible case where D = λ = 1 so that the skin and cancer cells are equally motileand equally proliferative. Additional results are presented in Figures 4-5 for D (cid:54) = 1 and λ = 1, and for D = 1 and λ (cid:54) = 1, respectively. Results in Figure 4-5 are presented in theexact same format as in Figure 3 where we consider results for c > c < f) (h)(a) (b) (c)(g) (d)(h)(e) (f) Figure 4:
Travelling wave solutions for D (cid:54) = 1 and λ = 1. All partial differential equation solutions are obtained with L = 60, β = 1, α = 0 . s (0) = 30. Results in (a) correspond to D = 0 . κ u = 1 . κ v = 0 .
5. Results in (e) correspond to D = 2, κ u = 0 . κ v = 1 . c = 0 . c = − .
1. Solutions of Equations (8)-(9) in (a)and (e) show u ( x, t ) in solid yellow and v ( x, t ) in solid green, at t = 20 ,
30 and 40. Phase planes in (b) and (f), and (c) and (g) showthe trajectories corresponding to the U ( z ) and V ( z ) travelling waves, respectively. Relevant trajectories in (b)-(c) and (f)-(g) are shownin dashed lines upon which we superimpose the solid lines from the numerical solution of Equations (8)-(9) transformed into travellingwave coordinates. Results in (d) and (h) show U ( z ) and V ( z ) as a function of z where results from the phase plane are given in dashedlines superimposed upon the solutions of Equations (8)-(9) shifted so that the moving boundary is at z = 0. In the phase planes theequilibrium points are shown with a black disc and the intersection of the phase plane trajectory with the vertical axis is shown with apink disc. f)(a)(e) (b)(f) (c) (d)(g) (h) Figure 5:
Travelling wave solutions for D = 1 and λ (cid:54) = 1. All partial differential equation solutions are obtained with L = 60, β = 1, α = 0 . s (0) = 30. Results in (a) correspond to λ = 2, κ u = 1 . κ v = 0 .
5. Results in (e) correspond to λ = 0 . κ u = 0 . κ v = 1 . c = 0 . c = − .
1. Solutions of Equations (8)-(9) in (a)and (e) show u ( x, t ) in solid yellow and v ( x, t ) in solid green, at t = 20 ,
30 and 40. Phase planes in (b) and (f), and (c) and (g) showthe trajectories corresponding to the U ( z ) and V ( z ) travelling waves, respectively. Relevant trajectories in (b)-(c) and (f)-(g) are shownin dashed lines upon which we superimpose the solid lines from the numerical solution of Equations (8)-(9) transformed into travellingwave coordinates. Results in (d) and (h) show U ( z ) and V ( z ) as a function of z where results from the phase plane are given in dashedlines superimposed upon the solutions of Equations (8)-(9) shifted so that the moving boundary is at z = 0. In the phase planes theequilibrium points are shown with a black disc and the intersection of the phase plane trajectory with the vertical axis is shown with apink disc. .6 Perturbation solution for | c | (cid:28) All results in Figures 3–5 rely on numerical solutions of Equations (26)-(27) to exploretrajectories in the phase plane. We now provide additional insight by constructing ap-proximate perturbation solutions to complement these numerical explorations. First were-write Equations (26)-(27) asd X d V = − cX − λV (1 − V ) DX , (28)for which we seek a perturbation solution about c = 0. Substituting the expansion X ( V ) = X ( V ) + cX ( V ) + O ( c ) leads to ordinary differential equations for X ( V ) and X ( V ) that can be solved exactly. With the initial condition X (1) = 0, the two-termperturbation solution can be written as X ( V ) = ± (cid:115) λD (cid:18) − V + 2 V (cid:19) − c ( V − V ) / + √ D ( V − √ V + O ( c ) . (29)Retaining just the first term on the right of Equation (29) gives us an approximationthat we refer to as an O (1) perturbation solution whereas retaining both terms on theright of Equation (29) gives us an approximation that we refer to as an O ( c ) perturbationsolution. We will now explore both these solutions.Results in Figure 6(a)-(b) show the ( U ( z ) , W ( z )) and ( V ( z ) , X ( z )) phase planes for c = 0 .
05, respectively. The numerically-generated trajectories are compared with boththe O (1) and O ( c ) perturbation solutions. Here we see that the O (1) perturbationsolution is a teardrop-shaped homoclinic trajectory to (1 , O (1) perturbation solution is a reasonably accurate approximation of the numericaltrajectory in the fourth quadrant for U ( W ). Similarly, in Figure 6(b) we see that the O (1)perturbation solution is a reasonably accurate approximation of the numerical trajectoryin the first quadrant for V ( X ). We also superimpose the O ( c ) perturbation solutionin Figures 6(a)-(b) but it is difficult to visually distinguish between the O (1) and O ( c )solutions for c = 0 . U ( z ) , W ( z )) and ( V ( z ) , X ( z )) phase planes for c = 0 .
5, respectively. In both cases we see that the O (1) perturbation solutions do not18 a)(c) (d)(b) Figure 6:
Perturbation solution in the phase plane when c > and D = λ = 1 . Phase planes in (a)-(b), (c)-(d) compare numerical phase plane trajectories and pertur-bation solutions for c = 0 .
05 and c = 0 .
5, respectively. Numerical estimates of the U ( W )and V ( X ) trajectories are shown in dashed red and dashed green respectively. The O (1)and O ( c ) perturbation solutions are shown in solid yellow and solid purple, respectively.Equilibrium points are shown with black discs. The points at which the various solutionsintersect the vertical axis are shown with various coloured discs corresponding to thecolour of the particular trajectory. 19rovide an accurate approximation of the numerical trajectories, whereas the O ( c ) pertur-bation solutions compare very well with the physical part of the phase plane trajectoriesin both cases. The comparison of the numerical phase plane trajectories and the pertur-bation solutions in Figure 6 is given for the most fundamental case where D = λ = 1 and c >
0. Additional comparisons for other choices of D , λ and c are provided in AppendixB. The comparison of the numerical phase plane trajectories with the perturbation so-lutions in Figure 6 shows the shape of W ( U ) and X ( V ) in the phase plane. To explorehow these solutions compare in the z plane, we integrate both sides of Equation (29)with respect to z numerically using a forward Euler approximation with constant stepsize, d z = 1 × − . This numerical integration leads to estimates of the shape of thetravelling waves that can be compared with the shapes of the travelling wave obtainedfrom long-time numerical solutions of Equations (8)-(9). Figure 7 compares the shape ofboth V ( z ) and U ( z ) obtained from the O ( c ) perturbation solution with those obtainedfrom Equations (8)-(9), where we see that the shape of both the V ( z ) and U ( z ) profilescompare extremely well for c = ± .
05, as expected. It is also pleasing that the shapeof the profiles compare quite well even for much larger values, c = ± .
5. All results inFigure 7 correspond to the simplest case where D = λ = 1 and additional comparisonsfor other choices of D and λ and are provided in Appendix B. All solutions in Figures 3-5 correspond particular choices of u ( x, v ( x, κ u and κ v thatlead to long time travelling wave solutions. However, we note that numerical simulationsand more rigorous analysis of the simpler single-phase Fisher-Stefan moving boundaryproblems gives rise to a spreading-vanishing dichotomy, whereby certain initial conditionsand parameter values lead to population extinction as t → ∞ [36, 38–45]. The mainfocus of our current work is to study travelling wave solutions since we are interested insituations where both populations are present, such as the images in Figure 1(a)-(b). Tocomplement these solutions in Figure 3–7 we now briefly consider additional numericalsolutions of Equations (8)-(9) where similar extinction behaviour occurs in the two-phase20 b)(a) (d)(c) Figure 7:
Perturbation solution for the shape of the travelling waves for D = λ = 1 . Comparison of U ( z ) and V ( z ) from the O ( c ) perturbation solution (purple solid)with numerical estimates of the travelling wave obtained by solving Equations (8)-(9)and shifting the profiles so that U (0) = V (0) = 0. Numerical estimates of U ( z ) and V ( z )are shown in dashed yellow and dashed green lines, respectively. Results are shown for:(a)–(b) c = ± .
05, and (c)–(d) ± c = − .
5. 21roblem.Figure 8 shows various results when we consider vary the initial condition and/orvalues of κ u and κ v . The first set of results in Figure 8(a)–(c) shows a case in which s (0) = 1. Here we see the solution evolving to travelling wave profile with positive speedof the type we have discussed in some detail. An important point to make here is that s (0) < π/
2, which is a critical length in the corresponding one-phase problem [36, 38].Based on the previously reported studies of the one-phase problem, our interpretation ofthe solution in Figure 8(a)–(c) is that even though s (0) < π/
2, travelling wave solutionsare still possible provided the initial mass (cid:82) s (0)0 u ( x,
0) d x is sufficient to overcome masslost at the moving boundary. On the other hand, in Figure 8(d)–(f) the solution has thesame parameter values as in Figure 8(a)–(c), except that κ u and κ v are now reduced.In this case the moving boundary x = s ( t ) moves to the right and approaches a steadystate value which is less than the critical length π/
2, while the left population u ( x, t )goes extinct as t → ∞ . The extinction is caused by the fact that the rate at which massassociated with the u ( x, t ) population is lost at x = s ( t ) exceeds the rate at which themass of u ( x, t ) is gained by proliferation. These two examples are consistent with thespreading-vanishing dichotomy in the one-phase problem [36, 38].22 a) (b) (c)(e) (f)(d)(g) (h) (i)(k)(j) (l) Figure 8:
Additional solutions with qualitatively different long time behaviour.
Four additional numerical solutions of Equations (8)–(9). Each solution corresponds to D = 1, λ = 1, with L = 30, s (0) = 1, β = 0 . α = 0 .
5. Results in each row correspond todifferent values of κ u and κ v : (a)–(c) corresponds to κ u = 2 . κ v = 0 . κ u = 0 . κ v = 0; (g)–(i) corresponds to κ u = 0 . κ v = 0 . κ u = 0 . κ v = 0 . t = 0 , , ,
12 and 16. The density profiles in (i)-(j) are shownat t = 0 , , , < x <
60; profiles in the middle column show the details of these solutions on0 < x <
10, and profiles in the right column show the evolution of s ( t ).23dditional results in Figures 8(g)–(i) and (j)–(l) show two further solutions withdifferent choices of κ u and κ v . In both these cases we see that the u ( x, t ) profile eventuallybecomes extinct whereas the v ( x, t ) profile eventually forms a travelling wave solutionwith c = − .
1. Subtle differences, highlighted in Figure 8(i) and (l), show the temporalbehaviour in terms of the movement of the interface, s ( t ). The case in Figure 8(l) leadsto a monotonically decreasing s ( t ), whereas the case in Figure 8(i) leads to s ( t ) that isinitially increasing before eventually decreasing at later time. This kind of nonmonotonebehaviour of s ( t ) is very interesting because the standard single phase Fisher-Stefan modelappears to only lead to monotone s ( t ), whereas our two-phase analogue leads to moreinteresting and nuanced behaviours. In this work we consider a novel mathematical model of cell invasion which takes theform of a two-phase moving boundary problem. This modelling strategy is both biologi-cally relevant and mathematically novel. The moving boundary model leads to travellingwave solutions with a clearly defined moving front. This is advantageous over the classi-cal Fisher–Kolmogorov model and extensions because travelling wave solutions of thosemodels do not have this property. From a biological point of view, our model describes themigration and proliferation of two populations of cells, u ( x, t ) and v ( x, t ), and this allowsus to model a population of cancer cells, u ( x, t ), invading into a population of surroundingcells, v ( x, t ). This scenario is relevant to melanoma cells invading into surrounding skincells, as shown in Figure 1(a)-(b). Interestingly, the moving boundary model leads totravelling wave solutions that move in either the positive or negative direction, meaningthat we can simulate malignant invasion as well as malignant retreat. This is very dif-ferent to travelling wave solutions of the Fisher-Kolmogorov and Porous-Fisher modelsbecause those models only ever predict malignant advance and never predict malignantretreat.The two-phase moving boundary model is also very interesting mathematically. Inthis work we analyse travelling wave solutions where we show that the U ( z ) = u ( x − ct )and V ( z ) = v ( x − ct ) travelling waves can be analysed in two separate phase planes. These24wo phase planes are identical to the phase plane that arises in the classical analysis oftravelling wave solutions of the Fisher-Kolmogorov model. This phase plane contains twoequilibria: (i) (1 ,
0) is a saddle for all c ; and (ii) (0 ,
0) is a stable node if c ≥ √ λD ora stable spiral for c < √ λD . Normally, in the case of travelling waves solutions of theFisher–Kolmogorov model we are interested in a heteroclinic trajectory between thesetwo equilibria and so we require c ≥ √ λD to avoid the nonphysical negative densitiesthat arise from spirals in the phase plane. In contrast, travelling wave solutions of ourmoving boundary model have c < √ λD and so these normally-discarded trajectoriesturn out to be very useful.For our two-phase moving boundary model we use numerical simulations and per-turbation methods to confirm that the travelling wave solutions for U ( z ) and V ( z ) areassociated with trajectories in the classical Fisher–Kolmogorov phase plane that are nor-mally disregarded as being nonphysical. In the cases we consider with c >
0, the U ( z )travelling wave is associated with the heteroclinic trajectory that leaves (1 ,
0) along theunstable manifold and eventually spiralling into (0 , U ( z ) >
0. Similarly, the V ( z ) travelling wave is associated with the trajectory thatapproaches (1 ,
0) along the stable manifold. Here we have the restriction that the trav-elling wave is associated with that part of the trajectory near (1 ,
0) where V ( z ) > c < U ( z ) travellingwave is associated with the trajectory that eventually moves into (1 ,
0) from infinity,whereas the V ( z ) travelling wave is associated with the trajectory that eventually spiralsinto (0 , t → ∞ with care todetermine how quickly travelling wave solutions develop. Acknowledgements . This work is supported by the Australian Research Council (DP200100177,DP190103757). We thank Dr Wang Jin for advice on the numerical algorithms, and wethank Emeritus Professor Sean McElwain for his feedback.
Contributions . All authors conceived and designed the study; M.El-H. performed numer-ical calculations. All authors drafted the article and gave final approval for publication.
Competing Interests . We have no competing interests.26
Appendix A: Numerical Methods
Liberally commented MATLAB implementations of all numerical algorithms used to gen-erate the solutions of the differential equations in this work are available on GitHub.
As we explained in the main document, the partial differential equation models are trans-formed to a fixed domain, Equations (16)-(17) on 0 ≤ ξ ≤ ≤ η ≤
2, respectively.To solve these transformed partial differential equations we discretize the ξ and η domainsuniformly. In principle we use m equally-spaced mesh points for ξ , m = 1 / ∆ ξ + 1, and n equally-spaced mesh points for η , n = 1 / ∆ η + 1. In practice we usually implementthe numerical solution with m = n by setting ∆ ξ = ∆ η . This is convenient, but notnecessary.Using a central difference approximation for the transformed spatial variable and animplicit Euler approximation for the temporal derivatives [60, 61], at the central nodeson both meshes we have u j +1 i − u ji ∆ t = (cid:0) u j +1 i +1 − u j +1 i + u j +1 i − (cid:1) ( s j +1 ∆ ξ ) + ξ i ( s j +1 − s j ) (cid:0) u j +1 i +1 − u j +1 i − (cid:1) s j +1 ∆ t ∆ ξ + u j +1 i (cid:0) − u j +1 i (cid:1) , i = 2 , . . . , m − , (30) v j +1 i − v ji ∆ t = D (cid:0) v j +1 i +1 − v j +1 i + v j +1 i − (cid:1) (( L − s j +1 ) ∆ η ) + (2 − η i ) (cid:0) v j +1 i +1 − v j +1 i − (cid:1) ( s j +1 − s j )2∆ t ∆ η ( L − s j +1 )+ λv j +1 i (cid:0) − v j +1 i (cid:1) , i = 2 , . . . , n − , (31)where the subscript i denotes the mesh point and the superscript j denotes the time,where t = j ∆ t .To enforce the boundary conditions we set ∂u/∂ξ = 0 at ξ = 0 and ∂v/∂η = 0 at η = 2, further we set u = v = 0 at the moving boundary where ξ = η = 1, leading to u j +12 = u j +11 , v j +1 n = v j +1 n − , u j +1 m = v j +11 = 0 . (32)27o advance the discrete system from time t to t + ∆ t we solve the system of nonlinearalgebraic equations, Equations (30)-(32), using Newton-Raphson iteration [62]. Duringeach iteration of the Newton-Raphson algorithm we estimate the position of the movingboundary using the discretised Stefan condition, s j +1 − s j ∆ t = − κ u u j +1 m − u j +1 m − s j +1 ∆ ξ − κ v v j +12 − v j +11 ( L − s j +1 ) ∆ η . (33)Within each time step the Newton-Raphson iterations continue until the maximum chancein the dependent variables is less than the tolerance (cid:15) . All results in this work are obtainedby setting (cid:15) = 1 × − , ∆ ξ = ∆ η = 2 . × − and ∆ t = 1 × − , and we find that thesevalues are sufficient to produce grid-independent results. However, we recommend thatcare be taken when using the algorithms on GitHub for different choices of parameters,especially when considering larger values of κ u and κ v , which can require a much densermesh to give grid-independent results. To construct the phase planes we solve Equations (26)–(27) 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 [63] to calculatethe eigenvalues and eigenvectors for the particular choice of c , D and λ of interest. Theflow of the dynamical system are plotted on the phase planes using the MATLAB quiver function [64]. 28 Appendix B: Additional results
Results in Figure 6 are presented for D = λ = 1 and c > D = λ = 1 only. Here, in Section 5.1 we present additionalphase plane results where D (cid:54) = 1, λ (cid:54) = 1 and c <
0. Similarly, here in Section 5.2 wepresent additional results where we plot U ( z ) and V ( z ) where D (cid:54) = 1 and λ (cid:54) = 1. In allcases we have a good match between the perturbation solutions and numerical solutionsprovided that the wavespeed is sufficiently close to zero, as expected.29 .1 Additional perturbation results in the phase plane (a)(c) (d)(b) Figure 9:
Perturbation solution for the phase plane trajectories when c < and D = λ = 1 . Phase planes in (a)-(b), (c)-(d) compare numerical phase plane trajec-tories and perturbation solutions for c = − .
05 and c = − .
5, respectively. Numericalestimates of the U ( W ) and V ( X ) trajectories are shown in dashed red and dashed greenrespectively. The O (1) and O ( c ) perturbation solutions are shown in solid yellow andsolid purple, respectively. Equilibrium points are shown with black discs. The points atwhich the various solutions intersect the vertical axis are shown with various coloureddiscs corresponding to the colour of the particular trajectory.30 a) (b)(c) (d) Figure 10:
Perturbation solution for the phase plane trajectories when c > , D = 0 . and λ = 1 . Phase planes in (a)-(b), (c)-(d) compare numerical phase planetrajectories and perturbation solutions for c = 0 .
05 and c = 0 .
5, respectively. Numericalestimates of the U ( W ) and V ( X ) trajectories are shown in dashed red and dashed greenrespectively. The O (1) and O ( c ) perturbation solutions are shown in solid yellow andsolid purple, respectively. Equilibrium points are shown with black discs. The points atwhich the various solutions intersect the vertical axis are shown with various coloureddiscs corresponding to the colour of the particular trajectory.31 a) (b)(c) (d) Figure 11:
Perturbation solution for the phase plane trajectories when c < , D = 0 . and λ = 1 . Phase planes in (a)-(b), (c)-(d) compare numerical phaseplane trajectories and perturbation solutions for c = − .
05 and c = − .
5, respectively.Numerical estimates of the U ( W ) and V ( X ) trajectories are shown in dashed red anddashed green respectively. The O (1) and O ( c ) perturbation solutions are shown in solidyellow and solid purple, respectively. Equilibrium points are shown with black discs. Thepoints at which the various solutions intersect the vertical axis are shown with variouscoloured discs corresponding to the colour of the particular trajectory.32 a) (b)(c) (d) Figure 12:
Perturbation solution for the phase plane trajectories when c > , D = 2 and λ = 1 . Phase planes in (a)-(b), (c)-(d) compare numerical phase planetrajectories and perturbation solutions for c = 0 .
05 and c = 0 .
5, respectively. Numericalestimates of the U ( W ) and V ( X ) trajectories are shown in dashed red and dashed greenrespectively. The O (1) and O ( c ) perturbation solutions are shown in solid yellow andsolid purple, respectively. Equilibrium points are shown with black discs. The points atwhich the various solutions intersect the vertical axis are shown with various coloureddiscs corresponding to the colour of the particular trajectory.33 a) (b)(c) (d) Figure 13:
Perturbation solution for the phase plane trajectories when c < , D = 2 and λ = 1 . Phase planes in (a)-(b), (c)-(d) compare numerical phase plane tra-jectories and perturbation solutions for c = − .
05 and c = − .
5, respectively. Numericalestimates of the U ( W ) and V ( X ) trajectories are shown in dashed red and dashed greenrespectively. The O (1) and O ( c ) perturbation solutions are shown in solid yellow andsolid purple, respectively. Equilibrium points are shown with black discs. The points atwhich the various solutions intersect the vertical axis are shown with various coloureddiscs corresponding to the colour of the particular trajectory.34 a) (b)(c) (d) Figure 14:
Perturbation solution for the phase plane trajectories when c > , D = 1 and λ = 0 . . Phase planes in (a)-(b), (c)-(d) compare numerical phase planetrajectories and perturbation solutions for c = 0 .
05 and c = 0 .
5, respectively. Numericalestimates of the U ( W ) and V ( X ) trajectories are shown in dashed red and dashed greenrespectively. The O (1) and O ( c ) perturbation solutions are shown in solid yellow andsolid purple, respectively. Equilibrium points are shown with black discs. The points atwhich the various solutions intersect the vertical axis are shown with various coloureddiscs corresponding to the colour of the particular trajectory.35 a) (b)(c) (d) Figure 15:
Perturbation solution for the phase plane trajectories when c < , D = 1 and λ = 0 . . Phase planes in (a)-(b), (c)-(d) compare numerical phaseplane trajectories and perturbation solutions for c = − .
05 and c = − .
5, respectively.Numerical estimates of the U ( W ) and V ( X ) trajectories are shown in dashed red anddashed green respectively. The O (1) and O ( c ) perturbation solutions are shown in solidyellow and solid purple, respectively. Equilibrium points are shown with black discs. Thepoints at which the various solutions intersect the vertical axis are shown with variouscoloured discs corresponding to the colour of the particular trajectory.36 a)(c) (b)(d) Figure 16:
Perturbation solution for the phase plane trajectories when c > , D = 1 and λ = 2 . Phase planes in (a)-(b), (c)-(d) compare numerical phase planetrajectories and perturbation solutions for c = 0 .
05 and c = 0 .
5, respectively. Numericalestimates of the U ( W ) and V ( X ) trajectories are shown in dashed red and dashed greenrespectively. The O (1) and O ( c ) perturbation solutions are shown in solid yellow andsolid purple, respectively. Equilibrium points are shown with black discs. The points atwhich the various solutions intersect the vertical axis are shown with various coloureddiscs corresponding to the colour of the particular trajectory.37 a)(c) (b)(d) Figure 17:
Perturbation solution for the phase plane trajectories when c < , D = 1 and λ = 2 . Phase planes in (a)-(b), (c)-(d) compare numerical phase plane tra-jectories and perturbation solutions for c = − .
05 and c = − .
5, respectively. Numericalestimates of the U ( W ) and V ( X ) trajectories are shown in dashed red and dashed greenrespectively. The O (1) and O ( c ) perturbation solutions are shown in solid yellow andsolid purple, respectively. Equilibrium points are shown with black discs. The points atwhich the various solutions intersect the vertical axis are shown with various coloureddiscs corresponding to the colour of the particular trajectory.38 .2 Additional perturbation results presented in the z coordi-nate (b)(a) (d)(c) Figure 18:
Perturbation solution for the shape of the travelling waves when D = 0 . and λ = 1 . Comparison of U ( z ) and V ( z ) from the O ( c ) perturbation solution(purple solid) with numerical estimates obtained by solving Equations (16)-(17) that areshifted so that U (0) = V (0) = 0. Numerical estimates of U ( z ) and V ( z ) are shown indashed yellow and dashed green lines, respectively. Results are shown for: (a) c = 0 . c = − .
05; (c) c = 0 .
5; and (d) c = − . b)(a) (d)(c) Figure 19:
Perturbation solution for the shape of the travelling waves when D = 2 and λ = 1 . Comparison of U ( z ) and V ( z ) from the O ( c ) perturbation solution(purple solid) with numerical estimates obtained by solving Equations (16)-(17) that areshifted so that U (0) = V (0) = 0. Numerical estimates of U ( z ) and V ( z ) are shown indashed yellow and dashed green lines, respectively. Results are shown for: (a) c = 0 . c = − .
05; (c) c = 0 .
5; and (d) c = − . b)(a) (d)(c) Figure 20:
Perturbation solution for the shape of the travelling waves when D = 1 and λ = 0 . . Comparison of U ( z ) and V ( z ) from the O ( c ) perturbation solution(purple solid) with numerical estimates obtained by solving Equations (16)-(17) that areshifted so that U (0) = V (0) = 0. Numerical estimates of U ( z ) and V ( z ) are shown indashed yellow and dashed green lines, respectively. Results are shown for: (a) c = 0 . c = − .
05; (c) c = 0 .
5; and (d) c = − . b)(a) (d)(c) Figure 21:
Perturbation solution for the shape of the travelling waves when D = 1 and λ = 2 . Comparison of U ( z ) and V ( z ) from the O ( c ) perturbation solution(purple solid) with numerical estimates obtained by solving Equations (16)-(17) that areshifted so that U (0) = V (0) = 0. Numerical estimates of U ( z ) and V ( z ) are shown indashed yellow and dashed green lines, respectively. Results are shown for: (a) c = 0 . c = − .
05; (c) c = 0 .
5; and (d) c = − . eferences [1] K.R. Swanson, C. Bridge, J.D. Murray, E.C. Alvord Jr. Virtual and real brain tu-mors: using mathematical modeling to quantify glioma growth and invasion. Journalof Neurological Sciences. 216 (2003) 1–10.[2] R.A. Gatenby, E.T. Gawlinski. A reaction–diffusion model of cancer invasion. CancerResearch. 56 (1996) 5745–5753.[3] T. Roose, S.J. Chapman, P.K. Maini. Mathematical models of avascular tumorgrowth. SIAM Review. 49 (2007) 179–208.[4] H.M. Byrne. Dissecting cancer through mathematics: from the cell to the animalmodel. Nature Reviews Cancer. 10 (2010) 221–230.[5] P.K. Maini, D.L.S. McElwain, D.I. Leavesley. Traveling wave model to interpret awound-healing cell migration assay for human peritoneal mesothelial cells. TissueEngineering. 10 (2004) 475–482.[6] P.K. Maini, D.L.S. McElwain, D. Leavesley. Traveling waves in a wound healingassay. Applied Mathematics Letters. 17 (2004) 575–580.[7] M.J. Simpson, K.K. Treloar, B.J. Binder, P. Haridas, K.J. Manton, D.I. Leavesley,D.L.S. McElwain, R.E. Baker. Quantifying the roles of motility and proliferation ina circular barrier assay. Journal of the Royal Society Interface. 10 (2013) 20130007.[8] M.J. Simpson, D.C. Zhang, M. Mariani, K.A. Landman, D.F. Newgreen. Cell prolif-eration drives neural crest cell invasion of the intestine. Developmental Biology. 302(2007) 553–568.[9] B.G. Sengers, C.P. Please, R.O.C. Oreffo. Experimental characterization and compu-tational modelling of two-dimensional cell spreading for skeletal regeneration. Jour-nal of the Royal Society Interface. 4 (2007) 1107–1117.[10] J.G. Skellam. Random dispersal in theoretical populations. Biometrika. 38 (1951)196–218.[11] N. Shigesada, K. Kawasaki, Y. Takeda. Modeling stratified diffusion in biologicalinvasions. American Naturalist. 146 (1995) 229–251.4312] P. Broadbridge, B.H. Bradshaw, G.R. Fulford, G.K. Aldis. Huxley and Fisher equa-tions for gene propagation: An exact solution. ANZIAM Journal. 44 (2002) 11–20.[13] B.H. Bradshaw-Hajek, P. Broadbridge. A robust cubic reaction-diffusion system forgene propagation. Mathematical and Computer Modelling. 39 (2004) 1151–1163.[14] R.A. Fisher. The wave of advance of advantageous genes. Annals of Eugenics. 7(1937) 355–369.[15] A.N. Kolmogorov, P.G. Petrovskii, N.S. Piskunov. A study of the diffusion equationwith increase in the amount of substance, and its application to a biological problem.Moscow University Mathematics Bulletin. 1 (1937) 1–26.[16] J. Canosa. On a nonlinear diffusion equation describing population growth. IBMJournal of Research and Development. 17 (1973) 307–313.[17] J.A. Sherratt, J.D. Murray. Models of epidermal wound healing. Proceedings of theRoyal Society of London: Series B. 241 (1990) 29–36.[18] J.D. Murray. Mathematical Biology I: An Introduction. Third edition, Springer, NewYork, (2002).[19] P. Haridas J.A. McGovern, S.D.L. McElwain, M.J. Simpson. Quantitative compari-son of the spreading and invasion of radial growth phase and metastatic melanomacells in a three-dimensional human skin equivalent model. PeerJ. 5 (2017) e3754.[20] P. Haridas, A.P. Browning, J.A. McGovern, D.L.S. McElwain, M.J. Simpson. Three-dimensional experiments and individual based simulations show that cell prolifera-tion drives melanoma nest formation in human skin tissue. BMC Systems Biology.12 (2018) 34.[21] D. Hanahan, R.A.Weinberg. The hallmarks of cancer. Cell. 100 (2000) 57–70.[22] F. S´anchez Garduno, P.K. Maini. An approximation to a sharp type solution ofa density-dependent reaction–diffusion equation. Applied Mathematics Letters. 7(1994) 47–51.[23] F. S´anchez Garduno, P.K. Maini. Traveling wave phenomena in some degeneratereaction–diffusion equations. Journal of Differential Equations. 117 (1994) 281–319.4424] T.P. Witelski. An asymptotic solution for traveling waves of a nonlinear–diffusionFisher’s equation. Journal of Mathematical Biology. 33 (1994) 1–16.[25] T.P. Witelski. Merging traveling waves for the porous–Fisher’s equation. AppliedMathematics Letters. 8 (1995) 57–62.[26] S. Harris. Fisher equation with density-dependent diffusion: special solutions. Jour-nal of Physics A: Mathematical and Theoretical. 37 (2004) 6267.[27] J.A. Sherratt, B.P. Marchant. Nonsharp travelling wave fronts in the Fisher equationwith degenerate nonlinear diffusion. Applied Mathematics Letters. 9 (1996) 33–38.[28] W. Jin, E.T. Shah, C.J. Penington, S.W. McCue, L.K. Chopin, M.J. Simpson. Re-producibility of scratch assays is affected by the initial degree of confluence: exper-iments, modelling and model selection. Journal of Theoretical Biology. 390 (2016)136–145.[29] D.J. Warne, R.E. Baker, M.J. Simpson. Using experimental data and informationcriteria to guide model selection for reaction–diffusion problems in mathematicalbiology. Bulletin of Mathematical Biology. 81 (2019) 1760–1804.[30] S.W. McCue, W. Jin, T.J. Moroney, K-Y. Lo, S.E. Chou, M.J. Simpson. Hole-closingmodel reveals exponents for nonlinear degenerate diffusivity functions in cell biology.Physica D: Nonlinear Phenomena. 398 (2019) 130–140.[31] M.J. Simpson, R.E. Baker, S.W. McCue. Models of collective cell spreading withvariable cell aspect ratio: A motivation for degenerate diffusion models. PhysicalReview E. 83 (2011) 021901.[32] J. Crank. Free and moving boundary problems. Oxford University Press, Oxford(1987).[33] J.M. Hill. One–dimensional Stefan problems: an introduction. Longman Scientific &Technical, Harlow (1987).[34] S.C. Gupta. The classical Stefan problem. Basic concepts, modelling and analysiswith quasi-analytical solutions and methods. Second edition. Elsevier, Amsterdam(2017). 4535] S.W. McCue, B. Wu, J.M. Hill. Classical two-phase Stefan problem for spheres. Pro-ceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences.464 (2008) 2055–2076.[36] M. El-Hachem, S. W. McCue, W. Jin, Y. Du, M. J. Simpson. Revisiting the Fisher–Kolmogorov–Petrovsky–Piskunov equation to interpret the spreading–extinction di-chotomy. Proceedings of the Royal Society A: Mathematical, Physical and Engineer-ing Sciences. 475 (2019) 20190378.[37] N.T. Fadai, M. J. Simpson. New travelling wave solutions of the Porous–Fisher modelwith a moving boundary. Journal of Physics A: Mathematical and Theoretical. 53(2020), 095601.[38] M.J. Simpson. Critical length for the spreading-vanishing dichotomy in higher di-mensions. ANZIAM Journal. In press, (2020) arXiv preprint[39] Y. Du, Z. Lin. Spreading–vanishing dichotomy in the diffusive logistic model with afree boundary. SIAM Journal on Mathematical Analysis. 42 (2010) 377–405.[40] Y. Du, Z. Guo. Spreading–vanishing dichotomy in a diffusive logistic model with afree boundary, II. Journal of Differential Equations. 250 (2011) 4336–4366.[41] G. Bunting, Y. Du, K. Krakowski. Spreading speed revisited: Analysis of a freeboundary model. Networks & Heterogeneous Media. 7 (2012) 583–603.[42] Y. Du, Z. Guo. The Stefan problem for the Fisher-KPP equation. Journal of Differ-ential Equations. 253 (2012) 996–1035.[43] Y. Du, H. Matano, K. Wang. Regularity and asymptotic behavior of nonlinear Stefanproblems. Archive for Rational Mechanics and Analysis. 212 (2014) 957–1010.[44] Y. Du, H. Matsuzawa, M. Zhou. Sharp estimate of the spreading speed determinedby nonlinear free boundary problems. SIAM Journal on Mathematical Analysis. 46(2014) 375–396.[45] Y. Du, B. Lou. Spreading and vanishing in nonlinear diffusion problems with freeboundaries. Journal of the European Mathematical Society. 17 (2015) 2673–2724.4646] F. Font, S.L. Mitchell, T.G. Myers. One–dimensional solidification of supercooledmelts. International Journal of Heat and Mass Transfer. 62 (2013) 411–421.[47] S.L. Mitchell and S.B.G. O’Brien. Asymptotic and numerical solutions of a freeboundary problem for the sorption of a finite amount of solvent into a glassy polymer.SIAM Journal on Applied Mathematics. 74 (2014) 697–723.[48] M.P. Dalwadi, S.L. Waters, H.M. Byrne, I.J. Hewitt. A Mathematical frameworkfor developing freezing protocols in the cryopreservation of cells. SIAM Journal onApplied Mathematics. 80 (2020) 657–689.[49] K.A. Landman, G.J. Pettet. Modelling the action of proteinase and inhibitor intissue invasion. Mathematical Biosciences. 154 (1998) 23—37.[50] A.J. Perumpanani, J.A. Sherratt, J. Norbury, H.M. Byrne, A two parameter familyof travelling waves with a singular barrier arising from the modelling of extracellularmatrix mediated cellular invasion. Physica D: Nonlinear Phenomena. 126 (1999)145–15.[51] K.J. Painter, J.A. Sherratt. Modelling the movement of interacting cell populations.Journal of Theoretical Biology. 225 (2003) 327–339.[52] M.J. Simpson, K.A. Landman, B.D. Hughes, D.F. Newgreen. Looking inside aninvasion wave of cells using continuum models: Proliferation is the key. Journal ofTheoretical Biology. 243 (2006) 343–360.[53] A.P. Browning, P. Haridas, M.J. Simpson. A Bayesian sequential learning frameworkto parameterise continuum models of melanoma invasion into human skin. Bulletinof Mathematical Biology. 81 (2019) 676–698.[54] S.L. Mitchell, M. Vynnycky. On the numerical solution of two–phase Stefan prob-lems with heat-flux boundary conditions, Journal of Computational and AppliedMathematics. 264 (2014) 49–64.[55] S.L. Mitchell. Applying the combined integral method to two-phase Stefan prob-lems with delayed onset of phase change. Journal of Computational and AppliedMathematics. 281 (2015)58–73. 4756] S.L. Mitchell, M. Vynnycky. On the accurate numerical solution of a two–phaseStefan problem with phase formation and depletion, Journal of Computational andApplied Mathematics. 300 (2016) 259–274.[57] K.K. Treloar, M.J. Simpson, D.L.S. McElwain, R.E. Baker. Are in vitro estimatesof cell diffusivity and cell proliferation rate sensitive to assay geometry? Journal ofTheoretical Biology. 356 (2014) 71–84.[58] K. Smallbone, D.J. Gavaghan, R.A. Gatenby, P.K. Maini. The role of acidity in solidtumour growth and invasion. Journal of Theoretical Biology. 235 (2005) 476–484.[59] A.B. Holder, M.R. Rodrigo, M.A. Herrero. A model for acid-mediated tumour growthwith nonlinear acid production term. Applied Mathematics and Computation. 227(2014) 176–198.[60] M.J. Simpson, K.A. Landman, T.P. Clement. Assessment of a non-traditional opera-tor split algorithm for simulation of reactive transport. Mathematics and Computersin Simulation. 70 (2005) 44–60.[61] M.J. Simpson, K.A. Landman, K. Bhaganagarapu. Coalescence of interacting cellpopulations. Journal of Theoretial Biology. 247 (2007) 525–543.[62] R.L. Burden, J.D. Faires. Numerical analysis. Ninth edition, Brooks/Coles, Boston(2011).[63] MathWorks eig . Retrieved May 2020 https://au.mathworks.com.[64] MathWorks quiverquiver