ON EXPLICIT LOWER BOUNDS AND BLOW-UP TIMES IN A MODEL OF CHEMOTAXIS

. This paper is concerned with a parabolic Keller-Segel system in R n , with n = 2 or 3, under Neumann boundary conditions. First, important theoretical and general results dealing with lower bounds for blow-up time estimates are summarized and analyzed. Next, a resolution method is proposed and used to both compute the real blow-up times of such unbounded solutions and analyze and discuss some of their properties.


1.
Introduction. In 1970 Keller and Segel proposed a mathematical model of chemotaxis phenomena; it consists of a system of parabolic equations (see [9]), whose solution describes the movement of cells in a bounded domain in response to the presence of a chemical substance which is inhomogeneously distributed in space. The mathematical formulation of this problem is prescribed by the following system of PDEs: (1) In system (1) the spatial distribution of the cells/bacteria is identified with the function u and that of the chemical substance/chemoattractant with v. The chemoattractant spreads diffusively and decays with rate k 2 ; it is also produced by the bacteria with rate 1. The bacteria diffuse with mobility µ and also drift in the direction of the gradient of the concentration of the chemoattractant with velocity χ|∇v|; χ is called chemosensitivity. Furthermore, the parameter τ ≥ 0 is a measure of the ratio of the time scales of the cell movement and the distribution of the chemical. Once an initial distribution for both the bacteria and the chemical is given, the homogeneous Neumann boundary conditions represent the most natural and real situation of system (1); no net flux of the distributions with the external domain is permitted. By rescaling the equations, we can always rewrite (1) with τ = 0 or τ = 1; in the former case the model is called parabolic-elliptic and in the latter parabolic-parabolic. In particular, for the parabolic-elliptic the spread rates degenerate to infinite, i.e. the chemical diffuses much faster than the cells move. We refer to paper [8] for a series of results dealing with some properties of the solutions of (1) in the parabolic-elliptic case.
In this manuscript we are interested in chemotactic collapse, a phenomenon connected to the same model (1); it is experimentally observed that the bacteria may concentrate in one or more points (see [7]). On the other hand, as discussed in [3] and [4], the so called isothermal collapse (related to relativistic or ultra-cold gases) represents another possible singular scenario.
From a mathematical point of view, the collapse corresponds to the blow-up of the solution at some finite time t . In this sense, the blow-up phenomena of solutions to various problems, particularly for nonlinear parabolic systems, have received considerable attention (see, for instance, [11]- [16]).
Specifically, this work is concerned with the behavior of nonnegative and classical solutions (u, v) = (u(x, t), v(x, t)), with x = (x 1 , . . . , x n ) for n = 2 or 3, of the following system where Ω is a regular domain of R n , whose boundary ∂Ω is sufficiently smooth, ν = (ν 1 , . . . , ν n ) stands for the unit outer normal to ∂Ω, k I (I = 1, 2, 3, 4) are positive constants, and u 0 (x) and v 0 (x) are nonnegative functions satisfying the corresponding compatibility conditions on ∂Ω.
We focus our attention on lower bounds and blow-up times of unbounded solutions of problem (2). Precisely, in §2 we recall the main results obtained by Payne and Song in [17]; in this work the authors derive a first order time dependent differential inequality by means of an appropriate auxiliary function (herein called W -norm) of the solution of (2), depending on u and ∆v. Subsequently, they give an explicit lower bound of t in the sense of this norm, for both n = 2 and n = 3, by means of analytical techniques. This lower bound estimate depends on the geometry of Ω, on the data k I (I = 1, 2, 3, 4), and on u 0 (x) and v 0 (x).
Furthermore, there exist numerous papers devoted to the quantitative analysis of blowing up solutions of problems defined on bounded or unbounded domains (see for example [1] and [2]). In this sense, §3 deals with the resolution of problem (2); starting from its weak formulation, we propose an algorithm based on a mixed Finite Element Method in space and Euler Method in time (see [10]) capable of numerically solving the above system. This resolution approach is implemented in the 2D case to analyze the behaviors of the aforementioned W -norm as well of the infinity norm (indicated with L-norm) of some blowing up solutions of (2); more precisely, the values of the corresponding blow-up times in these two norms and evolutions of the blow-up spatial point are studied.
Finally, the paper is complemented by some conclusions ( §4).
2. An explicit lower bound: main results. To derive a lower bound of the blow-up time t of a solution (u, v) of (2), let us introduce the auxiliary function for some positive constant α to be specified, and also give the following As we said, the details of these two theorems are available in paper [17]; precisely, by using the same nomenclature we have: (Lower bound for t in the 2D case; n = 2). Let (u, v) be a classical solution of (2), Ω being a bounded and convex domain in R 2 , with the origin inside. Assume α = k 2 4 k2 . If (u, v) becomes unbounded in W -norm (3) at some finite time t , then an estimate from below for t is given by where  k2 . If (u, v) becomes unbounded in W -norm (3) at some finite time t , then an estimate from below for t is given by where W 0 = W (0) = α Ω u 2 0 dx + Ω (∆v 0 ) 2 dx, A and B being positive constants depending on the geometry of Ω and on the constant k I (I = 1, 2, 3, 4).

