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 this process involves reaction–diffusion (R–D) partial differential equations (PDEs) [2-4]. A key feature of many such models of tumour invasion is the inclusion of degenerate, cross-dependent diffusion. The aim of this paper is to study this common characteristic by proposing a minimal model which captures the main components of the tumour invasion process and is analytically tractable. We seek two types of constant profile, constant speed travelling-wave solutions (TWS) for our model. Both types represent invasive fronts of tumour tissue into extracellular matrix (ECM), which represents healthy tissue, but they differ according to whether the density of ECM far ahead of the wave front is maximal or not. For the former, we prove the existence and uniqueness of TWS for all positive propagation speeds using the shooting argument developed by Gallay & Mascia [5]. For the latter, we expand this shooting argument to prove the existence and uniqueness of TWS for propagation speeds greater than or equal to a strictly positive minimal value. Finally, we characterize this minimal wave speed using a conjecture motivated by a combination of analytical results and numerical simulations.

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 , and , respectively, for , their model takes the form Here, it is assumed that healthy cells do not move, while tumour cells can invade in a density-dependent 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 , 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-16].

Nonlinear, degenerate diffusion: from scalar to multi-dimensional analysis

A key common component of the Gatenby–Gawlinksi model and its variations is the degenerate, cross-diffusion term in the equation for the tumour cell density. For scalar R–D equations with nonlinear, degenerate diffusion, TWS have been extensively studied; see for instance [17-25]. In general, if the dimensionless equation has a reaction term, , of Fisher–Kolmogoroff–Petrovsky–Piscounoff (Fisher–KPP) type, i.e. with and , then TWS exist and are unique if and only if the wave speed is greater than or equal to a minimal speed, , defined as the threshold speed below which no TWS exist. Further, if , then the TWS is of sharp type (that is, there is a discontinuity in the spatial derivative at the front) and, for each , 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, . 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-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-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.

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 and , respectively, for , we propose the following system of PDEs:

We assume that the tumour grows logistically, with maximum growth rate and carrying capacity . 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 and the ECM density that inhibits all tumour cell movement is denoted by . 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 . 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: where .

Here, we note that system (1.3) is similar to a reduced version of the model (1.1) from Moschetta & Simeoni [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, , 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 , 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].

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, is the trivial state, is a state in which the tumour has successfully invaded and degraded all ECM and and are, respectively, a tumour-free state at maximum ECM density and a continuum of tumour-free states. We distinguish from because of the degeneracy at 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 to and those connecting to . 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 to for any positive wave speed. We then show that, for each , system (1.3) has a unique TWS connecting to 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 on and , 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

Preliminaries

We seek constant profile, constant speed TWS of system (1.3) by introducing the travelling-wave coordinate . We require the wave speed so that the tumour invades the ECM from left to right in the spatial domain. Substituting the ansatz and into system (1.3), we deduce that TWS must satisfy the following ODE system: and

The TWS we seek connect spatially homogeneous steady states of system (1.3) and, equivalently, steady states of system (2.1a,). Thus, we require one of the following sets of asymptotic conditions to be satisfied: or 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. ) or at any value . As noted previously, the first equation in system (1.3) is a degenerate parabolic equation since the cross-diffusion coefficient is zero when . The existence of global classical solutions of this PDE system and the corresponding ODE system (2.1a,) is therefore unclear in cases where or, correspondingly, where . We therefore define a weak TWS in a similar way to the definition of a propagation front in [5].

Definition 2.1.

The triple is called a weak TWS for system (1.3) if

and ;

is a weak solution of (2.1a,), i.e. for all with compact support and

one of the pairs of asymptotic conditions given by (2.2) and (2.3), respectively, are satisfied.

We refer to as the travelling-wave profile and as the propagation speed.

Note.

Hence, unless otherwise stated, we refer to weak TWS in the sense of definition 2.1 as TWS.

If is a TWS for system (1.3), then we can show that and 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 is a TWS for system (1.3), then and are non-negative and bounded and, thus, the TWS is biologically realistic.

