Propagation fronts in a simplified model of tumor growth with degenerate cross-dependent self-diffusivity

https://doi.org/10.1016/j.nonrwa.2021.103387Get rights and content

Abstract

Motivated by tumor growth in Cancer Biology, we provide a complete analysis of existence and non-existence of invasive fronts for the reduced Gatenby–Gawlinski model tU=Uf(U)dV,tV=xf(U)xV+rVf(V),where f(u)=1u and the parameters d,r are positive. Denoting by (U,V) the traveling wave profile and by (U±,V±) its asymptotic states at ±, we investigate existence in the regimes d>1:(U,V)=(0,1)and(U+,V+)=(1,0),d<1:(U,V)=(1d,1)and(U+,V+)=(1,0),which are called, respectively, homogeneous invasion and heterogeneous invasion. In both cases, we prove that a propagating front exists whenever the speed parameter c is strictly positive. We also derive an accurate approximation of the front profile in the singular limit c0.

Introduction

Biological invasion is one of the basic features of Nature and its potentiality to modify its structure and its inherent vitality. Sometimes invasion of a new species can be regarded as a positive event, sometimes as a negative one depending on the property of the intruder and the invaded, see [1]. Here, motivated by Cancer Biology, we focus on the appearance of invasion in the form of propagating fronts in tumor growth. Precisely, we present a rigorous mathematical analysis of a reaction–diffusion system composed by two differential equations, for which we prove the existence of traveling wave solutions that can be interpreted as invasion fronts of a cancerous tissue into a healthy one. We urge the reader to pay attention to the specific form of the nonlinear diffusion term in our system and the consequences it has on the set of admissible propagation speeds.

The original motivation is the analysis of the so-called acid-mediated tumor growth, proposed by Otto Warburg as a mechanism responsible for tumor increase [2]. Precisely, the so-called Warburg effect refers to the observation that – even in aerobic conditions – cancer cells tend to favor metabolism via glycolysis rather than the more efficient oxidative phosphorylation pathway, usually preferred by most other cells of the body [3]. A simplified mathematical description for such a mechanism has been proposed in [4]. After an appropriate rescaling and using the notation f(s)=1s, the one-dimensional version of the model reads as tU=U{f(U)dW},tV=x{f(U)xV}+rVf(V),tW=ax2W+b(VW),where U=U(x,t) represents the (normalized) population of healthy cells, V=V(x,t) is the (normalized) population of tumor cells, and W=W(x,t) is the concentration of lactic acid. The reaction–diffusion system (1.1) is a sound description of the acid-mediated tumor growth mechanism. In what follows, we refer to (1.1) as the complete Gatenby–Gawlinski model to distinguish it from a corresponding reduced version to be introduced later.

The basics of such a modeling is rather clear. Firstly, healthy cells – denoted by U – have a certain reproduction level (supposed logistic with rate 1, for simplicity) and are deteriorated by the acid, following the standard mass action law with kinetic constant d. Secondly, the tumor cells V have the capability of spreading at a rate that depends on the quantity of healthy cells U, and they also reproduce according to a logistic law, with a different rate denoted by r. A rough justification of the dependence on U in the coefficient of xV is that tumor cells – possessing a high-degree of invasiveness – can hardly move when the density of healthy cells is high. Specifically, the coefficient is null (no motion of cancerous tissue) if the healthy cells are at carrying capacity. Finally the concentration W of lactic acid undergoes diffusion at constant rate a, and is increased proportionally to the unknown V, with kinetic constant b, until it reaches the saturation level W=V. Let us stress that the third unknown W has no direct effect on the dynamics of the variable V. The system is meaningful when the parameters a,b,d,r are all positive. Modifications of the original model have also been considered by many authors: among others, we quote here [5] (generalized Gatenby–Gawlinski model), [6] (linear diffusion in the tumor variable), [7] (effect of chemotherapy), [8] (stromal interaction), [9] (distinction between intracellular and extracellular proton dynamics). An interesting feature of the Gatenby–Gawlinski model (1.1) is the numerical evidence of existence of invasive propagation fronts, i.e. special solutions describing the invasion of the cancer cells into the healthy tissue. To our knowledge, no rigorous proof of the existence of such front is available so far, except in some limiting parameter regimes.

