# Travelling wave analysis of cellular invasion into surrounding tissues.

Maud El-Hachem  Scott W McCue  Matthew J Simpson  School of Mathematical
Sciences, Queensland University of Technology, Brisbane, Australia.

###### Abstract

Single-species reaction-diffusion equations, such as the Fisher-KPP and
Porous-Fisher equations, support travelling wave solutions that are often
interpreted as simple mathematical models of biological invasion. Such
travelling wave solutions are thought to play a role in various applications
including development, wound healing and malignant invasion. One criticism of
these single-species equations is that they do not explicitly describe
interactions between the invading population and the surrounding environment.
In this work we study a reaction-diffusion equation that describes malignant
invasion which has been used to interpret experimental measurements describing
the invasion of malignant melanoma cells into surrounding human skin tissues
[1]. This model explicitly describes how the population of cancer cells
degrade the surrounding tissues, thereby creating free space into which the
cancer cells migrate and proliferate to form an invasion wave of malignant
tissue that is coupled to a retreating wave of skin tissue. We analyse
travelling wave solutions of this model using a combination of numerical
simulation, phase plane analysis and perturbation techniques. Our analysis
shows that the travelling wave solutions involve a range of very interesting
properties that resemble certain well-established features of both the Fisher-
KPP and Porous-Fisher equations, as well as a range of novel properties that
can be thought of as extensions of these well-studied single-species
equations. Of particular interest is that travelling wave solutions of the
invasion model are very well approximated by trajectories in the Fisher-KPP
phase plane that are normally disregarded. This observation establishes a
previously unnoticed link between coupled multi-species reaction diffusion
models of invasion and a different class of models of invasion that involve
moving boundary problems. Matlab software to implement all calculations is
available at GitHub.

###### keywords:

Invasion; Cell invasion; Reaction-diffusion; Partial differential equation;
Cancer; Phase plane.

##  1 Introduction

The Fisher-KPP model [2, 3] is a very simple prototype mathematical model of
biological invasion that describes the spatial and temporal evolution of a
population where individuals undergo migration by linear diffusion and
proliferation via logistic growth. The Porous-Fisher model is an extension of
the Fisher-KPP model where the linear diffusion term is generalised to a
nonlinear degenerate diffusion term with a power law diffusivity [4, 5, 6, 7,
8]. Such single species partial differential equation (PDE) models have had a
major influence on the study of biological populations, both in cell biology
[9, 11, 10, 12, 13, 14, 15, 16] and in ecology [17, 18, 19, 20, 21, 22, 23],
since these models gives rise to travelling wave solutions that are thought to
represent invasion waves [24, 4]. While influential, an obvious limitation of
such single species models is that they focus on the properties of the
invading population alone, and neglect interactions between the invasive
population and the environment, or interactions between the invasive
population and other populations of interest. To overcome this limitation, a
number of extended multi-species models have been presented that explicitly
describe coupling between the invasive population and the environment or other
populations of interest, for example [25, 26, 27].

Figure 1: Experimental motivation and mathematical model schematic. (a) Cross
section through human skin tissues showing the invasion of melanoma cells
(dark brown) into surrounding skin tissue (light brown). The direction of
invasion is shown with the arrow. This image is reproduced from Haridas [28]
with permission. (b) Schematic solution of a one-dimensional PDE model showing
the spatial distribution of cell populations during invasion. To be consistent
with the experiments in (a), the density of cancer cells (brown) moves in the
direction of the black arrow, and this invasion is associated with the retreat
of the density of surrounding skin tissues (blue). The overlap region involves
a visually-distinct region where both populations are present at the same
location. (c) Shows travelling wave solutions for the cancer cell density
$U(z)$, and the skin cell density $V(z)$, on $-\infty<z<\infty$. In this work
we will refer to three different regions of the domain: (i) region I is the
invaded region where $U(z)\rightarrow 1$ and $V(z)\rightarrow 0$ as
$z\to-\infty$; (ii) region II is the overlap region where $U(z)>0$ and
$V(z)>0$; and, (iii) region III is the uninvaded region where $U(z)\rightarrow
0$ and $V(z)\rightarrow\mathcal{V}$ as $z\to\infty$.

One area of research where coupled mathematical models of invasion have played
an important role is in the study of cancer biology where the invasion of a
population of malignant cells is tightly coupled to the degradation of
surrounding tissues. The experimental image in Figure 1(a) shows a population
of highly aggressive melanoma cells growing within and invading into human
skin tissues [28, 29]. During these experiments, melanoma cells simultaneously
migrate and proliferate to form an invading front, and this invasion is
tightly coupled to the biochemical breakdown of the surrounding tissues by
proteases released by the melanoma cells. From a mathematical modelling point
of view, this kind of process leads to conceptual models like that shown in
Figure 1(b) where we think of a density profile of invading cancer cells that
is coupled to the retreat of a population of surrounding skin cells. In this
schematic we identify three regions: region I is the invading region; region
II is the overlap region; and, region III is the uninvaded region. Throughout
this work we refer to the density of invading cells as $\hat{u}$ and the
density of surrounding tissues as $\hat{v}$.

One of the first mathematical models of cellular invasion was proposed by
Gatenby and Gawlinski [30], who present a coupled system of PDEs that
described the spatial and temporal development of tumour tissue, normal tissue
and excess hydrogen ion concentration. Numerical exploration and phase plane
analysis suggested the formation of a pH gradient across the tumour-host
interface, leading to a hypocellular interstitial gap at the tumour-host
interface. This gap was then verified using both in vitro and in vivo data. A
similar model of cellular migration coupled to the degradation of surrounding
tissues in the context of developmental biology was later proposed by Landman
and Pettet [31]. In this case the simpler mathematical model was solved
exactly to reveal details of the coupling between the invading population and
the degradation of surrounding tissues. Since these first two mathematical
models were proposed and analysed in the late 1990s, a range of more detailed
mathematical models have since studied to examine different aspects of
cellular invasion [32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42].

Despite the fact that mathematical models of cellular invasion into
surrounding tissues had been considered for over twenty years, it was not
until 2019 that one of these models was first quantitatively calibrated to
experimental data. In 2019, Browning and co-workers [1] examined a simplified
model based on Gatenby and Gawlinski’s earlier, more general modelling
framework [30]. In this work, Browning took experimental data describing time-
series measurements of melanoma invasion into human tissues (Figure 1(a)) and
used a Bayesian sequential learning approach to estimate the diffusivity of
the melanoma cells, the proliferation rate of the melanoma cells and the rate
at which melanoma cells degraded the surrounding skin tissues [1]. This 2019
study was different to many previous mathematical studies of cellular invasion
since Browning did not consider any travelling wave solutions or travelling
wave analysis. In the present study we re-visit the model proposed by Browning
and explore various travelling wave solutions of that model. Using a
combination of numerical methods to solve the full time-dependent PDE model,
phase plane analysis and perturbation methods, we reveal several novel
features of the travelling wave solutions of this model. In particular we
unearth many important parallels and differences between the travelling wave
solutions of the invasion model and the very well-studied Fisher-KPP model.
There are two aspects of our analysis that are of particular interest. First,
we show that often disregarded phase plane analysis for the Fisher-KPP model
plays a role here in the travelling wave analysis of Browning’s model of
cellular invasion. Second, we show that Browning’s model of invasion leads to
travelling wave solutions that are reminiscent of travelling wave solutions of
a moving boundary type model [43, 44, 45].

##  2 Mathematical model and preliminary simulations

###  2.1 Browning’s model of cellular invasion

In 2019, Browning and co-workers [1] proposed the following simple dimensional
model to describe the invasion of melanoma cells into surrounding skin tissue,

$\displaystyle\dfrac{\partial\hat{u}}{\partial\hat{t}}=\hat{D}\dfrac{\partial}{\partial\hat{x}}\left[\left(1-\dfrac{\hat{v}}{\hat{K}}\right)\dfrac{\partial\hat{u}}{\partial\hat{x}}\right]+\hat{\lambda}\hat{u}\left(1-\dfrac{\hat{u}+\hat{v}}{\hat{K}}\right),$
$\displaystyle 0<\hat{x}<\infty,$ (1)
$\displaystyle\dfrac{\partial\hat{v}}{\partial\hat{t}}=-\hat{\gamma}\hat{u}\hat{v},$
$\displaystyle 0<\hat{x}<\infty,$ (2)

where $\hat{u}(\hat{x},\hat{t})>0$ is the density of melanoma cells, and
$\hat{v}(\hat{x},\hat{t})>0$ is the density of skin cells. Throughout this
work we use a circumflex to indicate dimensional parameters and variables.

Equation (1) governs the evolution of the cell density, and we see that the
melanoma cells move according to a nonlinear diffusion term, with diffusivity
$\hat{D}>0$ $[\mathrm{\mu m}^{2}\ \mathrm{h}^{-1}]$. The nonlinear diffusivity
function decreases linearly with the skin density such that the diffusion of
melanoma cells vanishes if the skin density reaches the carrying capacity
density, $\hat{v}(\hat{x},\hat{t})=\hat{K}>0$. Further, equation (1) specifies
that the melanoma cells grow logistically, with rate $\hat{\lambda}>0$
$[\mathrm{h}^{-1}]$, such that the net proliferation rate is a linearly
decreasing function of the total cell density,
$\hat{u}(\hat{x},\hat{t})+\hat{v}(\hat{x},\hat{t})$. Equation (2) governs the
evolution of the skin density such that melanoma cells degrade skin cells at a
rate governed by $\hat{\gamma}\geq 0$ $[\mathrm{\mu
m}^{2}\mathrm{cells}^{-1}\mathrm{h}^{-1}]$. This two-species PDE model is a
simplification of a three-species PDE extension that is fully described in the
Supplementary Material document reported by Browning et al. [1].