Lemma 2.2.

If is a weak TWS, in the sense of definition 2.1, that satisfies the asymptotic conditions (2.3) for , then there exists a unique point such that

and for ;

if , then and for all

Remark 2.3.

The case of TWS that satisfy the asymptotic conditions (2.3) for is not considered in lemma 2.2. By definition, such solutions satisfy for , which is only possible if on since is increasing for In this case, system (2.1a,) reduces to the Fisher–KPP equation, which has been extensively studied [37-39]. It is known that the Fisher–KPP equation admits classical TWS that satisfy the asymptotic conditions , and for all . This result, therefore, holds for TWS of (2.1a,) satisfying the asymptotic conditions (2.3) for .

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,) 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.

Desingularization of the ODE system

Definition 2.1 describes two types of TWS of system (2.1a,), which differ in the asymptotic conditions they satisfy at infinity. One type of solution converges to at infinity. Therefore, we need to elucidate the behaviour of solutions as they approach , which is precisely when system (2.1a,) is singular. A common approach to simplify the analysis is to remove this singularity by re-parametrizing the system. Given a solution of system (2.1a,) satisfying either (2.2) or (2.3), we introduce a new independent variable defined such that Further introducing the following dependent variables: we can apply the chain rule and use (2.6) to find that, for , the trajectories satisfy the following ODE system, for : and

In line with the asymptotic conditions (2.2) and (2.3), we require one of the following to hold: or

Importantly, system (2.8a,) is topologically equivalent to system (2.1a,) for . This follows from the fact that (2.7) defines a homeomorphism that maps the orbits of (2.1a,) onto the orbits of (2.8a,), while preserving their orientation—(2.6) implies that is an increasing function of for all . We also observe that, in contrast to system (2.1a,), system (2.8a,) has an additional continuum of steady states of the form , . 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.

We finally obtain a system of three first-order ODEs by introducing the additional variable and, using primes to denote derivatives with respect to , we have

In the following section, we set up a framework, first proposed in [5] for a different system, to study two distinct types of solutions of (2.11a–c). First, those that remain in the region , defined as and that satisfy , . Second, for , those that remain in the region , defined similarly to (2.12) as and that satisfy , .

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 , . We find that the results of [5] for TWS connecting the equilibrium points and 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.

Local analysis of the equilibrium point : defining the shooting parameter

The TWS of interest satisfy . We therefore study the behaviour of solutions of (2.11a–c) in a neighbourhood of the equilibrium point by performing a linear stability analysis. The Jacobian matrix at reduces to and it has the following eigenvalues and eigenvectors: and Since is negative and and are positive, is a three-dimensional hyperbolic saddle point with a two-dimensional unstable manifold, which locally is a plane through generated by the eigenvectors and . There is also a one-dimensional stable manifold which locally is a straight line spanned by the eigenvector . Trajectories defined by (2.11a–c) that leave must do so via the two-dimensional unstable manifold at . We therefore compute asymptotic expansions of all solutions of (2.11a–c) in a neighbourhood of that lie on the unstable manifold. Requiring that and , so that solutions leaving remain in , we obtain the following result.

Lemma 3.1.

Fix For any , the system (2.11a–c) has a unique solution such that, as , where and are given by (3.2) and

Remark 3.2.

The free parameter, , arises because the form taken by the unstable manifold at does not impose any condition on . In a sense, the choice of is a choice of how fast increases from and, accordingly, will influence the value that attains at . 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 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 , and and satisfies the asymptotic conditions (2.10) for .

Figure 1

Solutions of (2.11a–c) subject to the asymptotic conditions (3.3) for different values of the shooting parameter , and (a) or (b). The purple lines are the two continua of steady states of the system (2.11a–c), given by , , and , , respectively. Since , , are not spatially homogeneous steady states of (1.3), the dashed curves represent solutions that are not TWS of system (1.3). The dotted curves represent physically unrealistic solutions for which the -component becomes negative. The values and attain at infinity appear to increase monotonically (between 0 and 1) with . (Online version in colour.)

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 , which is the wave speed. In particular, we denote by 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:

