{'fraction_non_alphanumeric': 0.04772894157138362, 'fraction_numerical': 0.01542011958460086, 'mean_word_length': 4.110426159206647, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scalar auxiliary variable finite element scheme for the parabolic-parabolic Keller-Segel model

We describe and analyze a finite element numerical scheme for the parabolic-parabolic Keller-Segel model. The scalar auxiliary variable method is used to retrieve the monotonic decay of the energy associated with the system at the discrete level. This method relies on the interpretation of the Keller-Segel model as a gradient flow. The resulting numerical scheme is efficient and easy to implement. We show the existence of a unique non-negative solution and that a modified discrete energy is obtained due to the use of the SAV method. We also prove the convergence of the discrete solutions to the ones of the weak form of the continuous Keller-Segel model.


Introduction
Since chemotaxis is observed very widely in various areas of biology and medicine, it becomes a prolific subject in mathematical biology throughout the past decades. Among the different mathematical models used to represent chemotaxis of living organisms, the Keller-Segel equation is one of the most recognized. It has been introduced by Keller and Segel [14] to depict the movement of the Dictyostelium discoideum toward the location of high concentration of adenosine 3', 5'-cyclic monophosphate. The parabolic-parabolic Keller-Segel model (KS in short) is often set in a bounded domain Ω ⊂ R d , d = 1, 2, 3 with a Lipschitz boundary ∂Ω and reads ∂ t u = ∇ · (D u ∇u − χ c ϕ(u)∇c) in Ω × (0, +∞), where ν is the outward normal vector to the boundary. We assume in the following that the initial condition satisfies {u(0, x), c(0, x)} = {u 0 , c 0 } ∈ H 1 (Ω) × H 1 (Ω), and 0 ≤ u 0 ≤ 1 a.e. in Ω, 0 ≤ c 0 ≤ C a.e. in Ω, (1.4) where C is a positive finite constant. In the model (1.1)-(1.4), the cell density u(t, x) is attracted by the chemo-attractant given by c(t, x), its concentration. Cells can move randomly by diffusion with a coefficient of diffusion given by D u and by chemotaxis with χ c , a coefficient used to represent the strength of this movement. A small parameter τ > 0 is used to denote how fast the chemo-attractant is diffusing compared to the cells. Without a loss a generality, we will assume in the following that τ = 1.
ϕ(u) is the chemosensitivity and is given by This particular form of chemosensitivity prevents the unrealistic scenario of overcrowding of cells and therefore the blow-up of the solution. Due to this possible behavior of solutions, the Keller-segel system exhibits very interesting mathematical structure and the interested reader can refer to the review [13] and the work of Blanchet et al. [4]. The volume filling strategy was proposed in [15]  where B = D u /χ c and we define the integral of the free energy density Here, C is a positive constant such that F (u) > 0, ∀u ∈ [0, 1]. For latter convenience, we denote F ′ (u) = g(u), and we remark that g ′ (u) = 1 ϕ(u) . Thus, we can express the Keller-Segel model using its gradient flow structure [5] where the variational derivatives of the energy functional with respect to u and c are given respectively by Generally, a numerical scheme for gradient flow model is evaluated by several aspects: i) its capacity to keep the energy dissipation; ii) if it is convergent, and if error bounds can be established; iii) its efficiency; iv) its implementation simplicity. For a large a class of gradient flows, the Scalar Auxiliary Variable (SAV in short) [24] has shown to meet all the previous points. Applying this method to the Keller-Segel model is only possible starting from its gradient flow formulation (1.8)-(1.9) and gives what we call the SAV Keller-Segel model (1.10) and we define the scalar unknown In this article, we propose to study a finite element scheme to simulate the system (1.10)-(1.15) that is known to preserve the energy at the discrete level. Throughout the past decades, the Keller-Segel model as been at the center of many pieces of research. The analytical properties of the Keller-Segel model without volume filling have been extensively studied. One of the most important result was to show that the solution of the model blows up in finite time if a certain constraint on the initial mass is not satisfied. For the reader interested into the analytical results about this model without volume filling, we refer to the review paper [25]. The volume filling approach prevents this blow-up of the solution in finite time for any initial condition satisfying (1.4). Moreover, it seems to be more biologically relevant since it takes into account the finite size of the cells. A more general form of the Keller-Segel model is where the random movement of the cells is given by D u β(u) (that can be non-linear) and the chemosensitivity is given by χ c µ(u). Particular assumptions on both β(u) and µ(u) can be made to prevent the blow-up of solutions in finite time. The introduction of the parabolic-parabolic KS with volume filling and quorum-sensing is presented in the work of Painter and Hillen [15]. They described a discrete lattice model where the probability for cells to jump to a different location is dependent on the local density and on the concentration of the chemotactic agent. From this discrete model, they derived the continuous limit model and give the following conditions for β(u) and µ(u) where ψ(u) is a monotonically decreasing function. Under the assumptions ψ(0) > 0, ψ(u) > 0 for 0 < u < u, and ψ(u) = 0, the same authors proved the global existence of classical solutions in [12]. They also presented some numerical simulations where they were able to make observation of the behavior of the solution for longer times since the blow-up of the solution is prevented by the model. Many other variations of the Keller-Segel model have been proposed to take into account the effect of volume filling. For example, more recently, Bubba et al. [8] proposed to take where u max represents the density at which cells are too overcrowded. Numerical methods for the Keller-Segel model are numerous. Considering zero-flux boundary conditions, the conservation of the total mass of the cells, the non-negativity of the solution and the capacity to retrieve the energy at the discrete level are the key properties expected from a numerical scheme for this equation. For the parabolic-elliptic Keller-Segel equation where the equation for the chemo-attractant is given by −∆c = δu − αc, Saito and Suzuki proposed a conservative finite-difference scheme [20]. For the parabolic-parabolic version, Saito proposed and performed an error analysis for an upwind finite element scheme [17,18,19] using Baba and Tabata's method [3]. The finite volume method has also been applied for this problem: we can cite the work of Filbet [10] that deals with the classical Patlak-Keller-Segel model (without volume filling) and the work of Almeida et al. [2]. In the latter, the parabolic-elliptic model is used and the authors were able to prove the preservation of the important properties for two finite volume schemes. The difference between the two is that one uses the gradient flow structure of the model while the user uses an exponential rewriting inspired by the Scharfetter-Gummel discretization. The scheme we propose in this article follows the same idea. The Keller-Segel model has a gradient flow structure that can be useful for its numerical simulation. A recent numerical method to simulate gradient flows that ensures that the energy is preserved at the discrete level is the SAV method [23,24]. This method provides a robust framework to simulate gradient flows in an efficient way. In fact, the computation of the solution of any gradient flow model requires only the solving of two decoupled linear systems at each time step. This method has shown very interesting results for the simulation of the Cahn-Hilliard equation [21] for which the properties concerning the discrete energy and the conservation of the total mass are of main importance. We must stress that the energy recovered by the SAV method is a modified version of the energy of the real system. This is due to the discretization of the equation for the scalar variable. In a recent work of Bouchriti et al. [6], the authors showed that the use of the SAV method for the damped wave equation and the Cahn-Hilliard equation leads to the convergence to modified steady states as well.
In a recent work Shen and Xu [22] proposed an unconditionally energy stable method that is able to preserve the nonnegativity of the solution. This method relies on the use of the gradient flow structure of the Keller-Segel model and requires the solving of a nonlinear convex system.
To the best of our knowledge the SAV method has never been applied to the Keller-Segel model. The principal difference with previous works on the SAV method is that the mobility in the first equation of the Keller-Segel system is not constant through time, leading to the necessity to compute at each time step the associated matrix.
Therefore, in this article, we propose to use it to obtain a new model that we discretize in space using the finite element method. Altogether, we obtain a new way to simulate the parabolic-parabolic Keller-Segel equation with the certitude to be able to retrieve the energy associated with the model at the discrete level. First, we describe the method and explain the strategy to solve the resulting equations. Then, the well-posedness of the scheme is studied. We show the existence of a unique pair of solution that is non-negative and retrieve the expected L ∞ norm under some constraints on the spatio-temporal mesh. We also show that the initial mass of the cells is conserved. We prove that a modified energy is retrieved at the discrete level which is an inherent property of the SAV method. Lastly, a convergence analysis is performed. The convergence of subsequences in Bochner spaces can be proved and the system solved by the limit solution is the weak form of the Keller-Segel equation.

