An unconditionally energy stable and positive upwind DG scheme for the Keller-Segel model

The well-suited discretization of the Keller-Segel equations for chemotaxis has become a very challenging problem due to the convective nature inherent to them. This paper aims to introduce a new upwind, mass-conservative, positive and energy-dissipative discontinuous Galerkin scheme for the Keller-Segel model. This approach is based on the gradient-flow structure of the equations. In addition, we show some numerical experiments in accordance with the aforementioned properties of the discretization. The numerical results obtained emphasize the really good behaviour of the approximation in the case of chemotactic collapse, where very steep gradients appear.


Introduction
Since the introduction in the 70's of the biological Keller-Segel model for chemotaxis phenomena [30,31], it and many related variants have attracted a great deal of interest in the mathematical community.Chemotaxis, a biological process through which organisms (e.g.cells) migrate in response to a chemical stimulus, is modelled by means of nonlinear systems of partial differential equations (PDE).The classical one can be written as follows: find two real valued functions, u = u(x, t) and v = v(x, t), defined in Ω × [0, T ] such that: in Ω × (0, T ), (1a) Herein Ω is a bounded and smooth domain of R d , with d ∈ N the spatial dimension, and the parameters are k i > 0 for i ∈ {0, 1, 2, 3, 4}.The mathematical formulation of (1) can be interpreted in biological terms as follows: u denotes a certain cell distribution (or population of organisms, in general) at the position x ∈ Ω and time t ∈ [0, T ], whereas v stands for the concentration of chemoattractant (i.e. a chemical signal towards which cells are induced to migrate).Both cells and chemoattractant experiment some diffusion in the spatial domain.This auto-diffusion phenomena (experimented by cells and chemoattractant) are designed by the terms −k 0 ∆u and −k 2 ∆v, while the migration mechanism is modeled by the nonlinear crossdiffusion term −k 1 ∇ • (u∇v).This term is the major difficulty for theoretical analysis and also for numerical modelling of system 1.Further, the degradation and production of chemoattractant are associated with the terms −k 3 v and k 4 u, respectively.Note that the production of chemoattractant by the cells, to which cells are attracted, may eventually result in a chemotactic collapse, a phenomenon in which uncontrolled aggregation for u give rise to blowing up or exploding in finite time.This feature is well known and constitutes one of the outstanding characteristics of classical Keller-Segel model, and also one of its main challenges, specifically for numerical methods.Finally, the coefficient τ ∈ {0, 1} is considered to write at the same time the parabolic system when τ = 1, or the parabolic-elliptic for τ = 0.
Regarding the mathematical analysis for the system (1): some results on sufficient conditions to ensure global existence and boundedness of solutions along time can be shown (see e.g. the review of Bellomo et al [7] and the references therein).They are based on mass conservation for u and on an energy dissipation law for this model (see Section 3).For dimension d ≥ 2, these results require the initial density of cells, Ω u 0 , to be bigger than certain threshold.On the other hand, considerable research has been done in the direction of finding cases where chemotactic collapse arise.Among them, it is worth mentioning the first result in this direction, due to Herrero and Velázquez [26], where radially symmetric two-dimensional solutions which finite-time blow up are found.Other authors shed light on more general cases, for instance Horstmann and Wang [27] (non symmetric blow-up solutions) or Winkler [40] (higher dimensional case).
In the last decades, a lot of papers have been published dealing with these kinds of theoretical issues both for the classical model (1) and for other models based on some extensions or generalizations.In general, they start from the Keller-Segel classical equations and modify them with the purpose of avoiding the non-physical blow up of solutions, producing solutions which are closer to the "real chemotaxis" phenomena observed in biology.Models include logistic, non-linear diffusion or production terms, chemo-repulsion effects or coupling with fluid equations [11,23,37,38,41].See e.g.[5,7] for more examples.Thus, it is hoped that understanding the classical Keller-Segel equations may open new insights for dealing in depth with those other chemotaxis models.
On the other hand, taking into account the considerable efforts of the mathematical community in the theoretical analysis of Keller-Segel models, the number of papers dedicated to numerical analysis and simulation of chemotaxis equations is much lower, in relative terms.The main difficulty is the numerical approximation of the cross-diffusion term.Not only for its non-linearity, but due to its convective nature, which makes particularly difficult to deal with using the finite element (FE) method (see, for instance, [9,20] for more details about this method).This difficulty is specifically significant in steep-gradient regions for v, which are precisely relevant in blow-up settings.Furthermore, preserving the physical properties of the continuous model (mass conservation, positivity and energy dissipation) in the discrete case adds an extra level of complication when it comes to designing a well-suited approximation.
Despite that, many interesting works have been published on numerical simulation of chemotaxis equations using different kinds of approaches.For instance, the work by Saito [35] uses FE with upwind stabilization for the parabolic-elliptic Keller-Segel model (τ = 0 in (1)) showing mass conservation, positivity and error estimates.Also, Gutiérrez-Santacreu and Rodríguez-Galván demonstrate positivity, an energy law and a priori bounds for their FE scheme on acute meshes in [25].In addition, other sophisticated techniques have been applied to this problem.This is the case of the finite volume (FV) method (we recommend [32] on this topic) which has become a very popular and successful approach as we can observe in papers like [13] by Chertock and Kurganov in which the authors devised positive preserving methods for the parabolic-parabolic formulation (τ = 1 in (1)), with demonstrated high accuracy and robustness, specifically on chemotactic collapse.It might also be pointed out the works of Saad and others, for instance in [29], where a volume finite element scheme for the capture of spatial patterns for a volume-filling chemotaxis is analyzed.
In this sense, it is also worth mentioning the very recent works [6,12,28,36] that show and analyze different techniques to approximate the solution of (1) while achieving mass conservation, positivity and energy-dissipation for a strictly positive initial cell condition.In the paper by Badía et al. [6] a discrete scheme using stabilized FE with a graph-Laplacian operator and a shock detector is proposed.This discretization satisfies, in general, both the mass conservation and positivity properties, and, in the case of acute meshes, it is also energy-dissipative.Otherwise, Huang and Shen develop in [28] a time-discrete approximation, admitting any spatial discretization, that preserves the positivity for the cell distribution u and is energy stable for a modified energy.This latest approach uses a suitable transformation of the solution for the positivity and the scalar auxiliary variable (SAV) technique for the energy stability.Alternatively, Shen and Xu introduced in [36] another general, positive (for the cell distribution u) and energy-dissipative, time-discrete scheme based on the gradient flow structure of the continuous model that admits different kinds of spatial discretization such as FE, spectral methods or even finite-differences.The order of the latest approach and the blow-up phenomenon of the discrete solution, under a CFL condition, is studied in [12] by Chen et al.
Furthermore, discontinuous Galerkin (DG) methods (we refer the reader to [15,17,34] for a further insight) aroused the interest of researchers in recent years due to their flexibility for approximating, using standard meshes and computer libraries, different types of PDEs: elliptic, parabolic, hyperbolic.In the chemotaxis context, it is worth mentioning the paper of Y. Epshteyn and A. Kurganov [19], where the FV scheme given in [13] for the 2D Keller-Segel model ( 1) is extended to a DG scheme on cartesian meshes, obtaining good approximations even on blow-up regimes.Y. Epshteyn introduced two other related schemes in [18].In all cases, different discontinuous Interior Penalty (IP) methods and upwinding techniques were considered for defining DG approximations.The schemes are even applied to the simulation of an haptotaxis model of tumor invasion into healthy tissue.Error estimates are shown but no energy property or maximum principle for the schemes is proven.In fact, spurious oscillations and negative values in the solution are reported.
More recent works make further progress in this direction, for instance in [42], where the classical equations (1) are approximated by a positivity-preserving DG method with strong stability preserving (SSP) high order time discretizations.Error order estimates as well as positivity are shown in this work.Finally, in [24,33], the local discontinuous Galerkin method is applied, showing respectively positivity and energy dissipation.
In this work, we propose a new upwind DG scheme for the Keller-Segel model (1) that preserves the mass-conservation, positivity and energy stability properties of the continuous problem.As in [36], the proposed discretization takes advantage of the gradient flow structure of the model.First, section 2 sets the notation that we are going to consider throughout the paper.In section 3 we discuss the physical properties of the continuous model.Section 4 is the main part of the paper, in which we introduce the upwind DG scheme (8).In particular, in section 4.1, we define the upwind approximation based on the ideas introduced in [1] along with some geometrical considerations for the mesh family T h , and we discuss the properties of the scheme in section 4.2.Finally, in section 5, we show several numerical tests in which we reproduce some blow-up results shown in the literature with one steady peak, [13], and a peak moving towards the corner of the domain, [35], as well as pattern formation results with several peaks, [4,10,21,39].These numerical experiments endorse the good behaviour of the approximation obtained with the new scheme, which allows us to capture peaks reaching values up to the order of 10 7 .These kinds of numerical results are rare in the literature due to the steep-gradients inherent to such sort of tests.

