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

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:

where $\Omega\subset\mathbb{R}^{N}$ is a bounded domain with a sufficiently smooth boundary $\partial\Omega$ and outer unit normal $\eta$. Equation (1.1a) is doubly nonlinear, since we apply the $p$-Laplacian diffusion operator, where we assume $2\leq p<\infty$, to the integrated diffusion function $A(u):=\int_{0}^{u}a(s)\,ds$, where $a(\cdot)$ is a non-negative integrable function with support on the interval $[0,1]$.

In the biological phenomenon described by (1.1), the quantity $u=u(x,t)$ 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 $\nabla v$ of the concentration $v=v(x,t)$ of a substance called chemoattractant. The movement in response to $\nabla v$ also involves the density-dependent probability $f(u(x,t))$ for a cell located at $(x,t)$ to find space in a neighboring location, and a constant $\chi$ describing chemotactic sensitivity. On the other hand, the PDE (1.1b) describes the diffusion of the chemoattractant, where $d>0$ is a diffusion constant and the function $g(u,v)$ describes the rates of production and degradation of the chemoattractant; we here adopt the common choice

We assume that there exists a maximal population density of cells $u_{\mathrm{m}}$ such that $f(u_{\mathrm{m}})=0$. 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 $\Omega$ after their density attains a certain threshold value, and the chemotactic cross-diffusion term $\chi uf(u)$ vanishes identically when $u\geq u_{\mathrm{m}}$. We also assume that the diffusion coefficient $a(u)$ vanishes at $0$ and $u_{\mathrm{m}}$, so that (1.1a) degenerates for $u=0$ and $u=u_{\mathrm{m}}$, while $a(u)>0$ for $0<u<u_{\mathrm{m}}$. A typical example is $a(u)=\epsilon u(1-u_{\mathrm{m}})$, $\epsilon>0$. Normalizing variables by $\tilde{u}=u/u_{\mathrm{m}}$, $\tilde{v}=v$ and $\tilde{f}(\tilde{u})=f(\tilde{u}u_{\mathrm{m}})$, we have $\tilde{u}_{\mathrm{m}}=1$; 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 $p$-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 $a(u)$ that decays like a power at both degeneracy points (with $p=2$) with the technique by Porzio and Vespri [8] to study the $p$-Laplacian, with $a(u)$ degenerating at only one side. We recover both results as particular cases of the one studied here. To our knowledge, the $p$-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 $p>2$, compared with solutions to the standard case $p=2$, but including nonlinear diffusion.

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 $u=0$ 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 $\mathbb{R}^{N}$ 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 $p$-Laplacian equation $u_{t}=\mathrm{div}\,(|\nabla u|^{p-2}\nabla u)$, $p>1$, 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 $p=2$. 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 $\Delta v$ in $L^{\infty}(Q_{T})$ due to the parabolic nature of (1.1b).

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 $p^{\prime}$ the conjugate exponent of $p$ (we will restrict ourselves to the degenerate case $p\geq 2$): $\frac{1}{p}+\frac{1}{p^{\prime}}=1$. Moreover, $C_{w}(0,T,L^{2}(\Omega))$ denotes the space of continuous functions with values in (a closed ball of) $L^{2}(\Omega)$ endowed with the weak topology, and $\left\langle\cdot,\cdot\right\rangle$ is the duality pairing between $W^{1,p}(\Omega)$ and its dual $(W^{1,p}(\Omega))^{\prime}$.

To ensure, in particular, that all terms and coefficients are sufficiently smooth for this definition to make sense, we require that $f\in C^{1}[0,1]$ and $f(1)=0$, and assume that the diffusion coefficient $a(\cdot)$ has the following properties: $a\in C^{1}[0,1]$, $a(0)=a(1)=0$, and $a(s)>0$ for $0<s<1$. Moreover, we assume that there exist constants $\delta\in(0,1/2)$ and $\gamma_{2}\geq\gamma_{1}>1$ such that

where we define the functions $\smash{\phi(s):=s^{\beta_{1}/(p-1)}}$ and $\smash{\psi(s):=s^{\beta_{2}/(p-1)}}$ for $\beta_{2}>\beta_{1}>0$.

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

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 $a(u)$ by the regularized, strictly positive diffusion coefficient $a_{\varepsilon}(u):=a(u)+\varepsilon$, where $\varepsilon>0$ is the regularization parameter. Once the regularized problem is solved, we send the regularization parameter $\varepsilon$ 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 $\partial_{t}Q_{T}$ the parabolic boundary of $Q_{T}$, define $\smash{\tilde{M}:=\|u\|_{\infty,Q_{T}}}$, and recall the definition of the intrinsic parabolic $p$-distance from a compact set $K\subset Q_{T}$ to $\partial_{t}Q_{T}$ as

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

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 $p$-Laplacian structure of the diffusion term, as well as the fact that the diffusion coefficient $a(u)$ vanishes at $u=0$ and $u=1$. 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 $u$; that of $v$ 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].

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 $p$-Laplacian term, and in the Appendix we give further details about the numerical method used to treat the examples.

