Nonlinear stability of planar traveling wavesin a chemotaxis model of tumor angiogenesis with chemical diffusion

# Nonlinear stability of planar traveling waves in a chemotaxis model of tumor angiogenesis with chemical diffusion

Myeongju Chae Department of Mathematics, Hankyong University, Anseong-si, Gyeonggi-do, Republic of Korea  and  Kyudong Choi Department of Mathematical Sciences, Ulsan National Institute of Science and Technology, Ulsan, Republic of Korea
July 1, 2019
###### Abstract.

We consider a simplified chemotaxis model of tumor angiogenesis, described by a Keller-Segel system on the two dimensional infinite cylindrical domain , where is the circle of perimeter . The domain models a virtual channel where newly generated blood vessels toward the vascular endothelial growth factor will be located. The system is known to allow planar traveling wave solutions of an invading type. In this paper, we establish the nonlinear stability of these traveling invading waves when chemical diffusion is present if is sufficiently small. The same result for the corresponding system in one-dimension was obtained by Li-Li-Wang (2014) [16]. Our result solves the problem remained open in [3] at which only linear stability of the waves was obtained under certain artificial assumption.

###### Key words and phrases:
tumor, angiogenesis, chemotaxis, Keller-Segel, nonlinear stability, 2D infinite cylinder, planar traveling wave, infinite strip, chemical diffusion, Cole-Hopf transform
92B05, 35K45

## 1. Introduction

### 1.1. A Keller-Segel system

The formation of new blood vessels from pre-existing vessels, which is so-called angiogenesis, is the essential mechanism for tumour progression and metastasis. Focusing on the interaction between endothelial cells and growth factor, a simplified model of tumor angiogenesis can be described by the following Keller-Segel system [7, 14, 22]:

 (1.1) ∂tn−Δn=−∇⋅(nχ(c)∇c)∂tc−ϵΔc=−cmn.

We consider the above system in two-dimension with a front boundary condition in and a periodic condition in , both specified later, with and . In a general Keller-Segel context, the unknown is the bacterial density while the unknown is the concentration of chemical nutrient consumed by bacteria at position and time . Considering formation of new blood vessels, denotes the density of endothelial cells while does the concentration of the protein known as the vascular endothelial growth factor (VEGF). The chemosensitivity function is a given decreasing function, reflecting that the chemosensitivity gets lower as the concentration of the chemical gets higher. The positive constant is the diffusion rate constant for the chemical substance while indicates the consumption rate of nutrient .

When we model endothelial angiogenesis, we interpret that the endothelial cells behave as an invasive species, responding to signals produced by the hypoxic tissue. Accordingly, we choose the -axis by the propagating direction and the system (1.1) is given the front condition at left-right ends such that

 (1.2) limx→−∞n(x,y,t)=n−>0,limx→∞n(x,y,t)=0, (1.3) limx→−∞c(x,y,t)=0,limx→∞c(x,y,t)=c+>0.

To all functions in this paper, we impose the periodic condition in -variable of period .

A planar traveling wave solution of (1.1) is a traveling wave solution independent of the transversal direction :

 (1.4) n(x,y,t)=N(x−st),c(x,y,t)=C(x−st)

with a given wave speed which we always assume positive in this paper without loss of generality. We consider only waves satisfying the boundary conditions (1.2) and (1.3) which means

 (1.5) limx→−∞N(z)=n−>0,limx→+∞C(z)=c+>0,limx→+∞N(z)=limx→−∞C(z)=0.

We also assume that

 (1.6) limx→±∞N′(z)=limx→±∞C′(z)=0.

To have a traveling wave, it is known that that the chemosensitivity function needs to be singular near (e.g. see [13, 25]). In the paper [13], , which yields the logarithmic singularity , is assumed, which choice of is then adopted on modeling the formation of the vascular network toward cancerous cells (e.g. see [7, 14, 22]). The existence of traveling wave solution with an invading front might be an evidence of the tumor encapsulation (e.g. see [1, 2, 26]).

In this paper, we consider only the case and of (1.1):

 (1.7) ∂tn−Δn=−∇⋅(n∇cc),∂tc−ϵΔc=−cn,(x,y,t)∈R×Sλ×R+