Lemma 3.3.

If the solution is defined on some interval , with , and satisfies for all , then for all

Given lemma 3.3, we introduce the following variable for any and : Then, lemma 3.3 implies that only one of the following holds:

, so and . In this case, becomes negative for some and does not represent a valid TWS; we disregard these values of the shooting parameter .

, which means that we have a global solution which stays in for all . We are interested in finding TWS for these values of .

Remark 3.4.

Given , lemma 3.3 provides a condition under which solutions of (2.11a–c) that satisfy (3.3) remain in , but not necessarily in . In particular, even if for all , a solution can leave . In that case, for a solution of (2.11a–c) that satisfies (3.3) to converge to as , we must have for some values of (since is increasing for positive ). Therefore, searching for solutions that satisfy is a necessary condition for the existence of physically realistic TWS that converge to as , but not a sufficient one ( may not attain the value for positive ).

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 . This result, which can be proved following the approach in [5], can be formulated as follows:

Lemma 3.5.

Fix . If , then and the solutions of (2.11a–c) defined by (3.3) satisfy for all .

Lemma 3.5 shows that, for fixed , is an increasing function of . Since we seek TWS that satisfy , we define the following critical value of , which depends on : We then characterize as a function of in the following lemma, whose proof follows similarly to that of lemma 2.7 in [5]:

Lemma 3.6.

If , then . If , then .

Lemma 3.6 ensures that for all there exist some values of for which , and thus for all there exist some TWS. We still need to elucidate the behaviour of solutions at infinity to determine which TWS exist for each .

Remark 3.7.

The proof of lemma 3.6 relies on showing that, for any , we can choose sufficiently large such that there exists a solution of system (2.11a–c) that satisfies (3.3), remains in region and converges to as , with . Such solutions are not TWS, but their existence will be crucial in proving the existence of the TWS we seek.

Behaviour of solutions at infinity

By lemma 3.6, we know that, for any , there exist solutions of system (2.11a–c) that satisfy (3.3) and remain in region for all . It remains to characterize the behaviour of these solutions as , 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.

If , then the following limits exist: Moreover, if , then , and, if , then .

Lemma 3.8 defines the possible limits of solutions of (2.11a–c) that satisfy (3.3) and remain in region . We must now determine for which values of we can find such that the corresponding solution :

remains in and converges to as , or,

for each , remains in and converges to as .

We consider these cases separately in the two sections that follow.

Solutions converging to the equilibrium point

In this section, we show that, for each , there exists a unique value of such that the solution of system (2.11a–c) satisfying (3.3) remains in region and converges to the equilibrium point as . This then allows us to draw conclusions on the existence and uniqueness of TWS that satisfy the asymptotic conditions (2.2).

By remark 3.7, we have that, for any , there exists sufficiently large such that the solution of system (2.11a–c) satisfying (3.3) remains in region and satisfies . We can therefore define and prove the following result:

Lemma 3.9.

For any , we have .

Proof.

Fix . By remark 3.7, we know that there exists large enough such that and, hence, . If , then we know by lemma 3.6 that and, therefore, by the definition of , we must have .

If , then suppose, for a contradiction, that . A linear stability analysis about the equilibrium point with shows that it is non-hyperbolic with two negative eigenvalues , and one zero eigenvalue , Therefore, by the Centre Manifold Theorem, in a small, open neighbourhood of with , there exists a two-dimensional stable manifold spanned by the eigenvectors corresponding to . In this neighbourhood, there also exists a one-dimensional centre manifold spanned by , which comprises the family of equilibria with sufficiently close to . Therefore, for fixed , we can find a neighbourhood of that is foliated by two-dimensional stable leaves over a one-dimensional centre manifold, composed of points of the form with . Then, any solution that enters converges to as for some that satisfies .

