Travelling wave analysis of cellular invasion into surrounding tissues.

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

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 [50]. 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

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,

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 [32].

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 [32, 33]. Once we have characterised this first family of travelling wave solutions, we will then study another family of travelling wave solutions where $u(x,0)$ decays to zero as $x\to\infty$. While this second family of travelling wave solutions are mathematically interesting, the fact that these travelling wave solutions that arise from exponentially decaying initial conditions means that this second set of travelling waves are more difficult to interpret biologically since this kind of special initial condition is not relevant in practice. To study this second family of travelling wave solutions we work with an initial condition given by

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],

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],

We will now begin to explore these relationships.

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$.

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 since they are differentiable everywhere on the domain, and they do not have compact support since they decay to zero in the far field, $u(x,t)\to 0^{+}$ as $x\rightarrow\infty$. In this regard, the shape of these travelling wave solutions 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. These solutions are non-smooth since they contain a well defined front and $u(x,t)$ is not differentiable at the location of the front, $x=X$ (sometimes called the contact point). These solutions have compact support since $u(x,t)>0$ for $x<X$, and $u(x,t)=0$ for $x\geq X$. 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. The relationship between $c$, $\gamma$ and $\mathcal{V}$ is summarised 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 Equations (5)–(6. As we will establish later in Section 5, this means that these estimates of $c$ correspond to the minimum wave speed, $c_{\textrm{min}}$. Remembering that the Fisher-KPP model has a very simple minimum wave speed, $c_{\textrm{min}}=2$ [4, 24], here we see that the minimum wave speed for the invasion model is a relatively complicated function of $\mathcal{V}$ and $\gamma$ and we will devote much of this work to understanding this relationship in different limiting cases. 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$. Interestingly, for all choices of $\gamma$ and $\mathcal{V}>0$ we observe travelling wave solutions with a minimum wave speed that is less than the minimum travelling wave speed for the Fisher-KPP model, and we will return to this point later. 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.

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

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,

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

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.

Setting $\mathcal{V}<1$ leads to smooth travelling wave solutions that we will analyse in phase space. Since the travelling wave solutions are smooth the phase plane is non-singular. We refer to this as a traditional phase space since we do not have to consider any nonsingularities, and 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$ [51]. 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. Alternatively, simply set $\gamma=0$ in Equations (3)–(4) uncouples the system. Seeking travelling wave solutions of this uncoupled system leads to Fisher-KPP-like phase plane analysis, giving $c_{\textrm{min}}=2(1-\mathcal{V})$ which corroborates our previous result. Unfortunately, this simpler approach provides no insight when $\gamma>0$.

Biologically, setting $\mathcal{V}<1$ corresponds to the situation where the density of the skin tissue ahead of the invading front of is lower than the carrying capacity of the invading population. Intuitively, we might anticipate that the speed of the invading front would decrease with $\mathcal{V}$ and increase with $\gamma$. While the first of these intuitive expectations is consistent with our analysis and numerical observations, the second point highlights the value of our mathematical analysis and numerical explorations since our finding that the minimum wave speed is independent of $\gamma$, for sufficiently small $\gamma<\gamma_{\textrm{c}}$, is not at all obvious. The implication of this finding is that interventions seeking to reduce the decay rate to zero would not stop the invasion since $c_{\textrm{min}}>0$ when $\gamma=0$.

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

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

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$.

Biologically, setting $\mathcal{V}=1$ corresponds to the situation where the density of the skin tissue ahead of the invading front is identical to the carrying capacity of the invading population. In this situation we see that $c_{\textrm{min}}\to 0^{+}$ as $\gamma\to 0^{+}$, such that $\gamma_{\textrm{c}}=0$ and there is no interval of $\gamma$ for which the minimum wave speed is independent of $\gamma$. The implication of this finding is that interventions seeking to reduce the decay rate to $\gamma=0$ would eventually stop the invasion since $c_{\textrm{min}}\to 0^{+}$ as $\gamma\to 0^{+}$.

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,

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

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.

In this section we demonstrate a relationship between the three-dimensional dynamical system and phase space for the invasion model with the far simpler two-dimensional dynamical system and phase plane for the simpler Fisher-KPP model. . To motivate this 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.

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 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], but here we find that these 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.

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.

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

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

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

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

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.