A corrected decoupled scheme for chemotaxis models

02/26/2020 ∙ by M. Akhmouch, et al. ∙ 0

The main purpose of this paper is to present a new corrected decoupled scheme combined with a spatial finite volume method for chemotaxis models. First, we derive the scheme for a parabolic-elliptic chemotaxis model arising in embryology. We then establish the existence and uniqueness of the numerical solution, and we prove that it converges to a corresponding weak solution for the studied model. In the last section, several numerical tests are presented by applying our approach to a number of chemotaxis systems. The obtained numerical results demonstrate the efficiency of the proposed scheme and its effectiveness to capture different forms of spatial patterns.



There are no comments yet.


page 17

page 21

page 22

page 23

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

Chemotaxis refers to a phenomenon that enables cells (or organisms) to migrate in response to a chemical signal. This process has sparked the interest of many scientists since it is encountered in several medical and biological applications, such as bacteria aggregation, tumour growth, integumental patterns in animals etc.

In Oster1 , Oster and Murray discussed a cell-chemotaxis model involving motile cells that respond to a chemoattractant secreted by the cells themselves. In its dimensionless form, the model reads


where and are positive constants, is the cell density and is the concentration of chemoattractant.

The above system is based on the Keller-Segel model Keller1 , which is the most popular model for chemotaxis. The migration of cells is assumed to be governed by Fickian diffusion and chemotaxis, and the mass of cells is conserved. The chemoattractant is assumed also to diffuse, but it increases with cell density in Michaelis-Menten way and undergoes decay through simple degradation.

In Murray1 , Murray et al. suggest that the presented cell-chemotaxis model is an appropriate mechanism for the formation of stripe patterns on the dorsal integument of embryonic and hatchling alligators (Alligator mississippiensis). These skin pigment patterns is associated with the density of melanocyte cells: Melanocytes are abundant in the regions where the black stripes appear, and are insufficient in the regions of the white stripes. The formation of these stripes is a result of a chemoattractant secretion. The system (1.1) is known to produce propagating pattern of standings peaks and troughs in cell density in the case of one-dimensional space. This patterning process was numerically and analytically investigated by Myerscough and Murray in Myerscough1 .

In this paper, we will first focus on the following parabolic-elliptic system


with the homogeneous Neumann boundary and initial conditions


where , is an open bounded polygonal subset, is a fixed time and denotes the outward unit-normal on the boundary . The system (1.2) is a simpler version of the original model (1.1) in which the second equation of the system is elliptic, using the reasonable assumption that the chemoattractant diffuses much faster than cells.

In this work, we develop a decoupled finite volume scheme which can be applied to a class of chemotaxis models. For the convection-diffusion term, the approximation used is quite similar to the hybrid scheme of Spalding Spalding1 . Concerning the time discretization, which is the main aim of this paper, it is developed such that the scheme only requires to solve decoupled systems, which excludes fully implicit discretizations. We require also that the scheme converges without needing to fulfill any CFL condition, which is not the case of the fully explicit schemes. In the literature, a number of decoupled methods for the Keller-Segel model and its variants have been proposed (see, e.g., Saito2 ; Saito1 ; Strehl1 ; Andreianov1 ; Ibrahim1 ; Chamoun1 ; Akhmouch1 ; Akhmouch2 ). In all these works, the time discretization is based on the classical backward Euler scheme with an explicit approximation of some terms to avoid coupling of the system. However, it is well known that the main drawback of this strategy is its lack of accuracy. A more efficient approach will be presented in this work.

This paper provides also a convergence analysis of the proposed scheme applied to the system (1.2)–(1.3). It is proved that the convergence of the approximate solution can be obtained for any nonnegative initial cell density . Our proof uses some techniques from Filbet1 , where a fully implicit upwind finite volume scheme is studied for the classical Keller-Segel model.

The outline of this paper is as follows. In the next section, we present our corrected decoupled finite volume scheme to approximate the solution of (1.2)–(1.3

). In section 3, we prove the existence and uniqueness of the solution of the proposed scheme. Positivity preservation and mass conservation are also shown in this section. A priori estimates are given in Section 4. In Section 5, we use these estimates to prove that the approximate solution converges to a weak solution of the studied model. In Section 6, we present some numerical tests and compare the accuracy of our approach with that of more usual decoupled schemes. The paper ends with a conclusion.

2 Presentation of the numerical scheme

2.1 Spatial discretization of , definitions and preliminaries

We assume that is an open bounded polygonal subset. Following Definition 9.1 in Eymard1 , we consider an admissible finite volume mesh of , denoted by . This mesh is given by:

  • A family of control volumes which is commonly denoted by the same notation of the mesh . All control volumes are open and convex polygons.

  • A family of edges , where the set of edges of any control volume is denoted by . We denote also by the edge between and ( and ).

  • A family of points such that (for all ). The straight line going through and must be orthogonal to .

For all , we denote by the set of control volumes which have a common edge with , and by m the Lebesgue measure in or .

For all , we define

where d is the Euclidean distance, and we denote by the transmissibility coefficient given by:

The time discretization of is given by a uniform partition: with and for .

We denote by the maximal size (diameter) of the control volumes included in , and we define

In Section 4, the following time-step condition will be used: there exists such that


We will also need this additional constraint on the mesh: there exists such that


which is specially needed to apply the discrete Gagliardo-Nirenberg-Sobolev inequality (see Lemma 2.1).