Now, by remark 3.2, implies that . Thus, we can find large enough such that for all . By continuity of solutions with respect to , we can find such that for any . This implies that, for any such , converges to as . Since by our choice of , there exists such that . However, since , we must have for all and we have reached the desired contradiction.

Lemma 3.9 ensures that, for any and , the solution of (2.11a–c), subject to the asymptotic conditions (3.3), stays in region and satisfies , where . We would now like to show that, for any , there exists a unique such that .

For the rest of this section, we suppose that . A linear stability analysis at the equilibrium point shows that is non-hyperbolic, with one negative eigenvalue and two zero eigenvalues. Therefore, at , we have a one-dimensional stable manifold, , generated by the eigenvector associated with . We also have a two-dimensional centre manifold, , which is tangent at to the subspace spanned by the eigenvectors and associated with . Solutions of (2.11a–c) that satisfy (3.3) and remain in a small enough neighbourhood of for all sufficiently large converge to . Therefore, in order to study the dynamics around , 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 such that if and only if , 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 that converges to the origin as . We are interested in solutions of (2.11a–c) that satisfy (3.3) and remain in region for all . Equivalently, we seek solutions of (3.11a–c) that satisfy (3.3) and lie on a manifold , where

The following lemma characterizes such solutions that converge to the origin as (the proof corresponds, with minor modifications, to that of lemma 2.12 in [5]).

Lemma 3.10.

Up to translations in the variable , there exists a unique solution of (3.11a–c) that satisfies the asymptotic conditions (3.3), lies on the centre manifold and whose components converge to zero as , such that

Lemma 3.10 establishes the existence of at least one solution of (2.11a–c) that satisfies (3.3), stays in region and converges to as . Furthermore, this solution is uniquely determined on the centre manifold . Given that any solution of (2.11a–c) that satisfies (3.3), stays in region and converges to as must do so via and, given the monotonicity result of lemma 3.5, it is easy to prove the following lemma as in [5].

Lemma 3.11.

Given any , there exists at most one value of of the shooting parameter such that the solution of (2.11a–c) satisfying the asymptotic properties in lemma 3.1 converges to as .

Exploiting the continuity of solutions with respect to the shooting parameter, , we can extend lemma 3.11 to determine the unique value of , given , for which the solution of (2.11a–c) that satisfies (3.3) converges to as . For the proof of the following result, we refer the reader to the proof of lemma 2.14 in [5].

Lemma 3.12.

Given any , if , then the solution of (2.11a–c) satisfying the asymptotic properties in lemma 3.1 converges to as .

Using lemma 3.12 and reversing the change of variables (2.6), it is straightforward to construct a unique (up to translation) solution of system (2.1a,) that satisfies the asymptotic conditions (2.2). This leads to our first main result, which can be proved following the approach in [5]:

Theorem 3.13.

Fix . For any , system (1.3) has a smooth weak TWS connecting and . This solution is unique (up to translation), and and are monotonically strictly decreasing and increasing functions of , respectively.

Solutions converging to the equilibrium point with

In this section, we consider solutions of system (2.11a–c) subject to (3.3) that stay in region and converge to for as . Using arguments similar to those for the previous case, we can show that, for all , 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 , we denote the minimal wave speed by for each . This will enable us to draw some conclusions on the existence and uniqueness of TWS that satisfy the asymptotic conditions (2.3).

At this stage, we have no information about the possible values of for and . More specifically, given , we currently have for each . For , by remark 2.3, it is straightforward to show that for all . To characterize the minimal wave speed for , we begin by proving a non-existence result.

Lemma 3.14.

Fix and . If , then there is no such that the solution of (2.11a–c) that satisfies the asymptotic properties in lemma 3.1 converges to as .

Fix and and suppose that . We suppose for a contradiction that there exists such that the solution of (2.11a–c) that satisfies (3.3) converges to as . By the definition of , this implies that stays in region for all . Now, we can choose small enough such that and we can also find sufficiently large such that and for all . Solutions of the constant coefficient second-order ODE have infinitely many zeros in (since its characteristic equation has complex roots). Since for all , Sturm’s Comparison Theorem implies that must also have infinitely many zeros in . Therefore, exits region (and ), contradicting the assumption that .

