# Analysis of a fully discrete approximation for the classical Keller–Segel model: lower and a priori bounds

This paper is devoted to constructing approximate solutions for the classical Keller–Segel model governing chemotaxis. It consists of a system of nonlinear parabolic equations, where the unknowns are the average density of cells (or organisms), which is a conserved variable, and the average density of chemioattranct. The numerical proposal is made up of a crude finite element method together with a mass lumping technique and a semi-implicit Euler time integration. The resulting scheme turns out to be linear and decouples the computation of variables. The approximate solutions keep lower bounds – positivity for the cell density and nonnegativity for the chemioattranct density –, are mass conservative, satisfy a discrete energy law, and have a priori energy estimates. The latter is achieved by means of a discrete Moser–Trudinger inequality. As far as we know, our numerical method is the first one that can be encountered in the literature dealing with all of the previously mentioned properties at the same time.

There are no comments yet.

## Authors

• 1 publication
• 1 publication
10/13/2019

### A Conservative Finite Element Method for the Incompressible Euler Equations with Variable Density

We construct a finite element discretization and time-stepping scheme fo...
09/22/2021

### Numerical analysis of a finite element formulation of the P2D model for Lithium-ion cells

The mathematical P2D model is a system of strongly coupled nonlinear par...
11/07/2020

### Bresse-Timoshenko type systems with thermodiffusion effects: Well-possedness, stability and numerical results

Bresse-Timoshenko beam model with thermal, mass diffusion and theormoela...
05/31/2021

### Energy-preserving fully-discrete schemes for nonlinear stochastic wave equations with multiplicative noise

In this paper, we focus on constructing numerical schemes preserving the...
08/09/2019

### Numerical analysis for a chemotaxis-Navier-Stokes system

In this paper we develop a numerical scheme for approximating a d-dimens...
12/07/2021

### Finite Element numerical schemes for a chemo-attraction and consumption model

This work is devoted to design and study efficient and accurate numerica...
07/05/2017

### On the Fusion of Compton Scatter and Attenuation Data for Limited-view X-ray Tomographic Applications

In this paper we demonstrate the utility of fusing energy-resolved obser...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

### 1.1. Aims

In 1970/71 Keller and Segel [15, 16] attempted to derive a set of equations for modeling chemotaxis – a biological process through which an organism (or a cell) migrates in response to a chemical stimulus being attractant or repellent. It is nowadays well-known that the work of Keller and Segel turned out to be somehow biologically inaccurate since their equations provide unrealistic solutions; a little more precisely, solutions that blow up in finite time. Such a phenomenon does not occur in nature. Even though the original Keller–Segel equations are less relevant from a biological point of view, they are mathematically of great interest.

Much of work for the Keller–Segel equations has been carried out in developing purely analytical results, whereas there is very few numerical results in the literature. This is due to the fact that solving numerically the Keller–Segel equations is a challenging task because their solutions exhibit many interesting mathematical properties which are not easily adapted to a discrete framework. For instance, solutions to the Keller–Segel equations satisfy lower bounds (positivity and non-negativity) and enjoy an energy law, which is obtained by testing the equations against non linear functions. Cross-diffusion mechanisms governing the chemotactic phenomena are the responsible for the Keller–Segel equations to be so difficult not only theoretically but also numerically.

In spite of being a limited model, it is hoped that developing and analyzing numerical methods for the classical Keller–Segel equations may open new roads to deeper insights and better understandings for dealing with the numerical approximation of other chemotaxis models – biologically more realistic –, but which are, however, inspired on the original Keller–Segel formulation. In a nutshell, these other chemotaxis models are modifications of the Keller–Segel equations in order to avoid the non-physical blow up of solutions and hence produce solution being closer to chemotaxis phenomena. For these other chemotaxis models, it is recommended the excellent surveys of Hillen and Painter [11], Horstamann [12, 13], and, more recently, Bellomo, Bellouquid, Tao, and Winkler [1]. In these surveys the authors reviewed to date as to when they were written the state of art of modeling and mathematical analysis for the Keller–Segel equations and their variants.

