Travelling-wave analysis of a model of tumour invasion with degenerate, cross-dependent diffusion

In this paper, we carry out a travelling-wave analysis of a model of tumour invasion with degenerate, cross-dependent diffusion. We consider two types of invasive fronts of tumour tissue into extracellular matrix (ECM), which represents healthy tissue. These types differ according to whether the density of ECM far ahead of the wave front is maximal or not. In the former case, we use a shooting argument to prove that there exists a unique travelling-wave solution for any positive propagation speed. In the latter case, we further develop this argument to prove that there exists a unique travelling-wave solution for any propagation speed greater than or equal to a strictly positive minimal wave speed. Using a combination of analytical and numerical results, we conjecture that the minimal wave speed depends monotonically on the degradation rate of ECM by tumour cells and the ECM density far ahead of the front.


Introduction
Tissue invasion is a hallmark of malignant tumours [1] and a classical mathematical approach to study 2021 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

(a) Reaction-diffusion PDE models of tumour invasion
To invade the surrounding healthy tissue, a tumour must overcome the defences developed by the body to maintain homeostatic control. An important barrier to tumour invasion is the ECM, which is a strong scaffold of proteins that holds tissue cells in place and initiates signalling pathways for cellular processes such as migration, differentiation and proliferation [6,7]. The healthy cells encased by the ECM form another barrier to invasion by creating a competitive environment for the tumour cells. However, tumour cells have developed mechanisms to overcome both of these barriers. First, they can remodel or degrade the ECM by producing specific matrix-degrading enzymes, which act in close proximity to the cells producing them [8,9]. Second, by favouring glycolytic metabolism even in aerobic conditions (i.e. the 'Warburg effect'), tumour cells may acidify the tissue microenvironment, resulting in healthy cell death [10,11]. Matrix remodelling is a very localized process, in contrast to the diffusion of lactic acid, which occurs on a longer spatial range.
The pioneering model by Gatenby & Gawlinksi [2] describes the spatio-temporal dynamics of acid-mediated tumour invasion by considering the interactions of healthy tissue, tumour tissue and the lactic acid produced by the tumour cells. Denoting the dimensionless tumour and healthy tissue densities and acid concentration by N(x, t), M(x, t) and L(x, t), respectively, for (x, t) ∈ R × (0, +∞), their model takes the form and Here, it is assumed that healthy cells do not move, while tumour cells can invade in a densitydependent manner. Depending on the value of α, the model describes the total or partial destruction of normal tissue following tumour invasion. We refer the reader to the original paper for full details of the model. A numerical study of the TWS of system (1.1), with 0 < D N 1, showed the existence of an interstitial gap, i.e. a region devoid of cells, formed locally ahead of the invading tumour front, for large values of α [2]. Experimental evidence has confirmed that such an interstitial gap can exist and, in this way, the model has led to novel and accurate predictions regarding tumour invasion. This is one of the reasons why this model and its variations have been widely investigated [3,[12][13][14][15][16]. nonlinear, degenerate diffusion, TWS have been extensively studied; see for instance [17][18][19][20][21][22][23][24][25]. In general, if the dimensionless equation has a reaction term, f , of Fisher-Kolmogoroff-Petrovsky-Piscounoff (Fisher-KPP) type, i.e. f ∈ C[0, 1] with f (0) = f (1) = 0 and f (n) > 0 ∀ n ∈ (0, 1), then TWS exist and are unique if and only if the wave speed is greater than or equal to a minimal speed, c * > 0, defined as the threshold speed below which no TWS exist. Further, if c = c * , then the TWS is of sharp type (that is, there is a discontinuity in the spatial derivative at the front) and, for each c > c * , there exists a TWS of front type (that is, smooth). It is non-trivial to extend such an existence result to R-D systems with multiple equations because of the added complexity of studying trajectories in a phase space rather than a phase plane. Kawasaki et al. [26] do so for an R-D system with cross-dependent diffusion developed to describe spatio-temporal pattern formation in colonies of bacteria. More specifically, numerical and analytical investigations [25,27] have shown the existence of TWS for wave speeds above or equal to a critical value, c * > 0. Until recently, most comprehensive results on the existence of TWS for spatially resolved models of tumour invasion focused on models in which invasion is driven by haptotaxis or chemotaxis [28][29][30][31]. In particular, the existence of TWS for the Gatenby-Gawlinski model has been largely supported by a combination of numerical and analytical results [12][13][14][15]32,33]. This also holds for a simplified model of invasion by Browning et al. [34], as seen in [35]. However, key existence results were recently proved by Gallay & Mascia [5] for a reduced version of the Gatenby-Gawlinski model.