Since we are interested in travelling wave solutions we pose Equations (1)–(2)
on $0<\hat{x}<\infty$, however when we generate numerical solutions we take
the usual approach and examine solutions on a truncated domain,
$0<\hat{x}<\hat{L}$. For all numerical solutions, we specify no-flux boundary
conditions for $\hat{u}(\hat{x},\hat{t})$, so that
$\partial\hat{u}(0,\hat{t})/\partial\hat{x}=\partial\hat{u}(\hat{L},\hat{t})/\partial\hat{x}=0$.
Note that since Equation (2) does not involve any spatial derivatives, we do
not need to specify any boundary conditions for $\hat{v}(\hat{x},\hat{t})$.
The choice of $\hat{L}$ is unimportant provided that it is chosen to be
sufficiently large [46]. Matlab software on GitHub can be used to explore
different choices of $\hat{L}$ for all problems that we consider, and full
details of the numerical algorithms are outlined in the Appendix.

The mathematical model is nondimensionalised by introducing
$u=\hat{u}/\hat{K}$, $v=\hat{v}/\hat{K}$,
$x=\hat{x}\sqrt{\hat{\lambda}/\hat{D}}$, $t=\hat{\lambda}\hat{t}$, and
$\gamma=\hat{K}\hat{\gamma}/\hat{\lambda}$, which gives

$\displaystyle\dfrac{\partial u}{\partial t}=\dfrac{\partial}{\partial
x}\left[(1-v)\dfrac{\partial u}{\partial x}\right]+u(1-u-v),$ $\displaystyle
0<x<\infty$ (3) $\displaystyle\dfrac{\partial v}{\partial t}=-\gamma uv,$
$\displaystyle 0<x<\infty,$ (4)

which means that the nondimensional PDE model requires the specification of
just one model parameter, $\gamma\geq 0$.

For the first part of the study, we specify initial conditions such that the
initial distribution of skin density is a constant, and the initial
distribution of melanoma cells has compact support,

$\displaystyle u(x,0)$ $\displaystyle=\alpha\left[1-H(\beta)\right],$ (5)
$\displaystyle v(x,0)$ $\displaystyle=\mathcal{V},$ (6)

where $H(x)$ is the usual Heaviside function and $\beta>0$ is a constant so
that we have $u(x,0)=\alpha$ for $x<\beta$ and $u(x,0)=0$ for $x>\beta$. The
initial density of skin, $0\leq\mathcal{V}\leq 1$, is set to be a constant
that does not depend upon position. In summary, the non-dimensional invasion
model involves one model parameter, $\gamma\geq 0$, and when we specify the
initial conditions we introduce another parameter, $0\leq\mathcal{V}\leq 1$,
that we will study. In general, we interpret $\gamma$ as the rate at which
cancer cells degrade skin tissues, and $\mathcal{V}$ is the density of skin
tissues in the far field ahead of the invading front. Most of our analysis
will be concerned with how varying $\gamma$ and $\mathcal{V}$ influence the
shape and speed of the resulting travelling wave solutions. Setting the upper
value of $\mathcal{V}=1$ implies that the maximum density of skin cells is the
same as the maximum density of cancer cells. This assumption may be
biologically realistic as skin cells and cancer cells are often similar in
shape and size [28].

The first part of this work focuses on travelling wave solutions of the
invasion model that arise from initial conditions given by Equations (5)–(6)
since this form of initial condition with compact support is biologically-
relevant [28, 29]. Once we have characterised these travelling wave solutions
we will then study travelling wave solutions of the invasion model where
$u(x,0)$ decays to zero by setting

$\displaystyle u(x,0)$ $\displaystyle=\begin{cases}\alpha,&x<\beta,\\\
\alpha\,\textrm{exp}[-a(x-\beta)],&x\geq\beta,\end{cases}$ (7) $\displaystyle
v(x,0)$ $\displaystyle=\mathcal{V},$ (8)

where $a>0$ is a constant that we vary so that we can understand how
travelling wave solutions of the invasion model depend upon the initial
spatial decay of $u$.

One of the main themes of our work is to highlight surprising similarities and
differences between the invasion model (3)–(4) and some well-studied single-
species models. For example, setting $\mathcal{V}=0$ reduces the invasion
model to the well-known Fisher-KPP model [2, 3, 4],

$\dfrac{\partial u}{\partial t}=\dfrac{\partial^{2}u}{\partial x^{2}}+u(1-u).$
(9)

As we pointed out in Section 1, the Fisher-KPP model is a simplified
mathematical model of biological invasion that ignores any interaction between
the invading population and the local environment. In this work, we study
travelling wave solutions of Equations (3)–(4) and we find there are important
similarities and differences with travelling wave solutions of the Fisher-KPP
model, so it is useful to explicitly note the connection between these two
mathematical models at the outset. There are also interesting, but perhaps
less obvious, connections between the invasion model and the Porous-Fisher
model [4, 5, 6, 7, 8],

$\dfrac{\partial u}{\partial t}=\dfrac{\partial}{\partial
x}\left[u\dfrac{\partial u}{\partial x}\right]+u(1-u).$ (10)

We will now begin to explore these relationships.

###  2.2 Time-dependant solutions

We begin our study of the invasion model by computationally exploring various
travelling wave solutions of (3)–(4) to develop an intuitive understanding of
how their shape and speed depend upon the choice of $\gamma$ and
$\mathcal{V}$. Full details of the numerical method used to solve the PDE
model is given in the Appendix, and MATLAB software to study time-dependent
solutions of (3)–(4) is available on GitHub.

In the first instance we focus on travelling wave solutions with initial
condition given by Equations (5)–(6) so that the initial cancer density
profile has compact support. Time-dependent solutions are shown in Figure
2(a)–(d) for $\mathcal{V}=0.5$ for a range of $\gamma$. Additional time-
dependent solutions are shown in Figure 2(e)–(h) for $\mathcal{V}=1$ for a
range of $\gamma$. In all cases we see that the time-dependent solutions of
Equations (3)–(4) evolve to constant-speed travelling wave solutions, and in
each subfigure we give a numerical estimate of the travelling wave speed, $c$.

Figure 2: Time-dependant solutions of Equations (3)–(4) for different $\gamma$
and $\mathcal{V}$. Results in (a)–(d) correspond to $\mathcal{V}=0.5$, while
results in (e)–(h) correspond to $\mathcal{V}=1$. Solutions in (a) and (e),
(b) and (f), (c) and (g), and (d) and (h) correspond to $\gamma=0.1,1,10$ and
$100$, respectively. Density profiles for $u(x,t)$ are shown in brown, and
profiles of $v(x,t)$ are shown in blue. Results in (a)–(d) are shown at
$t=0,100,110,120$ and the invasion front moves in the positive $x$-direction,
whereas results in (e)–(h) are shown at $t=0,40,80,120$. For each solution we
show the initial condition in dashed lines, and the subsequent solutions for
$t>0$ are shown in solid lines. All PDE solutions are obtained using $\Delta
x=\Delta t=1\times 10^{-3}$.

Results in Figure 2(a)–(d) highlight several interesting properties of these
travelling wave solutions, especially when we compare them with the well-known
travelling wave solutions of the Fisher-KPP model (9). Each travelling wave
solution in Figure 2(a)–(d) leads to profiles for $u(x,t)$ that are smooth,
and without compact support since the travelling wave density profile decays
to zero for large $x$. In this regard, the shape of the travelling wave
solutions of Equations (3)–(4) is similar to those for the Fisher-KPP model
[4, 24]. However, in the invasion model we see that the speed of the
travelling wave depends upon the decay rate, $\gamma$, in a rather unexpected
way. First, results in Figure 2(a)–(b) suggest that for sufficiently small
$\gamma<\gamma_{\textrm{c}}$, the speed of the travelling wave is independent
of $\gamma$. Second, results in Figure 2(c)–(d) indicate that for sufficiently
large $\gamma>\gamma_{\textrm{c}}$, the travelling wave speed increases with
$\gamma$. It is of interest to note that all values of $\gamma$ considered in
Figure 2(a)–(d) lead to travelling wave solutions with $c<2$. This is very
different to travelling wave solutions of the dimensionless Fisher-KPP model
that evolve from initial conditions with compact support since these
travelling waves always have $c=2$ [4, 24].

Results in Figure 2(e)–(h) for $\mathcal{V}=1$ show that the invasion model
leads to non-smooth travelling wave solutions that have compact support. From
this point of view, the shape of the travelling wave solutions in Figure
2(e)–(h) is reminiscent of the shape of travelling wave solutions of the
Porous-Fisher model (10). Comparing the speeds of the travelling wave
solutions in Figure 2(e)–(h) indicates that $c$ increases with $\gamma$ but,
unlike the results in Figure 2(a)–(d), we do not see any evidence that the
wave speed is independent of $\gamma$ for any of the values considered. As
with Figure 2(a)–(d), all values of $\gamma$ explored in Figure 2(e)–(h) lead
to travelling wave solutions with $c<2$, which, again, is very different to
travelling wave solutions of the dimensionless Fisher-KPP equation (9).

In summary, preliminary numerical explorations reveal that we obtain smooth
travelling wave solutions of the invasion model for $\mathcal{V}<1$ whereas we
obtain non-smooth travelling waves with compact support when $\mathcal{V}=1$.
While the former are similar in shape to the travelling wave solutions of the
Fisher-KPP model and the latter are similar in shape to the travelling wave
solutions of the Porous-Fisher model, the speed of the travelling waves of the
invasion model are very different to the well-known bounds on travelling wave
solutions of the Fisher-KPP and Porous-Fisher models. We summarise the
relationship between $c$, $\gamma$ and $\mathcal{V}$ in Figure 3 where we
estimate the long time travelling wave speed $c$ from a large number of
simulations where $u(x,0)$ has compact support, given by Equation (5). As
noted in Section 2.1, setting $\mathcal{V}=0$ means that the invasion model
simplifies to the Fisher-KPP model and so we obtain $c=c_{\textrm{min}}=2$
regardless of $\gamma$. For intermediate values of $0<\mathcal{V}<1$ we see
that $c_{\textrm{min}}$ is a decreasing function of $\mathcal{V}$, but
independent of $\gamma$ provided that $\gamma$ is sufficiently small. In
contrast for intermediate values of $0<\mathcal{V}<1$ we see that
$c_{\textrm{min}}$ increases with $\gamma$, and approaches
$c_{\textrm{min}}=2$ for large $\gamma$. Finally, for $\mathcal{V}=1$, there
appears to be no positive value of $\gamma$ where the wave speed is
independent of $\gamma$, meaning that we have with $c_{\textrm{min}}\to 0^{+}$
as $\gamma\to 0^{+}$, and $c_{\textrm{min}}\to 2^{-}$ as $\gamma\to\infty$.
The key features of the travelling wave solutions of the invasion model, and
their relationship to the well-studied travelling wave solutions of the
Fisher-KPP and Porous-Fisher models are summarised in Table 1.