Given and , if the minimal wave speed, , exists, then lemma 3.14 yields a lower bound for . More specifically, for all and , .

Lemma 3.15.

Fix . If , then, for all , there exists a unique such that the solution of (2.11a–c) satisfying the asymptotic properties in lemma 3.1 converges to as .

Fix and suppose that . By lemmas 3.6 and 3.9, we know that . Then, by the definition of , we have that, for any , , and the solution of (2.11a–c) satisfying (3.3) stays in the region for all by lemma 3.3. Then, by lemma 3.8, we know that, for every , the limits , , exist. In addition, by monotonicity of and with respect to (see lemma 3.5) and the fact that by lemma 3.12, we must have and We recall that and are the unique values of the shooting parameter for which the solution of (2.11a–c) that satisfies (3.3) remains in region and converges to and as , respectively. Therefore, we find that, for every , the limits for as must satisfy

We will now prove that the mapping is continuous and strictly increasing on . Choose , such that . Suppose, for a contradiction, that Irrespective of the asymptotic conditions (3.3), we can solve equation (2.11c) for and impose to obtain Any solution for which must therefore take the form (3.16). Thus, and take the form (3.16), with replaced by and , respectively. Now, by lemma 3.5, we know that for all since . We therefore have, for any , Since for all by lemma 3.5, the inequality (3.18) cannot hold and we have reached a contradiction. Since we have that , by monotonicity of solutions with respect to , and that , by the above argument, we conclude that . This proves that the mapping is strictly increasing on . Using the fact that and are, respectively, the unique values of the shooting parameter for which the solution of (2.11a–c) given by lemma 3.1 converges to and as , we have that is strictly increasing on .

We now prove that the mapping is continuous on . For fixed , (3.15) implies that and . In the proof of lemma 3.9, we performed a linear stability analysis about the equilibrium point for . We showed that, for fixed , we can find a neighbourhood of that is foliated by two-dimensional stable leaves over a one-dimensional centre manifold, which comprises equilibria of the form for . Then, any solution that enters converges to as for some that satisfies . Since converges to as , we can find large enough such that for all . By continuity of solutions with respect to , we can find such that for any such that . This implies that, for any such , converges to as , for some (since is strictly increasing with ). By our choice of , , i.e. for any such that . This proves continuity of the mapping on .

We finally show continuity at . We fix and note that, since , we can find large enough such that for all . By continuity of solutions with respect to , we can find such that for any satisfying , i.e. for any . Therefore, we have that for any . Moreover, for any , the function is strictly increasing for all and bounded above by , so . In particular, for any , we have . This proves continuity of the mapping at .

We have now shown that the mapping is strictly increasing and continuous on . Since and , application of the Intermediate Value Theorem enables us to conclude that, for any , there exists a unique such that .

Remark 3.16.

Using a similar proof to the above, we can generalize lemma 3.15 to obtain the following result. Given , suppose that there exists a unique value of the shooting parameter, , such that the solution of (2.11a–c) satisfies (3.3) and converges to as for some . Then, for all , there exists a unique value of the shooting parameter, , for which the solution of (2.11a–c) that satisfies (3.3) stays in and converges to as .

Lemma 3.15 implies that, for all , the minimal wave speed, , exists and is bounded above by . Then, given , for any , we can define

We now improve the upper bound on for by formulating a conjecture. We consider the following generalized Fisher–KPP equation with reaction term, , of Fisher–KPP type, One typically seeks TWS such that is monotonically decreasing, in which case we can invert to obtain a function . Considering the new variable , we obtain the following first-order boundary value problem (BVP): subject to . Studying TWS of (3.20) and solutions of (3.21), subject to their respective asymptotic and boundary conditions, is equivalent [22]. Moreover, it is known that if , then (3.21) subject to has a unique solution if [40-42]. Therefore, TWS of (3.20) exist and are unique if .