(c) The mathematical model
We now present a minimal model of tumour invasion. There is increasing evidence that phenotypically heterogeneous tumours can contain sub-populations of cells with different traits, e.g. matrix-degrading cells and acid-producing cells [16]. Therefore, we make the simplifying assumption that the healthy tissue compartment solely comprises ECM, disregarding healthy cells, and we focus on the interactions of ECM-degrading tumour cells and ECM. Using a standard law for conservation of mass and denoting the tumour cell and ECM densities by N(x, t) and M(x, t), respectively, for (x, t) ∈ R × (0, +∞), we propose the following system of PDEs: We assume that the tumour grows logistically, with maximum growth rate ρ and carrying capacity K. Further, the ECM acts as a physical barrier that inhibits tumour cell movement, but not proliferation. Thus, following Gatenby & Gawlinski [2] and others [3,16,36], we define the diffusivity of tumour cells as a monotonically decreasing function of the ECM density to model the obstruction of movement by the ECM. The diffusivity of tumour cells in the absence of ECM is denoted by D N and the ECM density that inhibits all tumour cell movement is denoted by M Max . Finally, we assume that the ECM does not grow and is degraded at a rate that is proportional to the local tumour cell density, with a per cell degradation rate of k. We use a mass-action term to reflect the localized nature of matrix degradation.
To reduce the number of free parameters in the system and facilitate the analysis that follows, we non-dimensionalize equations (1.2) and, retaining the same dimensional state variables for notational convenience, we obtain the following system:  [15], studied by Gallay & Mascia [5], and a reduced model of melanoma invasion from Browning et al. [34], studied by El Hachem et al. [35]. In particular, while the models differ according to the reaction terms included for the healthy and tumour cell densities, each model retains the same degeneracy in the cross-diffusion term, which is the key focus of this paper.
Gallay & Mascia [5] rigorously proved the existence of a weak form of TWS for any positive wave speed, c > 0, for the model in [15]. These TWS represent the invasion of tumour tissue into healthy tissue, where the density of healthy tissue ahead of the wave front is at carrying capacity. El Hachem et al. [35] performed a numerical study which suggests that such TWS also exist for the model in [34] for any positive wave speed. In addition, their numerical results indicated the existence of another type of TWS for the model in [34] for wave speeds above a strictly positive minimal value. These TWS differ from the former in that the density of healthy tissue ahead of the wave front is below carrying capacity. Finally, El Hachem et al. [35] described the dependence of the minimal wave speed of this second type of TWS on the rescaled degradation rate of healthy tissue, which we denote by κ: it remains constant provided κ is below some threshold value, κ * , which is yet to be determined, and then increases with κ for κ ≥ κ * .
The key contribution of the present paper is to rigorously prove the existence of both aforementioned types of TWS for system (1.3), which we achieve by applying and expanding the shooting argument developed by Gallay & Mascia [5]. Similarly to [5,35], we will find that the first type of TWS exists for all c > 0, whereas there is a strictly positive minimal wave speed for the second type of TWS. We will see that this minimal wave speed for TWS of system (1.3) can be qualitatively characterized in the same way as that for equivalent TWS of the model in [34]. However, given κ, the value of this minimal wave speed for TWS of system (1.3) and of the system in [34] is not the same. A final contribution of our work compared with that in [35] is that we propose an expression for κ * for system (1.3), not reported in [35] for the model in [34].
(d) Structure of the paper We will seek constant profile, constant speed TWS for (1.3), which are heteroclinic trajectories of a three-dimensional dynamical system connecting two of its steady states. These correspond to spatially homogeneous, steady-state solutions of (1.3), which are given by Here, (N * 0 , M * 0 ) is the trivial state, (N * 1 , M * 1 ) is a state in which the tumour has successfully invaded and degraded all ECM and (N * 2 , M * 2 ) and (N * 3 , M * 3 ) are, respectively, a tumour-free state at maximum ECM density and a continuum of tumour-free states. We distinguish (N * 2 , M * 2 ) from (N * 3 , M * 3 ) because of the degeneracy at M = 1 in system (1.3). Since we are interested in studying the existence of TWS that describe the invasion of a tumour into healthy tissue, we will look for two types of heteroclinic trajectories: those connecting (N * 1 , M * 1 ) to (N * 2 , M * 2 ) and those connecting (N * 1 , M * 1 ) to (N * 3 , M * 3 ). In §2, we define the TWS we seek, prove preliminary results and derive the ordinary differential equation (ODE) system they must satisfy. In §3, we use the shooting argument developed by Gallay & Mascia [5] to show that system (1.3) has a unique TWS connecting (N * 1 , M * 1 ) to (N * 2 , M * 2 ) for any positive wave speed. We then show that, for eachM ∈ [0, 1), system (1.3) has a unique TWS connecting (N * 1 , M * 1 ) to (N * 3 , M * 3 ) for any wave speed greater than or equal to a strictly positive minimum value. Motivated by our numerical simulations and partial analytical results, we make a conjecture about the dependence of the minimal wave speed onM ∈ [0, 1) and κ > 0, the rescaled degradation rate of the ECM. In §4, we present numerical simulations of system (1.3) which support and complement the preceding analytical results. We conclude the paper in §5, where we discuss our results alongside future research perspectives.

