# A corrected decoupled scheme for chemotaxis models

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.

## Authors

• 1 publication
• 2 publications
• ### Finite Volume approximation of a two-phase two fluxes degenerate Cahn-Hilliard model

We study a time implicit Finite Volume scheme for degenerate Cahn-Hillia...
05/04/2020 ∙ by Clément Cancès, et al. ∙ 0

• ### Discontinuous Galerkin methods for semilinear elliptic boundary value problem

A discontinuous Galerkin (DG) scheme for solving semilinear elliptic pro...
01/26/2021 ∙ by Jiajun Zhan, et al. ∙ 0

• ### A convergent finite volume scheme for the stochastic barotropic compressible Euler equations

In this paper, we analyze a semi-discrete finite volume scheme for the t...
08/27/2021 ∙ by Abhishek Chaudhary, et al. ∙ 0

• ### A Positive and Energy Stable Numerical Scheme for the Poisson-Nernst-Planck-Cahn-Hilliard Equations with Steric Interactions

We consider numerical methods for the Poisson-Nernst-Planck-Cahn-Hilliar...
02/22/2020 ∙ by Yiran Qian, et al. ∙ 0

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

Bresse-Timoshenko beam model with thermal, mass diffusion and theormoela...
11/07/2020 ∙ by Mohammad Elhindi, et al. ∙ 0

• ### A computationally and cognitively plausible model of supervised and unsupervised learning

Both empirical and mathematical demonstrations of the importance of chan...
10/11/2020 ∙ by David M. W. Powers, et al. ∙ 0

• ### Existence, Uniqueness and Numerical Modeling of Wine Fermentation Based on Integro-Differential Equations

Predictive modeling is the key factor for saving time and resources with...
04/14/2021 ∙ by Christina Schenk, et al. ∙ 0

##### 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

 ⎧⎨⎩∂tu=μΔu−a∇⋅(u∇c),∂tc=Δc+uu+1−c, (1.1)

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

 (1.2)