Remark 1.
By analyzing the proofs of Theorems 2.2 and 2.3, it is seen that since the constants A 1 , B 1 , A and B are not uniquely determined, neither estimate in (4) and in (5) is. In this sense they represent one of the possible lower bounds of t in the W -norm. On the other hand, the integrals in (4) and in (5) can be explicitly computed: and First of all, let us underline that the expression A (7) is well defined. Moreover, (6) and (7) decrease with W 0 increasing (actually also t does, as shown in examples of §3.2.1); as a consequence, if W 0 is large enough, both estimates are very close to zero. Therefore, since the existence of blowing up solutions of (2) is also linked to the magnitudes of its initial data (i.e., W 0 ), we expect that the real value of t is in general considerably larger than the estimates given in the two theorems (see some examples in [19]).
We conclude this section by mentioning that in [15] the authors of this manuscript investigate system (2) under a more general assumption, consisting of choosing k I = k I (t) (I = 1, 2, 3, 4) and non constant; in fact, the time dependence of these coefficients may affect the behavior of the solution (u, v).
3. Numerical resolution method. Let us give a resolution method for system (2). We propose a mixed Finite Element Method in space and Euler Method in time algorithm. Let us point out that even though Euler's approximation produces results with first order precision, it is totally appropriate to the aims of our research (see [19]). (10) and (11) (um, vm) computing norm at step m by (13) ) ) (10) and (11) O O Solution: t ≈ m∆t Table 1. Computation of the blow-up time t . The necessary input data are the threshold ε 0 , the time step ∆t and the initial datum (u 0 , v 0 ); subsequently, it is possible to calculate the sequences (u m , v m ) and ||(u, v)|| m and, consequently, to compute t . Let us clarify that if ||(u, v)|| m = W m (respectively, L m ) t represents t W (respectively, t L ) 3.1. Finite element method: semi-discretization in space. If a mesh of Ω ⊂ R n (n = 2 and 3) is fixed and N represents the total number of nodes of Ω, let (U, V) be the numerical approximation of the solution (u, v) of (2): therefore, by separating variables where ϕ i (x) is the standard quadratic basis function at the vertex x i , for i = 1, ..., N . Thanks to the divergence theorem and the homogeneous boundary conditions of system (2), by multiplying its first two equations by a generic test function ϕ j (x), the following variational form in space is achieved: where (·, ·) denotes the usual L 2 inner product. To compute the time evolutions of both coefficients u i and v i appearing in (8), let ∆t = t m+1 −t m be a given time step, with m = 0, 1, 2, . . . (t 0 = 0), and (U m , V m ) the approximation of (U(x, t), V(x, t)) at time t m . By applying an implicit Euler finite difference approximation to system (9), it is seen that i.e., taking into account (8), and with M ∈ R N ×N (mass matrix), K ∈ R N ×N (stiffness matrix) and In this way the continuous solution of the nonlinear system (2) is identified to the discrete solution of the linear system composed by (10) and (11) 3.2. Numerical algorithm: computing the blow-up times. Theorems 2.2 and 2.3 return an estimate of a lower bound for the blow-up time of positive solutions of (2) in terms of W -norm (3). In this section we will analyze the blow-up phenomena also by considering the natural infinity norm (indicated with L-norm) of a solution (u, v) of the same system (see [7]): By means of the numerical approach presented in §3.1, we want to compare the blow-up times in the sense of the W -norm and the L-norm. We indicate these values with t W and t L ; of course, they depend on the temporal evolutions of the W and L norms, respectively.
To this end, first of all, we need to numerically quantify the values of the W -norm and the L-norm at step m, i.e. W m = W (t m ) and L m = L(t m ); according to (3) and (12), we can define with α as in Theorems 2.2 and 2.3. Therefore, we propose the following algorithm to calculate t W (i.e. t L ) for unbounded solutions of (2). Let ε 0 be a fixed threshold: once the initial data (u 0 , v 0 ) is given, u 1 and v 1 are computed from (10) and (11), respectively. Subsequently, (u 1 , v 1 ) is used to update (u 2 , v 2 ), and so on; we exit the loop when the W -norm (i.e. the L-norm) of the solution at step m, W m (i.e. L m ), is larger than the initial threshold ε 0 (Stopping Criterion); consequently, t W ≈ m∆t (i.e. t L ≈ m∆t). The scheme in Table 1 summarizes this algorithm.
3.2.1. Numerical tests for n = 2. Let us consider the domain Q = Ω × R + 0 , being Ω = [−2, 2] × [−2, 2] the square with center in the origin of the axes O = (0, 0) and length 4. We take a uniform mesh, obtained by dividing each side of the square into 200 equal parts (i.e. 40401 vertexes and 80000 triangles). We also choose v 0 (x) = 0.55e −(x 2 +y 2 ) (4 − x 2 ) 2 (4 − y 2 ) 2 as initial condition (the chemical signal at time t 0 = 0), and ε 0 = 10 6 as the threshold and ∆t = 10 −4 as the integration step. Moreover, we fix the following values for k I : k 1 = 0.2, k 2 = k 4 = 1 and k 3 = 0.1. Once these data are set, we solve (2) in two cases, depending on the choice of the other initial condition u 0 (x) (the bacteria at time t 0 = 0): We observe that both v 0 (x) and u 0 (x) chosen in tests 1 and 2 verify the compatibility conditions associated to system (2). The results obtained by running the algorithm of §3.2 are presented in Table 2; it points out that the blow-up time in the sense of Definition 2.1 is smaller than the value computed in the sense of Definition 3.1. This is justified by the different behaviors of the two norms, W and L; in fact, Figures 1(e) and 3(a) show that the W -norm grows considerably faster with respect the L-norm. In the same table we can observe that the blow-up time decreases with increasing the corresponding norm of the initial data (see Remark 1). Furthermore, as we briefly discussed in §1, if the initial cell  Table 2. Computing the blow-up times in the sense of the W -norm and the L-norm density is sufficiently large the phenomenon of cell aggregation is produced; the nonlocal chemical interaction dominates diffusion and results in a blow-up of the bacteria distribution, represented by function u. Both examples we are proposing obey this phenomenon; in either case, the contribution of the part of both norms W and L corresponding to u dominates that of v. In particular, for Test 1, this is observed in Figures 1(a), 1(b), 1(c) and 1(d).
On the other hand, let us also specify that the initial function u 0 of Test 1 has one peak in the origin O. The effect of the blow-up is to produce a larger and sharper peak at this same point (see Figure 2); Figure 1(f) exhibits the temporal evolution of u(O), that in this case coincides with that of the L-norm. On the contrary, the initial function u 0 of Test 2 has four peaks in the points P 1 = (x 0 , y 0 ), P 2 = (x 0 , −y 0 ), P 3 = (−x 0 , y 0 ) and P 4 = (−x 0 , −y 0 ) (with x 0 = y 0 = 4 − √ 10), where it achieves the same value. As shown in Figure 4 they coalesce and produce a single peak in the origin O. More exactly, the values of u in these four points decrease while its value at the origin O, which initially represents a local minimum of u 0 , increases and eventually blows up (Figures 3(b), 3(c) and 3(d)).