Figure 3: Relationship between the minimum travelling wave speed
$c_{\textrm{min}}$, $\gamma$ and $\mathcal{V}$. Numerical estimates of
$c_{\textrm{min}}$ are obtained by solving Equations (3)–(4) with initial
conditions (5)–(6), with $\alpha=1$ and $\beta=10$. Numerical solutions of the
PDE model are obtained with $\Delta x=\Delta t=1\times 10^{-2}$, for $1\times
10^{-3}\leq\gamma\leq 1\times 10^{6}$ and
$\mathcal{V}=1/8,2/8,3/8,4/8,5/8,6/8,7/8$ and $1$. Table 1: Comparing key
features of travelling wave solutions of the invasion model with travelling
wave solutions of the Fisher-KPP model and the Porous-Fisher model.

Solution features | Fisher-KPP (9) | Porous-Fisher (10) | Invasion model
---|---|---|---
$\mathcal{V}<1$ | $\mathcal{V}=1$
Travelling wave shape | Smooth | Sharp-fronted | Smooth | Sharp-fronted
No compact support | Compact support | No compact support | Compact support
Minimum wave speed | $c_{\textrm{min}}=2$ | $c_{\textrm{min}}=\dfrac{1}{\sqrt{2}}$ |  $c_{\textrm{min}}=2(1-\mathcal{V})$, $\gamma<\gamma_{\textrm{c}}$ | $\displaystyle{\lim_{\gamma\to 0^{+}}c_{\textrm{min}}=0^{+}}$
$\displaystyle{\lim_{\gamma\to\infty}c_{\textrm{min}}=2^{-}}$ | $\displaystyle{\lim_{\gamma\to\infty}c_{\textrm{min}}=2^{-}}$

##  3 Travelling wave analysis

We now attempt to formalise these preliminary numerical results by analysing
travelling wave solutions of Equations (3)–(4) in the phase-space.

###  3.1 Preamble

To study the travelling wave solutions of Equations (3)–(4) we seek solutions
in the form $u(x,t)=U(z)$ and $V(x,t)=V(z)$, where $z$ is the usual travelling
wave variable, $z=x-ct$. Re-writing the governing equations in terms of the
travelling wave coordinate gives

$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}z}\left[(1-V)\dfrac{\mathrm{d}U}{\mathrm{d}z}\right]+c\dfrac{\mathrm{d}U}{\mathrm{d}z}+U(1-U-V)$
$\displaystyle=0,$ $\displaystyle-\infty<z<\infty,$ (11) $\displaystyle
c\dfrac{\mathrm{d}V}{\mathrm{d}z}-\gamma UV$ $\displaystyle=0,$
$\displaystyle-\infty<z<\infty,$ (12)

with boundary conditions $U(z)\rightarrow 1$ and $V(z)\rightarrow 0$ as
$z\rightarrow-\infty$, and $U(z)\rightarrow 0$ and
$V(z)\rightarrow\mathcal{V}$ as $z\rightarrow\infty$. At this point it is
worthwhile to observe that if the solution $U(z)$ is known, we can solve
Equation (12) to give,

$V(z)=\mathcal{V}\exp{\left(\dfrac{\gamma}{c}\int_{z}^{\infty}U(\xi)\
\mathrm{d}\xi\right)}.$ (13)

To study this boundary value problem in phase space we re-write Equations
(11)–(12) as a first order system

$\displaystyle\dfrac{\mathrm{d}U}{\mathrm{d}z}$ $\displaystyle=W,$ (14)
$\displaystyle\dfrac{\mathrm{d}V}{\mathrm{d}z}$ $\displaystyle=\dfrac{\gamma
UV}{c},$ (15) $\displaystyle\dfrac{\mathrm{d}W}{\mathrm{d}z}$
$\displaystyle=\dfrac{1}{(1-V)}\left[\dfrac{\gamma
UVW}{c}-cW-U(1-U-V)\right],$ (16)

with boundary conditions $U\rightarrow 1,V\rightarrow 0$ and $W\rightarrow 0$
as $z\rightarrow-\infty$, and $U\rightarrow 0,V\rightarrow\mathcal{V}$ and
$W\rightarrow 0$ as $z\rightarrow\infty$. There are two equilibrium points of
this dynamical system: (i) $(\bar{U},\bar{V},\bar{W})=(1,0,0)$ corresponding
to $z\to-\infty$; and, (ii) $(\bar{U},\bar{V},\bar{W})=(0,\mathcal{V},0)$ for
$\mathcal{V}<1$, corresponding to $z\to\infty$. Before we proceed it is useful
to remark that the dynamical system is singular when $\mathcal{V}=1$, whereas
there is no such singularity for $\mathcal{V}<1$. Therefore, just like we did
in Section 2, we will treat these two cases separately.

###  3.2 $\mathcal{V}<1$

Setting $\mathcal{V}<1$ leads to smooth travelling wave solutions and a
regular phase space so we proceed by analysing Equations (14)–(16) directly.
Eigenvalues of the linearised dynamical system about
$(\bar{U},\bar{V},\bar{W})=(1,0,0)$ are given by the roots of
$\lambda^{3}-(\gamma-c^{2})\lambda^{2}/c-(1+\gamma)\lambda+\gamma/c=0$, so
that we have $\lambda_{1}=\gamma/c$ and
$\lambda_{2,3}=(-c\pm\sqrt{c^{2}+4})/2)$. These eigenvalues are all real with
$\lambda_{1}>0$, $\lambda_{2}>0$ and $\lambda_{3}<0$, which means that the
equilibrium point $(\bar{U},\bar{V},\bar{W})=(1,0,0)$ is a three-dimensional
saddle for all values of $c$ and $\gamma$ [47]. Note that the analogous phase
plane analysis for travelling wave solutions of the Fisher-KPP model (9)
involves an equilibrium point corresponding to the invaded region that is a
two-dimensional saddle for all $c$ [4].

Eigenvalues of the linearised dynamical system about
$(\bar{U},\bar{V},\bar{W})=(0,\mathcal{V},0)$ are given by the roots of
$\lambda^{3}+c\lambda^{2}/(\mathcal{V}-1)+\lambda=0$, so that we have
$\lambda_{1}=0$ and
$\lambda_{2,3}=(-c\pm\sqrt{c^{2}-4(1-\mathcal{V})^{2}})/[2(1-\mathcal{V})]$.
In this case $\lambda_{2}$ and $\lambda_{3}$ are real negative numbers when
$c\geq 2(1-\mathcal{V})$, whereas they are complex conjugates when
$c<2(1-\mathcal{V})$. This means that the equilibrium point associated with
the uninvaded region is a non-hyperbolic stable node when $c\geq
2(1-\mathcal{V})$ and a non-hyperbolic stable spiral when
$c<2(1-\mathcal{V})$. Since we are interested in travelling waves with
$U(z)>0$, this condition defines a minimum wave speed,
$c_{\textrm{min}}=2(1-\mathcal{V})$. Again, note that the analogous phase
plane analysis of the Fisher-KPP model also involves the equilibrium point
corresponding to the uninvaded region bifurcating from a stable node to a
stable spiral, and this defines an analogous minimum wave speed,
$c_{\textrm{min}}=2$ [4].

The phase space analysis leading to a minimum wave speed condition,
$c_{\textrm{min}}=2(1-\mathcal{V})$ is consistent with our numerical results
in Figure 3 computed with initial conditions with compact support. In
particular, for sufficiently small $\gamma<\gamma_{\textrm{c}}$ we have
travelling wave solutions that move with the minimum wave speed,
$c_{\textrm{min}}=2(1-\mathcal{V})$, and this speed is independent of
$\gamma$. For large $\gamma>\gamma_{\textrm{c}}$ our numerical results lead to
travelling wave solutions with $c>c_{\textrm{min}}=2(1-\mathcal{V})$, which is
again consistent with the phase space analysis. Unfortunately, this analysis
provides no insight into the behaviour of the wave speed for sufficiently
large $\gamma$, nor the critical value, $\gamma_{\textrm{c}}$, where the wave
speed dependence upon $\gamma$ changes.

###  3.3 $\mathcal{V}=1$

Setting $\mathcal{V}=1$ leads to nonsmooth travelling wave solutions and a
singularity in Equations (14)–(16). We proceed by defining a new independent
variable $\zeta$ ,
$(1-V)\mathrm{d}(\cdot)/\mathrm{d}z=\mathrm{d}(\cdot)/\mathrm{d}\zeta$ [4], so
that the desingularised dynamical system is given by

$\displaystyle\dfrac{\mathrm{d}U}{\mathrm{d}\zeta}$ $\displaystyle=W(1-V),$
(17) $\displaystyle\dfrac{\mathrm{d}V}{\mathrm{d}\zeta}$
$\displaystyle=\left(\dfrac{\gamma UV}{c}\right)(1-V),$ (18)
$\displaystyle\dfrac{\mathrm{d}W}{\mathrm{d}\zeta}$
$\displaystyle=\left(\dfrac{\gamma UV-c^{2}}{c}\right)W-U(1-U-V),$ (19)