where is the circle of perimeter This 2D cylindrical domain would be understood as a virtual channel where newly generated blood vessels toward the chemical (VEGF) will be located. We focus on establishing the time asymptotic stability of a planar traveling wave solution of (1.7). The restriction on is required for treating the singularity of by the Cole-Hopf transformation

 (1.8) p:=−∇lnc=−∇cc=−(∂xcc,∂ycc).

A well-written explanation of the system including biological interpretation can be found in [23, 24] (also refer to [20] and the references therein).

### 1.2. A parabolic system of conservation laws

By the Cole-Hopf transform, we translate the singular Keller-Segel system (1.7) into the following system of without singularity:

 (1.9) ∂tn−Δn=∇⋅(np),∂tp−ϵΔp=−2ϵ(p⋅∇)p+∇n,(x,y,t)∈R×Sλ×R+

with the notation .

By denoting

 (1.10) P:=−C′/Cand P:=(P,0),

we have a planar traveling wave solution of (1.9) of speed with the boundary conditions inherited from those of . The existence and some properties of those waves and can be found in [30, 17]. We put some of the results on the waves we need in Subsection 2.1.

The study on the existence of traveling wave solutions of a Keller-Segel model was initiated by the paper [13] then many works followed (see [9] and the references therein). We also refer to the survey paper [31] which is an excellent exposition of the topic. The existence of traveling waves with the front conditions (1.2) and (1.3) can be found in [32] for , and [17, 30] for . When considering the one dimensional system (i.e. no -dependency in (1.7)), the nonlinear stability results were shown in a weighted Sobolev space in [11] for , and [16] when is small (also see [20]). The weighted Sobolev space has commonly appeared when studying nonlinear stability of viscous shocks of conservations laws since [12] (also see [19]).

The study of higher dimensional traveling waves is a very interesting topic and remains open for many questions including existence and stability of such waves as indicated in [31]. As a special case in 2D, planar waves for an infinite cylinder was considered by [3] following the spirit of the nonlinear energy estimates developed in [11] for the whole line case. In angiogenesis, one may consider that a blood vessel in our body has a 2D cylinder structure.

The previous result [3] mainly proved two things: one is the nonlinear stability for and the other is the stability of the linearized equation for small under the additional mean-zero assumption in transversal direction for some technical reason. In addition to these two results, Theorem 1.6 in [3] gives a hint why studying planar waves is natural instead of doing general 2D traveling waves by showing that the -derivative of any smooth solution decays to zero in -sense under certain additional assumption.

In this paper, we show that traveling wave solutions of the nonlinear system (1.9) are globally stable under the smallness assumption on the parameters and without the artificial mean-zero assumption in transversal direction , which was needed in [3] even for the corresponding linearized system. Indeed, the main estimate (2.10) holds uniformly for small when the antiderivative of a perturbation of the form is sufficiently small in a weighted Sobolev space (see (2.8) and (2.9)). Our result can be considered as an extension of [3] into case and an extension of [16] into 2D case. See Theorem 2.9 and Subsection 2.2 for the precise set-up. We state the stability result in terms of the perturbation of in Theorem 2.9, then explain the implication of the theorem for the perturbation of to in Remark 2.10.

At first glance, the transformed -system (1.9) seems simpler than the -system (1.7) to analyze since this parabolic system (1.9) of conservation laws does not have the logarithmic singularity. As a price for this, however, the quadratic nonlinear term appears, and it is not clear at all if the linear term in the main perturbation equation (2.13) produced by the nonlinear term in (1.9) can be controlled by the diffusion term in (2.13) produced by the diffusion term in (1.9).

In this regard, the main obstacle is to handle the quantity in (3.3), which is the time integral of a localized -norm of multiplied by . We overcome the difficulty thanks to certain dissipations of a localized -norm of (not of ) together with a careful manipulation done in Lemma 3.3 (see Remark 3.2). In doing so, we need the smallness assumption on . This idea was first used in [16] for the one-dimensional system while for our two dimensional system, it becomes more delicate due to the non-symmetric nature of the main perturbation equation (2.13) on the propagating direction and the transversal direction . For instance, when we denote , we see the non-symmetric term in (3.5). The smallness condition on the chemical diffusion constant might be understood in the sense that the chemical in angiogenesis often diffuses in the dense network of extracellular matrix and tissues which are almost static as mentioned in [31].

