Mathematical and Numerical Analysis of a Modified Keller-Segel Model with General Diffusive Tensors

This paper is devoted to the mathematical analysis of a model arising from biology, consisting of diffusion and chemotaxis with volume filling effect. Motivated by numerical and modeling issues, the global existence in time and the uniqueness of weak solutions to this model is investigated. The novelty with respect to other related papers lies in the presence of a two-sidedly nonlinear degenerate diffusion and anisotropic heterogeneous diffusion tensors, where we prove global existence and uniqueness under further assumptions. Moreover, we introduce and we study the convergence analysis of the combined scheme applied to this anisotropic Keller-Segel model with general tensors. Finally, a numerical test is given to prove the effectiveness of the combined scheme. Keywords-Degenerate parabolic equation ; heterogeneous and anisotropic diffusion; global existence of solutions; combined scheme.


I. INTRODUCTION
Chemotaxis, the directed movement of cells in response to chemical gradients, plays an important role in many biological fields, such as embrogenesis, immunology, cancer growth and wound healing.This behavior enables many living organisms to locate nutrients, avoid predators or find animals of the same species.For example, bacteria often swim towards higher concentration of oxygen to survive (see [7]).In the following, we investigate a system consisting of the parabolic Keller-Segel equations with general tensors, where Ω is a bounded domain in R d ; d ∈ N * with boundary ∂Ω.This system of equations is supplemented by the following boundary conditions on Σ t = ∂Ω × (0, T ), where ν is the exterior unit normal to ∂Ω.The initial conditions on Ω, are given by, N (x, 0) = N 0 (x), C(x, 0) = C 0 (x).
Here N and C denote respectively the cell density and the concentration of a chemical.Moreover, a(N ) denotes the density-dependent diffusion coefficient and χ(N ) is usually written in the form χ(N ) = N h(N ) where h is commonly referred as a chemotactic sensitivity function.S(x) and M (x) are considered as general tensors which may be anisotropic and heterogeneous.
In the model (1), the cell density N diffuses and swims upwards chemical gradients.In addition to that, the chemical C also diffuses where g(N, C) is the rate of production and consumption of the chemical.
space dimension, in the presence of anisotropic and heterogeneous tensors, two-sidedly nonlinear degenerate diffusion and modified chemotactic sensitivity χ.The techniques of the proof of global existence are essentially those designed by [4] and [5].Under further assumptions and in the spirit of the duality method used in [8], the uniqueness of these weak solutions is guaranteed.Moreover, in order to discretize our model (1), it is wellknown that classical finite volume methods not permit to handle anisotropic tensors in the diffusive terms but it is very useful, especially the technique "upwind", to discretize the convective term since it does not allow instabilities in the numerical solution.The intuitive idea is hence to combine two numerical methods (see [2], [3]) by discretizing the diffusive terms with the finite element method enabling the discretization of anisotropic diffusion tensors without any restrictions on the meshes and the other terms with the finite volume method since we avoid numerical instabilities in the convectiondominated regime.Finally, a numerical test will be given to illustrate the effectiveness of our combined scheme applied to the anisotropic Keller-Segel model (1).