Returning to our original problem, by introducing and , we view the system (2.11a–c) subject to the conditions (2.3) as the following BVP: subject to the additional conditions In electronic supplementary material, S1, we show that is of Fisher–KPP type for and that if , where We conjecture that, if , then . By the preceding result for the generalized Fisher–KPP equation, this would imply that, given , the system (2.11a–c) subject to the conditions (2.3) has unique TWS for . Now, given , we let . Noting that if and only if , we formulate the following conjecture.

Conjecture 3.17.

Fix and . Given , there exists a unique such that the solution of (2.11a–c) that satisfies the asymptotic properties in lemma 3.1 converges to as . In particular, if , then .

Conjecture 3.17 implies that, given , there are values of such that the solutions of (2.11a–c) that satisfy the asymptotic conditions (2.3) behave similarly to solutions of a generalized Fisher–KPP equation with reaction term . In particular, the minimal wave speed for these TWS is defined similarly to that of a generalized Fisher–KPP equation, i.e. it is the smallest value of such that is a stable node, and not a stable spiral, for system (2.11a–c). In addition, using lemmas 3.14 and 3.15 and conjecture 3.17, we make the hypothesis that, if or, equivalently, if , then the minimal wave speed for TWS that converge to as should satisfy . In other words, in these cases, we expect that there is another mechanism that can lead to , even if is a stable node for the system (2.11a–c).

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) with evolve into travelling waves with constant propagation speed (see electronic supplementary material, S2 for corresponding travelling-wave profiles). We observe that, for , this speed is independent of , and, calculating the slopes of these lines, we find that it is approximately equal to . Additionally, when , 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 and the respective values of the propagation speed estimated using the solutions of the PDE system (results not shown). We observed that, given , the wave speed selected by the PDE appears to correspond to the smallest wave speed such that the solution of the system (2.11a–c), subject to (3.3), satisfies and converges to , .

Figure 2

We numerically solve system (1.3) on the one-dimensional spatial domain, , and impose the initial conditions (4.1) with (a) and (b). Each plot represents such that for when , with defined by (3.24), and when . We see that the front travels with a constant propagation speed that increases monotonically with . (Online version in colour.)

Now, suppose conjecture 3.17 is true. Then, given , for each and , 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 . If , then, for all , there exists a unique such that the solution of (2.11a–c) that satisfies the asymptotic properties in lemma 3.1 converges to as .

This lemma allows us to obtain a sharper upper bound on the minimal wave speed for solutions of (2.11a–c) subject to (3.3) that converge to as , where . We now summarize what we can conclude about the minimal wave speed .

Lemma 3.19.

Suppose conjecture 3.17 is true. Given , the minimal wave speed is a monotonically decreasing function on , such that

Fix . Suppose, for a contradiction, that is not a monotonically decreasing function of on . Then, we can find such that . Now, choose . Then, there exists a solution of (2.11a–c) that satisfies the asymptotic conditions (3.3), stays in region and converges to as , but there does not exist a solution of (2.11a–c) that satisfies the asymptotic conditions (3.3), stays in region and converges to as . As , remark 3.16 gives us a contradiction, hence is a decreasing function of on .

From lemma 3.14 and conjecture 3.17, we know that the minimal wave speed for all is . Since is a decreasing function of on , we must have for any . Finally, by lemma 3.14, we know that for any . This completes the proof of lemma 3.19.

While we do not have a complete characterization of the minimal wave speed for all and , 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.

Theorem 3.20.

Suppose conjecture 3.17 is true. Given , for any , there exists a minimal wave speed defined by (3.25) such that:

For , system (1.3) has no weak TWS connecting and .

For , system (1.3) has a smooth weak TWS connecting and . Moreover, this solution is unique (up to translation) and are monotonically strictly decreasing and increasing functions of , respectively.

Numerical solutions of the PDE model

