Traveling wave solutions in a model for tumor invasion with the acid-mediation hypothesis
P. N. Davis P. van Heijster R. Marangell M. R. Rodrigo































label1,label2]Paige N. Davis 
label1,label3]Peter van Heijster
label4]Robert Marangell
label5]Marianito R. Rodrigo

[label1]School of Mathematical Sciences, Queensland University of Technology, Brisbane, QLD 4000, Australia
[label2]Charles University, Faculty of Mathematics and Physics, Mathematical Institute Sokolovsk´a 83, 186 75 Prague 8, Czech Republic
[label3]Mathematical and Statistical Methods - Biometris, Wageningen University & Research, Wageningen 6708 PB, Netherlands
[label4]School of Mathematics and Statistics, University of Sydney, Sydney, NSW 2006, Australia
[label5]School of Mathematics and Applied Statistics, University of Wollongong, Wollongong, NSW 2522, Australia


In this manuscript, we prove the existence of slow and fast traveling wave solutions in the original Gatenby–Gawlinski model. We prove the existence of a slow traveling wave solution with an interstitial gap. This interstitial gap has previously been observed experimentally, and here we derive its origin from a mathematical perspective. We give a geometric interpretation of the formal asymptotic analysis of the interstitial gap and show that it is determined by the distance between a layer transition of the tumor and a dynamical transcritical bifurcation of two components of the critical manifold. This distance depends, in a nonlinear fashion, on the destructive influence of the acid and the rate at which the acid is being pumped.





Warburg effect acid-mediation hypothesis Gatenby–Gawlinski model interstitial gap geometric singular perturbation theory dynamical transcritical bifurcation.














§ INTRODUCTION


Altered energy metabolism is a characteristic feature of many solid cancer tumors and it has been recognized as a possible phenotypic hallmark <cit.>. The discovery of this altered metabolism feature dates back to the seminal work of Warburg <cit.>, who observed that certain carcinomas undergo glucose metabolism by glycolysis and not by mitochondrial oxidative phosphorylation (MOP), as normal cells do. MOP produces lactic acid as a toxic byproduct and is usually reserved for conditions of hypoxia. Paradoxically, cancer cells maintain the glycolytic phenotype even in the presence of sufficient oxygen to undergo MOP. This phenomenon is known as aerobic glycolysis or the Warburg effect.
The underlying causes of the Warburg effect still remain largely unknown. One explanation for this phenomenon is
the so-called acid-mediation hypothesis, that is, the hypothesis that tumor progression is facilitated by the acidification of the region around the tumor-host interface. This leads to a comparative advantage for tumor cells since they are more adapted to low pH environmental conditions than healthy cells. The resulting tissue degradation facilitates tumor invasion of the tissue microenvironment <cit.>.



 §.§ The Gatenby-Gawlinski model and extensions