The travelling-wave problem (a) Preliminaries
We seek constant profile, constant speed TWS of system (1.3) by introducing the travelling-wave coordinate ξ = x − ct. We require the wave speed c > 0 so that the tumour invades the ECM from left to right in the spatial domain. Substituting the ansatz N(x, t) = N (ξ ) and M(x, t) = M(ξ ) into system (1.3), we deduce that TWS must satisfy the following ODE system: The TWS we seek connect spatially homogeneous steady states of system (1.3) and, equivalently, steady states of system (2.1a,b). Thus, we require one of the following sets of asymptotic conditions to be satisfied: In other words, far behind the wave, the tumour density is at carrying capacity and the ECM has been fully degraded, whereas, far ahead of the wave, the tumour density is zero and the ECM density is either at carrying capacity (i.e. M = 1) or at any value M ∈ [0, 1). As noted previously, the first Hence, unless otherwise stated, we refer to weak TWS in the sense of definition 2.1 as TWS. If (N , M; c) is a TWS for system (1.3), then we can show that N (1 − N ) ∈ L 1 (R) and c > 0 using a proof identical to that of lemma 2.1 in [5] and, thus, we omit it. The following lemma, whose proof follows as in [5], states that if (N , M; c) is a TWS for system (1.3), then N and M are non-negative and bounded and, thus, the TWS is biologically realistic.   In this case, system (2.1a,b) reduces to the Fisher-KPP equation, which has been extensively studied [37][38][39]. It is known that the Fisher-KPP equation admits classical TWS that satisfy the asymptotic conditions lim ξ →−∞ N (ξ ) = 1, lim y→+∞ N (ξ ) = 0 and lim ξ →±∞ (dN /dξ )(ξ ) = 0 for all c ≥ 2. This result, therefore, holds for TWS of (2.1a,b) satisfying the asymptotic conditions (2.3) forM = 0.
A version of lemma 2.2 for TWS that satisfy the asymptotic conditions (2.2) follows similarly [5]. These results highlight that the solutions we seek are classical solutions of system (2.1a,b) on intervals of the form (−∞,ξ ). Further, ifξ = +∞, then the TWS are here called smooth. In contrast, ifξ < +∞, then lemma 2.2 implies that we have a corner point atξ and the TWS are here called sharp.
(b) Desingularization of the ODE system Definition 2.1 describes two types of TWS of system (2.1a,b), which differ in the asymptotic conditions they satisfy at infinity. One type of solution converges to (N , M) = (0, 1) at infinity. Therefore, we need to elucidate the behaviour of solutions as they approach M = 1, which is precisely when system (2.1a,b) is singular. A common approach to simplify the analysis is to remove this singularity by re-parametrizing the system. Given a solution (N , M) of system (2.1a,b) satisfying either (2.2) or (2.3), we introduce a new independent variable y = Φ(ξ ) defined such that Further introducing the following dependent variables: we can apply the chain rule and use (2.6) to find that, for 0 ≤ m ≤ 1, the trajectories satisfy the following ODE system, for y ∈ R: In line with the asymptotic conditions (2.2) and (2.3), we require one of the following to hold: (2.10) Importantly, system (2.8a,b) is topologically equivalent to system (2.1a,b) for (N , M) ∈ (0, 1) 2 . This follows from the fact that (2.7) defines a homeomorphism that maps the orbits of (2.1a,b) onto the orbits of (2.8a,b), while preserving their orientation-(2.6) implies that y is an increasing function of ξ for all 0 ≤ M < 1. We also observe that, in contrast to system (2.1a,b), system (2.8a,b) has an additional continuum of steady states of the form (n, m) = (n, 1),n ∈ (0, 1]. These are not spatially homogeneous steady states of the original PDE system (1.3), so we do not consider them as asymptotic conditions in the context of TWS.