In this section, we present numerical solutions of the PDE model (1.3), which complement our travelling-wave analysis. We solve (1.3) on the one-dimensional spatial domain , where , using the method of lines. A detailed description of the numerical methods employed is provided in electronic supplementary material, S2. Similarly to [16], we assume that the tumour has already spread to a position in the tissue and we impose initial conditions that satisfy, for , Here, 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 at and . We set , and for our simulations.

Remark 4.1.

Initial conditions for 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 , i.e. initial conditions of the type of (4.1) (see electronic supplementary material, S2).

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, [37-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 . 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) for 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.

Figure 3

We solve system (1.3) on the one-dimensional spatial domain, , and impose the initial condition (4.1) for . In (a), we plot such that for in the cases where and , and we see that the front travels with a strictly positive, constant propagation speed that increases monotonically with . In (b), we plot the speed of the travelling front that emerges for and observe that this speed is monotonically decreasing with , given . (Online version in colour.)

Given different values of , we calculated the speed of travelling fronts that emerge for solutions of (1.3) subject to the initial conditions (4.1) with . Our numerical simulations suggest that the wave speed selected by the PDE model is a continuous, decreasing function of , as illustrated in figure 3b, which represents this wave speed as a function of for . 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) with , appears to be equal to the minimal wave speed, , defined by (3.25). This result is interesting because the speed selected by the PDE model appears to be left-continuous at , 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, , 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 and , , we observed that the impact of decreasing the discretization step size becomes significant as approaches (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.

Comparing trajectories of the PDE and ODE models

From theorems 3.13 and 3.20, we know that system (1.3) has TWS connecting and , , for all if and for all , defined by (3.25), otherwise. Furthermore, we saw that solutions of (1.3) subject to the initial conditions (4.1) for evolve towards travelling waves and, in particular, that the wave speed is approximately equal to for . We should therefore be able to find agreement between the wave profiles of the solutions of the PDE system (1.3), subject to the initial conditions (4.1) for , and that of the ODE system (2.1a,), subject to the asymptotic conditions (2.2), if , and (2.3) otherwise, where we set to be the wave speed selected by the numerical solution of the PDE system to numerically solve the ODE system. We find good agreement between the wave profiles of the PDE and ODE solutions, and a couple of illustrative examples are shown in figure 4.

Figure 4

We compare solutions of (1.3), subject to the initial conditions (4.1) with (a) and (b), with solutions of (2.1a,), subject to the asymptotic conditions (2.2) (a) and (2.3) with (b). We use the wave speed estimated from the numerical solution of the PDE model to solve the ODE model and set . Solutions of the PDE and ODE models agree in both cases. (Online version in colour.)

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 and , for . Owing to the degeneracy in the first equation of (1.3) for , we distinguished the cases where and where .

In the first case, we proved the existence of smooth TWS for any positive wave speed, . 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 and . 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) with . 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 and 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, , defined by (3.25) for . Given , this minimal speed appears to be a monotonically decreasing, continuous function of . In particular, we conjectured that, given and , we can precisely define for . Similarly to the generalized Fisher–KPP equation, this minimal wave speed is the smallest such that the equilibrium , with , of system (2.11a–c) is a stable node and not a stable spiral. For , numerical simulations suggested that the wave speed selected by the PDE is strictly greater than , which is consistent with (3.25). The fact that the equilibrium of system (2.11a–c) is a stable node is then no longer a sufficient condition to ensure the positivity of the -component of the TWS in the desingularized variables and thus of the -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 and for TWS of (1.3) connecting and , 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 if . Yet, as , , and it is clear that, as increases, we can precisely describe the minimal wave speed for a decreasing range of values of . We would therefore like to provide a complete characterization of for all and . Now, we observed in §4a that the solution of system (1.3) subject to initial conditions (4.1) with evolves towards a travelling front with a - and -dependent wave speed. Importantly, given , it appears that this numerical wave speed is a continuous function of in , is equal to for all and is strictly greater than for all . We note that we have included in our preceding observations, which highlights our hypothesis that elucidating the minimal wave speed for (1.3) in the case could perhaps help us elucidate the minimal wave speed for (1.3) in the case , 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.