with boundary conditions $U\rightarrow 1,W\rightarrow 0$ and $V\rightarrow 0$
as $\zeta\rightarrow-\infty$, and $U\rightarrow 0,W\rightarrow 0$ and
$V\rightarrow\ \mathcal{V}$ as $\zeta\rightarrow\infty$. Similar to before,
given $U(\zeta)$, the solution for $V(\zeta)$ can be written as

$V(\zeta)=\dfrac{1}{\displaystyle{A\exp{\left(-\dfrac{\gamma}{c}\int_{\zeta}^{\infty}U(\eta)\
\mathrm{d}\eta\right)}}+1},$ (20)

where $A$ is an integration constant.

Eigenvalues of the linearised dynamical system around
$(\bar{U},\bar{V},\bar{W})=(1,0,0)$ are roots of
$\lambda^{3}-(\gamma-c^{2})\lambda^{2}/c-(1+\gamma)\lambda+\gamma/c=0$,
leading to $\lambda_{1}=\gamma/c$, $\lambda_{2,3}=(-c\pm\sqrt{c^{2}+4})/2)$.
This means that $\lambda_{2}$ and $\lambda_{3}$ are real numbers with opposite
sign, and so $(\bar{U},\bar{V},\bar{W})=(1,0,0)$ is a saddle. The eigenvalues
of the linearised dynamical system around $(\bar{U},\bar{V},\bar{W})=(0,1,0)$
are roots of $c\lambda^{2}+\lambda^{3}=0$, leading to $\lambda_{1}=-c$,
$\lambda_{2}=\lambda_{3}=0$, which means that the uninvaded equilibrium point
is always a degenerate stable node. Unlike the previous case where
$\mathcal{V}<1$, here there is no restriction on a minimum wave speed to
ensure $U(\zeta)>0$. One way of interpreting this result is that
$c_{\textrm{min}}=0$ which is consistent with the numerical results in Figure
3 where we have $c>c_{\textrm{min}}=0$ when $\mathcal{V}=1$.

##  4 Limiting cases

As we pointed out in Section 2.1, we are primarily interested in understanding
how travelling wave solutions of the invasion model (3)–(4) depend upon
choices of $\gamma$ and $\mathcal{V}$. We will start by considering limits of
fast and slow decay, $\gamma\gg 1$ and $\gamma\ll 1$, respectively, and
consider differences between $\mathcal{V}=1$ and $\mathcal{V}<1$ as
appropriate.

###  4.1 Fast decay: $\gamma\gg 1$

Preliminary numerical results in Figure 2 indicate that the width of the
overlap region (region II) decreases with $\gamma$. This trend is evident in
Figure 2(a)–(d) for $\mathcal{V}<1$ as well as in Figure 2(e)–(h) where
$\mathcal{V}=1$. One way to interpret this observation is that the overlap
between the $U(z)$ and $V(z)$ profiles becomes negligible as $\gamma$
increases.

Given our previous discussion in Section 2.1 where we observed that setting
$\mathcal{V}=0$ means that the evolution equation for $u(x,t)$ simplifies to
the Fisher-KPP model (9), we anticipate that the solution of the dynamical
system associated with travelling wave solutions of the invasion model for
$\gamma\gg 1$ can be approximated by the solution of the dynamical system
associated with travelling wave solutions of the Fisher-KPP model,

$\dfrac{\mathrm{d}^{2}U}{\mathrm{d}z^{2}}+c\dfrac{\mathrm{d}U}{\mathrm{d}z}+U(1-U)=0,$
(21)

with boundary conditions $U\rightarrow 1$ as $z\rightarrow-\infty$, and
$U\rightarrow 0$ as $z\rightarrow\infty$. In the usual way, this second-order
boundary value problem can be re-written in terms of a first order system

$\displaystyle\dfrac{\mathrm{d}U}{\mathrm{d}z}$ $\displaystyle=W,$ (22)
$\displaystyle\dfrac{\mathrm{d}W}{\mathrm{d}z}$ $\displaystyle=-cW-U(1-U),$
(23)

with boundary conditions $U\rightarrow 1$ and $W\rightarrow 0$ as
$z\rightarrow-\infty$, and $U\rightarrow 0$ and $W\rightarrow 0$ as
$z\rightarrow-\infty$. There are two equilibrium points: (i)
$(\bar{U},\bar{W})=(1,0)$ that is associated with the invaded region; and,
(ii) $(\bar{U},\bar{W})=(0,0)$ that is associated with the uninvaded region.

To demonstrate the relationship between the three-dimensional dynamical system
and phase space for the invasion model and the two-dimensional dynamical
system and phase plane for the simpler Fisher-KPP model, we compare various
solutions for $\gamma=10$ and $1000$ with $\mathcal{V}=0.5$ in Figure 4, and a
separate set of comparisons are made for $\gamma=10$ and $1000$ with
$\mathcal{V}=1$ in Figure 5.

Figure 4: Travelling wave solutions of the invasion model with $\gamma\gg 1$
and $\mathcal{V}<1$. (a) and (d) show late-time numerical solutions of
Equations (3)–(4) with $\gamma=10$ and $\gamma=1000$, respectively. The
solution for $U(z)$ is shown in solid brown and the solution for $V(z)$ is
shown in solid blue. The approximate solution for $U(z)$ obtained from the
Fisher-KPP phase plane is shown in dashed green. (b) and (c) show a projection
of the $(U,V,W)$ phase space for the invasion model projected onto the $(U,W)$
plane together with the projection of the three-dimensional heteroclinic orbit
in solid brown. The trajectory from Fisher-KPP phase plane is shown in dashed
green. (c) and (f) show magnified regions near the origin in (b) and (e),
respectively.

Results in Figure 4(a) and (d) show the long-time numerical solutions of the
invasion model (3)–(4) with $\gamma=10$ and $\gamma=1000$, respectively.
Comparing the shapes of these two travelling wave solutions confirms that the
width of region II decreases with $\gamma$. In particular, the travelling wave
profile in Figure 4(d) confirms that the overlap between the invading cancer
density and the retreating skin density is barely noticeable at this scale.
These travelling wave profiles in Figure 4(a) and (d) are first generated and
plotted in the three-dimensional $(U,V,W)$ phase space, and a projection of
this phase space and the trajectory is plotted in the $(U,W)$ plane. This
projection looks like a two-dimensional heteroclinic orbit between
$(\bar{U},\bar{W})=(1,0)$ and $(\bar{U},\bar{W})=(0,0)$. To make the
connection with the simpler Fisher-KPP model explicit, we superimpose a
numerical trajectory obtained from the Fisher-KPP phase plane (22)–(23) for
the appropriate value of $c$ obtained from the long-time numerical solutions
of (3)–(4). In both cases we see that the projection of the trajectory
associated with the invasion model and the trajectory associated with the
simpler Fisher-KPP model compare very well, particularly in Figure 4(e) where
$\gamma=1000$. The main discrepancy between the trajectories is near the
origin. Additional comparisons in Figure 4(c) and (f) to show details of the
trajectories near the origin where the differences are clear. Indeed, in both
cases we see that the trajectory for the simpler Fisher-KPP model spirals into
the origin, whereas the trajectory for the invasion model does not [4].

Overall, the results in Figure 4 reveal a very interesting and novel
application of the well-known phase plane for travelling wave solutions of the
Fisher-KPP model since we show that trajectories in this phase plane provide a
good approximation to projections of the three-dimensional trajectories
associated with travelling wave solutions of the more complicated invasion
model (3)–(4). This comparison is mathematically interesting since standard
phase plane analysis of travelling wave solutions to the Fisher-KPP model are
limited to $c\geq 2$, given, the heteroclinic trajectories for $c<2$ are
disregarded on physical grounds [4]. In contrast, here we find that these
normally disregarded trajectories provide a good approximation for the shape
of travelling wave solutions to the more complicated invasion model. To
highlight the value of this approximation, we take the $(U,W)$ trajectory from
the Fisher-KPP phase plane and plot the profile as a function of $z$ in Figure
4(a) and (d), where we see that the profile obtained from the simpler Fisher-
KPP model approximates the shape of travelling wave profile for the full
invasion model. This approximation is particularly accurate in Figure 4(d) for
$\gamma=1000$. In contrast, while the Fisher-KPP approximation in Figure 4(a)
is quite reasonable where $z<2$, it is relatively poor in the region $2<z<5$
because of the more pronounced oscillation.

Results in Figure 5 are presented in the exact same format as the results in
Figure 4, except that here we have $\mathcal{V}=1$ so we have sharp-fronted
travelling wave solutions. Comparing the shape of the travelling wave
solutions in Figure 5(a) and (d) again confirms that the width of the overlap
region decreases with $\gamma$. The projections of the three-dimensional phase
space onto the $(U,W)$ plane in Figure 5(b) and (e) take the form of a two-
dimensional heteroclinic orbit between $(\bar{U},\bar{W})=(1,0)$ and
$(\bar{U},\bar{W})=(0,0)$. Comparing the projections of the three-dimensional
trajectory with the numerical trajectory obtained from the Fisher-KPP model
(22)–(23), confirms that the simpler dynamical system provides a good
approximation to the three-dimensional dynamical system when $\gamma$ is
large, with only small discrepancies near the origin, as shown in Figure 5(c)
and (f). In this case, comparing the $U(z)$ profiles in Figure 5(d) shows that
the entire shape of the travelling wave is approximated very well when
$\gamma=1000$, but we see a more clear discrepancy in Figure 5(a) for
$\gamma=10$ since this leads to $c=0.68$ and more pronounced oscillations
about $U=0$ at the front of the travelling wave.