To decrease the complexity of (1.1), we consider the reduced system tU=U{f(U)dV},tV=x{f(U)xV}+rVf(V),obtainable as a formal limit in the regime aconst. and b, that is replacing the last dynamical equation with the trivial constitutive identity W=V (see [10], [11], [12]). Such a model intends to describe the case in which the tumor cells act directly on the healthy tissue with no additional specific intermediate (at the level of ODE, an analogous system has been discussed in [13]). We refer to (1.2) as the reduced Gatenby–Gawlinski model. It can be regarded as a simplified version of the system proposed in [5], obtained by rescaling the space and choosing a2=d,r2=r,a1=d1=d2=0,D=1,c=b,where the parameters a1,a2,d1,d2,r2 and D are as in [5]. Incidentally, let us observe that there are other possible reductions of the same original complete model, which are obtained by considering different parameter regimes and could also be worth investigating.

The diffusion coefficient f(U)=1U in the V-equation of system (1.2) can be seen as an inherent defense process exerted by the healthy tissue in response to the presence of the tumor. If the healthy cells are at their carrying capacity, normalized to 1, no invasion can occur; in contrast, tumor starts growing whenever the density of healthy cells is lower than 1, and in absence of healthy tissue the tumor is free to permeate all the space. Incidentally, we observe that f(U) takes negative values when U>1, so that the diffusion equation for V becomes ill-posed in such a regime. In the subsequent traveling wave analysis, we restrict our attention to values U(0,1), so that f(U)>0.

Propagation fronts attracted the interest of many researchers because they provide the simplest mathematical framework describing the process of biological invasion. Rigorous results concerning existence and asymptotic stability of traveling wave solutions for scalar equations of the form tU=x{ϕ(U)xU}+f(U),with reaction term f and nonlinear diffusion ϕ, have been obtained under various assumptions. When ϕ is a strictly positive constant, the diffusion is linear and Eq. (1.4) is semilinear. For ϕ dependent on U and attaining strictly positive values, the diffusion is nonlinear and non-degenerate. Here, we are interested in the case in which the function ϕ is non-negative and null at some specific points, usually U=0, the archetypal example being the porous medium equation for which ϕ(s)sp for some p>0. Let us stress that considering a function ψ such that ψ=ϕ, Eq. (1.4) can be rewritten as tU=x2ψ(U)+f(U),corresponding to a porous medium equation with reaction. Existence of traveling wave solutions has been widely explored in such a context, see [14], [15], [16], [17], [18], [19], [20], [21], [22] for a single degeneration point, and [17], [23], [24], [25] for multiple degenerations. Different forms of the multiplier function ϕ, e.g. depending on xU, have also been considered (see [26], [27]).

In the scalar case, the typical existence statement – valid for linear, nonlinear non-degenerate or degenerate diffusion – can be rephrased as follows. Let the function f be of logistic type, i.e. it has two zeros (say 0 and 1) and is positive in between. Then the scalar equation (1.4) supports traveling waves U(x,t)=U(xct) satisfying the asymptotic conditions U()=1 and U(+)=0 if and only if cc, for some strictly positive c. Moreover, the solution to the initial value problem with Heaviside-like data, characterized by a sharp jump from 0 to 1, converges in an appropriate sense to the traveling wave connecting the states 0 and 1 and moving at critical speed c.

For traveling waves of linear diffusion equations, stability analysis is a classical subject dating back to the pioneering papers by Fisher [28] and Kolmogorov, Petrovskii and Piscounov [29]. Much less is known, however, in the case of degenerate diffusion equations. One of the few results is contained in [30], following the general method outlined in [31] for general reaction–diffusion systems.