It is our aim in this work to design a fully discrete algorithm for the classical Keller–Segel equations based on a finite element discretization whose discrete solutions satisfy lower and a priori bounds.

### 1.2. The Keller–Segel equations

Let be a bounded domain, with

being its outward-directed unit normal vector to

, and let be a time interval. Take and . Then the boundary-value problem for the Keller–Segel equations reads a follows. Find and satisfying

 (1) {∂tu−Δu=−∇⋅(u∇v) in Q,∂tv−Δv=u−v in Q,

subject to the initial conditions

 (2) u(0)=u0 and v(0)=v0 in Ω,

and the boundary conditions

 (3) ∇u⋅n=0 and ∇v⋅n=0 on Σ.

Here denotes is the average density of organisms (or cells), which is a conserved variable, and is the average density of chemioattranct, which is a nonconserved variable.

System (1) was motivated by Keller and Segel [15] describing the aggregation phenomena exhibited by the amoeba Dictyostelium discoideum due to an attractive chemical substance referred to as chemoattractant, which is generated by the own amoeba and is, nevertheless, degraded by living conditions. Moreover, diffusion is also presented in the motion of ameobae and chemoattractant.

The diffusion phenomena performed by cells and chemoattractant are modelled by the terms and , respectively, while the aggregation mechanism is described by the term . It is this nonlinear term that is the major difficulty in studying system (1). Further the production and degradation of chemoattractant are associated with the term .

Concerning the mathematical analysis for system (1), Nagai, Senba, and Yoshida [17] proved existence, uniqueness and regularity of solutions under the condition by means of one particular instance of Moser–Tridinguer’s inequality; in the particular case that be a ball, such a condition becomes , and Herrero and Velázquez [10] dealt with the first blow up framework by constructing some radially symmetric solutions which blow up within finite time. The next progress in this sense with being non-radial and simply connected was the work of Horstmann and Wang [14] who found some unbounded solutions provided that and . So far there is no supporting evidence as to whether such solutions may evolve to produce a blow-up phenomenon within finite time or whether, on the contrary, may increase to infinity with time. The main tool in proving blow-up solutions is the energy law which stems from system (1). Observe that an inadequate approximation of lower bounds can trigger oscillations of the variables, which can lead to spurious, blow-up solutions.

Concerning the numerical analysis for system (1), very little is said about numerical algorithms which keep lower bounds, are mass conservative and have a discrete energy law. Proper numerical treatment of these properties is made difficult by the fact that the non-linearity occurs in the highest order derivative. Numerical algorithms are mainly designed so as to keep lower bounds. We refer the reader to [18, 5, 22, 3, 21] which mainly deal with lower bounds and mass conservation. As far as we are concerned, there is no numerical method coping with a discrete energy law.

### 1.3. Notation

We collect here as a reference some standard notation used throughout the paper. For , we denote by the usual Lebesgue space, i.e.,

 Lp(Ω)={v:Ω→R:v Lebesgue-measurable,∫Ω|v(x)|pdx<∞}.

or

 L∞(Ω)={v:Ω→R:v Lebesgue-measurable,esssupx∈Ω|v(x)|<∞}.

This space is a Banach space endowed with the norm if or if . In particular, is a Hilbert space. We shall use for its inner product and for its norm.

Let be a multi-index with , and let be the differential operator such that

 ∂α=(∂∂x1)α1(∂∂x2)α2.

For and , we consider to be the Sobolev space of all functions whose derivatives are in , i.e.,

 Wm,p(Ω)={v∈Lp(Ω):∂kv∈L2(Ω) ∀ |k|≤m}

associated to the norm

 ∥f∥Wm,p(Ω) =⎛⎝∑|α|≤m∥∂αf∥pLp(Ω)⎞⎠1/p for 1≤p<∞, ∥f∥Wm,p(Ω) =max|α|≤m∥∂αf∥L∞(Ω), for p=∞.

For , we denote . Moreover, we make of use the space

 H2N(Ω)={v∈H2(Ω):∫Ωv(x)dx=0 and ∂nv=0 on ∂Ω},

for which is known that .

### 1.4. Outline

The remainder of this paper is organized in the following way. In the next section we state our finite element space and some tools. In particular, we prove a discrete version of a variant of Moser–Trudinger’s inequality. In section 4, we apply our ideas to discretize system (1) in space and time for defining our numerical method and formulate our main result. Next is section 5 dedicated to demonstrating lower bounds, a discrete energy law, and a priori bounds all of which are local in time for approximate solutions. This is accomplished in a series of lemmas where the final argument is an induction procedure on the time step so as to obtain the above mentioned properties globally in time.

## 2. Technical preliminaries

This section is mainly devoted to setting out the hypotheses and some auxiliary results concerning the finite element space that will use throughout this work.

### 2.1. Hypotheses

We consider the finite element approximation of (1) under the following assumptions on the domain, the mesh and the finite element space. Moreover, these assumptions will be mentioned on stating our main result.

1. Let be a convex, bounded domain of with a polygonal boundary, and let be the minimum interior angle at the vertices of .

2. Let be a family of acute, shape-regular, quasi-uniform triangulations of made up of triangles, so that , where , with being the diameter of . More precisely, we assume that

1. there exists , independent of , such that

 min{diamBT:T∈Th}≥αh,

where is the largest ball contained in , and

2. there exists such that every angle between two edges of a triangle is bounded by .

Further, let be the coordinates of the nodes of .

3. Associated with is the finite element space

 Xh={xh∈C0(¯¯¯¯Ω):xh|T∈P1(T), ∀T∈Th},

where is the set of linear polynomials on . Let be the standard basis functions for .

### 2.2. Auxiliary results

In the subsequent sections, we will make use of some technical results concerning the above-mentioned hypotheses.

A key tool for proving lower bounds for the finite element approximation of (1) is the acuteness of .

###### Proposition 2.1.

Let be a polygonal. Consider to be constructed over being acute. Then, for each with vertices , there exists a constant , depending on , but otherwise independent of and , such that

 (4) ∫T∇φai⋅∇φajdx≤−Cneg

for all with , and

 (5) ∫T∇φai⋅∇φaidx≥Cneg

for all .

A proof of (4) and (5) can be found in [9].

Both local and global finite element properties for

will be needed such as inverse estimates and bounds for the interpolation error. We first recall some local inverse estimates. See

[2, Lem. 4.5.3] or [6, Lem. 1.138] for a proof.

###### Proposition 2.2.

Let be polygonal. Consider to be constructed over being quasi-uniform. Then, for each and , there exists a constant , independent of and , such that, for all ,

 (6) ∥∇xh∥Lp(T)≤Cinvh−1∥xh∥Lp(T)

and

 (7) ∥∇xh∥L∞(T)≤Cinvh−2p∥∇xh∥Lp(T).

Concerning global inverse inequalities, we need the following.

###### Proposition 2.3.

Let be polygonal. Consider to be constructed over being quasi-uniform. Then for each , there exists a constant , independent of , such that, for all ,

 (8) ∥xh∥L∞(Ω)≤Cinvh−1∥xh∥,
 (9) ∥∇xh∥Lp(Ω)≤Cinvh−1∥xh∥Lp(Ω),
 (10) ∥∇xh∥Lp(Ω)≤Cinvh−2(12−1p)∥∇xh∥L2(Ω),

and

 (11) ∥∇xh∥L∞(Ω)≤Cinvh−2p∥∇xh∥Lp(Ω).

We introduce , the standard nodal interpolation operator, such that for all . Associated with , a discrete inner product on is defined by

 (xh,¯xh)h=∫ΩIh(xh(x)¯xh(x))dx.

We also introduce

 ∥xh∥h=(xh,xh)12h.

Local and global error bounds for are as follows (c.f. [2, Thm. 4.4.4] or [6, Thm. 1.103] for a proof).

###### Proposition 2.4.

Let be polygonal. Consider to be constructed over being quasi-uniform. Then, for each , there exists , independent of and , such that

 (12) ∥φ−Ihφ∥L1(T)≤Capph2∥∇2φ∥L1(T)∀φ∈W2,1(T).
###### Proposition 2.5.

Let be polygonal. Consider to be constructed over being quasi-uniform. Then it follows that there exists , independent of , such that

 (13) ∥∇(φ−Ihφ)∥L2(Ω)≤Capph∥Δφ∥L2(Ω)∀φ∈H2(Ω).
###### Corollary 2.6.

Let be polygonal. Consider to be constructed over being quasi-uniform. Let . Then it follows that there exist two positive constants and , independent of , such that

 (14) ∥xnh−Ih(xnh)∥L1(Ω)≤Cappn(n−1)h2∫Ω|xh(x)|n−2|∇xh(x)|2dx,
 (15) ∥xh¯¯¯xh−Ih(xh¯¯¯xh)∥L1(Ω)≤Ccomh∥xh∥L2(Ω)∥∇¯¯¯xh∥L2(Ω)

and

 (16) ∥xnh∥L1(Ω)≤∥Ih(xnh)∥L1(Ω)≤Csta∥xnh∥L1(Ω),

where depends on .

###### Proof.

Let and compute

 ∇2(xnh)=n(n−1)xn−2hd∑i,j=1∂ixh∂jxh on T.

Then, from (12) and the above identity, we have

 ∥xnh−Ih(xkh)∥L1(T)≤Capph2K∥∇2(xnh)∥L1(T)≤Cappn(n−1)h2∫T|xh(x)|n−2|∇xh(x)|2dx.

Summing over yields (14). The proof of (15) follows very closely the arguments of (14) for . The first part of assertion (16) is a simple application of Jensen’s inequality, while the second part follows from (14) on using Hölder’s inequality, (9) for and, later on, Minkowski’s inequality. ∎

The proof of the following proposition can be found in [4]. It is a generalization of a Moser–Trudinger-type inequality.

###### Proposition 2.7 (Moser-Trudinger).

Let be polygonal with being the minimum interior angle at the vertices of . Then there exists a constant depending on such that for all with and , it follows that

 (17) ∫Ωeα|u(x)|2dx≤CΩ

where .

###### Corollary 2.8.

Let be polygonal with being the minimum interior angle at the vertices of . Consider to be constructed over being quasi-uniform. Let with . Then it follows that there exists a constant , independent of , such that

 (18) ∫ΩIh(euh(x))dx≤CΩ(1+CMT∥∇uh∥2)e18θΩ∥∇uh∥2+1|Ω|∥uh∥L1(Ω).
###### Proof.

From (14), we have

 (19) ∫ΩIh(euh(x))dx=∫Ω(1+uh(x))dx+∞∑n=21n!∫ΩIh(unh(x))dx≤∞∑n=01n!∫Ωunh(x)dx+∞∑n=2Cappn(n−1)h2n!∫Ω|∇uh(x)|2un−2h(x)dx=∫Ω(1+Capph2|∇uh(x)|2)euh(x)dx.

Let with . Young’s inequality gives

 (20) uh=∥∇uh∥vh+m≤18θΩ∥∇uh∥2+2θΩ|vh|2+m.

Thus, combining (19) and (20) yields, on noting (7) for and (17), that

 ∫ΩIh(euh(x))dx≤e18θΩ∥∇uh∥2+m∫Ω(1+Capph2∥∇uh(x)∥2|∇vh|2)e2θΩ|vh(x)|2dx≤e18θΩ∥∇uh∥2+m(1+CappCinv∥∇uh(x)∥2)∫Ωe2θΩ|vh(x)|2dx≤CΩ(1+CappCinv∥∇uh(x)∥2)e18θΩ∥∇uh∥2+m.

An (average) interpolation operator into will be required in order to properly initialize our numerical method. We refer to [19, 7].

###### Proposition 2.9.

Let be polygonal. Consider to be constructed over being quasi-uniform. Then there exists an (average) interpolation operator from to such that

 (21) ∥Qhψ∥Ws,p(Ω)≤Csta∥ψ∥Ws,p(Ω)for s=0,1 and 1≤p≤∞,

and

 (22) ∥Qh(ψ)−ψ∥Ws,p(Ω)≤Capph1+m−s∥ψ∥Wm+1,p(Ω)for 0≤s≤m≤1.

Moreover, let be defined from to as

 (23) −(~Δhxh,¯xh)h=(∇xh,∇¯xh) for all ¯xh∈Xh,

and let be such that

 (24) {−Δx(h)=−~Δhxh in% Ω,∂nx(h)=0 on ∂Ω.

From elliptic regularity theory, the well-posedness of (24) is ensured by the convexity assumption stated in and

 ∥x(h)∥H2N(Ω)≤C∥−˜Δhxh∥.

See [8] for a proof.

###### Proposition 2.10.

Let be a convex polygon. Consider to be constructed over being quasi-uniform. Then there exists a constant , independent of , such that

 (25) ∥∇(x(h)−xh)∥L2(Ω)≤CLaph∥~Δhxh∥L2(Ω).
###### Proof.

We refer the reader to [9] for a proof which uses (13) and (15). ∎

###### Corollary 2.11.

Let be a convex polygon. Consider to be constructed over being quasi-uniform. Then, for each , there exists a constant , independent of , such that

 (26) ∥∇xh∥Lp(Ω)≤Csta∥−˜Δhxh∥.
###### Proof.

The triangle inequality gives

 ∥∇xh∥Lp(Ω)≤∥∇xh−∇Qhx(h)∥Lp(Ω)+∥∇Qhx(h)−∇x(h)∥Lp(Ω)+∥∇x(h)∥Lp(Ω)

and hence applying (10), (25), (21), (22), and Sobolev’s inequality yields (26).

## 3. Presentation of main result

We now define our numerical approximation of system (1). Assume that with and a. e. in .

We begin by approximating the initial data by as follows. Define

 (27) u0h=Qhu0,

which satisfies

 (28) u0h>0 a. e. in Ω,∥u0h∥L1(Ω)≤Csta∥u0∥L1(Ω) and ∥u0h∥≤Csta∥u0∥

and

 (29) v0h=Qhv0,

which satisfies

 (30) v0h≥0 a. e. in Ω,∥v0h∥H1(Ω)≤C∥v0∥H1(Ω) and ∥˜Δhv0h∥≤C∥Δv0∥.

Given , we let be a uniform partitioning of [0,T] with time step . Given , find such that

 (31) (δtun+1h,xh)h+(∇un+1h,∇xh)=(∇vnh,un+1h∇xh)

and

 (32) (δtvn+1h,xh)h+(∇vn+1h,∇xh)+(vn+1h,xh)h=(un+1h,xh)h.

It should be noted that scheme (31)-(32) combines a finite element method together a mass-lumping technique to treat some terms and a semi-implicit time integrator. The resulting scheme is linear and decouples the computation of .

In order to carry out our numerical analysis we must rewrite the chemotaxis term by using a barycentric quadrature rule as follows. Let and consider to be the barycentre of . Then let be the interpolation of into , with being the space of all piecewise constant functions over , defined by

 (33) ¯¯¯un+1h|T=un+1h(bT).

As a result, one has

 (34) (∇vnh,un+1h∇xh)=∑T∈K|T|∇vnh⋅∇xhun+1h(bT)=(∇vnh,¯¯¯un+1h∇xh).

From now on we will use to denote a generic constant independent of the discretization parameters .

Let us define

 (35) E0(uh,vh)=12∥vh∥2h+12∥∇vh∥2−(uh,vh)h+(loguh,uh)h,
 (36) E1(uh,vh)=∥uh∥2h+∥˜Δhvh∥2.

and, for each ,

 (37) Rε,δ0(u0h,v0h):=1δe+∥u0h∥L1(Ω)δ(CΩε+ε+(1+δ)|Ω|(∥v0h∥L1(Ω)+∥u0h∥L1(Ω))).

Associated with the above definitions, consider

 B0(u0h,v0h)=1δE0(u0h,v0h)+Rδ,ε0(u0h,v0h),
 B1(u0h,v0h)=(1+1δ)E0(u0h,v0h)+Rδ,ε0(u0h,v0h)+2|Ω|e,

and

 B2(u0h,v0h)=E0(u0h,v0h)+B0(u0h,v0h)+B1(u0h,v0h).

Finally, define

 F(u0h,v0h)=eB2(u0h,v0h)+T12B122(u0h,v0h)(E0(u0h,v0h)+CTB31(u0h,v0h)+CT∥u0h∥L1(Ω)).

The definition of the above quantities will be apparent later.

We are now prepared to state the main result of this paper.

###### Theorem 3.1.

Assume that hypotheses are satisfied. Let with and , and take and defined by (27) and (29), respectively. Assume that fulfill

 (38) Ckh2F(u0h,v0h)<12

and

 (39) −Cneg+Ch1−2pF12(u0h,v0h)<0.

Then the sequence computed via (31) and (32) satisfies the following properties, for all :

• Lower bounds:

 (40) umh(x)>0

and

 (41) vmh(x)≥0

for all ,

• -bounds:

 (42) ∥umh∥L1(Ω)=∥u0h∥L1(Ω)

and

 (43) ∥vmh∥L1(Ω)≤∥v0h∥L1(Ω)+∥u0h∥L1(Ω).
• A discrete energy law:

 (44) E0(umh,vmh)+km∑r=1(∥δtvrh∥2h+k∥A−12h(urh)∇urh−A12h(urh)∇vr−1h∥2)≤E0(u0h,v0h),

where is defined in (60).

Moreover, if we are given such that

 (45) Ch1−2pE1(u0h,v0h)≤512,

it follows that

 (46) E1(umh,vmh)+k2m∑r=1(∥∇urh∥2+∥∇˜Δhvrh∥2)≤F(u0h,v0h).

As system (31)-(32) is linear, existence follows from uniqueness. The latter is an immediate outcome of a priori bounds for .

## 4. Proof of main result

In this section we address the proof of Theorem 3.1. Rather than prove en masse the estimates in Theorem 3.1, because all of them are connected, we have divided the proof into various subsections for the sake of clarity. The final argument will be an induction procedure on relied on the semi-explicit time discretization employed in (31).

### 4.1. Lower bounds and a discrete energy law

We first demonstrate lower bounds for and, as a consequence of this, a discrete local-in-time energy law is established.

###### Lemma 4.1 (Lower bounds).

Assume that are satisfied. Let and and let

 (47) ∥˜Δhvnh∥2≤F(u0h,v0h).

Then if one chooses satisfying (38) and (39), it follows that the solution computed via (31) and (32) are lower bounded, i.e, for all ,

 (48) un+1h(x)>0

and

 (49) vn+1h(x)≥0.
###### Proof.

Since and are piecewise linear polynomial functions, it will suffice to prove that (48) and (49) hold at the nodes. To do this, let be a fixed triangle with vertices , and choose two of them, i.e. with . Then, from (6), (7), (26), and (47), we have on noting (34) that

 (50) ∫T¯¯¯¯φai∇vnh⋅∇φajdx=∫Tφai(bT)∇vnh⋅∇φajdx≤|T|∥φai∥L∞(T)∥∇vn