ASYMPTOTIC BEHAVIOR OF SOLUTIONS OF A NONLINEAR DEGENERATE CHEMOTAXIS MODEL

. Pattern formation in various biological systems has been attributed to Turing instabilities in systems of reaction-diﬀusion equations. In this paper, a rigorous mathematical description for the pattern dynamics of aggre-gating regions of biological individuals possessing the property of chemotaxis is presented. We identify a generalized nonlinear degenerate chemotaxis model where a destabilization mechanism may lead to spatially non homogeneous so- lutions. Given any general perturbation of the solution nearby an homogenous steady state, we prove that its nonlinear evolution is dominated by the cor- responding linear dynamics along the ﬁnite number of fastest growing modes. The theoretical results are tested against two diﬀerent numerical results in two dimensions showing an excellent qualitative agreement.


(Communicated by Pierre Magal)
Abstract. Pattern formation in various biological systems has been attributed to Turing instabilities in systems of reaction-diffusion equations. In this paper, a rigorous mathematical description for the pattern dynamics of aggregating regions of biological individuals possessing the property of chemotaxis is presented. We identify a generalized nonlinear degenerate chemotaxis model where a destabilization mechanism may lead to spatially non homogeneous solutions. Given any general perturbation of the solution nearby an homogenous steady state, we prove that its nonlinear evolution is dominated by the corresponding linear dynamics along the finite number of fastest growing modes. The theoretical results are tested against two different numerical results in two dimensions showing an excellent qualitative agreement.
1. Introduction. The term chemotaxis is commonly used to describe the movement of a species or organism which has a tendency to aggregate by moving towards higher concentrations of chemical substances that they themselves produce [18]. It naturally arises in a wide variety of interesting biological settings, including cancer metastasis and bacteria displacement [16]. The Keller-Segel model for chemotaxis is the most widely studied models. It was developed focusing on the cellular slime mold "Dictyostelium Discoideum" [7]. As mentioned in the pioneer work in [21], one cell emits a signal cAMP (cyclic Adenosine Monophosphate); hence, the other cells are chemotactically attracted to this signal and in turn begin to emit cAMP themselves. In mathematical field, asymptotic analysis is a method of describing limiting behavior of solutions. It is a key tool to explore systems of partial differential equations which arise in the mathematical modeling of real-world phenomena. In many parts of qualitative theory of dynamical systems deal with asymptotic properties of solutions (or trajectories) to know what happens with the system under small perturbations after a long period of time. If the system is stable then a small perturbation will result small perturbations in the solution. If not, it will result a motion with large amplitude.
Patterns are solutions of a reaction-diffusion system which are stable in time and stationary inhomogeneous in space, while pattern formation in mathematics is the process that, the spatially homogeneous steady states lose stability under perturbations and therefore inhomogeneous solutions arise (see [22]). The experiments addresses the following question about morphogenesis phenomenon: How an initially spherically symmetric system, such as an embryo in early stages of development, lead to a creature as spherically asymmetric as a person or an animal with spotted coat? One of the first attempts to theoretically explain morphogenesis is by Alan Turing. In [31], Turing proposed that an instability mechanism is the underlying cause of morphogenesis by showing that under certain conditions, such systems which are stable without diffusion are destabilized with diffusion (see [17,27]). In addition to that, precise estimates and behavior are given for pattern formation without an a priori assumption for the smallness of the perturbation later in time as in [29].
The destabilization of homogeneous stationary states may give rise to pattern formation which is a fascinating area of biology with many unanswered questions. Over time, pattern formation in various biological systems has been attributed to such Turing instabilities. Such a mechanism is already present in the well-known Keller-Segel model [21]. It relies on the destabilization mechanism via a chemotactic term of an otherwise stable uniform steady state. This mechanism is also analyzed for the case of the two-equation system in one spatial dimension in [32]. Recently, Guo and Hwang in [14] investigated nonlinear dynamics near an unstable constant equilibrium for the standard Keller-Segel model and their result introduced the mathematical characterization for pattern formation. Then, Shengmao and Cao in [10] treated the classical Keller-Segel model with a source term.
In this paper, we will study a larger class of nonlinear degenerate Keller-Segel systems including a nonlinear chemotaxis term in higher spatial dimensions where a destabilization mechanism near a stable constant equilibrium may lead to pattern formation. We develop and prove the results introduced in the short proceedings [4] where we present a detailed proof that the solution of the nonlinear system behaves asymptotically as spatially non-homogeneous solutions. These solutions prevent blow-up and the global existence has been already proved in [5] due to the maximum principle. Then, the L ∞ −norm of of the species unknown function remains bounded in either finite or infinite time. However, most of the works, including A. M. Turing himself, dealt with linear instability and the nonlinear effect of models in large time have been in question to date from a mathematical point of view. The asymptotic behavior of an attraction-repulsion linear Keller-Segel system has been studied in one dimension (see for instance [20]) relying on the use of Lyapunov function to establish the global existence and large-time behavior of solutions, also the existence of stationary solutions in certain parameter regimes of the same system were established via the method of energy estimates and the phase plane analysis [26,30]. However, the global-in-time existence and boundedness of solutions for a nonlinear Keller-Segel model are proved under some conditions on the nonlinear sensitivity functions and growth source function (see for example [33]).
In general, nonlinear instability is treated with great delicacy for non dissipative systems as in [2] and [12]. Starting in [13], nonlinear instability was established from linear instability for dissipative systems [14], based on a bootstrap lemma. Then, in [8], nonlinear instability was proved with matrix theoretical tools under a sufficient condition that the chemotactic feedback is sufficiently strong. Next, the authors in [15] considered the classical Turing instability over a time scale before blow-up of solutions. Our result has several advantages over an earlier different results along this direction. One of the key ideas is that the qualitative behavior of solutions under perturbations can be analyzed using the linearization near the equilibrium where the question of instability of the nonlinear system is dominated by the corresponding dynamics of its linearization. But, we cannot use the bootstrap lemma in our case due to lack of regularity of solutions of our model. On the other hand, our solutions prevent blow-up and the global existence has been already proved in [5] due to the maximum principle. Then, the L ∞ −norm of species function remains bounded in either finite or infinite time.
This paper is concerned with the asymptotic behavior of solutions of the following degenerate nonlinear parabolic Keller-Segel system, where Ω = T n = n i=1 ]0, π[ is a n-dimensional box for n = 1, 2 or 3. In addition, we assume no-flux conditions for U and V . Therefore, this system of equations is supplemented by the following homogeneous Neumann boundary conditions on where η is the exterior unit normal vector to ∂Ω outward to Ω. The initial conditions on Ω are given by, The species U represents the cell density and V stands for the chemical (external signal) concentration. This system is a modified Keller-Segel model with a nonlinear degenerate diffusion and cross-diffusion laws. In the cell density equation, the diffusive flux is modeling undirected (random) cell migration and the cross-diffusion flux with velocity dependent on the gradient of the signal, is modeling the contribution of chemotaxis. The diffusion coefficient of chemoattractant is denoted by d > 0. The function g(U, V ) = αU − βV describes the rates of production and degradation of the chemical signal (chemoattractant). We assume that χ(0) = 0 and there exists a maximum density of cells U m such that χ(U m ) = 0. The threshold condition has a clear biological interpretation; the cells stop to accumulate at a given point of Ω after their density attains certain threshold value U m , therefore the chemosensitivity χ (U ) vanishes when U tends to U m . This interpretation is sometimes called the volume-filling effect, or prevention of overcrowding. Next, we assume that the density-dependent diffusion coefficient a (U ) degenerates for U = 0 and U = U m . This means that the diffusion vanishes when U attains the threshold U m , and also in the absence of cell-population. In the sequel, we assume U m = 1 without any restriction.
To sum up, the main assumptions about the system (1) The global existence of a weak solution to the system (1)-(3) coupled with Navier-Stokes equation, in the sense of Definition 1.1 below, has been treated in details in [5].
for all ψ 1 and ψ 2 ∈ L 2 0, T ; H 1 (Ω) , where C w 0, T ; L 2 (Ω) denotes the space of continuous functions with values in (a closed ball of) L 2 (Ω) endowed with the weak topology. Theorem 1.2. Suppose that there exists a constant C 0 > 0 such that Then, the system (1)-(3) has a unique global weak solution (U, V ) in the sense of Definition 1.1.
Note that condition (6) may be expressed in an another way given in [24,Proposition 4].

2.
Growing modes in the Keller-Segel model. In this section, we summarize the classical linear instability criterion for the modified Keller-Segel model (1)-(3). When the diffusion terms are ignored, a uniform constant solution forms a nontrivial homogeneous steady state provided In this paper, we study the nonlinear evolution of a perturbation around which satisfies the equivalent system: Then, the system can be written into a matrix form: where, The corresponding linearized Keller-Segel system then takes the form: Let q = (q 1 , ..., q n ) ∈ N n and let . Then, {e q } q∈N n forms a Hilbert basis of the space of functions in L 2 (Ω).
We look for a normal mode to the linear Keller-Segel system (10) of the following form: where r q is a vector depending on q.
Plugging equation (11) into system (10) yields A nontrivial normal mode can be obtained by setting This leads to the following dispersion formula for λ q : Thus we deduce a linear instability by requiring: Due to a non negative discriminant of (12) and det (B) < 0, then one can deduce the existence of two eigenvalues of different sign (see e.g. [28,Chapter II]). Consequently, the characteristic polynomial (12) has at least one positive root λ q .
Obviously, an elementary computation of the discriminant yields: Therefore, we can denote two distinct real roots for all q by The corresponding linearly independent eigenvectors are given by Clearly, one has a U d q 2 + β − χ U α > 0, for q large in equation (13).
Hence, there are only finite numbers of q such that h q 2 < 0 and λ + q > 0. We therefore denote the largest eigenvalue by λ max > 0 and define Moreover, there is one q (possibly more than one) having λ + q = λ max , when we regard λ + q as a function of q 2 . We also denote ν > 0 to be the gap between the λ max and the rest, that is, where Q is the set of q ∈ N n such that h q 2 < 0.
Given an initial perturbation We know that the unique solution of system (10) is given by The aim of this section is the following proposition. Proposition 1. Assume the instability criterion (13) and consider the solution (17) of the linearized Keller-Segel system (10). Then, there exists a constant C 1 depending on α, β, d, a U , χ U such that Proof. We prove the proposition by distinguishing two cases: a) t ≥ 1. By analyzing the quadratic formula (12), we have and it follows also that: Consider equation (16), one has: Later on, we will always denote universal constants by E i (i = 1, 2, . . . ). It is not hard to verify that there exist positive constants E 1 and E 2 , such that: where E 2 = max{ 1 α (E 1 + 1), β α + 1}. It follows from equation (14) and inequality Consequently, From inequalities (19)- (20), and for t ≥ 1, we have For all q ∈ N n , there exists a constant E 4 > 0 such that Therefore, From equation (16), we have It follows from inequality (22) and equation (23) that Therefore, b) For t ≤ 1. It suffices to derive the standard energy estimate of system (10) in L 2 . From the Neumann boundary conditions, we can add the first equation multiplied by u and the second multiplied by κv where κ is a constant to determine. Consequently, We need now to find the constant κ in order to obtain a non-negative integrand in the second integral. For that, using Young inequality and for κ = (χ(U)) 2 da(U ) , one has: Using again Young inequality, it follows that: and the Gronwall inequality implies and hence one can deduce the proposition when t ≤ 1.
3. Main result. The transition from linear into nonlinear instability is very delicate. The key is to control the nonlinear growth for the perturbation by the linear growth rate along the dominant eigenfunction. A bootstrap argument was developed in [10] and [14] to prove the H 2 −regularity of solutions needed to prove nonlinear instability before the possible blow-up time. The latter bootstrap argument cannot be used for our degenerate nonlinear system (9) due to lack of regularity of solutions (L 2 instead of H s ). Moreover, we will not consider any critical time due to the global existence of solutions. Furthermore, many energy estimates, inspired from [5], will be constructed to prove that the dynamics of a general perturbation is characterized by the linear dynamics and that the linear fastest growing modes determine unstable patterns to the full Keller-Segel system (9).
The main result of this paper is given in the upcoming theorem. It gives a mathematical description of pattern formation to our degenerate system (1)-(3). Notice that each initial perturbation can be different from another, which gives rise to the richness of the pattern and the finite number maximal growing modes determines the common characteristics of the pattern.
Theorem 3.1. Assume the instability criterion (13) and let: If the initial perturbation of the steady state [U , V ] is W (x, 0) = W 0 , then its nonlinear evolution W (t, x) (i.e. the solution of system (8) with initial condition W 0 ) satisfies  . We can conclude that for t large enough, the first term e −νt tends to zero and the second is considered small enough for a small perturbation of the initial condition. Consequently, the behavior of the nonlinear solution is similar to a heterogeneous stationary solution, which gives a mathematical description of pattern formation.
Before we start the proof of our main theorem, we will demonstrate a similar version of the Gronwall lemma. then, for almost everywhere t ∈]0, T [.
and using the assumption (27), on has The above equation is equivalent to : and consequently by simple integration, one gets which establish the estimate (28).
Proof of theorem 3.1. We denote by Pw ∈ H 2 (Ω) ∩ L 2 0 (Ω) the unique solution to −∇ · (∇Pw) = w, in Ω, ∇Pw · η = 0, on ∂Ω, where L 2 0 (Ω) = {u ∈ L 2 (Ω) , Ω u dx = 0}. We denote by (−∆) −1 the operator associated to problem (29) defined by : Pw = (−∆) −1 w. We have Consider a nonlinear evolution W = (N 1 , C 1 ) written in term of the steady state and considered as a weak solution of the following problem: and a linear evolution W L = (N 2 , C 2 ) of a perturbation around a steady state as a weak solution of the following problem: The problems (31) and (32) have the same initial data, that is Let us fix T > 0 and consider for (x, t) ∈ Ω × [0, T ], We subtract systems (31) and (32) and therefore one gets, Due to the lack of regularity for solutions of our problem (34), we will construct our estimates for that allow as to estimate the norm Upon multiplying the second equation of system (34) by C with ε < d 2α , one obtains: (36) Next, we multiply the first equation of system (34) by PN , one has: It can be written as: where, Therefore, using assumption (6) with ε < 1 E0 , one has: Moreover, equation (29) implies that Using the linear instability and the boundedness of the density N , one has In addition to that, we have Therefore, Gathering inequalities (37), (38) and (39) for < min( d 2α , 1 E0 , d 2Eχ ), one deduces The last inequality can be written as: where, and Due to Lemma 3.2 and the positivity of k(s), one has Then, using estimate (35), one obtain Due to equality (33), then As we have e λmaxt ≤ e 2λmaxt for λ max t ≥ 0, then Therefore, where E = 2 Consequently, Using estimate (21), we deduce that We know that there is one (or two) q 2 satisfying λ + (q) = λ max . If there is only one q 2 , we denote it by q 2 max and if there are q 2 1 and q 2 2 , we can let q 2 max = max{q 2 1 , q 2 2 }. Therefore, Now, we consider B 2 , From (23) and (24), we get Therefore, From (40), (41) and (42), one can obtain where E = max{E , E + 1}, < min( d 2α , 1 E0 , d 2Eχ ) and K 1 = 2(α+1) ε + 2β. This end the proof of theorem 3.1.
4. Asymptotic behavior of solutions in two dimensions. In this section, we will perform some numerical simulations in order to investigate the asymptotic behavior of the modified Keller-Segel model (1)-(3) and to numerically validate our theoretical results about pattern formation. According to theorem 3.1 and remark 2, it is known that for any initial condition given by equation (25), the asymptotic behavior of the solution of system (1)-(3) is similar to some heterogenous solutions managed by the positive linearly independent eigenvectors given in relation (14).

4.1.
Numerical method. The numerical scheme adopts a Method of Finite Volume approach in which the equations are first discretized in space on an admissible mesh verifying the orthogonality condition (see [9]), and the subsequent system of ODEs is then discretized in time using forward Euler method [11]. Discretization of diffusion terms is performed with a cell-centered finite volume method using the well known TPFA (two points flux approximation) [3,6,19], while the advective term is discretized using an upwind differencing scheme. This choice of finite volume scheme together with the upwind technique is essential to have a scheme that maintains the monotonicity of solutions to our modified system (1)-(3) i.e. to ensure the discrete maximum principle (e.g. see [1]). On the other hand, it is well known that upwind technique preserves the local conservativity of the numerical fluxes [1,25], i.e. the numerical flux is conserved from one discretization cell to its neighbor. In what follows, we consider two different tests to investigate the pattern formation for system (1)-(3). Newton's algorithm is used to approach the solution U n+1 of the nonlinear system defined by the first equation of system (1), this algorithm is coupled with a bigradient method to solve the linear systems arising from the Newton algorithm as well as the linear system given by the second equation of system (1)-(3). Finally, for the upwind technique, we use a numerical flux F defined by Assume that the chemosensitivity χ(U ) given in system (1) involves a chemosensitivity coefficient ζ such that χ(U ) = ζϕ(U ) where ϕ is a certain function verifying assumption (A2). Then, system (1)-(3) depends on many parameters such as the chemosensitivity χ and others involved in the kinetic function (α and β) that play a significant role in production of bifurcations. According to these bifurcations, we determine the parameters of system (1)-(3) for which we build the initial condition (25) such case we can expect the pattern formation as an asymptotic behavior for our modified model. Hereafter, we detail our strategy and determine the relationship between the bifurcation parameters.

4.1.2.
Bifurcation with chemotactic sensitivity ζ, growth rate α, and death rate β. We give first the influence of the chemosensitivity parameter ζ on the pattern formation for system (1)-(3). Generally, there exists a critical value ζ c such that there is no pattern formation if ζ is below this critical value ζ c , while pattern formation can be expected if ζ is somewhere else above the critical value. To determine explicitly the critical value, we fix all parameters of system (1)-(3) except the chemotactic sensitivity ζ. And, from the previous analysis, we know that the bifurcation occurs when condition (13) is fulfilled. Therefore, the critical chemosensitivity ζ c is given by ζ c = a(U )β ϕ(U )α and pattern formation is expected when ζ > ζ c together with the other fixed parameters.
In the same manner, if we fix all parameters of system (1)-(3) except the growth rate α. Then according to condition (13), the critical growth rate α c is given by α c = a(U )β ϕ(U )ζ and pattern formation is expected when α > α c . However, for the bifurcation with the death rate β, we get from condition (13) and after fixing all other variables of system (1)-(3) that the pattern formation is expected only when β is below the critical death rate β c = χ(U )α a(U ) and that there is no pattern formation if β is above the critical value β c since in this case condition (13) becomes not satisfied. We perform all tests on an unstructured triangular mesh of the space domain Ω = (0, 1) × (0, 1) verifying the orthogonality condition (see Fig. 1) and unless stated otherwise, we assume zero-flux boundary conditions. For the steady state defined by equation (7), we fix U = 0.5 and ∆t = 0.01 (time step). satisfying condition (13) (see Fig. 2) is obtained from the zeros 0 and q * 2 of the function h( q 2 ) defined by equation (13), this range of wave numbers is defined by We choose ζ = 10 −3 > ζ c such a way we can expect the pattern formation, next we determine the set of eigenvalues (see Fig. 3) and eigenvectors defined in equation (14), and finally we fix the initial condition such as the sum defined in equation (25). Fig. 4 shows the plot of the initial condition of the function u (x, t) with a small perturbation around zero since for the initial function of U (x, t) solution to system (1)-(3), we consider a perturbation (with an order of magnitude equal to 8 × 10 −2 ) around the steady state U = 0.5. Spatial evolution of u (x, t) and the computed stationary state. Here, we give a comparison between the nonlinear evolution u (x, t) computed by solving numerically the full system (1)-(3) and the stationary state obtained by the nonlinear analysis and generated by the set of eigenvectors having the largest eigenvalue λ max . Fig. 5 shows the spatial evolution of the aforementioned functions, we remark at the beginning that the initial condition leads to some diffusions in the space which gives rise to some aggregations of densities that start the merging process in such  a way to generate spatial patterns which are identical and positioned at the same location of those described by the heterogeneous stationary state. In Fig. 6, we give the shape of the spatial patterns for both nonlinear evolution u (x, t) and the heterogeneous stationary state at t = 1000, we see that the numerical result coincides with the analytical result. Numerically evaluated in the L 2 -norm, the distance between the nonlinear approximation and the numerical solution of the system (1)-(3) is about 2 × 10 −2 . This is in agreement with the main result given in theorem 3.1.  Time evolution of the difference between solutions in L 2 norm. Here, we are interested in the time-tracking of the relative error in L 2 norm between the nonlinear evolution u (x, t) and the heterogeneous stationary state generated by the eigenvectors associated to the maximum eigenvalue λ max . Fig. 7 shows the time evolution of the quantity given by E (t) = W (x,t)−e λmax t q∈Qmax w + q r + q eq(x) .
Note that for numerical calculation, we compute the L 2 −error norm instead of the weakest norm estimated in Theorem 3.1. The graph reveals that the difference between solutions is following an exponential growth curve and that the nonlinear evolution u (x, t) is closely consistent with the heterogeneous stationary state when t ≤ 1100. The results revealed from the graph affirm the consistency of our theoretical result, the nonlinear evolution gives rise to patterns which are installed directly at the same location of the modes defined by the heterogeneous stationary  state and that finite number maximal growing modes determines these patterns. Above t = 1100, we remark that the graph increases exponentially with the time, this aspect looks normal due to the exponential term present in the stationary state obtained by the nonlinear analysis. In addition, when t is large enough the nonlinear evolution u (x, t) cannot exceed 1 (numerically and theoretically) due to the discrete maximum principle verified by the finite volume scheme used in our case. On the other hand, there is no control of evolution for the heterogeneous state which justifies the exponential increase of the difference of L 2 norm between u (x, t) and the heterogeneous state. Test 2. In this numerical test, we fix α = 2, . Plot of h( q 2 ) as a function of q 2 defined by equation (13). When the death rate β decreases below the critical value β c , h( q 2 ) becomes negative for a finite range of unstable wave numbers q 2 marked with rhombi, and pattern formation can be expected.
such a way we can expect the pattern formation (see Fig. 8), and we determine the set of eigenvalues (see Fig. 9) and eigenvectors defined in equation (14), and finally we fix the initial condition such as the sum defined in equation (25).  Spatial evolution of u (x, t) and the computed stationary state. Here, we show a comparison between the nonlinear evolution u (x, t) and the stationary state obtained by the nonlinear analysis and generated by the set of eigenvectors having the largest eigenvalue λ max . Fig. 11 shows the spatial evolution of the aforementioned functions. As for the previous test, we remark that the initial condition leads to some diffusions in the space which gives rise to some aggregations of densities that start the merging process in such a way to generate spatial patterns which are identical and positioned at the same location of those described by the heterogeneous stationary state. In Figure 11. First row from left to right. Nonlinear evolution of the function u (x, t) at t = 10, t = 70, and t = 750. Second row from left to right. Evolution of the heterogeneous stationary solutions at the same moments as for u (x, t). Figure 12. Similarities of patterns between the nonlinear evolution u((x, t) (to the left) and the heterogeneous state (to the right). Fig. 12, we give the shape of the spatial patterns for both nonlinear evolution u (x, t) and the heterogeneous stationary state, we see that the numerical result coincides with the analytical result. Numerically evaluated in the L 2 -norm, the distance between the nonlinear approximation and the numerical solution of the system (1)-(3) is about 1.4 × 10 −2 . This is again in agreement with the main result given in theorem 3.1.

5.
Conclusion. This paper represents a mathematical description of the pattern formation for the full nonlinear degenerate Keller-Segel model. For that, we have studied the asymptotic behavior of stationary solutions and we have proved that the nonlinear evolution of a perturbation around an homogeneous steady state is controlled by the finite number of linear growing modes of this perturbation. This is a clarified quantitative description in order to understand many spatial inhomogeneities and pattern formation processes in a variety of biological developmental situations. Moreover, we have presented two different numerical results in two dimensions to illustrate the effectiveness of our theoretical results presented in this paper and investigate the pattern formation.