As expected, for systems, the situation is less clear. First of all, it is necessary to agree on the terminology. Let us consider, for simplicity, a 2 × 2 reaction–diffusion system of the form tu=x{ϕ11(u,v)xu+ϕ12(u,v)xv}+F(u,v),tv=x{ϕ21(u,v)xu+ϕ22(u,v)xv}+G(u,v),for some diffusivities ϕij with i,j{1,2} and reaction terms F,G.

Many examples for (1.5) with a significant applied perspective can be provided. Among others, the celebrated Keller–Segel chemotaxis model –proposed as a description for the motion of bacteria towards some optimal environment [32]– fits into the class choosing ϕ11(u,v)=a,ϕ12(u,v)=0,ϕ21(u,v)=bχ(u),ϕ22(u,v)=μ(u),and F(u,v)=κ(u)v, G(u,v)=0 for some parameters a,b>0 and functions κ, μ, χ. The Keller–Segel system can be regarded as a prototype of “exotaxis” models because the gradient in the concentration of one species induces a flux of another species.

A second example is [33] – the progenitor of a long lineage – where the terminology cross-diffusion system has been used to denote a particular case of (1.5) characterized by the presence of the diffusion term x2ψi in place of x(ϕi1xu+ϕi2xv) for i=1,2. In other words we assume that there exist ψ1,ψ2 such that ϕ11=uψ1, ϕ12=vψ1, ϕ21=uψ2, ϕ22=vψ2, which is not the case in general.

A simplified version of (1.5) is obtained by assuming the terms ϕ12 and ϕ21 to be null, that is focusing on systems with cross-dependent self-diffusivities: tu=xϕ11(u,v)xu+F(u,v),tv=xϕ22(u,v)xv+G(u,v).Both species are submitted to self-diffusion with a diffusivity coefficient that, in general, may depend on the other variable. The reduced Gatenby–Gawlinski model (1.2) fits into (1.6) with the choices ϕ11(u,v)=0,ϕ22(u,v)=f(u),F(u,v)=u{f(u)dv},G(u,v)=rvf(v).Coming back to the topic of invasion fronts in the case of reaction–diffusion systems, a huge difference with respect to the scalar case arises already for linear self-diffusion because the dimension of the phase-space for the traveling wave ODE is strictly larger than two. Nonlinear non-degenerate self-diffusions make the analysis harder, but, in principle, still manageable with an adapted strategy. In contrast, when the diffusion operator degenerates at some values, the situation becomes more involved with a pivotal role played by an appropriate desingularization procedure, which will be explained below.

Focusing on the case of cross-dependent self-diffusivities, after the pioneering contribution of Aronson [14] devoted to a predator–prey system, the attention moved towards the model proposed by Kawasaki et al. [34], which attempts to provide a detailed description of the patterns generate by some colonies of bacteria, called Bacillus subtilis (see [35] for a comprehensive review on cooperative self-organization of micro-organisms). This model is composed by two coupled evolution equations for the population density b and the concentration of nutrient n. Degenerate cross-dependent self-diffusion appears in the equation for the bacteria b and is proportional to the product of the two unknowns, i.e. D(n,b)nb. Investigations on existence of propagating fronts in bacteria growth models – either from a purely analytical point of view or from a numerical perspective – have been performed in [36], [37], [38], [39]. In particular, in [38], [39], the existence result is very similar to the one valid for the scalar case, including the existence of a traveling wave for the critical speed.

Existence of propagation fronts for the complete Gatenby–Gawlinski model (and its modifications) is cogently supported by partial results and numerical calculations, see [5], [10], [12], [40], [41]. However, rigorous mathematical results are very limited if not completely missing. Again, a distinguished feature of the model is the presence of cross-dependent self-diffusion. In contrast with the bacteria models, the cross-diffusion term f(U) in the V-equation is a monotone decreasing function of the variable U, and our goal is to explore in detail the consequences of such kind of coupling, in the particular example of the reduced Gatenby–Gawlinski model (1.2).