Unfortunately, we also need the smallness assumption on the perimeter of a 2D infinite cylinder, and it appears due to a technical reason in our proof. In fact, with wave speed , we ask the product to be smaller than a given absolute constant (see (3.2)). This condition enables us to employ Poincaré inequality (3.1) in the transversal direction in order to control a non-localized -norm of (see (3.5) and (3.6))). In our opinion, it is very challenging to remove this technical smallness assumption on .

For the Cauchy problem of (1.1), we refer to [4, 5, 6, 7, 15], where [4, 5] prove the existence of a global weak solution, and [15] proves the existence of a global classical solution considering the zero chemical diffusion case in a multi-dimension. When a bounded domain is considered, a boundary layer may appear. We refer to [10] and [21] for the stability questions of the layer.

The remaining parts of the paper are organized as follows. In Subsection 2.1, we introduce background materials including the existence and some properties of traveling wave solutions and state the main result (Theorem 2.9) with its set-up in Subsection 2.2. In Subsection 2.3, we state the local existence of a perturbative solution and its a priori uniform-in-time estimate. In Section 3, we prove the uniform-in-time estimate. The zero-th and first order estimates (Subsection 3.1 and 3.2) are the essential part. Then, the higher order estimate (Lemma 3.9) can be obtained in a similar way. We present its proof for completeness in Subsection 3.3.

## 2. Main theorem and background materials

### 2.1. Existence and Properties of traveling wave solutions

We collect some results on traveling wave solutions and with the front conditions introduced in Section 1.

We first observe that a traveling wave solution defined by (1.4) solve the following ODE system by plugging the expression (1.4) into (1.7):

 (2.1) −sN′−N′′=−(C′CN)′,−sC′−ϵC′′=−CN.
###### Theorem 2.1 ([30] Lemma 3.2, Lemma 3.4).

A monotone solution of (2.1) for with the boundary conditions (1.5) and (1.6) exists if the relation

 n−=(1+ϵ)s2

holds. More precisely it holds that

1. and ,

2. , as , and

3. and .

In [30], the author used the results of the KPP-Fisher equation to establish the above theorem.

The relation gives the system

 (2.2) −sN′−N′′=(PN)′,−sP′−ϵP′′=(N−ϵP2)′.

We observe that is a traveling wave solution of (1.9). From (1.5) with the above theorem, the wave is given the boundary condition

 (2.3) N(−∞)=(1+ϵ)s2,N(+∞)=0,P(−∞)=−s,P(+∞)=0,N′(±∞)=P′(±∞)=0.

We abbreviate by for any function on . Moreover, the following theorem holds.

###### Theorem 2.2 ([17] Theorem 2.1, or [16]).

For , if is small, then there exists a monotone solution to (2.2) with the boundary condition (2.3). In particular, it satisfies with and with , and it is unique up to a translation.

The next lemma gives a uniform estimate of and for any small .

###### Lemma 2.3.

For , there exist constants and such that if is a solution of (2.2) and (2.3) for some given by Theorem 2.2, then

 |N(k)|≤L,|P(k)|≤L, for 0≤k≤2, and
###### Proof.

The estimate for in the first line was proved in [16] while the proof of the rest can be found in Lemma in [3]. ∎

Lastly, we need the following lemma which gives a point in contained both in the transition layer of and in that of .

###### Lemma 2.4.

For any , there exists a constant such that if is a solution of (2.2) and (2.3) for some given by Theorem 2.2, then there exists a point satisfying

 P(z0)=−s2 and N(z0)≥s24.
###### Proof.

Since is continuous on and , there exists a point such that . To show for sufficiently small , recall the equation (2.2). From , we have

 −sP−ϵP′=(N−ϵP2)

Assume that is smaller than in Lemma 2.3. Then for any , we get

where is the constant from Lemma 2.3. We take small enough to have . Then for any . ∎

###### Remark 2.5.

The lemma is due to the fact which means the transition layers of is overlapped with that of in some extent when is small enough.

###### Remark 2.6.

The first equation in (2.2) with (1.6) and (2.3) gives the simple relation between and :

 (2.4) −N′N=s+P.
###### Remark 2.7.

If we denote , then the above lemma implies

 (2.5) w′(z)w(z)≥s2%forz≥z0andw(z)≤4s2≤16s4Nfor z≤z0.