II. SETTING OF THE PROBLEM
Let us now state the assumptions on the data we will use in the sequel, together with the main results obtained in this paper.We assume that χ(0) = 0 and the chemotactical sensitivity χ(N ) vanishes when N ≥ N m .This threshold condition has a clear biological interpretation called the volume-filling effect.In fact, the effect of a threshold cell density or a volume filling effect has been also taken into account in the modeling of chemotaxis phenomenon in [9].Upon normalization, we can assume that the threshold density is N m = 1.The main assumptions are: A standard example for χ is, The positivity of χ means that the chemical attracts the cells; the repellent case is the one of negative χ.In addition to that, will assume that the rate g is linear, The permeabilities S, M : ) is the set of symmetric matrices d × d, verify: and there exist c , where C w (0, T ; L 2 (Ω)) denotes the space of continuous functions with values in (a closed ball of) L 2 (Ω) endowed with the weak topology.
Our first result is the following existence Theorem of weak solutions proved by using a technique of semidiscretization in time (see [4]) for the regularized nondegenerate problem associated to (1).Next, when the parameter of regularization tends to zero, a similar strategy as in [5] will be followed to achieve the proof.
Theorem 2.1: Our second result is the following Theorem.Theorem 2.2: Under the assumptions of Theorem 2.1, and assume that the functions A and χ are of class C 1 and if there exists a constant C > 0 such that, Biomath 2 (2013), 1312071, http://dx.doi.org/10.11145/j.biomath.2013.12.071 Then the system (1) has a global unique weak solution.
Outline of the proof: It relies on a classical duality technique.We introduce the subset L 2 0 (Ω) = {w ∈ L 2 (Ω), Ω w dx = 0} and we denote by N w the unique solution to This dual problem (10) and similar guidelines of [8] (subsection 2.2) are followed to apply Gronwall Lemma and to prove the uniqueness statement of Theorem 2.2.
In fact, we can give a useful sufficient condition which implies (9).Indeed, by a straight application of the mean theorem applied to χ and A which are C 1 -functions, we can rewrite (9) as follows, Then it amounts to check the following sufficient condition to prove (9), This sufficient condition is also mentioned in [8], Proposition 4 for diffusion coefficients with one point degeneracy.
Example.In the model (1), if a(u) = u(1 − u) and χ(u) = (u(1−u)) β then the weak solution of the system (1) is unique provided β ≥ 3 2 .A simple computation is left to the reader to check the sufficient condition (11).