Notation
In this section we introduce the notation used throughout the paper.First, we consider a shaperegular triangular mesh T h = {K} K∈T h of size h over a bounded polygonal domain Ω ⊂ R d .Moreover, we note the set of edges or faces of T h by E h , which can be split into the interior edges E i h and the boundary edges . Now, we fix the following orientation for the unit normal vector n e associated to and edge e ∈ E h of the mesh T h : h is shared by the elements K and L, i.e. e = ∂K ∩ ∂L, then n e is exterior to K pointing to L (see Figure 1).
• If e ∈ E b h , then n e points outwards of the domain Ω.
Moreover, we denote the baricenter of the triangle K ∈ T h by C K .Now, we define the approximation spaces of discontinuous, P disc k (T h ), and continuous, P cont k (T h ), finite element functions over T h whose restriction to K ∈ T h are polynomials of degree k ≥ 0. In addition, the average { {•} } and the jump [[•]] of a scalar function v on an edge e ∈ E h are defined as follows: Regarding the time discretization, we take an equispaced partition 0 Finally we set the following notation for the positive and negative parts of a scalar function v:
To the best knowledge of the authors, the existence of global solutions in time is still not clear in the literature.
Remark 3.2.By taking u = 1 in (2a) (or (4a)), any solution u conserves the mass, because , and adding the resulting expressions, one has that any solution (u, v) satisfies the following energy law where E : and 4 Fully discrete scheme First, we regularize the chemical potential of u, defined in (3), by for some ε > 0. Note that ε is a regularization parameter since log(u + ε) is regular for all u ≥ 0. We will take ε = ε(h, ∆t) such that ε(h, ∆t) → 0 if (h, ∆t) → 0, being h > 0 the mesh size and ∆t > 0 the time step.Then, we propose the following decoupled fully discrete first order in time and upwind DG in space scheme for the model (1): Let v m ∈ P cont 1 (T h ) and u m ∈ P disc 0 (T h ) such that v m ≥ 0 in the case τ > 0 and u m ≥ 0 be given.
Step 1: Find for all v ∈ P cont 1 (T h ).
In order to preserve the positivity of v m+1 , we have done mass lumping in the terms ∂ t v m+1 , v h and k 3 v m+1 , v h in (8a).
In section 4.2, we will provide a way of computing the solution of (8) enforcing the nonnegativity restriction u m+1 ≥ 0.