We define the new diffusion term $A_{\varepsilon}(s):=A(s)+\varepsilon s$, with $a_{\varepsilon}(s)=a(s)+\varepsilon$, and consider, for each fixed $\varepsilon>0$, the non-degenerate problem

With $\bar{u}\in\mathcal{K}$ fixed, let $v_{\varepsilon}$ be the unique solution of the problem

Given the function $v_{\varepsilon}$, let $u_{\varepsilon}$ be the unique solution of the following quasilinear parabolic problem:

Here $v_{0}$ and $u_{0}$ are functions satisfying the assumptions of Theorem 1.1.

Since for any fixed $\bar{u}\in\mathcal{K}$, (2.2a) is uniformly parabolic, standard theory for parabolic equations [22] immediately leads to the following lemma.

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

We define a map $\Theta:\mathcal{K}\to\mathcal{K}$ such that $\Theta(\bar{u})=u_{\varepsilon}$, where $u_{\varepsilon}$ solves (2.3), i.e., $\Theta$ is the solution operator of (2.3) associated with the coefficient $\bar{u}$ and the solution $v_{\varepsilon}$ coming from (2.2). By using the Schauder fixed-point theorem, we now prove that $\Theta$ has a fixed point. First, we need to show that $\Theta$ is continuous. Let $\smash{\{\bar{u}_{n}\}_{n\in\mathbb{N}}}$ be a sequence in $\mathcal{K}$ and $\bar{u}\in\mathcal{K}$ be such that $\bar{u}_{n}\to\bar{u}$ in $L^{p}(Q_{T})$ as $n\to\infty$. Define $u_{\varepsilon n}:=\Theta(\bar{u}_{n})$, i.e., $u_{\varepsilon n}$ is the solution of (2.3) associated with $\bar{u}_{n}$ and the solution $v_{\varepsilon n}$ of (2.2). To show that $u_{\varepsilon n}\to\Theta(\bar{u})$ in $L^{p}(Q_{T})$, we start with the following lemma.

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

Lemmas 2.2–2.4 imply that there exist $u_{\varepsilon}\in L^{p}(0,T;W^{1,p}(\Omega))$ and $v_{\varepsilon}\in L^{2}(0,T;H^{1}(\Omega))$ such that, up to extracting subsequences if necessary, $u_{\varepsilon n}\to u_{\varepsilon}$ strongly in $L^{p}(Q_{T})$ and $v_{\varepsilon n}\to v_{\varepsilon}$ strongly in $L^{2}(0,T;H^{1}(\Omega))$ as $n\to\infty$, so $\Theta$ is indeed continuous on $\mathcal{K}$. Moreover, due to Lemma 2.3, $\Theta(\mathcal{K})$ is bounded in the set

Similarly to the results of [24], it can be shown that $\mathcal{W}\hookrightarrow L^{p}(Q_{T})$ is compact, and thus $\Theta$ is compact. Now, by the Schauder fixed point theorem, the operator $\Theta$ has a fixed point $u_{\varepsilon}$ such that $\Theta(u_{\varepsilon})=u_{\varepsilon}$. This implies that there exists a solution ($u_{\varepsilon},v_{\varepsilon})$ of

We now pass to the limit $\varepsilon\to 0$ in solutions $(u_{\varepsilon},v_{\varepsilon})$ 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 implies that there exists a constant $C>0$, which does not depend on $\varepsilon$, such that

Notice that, from (2.6) and (2.7), the term $g(u_{\varepsilon},v_{\varepsilon})$ is bounded. Thus, in light of classical results on $L^{r}$ regularity, there exists another constant $C>0$, which is independent of $\varepsilon$, such that

Taking $\varphi=A_{\varepsilon}(u_{\varepsilon})$ as a test function in (2.5) yields

then, using (2.7), the uniform $L^{\infty}$ bound on $u_{\varepsilon}$, an application of Young’s inequality to treat the term $\mathrm{\nabla}v_{\varepsilon}\cdot\mathrm{\nabla}A_{\varepsilon}(u_{\varepsilon})$, and defining $\smash{\mathcal{A}_{\varepsilon}(s):=\int_{0}^{s}A_{\varepsilon}(r)\,dr}$, we obtain

for some constant $C>0$ independent of $\varepsilon$.

Let $\varphi\in L^{p}(0,T;W^{1,p}(\Omega))$. Using the weak formulation (2.5), (2.7) and (2.8), we may follow the reasoning in [4] to deduce the bound

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

To establish the second convergence in (2.10), we have applied the dominated convergence theorem to $\smash{u_{\varepsilon}=A_{\varepsilon}^{-1}(A_{\varepsilon}(u_{\varepsilon}))}$ (recall that $A$ is monotone) and the weak-$\star$ convergence of $u_{\varepsilon}$ to $u$ in $L^{\infty}(Q_{T})$. We also have the following lemma, see [4] for its proof.

Next, we identify $\Gamma_{1}$ as $|\mathrm{\nabla}A(u)|^{p-2}\mathrm{\nabla}A(u)$ when passing to the limit $\varepsilon\to 0$ 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).