with the homogeneous Neumann boundary and initial conditions

 ∇u⋅ν=∇c⋅ν=0on ∂Ω×(0,Tf),u(.,0)=u0in Ω, (1.3)

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

 dσ={\rm d(xK,σ)if % σ⊂∂Ω,\rm d(xK,xL)otherwise , σ=K|L,

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

 τσ=\rm m(σ)dσ,σ∈E.

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

 δ=max(Δt,h).

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

 1−2aΔt≥α. (2.1)

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

 \rm d(xK,σ)≥ξ\rm d(xK,xL), (2.2)

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 :

 ∫T0∫Ω(u∂tψ−μ∇u⋅∇ψ+au∇c⋅∇ψ)dxdt+∫Ωu0ψ(x,0)dx=0, (2.3) ∫T0∫Ω∇c⋅∇ψ dxdt=∫T0∫Ω(uu+1−c)ψdxdt. (2.4)

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

 ∥v∥p=(∑K∈T\rm m(K)|vK|p)1/p,

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

 |v|1,p,T=(∑σ∈E\rm m(σ)dp−1σ|Dσv|p)1/p,
 ∥v∥1,p,T=∥v∥p+|v|1,p,T,

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

 ∥v∥r≤Cξ(p−1)θ/p∥v∥θ1,p,T∥v∥1−θs, ∀v∈X(T),

where is defined by

 θ=1/s−1/r1/s+1/d−1/p.

### 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

 \rm m(K)un+1K−unKΔt−μ∑σ∈EKτσDun+1K,σ +a∑σ∈EKσ=K|Lτσ(S(Dcn+1K,σ)un+1K−S(−Dcn+1K,σ)un+1L)=0, (2.5) −∑σ∈EKτσDcn+1K,σ+\rm m(K)cn+1K=\rm m(K)unKunK+1, (2.6)

with the compatible initial condition

 u0K=1\rm m(K)∫Ku0(x)dx, (2.7)

and where for all

 DvnK,σ={0for σ⊂∂Ω,vnL−vnKotherwise , σ=K|L.

In the above scheme, the function is defined by

 S(x)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩0,if x<2(−μ+ε)/a,x2,if |x|≤2(μ−ε)/a,x,if ,x>2(μ−ε)/a, (2.8)

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

 −∑σ∈EKτσDcn+1K,σ+\rm m% (K)cn+1K=\rm m(K)un+1Kun+1K+1, (2.9)

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

 Tn+1K=\rm m(K)(un+1Kun+1K+1−unKunK+1),T0K=0. (2.10)

The equation (2.9) can then be written as

 −∑σ∈EKτσDcn+1K,σ+\rm m% (K)cn+1K=\rm m(K)unKunK+1+Tn+1K. (2.11)

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:

 −∑σ∈EKτσDcn+1K,σ+\rm m% (K)cn+1K=\rm m(K)unKunK+1+βnTnK, (2.12)

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:

 βn=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩1,if T∗n=∅,minK∈T∗nunKunK+1un−1Kun−1K+1−unKunK+1,otherwise , (2.13)

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:

 uδ(x,t)=un+1K,cδ(x,t)=cn+1K,x∈K,t∈[tn,tn+1). (2.14)

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

 dvδ(x,t)=\rm m(σ)\rm m(MK,σ)Dvn+1K,σνK,σ,x∈MK,σ, t∈(tn,tn+1),

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 :

 un+1K≥0andcn+1K≥0for all K∈T,n=0,...,N−1, (3.1) ∑K∈T\rm m(K)un+1K=∑K∈T\rm m(K)u0K=∥u0∥L1(Ω),% for all n=0,...,N−1. (3.2)
###### Proof.

For all

, we define the vectors

, , and , for which:

 Un+1K=un+1K,FnK=\rm m(K)unK/Δt, Cn+1K=cn+1K,GnK=\rm m(K)unKunK+1+βnTnK.

We also define the matrix and by

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

 AnUn+1=Fn, (3.3) BCn+1=Gn. (3.4)

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

 |BK,K|−∑L∈TL≠K|BK,L|=\rm m(K)>0,

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

 |AnK,K|−∑L∈TL≠K|AnL,K|=\rm m(K)/Δt>0.

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:

 ∑K∈T\rm m(K)(un+1K−unK)=0.

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

 ∑K∈T\rm m(K)up+1K=∑K∈T%m(K)u0K=∥u0∥L1(Ω),for all p=0,...,N−1,

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

 cn+1K≤2, (4.1) ∑K∈T∑σ∈EKτσ∣∣Dcn+1K,σ∣∣2≤C. (4.2)

where is a positive constant which only depends on .

###### Proof.

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

 −∑σ∈EKτσDcn+1K,σ(cn+1K−2)+ ≤\rm m(K)(2−cn+1K)(cn+1K−2)+,

it follows that

 −∑σ∈EKτσDcn+1K,σ(cn+1K−2)+≤0.

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

 −∑σ∈EKτσDcn+1K,σ(cn+1K−2)+≥0.

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

 12∑K∈T∑σ∈EKτσ∣∣Dcn+1K,σ∣∣2+∑K∈T\rm m(K)∣∣cn+1K∣∣2≤4\rm m(Ω),

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

 N∑n=0∑K∈TΔt \rm m(K)∣∣unK∣∣2≤C. (4.3)
###### Proof.

Performing the following changes of variable:

 ~unK=unK+1for all K∈T,n=0,...,N, (4.4)

and

 ~S(x)=S(x)+(μ−ε)/a, (4.5)

and using the identity

 S(x)−S(−x)=x, (4.6)

the scheme (2.5) becomes

 \rm m(K)~un+1K−~unKΔt−ε∑σ∈EKτσD~un+1K,σ +a∑σ∈EKσ=K|Lτσ(~S(Dcn+1K,σ)~un+1K−~S(−Dcn+1K,σ)~un+1L)−a∑σ∈EKτσDcnK,σ=0. (4.7)

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

 ∑σ∈EKτσDcnK,σ≤3\rm m% (K).

Using the above inequality in (4.7), we obtain

 \rm m(K)~un+1K−~unKΔt−ε∑σ∈EKτσD~un+1K,σ ≤−a∑σ∈EKσ=K|Lτσ(~S(Dcn+1K,σ)~un+1K−~S(−Dcn+1K,σ)~un+1L)+3a\rm m(K).

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

 E1+E2≤E3+E4, (4.8)

where

 E1=∑K∈T\rm m(K)(~un+1K−~unK)log(~un+1K), E2=−ε∑K∈T∑σ∈EKΔtτσD~un+1K,σlog(~un+1K), E3=−a∑σ∈EKσ=K|LΔtτσ(~S(Dcn+1K,σ)~un+1K−~S(−Dcn+1K,σ)~un+1L)log(~un+1K),

From the expression of , we can see that

 E1= ∑K∈T\rm m(K)(~un+1Klog(~un+1K)−~unKlog(~unK)) −∑K∈T\rm m(K)~unK(log(~un+1K)−log(~unK)),

and since we have from a Taylor expansion of

 ~unK(log(~un+1K)−log(~unK))≤~un+1K−~unK,

we deduce using the mass conservation property (3.2) that

 E1≥∑K∈T\rm m(K)(~un+1Klog(~un+1K)−~unKlog(~unK)).

By an integration by parts on , we get

 E2=ε2∑K∈T∑σ∈EKΔtτσD~un+1K,σ(log(~un+1L)−log(~un+1K)).

Then, by a Taylor expansion of we get

 E2=ε2∑K∈T∑σ∈EKΔtτσ∣∣ ∣∣D~un+1K,σ√θn+1σ∣∣ ∣∣2,

where with . Hence, using the fact that

 ∣∣D~un+1K,σ∣∣√θn+1σ=√~un+1K+√~un+1L√θn+1σ∣∣∣D(√~un+1)K,σ∣∣∣≥∣∣∣D(√~un+1)K,σ∣∣∣,

we infer that

 E2≥ε2∑K∈T∑σ∈EKΔtτσ∣∣∣D(√~un+1)K,σ∣∣∣2.

Using a summation by parts, we obtain

 En+13=a2∑K∈T∑σ∈EKσ=K|LΔt τσ(~S(Dcn+1K,σ)~un+1K−~S(−Dcn+1K,σ)~un+1L) ×(log(~un+1L)−log(~un+1K)).

Now, we define

 ~En+13=a2∑K∈T∑σ∈EKσ=K|LΔt τσθn+1σDcn+1K,σ(log(~un+1L)−log(~un+1K)).

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

 En+13−~En+13= a2∑K∈T∑σ∈EKσ=K|LΔtτσ~S(Dcn+1K,σ)(~un+1K−θn+1σ) ×(log(~un+1L)−log(~un+1K))+a2∑K∈T∑σ∈EKσ=K|LΔtτσ~S(−Dcn+1K,σ) = a(1−tσ)2∑K∈T∑σ∈EKσ=K|LΔtτσ~S(Dcn+1K,σ)(~un+1K−~un+1L) ×(log(~un+1L)−log(~un+1K))+atσ2∑K∈T∑σ∈EKσ=K|LΔtτσ~S(−Dcn+1K,σ) ×(~un+1K−~un+1L)(log(~un+1L)−log(~un+1K)).

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

 En+13≤a2∑K∈T∑σ∈EKΔtτσDcn+1K,σD~un+1K,σ.

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

 En+13≤2aΔt∥u0+1∥L1(Ω).

Finally for , we have

 E4≤3aΔt∥u0+