Definition of upwind bilinear form a
Now, we are going to define the upwind bilinear form a upw h (•; •, •), introduced in the scheme (8).In order to achieve the energy stability with the scheme (8), we must consider the following hypothesis that will let us approximate the flux −∇µ accordingly: Hypothesis 1.The mesh T h of Ω is structured in the sense that the line between the baricenters of the triangles K and L is orthogonal to the interface e = K ∩ L ∈ E i h .Then, we define the following upwind bilinear form to be applied to the flux −∇µ which may be discontinuous over E i h : where, for every e ∈ E i h with e = K ∩ L, with Π 0 being the projection on P disc 0 (T h ) and D e (T h ) the distance between the baricenters of the triangles K and L of the mesh T h that share e ∈ E i h , denoted by C K and C L , respectively.This way, we can rewrite (9) as Remark 4.1.Since the quadrature formula of the baricenter (or centroid) is exact for polynomials of order 1, if µ ∈ P disc 1 (T h ), we have that where C K is the baricenter of K ∈ T h .Hence, if µ ∈ P disc 1 (T h ), the expression (10) is the slope of the line between the points (C K , µ(C K )) and (C L , µ(C L )), which, under the Hypothesis 1, this line is parallel to the vector n e , with e = K ∩ L ∈ E i h .This expression is considered as an approximation of the discontinuous numerical normal flux ∇µ • n e for µ ∈ P disc 1 (T h ).In addition, observe that, since the baricenters are located 1/3 of the median from the side and 2/3 of the median from the vertex of the triangle, the expression (10) does not degenerate when h → 0, i.e., D e (T h ) > 0 for every e ∈ E i h and h > 0. A visual representation of the regular polygonal structure given by the lines between the adjacent baricenters is given in Figure 2.
Furthermore, in order to preserve the positivity of the variable v in our fully discrete scheme, we assume the following hypothesis.
Hypothesis 2. The mesh T h is acute, i.e., the angles of the triangles of T h are less than or equal to π/2.
We give some examples, the meshes represented in Figure 3, which satisfy both Hypotheses 1 and 2. Theorem 4.2.Given the meshes represented in Figure 3 we can define D e (T h ) for these meshes as follows a) Mesh 1: where l is the length of the side of the highlighted squares of the mesh.
Proof.We will only prove the case a) Mesh 1 since the case b) Mesh 2 is analogous.
Observe Mesh 1 in Figure 3.The baricenter of the triangles △OAB, △OBC and △OCD are, respectively O+A+B