A mathematical model, aiming to describe how a population of melanoma cells invades into human skin and possessing a structure similar to the one explored here, has been recently proposed in [42] with a corresponding propagation front analysis considered in [43] (see also [44]).

A traveling wave for (1.2) is a solution of the form U(x,t)=U(xct),V(x,t)=V(xct),where the parameter cR is the propagation speed. If ξxct denotes the space variable in a comoving frame, the traveling wave profile formally satisfies the ODE system cdUdξ+U{f(U)dV}=0,ddξf(U)dVdξ+cdVdξ+rVf(V)=0.A propagation front is a special type of traveling wave, enjoying the asymptotic conditions limξ±(U,V)(ξ)=(U±,V±).The asymptotic values (U±,V±) are forced to be constant equilibria to system (1.2). While numerical evidence of existence of traveling waves has been provided in [11], [12], no rigorous result was obtained so far. In what follows, we will reconsider the above definition of propagation front in order to incorporate the presence of possible degeneracies.

When d1, system (1.2) has exactly four constant equilibria: the trivial state (Ū,V̄)=(0,0), the healthy state (1,0), the cancerous state (0,1) and the heterogeneous state (1d,1). Note that the last equilibrium is positive, hence biologically significant, if and only if d<1. In the limiting case d=1, the system has only three uniform equilibria.

In this paper, we mainly concentrate on the case d>1, which seems most relevant in cancerology (considering d as a measure of aggressiveness of the tumor), and we look for traveling wave solutions that describe the invasion of the healthy state by the infected state. In other words, we choose as asymptotic values (U,V)=(0,1)and(U+,V+)=(1,0).This situation is referred to as homogeneous invasion. We also study more succinctly the regime 0<d<1, in which homogeneous invasion is not possible. In that case, we focus on heterogeneous invasion which corresponds to the asymptotic states (U,V)=(1d,1)and(U+,V+)=(1,0).The particular case d=1 is non-generic, and will not be considered here.

It is important to keep in mind that the second equation in (1.2) is a degenerate parabolic equation, see [21], in the sense that the coefficient f(U) in front of the leading order term x2V vanishes when U=1. For that reason, it is not clear at all that system (1.2) has global classical solutions, and a similar caveat applies to the ODEs (1.7) satisfied by the traveling waves. Hence, we adopt here the following definition, which is adapted from [45] and based on the notion of weak solution.

Definition 1.1

The triple (U,V;c) is a propagation front for system (1.2) connecting the asymptotic states (U,V) and (U+,V+) if

  • (i)

    (U,V)C(R;[0,1])×C(R;[0,1]) and f(U)dVdξL2(R);

  • (ii)

    (U,V) is a weak solution to (1.7), i.e. for all (ϕ,ψ)C1(R)×C1(R) with compact support RUcdϕdξ[f(U)dV]ϕdξ=0,Rf(U)dVdξ+cVdψdξrVf(V)ψdξ=0;

  • (iii)

    the asymptotic conditions (1.8) are satisfied.

The couple (U,V) is the profile of the front and the value c is the speed of propagation.

The main result of this paper is the following.

Theorem 1.2

Assume that d>1 and r>0. For any c>0, the reduced Gatenby–Gawlinski system (1.2) has a propagation front (U,V;c) connecting (0,1) with (1,0). This solution is unique up to translations and both components U,V are strictly monotone functions of ξ=xct (see Fig. 1).

The proof also provides detailed information on the behavior of the front profile as ξ±. In particular, we can choose a translate of the wave such that U(ξ)=αeμξ+O(e(μ+η)ξ),V(ξ)=1eλξ+O(e(λ+η)ξ),asξ,for some α>0, where λ=12(c+c2+4r)>0,μ=d1c>0,η=min(λ,μ).Moreover, there exists β>0 such that U(ξ)=1βeγξ+O(e2γξ),V(ξ)=βr+1deγξ+O(e2γξ),as ξ+, where γ=r/c.