Finite element framework
Let L p (Ω), W m,p (Ω) with H m (Ω) = W m,2 (Ω), where 1 ≤ p ≤ +∞ and m ∈ N, be respectively the usual Lebesgue and Sobolev spaces. The corresponding norms are respectively || · || m,p,Ω , || · || m,Ω and semi-norms | · | m,p,Ω , | · | m,Ω . We denote L p (0, T ; V ) the Bochner spaces i.e. the spaces with values in Sobolev spaces [1]. The norm in these spaces is defined for all function η Bochner measurable by The standard L 2 inner product is denoted by (·, ·) Ω and the duality pairing between (H 1 (Ω)) ′ and Let Ω be a polyhedral domain and Since the mesh is assumed to be quasi-uniform, we know that it is shape-regular and it exists a positive constant C such that Since the domain is assumed to be polyhedral, the discrete domain Ω h exactly coïncides with the domain Ω. Hence, the closure of the domain can be written as the union of all the mesh elements We assume that the mesh is acute, i.e. for d = 2 the angles of the triangles can not exceed π 2 and for d = 3 the angle between two faces of the same tetrahedron can not exceed π 2 . We define by κ T the minimal perpendicular length of T and κ h = min T ∈T h κ T . We introduce the P-1 finite element space associated with the mesh T h where P 1 (T ) denotes the space of polynomials of order 1 on T . For latter convenience, we indicate the set of nodes of T h by J h and {x j } j=1,...,|J h | is the set of their coordinates. N h = |J h | stands for the total number of nodes. We denote by Λ i the set of nodes connected to the node x i by an edge ..,N h is the standard Lagrangian basis functions associated with the spatial mesh.
The standard interpolation operator is defined by For latter convenience, we state here some well-known results for the P-1 finite element method (see for e.g. [7], [16]) We define the standard mass M and stiffness K finite element matrices For the efficiency of the numerical scheme, it could be useful to use the lumped mass matrix which is a diagonal matrix with each term being the sum of the terms on the same row of the standard mass matrix.
From the hypothesis we made on the acuteness of the triangulation, we know that (see [11]) Therefore, we know that the non-diagonal entries of the stiffness matrix K and of the matrix A defined below by the equation (2.10) are non-positive.

