Doubly nonlinear chemotaxis model

# On a doubly nonlinear diffusion model of chemotaxis with prevention of overcrowding

Mostafa Bendahmane Raimund Bürger Ricardo Ruiz Baier  and  José Miguel Urbano
July 17, 2019
###### Abstract.

This paper addresses the existence and regularity of weak solutions for a fully parabolic model of chemotaxis, with prevention of overcrowding, that degenerates in a two-sided fashion, including an extra nonlinearity represented by a -Laplacian diffusion term. To prove the existence of weak solutions, a Schauder fixed-point argument is applied to a regularized problem and the compactness method is used to pass to the limit. The local Hölder regularity of weak solutions is established using the method of intrinsic scaling. The results are a contribution to showing, qualitatively, to what extent the properties of the classical Keller-Segel chemotaxis models are preserved in a more general setting. Some numerical examples illustrate the model.

###### Key words and phrases:
Chemotaxis, reaction-diffusion equations, degenerate PDE, parabolic -Laplacian, doubly nonlinear, intrinsic scaling
###### 1991 Mathematics Subject Classification:
AMS Subject Classification: 35K65, 92C17, 35B65
Departamento de Ingeniería Matemática, Universidad de Concepción, Casilla 160-C, Concepción, Chile.
E-mail: mostafab@ing-mat.udec.cl, rburger@ing-mat.udec.cl, rruiz@ing-mat.udec.cl
CMUC, Department of Mathematics, University of Coimbra, 3001-454 Coimbra, Portugal. E-mail: jmurb@mat.uc.pt

## 1. Introduction

### 1.1. Scope

It is the purpose of this paper to study the existence and regularity of weak solutions of the following parabolic system, which is a generalization of the well-known Keller-Segel model [1, 2, 3] of chemotaxis:

 ∂tu−div(|∇A(u)|p−2∇A(u))+div(χuf(u)∇v)=0in QT:=Ω×(0,T),T>0,Ω⊂RN, (1.1a) ∂tv−dΔv=g(u,v)in QT, (1.1b) |∇A(u)|p−2a(u)∂u∂η=0,∂v∂η=0on ΣT:=∂Ω×(0,T), (1.1c) u(x,0)=u0(x),v(x,0)=v0(x)on Ω, (1.1d)

where is a bounded domain with a sufficiently smooth boundary and outer unit normal . Equation (1.1a) is doubly nonlinear, since we apply the -Laplacian diffusion operator, where we assume , to the integrated diffusion function , where is a non-negative integrable function with support on the interval .

In the biological phenomenon described by (1.1), the quantity is the density of organisms, such as bacteria or cells. The conservation PDE (1.1a) incorporates two competing mechanisms, namely the density-dependent diffusive motion of the cells, described by the doubly nonlinear diffusion term, and a motion in response to and towards the gradient of the concentration of a substance called chemoattractant. The movement in response to also involves the density-dependent probability for a cell located at to find space in a neighboring location, and a constant describing chemotactic sensitivity. On the other hand, the PDE (1.1b) describes the diffusion of the chemoattractant, where is a diffusion constant and the function describes the rates of production and degradation of the chemoattractant; we here adopt the common choice

 g(u,v)=αu−βv,α,β≥0. (1.2)

We assume that there exists a maximal population density of cells  such that . This corresponds to a switch to repulsion at high densities, known as prevention of overcrowding, volume-filling effect or density control (see [4]). It means that cells stop to accumulate at a given point of after their density attains a certain threshold value, and the chemotactic cross-diffusion term vanishes identically when . We also assume that the diffusion coefficient vanishes at and , so that (1.1a) degenerates for and , while for . A typical example is , . Normalizing variables by , and , we have ; in the sequel we will omit tildes in the notation.