Travelling-wave analysis
In this section, we study the existence of TWS. To do so, we apply the shooting argument developed by Gallay & Mascia [5]. The crucial difference between Gallay & Mascia's model and system (1.3) is that the latter has an additional continuum of steady states of the form (0,M), M ∈ (0, 1). We find that the results of [5] for TWS connecting the equilibrium points (1, 0, 0) and (0, 0, 1) apply, with minor modifications, to the TWS of system (2.11a-c) that satisfy the same asymptotic conditions (2.2). Therefore, in what follows, we state the key results and present only those proofs which require a different approach. For TWS of system (2.11a-c) that satisfy the asymptotic conditions (2.3), we further develop the shooting argument to obtain new results.
(a) Local analysis of the equilibrium point (1, 0, 0): defining the shooting parameter The TWS of interest satisfy lim y→−∞ (n(y), p(y), m(y)) = (1, 0, 0). We therefore study the behaviour of solutions of (2.11a-c) in a neighbourhood of the equilibrium point P 1 := (1, 0, 0) by performing a linear stability analysis. The Jacobian matrix at P 1 reduces to and it has the following eigenvalues and eigenvectors: and Since λ 1 is negative and λ 2 and λ 3 are positive, P 1 is a three-dimensional hyperbolic saddle point with a two-dimensional unstable manifold, which locally is a plane through P 1 generated by the eigenvectors v 2 and v 3 . There is also a one-dimensional stable manifold which locally is a straight line spanned by the eigenvector v 1 . Trajectories defined by (2.11a-c) that leave P 1 must do so via the two-dimensional unstable manifold at P 1 . We therefore compute asymptotic expansions of all solutions of (2.11a-c) in a neighbourhood of P 1 that lie on the unstable manifold. Requiring that n ∈ (0, 1) and p < 0, so that solutions leaving P 1 remain in D 1 , we obtain the following result.
Remark 3.2. The free parameter, α, arises because the form taken by the unstable manifold at P 1 does not impose any condition on m. In a sense, the choice of α is a choice of how fast m increases from 0 and, accordingly, α will influence the value that m attains at y = +∞. We illustrate this in figure 1 and present some corresponding travelling-wave profiles in electronic supplementary material, S2. In addition, by remark 2.3, it is clear that α = 0 is the unique value of the shooting parameter such that the solution of (2.11a-c) that satisfies (3.3) stays in a region where n ∈ (0, 1), p < 0 and m = 0 and satisfies the asymptotic conditions (2.10) form = 0. Now, the idea is to view solutions of (2.11a-c) that satisfy (3.3) as functions of α, which we define to be our shooting parameter, and c, which is the wave speed. In particular, we denote by (n α,c , p α,c , m α,c ) the unique solution of (2.11a-c) satisfying (3.3). Our first result, which can be proved following the approach in [5], is the following: is defined on some interval J := (−∞, y 0 ), with y 0 ∈ R, and satisfies n α,c (y) > 0 for all y ∈ J, then (n α,c (y), p α,c (y), m α,c (y)) ∈ D 1 for all y ∈ J.
Given lemma 3.3, we introduce the following variable for any α > 0 and c > 0: Then, lemma 3.3 implies that only one of the following holds: -T(α, c) < +∞, so n α,c (T(α, c)) = 0 and p α,c (T(α, c)) < 0. In this case, n α,c (y) becomes negative for some y > T(α, c) and (n α,c , p α,c , m α,c ) does not represent a valid TWS; we disregard these values of the shooting parameter α. -T(α, c) = +∞, which means that we have a global solution which stays in D 1 for all y ∈ R. We are interested in finding TWS for these values of α.
Remark 3.4. Givenm ∈ (0, 1), lemma 3.3 provides a condition under which solutions of (2.11a-c) that satisfy (3.3) remain in D 1 , but not necessarily in Dm ⊂ D 1 . In particular, even if n α,c (y) > 0 for all y ∈ J, a solution can leave Dm. In that case, for a solution of (2.11a-c) that satisfies (3.3) to converge to (0, 0,m) as y → +∞, we must have n(y) < 0 for some values of y (since m is increasing for positive n). Therefore, searching for solutions that satisfy T(α, c) = +∞ is a necessary condition for the existence of physically realistic TWS that converge to (0, 0,m) as y → +∞, but not a sufficient one (m may not attain the valuem for positive n).