III. NUMERICAL SCHEME
This section is devoted to the formulation of a combined scheme for the anisotropic Keller-Segel model (1).First, we will describe the space and time discretizations, define the approximation spaces and then we will introduce the combined scheme.
Let Ω be an open bounded subset of R d with d = 2 or 3.In order to discretize the problem (1), we consider a family T h of meshes of the domain Ω, consisting of disjoint closed simplices such that Ω = ∪ K∈Th K and such that if K, L ∈ T h , K = L, then K ∩ L is either an empty set or a common face or edge of K and L. We denote by E h the set of all sides, by E int h the set of all interior sides, by E ext h the set of all exterior sides.The size of the mesh T h is defined by h:= size(T h )=max K∈Th diam(K), which has the sense of an upper bound for the maximum diameter of the control volumes in T h .We also define a geometrical factor, linked with the regularity of the mesh, by making the following shape regularity assumption on the family of triangulations {T h } h : There exists a positive constant k T ; We also use a dual partition D h of disjoint closed simplices called control volumes of Ω such that Ω = ∪ D∈Dh D. There is one dual element D associated with each side σ D = σ K,L ∈ E h .We construct it by connecting the barycenters of every K ∈ T h that contains σ D through the vertices of σ D .As for the primal mesh, we define F h , F int h and F ext h respectively as the set of all dual, interior and exterior mesh sides.For σ D ∈ F ext h , the contour of D is completed by the side σ D itself.We refer to the Figure 1 for the two-dimensional case.The time discretization is the sequence of discrete times t n = n∆t for n ∈ N, where ∆t > 0 is a given time-step.In the sequel, the exponent n denotes the approximation at time t n .Next, we define the following finite-dimensional spaces: The basis of X h is spanned by the shape functions ϕ D , D ∈ D h , such that ϕ D (P E ) = δ DE , E ∈ D h , δ being the Kronecker delta.We recall that the approximations in these spaces are nonconforming since X h ⊂ H 1 (Ω).Indeed, only the weak continuity of the solution is provided through the faces of the mesh and therefore the solution may be discontinuous on the faces.We equip X h with the seminorm, which becomes a norm on X 0 h .
The approximation of the flux S(x)∇C • η D,E on the interface σ D,E is denoted by δC D,E .Now, we have to approximate S(x)χ(N )∇C • η D,E by means of the values N D , N E and δC D,E that are available in the neighborhood of the interface σ D,E .To do this, we use a numerical flux function G(N D , N E , δC D,E ).One can find in [1] a way to construct this numerical flux G.
Finally, a combined finite volume-nonconforming finite element scheme for the discretization of the problem (1) is given by the following set of equations: For all D ∈ D h , and for all D ∈ D h , n ∈ {0, 1, ..., N }, (16) The matrix S (resp.M), of elements S D,E (resp.M D,E ) for D, E ∈ D int h , is the diffusion matrix which is the stiffness matrix of the nonconforming finite element method.So that: denotes the approximation of S(x)∇C • η D,E on the interface σ D,E : Definition 3.1: Using the values of N n+1 D , ∀D ∈ D h and n ∈ [0, N ], we will define the approximate finite volume solution Ñh,∆t defined as piecewise constant on the dual volumes in space and piecewise constant in time, such that: Now, we state the discrete maximum principle and a convergence result of the combined scheme under the assumption that all transmissibilities coefficients, defined in (17), are positive: Recall that 0 ≤ N 0 ≤ 1 and C 0 ≥ 0 a.e. in Ω.We have the following classical Proposition proved by a simple induction argument as in [3].
Outline of the proof.Adding to Proposition 3.1, we establish estimates on the discrete gradient of the function A( Ñh,∆t ) and Ñh,∆t and these discrete properties on the diffusive term allow to construct a priori estimate with respect to the norm X h defined in (13).Then, constructing estimates in time and in space imply that the sequence A( Ñh,∆t ) h,∆t satisfies the assumptions of the Kolmogorov's compactness criterion, and consequently A( Ñh,∆t ) h,∆t is relatively compact in L 2 (Q T ).This implies the existence of subsequences of A( Ñh,∆t ) h,∆t which converges strongly to some function Γ ∈ L 2 (Q T ).Using the monotonicity of A, we get Γ = A(N ).Moreover, due to the space translate estimate, [6], Theorem 3.10 gives that A(N ) ∈ L 2 (0, T ; H 1 (Ω)).As A −1 is well defined and continuous, applying the L ∞ bound on Ñh,∆t and the dominated convergence theorem of Lebesgue to Ñh,∆t = A −1 (A( Ñh,∆t )), there exists a subsequence of Ñh,∆t which converges strongly in L 2 (Q T ) and a.e. in Q T to the same function N .Then, we conclude by showing that the limit couple (N, C) is a weak solution of the continuous problem (1).IV.NUMERICAL TEST Through this numerical example, we would like to illustrate the effectiveness of the combined scheme for anisotropic Keller-Segel model (1).All the computations for this new numerical scheme have been implemented using the software package Fortran.The algorithm used to compute numerical solution of the discrete problem is the following: at each time step, we first calculate C n+1 solution of the linear system given by the equation of ( 16) and next we compute N n+1 as the solution of the nonlinear system defined by the first equation of (15).To this end, a Newton algorithm is implemented to approach the solution of nonlinear system and a bigradient conjugate method to solve linear systems arising from the Newton algorithm process.We will provide this test on the dual mesh D h of a refined initial admissible mesh T h , where the assumption (19) is satisfied.
In this test, we consider the following anisotropic diffusion tensors as, 3 ) with D = 0.005, χ(N ) = cN (1 − N ) 2 with c = 0.001.Finally, the diffusion coefficient of the chemo-attractant is D 1 = 10 −6 .The initial conditions are defined by region.The initial density is defined as N 0 (x, y) = 1 in the square (x, y) ∈ [0.45, 0.55] × [0.45, 0.55] and 0 otherwise.The initial chemoattractant is defined as C 0 (x, y) = 5 in the union of four squares (x, y) and 0 otherwise (see Figure 2).In the isotropic case (S = M = Id), the cells diffuse towards the four regions of the chemoattractant.In this test case and under the influence of the anisotropic diffusion, we observe in

V. CONCLUSION
In this article, we have explored the relevance of the Keller-Segel equations in the modeling of anisotropic chemotactic cell migration.Motivated by numerical simulations, we have proceeded to prove results pertaining to the existence and the uniqueness of global weak solutions of the anisotropic Keller-Segel model with general diffusive tensors and the convergence analysis of a new combined scheme introduced.This numerical scheme is a compromise between the nonconforming finite elements, enabling in particular the use of general meshes and the discretization of anisotropic diffusion tensors, and between the finite volumes enabling to avoid spurious oscillations in the convection-dominated regime.Finally, a numerical test was given to illustrate the combined scheme proposed.
Then, the system (1) has a global weak solution (N, C) in the sense of Definition 2.1.

Fig. 1 .
Fig. 1.Triangles K, L ∈ T h and dual volumes D, E ∈ D h associated with edges σD, σE ∈ E h .

Fig. 2 .
Fig. 2. Initial condition for the cell density N0 (left) and for the concentration of the chemo-attractant C0 (right) on the dual mesh D h

Figures 3
Figures 3 and 4 the movement of the cells only towards two chemoattractant regions during the stage of evolution in time of the cell density.