Most remarkably, Theorem 1.2 shows that there exists no minimal speed for the propagation fronts of system (1.2). This is in sharp contrast with what happens for scalar equations involving a degenerate diffusion, see [19], [20] and also [39] (for the Kawasaki system). To elaborate on that, we consider a further formal reduction of system (1.2), which seems reasonable at least for solutions that evolve slowly in time. In view of the first equation in (1.2), one can expect that the first component U should stay close to 1dV if dV1 and to zero if dV>1. Assuming this to be exactly true, we obtain a scalar evolution equation for the second component V  : tV=xϕ(V)xV+rVf(V),whereϕ(V)min(dV,1).Strictly speaking, the results of [20], [24] do not apply to (1.16) because the diffusion coefficient ϕ is only a Lipschitz function of V. Disregarding that technical issue, we expect nevertheless that Eq. (1.16) has monotone front solutions satisfying V()=1, V(+)=0 if and only if cc, for some minimal speed c=c(d,r)>0. Moreover, when c=c, the front profile is “sharp” in the sense that there exists ξ̄R such that V(ξ)=0 for all ξξ̄. Quite surprisingly, Theorem 1.2 shows that the PDE system (1.2) behaves differently: propagation fronts exist for all positive speeds c>0, no matter how small, and all front profiles are smooth and strictly monotone. Sharp fronts, which are typical for scalar equations with degenerate diffusion, do not exist in system (1.2).

Of course, this discrepancy means that the formal reduction leading to (1.16) is not justified. In fact, if (U,V) is the front profile given by Theorem 1.2, for some values of the parameters d,r,c, we can introduce the effective diffusion coefficient ϕ:(0,1)(0,1) defined by 1U(ξ)=ϕ(V(ξ)),ξR.By construction, the second component V of the front profile is a traveling wave solution of the scalar equation (1.16) with ϕ given by (1.17). In particular we must have cc(ϕ,r), where c is the minimal speed for the scalar equation. Numerical calculations shows that the shape of the effective diffusion coefficient ϕ depends strongly on the values of the parameters d,r,c, and is often very different from the naive guess ϕ(V)=min(dV,1). This is especially true when c is small, in which case ϕ(V) is found to be extremely flat near the origin V=0, see the discussion at the end of Section 3.

Going back to (1.2), in the case where 0<d<1, the system has an additional, biologically significant, equilibrium (Ū,V̄)=(1d,1) in which healthy and cancerous cells coexist. In that case, it is natural to consider traveling waves that connect the coexistence state to the healthy state given by (1.10). We have the following analogue of Theorem 1.2.

Theorem 1.3

Assume that 0<d<1 and r>0. For any c>0, the reduced Gatenby–Gawlinski system (1.2) has a propagation front (U,V;c) connecting (1d,1) with (1,0). This solution is unique up to translations and both components U and V are strictly monotone. Moreover, there is no propagation front connecting (0,1) with (1,0) in that case (see Fig. 2).

The rest of this article is organized as follows. Section 2 is mainly devoted to the proof of Theorem 1.2. Our starting point is the desingularization of the ODE system (1.7), using a standard procedure that was already known for scalar equations with degenerate diffusion. Propagation fronts are then constructed as heteroclinic connections between two equilibria of the desingularized system. The unstable manifold of the infected state is two-dimensional in that setting, which forces us to introduce an additional shooting parameter, and an important part of our analysis relies on monotonicity properties with respect to that shooting parameter. The proof of Theorem 1.3 goes along the same lines, and is briefly presented in Section 2.8. In Section 3, we explore the limiting regime where c0, and we derive an asymptotic expansion of the front profile that is remarkably accurate even at moderately small speeds. Finally, we draw some conclusions in Section 4 and we outline a few perspectives.