Remark 2.
In another example we have chosen as initial data a compatible function u 0 (x) presenting four peaks, none at the origin O and one of them considerably larger than the   ; also in this case the values at these points decrease contrary to that at the origin which increases; finally, they coalesce and produce a single peak in O.  Remark 3. These numerical simulations have been obtained by means of the software package FreeFem++ (see [6]). This is a free programming language based in the Finite Element Method, focused on solving partial differential equations. FreeFem++ is implemented in terms of the variational forms of the corresponding problems, so that it is straightforward to address problems involving PDEs (see, for instance, [5] and [18]).

Conclusions.
In this work we have discussed some aspects of a two and three dimensional Keller-Segel system connected to chemotaxis phenomena. Specifically, the blowing up behaviors of its positive solutions have been analyzed. We have started by reviewing theoretical results concerning lower bounds of unbounded solutions in terms of an appropriate norm (the W -norm). In parallel, we have also proposed a numerical resolution algorithm for the same system; it consists of a mixed finite-difference and finite-element method. We have applied this algorithm and given some numerical results for the 2D case; it is shown that with a fixed integration step highly concentrated blowing up solutions can be successfully captured. In this sense, we have computed the blow-up time t W in terms of the W -norm and also t L , i.e. that obtained in terms of the more natural infinity norm (the L-norm). Moreover, an analysis of the temporal evolution of these norms has been carried out; we observe that even though the W -norm increases considerably faster with respect to the L-norm, t W is not so much smaller than t L . Moreover, the analysis of the proposed tests also allows us to conclude that the method is coherent with the expected results: in fact, for instance, the solution depends on that of the corresponding initial function and the more the initial norm (as much W as well L) increases the more the associated blow-up time decreases.