Solvability of a Keller-Segel system with signal-dependent sensitivity and essentially sublinear production
Giuseppe Viglialoro Thomas E. Woolley

Spectral analysis of sheared nanoribbons

    Philippe Briet,^a  Hamza Abdou Soimadou^a  
and  David Krejčiřík^b
    



  a)
Aix-Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France;
briet@univ-tln.fr, abdou-soimadou-hamza@etud.univ-tln.fr

  b)
Department of Mathematics, Faculty of Nuclear Sciences and 
Physical Engineering, Czech Technical University in Prague, 
Trojanova 13, 12000 Prague 2, Czechia;
david.krejcirik@fjfi.cvut.cz




26 July 2018
=================================================================================================================================================================================================================================================================================================================================================================

Spectral analysis of sheared nanoribbons

    Philippe Briet,^a  Hamza Abdou Soimadou^a  
and  David Krejčiřík^b
    



  a)
Aix-Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France;
briet@univ-tln.fr, abdou-soimadou-hamza@etud.univ-tln.fr

  b)
Department of Mathematics, Faculty of Nuclear Sciences and 
Physical Engineering, Czech Technical University in Prague, 
Trojanova 13, 12000 Prague 2, Czechia;
david.krejcirik@fjfi.cvut.cz




26 July 2018
=================================================================================================================================================================================================================================================================================================================================================================



Giuseppe Viglialoro^1,♯ Thomas E. Woolley^2
^1Dipartimento di Matematica e Informatica
 Università di Cagliari
 Viale Merello 92, 09123. Cagliari (Italy)
 ^2Cardiff School of Mathematics
   Cardiff University
   Senghennydd Road,  Cardiff, CF24 4AG (United Kingdom)
 


In this paper we consider the zero-flux chemotaxis-system 

    u_t=Δ u-∇· (u χ(v)∇ v)    in  Ω× (0,∞), 
    
    0=Δ v-v+g(u)    in  Ω× (0,∞),

in a smooth and bounded domain Ω of ℝ^2. The chemotactic sensitivity χ is a general nonnegative function from C^1((0,∞)) whilst g, the production of the chemical signal v, belongs to C^1([0,∞)) and satisfies λ_1≤ g(s)≤λ_2(1+s)^β, for all s≥ 0, 0≤β≤1/2 and 0<λ_1≤λ_2. It is established that no chemotactic collapse for the cell distribution u occurs in the sense that any arbitrary nonnegative and sufficiently regular initial data u(x,0) emanates a unique pair of global and uniformly bounded functions  (u,v) which classically solve the corresponding initial-boundary value problem.  Finally, we illustrate the range of dynamics present within the chemotaxis system by means of numerical simulations.




§ INTRODUCTION, MOTIVATIONS AND MAIN RESULT



 §.§ Chemotaxis models: general overview and results