The main intention of the present work is to address the question of the regularity of weak solutions, which is a delicate analytical issue since the structure of equation (1.1a) combines a degeneracy of -Laplacian type with a two-sided point degeneracy in the diffusive term. We prove the local Hölder continuity of the weak solutions of (1.1) using the method of intrinsic scaling (see [5, 6]). The novelty lies in tackling the two types of degeneracy simultaneously and finding the right geometric setting for the concrete structure of the PDE. The resulting analysis combines the technique used by Urbano [7] to study the case of a diffusion coefficient that decays like a power at both degeneracy points (with ) with the technique by Porzio and Vespri [8] to study the -Laplacian, with degenerating at only one side. We recover both results as particular cases of the one studied here. To our knowledge, the -Laplacian is a new ingredient in chemotaxis models, so we also include a few numerical examples that illustrate the behavior of solutions of (1.1) for , compared with solutions to the standard case , but including nonlinear diffusion.

### 1.2. Related work

To put this paper in the proper perspective, we recall that the Keller-Segel model is a widely studied topic, see e.g. Murray [3] for a general background and Horstmann [1] for a fairly complete survey on the Keller-Segel model and the variants that have been proposed. Nonlinear diffusion equations for biological populations that degenerate at least for were proposed in the 1970s by Gurney and Nisbet [9] and Gurtin and McCamy [10]; more recent works include those by Witelski [11], Dkhil [12], Burger et al. [13] and Bendahmane et al. [4]. Furthermore, well-posedness results for these kinds of models include, for example, the existence of radial solutions exhibiting chemotactic collapse [14], the local-in-time existence, uniqueness and positivity of classical solutions, and results on their blow-up behavior [15], and existence and uniqueness using the abstract theory developed in [16], see [17]. Burger et al. [13] prove the global existence and uniqueness of the Cauchy problem in for linear and nonlinear diffusion with prevention of overcrowding. The model proposed herein exhibits an even higher degree of nonlinearity, and offers further possibilities to describe chemotactic movement; for example, one could imagine that the cells or bacteria are actually placed in a medium with a non-Newtonian rheology. In fact, the evolution -Laplacian equation , , is also called non-Newtonian filtration equation, see [18] and [19, Chapter 2] for surveys. Coming back to the Keller-Segel model, we also mention that another effort to endow this model with a more general diffusion mechanism has recently been made by Biler and Wu [20], who consider fractional diffusion.

Various results on the Hölder regularity of weak solutions to quasilinear parabolic systems are based on the work of DiBenedetto [5]; the present article also contributes to this direction. Specifically for a chemotaxis model, Bendahmane, Karlsen, and Urbano [4] proved the existence and Hölder regularity of weak solutions for a version of (1.1) for . For a detailed description of the intrinsic scaling method and some applications we refer to the books [5, 6].

Concerning uniqueness of solution, the presence of a nonlinear degenerate diffusion term and a nonlinear transport term represents a disadvantage and we could not obtain the uniqueness of a weak solution. This contrasts with the results by Burger et al. [13], where the authors prove uniqueness of solutions for a degenerate parabolic-elliptic system set in an unbounded domain, using a method which relies on a continuous dependence estimate from [21], that does not apply to our problem because it is difficult to bound in due to the parabolic nature of (1.1b).

### 1.3. Weak solutions and statement of main results

Before stating our main results, we give the definition of a weak solution to (1.1), and recall the notion of certain functional spaces. We denote by the conjugate exponent of (we will restrict ourselves to the degenerate case ): . Moreover, denotes the space of continuous functions with values in (a closed ball of) endowed with the weak topology, and is the duality pairing between and its dual .

###### Definition 1.1.

A weak solution of (1.1) is a pair of functions satisfying the following conditions:

 0≤u(x,t)≤1 and v(x,t)≥0 for a.e.~{}(x,t)∈QT, u∈Cw(0,T,L2(Ω)),∂tu∈Lp′(0,T;(W1,p(Ω))′),u(0)=u0, A(u)=∫u0a(s)ds∈Lp(0,T;W1,p(Ω)),

and, for all and ,

 ∫T0⟨∂tu,φ⟩dt+∬QT{|∇A(u)|p−2∇A(u)−χuf(u)∇v}⋅∇φdxdt=0, ∫T0⟨∂tv,ψ⟩dt+d∬QT∇v⋅∇ψdxdt=∬QTg(u,v)ψdxdt.