Indeed, for , we have

 w′w=(1/N)′1/N=−N′/N21/N=−N′N=s+P≥s+P(z0)=s2

thanks to (2.4) and . For , we have

 w=1N≤1N(z0)≤4s2=16s4⋅s24≤16s4N(z0)≤16s4N(z)

due to . We will use (2.5) in the proof of Lemma 3.7.

Figure describes the above discussions including monotonicity of waves.

### 2.2. Main theorem

We recall (1.9):

 (2.6) ∂tn−Δn=∇⋅(np)∂tp−ϵΔp=−2ϵ(p⋅∇)p+∇n,(x,y,t)∈R×Sλ×R+.

Let be a traveling wave solution of (2.6) with (2.3). In the below we introduce a weighted Sobolev space where our perturbative functions are constructed. We use the weight function (only in the horizontal direction) defined by

 w(z)=1N(z),z∈R

where this unbounded weight was essentially introduced in [16] to handle the difficulty coming from the vacuum state . Note that is monotonically increasing, and when by (2.3) and Theorem .
For an integer and for any , we define the Sobolev spaces and a weighted Sobolev space for functions periodic in with period as follows;

 Hk :={φ∈Hkloc(R2)|∑i+j≤k∑n∈Z∫Rn2j|∂izφn(z)|2dz<∞,φ(⋅z,⋅y+λ)=φ(⋅z,⋅y)},
 Hkw :={φ∈Hkloc(R2))|∑i+j≤k∑n∈Z∫Rn2j|∂izφn(z)|2w(z)dz<∞,φ(⋅z,⋅y+λ)=φ(⋅z,⋅y)}

where for each and for each , is the th Fourier coefficient of the (-) periodic (in ) function .
We define the norms by 111 The two quantities used to define are equivalent up to the transversal length scale . In this paper, we do not pursue any estimate which needs to hold uniformly on .

 ∥φ∥2Hk:=∑i+j≤k∫R×[0,λ]|∂iz∂jyφ(z,y)|2dzdy, ∥φ∥2Hkw:=∑i+j≤k∫R×[0,λ]|∂iz∂jyφ(z,y)|2w(z)dzdy.

Note that for any , we know

 (2.7) ∥f∥2Hk≤(1+ϵ)s2∥f∥2Hkw

due to .

We perturb the equation (2.6) around the wave

 (2.8) n(x,y,t)=N(x−st)+∇⋅φ(x−st,y,t)and p(x,y,t)=P(x−st)+∇ψ(x−st,y,t).

With in the moving frame, we expect that for each time , the perturbation lies on the following function class:

 (2.9) φ=(φ1,φ2)∈(H3w)2 and ψ∈H3 with ∇ψ∈H2w.
###### Remark 2.8.

Such a one-sided decaying function (in the weighted Sobolev space) appears typically with respect to the solvability of in the infinite cylinder (e.g. see [28]). An explanation relevant to the perturbation (2.8) is given in Remark 1.2 in [3].

Now we state the main theorem:

###### Theorem 2.9.

For any and for any such that the product is sufficiently small, there exist constants , , and such that if is a solution of (2.2) for some with (2.3) given by Theorem 2.2, then for any initial data of (2.6) in the form of and satisfying

 M0:=(∥φ0∥2H3w+∥ψ0∥2H3+∥∇ψ0∥2H2w)≤K0,

there exists a unique global solution of (2.6) in the form of

 n(x,y,t)=N(x−st)+∇⋅φ(x−st,y,t),p(x,y,t)=P(x−st)+∇ψ(x−st,y,t),

where and , and satisfies the following inequality:

 (2.10) supt∈[0,∞)(∥φ∥2H3w+∥ψ∥2H3+∥∇ψ∥2H2w)(t)+∫∞0(∥∇φ∥2H3w+∥∇ψ∥2H2w+ϵ∥∇4ψ∥2w)(t)dt≤C0M0.
###### Remark 2.10.

From (2.8), (1.10) and (1.8), we have . Together with , the above theorem implies

 supt∈[0,∞)(∥n(⋅z+st,⋅y,t)−N(⋅z)∥2H2w+∥∇(logc(⋅z+st,⋅y,t)−logC(⋅z))∥2H2w)