In this paper we  focus our attention on one of the innumerable variants of the landmark models proposed by Keller and Segel (<cit.> and <cit.>), which idealize chemotaxis phenomena (largely spread in the physical and biological sciences), describing situations where the motion of a certain individual cells u=u(x,t) is strongly influenced by the presence of a chemical signal v=v(x,t). Precisely a very general mathematical formulation of these models involves two coupled  partial differential equations and reads

    {[       u_t=  ∇· (A(u,v)∇ u-B(u,v)∇ v)+C(u,v)                               in  Ω× (0,∞),;                           τ v_t= Δ v+E(u,v)                               in  Ω× (0,∞),;                             ∂ u/∂ν=∂ v/∂ν=0                                in  ∂Ω, t>0,; u(x,0)=u_0(x) ≥ 0   and   v(x,0)=v_0(x) ≥ 0                                        x∈Ω, ].

where Ω⊂^n, with n≥ 1, is a bounded domain with smooth boundary, τ∈{0,1} and A, B, C and E are sufficiently regular functions of their arguments.  Additionally, u_0(x) and  v_0(x) are the initial cell and chemical distributions, ∂/∂ν stands for the outward normal derivative on ∂Ω, so that the zero-flux boundary conditions on both u and v indicate the total insulation of the domain.

Strong numerical methods support real experiments and observations indicating that the aforementioned movement may eventually degenerate to aggregation processes, where an uncontrolled gathering of cells at certain spatial locations is perceived as time evolves: the so called chemotactic collapse. In particular, such a coalescence phenomena can be intuitively justified in this way: for B(u,v)>0 the first equation in (<ref>) indicates that while the cell population naturally diffusing (through the law of A) and growing/dying (through the law of C), it is additionally driven by the concentration gradient of the chemical signal in the opposite direction of the diffusion (due to the positivity of B).

By the mathematical point of view, the chemotactic collapse implies that possibly u, in a particular  instant (blow-up time), becomes unbounded in one or more points of its domain. The appearance of this instability may be tied to the size of the initial data (u_0,v_0), the growth rate of the cell distribution u induced by the source term C(u,v), the mutual interplay between the diffusion A(u,v) and the chemotactic sensitivity B(u,v) (also in terms of the space dimension) as well as the production/degradation rate of the chemical signal v given by E(u,v).  

For the linear diffusion case A(u,v)≡ 1, and B(u,v)=χ u, χ>0, in absence of source C(u,v) and with production E(u,v)=-v+u, for τ=1 (parabolic-parabolic case), in <cit.> it is shown that in one-dimensional domains, all the solutions of  (<ref>) are global and uniformly bounded in time, while in the n-dimensional context, with n≥ 2, unbounded solutions to the same problem have been discovered (see, for instance, <cit.> and <cit.>). Similarly, for the parabolic-elliptic case (τ=0), in <cit.> for radial solutions and in <cit.> for non-radial, the authors prove that for n=2 a certain threshold value given by the  product between the chemosensitivity χ and the initial mass ∫_Ω u_0 decides whether the solution can blow up in finite time  or exists for all time t>0. 

By virtue of this, let us observe that for positive chemical and cell distributions, the expression E(u,v)=-v+u  in (<ref>) shows how an increasing of the cells favours a secretion of the signal, which depending on the expression of the chemosensitivity B(u,v) might strongly contrast the smoothing and equilibrating effect of the  diffusion A(u,v). In this sense an abundant literature concerning existence and properties of global, uniformly bounded or  blow-up solutions, is available; for a complete picture, we suggest the largely cited contributions <cit.> and <cit.> where, inter alia, reviews of various models about Keller-Segel-type systems are presented and analysed.

Finally, exactly in order to provide a more comprehensive picture of the general situation tied to this balance between the destabilizing and stabilizing effects (A versus B), the presence of an absorptive  logistic-type source of the type C(u,v)≃ ku-μ u^δ,  for k ∈, μ>0 and δ>1 in (<ref>)  may, or may not, have a certain relevance on the dynamics of the system. We mention the papers  <cit.> and <cit.> (and references therein), which deal with existence, blow-up and properties of solutions to  (<ref>), for both fully parabolic and parabolic-elliptic versions, and as above for the linear diffusion case A(u,v)≡ 1, and B(u,v)=χ u, χ>0.


 §.§ Some inspiring results related to our research: presentation of the main theorem

In this investigation the source C(u,v) is taken nil, whereas the sensitivity B(u,v) = uχ(v) is such that χ belongs to a general family of positive functions. In particular,

    B(u,v)=χ_0/vu,   v>0,  with some  χ_0>0,
 
 employed in the so-called Weber-Fechner law expressing the relation between the actual change in the stimulus and the perceived change (<cit.>), is a prototype of such functions, and it presents a singularity at v=0 which is the main source of both technical and numerical difficulties. Moreover, and as far as we know this is the novelty of our contribution  (at least in this context), we take the function E(u,v) in such a way that, basically, the corresponding reproduction rate implies a lower increasing of the signal  than that supplied by the linear one E(u,v)=-v+u. As a consequence of this, it could be expected that lessening the impact of high values of the cell distribution on the production of the chemical signal may enforce global existence of solutions. To be precise, what could we expect, for instance,  if the chemical were secreted with a sublinear rate depending on the bacterial density? 

Before giving a first answer to this question,  we want to motivate this work starting with an overview of previous achievements regarding some variants of systems like (<ref>), defined in n-dimensional domains. In particular, we believe that the coming contributions, all characterised by the presence of singular chemo-sensitivity and expression for E given by E(u,v)=-v+u, deserve to be discussed since they put into perspective and also inspire our current investigation: 


  * for τ=1, A(u,v)=1, B(u,v)=uχ(v), with χ(v)=χ_0/v, and C(u,v)=0, global existence of weak solutions under the assumption 0<χ_0<n/2 is proved (<cit.>);

  * for τ=0, A(u,v)=1, B(u,v)=uχ(v), with  χ(v)=χ_0/v, and C(u,v)=0, uniform boundedness and blow-up of radial solutions are positive addressed in <cit.>; more exactly, solutions are global and remain bounded when either n≥ 3 and 0<χ_0<n/n-2 or n=2 and χ_0>0 is arbitrary, whilst for n≥ 3, 0<χ_0<n/n-2 and ∫_Ω u_0|x| is sufficiently small, the solution blows up in finite time; 

  * for τ=0, A(u,v)=1, B(u,v)=uχ(v), for 0<χ(v)≤χ_0/v^k, with k≥ 1 and v>0, and C(u,v)=0,  global existence and uniform boundedness of classical non-radial solutions are discussed in <cit.>, where it is shown that the system possesses a unique global classical solution, uniformly bounded if 0<χ_0<2/n (k=1) and 0<χ_0<2/nk^k/(k-1)^k-1γ^k-1 (k>1), γ>0 being a constant depending on the data.

  * for τ=0, A(u,v)=1, B(u,v)=uχ(v), with  χ(v)=χ/v, and C(u,v)=ru-μ u^2, r∈, χ,μ>0, it is proved that in two-dimensional domains the logistic kinetics ensures global existence of classical solutions even for arbitrary large χ and any μ>0. Additionally, it is shown that if r is larger with respect to some expression of χ such solutions are also bounded (<cit.>); 

  * for τ=1, A(u,v)=1, B(u,v)=uχ(v), with χ(v)=χ/v, and C(u,v)=0, in <cit.> uniform boundedness of global classical solutions is shown in the two-dimensional setting and for χ∈(0,χ_0), for some χ_0>1. 

  * for τ=0, A(u,v)=(u+1)^m-1, B(u,v)=u(u+1)^α-1χ(v), with 0<χ(v)≤χ_0/v^k, k≥ 1 and v>0, and C(u,v)=0, for m,α∈ such that α≤max{m,m+1/2}, global existence and uniform boundedness of classical solutions are proved, provided some smallness assumptions on χ_0 are satisfied (<cit.>).
  
In accordance with these premises, we intend to enhance the knowledge of the mathematical analysis of general chemotaxis-systems, by studying this problem

    u_ t=Δ u-∇· (uχ(v)∇ v)    in  Ω× (0,∞), 
    
    0=Δ v-v+g(u)    in  Ω× (0,∞),
    ∂ u/∂ν=∂ v/∂ν=0    on  ∂Ω× (0,∞), 
    
    u(x,0)=u_0(x) ≥ 0    x∈Ω,

where Ω is a smooth and bounded domain of ℝ^2 and u_0(x)=u(x,0) is a nonnegative function such that 0≢u_0∈ C^0(Ω̅). 


Additionally, 0<χ∈ C^1((0,∞)) and g∈ C^1([0,∞)) satisfies this essentially sublinear growth (see Remark <ref> below):

    0<λ_1 ≤ g(s)≤λ_2(1+s)^β  for   s≥ 0,    0≤β≤1/2  and   0<λ_1≤λ_2.

We prove the existence and uniqueness of a global and bounded classical solution to problem (<ref>), and precisely we show that  the high sublinear action induced on v by g exerts a certain smoothing effect on u and it is sufficient to prevent  δ-singularities formation for the cell distribution u, even for widely large initial cells' density u_0(x) or strong sensitivity effects.

This conclusion is mathematically formulated in our main theorem:


Let Ω be a smooth and bounded domain of ℝ^2, 0<χ∈ C^1((0,∞)) and g∈ C^1([0,∞)) a function satisfying  (<ref>). Then for any nonnegative initial data 0≢u_0∈ C^0(Ω̅), problem (<ref>) admits a unique global classical solution (u, v). Moreover, both u and v are bounded in Ω× (0,∞).


For the chemotaxis model (<ref>) with linear production, i.e. g(u)=u, it is seen from the conservation of the total cell mass (see Lemma <ref>), that is ∫_Ω u=∫_Ω u_0, that also the same property for the chemical v holds; indeed, by integrating over Ω the second equation we have ‖ v‖_L^1(Ω)=∫_Ω g(u)=∫_Ω u=∫_Ω u_0>0 throughout the time. Conversely,  for chemotaxis models with sublinear production, g(u)=u^ϑ with 0<ϑ<1,  this is no longer true. In fact, Hlder's inequality implies ‖ v‖_L^1(Ω)=∫_Ω g(u)=∫_Ω u^ϑ≤ (∫_Ω u_0)^ϑ |Ω|^1-ϑ, which does not exclude the possibility of vanishing for ∫_Ω v at some time. Since in order to save the chemosensitivity χ=χ(v) from  singularities  we have to avoid this last scenario (see Lemma <ref>),  we assume   g(s)≥λ_1>0, for s>0. Subsequently, for the production source g restricted to grow no faster than u^β at infinity (essentially sublinear growth), i.e. g(s) ≤λ_2(1+ s)^β for all s>0 and some λ_2>0,  the assumption λ_2≥λ_1>0 makes consistent the lower and upper bounds in (<ref>).

The rest of the paper is structured as follows. First, in <ref>, we collect some necessary and preparatory material, then, in <ref>, we prove the local existence and uniqueness of a classical solution to (<ref>) and some of its properties. Successively, in <ref>, we establish how to ensure globability and boundedness of local solutions using their L^p-boundedness. Such a bound is derived in <ref>, which represents the main part of this report and that concludes with the proof of Theorem <ref>. Finally, the theoretical results presented here are investigated numerically in <ref>, where simulations are used to detect critical exponents for β which
delineate regions where different asymptotic behaviours of solutions to the same system (<ref>) may manifest. 



§ PRELIMINARIES AND AUXILIARY TOOLS

The coming results are supportive in the proof of the main theorem of this paper. To be precise, we mainly summarize and derive some general functional inequalities, also tied to elliptic regularity theory.

Let us first recall a special case of the well-known Gagliardo-Nirenberg inequality which will be used through the paper to prove the main theorem. 
 (Gagliardo-Nirenberg inequality)
Let Ω be a smooth and bounded domain of ^2. Then there is a constant C_GN>0 such that the following inequality holds: 
With 𝔮,𝔰∈[1,2], 𝔭∈ [2,4] and θ=1-𝔮/𝔭∈[0,1),  

    f _L^𝔭(Ω)≤ C_GN ( ∇ f _L^2(Ω)^θ f _L^𝔮(Ω)^1 - θ+   f _L^𝔰(Ω))

is satisfied for all f∈ L^𝔮(Ω) with ∇ f∈ L^2(Ω),  

See <cit.>.


The following lemma is fundamental in our computations and its validity is restricted to two-dimensional settings, which are those where our main problem is studied. Its proof is, fundamentally, a reformulation of <cit.> which we adapted to our presentation in order to make the present article more self-contained.

Let Ω be a smooth and bounded domain of ^2.  Then there exists a positive constant Ĉ such that for all p∈(1,2) and f∈ C^2(Ω̅), with ∂ f/∂ν=0 on ∂Ω, holds that

    ‖∇ f ‖_L^2p+2(Ω)≤Ĉ‖ -Δ f + f‖^1/2_L^p+1(Ω)‖ f ‖^1/2_L^2(Ω).


Given ψ∈ C^2(Ω̅), we apply (<ref>) with 𝔭=4 and 𝔮=2 to obtain

    ψ_L^4(Ω)^4 ≤ c_1  ∇ψ_L^2(Ω)^2ψ_L^2(Ω)^2+  c_1ψ_L^2(Ω)^4,

where c_1=(2C_GN)^4, having also made use of  

    (a+b)^α≤ 2^α(a^α+b^α)  for any   a,b≥ 0, α>0.

In addition, for any φ∈ C^1(Ω̅) we have that  for all p∈ (1,2) the function ψ=|φ|^p+1/2 belongs to W^1,2(Ω) and is such that |∇ψ|=p+1/2|φ|^p-1/2|∇φ|; consequently, inequality (<ref>)  explicitly reads

    ∫_Ω |φ|^2(p+1)≤  c_1(p+1)^2/4(∫_Ω |φ|^p-1|∇φ|^2)(∫_Ω |φ|^p+1)+c_1(∫_Ω |φ|^p+1)^2.

Now, the Hlder inequality enables us to get

    ∫_Ω |φ|^p-1|∇φ|^2 ≤(∫_Ω |∇φ|^p+1)^2/p+1(∫_Ω |φ|^p+1)^p-1/p+1,

and

    (∫_Ω |φ|^p+1)^2p/p+1 =(∫_Ω |φ|^p^2-1/p|φ|^p+1/p)^2p/p+1≤(∫_Ω | φ|^2(p+1))^p-1/p+1∫_Ω |φ|^2.

Further, by inserting these last two relations into (<ref>), a proper decomposition and the inequality (p+1)^2>4, valid for all 1<p<2, infer

    ∫_Ω |φ|^2(p+1)   ≤  c_1(p+1)^2/4[(∫_Ω |∇φ|^p+1)^2/p+1+(∫_Ω |φ|^p+1)^2/p+1]  
         ×(∫_Ω |φ|^p+1)^2p/p+1,

so that algebraic manipulations yield 

    ‖φ‖_L^2(p+1)(Ω)^4   =(∫_Ω |φ|^2(p+1))^2/p+1 =(∫_Ω |φ|^2(p+1))^1-p-1/p+1
       ≤  c_1(p+1)^2/4[(∫_Ω |φ|^p+1)^2/p+1
     +(∫_Ω |∇φ|^p+1)^2/p+1]∫_Ω |φ|^2
       ≤ c_2[(∫_Ω |φ|^p+1+∫_Ω |∇φ|^p+1)]^2/p+1∫_Ω |φ|^2
       
    =c_2‖φ‖_W^1,p+1(Ω)^2‖φ‖_L^2(Ω)^2,

where we used, in the last step, a^γ+b^γ≤ 2^1-γ(a+b)^γ, for all a,b≥ 0 and 0<γ≤ 1  and set c_2= c_1(p+1)^2/42^p-1/p+1.

On the other hand, let C>0 the constant from <cit.>; for f as in our hypothesis and p∈(1,2) we can rely on such inequality and obtain

    ‖ f ‖_W^2,p+1(Ω)≤ C ‖ -Δ f + f ‖_L^p+1(Ω).

In addition, the same properties of f also allow us to set φ=|∇ f| in (<ref>) so to have for c_3=c_2^1/4

    ‖∇ f ‖_L^2(p+1)(Ω)≤ c_3 ‖ f ‖_W^2,p+1(Ω)^1/2‖∇ f ‖_L^2(Ω)^1/2,

where we considered (<ref>) and the fact that 

    ‖∇ f‖_W^1,p+1(Ω)^p+1   =
    ∫_Ω |∇ f|^p+1+∫_Ω|∇  |∇ f||^p+1
       ≤∫_Ω |f|+∫_Ω |∇ f|^p+1+∫_Ω|∇  |∇ f||^p+1≤‖  f‖_W^2,p+1(Ω)^p+1.

We conclude the proof by combining (<ref>) and (<ref>), being Ĉ=c_3√(C).


Conforming to the comments in Remark <ref>, in the next result we will establish a lower bound for the second component of solutions to the parabolic-elliptic Keller-Segel system (<ref>). In particular we derive a quantitative estimate on positivity of solutions to the Neumann problem for the Helmholtz equation with nonnegative inhomogeneity having given norm in L^1(Ω).

For any n∈, let Ω be a smooth and bounded domain of ℝ^n and 0 ≢w∈ C^0(Ω̅) a nonnegative function.  If f∈ C ^2(Ω̅) is a solution of 

    -Δ f + f=w    in  Ω,
    ∂ f/∂ν=0   on  ∂Ω,

then there exists some positive constant η such that 

    f≥η∫_Ω w>0  in  Ω.


The proof is a consequence of the the positivity of the Green function to the Helmholtz equation.


We will also make use of the following elementary proof:

Let y be a positive real number verifying y ≤ k(y^l+1) for some k>0 and 0<l<1. Then y≤max{1,(2k)^1/1-l}.

 Since for y≤ 1 there is nothing left to show, we suppose y≥ 1. Then y^l ≥ 1 so that y ≤ k(y^l+1) ≤ 2k y^l and hence y≤  (2k)^1/1-l.




§ EXISTENCE OF LOCAL-IN-TIME SOLUTIONS AND MAIN PROPERTIES

We open this section with a lemma concerning the local-in-time existence of classical solutions (u,v) to system (<ref>). The proof is developed by adapting well-established methods involving an appropriate fixed point framework and standard parabolic and elliptic regularity results. Through the same lemma we also are able to achieve an important conservation of mass property and, relying on Lemma <ref>, a crucial uniform-in-time estimate for the component v of the solution, which ensures the uniform bound of any signal-dependent sensitivity χ=χ(v) taken from C^1((0,∞)).
  
For any n∈, let Ω be a smooth and bounded domain of ℝ^n, 0<χ∈ C^1((0,∞)) and and g∈ C^1([0,∞)) a function satisfying  (<ref>). Then for any nonnegative initial data 0≢u_0∈ C^0(Ω̅), problem (<ref>) admits a unique local-in-time classical solution 

    (u,v)∈ (C^0([0,T_max);C^0(Ω))∩ C^2,1(Ω̅× (0,T_max)))^2,

where T_max∈(0,∞], denoting the maximal existence time, is such that if T_max<∞ necessarily

    lim sup_t↗ T_max‖ u (·,t)‖_L^∞(Ω)=∞.

Moreover, for some γ,Γ>0 we have for all (x,t) in Ω× (0,T_max)

    u≥ 0,   γ≤ v ≤Γ  and  ‖χ(v) ‖_L^∞(Ω× (0,T_max))≤‖χ‖_L^∞([γ,Γ]),

and also

    ∫_Ω u (· ,t)  = m=∫_Ω u_0>0  for all    t∈ (0,T_max).


Existence. For any T∈ (0,1), 0≢u_0∈ C^0(Ω̅) nonnegative,  and R:=‖ u_0‖_L^∞(Ω)+1, let us consider the Banach space X:=C^0(Ω̅× [0,T]) and its closed subset

    S:={0≤ u ∈ X |‖ u (·, t)‖_L^∞(Ω)≤ R   for all    t∈[0,T]
    }.
  
For û∈ S, let v be the solution of 

    -Δ v+v=g(û)    in  Ω× (0,T),
    ∂ v/∂ν=0    on  ∂Ω× (0,T),

and, in turn, let u be the solution of 

    u_t-Δ u=∇· (u χ(v)∇ v)    in  Ω× (0,T), 
    ∂ u/∂ν=0    in  ∂Ω× (0,T), 
    
    u(x,0)=u_0(x) ≥ 0     x∈Ω.

In agreement with these statements, we shall show that for appropriate small T,  Φ: S→ S defined by Φ(û)=u is a compact map such that Φ(S)⊂ S; subsequently, due to the convexity of S, the Schauder fixed point theorem ensures the existence of  u∈ S such that Φ(u)=u.

First, we observe that for a certain fixed û∈ S well known elliptic regularity results, in conjunction with Morrey's theorem (<cit.>), infer a unique solution v(·,t) to problem (<ref>) in the space C^1,δ(Ω), for all δ∈ (0,1); this, in particular, implies that  ∇ v ∈ L^∞(Ω) for all t∈(0,T). Again continuing on the property of the solution v, because û∈ S, û is nonnegative so that from (<ref>) we have that g(û) is well defined and moreover g(û)≥λ_1>0. In this way, an application of Lemma <ref> to  problem (<ref>), together with the definition of S, leads to

    v(x, t) ≥η·∫_Ω g(û)=γ:=ηλ_1|Ω|>0   (x,t)∈Ω× (0,T).

Moreover, besides v(x,t)≥γ, the elliptic maximum principle, and again (<ref>) and the definition of S, provide

    v(x,t)≤sup_Ω g(û)≤Γ :=λ_2(1+ R)^β   (x,t)∈Ω× (0,T),

so that  on account of χ∈ C^1((0,∞)) and γ≤ s:=v(x,t) ≤Γ, with (x,t)∈Ω× (0,T), χ(s) is also from L^∞([γ,Γ]) and in the specific we have

    ‖χ(v) ‖_L^∞(Ω× (0,T))≤‖χ‖_L^∞([γ,Γ]).

Subsequently, for some positive constant c (which until the end of this proof might change line by line), using û∈ S and the gained bounds for ∇ v and χ(v), <cit.> applied to problem (<ref>) implies that u∈ C^δ_1,δ_1/2(Ω× (0,T)), for some δ_1 ∈ (0,1). Hence, 

    |u(x_1, t_1)-u(x_2,t_2)|≤ c (|x_1-x_2|^δ_1+|t_1-t_2|^δ_1/2)  for all   x_1,x_2∈Ω,   t_1,t_2∈ (0,T),

that is

    u(·, t)≤ u_0(·)+ct^δ_1/2  for all    t∈ (0,T).

Thereafter

    max_t∈[0,T]‖ u(·,t)‖_L^∞(Ω)≤‖ u_0‖_L^∞(Ω)+cT^δ_1/2,

and subsequently for T< c^-δ_1/2 we also deduce that 

    ‖ u(·,t)‖_L^∞(Ω)≤‖ u_0‖_L^∞(Ω)+1=R   for all    t∈ (0,T).

Additionally, u≡ 0 is a subsolution of the first equation in (<ref>) so that the parabolic comparison principle warrants the nonnegativity of u; hence Φ maps S into itself, compactly since C^δ_1,δ_1/2(Ω× (0,T)) ↪ X. Let u be a fixed point of Φ; by employing the elliptic and parabolic regularity theory to problems  (<ref>)  and (<ref>) (explicitly <cit.> and <cit.>, respectively), we have v∈ C^2+δ_1,δ_1/2(Ω̅,[τ,T]) and hence u∈ C^2+δ_1,1+δ_1/2(Ω̅× [τ,T]), for any τ∈ (0,T). Moreover, by standard bootstrap arguments the solution may be prolonged in the interval  [0, T_max), with T_max≤∞, T_max being finite if and only if (<ref>) holds. In this way, the gained nonnegativity of u(·,t) in [0,T], lower and upper estimates for v(·,t) in [0,T] (see  (<ref>) and (<ref>)) and bound for χ(v(·,t)) in [0,T] (see (<ref>)) remain preserved up to T_max, exactly as claimed in (<ref>). Finally, an integration over Ω in the first equation of (<ref>) and the no-flux boundary conditions on both u and v give ∫_Ω u(·,t)=∫_ Ω u_0 for all t∈ (0,T_max), and (<ref>) is also justified.

Uniqueness.  By absurdity let  (u_1,v_1) and (u_2,v_2) be two nonnegative different classical solutions of (<ref>) in Ω× (0,T_max) with the same initial data u_1(·,0)=u_2(·,0). In such circumstances, using the equation for u in (<ref>), we have

    1/2d/d t∫_Ω (u_1-u_2)^2   +∫_Ω |∇ (u_1-u_2)|^2
       
    =∫_Ω (u_1χ(v_1)∇ v_1-u_2χ(v_2)∇ v_2)·∇ (u_1-u_2)
       ≤1/2∫_Ω |(u_1χ(v_1)∇ v_1-u_2χ(v_2)∇ v_2)|^2
          +1/2∫_Ω |∇ (u_1-u_2)|^2    for all   t∈(0,T_max),

and of course 

    ∫_Ω (u_1(·,0)-u_2(·,0))^2=0.

Now, for all t∈ (0,T_0) with T_0<T_max we set

    s_1=s_1(T_0):=min{‖ u_1 ‖_L^∞(Ω× (0,T_0)),‖ u_2 ‖_L^∞(Ω× (0,T_0))}, 
       s_2=s_2(T_0):=max{‖ u_1 ‖_L^∞(Ω× (0,T_0)),‖ u_2 ‖_L^∞(Ω× (0,T_0))},

and the Mean Value Theorem applied to the function s↦ g(s) in the interval [s_1,s_2] infers g(s_1)-g(s_2)=g'(s̅)(s_1-s_2) for some s̅∈ (s_1,s_2). 
In light of this, through Young's inequality, the second equation of (<ref>) provides some positive C_1 depending on T_0 such that on (0,T_0)

    ∫_Ω|∇ (v_1-v_2)|^2   =-∫_Ω (v_1-v_2)^2+∫_Ω (g(u_1)-g(u_2))(v_1-v_2) 
       
    =-∫_Ω (v_1-v_2)^2+C_1∫_Ω (u_1-u_2)(v_1-v_2) 
       ≤C_1^2/2∫_Ω (u_1-u_2)^2-1/2∫_Ω (v_1-v_2)^2.

Additionally, the same elliptic and parabolic regularity results previously used allow us to find some C_2=C_2(T_0)>0 such that |∇ v_1|≤ C_2, |∇ v_2|≤ C_2 and u_2≤ C_2 on Ω× (0,T_0); subsequently, by virtue of the boundedness of χ in (<ref>) and Hlder's inequality,  some manipulations lead to

    1/2∫_Ω   | u_1χ(v_1)∇ v_1-u_2χ(v_2)∇ v_2 |^2≤3/2∫_Ω| (u_1-u_2)^2χ(v_1)^2|∇ v_1 |^2
          +3/2∫_Ω u_2^2χ(v_2)^2|∇ (v_1-v_2) |^2 +
    3/2∫_Ω u_2^2(χ(v_1)-χ(v_2))^2(v_1-v_2)^2|∇ (v_2) |^2
       
    
    ≤3/2‖χ‖_L^∞([γ,R])^2C_2^2 ∫_Ω (u_1-u_2)^2+3/2‖χ‖_L^∞([γ,R])^2C_2^2∫_Ω|∇ (v_1-v_2)|^2
          +3/2  C_2^4‖χ' ‖_L^∞([γ,R])^2∫_Ω  (v_1-v_2)^2  for all   t∈ (0,T_0),
 
where we also applied the Mean Value theorem to the function s↦χ(s), with s∈[s_1,s_2], and used the  inequality  (a+b+c)^2≤  3 (a^2+b^2+c^2), valid for all a,b,c ∈.
Hence, we can write for all t∈(0,T_max) 

    1/2∫_Ω   |(u_1χ(v_1)∇ v_1-u_2χ(v_2)∇ v_2)|^2≤
          C_3 (1/2∫_Ω (u_1-u_2)^2+∫_Ω|∇ (v_1-v_2)|^2 +1/2∫_Ω  (v_1-v_2)^2),

where C_3=C_3(T_0)= 3C_2^2max{‖χ‖_L^∞([γ,R])^2,‖χ' ‖_L^∞([γ,R])^2C_2^2}.












Finally, coming back to (<ref>), and plugging in it (<ref>)  and (<ref>) we arrive at



this initial problem 

    d/d tℱ≤ (C_3+C_3C_1^2) ℱ   t∈ (0,T_0),   ℱ(0)=0,

where ℱ(t):=∫_Ω (u_1-u_2)^2 and where (<ref>) has been also considered. Since (<ref>) admits the unique solution ℱ≡ 0 on (0,T_0), due to the arbitrary of T_0, we attain u_1=u_2 on (0,T_max) and hence, by using again (<ref>), also v_1=v_2  on (0,T_max).




§ FROM LOCAL TO GLOBAL-IN-TIME AND BOUNDED SOLUTIONS

The forthcoming important result shows how to achieve  uniform-in-time boundedness of solutions from their L^p-boundedness, for some suitable p>1. As a consequence, in order to establish our main theorem, it will be successively sufficient the derivation of such a L^p bound. 

Under the assumptions of Lemma <ref>, let  (u,v)  be the local-in-time classical solution  of problem (<ref>). If for some n/2<p<n the u-component belongs  to L^∞((0,T_max);L^p(Ω)), then (u,v) is global in time, i.e. T_max=∞, and moreover both u and v are bounded in Ω× (0,∞). 

For u∈ L^∞((0,T_max);L^p(Ω)), taking into consideration assumption (<ref>) on g, we have that 

    ∫_Ω g(u)^p≤λ_2^p∫_Ω (1+u)^β p≤λ_2^p ∫_Ω (1+u)^p  for all   t∈(0,T_max),

so that also g(u)∈ L^∞((0,T_max);L^p(Ω)). Consequently, standard elliptic regularity results applied to the second equation of (<ref>) warrant v ∈ L^∞((0,T_max);W^2,p(Ω)) and hence ∇ v∈  L^∞((0,T_max);W^1,p(Ω)), and finally Sobolev embedding theorems give v∈  L^∞((0,T_max);C^[2-(n/p)](Ω̅)) and ∇ v∈  L^∞((0,T_max);L^q(Ω)) for all n<q<p^*:=np/n-p. In particular, for some positive constant C_q we have that 

    ‖  v (·, t)‖_L^q(Ω)+‖∇ v (·, t)‖_L^q(Ω)≤ C_q  for all   t∈(0,T_max),

and in addition, because q>n implies W^1,q(Ω)  L^∞(Ω), we also attain v∈ L^∞(Ω× (0,T_max)).

As far as u is concerned, for any x∈Ω and t ∈ (0,T_max) we set t_0:=max{0,t-1} so that  the representation formula for u yields

    u (·,t)    ≤ e^(t-t_0)Δu(·,t_0)-∫_t_0^t e^(t-s)Δ∇· (u (·,s) χ(·,s)∇ v (·,s))ds
        
    =: u_1(·,t)+u_2(·,t).

Here, we invoke known smoothing estimates for the Neumann heat semigroup (see  <cit.>) which warrant the existence of positive constants C_s and Ĉ_S  such that for all t>0 and f∈ L^1(Ω)
 
    ‖  e^tΔ f  ‖_L^∞(Ω)≤ C_S (1+t^-n/2) ‖ f‖_L^1(Ω),
 
 and for all t>0, r>1 and f ∈ C^1(Ω̅) with ∂ f/∂ν=0 on ∂Ω
 
    ‖  e^tΔ∇·∇ f  ‖_L^∞(Ω)≤Ĉ_S (1+t^-1/2-n/2r)e^-λ_1 t‖∇ f‖_L^r(Ω),

 λ_1>0 denoting the first nonzero eigenvalue of -Δ in Ω under Neumann boundary conditions. 

Subsequently, if  t≤ 1 and hence t_0=0 we achieve from the parabolic maximum principle 

    ‖ u_1(·,t)  ‖_L^∞(Ω)≤‖ u_0 ‖_L^∞(Ω)  for all   t≤ 1.

Conversely, for all t>1 and hence t-t_0=1, an application of  (<ref>) with f=u(·,t_0) and bound (<ref>) infer for all 1< t < T_max this estimate

    ‖ u_ 1(·,t)  ‖_L^∞(Ω)≤ C_S (1+(t-t_0)^-n/2)‖ u(·,t_0) ‖_L^1(Ω)≤ 2 mC_S.

Furthermore, for any n<r<q we apply (<ref>) with ∇ f=u (·,t) χ(·,t)∇ v (·,t) and arrive for t∈ (0,T_max) at

    ‖ u_2(·,t) ‖_L^∞(Ω)≤∫_t_0^t ‖ e^(t-s)Δ∇· (u (·,s) χ(·,s)∇ v (·,s)) ‖_L^∞(Ω)ds
         ≤‖χ‖_Ĉ_S ∫_t_0^t (1+(t-s))^-1/2-n/2r‖ u (·,s) ∇ v (·,s) ‖_L^r(Ω) ds,

where once we invoked the boundedness of χ given in (<ref>). 

Now, for any given t'∈ (0,T_max) we consider the function defined by

    A(t'):=sup_t∈ (0,t')‖ u (·, t)‖_L^∞(Ω),

which is bounded in view of the properties of u. Hence, an interpolation inequality, the assumption u∈ L^∞((0,T_max);L^p(Ω)) (i.e. ∫_Ω u^p≤ C_p for some C_p>0 and all t∈(0,T_max)) and (<ref>) entail for all s ∈ (0,t')

    ‖ u (·,s) ∇ v (·,s) ‖_L^r(Ω)   ≤‖ u (·,s) ‖_L^rq/q-r(Ω)‖∇ v (·,s) ‖_L^q(Ω)
       ≤‖ u (·,s) ‖_L^∞(Ω)^1-p(q-r)/rq‖ u (·,s) ‖_L^p(Ω)^p(q-r)/rq‖∇ v (·,s) ‖_L^q(Ω)
       ≤ A^1-p(q-r)/rq (t')C_p^q-r/rqC_q=c_4A^l (t'),

with c_4=C_p^q-r/rqC_q and 0<l=1-p(q-r)/rq <1. Subsequently, in light of the estimate now gained for u (·,s) ∇ v (·,s) and the relation t-t_0≤ 1, bound (<ref>) reads

    ‖ u_2(·, t) ‖_L^∞(Ω)   ≤2r/r-n‖χ‖_Ĉ_S c_4 (2^1/2-n/2r-1)A^l (t')=:c_5A^l(t').

From expression (<ref>), by collecting  (<ref>)-(<ref>) and (<ref>) we infer

    ‖ u (·,t) ‖_L^∞(Ω)≤ c_6(A^l(t')+1)  for all   t ∈  (0,t'),

where c_6=max{max{||u_0||_L^∞(Ω),2mC_S},c_5}. Therefore recalling (<ref>) 

    sup_t∈ (0,t')‖ u (·, t)‖_L^∞(Ω)=:A(t')≤ c_6(A^l(t')+1)  for all   t'∈ (0,T_max),

which through Lemma <ref> yields this bound for u:

    ‖ u (·,t) ‖_L^∞(Ω)≤max{1,(2c_6)^1/1-l}=:ĉ  for all   t∈ (0,t').

Since ĉ>0 is time-independent and t'∈(0,T_max) is arbitrary, the above uniform bound for u holds up T_max. Hence the extensibility criterion (<ref>) of Lemma <ref> shows that T_max = ∞ and that both u and v are bounded in Ω× (0,∞). 




§ A PRIORI ESTIMATES AND PROOF OF THE MAIN RESULT

In this section we shall gain some uniform bound for u, by deriving an upper bound for ‖ u‖_L^p(Ω), with p sufficiently large and on the whole interval (0,T_max), which is given by a positive and time independent constant. This is attained by constructing an absorptive differential inequality for t↦∫_Ω u^p and using comparison principles, exactly as specified in this sequel of lemmas.
 
Under the assumptions of Lemma <ref>, let (u,v)  be the local-in-time classical solution of problem (<ref>). Then there exists a positive constant M such that 

    ∫_Ω|∇ v|^2≤ M  for all   t∈ (0,T_max).


Testing the second equation of  (<ref>) by v, we obtain

    -∫_Ω v Δ v=-∫_Ω v^2+∫_Ω g(u)v   for all   t∈ (0,T_max),

so that integrations by part, the mass conservation property (<ref>) and the Young inequality (recall also (<ref>)) lead for all t ∈ (0, T_max) to

    ∫_Ω |∇ v|^2   ≤ -∫_Ω v^2+∫_Ω v^2+λ_2^2/4∫_Ω (1+u)^2β≤λ_2^2/4(|Ω|+m)=:M.



  
Under the assumptions of Lemma <ref>, let (u,v) be the local-in-time classical solution  of problem (<ref>). Then, for any p>1

    d/dt∫_Ω u^p   +pp-1/2∫_Ω u^p-2|∇ u|^2 ≤
       
    +pp-1/2‖χ‖^2_L^∞([γ,Γ])∫_Ω u^p|∇ v|^2  on   (0,T_max).












For p>1, testing the first equation of problem (<ref>) by u^p-1 and using its boundary conditions  provide

    1/pd/dt∫_Ω u^p   =∫_Ω u^p-1u_t =-(p-1)∫_Ω u^p-2|∇ u|^2  
          +(p-1) ∫_Ω  u^p-1χ(v)∇ u ·∇ v   for all   t∈ (0,T_max).

As to the integral involving u^p-1χ(v)∇ u ·∇ v, the Young inequality and the uniform bound for χ derived in the third relation of (<ref>)  yield for all t∈(0,T_max)

    ∫_Ω|u^p-1χ(v)∇ u ·∇ v|   ≤1/2∫_Ω u^p-2χ^2(v)|∇ u|^2+ 1/2∫_Ω u^p|∇ v|^2 
       ≤1/2‖χ‖_^2∫_Ω u^p-2|∇ u|^2+1/2∫_Ω u^p|∇ v|^2,

so that equality (<ref>)  writes exactly as claimed.


 
For  n=2 and under the remaining assumptions of Lemma <ref>, let (u,v) be the local-in-time classical solution  of problem (<ref>). Then for any p∈(1,2) and all  t∈ (0,T_max) holds 

    d/d t∫_Ω u^p ≤ -p(p-1)/2∫_Ω u^p-2|∇ u|^2+(ε+c_8δ)∫_Ω u^p+1+c_10,

where c_8 and c_10 are computable positive constants, ε an arbitrary  positive number and δ another arbitrary nonnegative real such that 

    δ=0    if  β=0,
    δ> 0    if  β>0.
 

In order to estimate the term ∫_Ω u^p |∇ v|^2 appearing in the above Lemma <ref>, we observe that an application of the Hlder inequality infers for any p>1 

    ∫_Ω u^p |∇ v|^2 ≤(∫_Ω u^p+1)^p/p+1(∫_Ω |∇ v|^2p+2)^1/p+1=‖ u ‖^p_L^p+1(Ω)‖∇ v ‖^2_L^2p+2(Ω).
 
Now, in view of the fact that the v-component solves -Δ v+ v=g(u) in Ω× (0,T_max),  for p∈(1,2) we can invoke Lemma <ref> with f=v which, together with (<ref>),  implies

    (∫_Ω |∇ v|^2p+2)^1/p+1   =‖∇ v ‖^2_L^2p+2(Ω)≤Ĉ^2 ‖ -Δ v + v‖_L^p+1(Ω)‖∇ v‖_L^2(Ω)
       ≤Ĉ^2√(M)‖ g(u)‖_L^p+1(Ω)  for all   t∈ (0,T_max).

Moreover, assumption (<ref>) gives

    ‖ g(u)‖_L^p+1(Ω)≤λ_2 (∫_Ω (1+u)^β(p+1))^1/p+1  on   (0,T_max),

so that, by means of Young's inequality and the introduction of a positive real number ε, we deduce from (<ref>), (<ref>),  (<ref>) and (<ref>) that on  (0,T_max) 

    ∫_Ω u^p |∇ v|^2    ≤(∫_Ω u^p+1)^p/p+1((Ĉ^2√(M)λ_2)^p+1∫_Ω (1+u)^β(p+1))^1/p+1
       ≤ε/c_8∫_Ω u^p+1 +c_7/p+1(ε(p+1)/c_8p)^-p∫_Ω (1+u)^β(p+1)
       ≤ε/c_8∫_Ω u^p+1 +c_92^β(p+1)|Ω|+c_92^β(p+1)∫_Ω u^β(p+1),   β >0,
    ε/c_8∫_Ω u^p+1 +c_9|Ω|   β =0,
 
with c_7=(Ĉ^2√(M)λ_2)^p+1, c_8=p(p-1)/2‖χ‖_^2 and c_9=c_7/p+1(ε(p+1)/c_8p)^-p. In particular, for β >0, similar fashions allow us to show that for some δ>0, on (0,T_max) also holds 

    ∫_Ω   u^β(p+1)≤δ/c_9 2^β(p+1)∫_Ω u^p+1 + (1-β)(δ/c_9 β 2^β(p+1))^-β/1-β|Ω|.

Finally, taking into account (<ref>)  of Lemma <ref>, with relations  (<ref>) and (<ref>)  we conclude the proof setting 

    c_10=
    
    c_8c_9 |Ω|    for  β=0,
    
    c_8c_9 2^β(p+1)|Ω|(1+(1-β)(δ/c_9 β 2^β(p+1))^-β/1-β)   for  β>0.



As announced, the succeeding lemma is dedicated to establish the desired uniform-in-time bound for ‖ u‖_L^p(Ω), with some proper p>1. 

For  n=2 and under the remaining assumptions of Lemma <ref>, let (u,v) be the local-in-time classical solution  of problem (<ref>). Then for any p∈(1,2) there exists a positive constant C_p such that 

    ∫_Ω u^p ≤ C_p   for all   t ∈ (0,T_max).


We rely on the Gagliardo-Nirenberg inequality to estimate the term ∫_Ω u^p+1 appearing in (<ref>). Precisely, expression (<ref>) with f=u^p/2, 𝔧=0, 𝔪=1, 𝔯=n=2, 𝔭=2(p+1)/p and 𝔮=𝔰=2/p, in conjunction with (<ref>), provides  for  all t∈ (0,T_max) the relation

    ∫_Ω u^p+1   =|| u^p/2||_L^2(p+1)/p(Ω)^2(p+1)/p≤ c_11||∇ u^p/2||_L^2(Ω)^2(p+1)/pθ_1|| u^p/2||_L^2/p(Ω)^2(p+1)/p(1-θ_1) +c_11|| u^p/2||^2(p+1)/p_L^2/p(Ω),

being c_11= (2C_GN)^2(p+1)/p and 0<θ_1=p/p+1<1. In particular, this gained inequality and (<ref>) give

    (ε+c_8δ)∫_Ω u^p+1   ≤(ε+c_8δ) c_11 m ∫_Ω |∇ u^p/2|^2 
          + (ε+c_8δ)c_11m^p+1  for all   t ∈ (0,T_max),

so that (<ref>) is transformed in 
 
    d/d t∫_Ω u^p ≤(c_11m(ε+c_8δ)-2(p-1)/p)∫_Ω |∇ u^p/2|^2 + c_12,

with  c_12=(ε+c_8δ)c_11m^p+1+c_10. 

Successively we again use Lemma <ref>  and in the Gagliardo-Nirenberg inequality  (<ref>) we take f=u^p/2, 𝔮=𝔰=2/p, 𝔭=n=𝔯=2, 𝔧=0 and 𝔪=1; we infer through (<ref>) that on (0,T_max) one has

    ∫_Ω u^p   =|| u^p/2||_L^2(Ω)^2≤ c_13||∇ u^p/2||_L^2(Ω)^2θ_2|| u^p/2||_L^2/p(Ω)^2(1-θ_2) +c_13|| u^p/2||^2_L^2/p(Ω),

where 0<θ_2=p-1/p<1 and c_13= (2C_GN)^2. Considering once again bound (<ref>) and  introducing c_14=c_13max{m,m^p}, by making again use of (<ref>)  inequality (<ref>) can also be rewritten as 

    -∫_Ω|∇ u^p/2|^2≤ 1-(2c_14)^p/1-p(∫_Ω u^p)^p/p-1  for all   t ∈ (0,T_max).

We, then, choose 

    if  β=0,    (and so  δ=0)   ε=p-1/c_11mp>0,
    if  β >0,    ε=p-1/2c_11mp>0   and  δ=p-1/2c_8c_11mp>0,
 
and plug such values and  estimate (<ref>) into (<ref>); in this way we arrive at this initial problem

    Φ'(t)≤ c_15-c_16Φ(t)^p/p -1   t ∈ (0,T_max),
    Φ(0)=∫_Ω u_0^p,

where Φ(t):=∫_Ω u^p, c_15=c_12+p-1/p and c_16=p-1/p(2c_14)^p/1-p. Ultimately, we conclude the demonstration by an application of a comparison principle implying 

    Φ(t)≤max{Φ(0),(c_15/c_16)^p -1/p}=:C_p  for all   t∈(0,T_max).



As a consequence of all of the above preparations, we finally can prove our claimed statement:

 

Proof of Theorem <ref>. For n=2  and any nonnegative 0≢u_0∈ C^0(Ω̅), Lemma <ref> provides a unique local-in-time classical solution (u,v) to problem (<ref>). Thereafter, by virtue of Lemma <ref>, for all p∈(1,2) relation (<ref>) is achieved, then u∈ L^∞((0,T_max); L^p(Ω)) and we can conclude through Lemma <ref>.



§ NUMERICAL SIMULATIONS

In this section we numerically test the presented theoretical results by simulating system (<ref>) in two dimensions; for simplicity and with no possibility of confusion, the spatial variable x=(x_1,x_2) is indicated with (x,y). Further, we investigate whether the solutions are: stationary (time independent) and homogeneous (spatially uniform density); stationary and heterogeneous (spatially non-uniform density), but bounded; or suffer from chemotactic blow-up in finite time.

Specifically, we use finite element methods to simulate system (<ref>) with a variety of functions defining χ(v) and g(u). The solution algorithm is based on an adaptive, implicit Runge-Kutta finite element method (see <cit.>). The space Ω is defined to be the interior of the [0,0.1]×[0,0.1] square, making ∂Ω the boundary of the square. To some extent the size of the space is arbitrary since we can use a larger domain, rescale the system with respect to the size and, thus, reproduce the same dynamics. However, because of anticipating large and thin spike structures in the cells' density we choose to simulate a small domain, which allows us to resolve such small heterogeneous solutions more accurately, without increasing the overall simulation mesh resolution.

From system (<ref>) we impose to the unknowns u and v to satisfy Neumann conditions on the boundary. The initial condition for u was chosen to be uniform constant, u̅, plus noise, that is 

    u(x,y,0)=u_0=|u̅+ση(x,y)|,

where η:Ω→ [-1/2,1/2] is a continuous uniformly random variable and σ is a positive constant that allows the stochastic perturbation to be scaled with respect to the mean value u̅. On the other hand, even though v(x,y,0)=v_0 is not required in system (<ref>) (exactly because the equation for v is elliptic), in order to run the Runge-Kutta iterative method an initial condition for v has to be also assigned: in particular, we take v_0 as the solution of -Δ v_0+v_0=g(u_0) under Neumann boundary conditions. 

Since we are looking for differences between solutions that are bounded, versus those that suffer from blow-up, if a solution is found to increase indefinitely, then the same simulation was repeated with a finer discretisation to ensure that this outcome is the true numerical solution, rather than a numerical artefact. Specifically, whenever a solution was observed to be grow without bound, the grid was refined to have ten times as many elements as previously simulated, to ensure the outcome. As we will see later, system (<ref>) can support heterogeneous spike solutions in the densities of u. Critically, as parameters of interest are altered, the spikes densities become larger, whilst their support becomes smaller. As the spikes tend to the form of a δ-function, numerical solution convergence requires finer discretisations in order to resolve the spike's shape. However, limitations in computer memory and processing power limit this procedure of refinement. Such cases will be clearly highlighted.

Figure <ref> illustrates one of the results discussed in <ref>. To be precise, when χ(v)=χ>0, constant, and g(u)=u then solutions blow up if the initial mass ∫_Ω u_0 and χ are such that χ∫_Ω u_0 is  sufficiently large. Critically, when χ=10^3 and u̅=10 (right image of Figure <ref>), then u grows to over 10^10 in less than 10^-5 time units. The three simulations from the left of Figure <ref> show what happens when either (or both) the initial condition, u̅, or the sensitivity coefficient, χ, is reduced by a factor of 10, namely the solution rapidly converges to the homogeneous steady state,

    (u_s,v_s)=( 1/|Ω|∫_Ω u(x,y,0),g(u_s)).



In Figure <ref> we simulate system (<ref>) again, but this time with g(u)=(1+u)^1/2. Initially, we simulated this new system using the same parameters, u̅, σ and χ, as those specified in Figure <ref> (data not shown) and quickly realised that the simulations no longer succumbed to blow-up, rather each variable (u,v) converged to a finite stationary distribution. Figure <ref> demonstrates we are able to increase u̅ and χ without fear of blow-up. However, increasing these two parameters leads to the uniform steady stable being driven unstable, thus, the system evolves to a stable, bounded heterogeneous density.



 As specified above, increasing the initial condition, u̅, and/or the value of the sensitivity, χ, causes the spike solution to become sharper, giving problems with numerical convergence as the values are increased. Thus, although the simulated solutions begin to grow as u̅ and χ are enlarged, this is probably an issue of the numerical resolution, rather than an analytical singularity.


Next, in Figure <ref>, we show that the solutions are bounded even when χ(v) is chosen to be non-trivial. Further, for the same parameter values, changing between χ(v)=χ/v and χ(v)=χlog(v) also causes the simulations to alter between a homogeneous solution and a heterogeneous solution (compare figures <ref> and <ref>).


Initially, and in accordance to the theoretical results, our simulations corroborated that if g(u)=u the solution could become unbounded, whereas the global boundedness is guaranteed when g(u)=(1+u)^1/2. Thus, a natural question consists in analysing the possible existence of a  critical exponent, β in g(u)=(1+u)^β, under which system (<ref>) admits stationary and homogeneous, stationary and heterogeneous, but bounded, or (even) unbounded solutions: this is addressed in the examples included in Figure <ref>. Suppose χ(v)=χ>0 and parameters for the data u̅ and χ are chosen such that only the simulation tends to a homogeneous steady state: then a simple parameter sweep suggests that β≈ 1 is the bifurcation point (data not shown). However, if the parameters are chosen such that spike solutions exist, then the maximum values of the spikes grow as β increases, as shown in Figure <ref>. Here, we see that if β≤ 2/5 only the homogeneous solution is found, whereas, for β≥ 2/5 a heterogeneous solution appears. Notably, the spikes become numerically unstable for β≈ 0.76. Again, this may be due to the coarsity of the underlying mesh trying to resolve the large, but thin spikes (note that the vertical axes in Figure <ref> are logarithmic). These simulations can be compared with those seen in Figure <ref>, where χ(v)=χ/v. Specifically, although we see a discontinuous jump between the homogeneous and heterogeneous solution branches at around β=0.64 (suggesting a subcritical bifurcation), we note that the maximum value of u does not appear to increase much beyond 10^7. Thus, despite the system having a spike solution, the non-trivial sensitivity produces a sort of “controlling effect” on the growth of the cells' distribution, u, even for β>1.


Over all our simulations match the analytical insights, in the ranges that are numerically feasible. Specifically, we have shown that the general properties of solutions to system (<ref>) heavily depend on the form of g(u) and χ(v). Our numerical solutions also illustrate the need for the presented theoretical results, which confirm the existence and boundedness of solutions. Namely, since sharply peaked solutions are very difficult to numerically resolve, due to disparate scales within the system, the theoretical analysis is required to inform us as to the accuracy of our numerical schemes.





  §.§.§ Acknowledgments


GV is member of the Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Na­zio­na­le di Alta Matematica (INdAM) and is partially supported by the research project Integro-differential Equations and Non-Local Problems, funded by Fondazione di Sardegna (2017). 









10

Ascher-1997-IER
U. M. Ascher, S. J. Ruuth, and R. J. Spiteri.
Implicit-explicit Runge-Kutta methods for time-dependent partial
  differential equations.
Appl. Numer. Math., 25(2-3):151–167, 1997.

BellomoEtAl
N. Bellomo, A. Bellouquid, Y. Tao, and M. Winkler.
Toward a mathematical theory of Keller–Segel models of pattern
  formation in biological tissues.
Math. Models Methods Appl. Sci., 25(09):1663–1763, 2015.

BilerLocalAndGlobal
P. Biler.
Local and global solvability of some parabolic systems modelling
  chemotaxis.
Adv. Math. Sci. Appl., 8(2):715–743, 1998.

BrezisBook
H. Brezis.
Functional analysis, Sobolev spaces and partial differential
  equations.
Universitext. Springer, New York, 2011.

FujieWinklerTomomiParaEllip
K. Fujie, M. Winkler, and T. Yokota.
Blow-up prevention by logistic sources in a parabolic-elliptic
  Keller-Segel system with singular sensitivity.
Nonlinear Anal., 109:56–71, 2014.

FujieWinklerYokotaSignalDependent
K. Fujie, M. Winkler, and T. Yokota.
Boundedness of solutions to parabolic-elliptic Keller-Segel
  systems with signal-dependent sensitivity.
Math. Methods Appl. Sci., 38(6):1212–1224, 2015.

Hillen2009UGP
T. Hillen and K. J. Painter.
A user's guide to PDE models for chemotaxis.
J. Math. Biol., 58(1):183–217, 2009.

HorstWang
D. Horstmann and G. Wang.
Blow-up in a chemotaxis model without symmetry assumptions.
Eur. J. Appl. Math., 12(2):159–177, 2001.

JaLu
W. Jäger and S. Luckhaus.
On explosions of solutions to a system of partial differential
  equations modelling chemotaxis.
Trans. Amer. Math. Soc., 329(2):819–824, 1992.

K-S-1970
E. F. Keller and L. A. Segel.
Initiation of slime mold aggregation viewed as an instability.
J. Theor. Biol., 26(3):399–415, 1970.

Keller-1971-TBC
E. F. Keller and L. A. Segel.
Traveling bands of chemotactic bacteria: A theoretical analysis.
J. Theor. Biol., 30(2):235, 1971.

LSUBookInequality
O. A. Ladyženskaja, V. A. Solonnikov, and N. N. Ural'ceva.
Linear and Quasi-Linear Equations of Parabolic Type.
In Translations of Mathematical Monographs, volume 23.
  American Mathematical Society, 1988.

Lankeit
J. Lankeit.
Eventual smoothness and asymptotics in a three-dimensional chemotaxis
  system with logistic source.
J. Differential Equations, 258(4):1158–1191, 2015.

LankeitANewApproach
J. Lankeit.
A new approach toward boundedness in a two-dimensional parabolic
  chemotaxis system with singular sensitivity.
Math. Methods Appl. Sci., 39(3):394–404, 2016.

Nagai
T. Nagai.
Blowup of nonradial solutions to parabolic-elliptic systems modeling
  chemotaxis intwo-dimensional domains.
J. Inequal. Appl., 6(1):37–55, 2001.

NagaiSenba
T. Nagai and T. Senba.
Global existence and blow-up of radial solutions to a
  parabolic-elliptic system of chemotaxis.
Adv. Math. Sci. Appl., 8(1):145–156, 1998.

NirenbergOnEllipticPDE
L. Nirenberg.
On elliptic partial differential equations.
Ann. Scuola Norm. Sup. Pisa (3), 13:115–162, 1959.

OsYagUnidim
K. Osaki and A. Yagi.
Finite dimensional attractor for one-dimensional Keller-Segel
  equations.
Funkcial. Ekvacioj., 44(3):441–470, 2001.

SleemanLevine
B. D. Sleeman and H. A. Levine.
Partial differential equations of chemotaxis and angiogenesis.
Math. Methods Appl. Sci., 24(6):405–426, 2001.

OthmerStevens
A. Stevens and H. G. Othmer.
Aggregation, blowup, and collapse: The abc's of taxis in reinforced
  random walks.
SIAM J. Appl. Math., 57(4):1044–1081, 1997.

TelloWinkParEl
J. I. Tello and M. Winkler.
A chemotaxis system with logistic source.
Commun. Part. Diff. Eq., 32(6):849–877, 2007.

ViglialoroPreprintNonLinear
G. Viglialoro.
Global in time and bounded solutions to a parabolic-elliptic
  chemotaxis system with nonlinear diffusion and signal-dependent sensitivity.
Preprint.

ViglialoroVeryWeak
G. Viglialoro.
Very weak global solutions to a parabolic–parabolic
  chemotaxis-system with logistic source.
J. Math. Anal. Appl., 439(1):197–212, 2016.

ViglialoroBoundnessVeryWeak
G. Viglialoro.
Boundedness properties of very weak solutions to a fully parabolic
  chemotaxis-system with logistic source.
Nonlinear Anal. Real World Appl., 34:520–535, 2017.

ViglialoroWolleyDCDS
G. Viglialoro and T. Woolley.
Eventual smoothness and asymptotic behaviour of solutions to a
  chemotaxis system perturbed by a logistic growth.
Discrete Continuous Dyn. Syst. Ser. B., 22(5):1–23, 2017.

WinklAggre
M. Winkler.
Aggregation vs. global diffusive behavior in the higher-dimensional
  Keller-Segel model.
J. Differential Equations, 248(12):2889–2905, 2010.

WinDespiteLogistic
M. Winkler.
Blow-up in a higher-dimensional chemotaxis system despite logistic
  growth restriction.
J. Math. Anal. Appl., 384(2):261–272, 2011.

WinklerFiniteTeimeBlowUpHigher
M. Winkler.
Finite-time blow-up in the higher-dimensional parabolic-parabolic
  Keller-Segel system.
J. Math. Pures Appl., 100(5):748–767, 2013.

WinklerZAMPLogisticBlowUP
M. Winkler.
Finite-time blow-up in low-dimensional Keller–Segel systems with
  logistic-type superlinear degradation.
Z. Angew. Math. Phys., 69(2):69:40, 2018.