Section snippets

Existence of propagation fronts

This section is devoted to the proofs of Theorem 1.2, Theorem 1.3. We assume that the reader is familiar with center manifold theory for ODEs and normal form theory. All necessary material can be found in classical monographs such as [46], [47], [48], to which we shall refer when needed.

In what follows, we fix the parameters d>0 and r>0 in system (1.7). All quantities that appear in the proof depend on d,r, but for notational simplicity this dependence will not be indicated explicitly.

Asymptotic analysis of slowly propagating fronts

Perhaps the most striking aspect of Theorem 1.2, Theorem 1.3 is the absence of a minimal speed for the monotone traveling waves of system (1.2). To understand what happens in the singular limit c0, we compute in this section the leading term of a (formal) asymptotic expansion of the front profile. We do not feel the necessity of rigorous proofs at this stage, but we provide numerical illustrations supporting our arguments. We always assume that d>1, and we consider propagation fronts

Conclusions and perspectives

We conclude this paper with a list of possible questions that are, in our opinion, worth investigating in the future.

Acknowledgments

The research of the first named author is partially supported by the grant ISDEEC ANR-16-CE40-0013 of the French Ministry of Higher Education, Research and Innovation . This work benefited from mutual invitations, on the occasion of which the hospitality of Université Grenoble Alpes and Sapienza, Università di Roma is gratefully acknowledged.

References (54)

  • KawasakiK. et al.

    Modeling spatio-temporal patterns generated by Bacillus subtilis

    J. Theoret. Biol.

    (1997)
  • HolderA. et al.

    A model for acid-mediated tumour growth with nonlinear acid production term

    Appl. Math. Comput.

    (2014)
  • HilhorstD. et al.

    Interface dynamics of the Fisher equation with degenerate diffusion

    J. Differential Equations

    (2008)
  • FenichelN.

    Geometric singular perturbation theory for ordinary differential equations

    J. Differential Equations

    (1979)
  • LewisM. et al.
  • WarburgO. et al.

    The metabolism of tumors in the body

    J Gen. Phyisiol.

    (1927)
  • AlfaroukK. et al.

    Glycolysis, tumor metabolism, cancer growth and dissemination. A new pH-based etiopathogenic perspective and therapeutic approach to an old cancer question

    Oncoscience

    (2014)
  • GatenbyR. et al.

    A reaction-diffusion model of cancer invasion

    Cancer Res.

    (1996)
  • McGillenJ. et al.

    A general reaction-diffusion model of acidity in cancer invasion

    J. Math. Biol.

    (2014)
  • de AraujoA. et al.

    An analysis of a mathematical model describing acid-mediated tumor invasion

    Math. Methods Appl. Sci.

    (2019)
  • StinnerC. et al.

    A multiscale model for pH-tactic invasion with time-varying carrying capacities

    IMA J. Appl. Math.

    (2015)
  • MasciaC. et al.

    Numerical investigation of some reductions for the gatenby–gawlinski model

    (2021)
  • MoschettaP. et al.

    Numerical investigation of the Gatenby-Gawlinski model for acid-mediated tumour invasion

    Rend. Mat. Appl.

    (2019)
  • EnglerH.

    Relations between travelling wave solutions of quasilinear parabolic equations

    Proc. Amer. Math. Soc.

    (1985)
  • GildingB. et al.

    A Fisher/KPP-type equation with density-dependent diffusion and convection: traveling-wave solutions

    J. Phys. A: Math. Gen.

    (2005)
  • Sánchez-GardunoF. et al.

    Existence and uniqueness of a sharp travelling wave in degenerate non-linear diffusion Fisher–KPP equations

    J. Math. Biol.

    (1994)
  • Sánchez-GardunoF. et al.

    Travelling wave phenomena in non-linear diffusion degenerate Nagumo equations

    J. Math. Biol.

    (1997)
  • Cited by (0)

    View full text