Before closing this subsection we give a summary on notations used in the paper.

 Ω=R×[0,λ],
 w(z):=1N(z),
 M(t):=sups∈[0,t](∥φ(s)∥2H3w+∥ψ(s)∥2H3+∥∇ψ(s)∥2H2w),
 M0:=(∥φ0∥2H3w+∥ψ0∥2H3+∥∇ψ0∥2H2w)
 ∥f∥:=∥f∥L2(Ω),
 ∥f∥2k:=∥f∥2Hk=k∑|α|=0∫Ω|Dαf|2dzdy,
 ∥f∥2k,w:=∥f∥2Hkw=k∑|α|=0∫Ω|Dαf(z,y)|2w(z)dzdy,
 ∫f:=∫Ωf(z,y)dzdy,
 ∫t0g:=∫t0g(σ)dσ.

Here we use the notation to indicate certain norm in space only. For instance, when is time-dependent then means in the sequel.

### 2.3. Perturbation equation

In this subsection, we derive the system on first. Next we state the main propositions including results on the local existence and the uniform estimates of .

From (2.6) and (2.8), by setting

 (2.11) u=∇⋅φ and v=∇ψ

temporarily in (2.8), we obtain

 (2.12) ut−suz−Δu=∇⋅(Np+Pu+up),vt−svz−ϵΔv=−2ϵ(((P+v)⋅∇)(P+v)−(P⋅∇)P)+∇u.

Plugging the relation (2.11) in (2.12) and taking off derivatives, we find that the antiderivative of satisfies the system

 (2.13) φt−sφz−Δφ=N∇ψ+P∇⋅φ+∇⋅φ∇ψ,ψt−sψz−ϵΔψ=−2ϵP⋅∇ψ−ϵ|∇ψ|2+∇⋅φ

for . In doing so, we use the curl free property of and . Here the term means the vector . The multidimensional setting (2.11) was proposed in [3]. Looking for a perturbation in our system as an antiderivative follows the setting in one dimensional works [11], [16]. This method can be found in the study on the nonlinear stability of shock profiles of viscous conservation laws under the mean zero condition with a weight function since the papers [12] and [8]. Without the mean zero condition, we refer to [18], [29], [27] and references therein.

We obtain Theorem 2.9 immediately without any difficulty once we prove the proposition below.

###### Proposition 2.11.

For any and any such that the product is sufficiently small, there exist constants , and such that if is a solution of (2.2) for some with (2.3) given by Theorem 2.2, then we have the following:

For any initial data of (2.13) satisfying

 M0:=∥φ0∥2H3w+∥ψ0∥2H3+∥∇ψ0∥2H2w≤K0,

there exists a unique global solution of (2.13) where and , and satisfies the following inequality:

 supt∈[0,∞)(∥φ∥2H3w+∥ψ∥2H3+∥∇ψ∥2H2w)+∫∞0(∥∇φ∥2H3w+∥∇ψ∥2H2w+ϵ∥∇4ψ∥2w)dt≤C0M0.

Proposition 2.11 is a consequence of the following two propositions: Proposition 2.12 which gives a local-in-time existence result and Proposition 2.13 which shows an a priori uniform-in-time estimate.

###### Proposition 2.12.

Let and . For sufficiently small , if is a solution of (2.2) with (2.3) given by Theorem 2.2, then for any , there exists such that for any data with , the system (2.13) has a unique solution on satisfying

 φ∈L∞(0,T0;H3w),ψ∈L∞(0,T0;H3),∇ψ∈L∞(0,T0;H2w) with φ|t=0=φ0,ψ|t=0=ψ0

and

 supt∈[0,T0](∥φ∥2H3w+∥ψ∥2H3w+∥∇ψ∥2H2w)≤2M.

The local solution of (2.13) can be obtained by the usual contraction method and by a similar computation as in the proof of Proposition 2.13, for which we omit its proof (or see [3]).

The following proposition gives a uniform-in-time estimate, which is the main heart of this paper.

###### Proposition 2.13.

For any and any such that the product is sufficiently small, there exist constants , and such that if is a solution of (2.2) for some with (2.3) given by Theorem 2.2, then we have the following:
If be a local solution of (2.13) on for some with , then we have

 M(T)+∫T04∑l=1∥∇lφ∥2w+∫T03∑l=1∥∇lψ∥2w+ϵ∫T0∥∇4ψ∥2w≤C0M(0).