Fully discrete scheme
Given N T ∈ N * , let ∆t := T /N T be the constant time-step and t n := n∆t, for n = 0, . . . , N T − 1. We consider a partitioning of the time interval [0, T ] = N T −1 n=0 [t n , t n+1 ]. We approximate the continuous time derivative using a forward Euler The finite element numerical problem associated with the system (1.10)-(1.15) is: are respectively the finite element approximations of the cell density u and the concentration of the chemo-attractant c. We also have used the . We add to this system the following initial conditions (2.9)

Matrix formulation
Let us define A the finite element matrix associated with the right-hand side of (2.4) and the variable . (2.11) We denote in the following by capital letters the vectors associated with the quantities denoted by small letters in the finite element problem. Therefore, the system (2.4)-(2.8) can be rewritten into a matrix formulation

Linear system
Replacing (2.14) into (2.12) but also (2.15) into (2.13), we obtain the system where we have used the notation Multiplying equation (

Existence of a discrete non-negative solution
Theorem 1 (Existence of a unique non-negative discrete solution) Let d ≤ 3 and assume that κ h > 0, ∆t > 0 such that Given an initial condition {u 0 h , c 0 h } such that (1.4) and (2.9) are satisfied, there are two positive con-

2)
and where c is a positive and finite constant.

Proof.
Step and we know that f is a continuous function of both of its arguments. Therefore, from the equation h , c n+1 h } ∈ V h ×V h of the problem.
Step 2: Conservation of mass. To prove mass conservation, we use the identity Summing over the nodes, we get Using the symmetry of the matrix A, the property (3.6) and the fact that the mesh is acute, we obtain which implies mass conservation .
Step 3: Non-negativity and L ∞ bound for {u n+1 h , c n+1 h }. Using the equation (3.4), we find for all where the diffusion coefficient is given by