Figure 5: Travelling wave solutions of the invasion model with $\gamma\gg 1$
and $\mathcal{V}=1$. (a) and (d) show late-time numerical solutions of
Equations (3)–(4) with $\gamma=10$ and $\gamma=1000$, respectively. The
solution for $U(z)$ is shown in solid brown and the solution for $V(z)$ is
shown in solid blue. The approximate solution for $U(z)$ obtained from the
Fisher-KPP phase plane is shown in dashed green. (b) and (c) show a projection
of the $(U,V,W)$ phase space for the invasion model projected onto the $(U,W)$
plane together with the projection of the three-dimensional heteroclinic orbit
in solid brown. The trajectory from Fisher-KPP phase plane is shown in dashed
green. (c) and (f) show magnified regions near the origin in (b) and (e),
respectively.

Overall, results in Figures 4 and 5 show that the well-known phase plane
associated with travelling wave solutions of the Fisher-KPP model are related
to travelling wave solutions of Browning’s invasion model (3)–(4). This is
mathematically interesting because standard analysis of travelling wave
solutions of the Fisher-KPP model disregards the phase plane for $c<2$,
whereas here we find that this previously disregarded phase plane is closely
related to the model of invasion. While some previous studies have shown that
this typically disregarded phase plane plays a role in our understanding of
certain moving boundary problems [48, 49, 50, 51, 52], this is the first time
that this previously disregarded phase plane has been shown to play a role in
our understanding of mathematical models of cellular invasion.

###  4.2 Slow decay: $\gamma\ll 1$

We now turn our attention to the limit where cancer cells consume skin cells
very slowly, $\gamma\ll 1$. In this limit we find that it is necessary to
treat the cases $\mathcal{V}<1$ and $\mathcal{V}=1$ separately, as we will now
illustrate.

####  4.2.1 $\gamma\ll 1$ and $\mathcal{V}<1$

Preliminary numerical results in Figure 2(a)–(d) show that the width of region
II, the overlap region, increases as $\gamma\rightarrow 0$. To analyse this
behaviour we seek a perturbation solution by treating $\gamma$ as a small
parameter, and it is useful to recall from Section 2 that when $\mathcal{V}<1$
we have $c=c_{\textrm{min}}=2(1-\mathcal{V})$ for $\gamma<\gamma_{\textrm{c}}$
so that we take $c=c_{\textrm{min}}$ in our small $\gamma$ analysis. Re-
scaling the dependent variable $\tilde{z}=\gamma z$ gives

$\displaystyle\gamma^{2}\dfrac{\mathrm{d}}{\mathrm{d}\tilde{z}}\left[(1-V)\dfrac{\mathrm{d}U}{\mathrm{d}\tilde{z}}\right]+c_{\textrm{min}}\gamma\dfrac{\mathrm{d}U}{\mathrm{d}\tilde{z}}+U(1-U-V)$
$\displaystyle=0,$ $\displaystyle-\infty<\tilde{z}<\infty,$ (24)
$\displaystyle c_{\textrm{min}}\dfrac{\mathrm{d}V}{\mathrm{d}\tilde{z}}-UV$
$\displaystyle=0,$ $\displaystyle-\infty<\tilde{z}<\infty.$ (25)

We now seek perturbation solutions of the form

$\displaystyle U(\tilde{z})$ $\displaystyle=U_{0}(\tilde{z})+\gamma
U_{1}(\tilde{z})+\mathcal{O}(\gamma^{2}),$ (26) $\displaystyle V(\tilde{z})$
$\displaystyle=V_{0}(\tilde{z})+\gamma
V_{1}(\tilde{z})+\mathcal{O}(\gamma^{2}).$ (27)

Substituting these expansions into Equation (24) shows that we have
$U_{0}(1-U_{0}-V_{0})=0$ at leading order, so that
$V_{0}(\tilde{z})+U_{0}(\tilde{z})=1$. The differential equations governing
the terms in the perturbation solution are therefore given by

$\displaystyle
c_{\textrm{min}}\dfrac{\mathrm{d}U_{0}}{\mathrm{d}\tilde{z}}+U_{0}(1-U_{0})$
$\displaystyle=0,$ (28) $\displaystyle
c_{\textrm{min}}\dfrac{\mathrm{d}U_{0}}{\mathrm{d}\tilde{z}}-U_{0}(U_{1}+V_{1})$
$\displaystyle=0,$ (29) $\displaystyle
c_{\textrm{min}}\dfrac{\mathrm{d}V_{1}}{\mathrm{d}\tilde{z}}-V_{1}(2U_{0}-1)+V_{0}^{2}$
$\displaystyle=0.$ (30)

with boundary conditions $U_{0}\rightarrow 1,V_{0}\rightarrow
0,U_{1}\rightarrow 0,V_{1}\rightarrow 0$ as $\tilde{z}\rightarrow-\infty$,
$U_{0}\rightarrow 0,V_{0}\rightarrow\mathcal{V},U_{1}\rightarrow
0,V_{1}\rightarrow 0$ as $\tilde{z}\rightarrow\infty$.

We solve these differential equations using the following strategy. Equation
(28) can be solved for $U_{0}(\tilde{z})$ directly using separation of
variables, and from this we can evaluate
$V_{0}(\tilde{z})=1-U_{0}(\tilde{z})$. Given the $\mathcal{O}(1)$ solutions,
we simply rearrange Equation (29) to obtain $U_{1}(z)$ directly, and the
solution of Equation (29) is obtained using an integrating factor. In summary,
the solutions of these differential equations are

$\displaystyle U_{0}(\tilde{z})$
$\displaystyle=\dfrac{1}{1+\left(\dfrac{\mathcal{V}}{1-\mathcal{V}}\right)\exp{\left(\dfrac{\tilde{z}}{c_{\textrm{min}}}\right)}},\quad
V_{0}(\tilde{z})=\dfrac{1}{1+\left(\dfrac{1-\mathcal{V}}{\mathcal{V}}\right)\exp{\left(-\dfrac{\tilde{z}}{c_{\textrm{min}}}\right)}},$
(31) $\displaystyle U_{1}(\tilde{z})$
$\displaystyle=\dfrac{-\mathcal{V}\exp{\dfrac{\tilde{z}}{c_{\textrm{min}}}}}{\left(\mathcal{V}\left[\exp{\left(\dfrac{\tilde{z}}{c_{\textrm{min}}}\right)}-1\right]+1\right)^{2}},\quad
V_{1}(\tilde{z})=\dfrac{\mathcal{V}^{2}\left[1-\exp{\left(\dfrac{\tilde{z}}{c_{\textrm{min}}}\right)}\right]\exp{\left(-\dfrac{\tilde{z}}{c_{\textrm{min}}}\right)}}{\left[\exp{\left(-\dfrac{\tilde{z}}{c_{\textrm{min}}}\right)}(\mathcal{V}-1)-\mathcal{V}\right]^{2}}.$
(32)

where we have evaluated the constants of integration by setting
$V(0)=\mathcal{V}/2$.

Results in Figure 6 show long-time solutions of (3)–(4) in (a), (d) and (g)
for $\gamma=0.1,0.01$ and $0.005$, respectively. Comparing the shapes of these
travelling waves confirms that the width of the overlap region increases as
$\gamma$ decreases. Results in Figure 6(b), (e) and (h) compare the shape of
the late-time PDE solutions with the $\mathcal{O}(\gamma)$ perturbation
solutions, and we see that the accuracy of the approximate perturbation
solutions improves as $\gamma$ decreases, as expected. Results in Figure 6(c),
(f) and (i) compare the numerical solutions and the perturbation solutions
within the regions highlighted by the dashed rectangles in Figure 6(b), (e)
and (h). Again, we see the accuracy of the perturbation solution increases as
$\gamma$ decreases, and the perturbation solution captures the sharp
transition region reasonably accurately as $\gamma\to 0$.

Figure 6: Comparison of PDE solutions with perturbation solution when
$\gamma\ll 1$ and $\mathcal{V}=0.5$. Travelling wave solutions $U(z)$ and
$V(z)$ obtained from PDE with $\Delta x=\Delta t=1\times 10^{-2}$,
$\gamma=0.1$ in (a)–(c), $\gamma=0.01$ in (d)–(f) and $\gamma=0.005$ in
(g)–(i) are illustrated in solid brown and solid blue for $U(z)$ and $V(z)$
respectively. $\mathcal{O}(\varepsilon)$ perturbation solution is illustrated
in dashed yellow respectively for $U(z)$ and in dashed red for $V(z)$.
Solutions shown in (c),(f) and (i) are magnification of region of interest in
(b),(e) and (h) respectively.

###  4.3 $\gamma\ll 1$ and $\mathcal{V}=1$

To analyse the shape of the travelling wave for $\gamma\ll 1$ with
$\mathcal{V}=1$ we consider the desingularised system

$\displaystyle\dfrac{\mathrm{d}^{2}U}{\mathrm{d}\zeta^{2}}+c\dfrac{\mathrm{d}U}{\mathrm{d}\zeta}+U(1-V)(1-U-V)$
$\displaystyle=0,$ (33)
$\displaystyle\dfrac{\mathrm{d}V}{\mathrm{d}\zeta}-\dfrac{\gamma UV}{c}(1-V)$
$\displaystyle=0,$ (34)

with boundary conditions $U\rightarrow 1,V\rightarrow 0$ as
$\zeta\rightarrow-\infty$, $U\rightarrow 0,V\rightarrow 1$ as
$\zeta\rightarrow\infty$.

Numerical results in Figure 2 show that when $\mathcal{V}=1$ we have
$c\rightarrow 0$ as $\gamma\rightarrow 0$, so we write $c=\tilde{c}\gamma$ so
that $\tilde{c}$ is $\mathcal{O}(1)$. Like in the previous section, we re-
scale the independent variable $\tilde{\zeta}=\gamma\zeta$ to give

$\displaystyle\dfrac{\mathrm{d}^{2}U}{\mathrm{d}\tilde{\zeta}^{2}}+\tilde{c}\gamma\dfrac{\mathrm{d}U}{\mathrm{d}\tilde{\zeta}}+U(1-U-V)(1-V)$
$\displaystyle=0,$ $\displaystyle-\infty<\tilde{\zeta}<\infty$ (35)
$\displaystyle\tilde{c}\dfrac{\mathrm{d}V}{\mathrm{d}\tilde{\zeta}}-UV(1-V)$
$\displaystyle=0,$ $\displaystyle-\infty<\tilde{\zeta}<\infty,$ (36)