We define a weak solution of the system (1.2) with boundary and initial conditions (1.3) as follows:

Definition 2.1.

A weak solution of the initial-boundary value problem (1.2)–(1.3) is a pair of functions which verify the following identities for all test functions :


We denote by the set of functions from to which are constant over each control volume of the mesh. Let , if , the corresponding discrete norm reads

where for all and for all . We define also the discrete seminorm and the discrete norm:

where for all , if and otherwise, with .

We now recall the discrete Gagliardo-Nirenberg-Sobolev inequality (see chatard2 ), which will be useful to establish a priori estimates in Section 4.

Lemma 2.1.

Let be a an open bounded polyhedral domain of , . Let a mesh satisfying (2.2) and .

  • If , let

  • If , let

Then there exists a constant only depending on and such that

where is defined by

2.2 The corrected decoupled finite volume scheme

We begin this section by presenting a classical decoupled finite volume scheme for the problem (1.2)–(1.3):

for all and


with the compatible initial condition


and where for all

In the above scheme, the function is defined by


where is a small constant such that and .

The terms and denote respectively the approximations of the quantities and . As we can see, the proposed finite volume scheme is decoupled: at each time-step, we begin by solving (2.6) to compute and then, we compute from (2.5). The discretization used for is equivalent to the second order central difference scheme when and to the first order upwind scheme when or . When , the scheme is identical to that of Spalding Spalding1 (see also Patankar1 ).

It is clear that we can obtain a best accuracy if we replace (2.6) by the equation


however, it will be expensive in term of computational cost to find the solution of the scheme (2.5),(2.9) since we have to solve a large nonlinear system at each time step.

Now, for all and , we define


The equation (2.9) can then be written as


As we can see, the only difference between (2.11) and (2.6) is the term , so we conjecture that we can improve the accuracy of the decoupled scheme (2.5)–(2.6) if we add to the right hand side of (2.6) a correction term which approximates . Hence, we propose to replace (2.6) with the following scheme:


where . The purpose of is to ensure the nonnegativity of the right-hand side of (2.12), to do not affect the nonnegativity of . Then, by supposing that (this will be proved in Section 3), we define for as follows:


where for and . We can easily verify that . We mention however that, in practice, we will take , which seems the most natural choice. Indeed, several numerical tests are performed and, as expected, the right-hand side of (2.12) is always positive for this value unless the time-step size is extremely large.

We define , the finite volume approximations of and by:


We define also approximations of the gradients of and . To this end, we begin by defining , which is the cell formed from the vertices of , and if , and from the vertices of and if .

Following Hillairet1 , we define the discrete gradient which is the approximation of by

for all and , where denotes the unit normal on which is outward to .

3 Existence and uniqueness of a discrete solution

In this section, we prove existence and uniqueness of the solution of the proposed scheme. We show also that the scheme is mass conserving and positivity preserving.

Proposition 3.0.

Assume that . Then there exists a unique solution to the scheme (2.5),(2.12) which satisfies the following properties :


For all

, we define the vectors

, , and , for which:

We also define the matrix and by

The decoupled scheme (2.5),(2.12) can be then written equivalently as


To prove the desired results, we proceed by induction on . The proof in the case of is similar to that of the inductive step. We argue now that the vectors and are defined and nonnegative. We have

which means the matrix is strictly diagonally dominant by rows. Moreover, since the diagonal elements of are positive and its off-diagonal entries are nonpositive, we can conclude that is a nonsingular M-matrix, which implies the unique solvability of (3.4) with . Therefore, in view of the nonnegativity of (by the definition of (2.13)), we deduce that . Since for all , , we have for all

Then is strictly diagonal dominant by columns. Now, since for all (by the definition (2.8)), the matrix has nonpositive off-diagonal and positive diagonal entries, which implies that is a nonsingular M-matrix. Consequently, the existence and uniqueness of is proved, and since , it is clear that .

Now, summing (2.5) over , we have:

Then, summing over , with we get (3.2):

which ends the proof.

4 A priori estimates

In all this section, we assume that and . We assume also that is the solution of the scheme (2.5),(2.12).

Proposition 4.0.

For all and for all , we have


where is a positive constant which only depends on .


Let be the control volume which verifies . Multiplying the equation (2.12) by , and using the fact that , we get

it follows that

In view of the choice of , for all , which leads to

Hence, we deduce that . This establish (4.1).

Now, multiplying the equation (2.12) by , summing over , using a summation by parts and (4.1), we get

which gives (4.2).

Lemma 4.1.

Assume that (2.2) is fulfilled and that (see (2.8)). Then, there exists a constant only depending on , ,, , and such that


Performing the following changes of variable:




and using the identity


the scheme (2.5) becomes


From (2.12), using (4.1) and since , we can easily see that

Using the above inequality in (4.7), we obtain

Let us now multiply the above inequality by , and summing over , we find that



From the expression of , we can see that

and since we have from a Taylor expansion of

we deduce using the mass conservation property (3.2) that

By an integration by parts on , we get

Then, by a Taylor expansion of we get

where with . Hence, using the fact that

we infer that

Using a summation by parts, we obtain

Now, we define

Next, using the identity (4.6) and the last expression of , we can write

Since we have from (2.8) that , it follows that . Consequently, the above equality implies that , which yields by a Taylor expansion of log to

Multiplying (2.12) by and summing by parts we can easily see that

Finally for , we have