8)
Therefore, for each node x i ∈ J h , we have Because the spatio-temporal mesh satisfies the conditions (3.1) and one can easily find that 0 ≤ u n+1 h ≤ 1. The only difficulty with this latter is because it depends on r n+1 . However, as seen in equation (3.5), this term can be calculated from the solution {u n h , c n h }. Furthermore, this coefficient is bounded at all time t n since we assumed that r n is bounded and u n h ∈ [0, 1]. For the condition (3.10), since c n h is bounded and in V h , b n ∞ remains bounded from above at all time. Therefore, there are two positive constants C 1 and C 2 such that for n = 1, ..., N T , and we obtain the two conditions (3.2), (3.3). Under these time-dependent assumptions, we proved that the numerical scheme preserves the bounds From this result, the non-negativity and the existence of an upper bound c such that

Discrete energy a priori estimate
Since we are using the SAV method, we are preserving the energy at the discrete level. and The same can be done starting from equation (2.5) to obtain Therefore, summing the two previous equations, we obtain from which we conclude (3.13). Consequently, we already recover the monotonic decay of the discrete energy. To obtain the expression of the energy, we replace φ = u n+1 However, using the equation (2.8), we have Moreover, using the inequality a(a − b) ≥ 1 2 a 2 − b 2 , we get (3.14) Then, performing the same calculations starting from the equation (2.7), we obtain Summing equation (3.14) with (3.15), we obtain the inequality from which we deduce the definition and the decay of the discrete energy (3.12)-(3.13).
Remark 3 From the fact that both u n+1 h and c n+1 h are bounded (see theorem 1), the energy defined by (3.12) is bounded from below and can be used to obtain inequalities.

Notations
To prove the convergence of the discrete solutions, we need some further notations. We define the sequence of approximate solutions and each of them lies in the Cartesian product space To construct the sequences u h∆t , c h∆t , we define the linear operator S ih∆t : These two operators are inductively defined by the system

Preliminary results
Let us define the quantity Proposition 4 (Inverse inequality) Assuming that the mesh is quasi-uniform, the quantity (4.3) is finite and we have where C is a positive constant.
Proof. The proof is given in [9], corollary 1.141 on global inverse inequalities. Proof.

Stability of the scheme
Proof of the inequality (4.4). Starting from equation (4.1), taking φ = 2∆tu n+1 h and using the property However, using the coercivity of the operator a(t n+1 , ·, ·) = (∇·, ∇·), there is a positive constant α such that

From the inverse inequality (4.3) and Young's inequality, there is
Similarly, there is a constant 0 < κ 2 < 1 such that Altogether, we obtain the inequality (4.6) With the choice ϕ(s) = s(1 − s), we assume that κ 1 ,κ 2 and ∆t satisfy where the second condition is strongly related to (3.3). Hence, using the previous assumptions together with c p u n+1 h 0 ≤ u n+1 h 1 , it exists a positive constant C such that The right-hand side of the previous inequality is bounded by a constant that we denoted C 1 due to the energy inequality. Then, since we assumed that the conditions (4.7) holds, it exists a positive constant C 2 such that summing the equation (4.6) from n = 0 → N T − 1 gives where the right-hand side is bounded using the energy estimate (3.12)-(3.13) and we obtain (4.4).
Proof of the inequality (4.5). Starting from the equation (4.2) and taking φ = 2∆tc n+1 h , we have Furthermore, we know that it exists a positive real value α 2 such that . Therefore, again using Young's inequality, for 0 < κ 3 < 1, we have We assume that the condition 0 ≤ 2α 2 c p − κ 3 , is satisfied. Hence, summing from n = 0 → N T − 1 and using (4.4), we obtain Moreover, from the equation (4.8), we know that Summing the previous inequality from n = 0 → N T − 1, we have which gives (4.5).

Conclusion
We presented the application of scalar auxiliary variable method to the parabolic-parabolic Keller-Segel with volume filling using the gradient flow structure of the model. The resulting equations were approximated using a simple P − 1 finite element method. The system is composed of two linear decoupled equations that can be solved efficiently. We were able to prove for this system the existence of a unique non-negative solution and the preservation of the monotonic decay of the discrete energy. We must stress that from the use of the SAV method, the energy that we are able to recover is a modified version of the standard one. However, we were to prove that, in the limit of the discretization parameters, subsequences of solutions converge in Bochner spaces and the limit is the solution of the weak form of the Keller-Segel model.