with boundary conditions $U\rightarrow 1,V\rightarrow 0$ as
$\tilde{\zeta}\rightarrow-\infty$, $U\rightarrow 0,V\rightarrow 1$ as
$\tilde{\zeta}\rightarrow\infty$.

Seeking a perturbation solution of the form

$\displaystyle U(\tilde{\zeta})$ $\displaystyle=U_{0}(\tilde{\zeta})+\gamma
U_{1}(\tilde{\zeta})+\mathcal{O}(\gamma^{2}),$ (37) $\displaystyle
V(\tilde{\zeta})$ $\displaystyle=V_{0}(\tilde{\zeta})+\gamma
V_{1}(\tilde{\zeta})+\mathcal{O}(\gamma^{2}),$ (38)

leads to a system of coupled differential equations for the leading order
terms,

$\displaystyle\dfrac{\mathrm{d}^{2}U_{0}}{\mathrm{d}\tilde{\zeta}^{2}}+U_{0}(1-U_{0}-V_{0})(1-V_{0})$
$\displaystyle=0,$ (39)
$\displaystyle\tilde{c}\dfrac{\mathrm{d}V_{0}}{\mathrm{d}\tilde{\zeta}}-U_{0}V_{0}(1-V_{0})$
$\displaystyle=0,$ (40)

with boundary conditions $U_{0}\rightarrow 1,V_{0}\rightarrow 0$ as
$\tilde{\zeta}\rightarrow-\infty$, $U_{0}\rightarrow 0,V_{0}\rightarrow 1$ as
$\tilde{\zeta}\rightarrow\infty$. Unfortunately we are unable to obtain a
closed-form solution of Equations (39)–(40) and we do not proceed further with
this approach.

##  5 Fast travelling waves, $c>c_{\textrm{min}}$

Our focus so far has been on the long-time limit of the time-dependent
solutions of (3)–(4) with initial conditions given by Equations (5)–(6) so
that $u(x,0)$ has compact support. This leads to travelling wave solutions
with the minimum wave speed, $c_{\textrm{min}}$. We now examine travelling
wave solutions of the same model but with initial conditions given by
Equations (7)–(8), where $a>0$ controls the far-field decay rate of $u(x,0)$.
Results in Figure 7 summarise the numerically-observed travelling wave speed,
$c$, for $a=0.1,0.25,0.5$ and $1$, as a function of $\mathcal{V}$ and $\gamma$
as indicated. Comparing these results with those in Figure 3 for initial
conditions with compact support, we see that some general features of the
relationship between $c$, $\mathcal{V}$ and $\gamma$ are maintained, while
other features are different. In general we see that $c$ is a decreasing
function of $a$, and that all results suggest that $c$ is independent of
$\gamma$ for sufficiently small $\gamma$ when $\mathcal{V}<1$. As $\gamma$
increases, we see that $c$ increases with $\gamma$, but that the limiting
value of $c$ as $\gamma\to\infty$ depends upon the decay rate, $a$. Careful
comparison of the results in Figure 7(d) for $a=1$ shows that different
choices of $\gamma$ and $\mathcal{V}$ lead to the exact same travelling wave
speed as in Figure 3 for initial conditions with compact support, which can be
thought of as letting $a\to\infty$ in (7)–(8). We will now explain some of
these observed trends analytically.

Figure 7: Relationship between $c$, $\gamma$ and $\mathcal{V}$ for slowly
decaying initial $u(x,0)$. Numerical estimates of $c$ are obtained from long-
time solutions of Equations (3)–(4) with the initial condition given by
Equations (7)–(8). Results are presented with $\beta=10$ and $a=0.1,0.25,0.5$
and $1$, as indicated. Time-dependent PDE solutions are obtained using $\Delta
x=\Delta t=1\times 10^{-2}$, for $1\times 10^{-3}\leq\gamma\leq 1\times
10^{6}$ and $\mathcal{V}=1/8,2/8,3/8,4/8,5/8,6/8,7/8$ and $1$.

To understand the relationship between the decay rate of the initial condition
and the asymptotic wave speed, $c$, we examine the leading edge of the
travelling wave where $u\ll 1$, giving

$\displaystyle\dfrac{\partial\tilde{u}}{\partial
t}=(1-\tilde{v})\dfrac{\partial^{2}\tilde{u}}{\partial
x^{2}}+(1-\tilde{v})\tilde{u},$ (41)
$\displaystyle\dfrac{\partial\tilde{v}}{\partial
t}=-\gamma\tilde{u}\tilde{v}.$ (42)

Assuming the travelling wave solution takes the form
$\tilde{u}\sim\textrm{exp}[-a(x-ct)]$ for large $x$, substituting this into
Equation (41) gives

$c=(1-\mathcal{V})\left(a+\dfrac{1}{a}\right).$ (43)

This dispersion relationship is similar to the analogous result for the
Fisher-KPP model [4]. This simple relationship explains some of the
observations in Figure 7 where we see that $c$ is independent of $\gamma$ for
$\gamma<\gamma_{\textrm{c}}$. Indeed, this constant speed is given by Equation
(43). Unfortunately, this dispersion relationship does not provide any insight
into the relationship between $c$, $\mathcal{V}$ and $\gamma$ for
$\gamma>\gamma_{\textrm{c}}$, nor any insight into the shape of the resulting
travelling waves.

To address the shape of the travelling waves for $c>c_{\textrm{min}}$ we
extend the approach of Canosa [24] who noted that travelling wave solutions of
the Fisher-KPP model become increasingly wide as $c\to\infty$. Our numerical
simulation results suggest that travelling wave solutions of the invasion
model behave similarly, so we explore this behaviour by re-scaling the
independent variable, $\overline{z}=z/c$, giving

$\displaystyle\dfrac{\mathrm{d}U}{\mathrm{d}\overline{z}}+\frac{1}{c^{2}}\dfrac{\mathrm{d}}{\mathrm{d}\overline{z}}\left[(1-V)\dfrac{\mathrm{d}U}{\mathrm{d}\overline{z}}\right]+U(1-U-V)$
$\displaystyle=0,\quad-\infty<\overline{z}<\infty,$ (44)
$\displaystyle\dfrac{\mathrm{d}V}{\mathrm{d}\overline{z}}-\gamma UV$
$\displaystyle=0,\quad-\infty<\overline{z}<\infty,$ (45)

with boundary conditions $U(\overline{z})\rightarrow 1$ and
$V(\overline{z})\rightarrow 0$ as $\overline{z}\rightarrow-\infty$, and
$U(\overline{z})\rightarrow 0$ and $V(\overline{z})\rightarrow\mathcal{V}$ as
$\overline{z}\rightarrow\infty$. In this re-scaled coordinate we seek
perturbation solutions in terms of the small parameter $\varepsilon=1/c^{2}$,

$\displaystyle U(\overline{z})$ $\displaystyle=U_{0}(\overline{z})+\varepsilon
U_{1}(\overline{z})+\mathcal{O}(\varepsilon^{2}),$ (46) $\displaystyle
V(\overline{z})$ $\displaystyle=V_{0}(\overline{z})+\varepsilon
V_{1}(\overline{z})+\mathcal{O}(\varepsilon^{2}).$ (47)

Substituting Equations (46)–(47) into (44)–(45) leads to

$\displaystyle\dfrac{\mathrm{d}U_{0}}{\mathrm{d}\overline{z}}+U_{0}(1-U_{0}-V_{0})$
$\displaystyle=0,$ (48)
$\displaystyle\dfrac{\mathrm{d}V_{0}}{\mathrm{d}\overline{z}}-\gamma
U_{0}V_{0}$ $\displaystyle=0,$ (49)

with boundary conditions $U_{0}(\overline{z})\rightarrow 1$ and
$V_{0}(\overline{z})\rightarrow 0$ as $\overline{z}\rightarrow-\infty$,
$U_{0}(\overline{z})\rightarrow 0$ and
$V_{0}(\overline{z})\rightarrow\mathcal{V}$ as
$\overline{z}\rightarrow\infty$. Unlike Canosa [24], these differential
equations for the $\mathcal{O}(1)$ terms in the perturbation solution do not
have a closed-form solution. Nevertheless, we make progress by re-writing
Equations (48)–(49) as

$\dfrac{\mathrm{d}U_{0}}{\mathrm{d}V_{0}}=-\dfrac{1-U_{0}-V_{0}}{\gamma
V_{0}},$ (50)

with $U_{0}(\mathcal{V})=0$. The solution of this problem is given by

$U_{0}(V_{0})=\dfrac{1}{(\gamma-1)}\left[\left(V_{0}+\gamma-1\right)-({\mathcal{V}}+\gamma-1)\left(\dfrac{V_{0}}{\mathcal{V}}\right)^{\left(\dfrac{1}{\gamma}\right)}\right].$
(51)

For the special case $\gamma=1$, this solution can be written as

$U_{0}(V_{0})=1+\left(\dfrac{V_{0}}{\mathcal{V}}\right)\left(\mathcal{V}\log\left[\dfrac{V_{0}}{\mathcal{V}}\right]-1\right).$
(52)

Results in Figure 8 compare travelling wave solutions with various $c$ with
the $\mathcal{O}(1)$ perturbation solution (51). Results in the left column of
Figure 8 show $U$ as a function of $V$, as explicitly defined by Equation (51)
superimposed on curves obtained from long-time numerical solutions of
Equations (3)–(4) that are plotted in the same format. Results in the right
column of Figure 8 show the late time solutions of Equations (3)–(4) plotted
in the travelling wave coordinate, $z$, superimposed with the perturbation
solution. Results in Figure 8(a)–(d) show perturbation results for $c=5.1$ and
different choices of $\gamma$. At this scale the perturbation solution is
visually indistinguishable from the numerical solutions. Results in Figure
8(e)–(f) show that the perturbation solution performs well for $c=2$, but that
we can begin to see some discrepancy between the perturbation and numerical
results. Interestingly, the perturbation solution in Figure 8(g)–(h) leads to
reasonably accurate solutions for $c=1$, despite the fact that Equation (51)
is valid in the limit $c\to\infty$.