To ensure, in particular, that all terms and coefficients are sufficiently smooth for this definition to make sense, we require that and , and assume that the diffusion coefficient has the following properties: , , and for . Moreover, we assume that there exist constants and such that

 γ1ϕ(s)≤a(s)≤γ2ϕ(s)for s∈[0,δ],γ1ψ(1−s)≤a(s)≤γ2ψ(1−s)for s∈[1−δ,1], (1.3)

where we define the functions and for .

Our first main result is the following existence theorem for weak solutions.

###### Theorem 1.1.

If with and a.e. in , then there exists a weak solution to the degenerate system (1.1) in the sense of Definition 1.1.

In Section 2, we first prove the existence of solutions to a regularized version of (1.1) by applying the Schauder fixed-point theorem. The regularization basically consists in replacing the degenerate diffusion coefficient by the regularized, strictly positive diffusion coefficient , where is the regularization parameter. Once the regularized problem is solved, we send the regularization parameter  to zero to produce a weak solution of the original system (1.1) as the limit of a sequence of such approximate solutions. Convergence is proved by means of a priori estimates and compactness arguments.

We denote by the parabolic boundary of , define , and recall the definition of the intrinsic parabolic -distance from a compact set to as

 p-dist(K;∂tQT):=inf(x,t)∈K, (y,s)∈∂tQT(|x−y|+~M(p−2)/p|t−s|1/p).

Our second main result is the interior local Hölder regularity of weak solutions.

###### Theorem 1.2.

Let be a bounded local weak solution of (1.1) in the sense of Definition 1.1, and . Then is locally Hölder continuous in , i.e., there exist constants and , depending only on the data, such that, for every compact ,

 ∣∣u(x1,t1)−u(x2,t2)∣∣≤γ~M{|x1−x2|+~M(p−2)/p|t2−t1|1/pp-dist(K;∂tQT)}α,∀(x1,t1),(x2,t2)∈K.

In Section 3, we prove Theorem 1.2 using the method of intrinsic scaling. This technique is based on analyzing the underlying PDE in a geometry dictated by its own degenerate structure, that amounts, roughly speaking, to accommodate its degeneracies. This is achieved by rescaling the standard parabolic cylinders by a factor that depends on the particular form of the degeneracies and on the oscillation of the solution, and which allows for a recovery of homogeneity. The crucial point is the proper choice of the intrinsic geometry which, in the case studied here, needs to take into account the -Laplacian structure of the diffusion term, as well as the fact that the diffusion coefficient vanishes at and . At the core of the proof is the study of an alternative, now a standard type of argument [5]. In either case the conclusion is that when going from a rescaled cylinder into a smaller one, the oscillation of the solution decreases in a way that can be quantified.

In the statement of Theorem 1.2 and its proof, we focus on the interior regularity of ; that of  follows from classical theory of parabolic PDEs [22]. Moreover, standard adaptations of the method are sufficient to extend the results to the parabolic boundary, see [5, 23].

### 1.4. Outline

The remainder of the paper is organized as follows: Section 2 deals with the general proof of our first main result (Theorem 1.1). Section 2.1 is devoted to the detailed proof of existence of solutions to a non-degenerate problem; in Section 2.2 we state and prove a fixed-point-type lemma, and the conclusion of the proof of Theorem 1.1 is contained in Section 2.3. In Section 3 we use the method of intrinsic scaling to prove Theorem 1.2, establishing the Hölder continuity of weak solutions to (1.1). Finally, in Section 4 we present two numerical examples showing the effects of prevention of overcrowding and of including the -Laplacian term, and in the Appendix we give further details about the numerical method used to treat the examples.

## 2. Existence of solutions

We first prove the existence of solutions to a non-degenerate, regularized version of problem (1.1), using the Schauder fixed-point theorem, and our approach closely follows that of [4]. We define the following closed subset of the Banach space :

 K:={u∈Lp(QT):0≤u(x,t)≤1for a.e. (x,t)∈QT}.

### 2.1. Weak solution to a non-degenerate problem

We define the new diffusion term , with , and consider, for each fixed , the non-degenerate problem

 ∂tuε−div(|∇Aε(uε)|p−2∇Aε(uε))+div(χf(uε)∇vε)=0in QT, (2.1a) ∂tvε−dΔvε=g(uε,vε)in QT, (2.1b) |∇Aε(uε)|p−2aε(uε)∂uε∂η=0,∂vε∂η=0on ΣT, (2.1c) uε(x,0)=u0(x),vε(x,0)=v0(x)for x∈Ω. (2.1d)