Note that does not depend on .

###### Proof of Proposition 2.11 from Proposition 2.12 and Proposition 2.13.

We include the proof here for readers’ convenience even if this continuation argument is now standard (or see [3]). Let’s take and where and are the constants in Proposition 2.13. Due to , we know . Consider the initial data with . By using the constant to the local-existence result (Proposition 2.12), there exist , and there is the unique local solution on with . Due to , we can use the result of Proposition 2.13 to obtain , which implies . Hence we can extend the solution from the time up to the time by Proposition 2.12 and we obtain . Again by Proposition 2.13, it implies . Thus we can repeat this process of the extension to get for any . ∎

In the rest of the paper, we focus on proving Proposition 2.13.

## 3. Uniform-in-time estimate: Proof of Proposition 2.13

Let and . Recall the Poincaré inequality on intervals which says that there is a constant such that for any and for any , the inequality

 (3.1) ∥f−¯¯¯f∥L2(0,λ)≤λCp∥f′∥L2(0,λ)

holds. Here the mean value of is defined by . We assume that the product is small to have

 (3.2) s⋅λ⋅Cp≤116.

From now on, these values and are fixed until the end of the proof. Let’s assume and which will be taken sufficiently small later in the proof several times.

We suppose first that is sufficiently small so that any meets the assumption of Theorem 2.2. Let be a solution of (2.2) for some with (2.3) given by Theorem 2.2. Let be a local solution of (2.13) on for some with .

In the sequel, denotes a positive constant which may change from line to line, but which stays independent on ANY choice of and as long as the positive parameters and are sufficiently small.

### 3.1. Zero-th order estimate

###### Lemma 3.1.

If the positive constants are sufficiently small, then there exists a constant such that for any ,

 (3.3) ∥ψ∥2+∥φ∥2w+∫t0∥∇φ∥2w+ϵ∫t0∥∇ψ∥2+∫t0∫((N′)2N3|φ|2+P′N(φ1)2+PN′N2(φ2)2)≤C1⋅(∥ψ0∥2+∥φ0∥2w)+ϵ⋅C1∫t0∫P′|ψ|2+C1⋅M(t)⋅∫t0∥∇ψ∥2w.
###### Remark 3.2.

We note that the term

 ∫t0∫((N′)2N3|φ|2+P′N(φ1)2+PN′N2(φ2)2)

in the left-hand side of (3.3) plays a role of dissipation on the zero-th order. This is non-symmetric for and due to the non-symmetric structure of the main equation (2.13). In Lemma 3.3, these localized -norms of will be used to control

 ϵ∫t0∫P′|ψ|2

in the right-hand side of (3.3), which is a localized -norm of multiplied by .

###### Proof.

We multiply to the equation and to the equation:

 1Nφ⋅(φt−sφz−Δφ)+ψ(ψt−sψz−ϵΔψ)=1Nφ⋅(N∇ψ+P∇⋅φ+∇⋅φ∇ψ)+ψ(−2ϵP⋅∇ψ−ϵ|∇ψ|2+∇⋅φ).

Thus we get

 (1N|φ|2/2)t−s(1N|φ|2/2)z+s(1N)′|φ|2/2+1Nφ⋅(−Δφ)+(|ψ|2/2)t−s(|ψ|2/2)z+ψ(−ϵΔψ)=φ⋅∇ψ+PNφ1∇⋅φ+1N(∇⋅φ)φ⋅∇ψ−2ϵPψψz−ϵψ|∇ψ|2+ψ∇⋅φ.

By integrating in space , we have

 12ddt(∫|φ|2N+∫|ψ|2)+∫|∇φ|2N+ϵ∫|∇ψ|2+s2∫|φ|2(1N)′ =−∫φ⋅φz(1N)′+∫PNφ1∇⋅φ+∫φ⋅∇ψN∇⋅φ−2ϵ∫Pψzψ−ϵ∫|∇ψ|2ψ.

Here we use the notation

 |∇φ|2:=2∑i=1|∇φi|2.

Recall the Sobolev embedding which gives us a constant such that for any , the inequality

 ∥f∥L∞≤CSV∥f∥H2

holds. We control the cubic term:

 (3.4)

where we used due to (2.7) and .