Figure 8: Comparison of numerical and perturbation solutions for $c\geq
c_{\textrm{cmin}}$. The $\mathcal{O}(1)$ perturbation solution shown in dashed
green is compared to solution $U(V)$ from PDE shown in solid pink, in (a),
(c), (e) and (g). Solutions $U(z)$ and $V(z)$ from PDE are shown in solid
brown and in solid blue respectively in (b), (d), (f) and (h). Perturbation
solution is shown in dashed yellow as $U(z)$ by using $V_{0}=V(z)$ from PDE
solution. The initial conditions (7)–(8) with $\beta=10$ are used in (a)–(d),
where $a=0.1$, and in (e)–(f) where $a=0.27$. The initial conditions with
compact support are used in (g)–(h).

##  6 Conclusion and Outlook

In this work we study travelling wave solutions of a model of cellular
invasion, (3)–(4), where the migration and proliferation of the invasive
population is coupled to the degradation of surrounding skin tissues [1].
Time-dependent numerical solutions of the governing PDEs show that there is a
complicated relationship between the travelling wave speed $c$ and: (i)
$\gamma$, the rate at which the invasive population degrades the surrounding
skin tissues; and, (ii) $\mathcal{V}$, the far field density of surrounding
tissues. Numerical exploration shows that long-time travelling wave solutions
are smooth without compact support for $\mathcal{V}<1$, or sharp-fronted with
compact support with $\mathcal{V}=1$. The relationship between $c$ and the
parameters in the model are partially established in this work. Numerical
simulations and phase space analysis show that for initial conditions with
compact support, we have $c=2(1-\mathcal{V})$, which is independent of
$\gamma$ for $\gamma<\gamma_{\textrm{c}}$. Further numerical simulations show
that $c$ increases with $\gamma$ for $\gamma>\gamma_{\textrm{c}}$, with $c\to
2^{-}$ as $\gamma\to\infty$, but the precise details of this relationship are
not revealed through standard phase space analysis. Of great interest is the
fact that we always have $c<2$ for the invasion model for initial conditions
with compact support; this is very different to the standard Fisher-KPP model
where $c\geq 2$. We also show that time-dependent PDE solutions for initial
conditions without compact support lead to travelling wave solutions with
larger wave speed.

Our analysis of the invasion model for fast decay, $\gamma\gg 1$, indicates
that the width of the overlap region decreases with $\gamma$. This means that
the density of the invading population becomes uncoupled from the density of
the surrounding skin tissues, and suggests that the shape of the invading
density profile is related to the shape of the travelling wave solution of the
Fisher-KPP model. This is intriguing since the invasion model is associated
with travelling waves with $c<2$ whereas the Fisher-KPP model is associated
with $c\geq 2$. Indeed, our results show that the well-known phase plane
associated with travelling wave solutions of the Fisher-KPP model provides a
novel approximation to the shape of the travelling wave solution of the
invasion model for fast decay. This observation is mathematically interesting
because standard analysis of the Fisher-KPP model disregards the phase plane
for $c<2$ [4], whereas here we find that this previously disregarded phase
plane is closely related to our model of invasion.

The observation that the normally-disregarded phase plane of the Fisher-KPP
model for $c<2$ plays a role in the study of coupled reaction-diffusion models
of invasion provides a very interesting and previously unnoticed link with a
very different type of mathematical model of invasion, the Fisher-Stefan model
[48, 49, 50, 51, 52, 53, 46]. The Fisher-Stefan model involves studying the
Fisher-KPP model on a moving boundary, $0<x<L(t)$. In this model the density
vanishes on the moving boundary, $u(x,L(t))=0$, and the speed of the moving
boundary is driven by a one-phase Stefan condition,
$\textrm{d}L(t)/\textrm{d}t=-\kappa\partial u(L(t),t)/\partial x$. Here,
$\kappa>0$ is a constant that relates the speed of the moving boundary to the
spatial gradient of density at the moving boundary. For both the invasion
model and the Fisher-Stefan model it has been shown that time-dependent PDE
solutions eventually evolve to travelling wave solutions with $c<2$ [48, 49,
50, 51, 52, 53, 46] and the shape of the invading profiles in both cases is
given by the normally-disregarded phase plane of the well-known Fisher-KPP
model. This is very interesting because the normally-disregarded phase plane
trajectories imply $U(z)<0$ for certain intervals in $z$. However, here and in
the Fisher-Stefan model, the profile of interest is given by a truncated
trajectory in the phase plane where $U(z)>0$.

There are many opportunities to extend the work presented in this study. There
are several assumptions in the mathematical model (3)–(4), that could be
relaxed or varied. Such extensions could involve working with different
nonlinear diffusivity functions in (3), such as a power law [8]. Another
extension of interest would be to explore the impact of using a different
nonlinear source term in (3) to model the proliferation of cells [5]. One of
the limitations of our study is that we have not been able to arrive at a
mathematical expression for $\gamma_{\textrm{c}}$, and it would be of great
interest to arrive at some approximate expression for this critical decay
rate, or to place some bound on that value. Further extensions could involve
working in a different geometry since models of melanoma invasion in both two
and three-dimensions with radial symmetry are of great interest for studying
malignant invasion, e.g [43, 54, 55] .

##  Appendix A Numerical methods

We solve Equations (3)–(4) on $0<x<L$ by uniformly discretising the domain
with $N$ mesh points, with spacing $\Delta x$. We approximate the spatial
derivatives in Equations (3)–(4) using a central difference approximation, and
solve the resulting system of coupled ordinary differential equation through
time using an implicit Euler approximation, giving

$\displaystyle\dfrac{u_{i}^{j+1}-u_{i}^{j}}{\Delta t}=$
$\displaystyle\dfrac{1}{2\Delta
x^{2}}\bigg{(}\left[2-\left(v_{i+1}^{j+1}+v_{i}^{j+1}\right)\right]\left(u_{i+1}^{j+1}-u_{i}^{j+1}\right)-\left[2-\left(v_{i}^{j+1}+v_{i-1}^{j+1}\right)\right]\left(u_{i}^{j+1}-u_{i-1}^{j+1}\right)\bigg{)}$
$\displaystyle+u_{i}^{j+1}\left(1-u_{i}^{j+1}-v_{i}^{j+1}\right),$ (53)
$\displaystyle\dfrac{v_{i}^{j+1}-v_{i}^{j}}{\Delta t}$ $\displaystyle=-\gamma
u_{i}^{j+1}v_{i}^{j+1},$ (54)

where $\Delta t$ is the time step, $i$ is the spatial finite difference mesh
index and $j$ is the temporal index so that $u_{i}^{j}\approx u(x=(i-1)\Delta
x,j\Delta t)$. Discretising the boundary conditions for $u(x,t)$ leads to

$u_{2}^{j+1}-u_{1}^{j+1}=0,\quad u_{N}^{j+1}-u_{N-1}^{j+1}=0.$ (55)

Note that there are no boundary conditions for $v$, so the spatial index on
the discretisation for $v$ is $i=1,2,\ldots,N$. This discretisation leads to a
coupled system of nonlinear algebraic equations for $u_{i}^{j+1}$ and
$v_{i}^{j+1}$, which are solved sequentially to take advantage of the
tridiagonal structure of the discretised equations. The nonlinear equations
are solved using Newton-Raphson iterations that continue until the maximum
change in the dependent variables falls below some tolerance $\epsilon$ in
each time step. For all problems considered we always check that our choices
of $\Delta x$, $\Delta t$ and $\epsilon$ lead to grid-independent results.
Matlab software to implement these numerical solutions is available on GitHub.

To estimate the travelling wave speed from our time-dependent PDE solutions we
specify a contour value, $u(x,t)=U$. At the end of each time step in we use
linear interpolation to find the value of $X$ such that $u(X,t)=U$. At the end
of each time step we have estimates of both $X(t+\Delta t)$ and $X(t)$,
allowing us to estimate the speed at which the contour moves

$c=\dfrac{X(t+\Delta t)-X(t)}{\Delta t}.$ (56)

Evaluating Equation (56) at each time step leads to a time series of estimates
of $c$, and we find that these estimates asymptote to some positive constant
value as $t$ becomes sufficiently large. For all results presented we set
$U=0.5$, but we find that our estimates of $c$ are not particularly sensitive
to our choice of $U$ provided that our choice of $U$ is reasonably close to
the middle of the interval, $0<U<1$.

To construct phase planes for Fisher-KPP equation, we solve Equation (21)
numerically using Heun’s method with a constant step size $\textrm{d}z$. Since
we are interested in examining trajectories that leave the saddle $(1,0)$
along the unstable manifold we choose the initial location on the trajectory
to be on the appropriate unstable manifold and sufficiently close to $(1,0)$.
Matlab software to generate these phase planes and associated trajectories are
available on GitHub.

##### Acknowledgements

This work is supported by the Australian Research Council (DP200100177).

##### Contributions

All authors conceived and designed the study; M.El-H. performed numerical and
symbolic calculations. All authors drafted the article and gave final approval
for publication.

##### Competing Interests

We have no competing interests.