With fixed, let be the unique solution of the problem

 ∂tvε−dΔvε=g(¯u,vε)in QT, (2.2a) ∂vε∂η=0on ΣT,vε(x,0)=v0(x)for x∈Ω. (2.2b)

Given the function , let be the unique solution of the following quasilinear parabolic problem:

 ∂tuε−div(|∇Aε(uε)|p−2∇Aε(uε))+div(χuεf(uε)∇vε)=0in QT, (2.3a) |∇Aε(uε)|p−2aε(uε)∂uε∂η=0on ΣT,uε(x,0)=u0(x)for % x∈Ω. (2.3b)

Here and  are functions satisfying the assumptions of Theorem 1.1.

Since for any fixed , (2.2a) is uniformly parabolic, standard theory for parabolic equations [22] immediately leads to the following lemma.

###### Lemma 2.1.

If , then problem (2.2) has a unique weak solution , for all , satisfying in particular

 ∥vε∥L∞(QT)+∥vε∥L∞(0,T;L2(Ω))≤C,∥vε∥L2(0,T;H1(Ω))≤C,∥∂tvε∥L2(QT)≤C, (2.4)

where is a constant that depends only on , , , and .

The following lemma (see [22]) holds for the quasilinear problem (2.3).

###### Lemma 2.2.

If , then, for any , there exists a unique weak solution to problem (2.3).

### 2.2. The fixed-point method

We define a map such that , where solves (2.3), i.e., is the solution operator of (2.3) associated with the coefficient  and the solution coming from (2.2). By using the Schauder fixed-point theorem, we now prove that has a fixed point. First, we need to show that is continuous. Let be a sequence in and be such that in as . Define , i.e., is the solution of (2.3) associated with and the solution of (2.2). To show that in , we start with the following lemma.

###### Lemma 2.3.

The solutions to problem (2.3) satisfy

• for a.e. .

• The sequence is bounded in .

• The sequence is relatively compact in .

###### Proof.

The proof follows from that of Lemma 2.3 in [4] if we take into account that is uniformly bounded in . ∎

The following lemma contains a classical result (see [22]).

###### Lemma 2.4.

There exists a function such that the sequence converges strongly to in .

Lemmas 2.22.4 imply that there exist and such that, up to extracting subsequences if necessary, strongly in and strongly in as , so is indeed continuous on . Moreover, due to Lemma 2.3, is bounded in the set

 W:={u∈Lp(0,T;W1,p(Ω)):∂tu∈Lp′(0,T;(W1,p(Ω))′)}.

