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.
). 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 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.
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
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
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.
, we define the vectors, , and , for which:
We also define the matrix and by
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 .
4 A priori estimates
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).
which gives (4.2).
Performing the following changes of variable:
and using the identity
the scheme (2.5) becomes
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