(b) Monotonicity of solutions with respect to the shooting parameter
A key component of our analysis is that solutions of system (2.11a-c) satisfying (3.3) are monotonic functions of the shooting parameter, α, provided n > 0. This result, which can be proved following the approach in [5], can be formulated as follows: for all y ∈ (−∞, T(α 1 , c)).
Lemma 3.5 shows that, for fixed c > 0, T(α, c) is an increasing function of α. Since we seek TWS that satisfy T(α, c) = +∞, we define the following critical value of α, which depends on c: We then characterize α 0 as a function of c in the following lemma, whose proof follows similarly to that of lemma 2.7 in [5]: Lemma 3.6 ensures that for all c > 0 there exist some values of α for which T(α, c) = +∞, and thus for all c > 0 there exist some TWS. We still need to elucidate the behaviour of solutions at infinity to determine which TWS exist for each c > 0.
Remark 3.7. The proof of lemma 3.6 relies on showing that, for any 0 < c < 2, we can choose α(c) > 0 sufficiently large such that there exists a solution of system (2.11a-c) that satisfies (3.3), remains in region D 1 and converges to (n, 0, 1) as y → +∞, withn ∈ (0, 1). Such solutions are not TWS, but their existence will be crucial in proving the existence of the TWS we seek.

(c) Behaviour of solutions at infinity
By lemma 3.6, we know that, for any c > 0, there exist solutions of system (2.11a-c) that satisfy (3.3) and remain in region D 1 for all y ∈ R. It remains to characterize the behaviour of these solutions as y → +∞, and, in so doing, to establish whether they are TWS. Denoting the limits of components of the solution at infinity as we introduce the following lemma, which can be proved following the approach in [5].  Lemma 3.8 defines the possible limits of solutions (n α,c , p α,c , m α,c ) of (2.11a-c) that satisfy (3.3) and remain in region D 1 . We must now determine for which values of c > 0 we can find α(c) > 0 such that the corresponding solution (n α,c , p α,c , m α,c ): (i) remains in D 1 and converges to (0, 0, 1) as y → +∞, or, (ii) for eachm ∈ (0, 1), remains in Dm and converges to (0, 0,m) as y → +∞.
We consider these cases separately in the two sections that follow.
(d) Solutions converging to the equilibrium point (0, 0, 1) In this section, we show that, for each c > 0, there exists a unique value of α > 0 such that the solution (n α,c , p α,c , m α,c ) of system (2.11a-c) satisfying (3.3) remains in region D 1 and converges to the equilibrium point P 2 := (0, 0, 1) as y → +∞. This then allows us to draw conclusions on the existence and uniqueness of TWS that satisfy the asymptotic conditions (2.2).
For the rest of this section, we suppose that α ≥ α 1 (c). A linear stability analysis at the equilibrium point P 2 shows that P 2 is non-hyperbolic, with one negative eigenvalue and two zero eigenvalues. Therefore, at P 2 , we have a one-dimensional stable manifold, W S ⊂ R 3 , generated by the eigenvector v 1 = (1/c, 1, 0) associated with λ 1 = −c. We also have a two-dimensional centre manifold, W C ⊂ R 3 , which is tangent at P 2 to the subspace spanned by the eigenvectors v 2 = (1, 0, 0) and v 3 = (0, 0, 1) associated with λ 2 = λ 3 = 0. Solutions of (2.11a-c) that satisfy (3.3) and remain in a small enough neighbourhood of P 2 for all sufficiently large y > 0 converge to W C . Therefore, in order to study the dynamics around P 2 , we perform a nonlinear local stability analysis. We begin by transforming system (2.11a-c) into normal form by introducing the following variables:ñ which satisfy the following system: Then, we know that, in a neighbourhood of the origin, the centre manifold can be described by a function P(ñ,m) such that (ñ,p,m) ∈ W C if and only ifp = P(ñ,m), where Using this expression for the centre manifold in a neighbourhood of the origin, we must now prove that there is a solution of system (3.11a-c) converging to the centre manifold W C that converges to the origin as y → +∞. We are interested in solutions (n, p, m) of (2.11a-c) that satisfy (3.3) and remain in region D 1 for all y ∈ R. Equivalently, we seek solutions (ñ,p,m) of (3.11a-c) that satisfy (3.3) and lie on a manifold W + C ⊂ W C , where The following lemma characterizes such solutions that converge to the origin as y → +∞ (the proof corresponds, with minor modifications, to that of lemma 2.12 in [5]). Lemma 3.10. Up to translations in the variable y, there exists a unique solution of (3.11a-c) that satisfies the asymptotic conditions (3.3), lies on the centre manifold W + C and whose components converge to zero as y → +∞, such that (3.14) Lemma 3.10 establishes the existence of at least one solution of (2.11a-c) that satisfies (3.3), stays in region D 1 and converges to (1, 0, 0) as y → +∞. Furthermore, this solution is uniquely determined on the centre manifold W + C . Given that any solution of (2.11a-c) that satisfies (3.3), stays in region D 1 and converges to (0, 0, 1) as y → +∞ must do so via W + C and, given the monotonicity result of lemma 3.5, it is easy to prove the following lemma as in [5].
Using lemma 3.12 and reversing the change of variables (2.6), it is straightforward to construct a unique (up to translation) solution (N (ξ ), M(ξ )) = (n(Φ(ξ )), m(Φ(ξ ))) of system (2.1a,b) that satisfies the asymptotic conditions (2.2). This leads to our first main result, which can be proved following the approach in [5]: In this section, we consider solutions of system (2.11a-c) subject to (3.3) that stay in region Dm ⊂ D 1 and converge to (0, 0,m) form ∈ (0, 1) as y → +∞. Using arguments similar to those for the previous case, we can show that, for allm ∈ [0, 1), there exists a strictly positive, real-valued wave speed above which the solutions we seek exist and are unique. We will refer to this wave speed as the minimal wave speed and we will observe that it depends on κ, the rescaled degradation rate of ECM. In particular, given κ > 0, we denote the minimal wave speed by c * κ (m) for eachm ∈ [0, 1). This will enable us to draw some conclusions on the existence and uniqueness of TWS that satisfy the asymptotic conditions (2.3).
The preceding hypothesis and conjecture 3.17 are supported by numerical simulations of the PDE system (1.3) and ODE system (2.11a-c). In figure 2, we show that solutions of system (1.3) subject to the initial conditions (4.1) withM ∈ [0, 1) evolve into travelling waves with constant propagation speed (see electronic supplementary material, S2 for corresponding travelling-wave profiles). We observe that, for 0 < κ ≤ κ * (M), this speed is independent of κ, and, calculating the slopes of these lines, we find that it is approximately equal to 2 1 −M. Additionally, when κ > κ * (M), we observe that the wave speed selected by the PDE increases with κ. We also solved numerically the system (2.11a-c), subject to the asymptotic conditions (3.3), for the same values of κ > κ * (M) and the respective values of the propagation speed estimated using the solutions of the PDE system (results not shown). We observed that, given κ > κ * (M), the wave speed selected by the PDE appears to correspond to the smallest wave speed such that the solution (n, p, m) of the system (2.11a-c), subject to (3.3), satisfies n(y) > 0 ∀ y ∈ R and converges to (0, 0,m),m =M. Now, suppose conjecture 3.17 is true. Then, given κ > 0, for eachm ∈ [0, m * (κ)] and c ≥ 2 √ 1 −m, αm(c) as defined by (3.19) exists and is the unique α mentioned in the statement of conjecture 3.17. Using remark 3.16 and conjecture 3.17, the subsequent result follows naturally (we omit the proof for brevity). Lemma 3.18. Suppose conjecture 3.17 is true and fix κ > 0. If c ≥ 2 1 − m * (κ), then, for all m ∈ (m * (κ), 1), there exists a unique α ∈ (α m * (κ) (c), α 1 (c)) such that the solution of (2.11a-c) that satisfies the asymptotic properties in lemma 3.1 converges to (0, 0,m) as y → +∞.
While we do not have a complete characterization of the minimal wave speed for all κ > 0 and m ∈ (0, 1), we can now state our second main result. It can be proved similarly to theorem 3.13, so we refer the reader to [5] for its proof.  has already spread to a position x = σ < L in the tissue and we impose initial conditions that satisfy, forM (4.1) Here, 0 < ω < σ represents how sharp the initial boundary between the tumour and healthy tissue is. We complete the mathematical problem by imposing zero-flux boundary conditions for N at x = 0 and x = L. We set L = 200, σ = 2 and ω = 1 for our simulations.
Remark 4.1. Initial conditions for N with compact support, such as those given by (4.1), are biologically relevant. We verified that the travelling-wave profile and wave speed are preserved across different initial conditions with compact support for N, i.e. initial conditions of the type of (4.1) (see electronic supplementary material, S2).
(a) Elucidating the wave speed that emerges in the PDE model A characteristic feature of the well-studied Fisher-KPP model is that any non-negative initial condition with compact support will evolve towards a travelling front with speed equal to the minimal wave speed, c = 2 [37][38][39]. One may, therefore, question whether this result extends to more complex R-D systems that exhibit travelling waves. For our model, the results from §3e suggest that this does hold for solutions of (1.3) subject to the initial conditions (4.1) for M ∈ [0, 1). By contrast, the results from §3d show that there is no strictly positive minimal wave speed for TWS of (1.3) that satisfy the asymptotic conditions (2.2). Yet, the solution of (1.3) subject to the initial conditions (4.1) forM = 1 appears to evolve towards a travelling front with a strictly positive speed, as illustrated in figure 3a for different values of κ (see electronic supplementary material, S2 for a travelling-wave profile). In this way, the solutions of the PDE system preferentially select a wave speed in a way that the corresponding ODE system does not.
Given different values of κ > 0, we calculated the speed of travelling fronts that emerge for solutions of (1.   suggest that the wave speed selected by the PDE model is a continuous, decreasing function of M ∈ [0, 1], as illustrated in figure 3b, which represents this wave speed as a function ofM ∈ [0, 1] for κ ∈ {0, 1, 10, 100, 1000, 10 000}. This is consistent with lemma 3.19 and our observation that the speed of travelling fronts that emerge for solutions of (1.3), subject to the initial conditions (4.1) withM ∈ [0, 1), appears to be equal to the minimal wave speed, c * κ (M), defined by (3.25). This result is interesting because the speed selected by the PDE model appears to be leftcontinuous atM = 1, despite the fact that the minimal wave speed for the existence of TWS is not.

Remark 4.2.
Here, we investigated numerically the value of the propagation speed selected by the PDE model. Since the spatial domain, X , must be discretized to solve (1.3) using the method of lines, a natural question is whether the size of the discretization step influences the value of the wave speed. Given TWS of (1.3) that connect (1, 0) and (0,M),M ∈ [0, 1], we observed that the impact of decreasing the discretization step size becomes significant asM approaches 1 (see electronic supplementary material, S2). On the basis of the results illustrated in electronic supplementary material, S2, the discretization step size we used for our numerical simulations ensures that the numerical results obtained are weakly affected by numerical diffusion. In particular, our qualitative descriptions of the wave speed selected by the PDE model are unaffected.

Discussion and perspectives
Understanding the process of tumour invasion is at the forefront of cancer research. The seminal model of acid-mediated tumour invasion developed by Gatenby & Gawlinski [2] generated new biological insights and formed the basis for subsequent mathematical work on this topic. Owing to the model's complexity, most existing results in the literature on the existence of TWS of the model stem from numerical investigations, which are complemented by partial analytical results. In particular, obtaining a complete understanding of the existence of TWS has proven difficult and this has prompted the derivation of simpler models [15,34]. In this paper, we carried out a travelling-wave analysis for the simplified model (1.3). We found that system (1.3) can support a continuum of smooth TWS, which are here defined as TWS for whichξ , introduced in lemma 2.2, satisfiesξ = +∞. These TWS represent the invasion of healthy tissue, consisting of ECM, by tumour cells and differ according to the density of ECM far ahead of the wave front. More specifically, we characterized TWS connecting the two spatially homogeneous steady states (1, 0) and (0,M), forM ∈ [0, 1]. Owing to the degeneracy in the first equation of (1.3) for M = 1, we distinguished the cases whereM = 1 and whereM ∈ [0, 1).
In the first case, we proved the existence of smooth TWS for any positive wave speed, c > 0. This result is particularly interesting as it differs from previous results for degenerate diffusion in a scalar or multi-equation setting, where TWS exist if and only if the wave speed is greater than or equal to a strictly positive minimal wave speed [22,24,25]. It is important to note that this does not imply that a positive wave speed which is preferentially selected does not exist for solutions of (1.3) that connect (1, 0) and (0, 1). In fact, we saw in §4a that a strictly positive, κ-dependent wave speed appears to be selected by (1.3) subject to the initial conditions (4.1) withM = 1. It would, therefore, be interesting to study the stability of the TWS defined by theorem 3.13. We may gain insight on the minimal wave speed for solutions of (1.3) that connect (1, 0) and (0, 1) by determining parameter regimes in which solutions are unstable.
In the second case, we proved that smooth TWS exist if and only if the wave speed is greater than or equal to a strictly positive minimal wave speed, c * κ (M), defined by (3.25) for M ∈ [0, 1). Given κ > 0, this minimal speed appears to be a monotonically decreasing, continuous function ofM. In particular, we conjectured that, given κ > 0 and m * (κ) := 1/(κ + 1), we can precisely define c * κ (M) = 2 1 −M forM ∈ [0, m * (κ)]. Similarly to the generalized Fisher-KPP equation, this minimal wave speed is the smallest c > 0 such that the equilibrium (0, 0,m), with m =M, of system (2.11a-c) is a stable node and not a stable spiral. ForM ∈ (m * (κ), 1), numerical simulations suggested that the wave speed selected by the PDE is strictly greater than 2 1 −M, which is consistent with (3.25). The fact that the equilibrium (0, 0,m) of system (2.11a-c) is a stable node is then no longer a sufficient condition to ensure the positivity of the n-component of the TWS in the desingularized variables and thus of the N -component of the TWS in the original variables. This reflects the differences that can be observed in systems of equations compared with scalar equations, which can be attributed to the higher dimensionality of the problem.
Our results regarding the dependence of the minimal wave speed on the model parameters κ andM for TWS of (1.3) connecting (1, 0) and (0,M),M ∈ [0, 1) rely on a conjecture. Our aim is to rigorously prove this result in future work. In addition, we do not have an expression for the minimal wave speed ifM ∈ (m * (κ), 1). Yet, as κ → +∞, m * (κ) → 0, and it is clear that, as κ increases, we can precisely describe the minimal wave speed for a decreasing range of values of M ∈ [0, 1). We would therefore like to provide a complete characterization of c * κ (M) for all κ > 0 andM ∈ (m * (κ), 1). Now, we observed in §4a that the solution of system (1.3) subject to initial conditions (4.1) withM ∈ [0, 1] evolves towards a travelling front with a κ-andM-dependent wave speed. Importantly, given κ > 0, it appears that this numerical wave speed is a continuous function ofM in [0, 1], is equal to c * κ (M) = 2 1 −M for allM ∈ [0, m * (κ)] and is strictly greater than 2 1 −M for allM ∈ (m * (κ), 1]. We note that we have includedM = 1 in our preceding observations, which highlights our hypothesis that elucidating the minimal wave speed for (1. M ∈ (m * (κ), 1), or vice versa. It is, therefore, important to also study the stability of the travelling waves defined by theorem 3.20.
Finally, while it is of mathematical interest to obtain a comprehensive description of the minimal wave speed for all TWS of (1.3), it is also of biological interest. Our results indicate that the minimal wave speed is highly dependent on the value of κ, which is the rescaled ECM degradation rate. Since this parameter represents, in a sense, the aggressivity of the tumour cell population towards the ECM, it is significant from an oncological perspective. Hence, our results have the long-term potential of revealing promising targets for therapeutic intervention.
Data accessibility. Matlab codes to numerically solve the models are available in the electronic supplementary material.
Authors' contributions. F.S.-G. and P.K.M. conceived the study, all authors designed the study, C.C. and T.L.
executed the study and all authors contributed to the writing of the manuscript. P.K.M. and T.L. are joint senior authors.
Competing interests. We declare we have no competing interests. Funding. C.C. is supported by funding from the Engineering and Physical Sciences Research Council (EPSRC).