Similarly to the results of [24], it can be shown that is compact, and thus is compact. Now, by the Schauder fixed point theorem, the operator has a fixed point such that . This implies that there exists a solution ( of

 (2.5) ∀φ∈Lp(0,T;W1,p(Ω)) and ∀ψ∈L2(0,T;H1(Ω)).

### 2.3. Existence of weak solutions

We now pass to the limit in solutions to obtain weak solutions of the original system (1.1). From the previous lemmas and considering (2.1b), we obtain the following result.

###### Lemma 2.5.

For each fixed , the weak solution to (2.1) satisfies the maximum principle

 0≤uε(x,t)≤1\em andvε(x,t)≥0\em for a.e. (x,t)∈QT. (2.6)

Moreover, the first two estimates of (2.4) in Lemma 2.1 are independent of .

Lemma 2.5 implies that there exists a constant , which does not depend on , such that

 ∥vε∥L∞(QT)+∥vε∥L∞(0,T;L2(Ω))≤C, ∥vε∥L2(0,T;H1(Ω))≤C. (2.7)

Notice that, from (2.6) and (2.7), the term is bounded. Thus, in light of classical results on regularity, there exists another constant , which is independent of , such that

 ∥∂tvε∥Lr(QT)+∥vε∥Lr(0,T;W1,r(Ω))≤C for all r>1.

Taking as a test function in (2.5) yields

 ∫T0⟨∂tuε,A(uε)⟩dt+ε∫T0⟨∂tuε,uε⟩dt+∬QT|∇Aε(uε)|pdxdt−∬QTχf(uε)∇vε⋅∇Aε(uε)dxdt=0;

then, using (2.7), the uniform bound on , an application of Young’s inequality to treat the term , and defining , we obtain

 sup0≤t≤T∫ΩAε(uε)(x,t)dx+εsup0≤t≤T∫Ω|uε(x,t)|22dx+∬QT|∇Aε(uε)|pdxdt≤C (2.8)

for some constant independent of .

Let . Using the weak formulation (2.5), (2.7) and (2.8), we may follow the reasoning in [4] to deduce the bound

 ∥∂tuε∥Lp′(0,T;(W1,p(Ω))′)≤C. (2.9)

Therefore, from (2.7)–(2.9) and standard compactness results (see [24]), we can extract subsequences, which we do not relabel, such that, as ,

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Aε(uε)→A(u) % strongly in Lp(QT) and a.e.,uε→u strongly in Lq(QT) for all q≥1,vε→v strongly in L2(QT),∇vε→∇v weakly in L2(QT) and ∇Aε(uε)→∇A(u) weakly in Lp(QT),|∇Aε(uε)|p−2∇Aε(uε)→Γ1 weakly in Lp′(QT),vε→v weakly in L2(0,T;H1(Ω)),∂tuε→∂tu weakly in Lp′(0,T;(W1,p(Ω))′) and ∂tvε→∂tv weakly in L2(0,T;(H1(Ω))′). (2.10)

To establish the second convergence in (2.10), we have applied the dominated convergence theorem to (recall that is monotone) and the weak- convergence of to in . We also have the following lemma, see [4] for its proof.

###### Lemma 2.6.

The functions converge strongly to in as .

Next, we identify as when passing to the limit in (2.5). Due to this particular nonlinearity, we cannot employ the monotonicity argument used in [4]; rather, we will utilize a Minty-type argument [25] and make repeated use of the following “weak chain rule” (see e.g. [26] for a proof).

###### Lemma 2.7.

Let be Lipschitz continuous and nondecreasing. Assume is such that , , a.e. on , with . If we define , then

 −∫s0⟨∂tu,b(u)ϕ⟩dt=∫s0∫ΩB(u)∂tϕdxdt+∫ΩB(u0)ϕ(x,0)dx−∫ΩB(u(x,s))ϕ(x,s)dx

holds for all and for any .

###### Lemma 2.8.

There hold and strongly in .

###### Proof.

We define . The first step will be to show that

 ∭QT(Γ1−|∇σ|p−2∇σ)⋅(∇A(u)−∇σ)dxdsdt≥0,∀σ∈Lp(0,T;W1,p(Ω)). (2.11)

For all fixed , we have the decomposition

 I1:=∭QT|∇Aε(uε)|p−2∇Aε(uε)⋅(∇A(u)−∇Aε(uε))dxdsdt, I2:=∭QT(|∇Aε(uε)|p−2∇Aε(uε)−|∇σ|p−2∇σ)⋅(∇Aε(uε)−∇σ)dxdsdt, I3:=∭QT|∇σ|p−2∇σ⋅(∇Aε(uε)−∇A(u))dxdsdt.

Clearly, and from (2.10) we deduce that as . For , if we multiply (2.1a) by and integrate over , we obtain

Now, if we take and use Lemma 2.7, we obtain

 I1= −∫T0∫t0⟨∂tuε,A(u)⟩dsdt+∫T0∫t0⟨∂tuε,Aε(uε)⟩dsdt +∭QTχuεf(uε)∇vε⋅(∇A(u)−∇Aε(uε))dxdsdt = −∫T0∫t0⟨∂tuε,A(u)⟩dsdt+∬QTAε(uε)dxdt−T∫ΩAε(u0)dx +∭QTχuεf(uε)∇vε⋅(∇A(u)−∇Aε(uε))dxdsdt.

Therefore, using (2.10) and Lemma 2.6 and defining , we conclude that

 limε→0I1=−∫T0∫t0⟨∂tu,A(u)⟩dsdt+∫T0∫ΩA(u(x,t))dxdt−T∫ΩA(u0(x))dx,

and from Lemma 2.7, this yields as . Consequently, we have shown that

 limε→0∭QT(|∇Aε(uε)|p−2∇Aε(uε)−|∇σ|p−2∇σ)⋅(∇A(u)−∇σ)dxdsdt≥0,

which proves (2.11). Choosing with and and combining the two inequalities arising from and , we obtain the first assertion of the lemma. The second assertion directly follows from (2.11). ∎

With the above convergences we are now able to pass to the limit , and we can identify the limit as a (weak) solution of (1.1). In fact, if is a test function for (2.5), then by (2.10) it is now clear that

 ∬QT|∇Aε(uε)|p−2∇Aε(uε)⋅∇φdxdt→∬QT|∇A(u)|p−2∇A(u)⋅∇φdxdtas ε→0.

Since is bounded in and by Lemma 2.6, in , it follows that

 ∬QTχuεf(uε)∇vε⋅∇φdxdt→∬QTχuf(u)∇v⋅∇φdxdtas ε→0.

We have thus identified as the first component of a solution of (1.1). Using a similar argument, we can identify as the second component of a solution.

## 3. Hölder continuity of weak solutions

### 3.1. Preliminaries

We start by recasting Definition 1.1 in a form that involves the Steklov average, defined for a function and by

 wh:=⎧⎪⎨⎪⎩1h∫t+htw(⋅,τ)dτif t∈(0,T−h],0if t∈(T−h,T].
###### Definition 3.1.

A local weak solution for (1.1) is a measurable function such that, for every compact and for all ,

 ∫K×{t}{∂t(uh)φ+(|∇A(u)|p−2∇A(u))h⋅∇φ−(χuf(u)∇v)h⋅∇φ}dx=0,∀φ∈W1,p0(K). (3.1)

The following technical lemma on the geometric convergence of sequences (see e.g., [27, Lemma 4.2, Ch. I]) will be used later.

###### Lemma 3.1.

Let and , , be sequences of positive real numbers satisfying

where , , and are given constants. Then as provided that

 X0+Z1+κ0≤(2C)−(1+κ)/σb−(1+κ)/σ2,\em with σ=min{α,κ}.

### 3.2. The rescaled cylinders

Let denote the ball of radius centered at . Then, for a point , we denote the cylinder of radius and height by

 (x0,t0)+Q(τ,ρ):=Bρ(x0)×(t0−τ,t0).

Intrinsic scaling is based on measuring the oscillation of a solution in a family of nested and shrinking cylinders whose dimensions are related to the degeneracy of the underlying PDE. To implement this, we fix ; after a translation, we may assume that . Then let and let be small enough so that , and define

 μ+:=esssupQ(Rp−ε,2R)u,μ−:=essinfQ(Rp−ε,2R)u,ω:=essoscQ(Rp−ε,2R)u≡μ+−μ−.

Now construct the cylinder , where

 a0=(ω2)2−p1ϕ(ω/2m)p−1,

with to be chosen later. To ensure that , we assume that

 1a0=(ω2)p−2ϕ(ω2m)p−1>Rε, (3.2)

and therefore the relation

 essoscQ(a0Rp,R)u≤ω (3.3)

holds. Otherwise, the result is trivial as the oscillation is comparable to the radius. We mention that for small and for , the cylinder is long enough in the direction, so that we can accommodate the degeneracies of the problem. Without loss of generality, we will assume .

Consider now, inside , smaller subcylinders of the form

These are contained in if , which holds whenever and

 t∗∈((ω/2)2−pRpψ(ω/4)p−1−(ω/2)p−2Rpϕ(ω/2m)p−1,0).

These particular definitions of and of turn out to be the natural extensions to the case of their counterparts in [7]. Notice that for and , we recover the standard parabolic cylinders.

The structure of the proof will be based on the analysis of the following alternative: either there is a cylinder where is essentially away from its infimum, or such a cylinder can not be found and thus is essentially away from its supremum in all cylinders of that type. Both cases lead to the conclusion that the essential oscillation of within a smaller cylinder decreases by a factor that can be quantified, and which does not depend on .

###### Remark 3.1.

(See [8, Remark 4.2]) Let us introduce quantities of the type , where and