A chemorepulsion model with superlinear production: analysis of the continuous problem and two approximately positive and energy-stable schemes

We consider the following repulsive-productive chemotaxis model: find u ≥ 0, the cell density, and v ≥ 0, the chemical concentration, satisfying 1∂tu−Δu−∇⋅(u∇v)=0inΩ,t>0,∂tv−Δv+v=upinΩ,t>0,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left\{ \begin{array}{l} \partial_t u - {\Delta} u - \nabla\cdot (u\nabla v)=0 \ \ \text{ in}\ {\Omega},\ t>0,\\ \partial_t v - {\Delta} v + v = u^p \ \ { in}\ {\Omega},\ t>0, \end{array} \right. $$\end{document}with p ∈ (1, 2), Ω⊆ℝd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${\Omega }\subseteq \mathbb {R}^{d}$\end{document} a bounded domain (d = 1, 2, 3), endowed with non-flux boundary conditions. By using a regularization technique, we prove the existence of global in time weak solutions of (1) which is regular and unique for d = 1, 2. Moreover, we propose two fully discrete Finite Element (FE) nonlinear schemes, the first one defined in the variables (u,v) under structured meshes, and the second one by using the auxiliary variable σ = ∇v and defined in general meshes. We prove some unconditional properties for both schemes, such as mass-conservation, solvability, energy-stability and approximated positivity. Finally, we compare the behavior of these schemes with respect to the classical FE backward Euler scheme throughout several numerical simulations and give some conclusions.

with 1 2 , a bounded domain ( 1 2 3), endowed with non-flux boundary conditions. By using a regularization technique, we prove the existence of global in time weak solutions of (1) which is regular and unique for 1 2. Moreover, we propose two fully discrete Finite Element (FE) nonlinear schemes, the first one defined in the variables under structured meshes, and the second one by using the auxiliary variable and defined in general meshes. We prove some unconditional properties for both schemes, such as mass-conservation, solvability, energy-stability and approximated positivity. Finally, we compare the behavior of these schemes with respect to the classical FE backward Euler scheme throughout several numerical simulations and give some conclusions.

Introduction
Chemotaxis is the biological process of the movement of living organisms in response to a chemical stimulus, movement that can be addressed towards a higher (chemo-attraction) or lower (chemorepulsion) concentration of a chemical substance. At the same time, the presence of living organisms can produce or consume chemical substance. A repulsive-productive chemotaxis model can be given by the following parabolic PDE's system: where 0 and 0 denote, respectively, the cell density and the concentration of a repulsive chemical signal at position ( 1 2 3, being a bounded domain with boundary ) and at time 0. Moreover, 0 (if 0) is the production term. In this paper, we consider the particular case of superlinear signal production, that is, , with 1 2, and then we focus on the initial-boundary value problem: From the biological point of view, the nonlinear signal production considered in model (3) is justified and explains the saturation effects of chemotactic signal production at large (or short) densities of cells (see [32] and references therein). The theoretical analysis of chemorepulsion models has included the study of some qualitative properties of the solutions, such as existence, uniqueness, regularity and behavior at infinite time, among others [9,14,17,30,31]. In the case of linear ( 1) or quadratic ( 2) production term, problem (3) is well-posed (see [9,17] respectively) in the following sense: there exist global in time weak solutions in 3 domains, which are regular (and unique) for 1 and 2 domains. In [14], the uniqueness and global existence of solution for a chemorepulsion model with linear production and superlinear diffusion in domains (for 3) have been proved. In the context of Lotka-Volterra competition models, the effect of a chemorepulsive signal has been considered by Tello and Wrzosek in [31], proving the existence of global classical solution for the model in domains (for 1). A chemorepulsion model with nonlinear chemotactic sensitivity has been studied in [30], obtaining the existence of bounded classical solution and the convergence at infinite time to a constant steady state in domains (for 3). With respect to the study of chemotaxis models with nonlinear signal production ( ) the literature is scarce (we refer [23,32]). In [32], Winkler studied radially symmetric solutions of a parabolic-elliptic system, proving the existence of global bounded classical solution under some conditions on the power . Considering nonlinear chemotactic sensitivity, chemorepulsion and nonlinear production, in [23] Lai and Xiao analyzed the existence, uniform boundedness and asymptotic behavior of global classical solutions also for a parabolic-elliptic model. However, as far as we know, there are not works studying the parabolic-parabolic problem (3) with production (for 1 2). Therefore, the first aim of this work is to study the existence of global weak solutions of (3) (in the three-dimensional case) and global regularity (in the two and one-dimensional cases).
On the other hand, the second aim is to design numerical methods for model (3) conserving properties of the continuous problem such as: mass-conservation, energystability and positivity. It is important to mention that approaching chemorepulsion problems by using Finite Element (FE) approximations is not an easy task, because negative (discrete) solutions can be computed (see [17,19,20]). In those cases, some spurious oscillations may appear (see, for instance, [19] for a chemorepulsion model with quadratic production).
Energy-stable numerical approximations have also been studied in the chemotaxis context. A conditionally energy-stable FV scheme for a chemo-attraction model with an additional cross-diffusion term was analyzed by Bessemoulin and Jüngel [6]. Energy-stability of time-discrete numerical approximations and fully discrete FE schemes for a chemorepulsion model with quadratic production have been analyzed in [17] and [18,19], respectively; while, in the case of linear production, we refer [20]. However, as far as we know, for the chemorepulsion model with production term given in (3) there are not works studying energy-stable numerical schemes. Likewise, the positivity or approximated positivity properties have been studied on numerical schemes for chemotaxis models. In [8], Chamoun and collaborators proved a discrete maximum principle for a combined FV-FE scheme approaching a chemotaxis-fluid model. The positivity of only time-discrete schemes and approximated positivity of a fully discrete FE scheme associated with a chemorepulsion model with quadratic signal production were proved in [17] and [19], respectively; while, for the case of linear production, we refer to [20]. Positive numerical methods, using FE techniques, associated with a generalized Keller-Segel model were studied in [10]. In [34], a positive FV scheme for a parabolic-elliptic chemotaxis model was analyzed. However, there are not works studying positive (or approximately positive) FE schemes for model (3).
The idea here is to extend the analysis made in [20], although in this case, we need to use two matrix operators (see (51) and (52) below) in order to obtain energystability and approximated positivity. The first one is the operator defined in [2] (and used in [20]); while, the second one, is obtained by constructing regularized functions associated with the test function 1 . For the second operator, it was necessary to prove technical Lemmas (see Lemmas 4.4 and 4.10 below) which are requiered in order to obtain the desired properties for the numerical schemes.
Consequently, the main novelties in this paper are the following: The analysis of the existence of weak solutions of model (3) in the 3D case (which are regular and unique in the 2D and 1D cases) satisfying, in particular, an energy inequality (see (8) below). The introduction of a FE scheme (see scheme UV in Section 4.1 below) for model (3) which is energy-stable with respect to an energy in the primitive variables and approximately positive, under a right-angled constraint on the spatial triangulation (see hypothesis (H) in (46) below). The introduction of another FE scheme (see scheme US in Section 4.2 below) for model (3) which is unconditionally energy-stable with respect to a modified energy and approximately positive, without imposing the restriction (H) on the mesh.
The outline of this paper is as follows: In Section 2, we give the notation and some preliminary results. In Section 3, we prove the existence of weak-strong solutions of model (3) (in the sense of Definition 3.1 below) by using a regularization technique. In Section 4, we propose two fully discrete FE nonlinear approximations of problem (3), where the first one is defined in the variables , and the second one introduces as an auxiliary variable. We prove some unconditional properties such as mass-conservation, energy-stability, approximated positivity and solvability of the schemes. In Section 5, we compare the behavior of the schemes with respect to classical FE backward Euler scheme throughout several numerical simulations, including experimental convergence rates; and in Section 6, the main conclusions are summarized.

Notation and preliminary results
Along this paper, we will consider the usual Lebesgue spaces 1 with norm . In particular, the 2 -norm will be denoted by 0 . From now on, will denote the standard 2 -inner product over . The usual Sobolev spaces , for a multi-index , and 1, with norm denoted by will be also considered. If 0 is not integer, the space is a subspace of (where is the integer part of ) of functions with finite norm (see [26]): 1 .
In the case when 2, we denote 2 , with respective norm . Moreover, the following spaces are set n n 0 on (for 1 1 ) and the following equivalent norms in 1 and 1 , respectively, will be used (see [26] and [1, Corollary 3.5], respectively): Here rot denotes the well-known rotational operator (also called curl) which is a scalar operator for 2D domains and vectorial for 3D ones. In particular, (4) implies that, for all If is a general Banach space, its topological dual space will be denoted by . Moreover, the letters will denote different positive constants which may change from line to line. The following result will be used along this paper: Then, the problem in 0 When large time estimates will be treated, the following result will be used (see [21]):

Analysis of the continuous model
In this section, the existence of weak-strong solutions of problem (3) will be proved in the sense of the following definition.
Observe that any weak-strong solution of (3) is conservative in , because the total mass remains constant in time. In fact, by taking 1 in (6): In addition, integrating (7) in , one has .

Regularized problem
In order to prove the existence of weak-strong solution of problem (3) in the sense of Definition 3.1, we introduce the following regularized problem associated with model (3): Let 0 1 , find , with 0 a.e. in 0 , such that, for all 0, Notice that from (12), system (13) is satisfied a.e. in 0 . From now on in this section, we will denote solution of (14) only by . Observe that if is any solution of (13), then (10) and (11) are satisfied for . (12) and (13).

Theorem 3.2 There exists at least one solution of problem
Proof We will use the Leray-Schauder fixed point theorem. With this aim, we denote where and, in general, we denote max 0 . Then, is a solution of (13) iff is a fixed point of the operator defined in (16). Let us check every hypotheses of Leray-Schauder Theorem: 1.
is well defined. Observe that if , from the 2 and 3 -regularity of problem (14) (see [ Then, from (18) and (19) and using (5), we deduce the following estimates with respect to : is bounded in 0 Then, from (20) we conclude that is bounded in . Moreover, testing (17)  from which, taking into account (20) and using the Gronwall Lemma, we deduce that is bounded in . 3.
is compact. Let be a bounded sequence in . Then solves (16) (with and instead of and respectively). Therefore, analogously as in item 1, we obtain that and are bounded in ; and therefore, from Theorem 2.1 we conclude that is bounded in which is compactly embedded in , and thus is compact. Observe that the compactness embedding comes from the continuous embedding (using embeddings , see [24, Theorem 9.6]): , hence the compactness holds by applying the Aubin-Lions Lemma (see [29]). 4.
is continuous from into . Let be a sequence such that in as .
Therefore, is bounded in , and from item 3 we deduce that is bounded in . Then, there exist and a subsequence of still denoted by such that weakly in and strongly in .
Then, from (21) and (22), a standard procedure allows us to pass to the limit, as goes to , in (16) (with and instead of and respectively), and we deduce that . Therefore, we have proved that any convergent subsequence of converges to strong in , and from uniqueness of , we conclude that the whole sequence in . Thus, is continuous.
Therefore, the hypotheses of the Leray-Schauder fixed point theorem are satisfied and we conclude that the map has a fixed point , that is, , which is a solution of problem (12) and (13).

Theorem 3.3 There exists at least one weak-strong solution of problem (3).
Proof Observe that a variational problem associated with (13) Recall that is the unique solution of problem (14). From (18) we have that satisfies the following energy equality: Then, from (24) and using (19) we deduce the following estimates (independent of ) 2 is bounded in 0 Moreover, taking into account that from (25)  . Therefore, we deduce that is bounded in Notice that from (14) and (25) Then, from (25) Thus, from (30) and (31) and using that is bounded in Moreover, since strongly in 0 , we have that strongly in 1 0 Thus, taking to the limit when 0 in (23), and using (29), (32) and (33), we obtain that satisfies with 0 on . Notice that the limit function is nonnegative. In fact, it follows by testing (36) by and using that 0 0. Finally, we will prove that satisfies the energy inequality (8). Indeed, integrating (24) in time from 0 to 1 , with 1 0 0, and taking into account that Observe that, in 1 , 2 or 3 domains, condition (42) reads as for all 0 , for 1, and 3 2, respectively. Since the weakstrong regularity of given in Definition 3.1 only guarantees the boundedness of , therefore this regularity result can only be applied for 1 domains. On the other hand, in [9] it is proved that in 2 domains, assuming 2 0 2 and reasoning over the equation (1) 1 , then 0 for any . Consequently, since the weak-strong regularity guarantees that 2 0 2 and the -equation in our model and in [9] is the same, then one also has the existence and uniqueness of global in time regular solutions for (1) in 2 domains.

Fully discrete numerical schemes
In this section, two fully discrete numerical schemes associated with model (3) are proposed. Some unconditional properties such as mass-conservation, energystability, approximated positivity and solvability of the schemes are proved.

Scheme UV
In this section, in order to construct an energy-stable fully discrete scheme for model (3), we are going to make a regularization procedure, in which we will adapt the ideas of [2] (see also [16] Then, is obtained by integrating in (43) and imposing the conditions Observe that (at least formally) multiplying (45) 1 by , (45) 2 by , integrating over and adding, the chemotaxis and production terms cancel and we obtain the following energy law In particular, the modified energy 1 2 2 is decreasing in time. Thus, we are going to consider a fully discrete approximation of the regularized problem (45) using a FE discretization in space and the backward Euler discretization in time (considered for simplicity on a uniform partition of 0 with time step 0 ). Let be a polygonal domain. We consider a shape-regular and quasi-uniform family of triangulations of , denoted by 0 , with simplices , and max , so that . Further, let a denote the set of all the vertices of , and in this case we will assume the following hypothesis: (H) The triangulation is structured in the sense that all simplices have a right angle.
(46) We choose the following continuous FE spaces for and : 1 2 generated by 1 with 1.

Remark 4.2
The right-angled constraint (H) and the approximation of by 1continuous FE are necessary to obtain the relations (49) and (50) below, which are essential in order to obtain the energy-stability of the scheme UV (see Theorem 4.7 below).
We denote the Lagrange interpolation operator by , and we introduce the discrete semi-inner product on (which is an inner product in ) and its induced discrete seminorm (norm in ): (47)

Remark 4.3
In , the norms and 0 are equivalent uniformly with respect to [4].
We consider also the 2 -projection on , 2 and the classical 1 -projection 1 given by (48) . Moreover, owing to the right-angled constraint (H) and the choice of 1 -continuous FE for , following the ideas of [2] (see also [16]), for each 0 1 , we can construct two operators ( 1 2) such that are symmetric matrices and 1 is positive definite, for all and a.e. in , and satisfy 1 in (49) 2 1 in .
We emphasize that thanks to the choice of 0 made up of simplices (triangles in 2D and tetrahedra in 3D), and the fact that the gradient of a 1 -function is constant over each element of the triangular mesh, the operators ( 1 2) are constant by elements matrices such that (49) and (50) hold in each element . This condition is not satisfied when rectangular meshes are considered or approximation for 2. In the 1-dimensional case, are constructed as follows: For all and with vertices a 0 and a 1 , we set for some 1 2 . Following [2] (see also [16]), these constructions can be extended to dimensions 2 and 3, and from (51) the following estimate holds: The following result will be useful to prove the well-posedness of the scheme UV and we write its proof in the Appendix. where, in general, we denote 1 .

Remark 4.5 (Positivity of )
By using the mass-lumping technique in all terms of (56) 2 excepting the self-diffusion term , and approximating by 1continuous FE, we can prove that if 1 0 then 0. In fact, it follows testing (56) 2 by , where min 0 (see Remark 3.12 in [20]).

Mass-conservation, energy-stability and solvability
Since 1 and 1 , we deduce that the scheme UV is conservative in , that is,   , but is independent of and . Therefore, from the discrete energy law (60) and estimate (66), we have  The hypotheses of the Leray-Schauder fixed point theorem are satisfied as in Theorem 3.13 of [20], but applying in this case Lemma 4.4 in order to prove the continuity of the operator . Thus, we conclude that the map has a fixed point , that is , which is a solution of the scheme UV .

Approximated positivity of u n
In this subsection we are going to prove the property of approximated positivity for solution of the scheme UV , in the sense that 0 as 0 in the 2 norm, where min 0 0. With this aim, we will prove first a preliminary result.   Moreover, from (69)  Remark 4.13 From (71) 1 one has that, in order to guarantee the approximated positivity property for the scheme UV , it is necessary to choose such that 2 1 1 2 0 as 0.

Scheme US
In this section, we are going to construct another energy-stable fully discrete scheme for (3) considering the auxiliary variable an the regularized function 1 . We will also use the regularized functions , and This kind of formulation considering as auxiliary variable has been used in the construction of numerical schemes for other chemotaxis models (see for instance [18,20,33]). Once problem (74)  Observe that (at least formally) multiplying (74) 1 by , (74) 2 by , integrating over and adding both equations, the terms cancel, and we obtain the following energy law In particular, the modified energy 1 2 2 is decreasing in time. Then, we consider a fully discrete approximation of the regularized problem (74) using a FE discretization in space and the backward Euler discretization in time (considered for simplicity on a uniform partition of 0 with time step 0 ). Concerning the space discretization, we consider the triangulation as in the scheme UV , but in this case without imposing the constraint (H) related with the right-angled simplices. We choose the following continuous FE spaces for , , and : We recall that is the Lagrange interpolation operator, and the discrete semi-inner product was defined in (47). Once the scheme US is solved, given 1 , we can recover solving: Given and 1 , Lax-Milgram theorem implies that there exists a unique solution of (76). Moreover, notice that the result concerning to the positivity of solution of scheme UV established in Remark 4.5 remains true for in the scheme US .

Mass-conservation, energy-stability, solvability and approximated positivity
Observe that the scheme US is also conservative in (satisfying (57)), and we have the following behavior for : 1 .   , but is independent of and . Therefore, from the discrete energy law (78) and estimate (80), we have   Proof The proof follows as in Theorem 4.6 of [20], by using the Leray-Schauder fixed point theorem.

Numerical simulations
In this section, we will compare the results of several numerical simulations using the schemes derived through the paper. The spaces for , and have been generated by problem (3), which is given for the following first order in time, nonlinear and coupled scheme: Scheme UV: Initialization: Let 0 0 an approximation of 0 0 as 0. Time step n: Given 1 1 , compute by solving .

Remark 5.1
The scheme UV has not been analyzed in the previous sections because it is not clear how to prove neither its energy-stability nor its approximated positivity. In fact, observe that the scheme UV (which is the "closest" approximation to the scheme UV considered in this paper) differs from the scheme UV in the use of the regularized function and its derivatives (see Fig. 1) and in the approximation of the cross-diffusion and production terms, and respectively, which are crucial for the proof of the energy-stability of the scheme UV , and consequently for the approximated positivity. Note that a residual term 1 is considered. This term is required in order to improve the convergence of this iterative method. Indeed, since the self-diffusion term of the u-equation is rewritten in a nonlinear form, we have checked that this fact makes the convergence of the corresponding iterative method worse. (iii) Picard method to approach a solution of the scheme UV: Given , compute 1 1 solving the decoupled problems

Positivity of u n
In this subsection, the positivity of the variable in the three schemes is compared. We recall that for the two schemes studied in this paper, namely schemes UV and US , the positivity of the variable is not clear. However, it was proved that 0 as 0 (see Note that 0 0 0 in , min 0 0 1 1 0.01 and max 0 0 1 1 80.01. We obtain that: (i) All the schemes take negative values for the minimum of in different times 0, for the different values considered for and . However, in the case of the schemes UV and US , it is observed that these values are closer to 0 as 0 (see Figs. 3, 4, and 5). (ii) In the cases 1.1 y 1.5, the scheme US "preserves" better the positivity than the other schemes; while for 1.9, the scheme UV evidence "better" the positivity (see Figs. 3, 4, and 5).

Energy-stability
In this subsection, we compare numerically the stability of the schemes UV , US and UV with respect to the "exact" energy It was proved that the schemes UV and US are unconditionally energy-stables with respect to modified energies defined in terms of the variables of each scheme, and some energy inequalities are satisfied (see Theorems 4.7 and 4.15). However, it is not clear how to prove the energy-stability of these schemes with respect to the "exact" energy given in (81), which comes from the continuous problem (3) (see (8) and (9)). Therefore, it is interesting to compare numerically the schemes with respect to this energy , and to study the behavior of the following "residual" of the discrete energy law Then, we obtain that: (i) All the schemes UV , US and UV satisfy the energy decreasing in time property for the exact energy (see Fig. 7a), that is, (ii) All the schemes show 0 for some 0; being those corresponding to scheme UV that reach higher values, while the scheme US evidence the smallest values. Moreover, it is observed that the scheme UV introduces lower numerical source than the scheme UV, and lower numerical dissipation than the scheme US (see Fig. 7b).

Experimental convergence rates
In order to show the accuracy of the schemes proposed in this paper, we compare the schemes UV , US and UV against an exact solution and on several meshes. With this aim, in this experiment we consider the exact solution Note that 0 on . Moreover, we use a uniform partition with 1 nodes in each direction. We consider 0 1 2 and 10 6 . Numerical results of convergence rates in space are listed in Tables 1, 2, 3, 4, 5, 6, 7, 8 and 9 for 5 10 5 with respect to the final time 0.1. We denote the total errors by and . For the three schemes UV , US and UV, and different values of , we obtain optimal order of convergence in space, that is, second-order for in 2 -norm and first order in 2 1 -norm.

Conclusions
In this paper, the existence of global in time weak solutions for the chemorepulsion with -power production model (3) and satisfying the energy inequality (8) has been proved in the 3D case, which are regular and unique in the 2D and 1D cases. In addition, two new mass-conservative, unconditionally energy-stable and approximated positive fully discrete FE schemes for model (3), namely UV and US have been developed. From the theoretical point of view, the following statements have been deduced: (i) The solvability of both schemes. (ii) The scheme UV is energy-stable with respect to the modified energy (given in (59)), under the right-angled constraint (H); while the scheme US is unconditionally energy-stable with respect to the modified energy (given in (77)), without this restriction (H) on the mesh. (iii) It is not clear how to prove the energy-stability of the scheme UV (see Remark 5.1).
(iv) In the schemes UV and US there is a control for in 2 -norm, which tends to 0 as 0. This allows to conclude the non negativity of the solution in the limit as 0.
On the other hand, from the numerical simulations, the following deductions can be made: (i) The three schemes have decreasing in time energy , independently of . (ii) All the schemes show 0 for some 0; reaching highest values the scheme UV, and the smallest values the scheme US . Moreover, scheme UV introduces lower numerical source than the scheme UV, and lower numerical dissipation than scheme US . (iii) Both schemes UV and US satisfying that min 0 0 as 0.
(iv) All the schemes have optimal order of convergence in space, independent of the -values.