## References

  * [1] A.P. Browning, P. Haridas, M.J. Simpson. A Bayesian sequential learning framework to parameterise continuum models of melanoma invasion into human skin. Bulletin of Mathematical Biology, 81 (2019) 676–698.
  * [2] R.A. Fisher. The wave of advance of advantageous genes. Annals of Eugenics, 7 (1937) 355–369.
  * [3] A.N. Kolmogorov, P.G. Petrovskii, N.S. Piskunov. A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem. Moscow University Mathematics Bulletin, 1 (1937) 1–26.
  * [4] J.D. Murray. Mathematical Biology I: An Introduction. Third edition, Springer, New York, (2002).
  * [5] F. Sánchez Garduño, P.K. Maini. An approximation to a sharp type solution of a density-dependent reaction-diffusion equation. Applied Mathematics Letters, 7 (1994) 47–51.
  * [6] T.P. Witelski. An asymptotic solution for traveling waves of a nonlinear-diffusion Fisher’s equation. Journal of Mathematical Biology, 33 (1994) 1–16.
  * [7] T.P. Witelski. Merging traveling waves for the porous-Fisher’s equation. Applied Mathematics Letters, 8 (1995) 57–62.
  * [8] S.W. McCue, W. Jin, T.J. Moroney, K-Y. Lo, S-E. Chou, M.J. Simpson. Hole-closing model reveals exponents for nonlinear degenerate diffusivity functions in cell biology. Physica D: Nonlinear Phenomena, 398 (2019) 130–140.
  * [9] J.A. Sherratt, J.D. Murray. Models of epidermal wound healing. Proceedings of the Royal Society of London: Series B, 241 (1990) 29–36.
  * [10] K.R. Swanson, C. Bridge, J.D. Murray, E.C. Alvord Jr. Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion. Journal of Neurological Sciences, 216 (2003) 1–10.
  * [11] P.K. Maini, D.L.S. McElwain, D.I. Leavesley. Traveling wave model to interpret a wound-healing cell migration assay for human peritoneal mesothelial cells. Tissue Engineering, 10 (2004) 475–482.
  * [12] B.G. Sengers, C.P. Please, R.O.C. Oreffo. Experimental characterization and computational modelling of two-dimensional cell spreading for skeletal regeneration. Journal of the Royal Society Interface, 4 (2007) 1107–1117.
  * [13] P. Gerlee, S. Nerlander. The impact of phenotypic switching on glioblastoma growth and invasion. PLoS Computational Biology, 8(6) (2012) e1002556.
  * [14] W. Jin, K-Y. Lo, Y-S. Sun, YH Ting, M.J. Simpson. Quantifying the role of different surface coatings in experimental models of wound healing. Chemical Engineering Science 220 (2020) 115609.
  * [15] J.H. Lagergren, J.T. Nardini, G. Michael Lavigne, E.M. Rutter, K.B. Flores. Learning partial differential equations for biological transport models from noisy spatio-temporal data. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476 (2021) 20190800.
  * [16] J.H. Largergren, J.T. Nardini, R.E. Baker, M.J. Simpson, K.B. Flores. Biologically-informed neural networks guide mechanistic modelling from sparse experimental data. PLOS Compuational Biology, 16(12): e1008462.
  * [17] J.G. Skellam. Random dispersal in theoretical populations. Biometrika, 38 (1951) 196–218.
  * [18] M.A. Lewis, P. Kareiva. Allee dynamics and the spread of invading organisms. Theoretical Population Biology, 43 (1994) 141–158.
  * [19] E.E. Holmes, M.A. Lewis, J.E. Banks, R.R. Veit. Partial differential equations in ecology: spatial interactions and population dynamics. Ecology, 74 (1994) 17–29.
  * [20] N. Shigesada, K. Kawasaki, Y. Takeda. Modeling stratified diffusion in biological invasions. American Naturalist, 146 (1995) 229–251.
  * [21] J. Steele, J. Adams, T. Sluckin. Modelling paleoindian dispersals. World Archaeology, 30 (1998), 286–305.
  * [22] M. Kot. Elements of Mathematical Ecology. Cambridge University Press, Cambridge, (2003).
  * [23] S.A. Levin, H.C. Muller-Landau, R. Nathan, J. Chave. The ecology and evolution of seed dispersal: a theoretical perspective. Annual Review of Ecology, Evolution, and Systematics. 34 (2003) 575–604.
  * [24] J. Canosa. On a nonlinear diffusion equation describing population growth. IBM Journal of Research and Development, 17 (1973) 307–313.
  * [25] K.J. Painter, J.A. Sherratt. Modelling the movement of interacting cell populations. Journal of Theoretical Biology, 225 (2003) 327–339.
  * [26] H. Byrne, L. Preziosi. Modelling solid tumour growth using the theory of mixtures. Mathematical Medicine and Biology: A Journal of the IMA, 20 (2003) 341–366.
  * [27] H.M. Byrne, J.R. King, D.L.S. McElwain, L Preziosi. A two-phase model of solid tumour growth. Applied Mathematics Letters, 16 (2003) 567–573.
  * [28] P. Haridas, J.A. McGovern, D.L.S. McElwain, M.J. Simpson. Quantitative comparison of the spreading and invasion of radial growth phase and metastatic melanoma cells in a three-dimensional human skin equivalent model. PeerJ, 5 (2017) e3754.
  * [29] 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 proliferation drives melanoma nest formation in human skin tissue. BMC Systems Biology, 12 (2018) 34.
  * [30] R.A. Gatenby, E.T. Gawlinski. A reaction-diffusion model of cancer invasion. Cancer Research, 56 (1996) 5745–5753.
  * [31] K.A. Landman, G.J. Pettet. Modelling the action of proteinase and inhibitor in tissue invasion. Mathematical Biosciences, 154 (1998) 23–37.
  * [32] A.J. Perumpanani, J.A. Sherratt, J. Norbury, H.M. Byrne. A two parameter family of travelling waves with a singular barrier arising from the modelling of extracellular matrix mediated cellular invasion. Physica D: Nonlinear Phenomena, 126 (1999) 145–159.
  * [33] B.P. Marchant, J. Norbury, A.J. Perumpanani, Travelling shock waves arising in a model of malignant invasion. SIAM Journal on Applied Mathematics, 60 (2000) 463–476.
  * [34] K. Smallbone, D.J. Gavaghan, R.A. Gatenby, P.K. Maini. The role of acidity in solid tumour growth and invasion. Journal of Theoretical Biology, 235 (2005) 476–484.
  * [35] K.A. Landman, M.J. Simpson, G.J. Pettet, Tactically-driven nonmonotone travelling waves. Physica D: Nonlinear Phenomena, 237 (2008) 678–691.
  * [36] A.R.A. Anderson, V. Quaranta. Integrative mathematical oncology. Nature Reviews Cancer, 8 (2008) 227–234.
  * [37] S. Astanin, L. Preziosi. Mathematical modelling of the Warburg effect in tumour cords. Journal of Theoretical Biology, 258 (2009) 578–590.
  * [38] A. Fasano, M.A. Herrero, M.R. Rodrigo. Slow and fast invasion waves in a model of acid-mediated tumour growth. Mathematical Biosciences, 220 (2009) 45–56.
  * [39] H.M. Byrne. Dissecting cancer through mathematics: from the cell to the animal model. Nature Reviews Cancer, 10 (2010) 221–230.
  * [40] M.J. Tindall, L. Dyson, K. Smallbone, P.K. Maini. Modelling acidosis and the cell cycle in multicellular tumour spheroids. Journal of Theoretical Biology, 298 (2012) 107–115.
  * [41] A.B. Holder, M.R. Rodrigo, M.A. Herrero. A model for acid-mediated tumour growth with nonlinear acid production term. Applied Mathematics and Computation, 227 (2014) 176–198.
  * [42] A.B. Holder, M.R. Rodrigo. Model for acid-mediated tumour invasion with chemotherapy intervention II: Spatially heterogeneous populations. Mathematical Biosciences, 270 (2015) 10–29.
  * [43] J.P. Ward, J.R. King. Mathematical modelling of avascular tumour growth. Mathematical Medicine and Biology: A Journal of the IMA, 14 (1997) 39–69.
  * [44] H.M. Byrne, M.A.J. Chaplain. Free boundary value problems associated with the growth and development of multicellular spheroids. European Journal of Applied Mathematics, 8 (1997) 639-658.
  * [45] E.A. Gaffney, P.K. Maini, C.D. McCaig, M. Zhao, J.V. Forrester. Modelling corneal epithelial wound closure in the presence of physiological electric fields via a moving boundary formalism. Mathematical Medicine and Biology: A Journal of the IMA, 16 (1999) 369–393.
  * [46] M. El-Hachem, S.W. McCue, M.J. Simpson. Invading and receding sharp-fronted travelling waves. Bulletin of Mathematical Biology, 83 (2021) 35.
  * [47] S. Wiggins. Introduction to Applied Nonlinear Dynamical Systems and Chaos. Second edition, Springer, New York, (2003).
  * [48] Y. Du, Z. Lin. Spreading-vanishing dichotomy in the diffusive logistic model with a free boundary. SIAM Journal on Mathematical Analysis, 42 (2010) 377–405.
  * [49] Y. Du, Z. Guo. Spreading-vanishing dichotomy in a diffusive logistic model with a free boundary, II. Journal of Differential Equations, 250 (2011) 4336–4366.
  * [50] Y. Du, Z. Guo. The Stefan problem for the Fisher-KPP equation. Journal of Differential Equations, 253 (2012) 996–1035.
  * [51] Y. Du, H. Matsuzawa, M. Zhou. Sharp estimate of the spreading speed determined by nonlinear free boundary problems. SIAM Journal on Mathematical Analysis, 46 (2014) 375–396.
  * [52] 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 dichotomy. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences. 475 (2019) 20190378.
  * [53] M. El-Hachem, S.W. McCue, M.J. Simpson. A sharp-front moving boundary model for malignant invasion. Physica D: Nonlinear Phenomena, 412 (2020) 132639.
  * [54] J.P. Ward, J.R. King. Mathematical modelling of avascular-tumour growth II: modelling growth to saturation. Mathematical Medicine and Biology, 16 (1999) 171–211.
  * [55] W. Jin, L. Spoerri, N.K. Haass, M.J. Simpson. Mathematical model of tumour spheroid experiments with real-time cell cycle imaging. Bulletin of Mathematical Biology. 83 (2021) 44.