Original article Biomath 2 (2013), 1312071, 1–6 B f Volume ░, Number ░, 20░░ BIOMATH ISSN 1314-684X Editor–in–Chief: Roumen Anguelov B f BIOMATH h t t p : / / w w w . b i o m a t h f o r u m . o r g / b i o m a t h / i n d e x . p h p / b i o m a t h / Biomath Forum Mathematical and Numerical Analysis of a Modified Keller-Segel Model with General Diffusive Tensors Georges Chamoun∗,∗∗, Mazen Saad∗, Raafat Talhouk∗∗ *Ecole Centrale de Nantes, ED STIM, LMJL UMR6629 CNRS-UN 1, rue de la Noé, 44321 Nantes, France . Emails: georges.chamoun@ec-nantes.fr, mazen.saad@ec-nantes.fr **Université Libanaise, EDST et Faculté des sciences I, Laboratoire de Mathématiques, Hadath, Liban . Email: rtalhouk@ul.edu.lb Received: 18 September 2013, accepted: 7 December 2013, published: 23 December 2013 Abstract—This paper is devoted to the mathematical analysis of a model arising from biology, consisting of dif- fusion 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 de- generate 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 ; heteroge- neous and anisotropic diffusion; global existence of so- lutions; combined scheme. I. INTRODUCTION Chemotaxis, the directed movement of cells in re- sponse to chemical gradients, plays an important role in many biological fields, such as embrogenesis, immunol- ogy, cancer growth and wound healing. This behavior enables many living organisms to locate nutrients, avoid predators or find animals of the same species. For ex- ample, 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,{ ∂tN −∇· ( S(x) ( a(N)∇N −χ(N)∇C )) = 0, ∂tC −∇· (M(x)∇C) = g(N,C), (1) where Ω is a bounded domain in Rd; d ∈ N∗ with bound- ary ∂Ω. This system of equations is supplemented by the following boundary conditions on Σt= ∂Ω × (0,T), S(x)a(N)∇N ·η = 0, M(x)∇C ·η = 0, (2) where ν is the exterior unit normal to ∂Ω. The initial conditions on Ω, are given by, N(x, 0) = N0(x), C(x, 0) = C0(x). (3) 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) = Nh(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. This article is concerned with the global existence in time of weak solutions to the model (1), for any Citation: Georges Chamoun, Mazen Saad, Raafat Talhouk, Mathematical and Numerical Analysis of a Modified Keller-Segel Model with General Diffusive Tensors, Biomath 2 (2013), 1312071, http://dx.doi.org/10.11145/j.biomath.2013.12.071 Page 1 of 6 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2013.12.071 G Chamoun, Mathematical and numerical analysis of a modified Keller-Segel model... 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 well- known 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 convection- dominated 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 ≥ Nm. 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 Nm = 1. The main assumptions are: χ : [0, 1] 7−→ R is continuous and χ(0) = χ(1) = 0 , a : [0, 1] 7−→ R+ is continuous, a(0) = a(1) = 0 and a(s) > 0 for 0 < s < 1 . (4) A standard example for χ is, χ(N) = N(1 −N); N ∈ [0, 1] . (5) 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, g(N,C) = αN −βC; α, β ≥ 0 . (6) The permeabilities S, M: Ω −→ Md(R) where Md(R) is the set of symmetric matrices d×d, verify: Si,j ∈ L∞(Ω), Mi,j ∈ L∞(Ω), ∀i,j ∈{1, ..,d} , (7) and there exist cS ∈ R∗+ and cM ∈ R∗+ such that a.e x ∈ Ω, ∀ξ ∈ Rd, S(x)ξ · ξ ≥ cS|ξ|2, M(x)ξ · ξ ≥ cM|ξ|2 . (8) Definition 2.1: Assume that 0 ≤ N0 ≤ 1, C0 ≥ 0 and C0 ∈ L∞(Ω). A couple (N,C) is said to be a weak solution of (1) if, 0 ≤ N(x,t) ≤ 1, C(x,t) ≥ 0 a.e. in QT = Ω × [0,T], N ∈ Cw([0,T]; L2(Ω)), ∂tN ∈ L2(0,T; (H1(Ω)) ′ ) , A(N) := ∫ N 0 a(r)dr ∈ L2(0,T; H1(Ω)), C ∈ L∞(QT ) ∩L2(0,T; H1(Ω)) ∩C([0,T]; L2(Ω)) , ∂tC ∈ L2(0,T; (H1(Ω)) ′ ), and (N,C) satisfy,∫ T 0 < ∂tN,ψ1 >(H1)′,H1 dt+∫∫ QT ( S(x)a(N)∇N−S(x)χ(N)∇C ) ·∇ψ1 dxdt = 0 ,∫ T 0 < ∂tC,ψ2 >(H1)′,H1 dt+ ∫∫ QT M(x)∇C ·∇ψ2 dxdt = ∫∫ QT g(N,C)ψ2 dxdt, for all ψ1, ψ2 ∈ L2(0,T; H1(Ω)), where Cw(0,T; L 2(Ω)) denotes the space of continuous functions with values in (a closed ball of) L2(Ω) endowed with the weak topology. Our first result is the following existence Theorem of weak solutions proved by using a technique of semi- discretization in time (see [4]) for the regularized non- degenerate 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: Assume that 0 ≤ N0 ≤ 1 and C0 ≥ 0 a.e. in Ω, C0 ∈ L∞(Ω). Then, the system (1) has a global weak solution (N,C) in the sense of Definition 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 C1 and if there exists a constant C > 0 such that, (χ(N1) −χ(N2))2 ≤ C(N1 −N2)(A(N1) −A(N2)), ∀N1,N2 ∈ [0, 1] . (9) Biomath 2 (2013), 1312071, http://dx.doi.org/10.11145/j.biomath.2013.12.071 Page 2 of 6 http://dx.doi.org/10.11145/j.biomath.2013.12.071 G Chamoun, Mathematical and numerical analysis of a modified Keller-Segel model... 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 L20(Ω) = {w ∈ L2(Ω), ∫ Ω wdx = 0} and we denote by Nw the unique solution to { −∇· ( S(x)∇Nw ) = w S(x)∇Nw ·η = 0 . (10) 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 C1-functions, we can rewrite (9) as follows, (χ(N1)−χ(N2))2 = (χ′(ξ))2 a(ξ1) (N1−N2)(A(N1)−A(N2)) . Then it amounts to check the following sufficient condi- tion to prove (9), ∃C > 0; max ξ,ξ1∈]0,1[ (χ′(ξ))2 a(ξ1) ≤ C . (11) 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 Rd with d = 2 or 3. In order to discretize the problem (1), we consider a family Th of meshes of the domain Ω, consisting of disjoint closed simplices such that Ω̄ = ∪K∈ThK̄ and such that if K,L ∈Th, K 6= L, then K ∩L is either an empty set or a common face or edge of K and L. We denote by Eh the set of all sides, by Einth the set of all interior sides, by Eexth the set of all exterior sides. The size of the mesh Th is defined by h:= size(Th)=maxK∈Th diam(K), which has the sense of an upper bound for the maximum diameter of the control volumes in Th. 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 {Th}h: There exists a positive constant kT ; min K∈Th |K| (diam(K))d ≥ kT , ∀h > 0 . (12) We also use a dual partition Dh of disjoint closed simplices called control volumes of Ω such that Ω̄ = ∪D∈DhD̄. There is one dual element D associated with each side σD = σK,L ∈ Eh. We construct it by connecting the barycenters of every K ∈Th that contains σD through the vertices of σD. As for the primal mesh, we define Fh, Finth and F ext h respectively as the set of all dual, interior and exterior mesh sides. For σD ∈Fexth , the contour of D is completed by the side σD itself. We refer to the Figure 1 for the two-dimensional case. Fig. 1. Triangles K,L ∈ Th and dual volumes D, E ∈ Dh associated with edges σD, σE ∈Eh. We use the following notations: • |D|= meas(D)= d-dimensional Lebesgue measure of D and |σ| is the (d−1)-dimensional measure of σ. • PD is the barycenter of the side σD. • N(D) is the set of neighbors of the volume D. • Dinth and D ext h are respectively the set of all interior and boundary dual volumes. The time discretization is the sequence of discrete times tn = n∆t for n ∈ N, where ∆t > 0 is a given time-step. In the sequel, the exponent n denotes the approximation at time tn. Next, we define the following finite-dimensional spaces: Xh := {ϕh ∈ L2(Ω); ϕh|K is linear ∀K ∈Th, ϕh is continuous at the points PD, D ∈Dinth } , X0h := {ϕh ∈ Xh; ϕh(PD) = 0, ∀D ∈D ext h } . The basis of Xh is spanned by the shape functions ϕD, D ∈ Dh, such that ϕD(PE) = δDE, E ∈ Dh, δ being the Kronecker delta. We recall that the approximations Biomath 2 (2013), 1312071, http://dx.doi.org/10.11145/j.biomath.2013.12.071 Page 3 of 6 http://dx.doi.org/10.11145/j.biomath.2013.12.071 G Chamoun, Mathematical and numerical analysis of a modified Keller-Segel model... in these spaces are nonconforming since Xh 6⊂ H1(Ω). 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 Xh with the seminorm, ||Nh||2Xh := ∑ K∈Th ∫ K |∇Nh|2dx, (13) which becomes a norm on X0h. The approximation of the flux S(x)∇C · ηD,E on the interface σD,E is denoted by δCD,E. Now, we have to approximate S(x)χ(N)∇C · ηD,E by means of the values ND,NE and δCD,E that are available in the neighborhood of the interface σD,E. To do this, we use a numerical flux function G(ND,NE,δCD,E). One can find in [1] a way to construct this numerical flux G. Finally, a combined finite volume-nonconforming fi- nite element scheme for the discretization of the problem (1) is given by the following set of equations: For all D ∈ Dh, N0D = 1 |D| ∫ D N0(x)dx, C 0 D = 1 |D| ∫ D C0(x)dx, (14) and for all D ∈Dh, n ∈{0, 1, ...,N}, |D| Nn+1D −N n D ∆t − ∑ E∈Dh SD,EA(Nn+1E ) + ∑ E∈N(D) G(Nn+1D ,N n+1 E ; δC n+1 D,E ) = 0 , (15) |D| Cn+1D −C n D ∆t − ∑ E∈Dh MD,ECn+1E = |D|g(N n D,C n+1 D ) . (16) The matrix S (resp. M), of elements SD,E (resp. MD,E) for D,E ∈ Dinth , is the diffusion matrix which is the stiffness matrix of the nonconforming finite element method. So that: SD,E = − ∑ K∈Th (S(x)∇ϕE,∇ϕD)0,K . (17) Moreover, δCn+1D,E denotes the approximation of S(x)∇C ·ηD,E on the interface σD,E: δCn+1D,E = SD,E(C n+1 E −C n+1 D ) . (18) Definition 3.1: Using the values of Nn+1D ,∀D ∈ Dh 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: Ñh,∆t(x, 0) = N 0 D for x ∈ D, D ∈Dh , Ñh,∆t(x,t) = N n+1 D for x ∈ D, D ∈Dh, t ∈]tn, tn+1] . 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: SD,E ≥ 0 and MD,E ≥ 0, ∀D ∈Dh, E ∈N(D) . (19) Recall that 0 ≤ N0 ≤ 1 and C0 ≥ 0 a.e. in Ω. We have the following classical Proposition proved by a simple induction argument as in [3]. Proposition 3.1: (Discrete maximum principle) Under the assumption (19), one has, 0 ≤ Nn+1D ≤ 1, C n+1 D ≥ 0, ∀D ∈Dh, n ∈{0, 1, ..,N}. Theorem 3.1: (Convergence of the scheme) Assume (4) to (8). Consider C0 ∈ L∞(Ω), C0 ≥ 0 and 0 ≤ N0 ≤ 1 a.e. on Ω. Under the assumption (19), 1) There exists a solution (Nh,Ch) of the discrete system (15)-(16) with initial data (14). 2) Any sequence (hm)m decreasing to zero possesses a subsequence such that (Nhm,Chm ) converges a.e. on QT to a weak solution (N,C) of (1). Outline of the proof. Adding to Proposition 3.1, we establish estimates on the discrete gradient of the func- tion 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 Xh defined in (13). Then, con- structing estimates in time and in space imply that the se- quence ( 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 L2(QT ). This implies the existence of subsequences of ( A(Ñh,∆t) ) h,∆t which converges strongly to some function Γ ∈ L2(QT ). Using the monotonicity of A, we get Γ = A(N). Moreover, due to the space translate estimate, ( [6], Theorem 3.10 ) gives that A(N) ∈ L2(0,T; H1(Ω)). 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 L2(QT ) and a.e. in QT 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). Biomath 2 (2013), 1312071, http://dx.doi.org/10.11145/j.biomath.2013.12.071 Page 4 of 6 http://dx.doi.org/10.11145/j.biomath.2013.12.071 G Chamoun, Mathematical and numerical analysis of a modified Keller-Segel model... Fig. 2. Initial condition for the cell density N0 (left) and for the concentration of the chemo-attractant C0 (right) on the dual mesh Dh 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 Cn+1 solution of the linear system given by the equation of (16) and next we compute Nn+1 as the solution of the nonlinear system defined by the first equation of (15). To this end, a Newton algorithm is implemented to ap- proach 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 Dh of a refined initial admissible mesh Th, where the assumption (19) is satisfied. In this test, we consider the following anisotropic diffusion tensors as, S = M = [ 8 −7 −7 20 ] . Further, dt = 0.0005, α = 0.01, β = 0.05, A(N) = D(N 2 2 − N 3 3 ) with D = 0.005, χ(N) = cN(1 − N)2 with c = 0.001. Finally, the diffusion coefficient of the chemo-attractant is D1 = 10−6. The initial conditions are defined by region. The initial density is defined as N0(x,y) = 1 in the square (x,y) ∈ ( [0.45, 0.55] × [0.45, 0.55] ) and 0 otherwise. The initial chemoattrac- tant is defined as C0(x,y) = 5 in the union of four squares (x,y) ∈ ( [0.2, 0.3] × [0.7, 0.8] ) ∪ ( [0.2, 0.3] × [0.2, 0.3] ) ∪ ( [0.7, 0.8]×[0.2, 0.3] ) ∪ ( [0.7, 0.8]×[0.7, 0.8] ) 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 Fig. 3. Test 1- The cell density (N), at time t = 1 with 0 ≤ N ≤ 0.1752 (left), at time t = 5 with 0 ≤ N ≤ 0.08036 (right) . Fig. 4. Test 1- The cell density (N), at time t = 10 with 0 ≤ N ≤ 0.1456 (left), at time t = 15 with 0 ≤ N ≤ 0.1666, at time t = 20 with 0 ≤ N ≤ 0.1626 (right) . 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. 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 sim- ulations, 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. ACKNOWLEDGMENT The authors would like to thank the National Council for Scientific Research (Lebanon), the Ecole Centrale de Nantes and the Lebanese University for their support to this work. Biomath 2 (2013), 1312071, http://dx.doi.org/10.11145/j.biomath.2013.12.071 Page 5 of 6 http://dx.doi.org/10.11145/j.biomath.2013.12.071 G Chamoun, Mathematical and numerical analysis of a modified Keller-Segel model... REFERENCES [1] B. Andreianov, M. Bendahmane and M. Saad, Finite volume methods for degenerate chemotaxis model. Journal of com- putational and applied mathematics, 235: p. 4015-4031, 2011. http://dx.doi.org/10.1016/j.cam.2011.02.023 [2] P. Angot, V. Dolejsi, M. Feistauer and J. Felcman, Analysis of a combined barycentric finite volume-nonconforming finite element method for nonlinear convection-diffusion problems. Appl.Math., 43(4), p. 263-310, 1998. http://dx.doi.org/10.1023/A:1023217905340 [3] R. Eymard, D. Hilhorst and M. Vohralik, A combined finite volume-nonconforming/mixed hybrid finite element scheme for degenerate parabolic problems. Numer.Math., 105: p. 73-131, 2006. [4] M. Bendahmane, K. Karlsen and J.M. Urbano, On a two-sidedly degenerate chemotaxis model with volume-filling effect. Math. Methods Appl. Sci., 17(5): p. 783-804, 2007. [5] F. Marpeau and M. Saad, Mathematical analysis of radionuclides displacement in porous media with nonlinear adsorption. J. Differential Equations 228, p. 412-439, 2006. http://dx.doi.org/10.1016/j.jde.2006.03.023 [6] R. Eymard, T. Gallouet and R. Herbin, Finite volume methods, Handbook of numerical analysis. Handbook of numerical anal- ysis, Vol VII North-Holland, Amsterdam, p. 713-1020, 2000. [7] I. Tuval, L. Cisneros, C. Dombrowski, C.W. Wolgemuth, J.O.Kessler and R.E.Goldstein, Bacterial swimming and oxygen transport near contact lines. Proc. Natl. Acad. Sci. USA, 102, pp. 2277-2282, 2005. http://dx.doi.org/10.1073/pnas.0406724102 [8] P. Laurencot et D. Wrzosek, A Chemotaxis model with thresh- old density and degenerate diffusion. Nonlinear differential equations and their applications, vol. 64, Birkhauser, Boston, p. 273-290, 2005. [9] K. Painter and T. Hillen, Volume filling effect and quorum-sensing in models for chemosensitive movement. Canadian Appl. Math. Q. 10, p. 501-543, 2002. Biomath 2 (2013), 1312071, http://dx.doi.org/10.11145/j.biomath.2013.12.071 Page 6 of 6 http://dx.doi.org/10.1016/j.cam.2011.02.023 http://dx.doi.org/10.1023/A:1023217905340 http://dx.doi.org/10.1016/j.jde.2006.03.023 http://dx.doi.org/10.1073/pnas.0406724102 http://dx.doi.org/10.11145/j.biomath.2013.12.071 Introduction Setting of the problem Numerical Scheme Numerical test conclusion References