3
, O+B+C 3 and O+C+D

3
. Hence, if we denote by e 1 the edge between △OAB and △OBC and by e 2 the edge between △OBC and △OCD, then Now, since l/|e 1 | = 1 and l/|e 2 | = 1/ √ 2, we can define D e (T h ) for any e ∈ E i h as

Properties of the scheme
Finally, we discuss the different properties of the scheme (8) with the upwind bilinear form defined in (9) on meshes under the Hypotheses 1 and 2.
With the purpose of proving the existence of solution of the scheme (8) and providing a way of computing it enforcing the restriction u m+1 ≥ 0, we define the following auxiliary scheme where a cut-off operator is introduced in (8b): Step for all v ∈ P cont 1 (T h ).
Remark 4.3.In order to preserve the positivity of the solution u m+1 of (13), we have introduced a truncation of this function taking its positive part (u m+1 ) ⊕ in the upwind part of (13b), which is consistent as the solution of the continuous model (1) satisfies u ≥ 0. Since Theorem 4.5 below will guarantee that u m+1 ≥ 0, then log(u m+1 + ε), µ in (13b) is well-defined.
Proof.Proving that if u m ≥ 0 and v m ≥ 0 (when τ > 0) then v m+1 ≥ 0 using that the mesh is acute is a classic result which can be found, for example, in [14,22].
Moreover, if u m ≥ 0 and v m ≥ 0 (when τ > 0) it follows from the equation (13b) that u m+1 ≥ 0 using the same arguments that are shown to prove the positivity result in [1] since the proof given in [1] is independent of the flux ∇µ m+1 .Proposition 4.6.There is a unique solution v m+1 of the linear equation (13a).
Proof.Since we are dealing with a discrete linear problem, existence and unicity of the solution are equivalent.Hence, we just need to assume that there are two solutions of (8a), v 1 and v 2 , substract the expressions resulting of evaluating both solutions and test with v 1 − v 2 to prove unicity of the solution.
Proof.Consider the following well-known theorem: Theorem 4.8 (Leray-Schauder fixed point theorem).Let X be a Banach space and let T : X −→ X be a continuous and compact operator.If the set {x ∈ X : x = α T (x) for some 0 ≤ α ≤ 1} is bounded (uniformly with respect to α), then T has at least one fixed point.
Given u m ∈ P disc 0 (T h ) with u m ≥ 0 and the unique solution v m+1 of (13a), we define the map is the unique solution of the linear (and decoupled) problem: To check that T is well defined, it is straightforward to see that the solutions u of (14a) and µ of (14b) are unique, which involves their existence as P disc 0 (T h ) is a finite-dimensional space.
Secondly, we will check that T is continuous.Let { u j } j∈N , { µ j } j∈N ⊂ P disc 0 (T h ) be sequences such that lim j→∞ u j = u and lim j→∞ µ j = µ.Taking into account that all norms are equivalent in P disc 0 (T h ) since it is a finite-dimensional space, the convergences u j → u and µ j → µ are equivalent to the elementwise convergences ( u j ) K → u K and ( µ j ) K → µ K for every K ∈ T h (this may be seen, for instance, by using the norm ∥•∥ L ∞ (Ω) ).Taking limits when j → ∞ in ( 14) (with u . .= u j , µ . .= µ j and (u, µ) . .= T ( u j , µ j )), using the notion of elementwise convergence and the fact that log(( u) ⊕ + ε) is continuous, we get that hence T is continuous.In addition, T is compact since P disc 0 (T h ) have finite dimension.Finally, let us prove that the set is bounded (independent of α).The case α = 0 is trivial so we will assume that α ∈ (0, 1]. Now, testing (15) with u = 1, we get that and, as u m ≥ 0 and it can be proved that u ≥ 0 using the same arguments than in Theorem 4.5, we get that Moreover, since u ≥ 0, µ ∈ P disc 0 (T h ) is the solution of the equation Hence, Thus, taking into account that u is bounded in P disc 0 (T h ), we conclude that µ is bounded in P disc 0 (T h ).Since P disc 0 (T h ) is a finite-dimensional space where all the norms are equivalent, we have proved that B is bounded.
Since every solution of ( 13) is positive according to Theorem 4.5, the schemes ( 13) and ( 8) are equivalent in the sense that any solution of ( 13) is solution of (8) and vice versa.Therefore, using Propositions 4.6 and 4.7 the following result holds.Corollary 4.9.There is at least one solution of the decoupled non-truncated scheme (8).Moreover, v m+1 is unique.Remark 4.10.Obtaining the nonnegative solutions of (8) can be enforced by solving the scheme (13) including the cut-off operator (13b).In practice, the same solution was found in our numerical experiments using either the auxiliary truncated scheme (13) or the non-truncated scheme (8) without explicitly imposing the nonnegativity restriction u m+1 ≥ 0 (see Remark 5.2).
Remark 4.11.Showing uniqueness of solution of (8) is not straightforward and it might require using inverse inequalities that would probably involve some kind of restriction on the time step and mesh size, and this is beyond the scope of this work.Theorem 4.12.Any solution of the scheme (8) satisfies the following discrete energy law at the time step m + 1: where 8) and consider the equalities Then, adding the resulting expressions for (8a) and (8c) and substracting (8b), we obtain Now, using that Hence, using Proposition 4.4, Thus, taking into account Theorems 19 and 20, we obtain the discrete energy law (17).
Corollary 4.13.Given a solution of the scheme (8), the upwind bilinear form defined in (9) satisfies a upw h (µ m+1 ; u m+1 , µ m+1 ) ≥ 0. In consequence, the scheme (8) is unconditionally energy stable with respect to the approximated energy E ε , that is Proof.Since we know that the discrete energy satisfies (17), it suffices to prove that a upw h (µ m+1 ; u m+1 , µ m+1 ) ≥ 0 to show δ t E ε (u m+1 , v m+1 ) ≤ 0. Now, take u = µ m+1 and use the definition (11) of the upwind bilinear form to get the following: Remark 4.14.Notice that the energy stability of the scheme (8) is obtained thanks to the approximation of the flux −∇µ made in the upwind bilinear form a upw h (•; •, •).In addition, the approximation −∇ 0 ne µ on the edges e ∈ E i h requires the assumption of the Hypothesis 1 for the mesh T h as discussed in the Remark 4.1.
In the test 5.1 we reproduce the first numerical experiment shown in the paper [13] by A. Chertock and A. Kurganov.In this paper, they use a scheme that preserves the positivity of both variables u and v, although they do not show any energy related result.Hence, in our case, we can improve the results shown in the aforementioned paper assuring that our scheme preserves the energy law of the continuous Keller-Segel model.
Then, in the test 5.2 we simulate the qualitative behaviour of the solution in the numerical experiment made by N. Saito in [35].However, since the initial conditions used for the experiment in [35] are not specified we cannot reproduce the exact same test shown in this paper.In the case of our numerical test, the qualitative behaviour of the solution is similar to the one in [35] until the mesh is refined enough so that we capture the blow-up phenomenon.
Finally, in the test 5.3 we reproduce the qualitative behaviour of the results in [4,10,39] for different variations of chemotaxis equations and in [21] for the Keller-Segel equations, where pattern formations with multiple peaks are shown.
Remark 5.1.The scheme (8) preserves the positivity and conserves the mass of u m+1 which implies u m+1 ∈ L 1 (Ω).Therefore, we cannot expect an actual blow-up in the discrete case as it occurs in the continuous model.However, we observe in the numerical tests how the mass accumulates in some elements of T h leading to the formation of peaks.In fact, we are able to capture peaks that reach values up to the order of 10 7 using the approximation shown in (8).These kinds of numerical results are not usual in the literature due to the difficulties when approximating the steep gradients that this process involved.
The numerical results shown in this paper have been obtained using the Python library FEniCS, [3].In order to improve the efficiency of the code, these have been run in parallel using several CPUs.
For the sake of a better visualization of the results, a P cont 1 -projection of u is represented in 3D using Paraview, [2].
Remark 5.2.All the tests were carried out using both the non-truncated equation (8b) without the restriction u m+1 ≥ 0 and the truncated version (13b) to enforce the nonnegativity.The approximations of u and v obtained in every case using both versions of the scheme were identical in all the degrees of freedom.