Gatenby and Gawlinski <cit.> formulated the acid-mediation hypothesis in a reaction-diffusion framework. They proposed a reaction-diffusion system in which tumor cells produce an excess of H^+ ions due to aerobic glycolysis, which results in local acidification and thus destruction of the surrounding healthy tissue. After a suitable nondimensionalization <cit.>, the Gatenby–Gawlinski model can be written as the following system of singularly perturbed partial differential equations (PDEs) with nonlinear diffusion (in the V-component):

    {∂ U/∂τ    = U (1 - U - α W),
    ∂ V/∂τ    = β V (1 - V) + ∂/∂ x[(1 - U) ∂ V/∂ x],
    ∂ W/∂τ    = γ (V - W) + ∂^2 W/∂ x^2.
    .

Here, x ∈ and τ≥ 0 are the spatial and temporal variables, respectively. The quantities U(x,τ),  V(x,τ), and W(x,τ) represent nondimensionalized versions of the normal cell density, tumor cell density, and excess acid concentration, respectively. 
As in the quantitative discussions presented in <cit.>,  is assumed to be a small nonnegative parameter, i.e.  0 ≤≪ 1.
In addition, the constants α, β, and γ are all positive and strictly 𝒪(1) with respect to . 
The parameter α measures the destructive influence of H^+ ions on the normal tissue and therefore its value can be taken as an indicator of tumor aggressivity. For α≥ 1, solutions of (<ref>) model the situation in which total destruction of normal tissue occurs after the invasion of tumor tissue. On the other hand, for 0 < α < 1, solutions of (<ref>) correspond to the case where a residual concentration with value 1 - α of healthy tissue remains behind the spreading benign wave.

Gatenby and Gawlinski <cit.> investigated the traveling wave (TW) solutions that are compatible with (<ref>) and a number of interesting results were obtained. For instance, numerical simulations hinted at the existence of an interstitial gap (i.e. a region practically devoid of cells and located ahead of the invading tumor front) for large values of the parameter α. Subsequently, the existence of such a gap was verified experimentally; see Fig. 4 of <cit.>. 
In addition, arguments pointing toward comparatively faster invasive processes when α > 1 were provided in <cit.>.
Fasano, Herrero, and Rodrigo <cit.> further investigated the TW solutions that are compatible with (<ref>).
Using a nonstandard matched asymptotic analysis they showed that (<ref>) supports TW solutions that travel with speed 𝒪(1) and TW solutions that travel with speed 𝒪(^p) for 0 < p ≤ 1/2. 
They called the former TWs fast TW solutions and the latter TWs slow TW solutions, and the authors also obtained bounds for the wave speed in terms of the model parameters. Most notably, the authors identified slow TWs with an interstitial gap when α > 2 and the leading order width of this gap was estimated as

    z_+ = 1/√(γ)logα/2>0.

This interstitial gap ceases to exist when 0 < α≤ 2. Finally, the authors of <cit.> showed that TW solutions cannot be found when p > 1/2. See 
Fig. <ref> for a slow TW solution with an interstitial gap obtained by a numerical simulation of (<ref>). 


Different generalizations of the original Gatenby–Gawlinski model have also been proposed in the literature. 
For instance,
Holder, Rodrigo, and Herrero <cit.> included a cellular competition term in the U-equation and replaced the acid production term in the W-equation by a logistic-type reaction term. After nondimensionalization, this generalized Gatenby–Gawlinski model becomes

    {∂ U/∂τ    = U (1 - U - α (V+W)),
    ∂ V/∂τ    = β V (1 - V) + ∂/∂ x[(1 - U) ∂ V/∂ x],
    ∂ W/∂τ    = δ V(1-V) - γ W + ∂^2 W/∂ x^2.
    .

This generalization was motivated by the fact that tumors tend to present with very heterogeneous acid profiles and there is some experimental evidence of higher acid concentrations near the region of the interstitial gap. 
As a consequence of the addition of the nonlinear acid production term to the model, the profile of the excess acid concentration became pulse-like (instead of front-like in the original Gatenby–Gawlinski model; see, for instance, Fig. <ref>). The authors obtained results with regards to fast and slow TW solutions via matched asymptotic analysis similar to those in <cit.> and they also obtained estimates for the interstitial gap.


A different generalization of the Gatenby–Gawlinski model (<ref>) was given by McGillen et al. <cit.>. Here, the authors added cellular competition terms for both the U- and V-equations, as well as a term in the V-equation that incorporates acid-mediated tumor cell death. After nondimensionalization, this generalized Gatenby–Gawlinski model becomes

    {∂ U/∂τ    = U (1 - U - α_1 V - α_2 W),
    ∂ V/∂τ    = β V (1 - V) - δ_1 UV - δ_2 VW + ∂/∂ x[(1 - U) ∂ V/∂ x],
    ∂ W/∂τ    = γ(V-W) + ∂^2 W/∂ x^2,
    .

 and results analogous to those in <cit.> were derived.



 §.§ Results and outline

In this manuscript, we study the original nondimensionalized Gatenby–Gawlinski model (<ref>) and prove the formal results of <cit.> regarding the existence of fast and slow TW solutions[See the discussion in <ref> regarding using the techniques of this manuscript to analyze TW solutions found in (<ref>) and (<ref>).]. Moreover, we explain – from a mathematical perspective – the origin of the interstitial gap. We focus on the two critical cases p=0  (fast TW solutions) and p=1/2 (slow TW solutions).


We separate our results into two main theorems. 

For 0≤≪1, there exist traveling wave solutions (U_ F,V_ F ,W_ F) to (<ref>) which move with an 𝒪(1)-speed c. Upon introducing the traveling wave coordinate z = x -c τ, 
the profiles of these traveling wave solutions are, to leading order in , given by

(U_ F,V_ F ,W_ F)(x,τ) = (u_0, v_0, w_0)(z), with

    v_0(z)     = 1/1 + ^β z/c ,
    
    w_0(z)     = γ/ρ_+ - ρ_-(∫_z^∞^ρ_+ (z - ξ) v_0(ξ)  ξ̣+ ∫_-∞^z ^ρ_- (z - ξ) v_0(ξ)  ξ̣) ,
    
    u_0(z)     = c Φ_0(z)/∫_z^∞Φ_0(ξ)  ξ̣,   Φ_0(z) = ^-(1/c) ∫_0^z (1 - α w_0(ξ))  ξ̣ ,

where
ρ_± = (-c ±√(c^2 + 4 γ))/2 .

See Fig. <ref> for a fast TW solution obtained by directly simulating the Gatenby–Gawlinsky model (<ref>).





Let α∈ (0,∞)/{1,2}, then for 0≤≪1, there exist traveling wave solutions (U_ S, V_ S, W_ S) to (<ref>) which move with an 𝒪(√())-speed √()c.[The cases α =1 and α=2 are the border values for which the characteristics of the slow TW solution change, see Fig. <ref>. Therefore, they are excluded from Theorem <ref> as, for instance, for α=2 the layer transition now occurs at the same time as the transcritical bifurcation. This loss of normal hyperbolicity of the critical manifold at the layer transition complicates the proof of the theorem and is hence omitted, see <ref> for more details. That being said, we fully anticipate that the result also holds for α =1 and α=2. That is, for α=1 we expect that U=0 only in the limit x → -∞, while for α=2 the normal cell density is expected to start to grow at the tumor front.  ] 
Upon introducing the traveling wave coordinate z = x - √() c τ,
the profiles of these traveling wave solutions are, to leading order in , given by
(U_ S,V_ S ,W_ S)(x,τ) = (u, v, w)(z), with

    u(z)    = {( (1-α) + α/2^√(γ) z)_+ ,         z < 0 ,
    ( 1-α/2^-√(γ) z)_+ ,        z ≥ 0  ,
    .

where 

    (  ·  )_+ := max{ · , 0} ,

and

    w(z)    = {
    1 -1/2^√(γ) z ,         z <0 ,
    1/2^-√(γ) z ,        z ≥ 0 .
    .

The v-profile is, to leading order, given by the solution of 

    min{α/2,1}^̣2 v/ỵ^2 +c ṿ/ỵ   +β v(1-v) =0 ,


    ^̣2 v/ỵ^2 +c ṿ/ỵ   +β v(1-v) =0 ,

which connects 
v=1 as y → -∞ to v=0 as y →∞. Here, y=√() z.


In particular, these traveling wave solutions have an interstitial gap when α>2 and the leading order width of this gap is determined by (see (<ref>)) 

    1-α/2^-√(γ) z_+ = 0  z_+ = 1/√(γ)logα/2.

 
Depending on the magnitude of α, 
Theorem <ref> describes three different types of slow TW solutions, see Fig. <ref> and
note that (<ref>) is exactly the TW ODE associated to TWs in the classical Fisher–Kolmogorov–Petrovsky–Piskunov (Fisher-KPP) equation <cit.> 

    V_τ = β V (1 - V) + min{α/2,1} V_yy.






To prove Theorems <ref> and <ref> (and thus rigorously justify the asymptotic results from <cit.>), we rewrite the PDE model (<ref>) in its traveling wave framework upon introducing  (z, t) := (x - ^p c τ, τ) with p =0 or p=1/2 and with 𝒪(1)-wave speed c. TW solutions to (<ref>) now correspond to stationary solutions in this new framework and the problem reduces to studying heteroclinic orbits in an ordinary differential equation (ODE). Next, we use the multi-scale structure of (<ref>) to write this resulting ODE problem in a five-dimensional slow-fast system of first order ODEs <cit.>[Note that the slow and fast in slow-fast system is not related to the slow and fast in slow TW solution and fast TW solution. This terminology is standard in the GSPT literature and we decided not to change it.]. For the fast TW solutions there will be one fast component and four slow components, while the slow-fast splitting for the slow TW solutions is three fast components and two slow components. The details regarding the formulation of the slow-fast systems are given in <ref>.
 
We study these slow-fast systems for the fast TW solutions (see <ref>) and the slow TW solutions (see <ref>) using geometric singular perturbation theory (GSPT) <cit.>. In particular, we study the dynamics of the associated lower dimensional fast layer problems and slow reduced problems in the singular limits as 0. Next, we appropriately concatenate the dynamics of these lower dimensional systems to obtain information regarding the heteroclinic orbit – and thus fast and slow TW solutions to (<ref>) – in the singular limit as 0. Finally, we use Fenichel theory <cit.> to show that these solutions persist for positive but small . 
It turns out that for the fast TW solutions as discussed in Theorem <ref> – independent of the value of α – all the dynamics takes place on the attracting critical manifold of the slow reduced problem and the application of GSPT and Fenichel theory is straightforward. In essence, the model is a regularly perturbed problem for the fast TW solutions, and we will show that the asymptotic results of <cit.> are correct and persist for 0< ≪ 1, that is, we prove Theorem <ref>. 
See <ref> for the details.

In <ref> we prove the existence of slow TW solutions as discussed in Theorem <ref> and now the tumor aggressivity parameter α becomes important.
In particular, we have to distinguish between three cases: 0<α<1, 1<α<2, and α>2. In the first case, a slow TW solution in the singular limit 0 starts on one branch of the critical manifold (at z=-∞) and transitions through the fast layer problem (which we assume, without loss of generality, to happen at z=0) to a second branch of the critical manifold, and 
the layer dynamics will have a Fisher–KPP imprint <cit.>.
Again, we will show that such a slow TW solution persists for 0< ≪ 1 by applying GSPT and Fenichel Theory. 
In the latter two cases – 1<α<2 and α>2 – there is an additional complication related to a dynamical transcritical bifurcation of the two connected components on each branch of the critical manifold <cit.>. 
For 1<α<2, the transcritical bifurcation occurs before the fast transition through the layer problem (at z=0), while the bifurcation occurs after the transition for α>2, see Fig. <ref>. In particular, for 1<α<2 the transcritical bifurcation occurs (to leading order in ) when (1-α) + (α/2) ^√(γ) z_-=0, see (<ref>). That is, it occurs at 

    z_- = 1/√(γ)log2(α-1)/α<0 .

For α>2, the transcritical bifurcation occurs (to leading order in ) at z_+ (<ref>), see also <cit.>.
In other words, for α>2 the length of the interstitial gap is to leading order determined by the distance between the fast transition through the layer problem and the dynamical transcritical bifurcation. We conclude the manuscript with a summary and outlook regarding future projects.






§ SETUP OF THE SLOW-FAST SYSTEMS


Since we are looking for TW solutions supported by (<ref>), we introduce the traveling frame coordinates (z, t) := (x - ^p c τ, τ) for p ∈ℝ. Here, the speed c of the TW solution is assumed to be strictly 𝒪(1) with respect to . Moreover, as we are interested in waves of invasion, we assume, without loss of generality, that c > 0. A TW solution is stationary in this co-moving frame and will therefore satisfy the following system of ODEs: 

    {
    -^p c ụ/ẓ    = u (1 - u - α w),
    
    -^p c ṿ/ẓ    = β v (1 - v) + /ẓ[(1 - u) ṿ/ẓ],
    
    -^p c ẉ/ẓ    = γ (v - w) + ^̣2 w/ẓ^2,
    .

with asymptotic boundary conditions 
(u,v,w)  ((1 - α)_+,1,1) as z  -∞ and 
(u,v,w)  (1,0,0) as z ∞, see (<ref>) for the definition of (1 - α)_+.


Upon introducing the two new variables 
r := ^1 - p (1 - u) v_z+ c v (see Remark <ref>) and s := w_z, we can rewrite (<ref>) as an equivalent slow-fast system of five first order ODEs 

    {^p ụ/ẓ    = -1/c u (1 - u - α w),
    ^1 - pṿ/ẓ    = r - c v/1 - u,
    ^p ṛ/ẓ    = -β v (1 - v),
    ẉ/ẓ    = s,
    ṣ/ẓ    = -^p c s - γ (v - w).
    .

TW solutions of (<ref>) now correspond to heteroclinic orbits of (<ref>) connecting its two equilibrium points. That is,

    lim_z  -∞ (u,v,r,w,s)    =  
    ((1-α)_+,1,c,1,0) =: Z^- , 
    
    lim_z ∞ (u,v,r,w,s)    =   (1,0,0,0,0) =: Z^+ .

There are three critical p-values that balance the asymptotic scalings of (<ref>), namely, p=0, p=1/2, and p=1. In <cit.> it was shown that the case p=1 does not lead to TW solutions and we therefore do not consider this case in this manuscript (actually it was shown in <cit.> that there are no TWs for p>1/2). In addition, (<ref>) has three asymptotic scalings for 0<p<1/2. In this manuscript we consider only the cases p=0 – corresponding to fast TW solutions – and p=1/2 – corresponding to slow TW solutions.

When 0<p<1/2 the existence of slow TW solutions follows similarly to the proof for p=1/2. In <cit.> the solution profiles are obtained via asymptotic approximations on the outer (|z|≫1) and inner (z=0) regions and assuming solutions are sufficiently smooth in order to match the regions. Similar to the case when p=1/2 , it can be shown through GSPT that these solutions persist for 0<≪ 1. We refer the reader to <cit.> for more information on the procedure to apply when 0<p<1/2.

Equation (<ref>) is in its slow formulation[Recall that the slow in slow formulation is not related to the slow in slow TW solution, that is, (<ref>) is the slow formulation of the ODEs associated to both the slow TW solutions with p=1/2 and the fast TW solutions with p=0.] <cit.>.  Upon introducing the fast variable y := ^p-1z,  
the ODEs can be written in their fast formulation

    {ụ/ỵ    = -^1-2p/cu (1 - u - α w),
    ṿ/ỵ    =  r - c  v/1 -  u,
    ṛ/ỵ    = -^1-2pβ  v (1 -  v),
    ẉ/ỵ    = ^1-p s,
    ṣ/ỵ    = -^p(1-p) c s -  ^1-pγ (v - w).
    .

The slow problem (<ref>) and fast problem (<ref>) are equivalent for ≠ 0. However, they differ in the singular limit 0. 
In particular,
for the fast TW solutions, i.e. when p=0, the (u,r,w,s)-variables are slow variables and the v-variable is a fast variable. 
That is, for p=0 the slow problem (<ref>) in the singular limit 0 is a four-dimensional system of ODEs (in the slow variables) with one algebraic constraint (determined by the original equation for the fast variable). In contrast, the fast problem (<ref>) for p=0 in the singular limit 0 is a one-dimensional ODE (in the fast variable) with (up to) four additional parameters (coming from the slow equations). 
For the slow TW solutions, i.e. when p=1/2, only the (w,s)-variables are slow variables and the (u,v,r)-variables are fast variables. 



The scaling of the new variable r as
r := ^1 - p (1 - u) v_z+ c v is chosen such that -^p r_z is equal to the reaction term of the v-component in the original ODE model (<ref>). That is, -^p r_z = β v(1-v) (<ref>).
This particular scaling of r
 is inspired by a series of manuscripts <cit.> on TW solutions for chemotaxis-driven and haptotaxis-driven cell migration problems and it arises naturally when writing an extended version of (<ref>) as a singularly perturbed system of coupled balance laws. 






§ PROOF OF THEOREM <REF> ON THE EXISTENCE OF FAST TRAVELING WAVE SOLUTIONS







We start with studying the fast TW solutions supported by (<ref>) and show that the asymptotic results of <cit.> persist for 0 < ≪1. That is, we prove Theorem <ref> which states that a fast TW solution to (<ref>) is, to leading order in , given by (U_ F,V_ F ,W_ F)(x,τ) = (u_0, v_0, w_0)(z), with (u_0, v_0, w_0)(z) given in (<ref>). As eluded to above, a fast TW solution corresponds to a heteroclinic orbit in (<ref>)/(<ref>) with p=0 connecting Z^- to Z^+ (<ref>). Therefore, to prove the existence of fast TW solutions as stated in Theorem <ref> we first prove the existence of these heteroclinic orbits.

Equation (<ref>)/(<ref>) with p=0 supports a heteroclinic orbit connecting Z^- to Z^+. 


Taking p=0 in the fast system of ODEs (<ref>) and considering the singular limit 0 leads to the fast layer problem for the fast TW solutions[We rearranged the order of the equations in (<ref>) to emphasize the slow-fast structure of the problem.]

    {ṿ/ỵ    =  r - c  v/1 -  u,
    ụ/ỵ    = 0,
    ṛ/ỵ    =0,
    ẉ/ỵ    = 0,
    ṣ/ỵ    = 0.
    .

All of the variables except v are constant in (<ref>) and it can thus been seen as a single first order ODE with four additional parameters. 
It follows directly from (<ref>) that v= r/c is an equilibrium point.
Therefore, we define the four-dimensional critical manifold 

    S_ F^0 := {(u,v,r,w,s)  |    v = r/c} .

Since c>0 by assumption, we have that the critical manifold S_ F^0 is an attracting, normally hyperbolic manifold <cit.> for u<1. 
The critical manifold S_ F^0 loses normal hyperbolicity for u=1 and is repelling for u>1. 
As we will show, the u-component is always between 0 and 1 and only approaches 1 as z ∞; see (<ref>), (<ref>) and, in particular, Remark <ref>. Moreover, both asymptotic boundary conditions Z^±
(<ref>) lie on the critical manifold S_ F^0.

Taking p=0 in the slow system of ODEs (<ref>) and considering the singular limit 0 leads to the slow reduced problem for the fast TW solutions

    {
    0    = r - c v/1 - u,
    ụ/ẓ    = -1/c u (1 - u - α w),
    ṛ/ẓ    = -β v (1 - v),
    ẉ/ẓ    = s,
    ṣ/ẓ    = - c s - γ (v - w).
    .

Hence the reduced problem is a system of four first order ODEs restricted to the critical manifold S_ F^0 (<ref>). Upon imposing the algebraic constraint v=r/c, the system of four first order ODEs of (<ref>) can be written as  

    {ụ/ẓ   = -1/c u (1 - u - α w), 
    ṿ/ẓ   = -β/c v (1 - v), 
    ^̣2 w/ẓ^2 + c ẉ/ẓ - γ w    = -γ v .
    .

It was shown in <cit.> that this system – with boundary conditions as in (<ref>) – is solved by (<ref>). Hence, the u-component is strictly increasing and approaching one in the limit z ∞ <cit.>.

In the singular limit 0, the critical manifold S_ F^0 (<ref>) is normally hyperbolic and attracting in the fast direction for u<1, the asymptotic boundary conditions (<ref>) lie on S_ F^0, and the reduced problem (<ref>) restricted to the critical manifold supports the appropriate heteroclinic orbit (for which u(z)<1 for all z ∈ℝ). 
Therefore, by applying standard GSPT and Fenichel theory <cit.> (see Remark <ref>), we can conclude that this heteroclinic orbit persists in (<ref>)-(<ref>) – with p=0 – for 0 < ≪ 1. This completes the proof of Lemma <ref>. To complete the proof of Theorem <ref>, we observe that the persisting heteroclinic orbit is to leading order in  given by its singular limit. This heteroclinic orbit corresponds to the fast TWs of (<ref>) and the fast TWs are thus to leading order given by (<ref>).  This completes the proof of Theorem <ref>.



The slow problem (<ref>) and fast problem (<ref>) are – both for p=0 and p=1/2 – singular along  
{u=1}. However, u is always smaller than one, and it only approaches one in the limit z ∞, see, for instance, (<ref>) and (<ref>). A similar type of singularity is encountered in, for instance, a version of the generalized Gierer–Meinhardt model <cit.> and the Keller–Segel model <cit.>.
We refer to <cit.> for details on how GSPT and Fenichel theory can be extended to deal with this type of singularity at an asymptotic boundary condition.







§ PROOF OF THEOREM <REF> ON THE EXISTENCE OF SLOW TRAVELING WAVE SOLUTIONS


Next, we study the slow TW solutions (U_ S, V_ S, W_ S) supported by the Gatenby–Gawlinsky model (<ref>) and prove the formal asymptotic results of <cit.> and show their persistence for sufficiently small . That is, we prove Theorem <ref>.
A slow TW solution corresponds to a heteroclinic orbit in (<ref>)/(<ref>) with p=1/2 connecting Z^- to Z^+ (<ref>). Therefore, to prove the existence of slow TW solutions as stated in Theorem <ref> we first prove the existence of these heteroclinic orbits.


Equation (<ref>)/(<ref>) with p=1/2 supports a heteroclinic orbit connecting Z^- to Z^+. 





Proof.
Taking p=1/2 in the fast system of ODEs (<ref>) and considering the singular limit 0 leads to the fast layer problem for the slow TW solutions

    {ụ/ỵ    = -1/cu (1 - u - α w),
    ṿ/ỵ    =  r - c  v/1 -  u,
    ṛ/ỵ    = -β  v (1 -  v),
    ẉ/ỵ    = 0,
    ṣ/ỵ    = 0.
    .

The fast layer problem (<ref>) is again singular for u=1. However, as in the fast TW case, we will show that u-components associated to the heteroclinic orbits of interest stay smaller than one and only approach one in the limit z ∞. Therefore, this singularity does not lead to any significant complications, see Remark <ref>. 
Analysis of the equilibrium points of the layer problem (<ref>) yields a two-dimensional critical manifold S_ S^0 in ℝ^5. This critical manifold consists of two 
disjoint branches S_ S^A,B.
In turn, each of these branches consists of two connected components. In other words, the critical manifold S_ S^0 is the union of the
four two-dimensional manifolds S_ S^1,2,3,4. These four manifolds are parameterized by the slow variables (w,s) and are given by

    S_ S^A:        {
    S_ S^1   := {(u,v,r,w,s)   |    u = 0, v=0, r=0} ,
    
    S_ S^2   := {(u,v,r,w,s)   |    u = 1-α w, v=0, r=0} , 
    .
    
    
    S_ S^B:        {
    S_ S^3   := {(u,v,r,w,s)   |    u = 0, v=1, r=c} , 
    
    S_ S^4   := {(u,v,r,w,s)  |    u = 1-α w, v=1, r=c} . 
    .

The manifolds S_ S^1 and S_ S^2 intersect on S_ S^A along the line α w =1.
Similarly, S_ S^3 and S_ S^4 intersect on S_ S^B (which is disjoint from S_ S^A) along the line α w =1.
These intersections are nondegenerate in nature since α≠ 0, see Fig. <ref>. 

The three different types of slow TW solutions, see Fig. <ref>, can now be understood from the different pathways these TW solutions take through phase space along the four manifolds S_ S^1,2,3,4 in the singular limit:   



  * For 0<α < 1, the right asymptotic boundary condition Z^+ 
(<ref>) is located on S_ S^2 (as is the case for α>1), while the left asymptotic boundary condition Z^- 
(<ref>) is located on S_ S^4. Since both α and w are positive but less than 1, α w ≠ 1. As a result, the heteroclinic orbit associated to a slow TW solution starts at Z^- on S_ S^4 and transitions, via the layer dynamics, to S_ S^2. Subsequently,  it asymptotes to Z^+. 

  * For 1 < α < 2, the right asymptotic boundary condition Z^+ (<ref>) is located on S_ S^2, while the left asymptotic boundary condition Z^- (<ref>) is located on S_ S^3. The heteroclinic orbit associated to a slow TW solution thus starts at Z^- on S_ S^3, switches – via a dynamical transcritical bifurcation <cit.> – to S_ S^4 at z=z^- (<ref>) (i.e. when w(z^-)=1/α), before transitioning, via the layer dynamics, to S_ S^2.  Subsequently, it asymptotes to Z^+. 

  * For α > 2, the right asymptotic boundary condition Z^+ (<ref>) is located on S_ S^2, while the left asymptotic boundary condition Z^- (<ref>) is again located on S_ S^3. The heteroclinic orbit associated to a slow TW solution now starts at Z^- on S_ S^3, transitions, via the layer dynamics, to S_ S^1 and switches – via a dynamical transcritical bifurcation – to S_ S^2 at z=z_+ (<ref>) (i.e. when w(z_+)=1/α). Subsequently, it asymptotes to Z^+. In this case we expect to see an interstitial gap since both u and v are (to leading order) zero on S_ S^1. 

See also  Fig. <ref> for a schematic depiction of the four manifolds S_ S^1,2,3,4 (<ref>) and the three different heteroclinic orbits associated to the three different types of slow TW solutions. Finally, note that Z^- lies on the intersection of S_ S^3 and S_ S^4 for the boundary case α=1. Similarly, for α=2 the transition through the fast field occurs, in the singular limit, at the intersection of S_ S^3 and S_ S^4.







 §.§ The properties of the critical manifold

To understand the hyperbolic properties of the critical manifold S_ S^0, we compute Jacobian J of the fast equations of (<ref>)

    J = [ -1c(1-2u-α w)             0             0;   r-cv(1-u)^2         -c1-u          11-u;             0       β(2v-1)             0 ] .

The eigenvalues of the Jacobian J are given by 

    λ_1    = -1/c(1-2u-α w) ,   λ_2,3 = 1/2(1-u)(-c ±√(c^2+ 4β(2v-1)(1-u))),

with the associated eigenvectors 

    v⃗_1    = 
    (
    f(u,r,v;α, c,w), λ_1(r-cv), β(2v-1)(r-cv))^t ,       
    v⃗_2,3   = (
    0, λ_2,3, β(2v-1))^t ,

where 

    f(u,r,v;α, c,w) = (1-u)(λ_1(λ_1(1-u)+c) - β(2v-1)) .

The eigenvalues (<ref>) on the four manifolds S_ S^1,2,3,4 (<ref>) reduce to  

    S_ S^1:     λ_1^1 = -1/c(1-α w) ,       λ_2,3^1 = 1/2(-c ±√(c^2- 4β)) , 
    
    S_ S^2:     λ_1^2 = 1/c(1-α w) ,       λ_2,3^2 = 1/2α w(-c ±√(c^2- 4αβ  w)) ,
    
    S_ S^3:     λ_1^3 = -1/c(1-α w) ,       λ_2,3^3 = 1/2(-c ±√(c^2+ 4β)) ,
    
    S_ S^4:     λ_1^4 = 1/c(1-α w) ,       λ_2,3^4 = 1/2α w(-c ±√(c^2+ 4αβ w)) .

So, since the system parameters and the speed c are assumed to be positive, (λ_3^1,2,3,4)<0 on the associated manifolds. In addition, (λ_2^1,2)<0, while 
λ_2^3,4>0 (since β and αβ w are positive). The signs of the eigenvalues indicate that the fast transition, which is either from S_ S^4 to S_ S^2 or from S_ S^3 to S_ S^1, is always from a component of the manifold with two unstable eigenvalues to a component with only one unstable eigenvalue (since, as will follow from the upcoming analysis, λ_1^1,2,3,4 > 0 during the fast transition). Crucially, this latter unstable eigenvalue remains unchanged by the fast transition, i.e. λ_1^1 = λ_1^3 and λ_1^2 = λ_1^4.  
Furthermore, λ_1^1,2,3,4 have real part zero if, and only if, α w =1. Consequently, the critical manifold S_ S^0 loses normal hyperbolicity at w=1/α (i.e. where S_ S^1 coincides with S_ S^2 and  S_ S^3 coincides with S_ S^4) and this loss happens through the first eigenvalue. This loss of normal hyperbolicity is nondegenerate and transcritical in nature since α≠ 0, see Fig. <ref>. 
In other words, we have an exchange of stability between the two components on each of the two branches S_ S^A,B   
at w=1/α and the critical manifold S_ S^0 undergoes a dynamical transcritical bifurcation <cit.>.  
For α>2, this point (w=1/α) determines the rightmost point of the interstitial gap. 



We next study the slow reduced dynamics on the critical manifold S_ S^0. 
 Taking p=1/2 in the slow system of ODEs (<ref>) and considering the singular limit 0 leads to the slow reduced problem for the slow TW solutions

    {
    0    = -1/c u (1 - u - α w),
    
    0     = r - c v/1 - u,
    
    0     = -β v (1 - v),
    ẉ/ẓ    = s,
    ṣ/ẓ    =  - γ (v - w).
    .

So, the slow reduced dynamics on the four manifolds S_ S^1,2,3,4 is given by the linear equations

    ẉ/ẓ = s,   ṣ/ẓ = - γ (v^* - w),

where v^* = 0 on S_ S^1,2 and v^* = 1 on S_ S^3,4. 
These are solved by  

    w(z) = C_1^1,2^√(γ) z + C_2^1,2^-√(γ) z,   s(z) = C_1^1,2√(γ)^√(γ) z - C_2^1,2√(γ)^-√(γ) z

on S_ S^1,2, 
and 

    w(z) = 1 + C_1^3,4^√(γ) z + C_2^3,4^-√(γ) z,   s(z) = C_1^3,4√(γ)^√(γ) z - C_2^3,4√(γ)^-√(γ) z

on S_ S^3,4, for arbitrary constants C_1,2^1,2,3,4∈ℝ. These constants are determined by the asymptotic boundary conditions (<ref>) and by the dynamics of the layer problem (<ref>). Consequently, the constants are dependent on the specific α-value, see Fig. <ref>. 
To finalize the proof of Lemma <ref>, and thus Theorem <ref>, we distinguish between two different α-cases: 0<α<1 and α>1. Recall that in the former case Z^- ∈ S_ S^4, while in the latter case Z^- ∈ S_ S^3, see Fig. <ref>. 




 §.§ Proof of Lemma <ref> and Theorem <ref> for 0 < α < 1


To prove the existence of the slow TW solutions for 0<α<1, 
we first divide our spatial domain (in the slow variable z) into two slow fields I_ s^± – away from the layer dynamics – and one fast field I_ f – near the layer dynamics. In particular,

    I_ s^- :=  (-∞, -^3/8) ,   
    I_ f :=  [-^3/8,^3/8] ,   
    I_ s^+ :=  (^3/8,∞) ,

where we, without loss of generality, assumed that the layer dynamics is centered around zero. 
The asymptotic scaling ^3/8 of the boundaries of these fast and slow fields is chosen such that it is asymptotically small with respect to the slow variable z and asymptotically large with respect to the fast variable y:=^-1/2 z. In particular,  
^3/8≪ 1, while ^3/8-1/2≫ 1. 


As z  -∞ the heteroclinic orbit associated to the slow TW solution should approach Z^- (<ref>) and, hence, the critical manifold of interest is S_ S^4 for z ∈ I_ s^- (see the top left frame of Fig. <ref>). Consequently, the slow w and s components are given by (<ref>). To ensure that the solution has the correct asymptotic behavior as z  -∞ we must set C_2^4=0.  
Similarly, for z ∈ I_ s^+ the critical manifold of interest is S_ S^2 (see the bottom left frame of Fig. <ref>) and the slow w and s components are given by (<ref>) with C_1^2=0. 

During the transition through the fast field I_ f, the -dependent slow equations (w,s) are given by

    ẉ/ỵ =  s,   ṣ/ỵ = - ^1/4 c s - γ (v - w).

Therefore, and by the asymptotic scale of the fast field[^1/4≪^-(3/8-1/2).], the change of both w and s are, to leading order, constant during this transition. 
In other words, both w and s should match to leading order at zero. This determines the two remaining integration constants C_1^4 and C_2^2 and gives 

 
    w(z) = {
    1 -1/2^√(γ) z ,         z ∈ I_ s^- ,
    1/2^-√(γ) z ,        z ∈ I_ s^+ ,
    .     
    s(z) = {
    -1/2√(γ)^√(γ) z ,         z ∈ I_ s^- ,
    
     -1/2√(γ)^-√(γ) z ,        z ∈ I_ s^+ ,
    .

which coincides with (<ref>).
Hence, the fast transition always occurs at w=1/2 and the leading order profiles in the slow fields are now known (by combining (<ref>) and (<ref>)) for the five different components. In particular, 
 
    u(z) = {
    (1-α) + α/2^√(γ) z ,         z ∈ I_ s^- ,
    
     1-α/2^-√(γ) z ,        z ∈ I_ s^+ ,
    .

which coincides with (<ref>) for 0<α<1.

What remains is understanding the layer dynamics in the fast field I_ f. 
In this fast field the dynamics of the heteroclinic orbit is, to leading order, determined by (<ref>), and the orbit has to transition from S_ S^4 (where (λ_1,2^4)> 0 and (λ_3^4)<0) to S_ S^2 (where (λ_1^2)>0 and (λ_2,3^2)<0). Since w is to leading order constant in the fast field,
the u-equation of (<ref>) is of logistic-type and, by (<ref>), u=1-α w on both S_ S^2,4. 
Consequently, and since the logistic equation does not support pulse-type solutions, u is also constant during the fast transition. In particular, u=1-α w = 1-α/2 in I_ f, see (<ref>).  
The resulting (v,r)-equations (<ref>) – with u=1-α/2>0 – can be written as (<ref>) with the observation that min{α/2,1} = α/2 since 0<α<1.
This is exactly the TW ODE associated to TWs in the classical Fisher–KPP equation[This does not come as a surprise since the V-component of the original PDE (<ref>), in the fast variable y and for U=1-1/2α, is the Fisher–KPP equation
V_τ = β V (1 - V) + α/2 V_yy.]. Hence, there exists a heteroclinic connection between (v,r)=(1,0) and (v,r)=(0,0) in the fast field, see <cit.> and references therein. In addition, the (v,r)-components are nonnegative during this transition if, and only if, c ≥ c_min:=√(2 αβ)[The expression for c_min also arose from the formal analysis of <cit.>.] – the so-called minimum wave speed of the associated Fisher-KPP equation – see, for instance, <cit.> and references therein. The last observation also follows directly from the fact that λ_2,3^2 (<ref>) – with w=1/2 – are complex-valued for c < c_min. Moreover, observe that the first components of the eigenvectors v⃗_2,3 (<ref>) associated to λ_2,3 are zero, that is, the u-component indeed does not change during the fast transition. This completes the analysis of the layer problem, and hence the analysis of the heteroclinic orbits for 0 < α < 1, in the singular limit 0.


We show the persistence of the singular heteroclinic orbits for sufficiently small  in (<ref>)-(<ref>) (with p=1/2) and thus the existence of slow TW solutions in (<ref>). By (<ref>), a singular orbit only approaches u=1 in the limit z ∞ (see also Remark <ref>). Furthermore, as 0<α<1 and as w is given by (<ref>), we have that α w ≠ 1 along the singular orbit. Therefore, the critical manifold S_ S^0 does not lose normal hyperbolicity along the singulars orbit and each singular orbit is a heteroclinic connection between two normally hyperbolic components of the critical manifold. Fenichel's First Persistence Theorem <cit.> states that, for  small enough (and after appropriately compactifying S_ S^2 and S_ S^4), there exist locally invariant slow manifolds S_ S,^2 and S_ S,^4 in the full -dependent system (i.e. (<ref>)-(<ref>) with p=1/2) that are 𝒪(√())-close to S_ S^2 and S_ S^4, respectively. 
Observe that Z^± (<ref>) are independent of  and, hence, S_ S,^2,4 coincide with S_ S^2,4 in the asymptotic limits z ±∞.
Fenichel's Second Persistence Theorem <cit.> states that the full -dependent system also possesses locally invariant stable and unstable manifolds 𝒲^u(S_ S,^4) and 𝒲^s(S_ S,^2) which are 𝒪(√())-close to the stable and unstable manifolds 𝒲^u(S_ S^4) and 𝒲^s(S_ S^2), respectively. 
We also have the necessary property of the singular problem that the heteroclinic connections (singular orbits) are contained in the intersection 𝒲^u(S_ S^4)∩𝒲^s(S_ S^2) and it follows that the orbit persists (in the intersection of 𝒲^u(S_ S,^4)∩𝒲^s(S_ S,^2)) for 0<≪1 if the intersection 𝒲^u(S_ S^4)∩𝒲^s(S_ S^2) is transversal, see <cit.>.

 
The slow TW problem has three fast variables (u,v,r) and two slow variables (w,s). Moreover, for 0<α<1, (λ_1^2)>0 and (λ_2,3^2)<0, see (<ref>). Therefore, 
(𝒲^s(S_ S,^2)) = (𝒲^s(S_ S^2))=2+2=4.[The first “2” originates from the number of eigenvalues (<ref>) on S_ S^2 with negative real part (i.e the number of fast stable eigenvalues), while the second “2”  comes from the number of slow variables.]
Similarly, (λ_1,2^4)>0 and (λ_3^4)<0 and, consequently, (𝒲^u(S_ S,^4)) =(𝒲^u(S_ S^4))=2+2=4. Generically, two four-dimensional objects in a five-dimensional phase space intersect transversally. 
The transversality of the intersections is typically shown through a Melnikov integral <cit.>. However, for this specific system, we take advantage of the additional structures of the problem. We define the so-called take-off curve as the unstable direction from which the singular orbit leaves Z^- on S_ S^B, the jump point as the point on the take-off curve where a solution leaves the critical manifold to make the fast transition, and the touchdown curve as the union of points on S_ S^A a solution could land on after the fast transition. Due to the fact that u,w,s are, to leading order, constant across the fast transition, the touchdown curve is the projection of the take-off curve onto S_ S^A. The existence of an orbit relies on the fact that the touchdown curve intersects the stable direction of Z^+ and it is clear this intersection is transversal, see Fig. <ref>. The fact that this stable direction intersects the touchdown curve transversally is an indicator that the intersection 𝒲^u(S_ S^4)∩𝒲^s(S_ S^2) is also transversal.
Furthermore, during the fast transition, i.e. in the intersection 𝒲^s(S_ S^4)∩𝒲^u(S_ S^4), u is constant and the dynamics during this transition are controlled by a Fisher-KPP-type equation (<ref>) whose end state (in the two-dimensional state space (v,r)) has no unstable directions and supports a continuous family of TWs in c, implying the persistence of solutions under an  perturbation. We exploit these structures in order to prove the transversality of the intersection 𝒲^u(S_ S^4)∩𝒲^s(S_ S^2).

We first analyse the behaviour of the 4-dimensional stable subspace 𝒲^s(S_ S^2) and observe that the tangent space T𝒲^s(S_ S^2) at points in S_ S^2 is spanned by the four vectors (0 , λ_2,3^2 , -β , 0, 0)^T, ((1-α)_+ , 0 , 0 , 1 , 0)^T, (0 , 0 , 0 , 0 , 1)^T. 
The first three elements of the vectors (0 , λ_2,3^2 , -β , 0, 0)^T are the stable eigenvectors v⃗_2,3 respectively, see (<ref>), of the Jacobian evaluated on S_ S^2 appended with two 0 components representing w,s – components which remain constant across the fast transition. The latter vectors ((1-α)_+ , 0 , 0 , 1 , 0)^T, (0 , 0 , 0 , 0 , 1)^T are the span of the manifold S_ S^2. 
Of the vectors that span T 𝒲^s(S_ S^2) only (0 , λ_2,3^2 , -β , 0, 0)^T will change under the evolution along the layer fiber. This is because the layer transition is governed by a Fisher-KPP-type equation in v,r, and the other components are to leading order constant. 
Additionally, as the end state of the Fisher-KPP equation has no unstable directions the space spanned by these two vectors will always contain the space spanned by (0,1,0,0, 0)^T and (0,0,1,0,0)^T, i.e. the basis vectors of the (v,r) phase space.
Furthermore, v⃗_1∈𝒲^u(S_ S^4) and v⃗_1→ (f(1-α,1,c,α,c,1/2),0,0) as the orbit approaches S_ S^4 in backwards z. Thus, v⃗_1, appended with zeros for w,s, is in the tangent space T 𝒲^u(S_ s^4) and is proportional to (1,0,0,0,0)^T. This vector is linearly independent to the four vectors that span T 𝒲^s(S_ S^2). At any point along the layer fibre, the combined tangent spaces of 𝒲^s(S_ S^2) and 𝒲^u(S_ S^4) contain the full tangent space to ℝ^5.
From this, it follows directly that the intersection is transversal and the heteroclinic connection persists for 0<≪1  <cit.>. Consequently, (<ref>) supports slow TW solutions for 0 < α < 1 and for sufficiently small . This completes the proof of Lemma <ref> and Theorem <ref> for 0<α < 1.



 §.§ Proof of Lemma <ref> and Theorem <ref> for α > 1


For α>1 the situation is more involved since a dynamical transcritical bifurcation of critical manifolds is involved (when α w =1), see Fig. <ref>. This critical bifurcation occurs to the left of the layer transition (at z=0) for 1< α< 2, while it occurs to the right of the layer transition for α> 2. The latter case results in an interstitial gap only because part of the heteroclinic orbit is on S_ S^1 where both u, representing the normal cell density, and v, representing the tumor cell density, are zero to leading order. 
However, in both cases we can still use the same slow-fast splitting of the spatial domain (<ref>) in the singular limit 0. Furthermore, the layer problem still exhibits Fisher–KPP type behavior described by (<ref>). 

In more detail, since α>1 the heteroclinic orbit associated to the slow TW solution should approach Z^- ∈ S_ S^3, see (<ref>) and (<ref>), as z  -∞. Hence, the critical manifold of interest is S_ S^3 (<ref>) for -z ≫ 1. Consequently, the slow w and s components are given by (<ref>) and – to ensure that the solution has the correct asymptotic behavior – C_2^3=0.  That is, 

    w(z) = 1 + C_1^3^√(γ) z ,     s(z) = C_1^3√(γ)^√(γ) z ,   for   -z ≫ 1.

Similarly, for z ∈ I_ s^+ the critical manifold of interest is S_ S^2 (since Z^+ ∈ S_ S^2) and the slow w and s components are given by (<ref>) with C_1^2=0: 

    w(z) = C_2^2^-√(γ) z ,     s(z) = -C_2^2√(γ)^-√(γ) z ,    for   z ∈ I_ s^+.

The two critical manifolds S_ S^2,3 both undergo a (different) dynamical  trans­critical bifurcation at α w =1. If this bifurcation occurs at z=ž <0 (to the left of the layer transition at z=0) then the heteroclinic orbit passes from S_ S^3 onto S_ S^4. In contrast, if this bifurcation occurs at z=ẑ>0 (to the right of the layer transition) then the heteroclinic orbit transitions from S_ S^1 onto S_ S^2.

In the former case where the transition occurs at z=ž <0, we get that the slow w and s components after the transition are given by 

    w(z) = 1 + C_1^4^√(γ) z + C_2^4^-√(γ) z ,    s(z) = C_1^4√(γ)^√(γ) z - C_2^4√(γ)^-√(γ) z , 
      for   z ∈ I_ s^-   and   z>ž,

see (<ref>). However, by construction, the slow components should match as z approaches ž. So, from combining (<ref>) and (<ref>), we get  

    w(z) = 1 + C_1^3^√(γ) z ,    s(z) = C_1^3√(γ)^√(γ) z ,   for   z ∈ I_ s^-,

see Fig. <ref>.
Since the change of both w and s are, to leading order, constant during the transition through the fast field I_ f, see (<ref>), if follows that
(<ref>) and (<ref>) should match as z approaches zero. Thus, similar to the case α<1, the slow components are given by (<ref>)/(<ref>).

Hence, ž∈ I_ s^- such that α w(ž)=1 is given by
ž =γ^-1/2log(2(α-1)/α) =: z_- (<ref>), and ž is negative only for 1<α<2. That is, the dynamical transcritical bifurcation occurs only to the left of the layer transition, and the heteroclinic orbit transitions from S_ S^3 to S_ S^4, if 1<α<2. See also Fig. <ref> and Fig. <ref>.
As before, the leading order profiles in the slow fields are now known for all the components, and, in particular,
 
    u(z) = {
    0 ,         z < z_- ,
    
    (1-α) + α/2^√(γ) z ,          z > z_-  and  z∈ I_ s^- ,
    
     1-α/2^-√(γ) z ,         z∈ I_ s^+ ,
    .

which coincides with (<ref>) for 1<α<2.

 

We proceed in a similarly fashion in the case where the bifurcation occurs to the right of the layer transition at z=ẑ>0. Again, we obtain that the slow components in the slow fields are given by (<ref>)/(<ref>). Consequently, ẑ∈ I_ s^+ such that α w(ẑ)=1 is given by
ẑ = γ^-1/2log(α/2) =: z_+ (<ref>), and ẑ is positive only for α>2. That is, the dynamical transcritical bifurcation only occurs to the right of the layer transition and the heteroclinic orbit transitions from S_ S^1 to S_ S^2, if α>2, see Fig. <ref> and Fig. <ref>. For α>2 the positive value of ž corresponds to the existence of the interstitial gap. Furthermore, the value of z_+ indicates the width of the interstitial gap and thus confirms the estimate of the width from <cit.>. The leading order profiles in the slow fields are now known and the u-component is given by
 
    u(z) = {
    0 ,          z∈ I_ s^- ,
    
    0 ,         z < z_+  and  z∈ I_ s^- ,
    
     1-α/2^-√(γ) z ,          z>z_+,
    .

which coincides with (<ref>) for α>2.

For both 1<α<2 and α>2, the layer dynamics in the fast field I_ f is the same as for 0<α<1 in <ref>. That is, due to the logistic nature of the u-component in (<ref>) and the particulars of the critical manifolds involved, the fast u-component actually does not change during the transition through the fast field I_ f. Consequently, the layer transition is still associated to 
the Fisher–KPP equation (<ref>). The difference between 1<α<2 and α>0 arises from the fact that u=1-α/2 during the transition for 1<α<2, while u=0 during the transition for α>2, see (<ref>)) and (<ref>)).
Hence, in both cases there exists a heteroclinic connection between (v,r)=(1,0) and (v,r)=(0,0) in the fast field. The (v,r)-components are nonnegative for 1<α<2 if, and only if, c ≥ c_min:=√(2 αβ) (i.e. λ_2,3^2 (<ref>) are real-valued). In contrast, the (v,r)-components are nonnegative for α>2 if, and only if, c ≥c̅_min:=2√(β) (i.e. λ_2,3^1 (<ref>) are real-valued). This completes the analysis of the layer problem, and hence the analysis of the heteroclinic orbits in the singular limit 0, for α>1.







We show the persistence of the singular heteroclinic orbits for sufficiently small  in (<ref>)-(<ref>) (with p=1/2) and thus the existence of slow TW solutions in (<ref>). 
The added complexity – compared to the 0<α<1 case discussed in <ref> – is related to showing the persistence of the transcritical dynamical bifurcation structure around α w = 1 since the critical manifold S_ S^0 loses normal hyperbolicity here. In addition, as in the 0<α<1 case, the persistence of solutions across the fast transition will be shown.


The transcritical singularity results from the self-intersection of the critical manifold along the line α w=1. The persistence of the transcritical dynamical bifurcation structure around α w = 1 follows from the observation that u=0 is invariant for the full -dependent system ((<ref>) with p=1/2). Hence, we have u=0 on the perturbed manifolds S_ S,^1,3. Furthermore, away from α w=1 the perturbed manifolds S_ S,^2,4 are, to leading order, given by S_ S^2,4. Therefore, the intersection between S_ S^4 and S_ S^3 and the intersection between S_ S^2 and S_ S^1 must persist in the full -dependent system.

The persistence of singular orbits across the fast transition for 0<≪ 1 is shown by proving the transversality of the intersection 𝒲^u(S_ S^4)∩𝒲^s(S_ S^2) for 1 < α < 2, and the transversality of the intersection 𝒲^u(S_ S^3)∩𝒲^s(S_ S^1) for α > 2. The argument follows similarly to the 0<α<1 case. The fast transition is governed by a Fisher-KPP-type equation (<ref>) in each case and one can explicitly calculate the spanning vectors of the relevant tangent spaces in order to prove that the combined tangent spaces (of 𝒲^u(S_ S^4) and 𝒲^s(S_ S^2) for 1 < α < 2 and of 𝒲^u(S_ S^3) and 𝒲^s(S_ S^1) for α > 2) contain the full tangent space to ℝ^5. Hence, the intersection is transversal in each case and the heteroclinic connections persists for 1<α<2 and α>2  <cit.>. Consequently, (<ref>) supports slow TW solutions for 1<α<2 and α>2  for sufficiently small . 
This completes the proof of Lemma <ref> and Theorem <ref> for α > 1.




§ SUMMARY AND OUTLOOK


In this manuscript, we analyzed TW solutions supported by the nondimensionalized Gatenby–Gawlinski model (<ref>). This model was originally proposed by Gatenby and Gawlinski in <cit.> to investigate the acid-mediation hypothesis of the Warburg effect, also known as aerobic glycolysis <cit.>. 
This hypothesis postulates that this Warburg effect is caused by the fact that the progression of certain tumors is facilitated by the acidification of the region around the tumor-host TW interface and this leads to an advantage of the tumor cells <cit.>.
In the model, the acid-mediation hypothesis is characterized by an interstitial gap, a region in front of the invading TW interface devoid of cells, see also Fig. <ref> and Fig. 4 of <cit.>. 
The TW solutions of (<ref>) have been analyzed numerically in <cit.> and by using formal matched asymptotics in <cit.>. In particular, in <cit.> it was shown that the Gatenby–Gawlinski model (<ref>) supports slow and fast TW solutions. Here, “slow” and “fast” refer to the asymptotic scaling of the speed c of a TW solution with respect to the small parameter  (that measures the strength of the nonlinear diffusion of the tumor). 

In this manuscript, we embedded the TW problem associated to (<ref>) into a slow-fast[Here, slow-fast refers to the difference in asymptotic scaling of the (nonlinear) diffusion coefficient of (<ref>)] structure and use geometric singular perturbation techniques to prove Theorems <ref> and <ref> – and thus prove the formal results of <cit.> – in the critical cases (c ∼𝒪(1) and c ∼𝒪(√()) respectively). In particular, we showed that the interstitial gap is present only if the destructive influence of the acid, modeled by the parameter α in (<ref>), is strong enough. That is, the interstitial gap exists only for α>2, see also <cit.>. We showed that, from a geometric perspective, the interstitial gap can be understood as the distance between the TW interface – which has the characteristics of a 
Fisher–KPP wave – and a dynamical transcritical bifurcation of two parts of the critical manifold.  
For moderate strengths of the destructive influence of the acid, i.e. for 1<α<2, parts of the critical manifold involved still undergo a 
dynamical transcritical bifurcation, however, this now occurs behind the TW interface and no region devoid of cells is thus created, see, for instance, the middle panel of Fig. <ref>. 




The size of the interstitial gap 

    z_+ = 1/√(γ)logα/2

depends explicitly on the destructive influence α of the acid, scales with 1/√(γ) related to the relative absorption rate of the H^+ ions, but is independent of β which is the relative growth rate of the neoplastic tissue. A priori, this maybe comes as a surprise, but not after a closer heuristic inspection of the model equations (<ref>).
The interstitial gap is the absence of cells ahead of the invading tumor front – which is located at 0 by construction and without loss of generality – and 
the normal cell density is modelled by 

    ∂ U/∂τ  = U (1 - U - α W) .

The equilibrium points are given by U(x)=0 and U(x) = 1 - α W(x) and these exactly coincide, i.e. the transcritical bifurcation occurs, when x=x̅ is such that 1 = α W(x̅). 
Hence, the size of the interstitial gap x̅ is fully expected to explicitly depend on α, but also on the acid concentration W. 
The equation for the acid concentration is  

    ∂ W/∂τ  = γ (V - W) + ∂^2 W/∂ x^2,

and since the reaction term scales with γ we expect the interstitial gap to scale with 1/√(γ). While the equation for the acid concentration also depends on the tumor cell density V, the scale separation 
in the equation 
for the tumor cell density 

    ∂ V/∂τ  = β V (1 - V) + ε∂/∂ x[(1 - U) ∂ V/∂ x] ,

enforces that tumor cell density only plays a crucial dynamic role near the interface of the invading tumor front. Therefore, it does not influence the size of the interstitial gap, which is thus independent of β. 
Finally,
for an interstitial gap we need the normal cell density U to be 0 ahead of the tumor front V (located at the origin). That is, we require the transcritical bifurcation to occur at x=x̅>0. In other words, we want the solution of  
1 = α W(x)  to have a positive x-value. Since W is decreasing – the acid concentration is high behind the front put low ahead of the front – we expect that larger α-values lead to larger interstitial gaps and that there potentially is a lower-bound for which the interstitial gap does not exists. 

These heuristic arguments can of course also be used to predict the potential existence of the interstitial gap in other models – such as the generalized models (<ref>) and (<ref>) studied in <cit.>, respectively <cit.>. The Gatenby–Gawlinski model (<ref>) is amendable for analysis because the nonlinear diffusion term in the equation for the tumor cells acts as a regular perturbation to the normal diffusion term (as U is constant to leading order during the fast transition), and the underlying equation has a Fisher-KPP imprint. A simplified model, obtained via a quasi-steady state reduction <cit.> of the full model, is given by

    {∂ U/∂τ    = U (1 - U - α W),
    ∂ W/∂τ    = γ (H(-x) - W) + ∂^2 W/∂ x^2,
    .

where H(·) is the Heaviside step-function replacing the V-component of (<ref>).
This simplified model has similar characteristics to the full model (<ref>), and, crucially, 
still supports TW solutions with an interstitial gap of length z_+ for α>2.
Other dynamical expressions for the tumor cell density with the same properties are expected to yield the same, or at least similar, results. For instance, the tumor cell density in the generalized model (<ref>) <cit.> is given by 

    ∂ V/∂τ  = β V (1 - V) - δ_1 UV - δ_2 VW + ε∂/∂ x[(1 - U) ∂ V/∂ x].

Since we still have the scale separation in the model, it is anticipated that the tumor front still has a sharp front, while the normal cell density U and acid concentration W vary more gradually. Consequently, both the normal cell density U and acid concentration W are effectively constant near this front. Therefore, for δ_1 and δ_2 not too large, we still expect to see an invasive front with similar characteristics like an interstitial gap. However, behind the front the tumor density will not be 1, as for (<ref>), but 1 -  δ_1 U̅- δ_2 W̅, where U̅ and W̅ are the constant values of U and W near the front. So, while the details will be different, we fully expect the existence of an interstitial gap for this model (since also the characteristic of the normal cell density equation and acid concentration equation are not significantly different in (<ref>)). 


Similarly, an “inhibitor-type” acid concentration equation with a more general reaction term that increases as function of the tumor cell density and does not depend on the normal cell concentration, coupled with a normal cell concentration equation with a reaction term of the form U f(U,V,W) with two intersecting and interacting roots, is fully expected to have similar characteristics. Looking at the generalized model (<ref>) <cit.>, we observe that the equations for the normal cell density U, as well as the tumor cell density V, is still in the prescribed generalized form. However, the reaction term in the acid concentration equation 

    ∂ W/∂τ  = δ V(1-V) - γ W + ∂^2 W/∂ x^2

does not increase as function of the tumor cell density for V ≈ 1. Therefore, it is expected that some of the characteristics of (<ref>) will be different from (<ref>), as is observed in <cit.>. For instance, the acid concentration profile is not front-shaped but pulse-shaped. 



The results of this manuscript show that the Gatenby–Gawlinski model (<ref>) supports, even for a fixed parameter set, a myriad of TW solutions with different speeds. 
A logical next question to answer is related to wave speed selection. That is, given a specific parameter set and initial condition, what is – if any – the speed of the TW solution the initial condition converge to? 
Because of the Fisher-KPP imprint of the V-component of the model, it can be expected that a dispersion relation relating the asymptotic behavior of an initial condition around plus infinity and the linear spreading speed of the TW solution can be derived, see, for instance, <cit.>. However, a TW solution will not always travel with this linear spreading speed, see, for instance, <cit.>. It is also interesting to see if the observed wave speeds for the slow TW solutions equal the minimum wave speeds of the associated Fisher-KPP equations (c_min:=√(2 αβ) for 0 < α < 2 and c̅_min:=2√(β) for α>2, see <ref>). That is, are the observed slow TW solutions pushed or pulled fronts <cit.>? 
A first natural step to start tackling these questions is to study the stability properties of the slow and fast TW solutions, and
a potential approach is to combine the analytic approach used in <cit.> (to study the spectral stability of TW solutions in a Keller–Segel model) with the Ricatti Evans function approach developed in <cit.> to numerically compute eigenvalues.
This is part of future work. 






Finally, while we only rigorously establish the existence of slow and fast TW solutions to the original Gatenby–Gawlinski model (<ref>), the methodology of embedding the problem into a slow-fast structure and subsequently studying the dynamics of the reduced and layer problems can also be used to prove the existence of TW solutions in generalizations of the Gatenby–Gawlinski model (such as models (<ref>) and (<ref>) studied in <cit.>, respectively <cit.>).  
The argument for the persistence of solutions across the dynamical transcritical bifurcation for 0<≪1 follows from the invariance of u=0 in the full -dependent system (<ref>). 
A mathematically interesting question is whether this dynamical transcritical bifurcation also persists for similar systems where this invariance is broken, see <cit.>.






§ ACKNOWLEDGEMENTS

PD and PvH acknowledge support under the Australian Research Council's Discovery Early
Career Researcher Award funding scheme DE140100741. PD acknowledges the support of the grant No. 20-11027X financed by Czech Science Foundation (GAČR). PvH and RM acknowledge support under the Australian Research Council's Discovery Project DP200102130.
The authors would like to thank Martin Wechselberger, Hinke Osinga, and Bernd Krauskopf for their insightful and productive discussions. 
PvH and RM would also like to thank the University of Wollongong for their hospitality. 
This research was initiated during the first Joint Australia-Japan Workshop on Dynamical Systems with Applications in Life Sciences at Queensland University of Technology.

















99


P1 P. N. Davis, P. van Heijster, and R. Marangell (2017), Absolute instabilities of travelling wave solutions in a Keller–Segel model, Nonlinearity 30, pp. 4029–4061. 

P2 P. N. Davis, P. van Heijster, and R. Marangell (2018), Spectral stability of travelling wave solutions in a Keller–Segel model, Appl. Numer. Math. (in press).  


DGK A. Doelman, R. A. Gardner, and T. J. Kaper (2001), Large stable pulse solutions in reaction-diffusion equations, Indiana U. Math. J., 50, pp. 443–507.

FHR A. Fasano, M. A. Herrero, and M. R. Rodrigo (2009), Slow and fast invasion waves in a model of acid-mediated tumour growth, Math. Biosci., 220, pp. 45–56. 

F N. Fenichel (1979), Geometric singular perturbation theory for ordinary differential equations, J. Differ. Equations, 31, pp. 53–98.
	
KPP A. Kolmogorov, I. Petrovsky, and N. Piscounov (1937). Étude de l'equation de la diffusion avec croissance de lat quantité de matière et son application à un problèm biologique. Moscow Univ. Math. Bull., (1) pp. 1–25.
	
GG1 R. A. Gatenby and E. T. Gawlinski (1996), A reaction-diffusion model for cancer invasion, Cancer Res., 56, pp. 5745–5753.

GG2 R. A. Gatenby and R. J. Gillies (2004), Why do cancers have high aerobic glycolysis?, Nat. Rev. Cancer, 4, pp. 891–899.

HW D. Hanahan and R. A. Weinberg (2011), Hallmarks of cancer: the next generation, Cell, 144, pp. 646–674.

H1 K. Harley, P. van Heijster, R. Marangell, G. J. Pettet, and M. Wechselberger (2014), Existence of traveling wave solutions for a model of tumor invasion, SIAM J. Appl. Dyn. Syst., 13, pp. 366–396.

H2 K. Harley, P. van Heijster, R. Marangell, G. J. Pettet, and M. Wechselberger (2014), Novel solutions for a model of wound healing angiogenesis, Nonlinearity, 27, pp. 2975–3003.

H3K. Harley, P. van Heijster, R. Marangell, G. J. Pettet, and M. Wechselberger (2015), Numerical computation of an Evans function for travelling waves, Math. Biosci. 266, pp. 36–51.

HvHP K. Harley, P. van Heijster, and G. J. Pettet (2014), A geometric construction of traveling wave solutions to the Keller–Segel model, ANZIAM J., 55 (EMAC2013), pp. C399–C415.





H G. Hek (2010), Geometric singular perturbation theory in biological practice, J. Math. Biol., 60, pp. 347–386.

HRH A. B. Holder, M. R. Rodrigo, and M. A. Herrero (2014), A model for acid-mediated tumour growth with a nonlinear acid production term, Appl. Math. Comp., 227, pp. 176–198.

J C. K. R. T. Jones (1995), Geometric singular perturbation theory, in Dynamical systems (Montecatini Terme, 1994), volume 1609 of Lecture Notes in Math., pp. 44–118. Berlin: Springer.

K T. J. Kaper (1999), An introduction to geometric methods and dynamical systems theory for
singular perturbation problems, in Analyzing multiscale phenomena using singular perturbation methods (Baltimore, MD, 1998), volume 56 of Proc. Sympos. Appl. Math., pp. 85–131. Providence, RI: Amer. Math. Soc.

KS M. Krupa and P. Szmolyan (2001), Extending slow manifolds near transcritical and pitchfork singularities, Nonlinearity, 14, pp. 1473–1491.
 


K15 C. Kuehn (2015), Multiple Time Scale Dynamics (Vol. 191), Berlin: Springer.

LAR D. A. Larson (1978), Transient bounds and time-asymptotic behavior of solutions to nonlinear equations of
Fisher type, SIAM J. Appl. Math., 34, pp. 93–103.

LIU W. Liu, D. Xiao, and Y. Yi (2003), Relaxation oscillations in a class of predator–prey systems, J. Differ. Equations, 188, pp. 306–331.

MGMM J. B. McGillen, E. A. Gaffney, N. K. Martin, and P. K. Maini (2014), A general reaction-diffusion model of acidity in cancer invasion, J. Math. Biol., 68, pp. 1199–1224.

MCK H. P. McKean (1975), Application of Brownian motion to the equation of Kolmogorov–Petrovskii–Piskunov,
Comm. Pure Appl. Math., 28, pp. 323–331.

MUR J. D. Murray (2002),  Mathematical Biology. I, volume 17 of Interdisciplinary Applied Mathematics, Springer-Verlag, New York.

R C. Robinson (1983), Sustained resonance for a nonlinear system with slowly varying coefficients, SIAM J. Math. Anal., 5, pp. 847–860.

SAAR W. van Saarloos (2003), Front propagation into unstable states, Phys. Rep., 386, pp. 29–222.

S1 L. Sewalt, K. Harley, P. van Heijster, and S. Balasuriya (2016), Influences of Allee effects in the spreading of malignant tumours, J. Theor. Biol., 394, pp. 77–92.

S P. Szmolyan (1991), Transversal heteroclinic and homoclinic orbits in singular perturbation problems, J. Differ. Equations, 92, pp. 252–281.


W O. Warburg (1930), The Metabolism of Tumors, Arnold Constable, London.


WP M. Wechselberger and G. J. Pettet (2010), Folds, canards and shocks in advection-reaction-diffusion models, Nonlinearity, 23, pp. 1949–1969.

ZZ W. Zhang, V. Kirk, J. Sneyd, and M. Wechselberger (2011), Changes in the criticality of Hopf bifurcations due to certain model reduction techniques in systems with multiple timescales,  J. Math. Neurosci., 1.