One bulge of cells
First, we reproduce the results shown in [13].For this purpose, we consider the radially symmetric initial conditions u 0 = 1000e −100(x 2 +y 2 ) , v 0 = 500e −50(x 2 +y 2 ) , which are plotted in Figure 4.As stated in [13], the u and v components of the solution are expected to blow up in a finite time due to the initial conditions chosen.The result of the test using the scheme (8) with h ≈ 1.41 • 10 −3 and ∆t = 10 −6 is shown in Figures 5 and 6.In fact, we observe a blow-up phenomenon for a certain finite time in the range conjectured by A. Chertock and A. Kurganov in [13], t * ∈ (4.4 • 10 −5 , 10 −4 ), as our discrete approximation reaches values of order 10 6 in this time interval.
Moreover, the positivity is preserved for both u and v as stated in Theorem 4.5 and, unlike the scheme presented in [13], we are certain that the discrete energy decreases (both for E(•, •) and E(•, •)) using the scheme (8) as proved in Theorem 4.12.See Figures 7 and 8.
Remark 5.3.This test have been computed with greater and lower values of ε including the limiting case ε = 0.The difference in norms L 2 and L ∞ are shown in Table 1.From a qualitative point of view, the solutions are indistinguishable.
In this case, the numerical approximation works with ε = 0 since the minimum of u remains strictly positive and does not tend to 0, it takes values around 10 −19 during all the iterations computed.Below, in Remark 5.4, we show a different test where we do have to take ε > 0 to ensure convergence of the scheme.An accuracy test in space (with ε = 10 −10 ) has been also carried out where the solution obtained with h ≈ 7.071 • 10 −4 and ∆t = 10 −6 has been taken as reference solution.The results shown in Tables 2 and 3 suggest first order of convergence in space in norm L 2 both for u and v.
It is remarkable to notice that, although the L 2 errors of the approximation of u may seem huge at first, particularly as it approaches the blow-up time, they are not that big in relative terms.As it can be observed in Figure 5, the maximum value reached by u is around 10 3 bigger than the L 2 errors shown in Table 2 at each time step.These errors will tend to vanish as the mesh is refined so that the spiky bulge in the middle of the domain is more accurately approximated.
Also, we would like to emphasize the difficulty of achieving such results as obtaining a reference solution in a blow-up situation where the exact solution tends to degenerate and huge gradients appear require a significant computational effort.In this regard, the reference solution has been computed in parallel using a domain decomposition technique.

Three bulges of cells
Now, we show the results for a similar test to the one that appears in [35].In this case, we take the parabolic-elliptic case (τ = 0) so that the characteristic speed of v is much faster than the characteristic speed of u.Moreover, we consider the initial condition u 0 = 900e −100((x−0.2) 2 +y 2 ) + 800e −100(x 2 +(y−0.2) 2 ) + 1000e −100((x−0.3) 2 +(y−0.3) 2 ) , which is plotted in Figure 9.As stated in [35] and the references therein, since ∥u 0 ∥ L 1 (Ω) > 8π, the solution is expected to blow-up in finite time.In Figures 10 and 11 we can observe the result of the test with h ≈ 2.83•10 −2 and ∆t = 10 −5 .In this case, the qualitative behavior of the solution is similar to the one shown in [35], with the peak of cells moving towards a corner of the domain.However, the qualitative behavior of the solution is different if we take h ≈ 7.07 • 10 −3 and ∆t = 10 −5 .Now, as represented in Figures 12 and 13, a blow-up phenomenon seems to occur in finite time and the peak of cells remains motionless far away from the corners of the domain.
In both cases, the positivity is preserved and the energy decreases in the discrete case.See Figures 14 and 16

Pattern formation with multiple peaks
Finally, we show the results for a test in which we obtain a numerical solution describing a pattern with multiple peaks as it occurs, for instance, in the cases that appear in [4,10,39] for different variations of chemotaxis equations and in [21] for the Keller-Segel equations.For this purpose, we

Figure 1 :
Figure 1: Orientation of the unit normal vector n e .

Figure 3 :
Figure 3: Representation of D e (T h ).

u 0 v 0 Figure 4 :
Figure 4: Initial conditions for blow-up as in [13] (different scales are used for u and v). u

Figure 7 :
Figure 7: Minimum and maximum of u and v over time in the case shown in [13].

Table 2 :
Accuracy test in norm L 2 for u (test in subsection 5.1).