key: cord-0641346-jueoz0ok
authors: Calvez, Vincent; Hoffmann, Franca
title: Nonlinear stability of chemotactic clustering with discontinuous advection
date: 2020-09-23
journal: nan
DOI: nan
sha: 3afeb0209220f4952e6c158b864b43634a8a9510
doc_id: 641346
cord_uid: jueoz0ok

We perform the nonlinear stability analysis of a chemotaxis model of bacterial self-organization, assuming that bacteria respond sharply to chemical signals. The resulting discontinuous advection speed represents the key challenge for the stability analysis. We follow a perturbative approach, where the shape of the cellular profile is clearly separated from its global motion, allowing us to circumvent the discontinuity issue. Further, the homogeneity of the problem leads to two conservation laws, which express themselves in differently weighted functional spaces. This discrepancy between the weights represents another key methodological challenge. We derive an improved Poincar'e inequality that allows to transfer the information encoded in the conservation laws to the appropriately weighted spaces. As a result, we obtain exponential relaxation to equilibrium with an explicit rate. A numerical investigation illustrates our results.

This work is devoted to the stability analysis of stationary clusters of bacterial cells under the effect of chemotaxis, as reported in the biophysical literature, see e.g. [16] . Certain type of bacteria such as Escherichia coli are sensible to chemical gradients and can communicate with each other by means of chemoattractants. In the absence of nutrients, the long-time asymptotics of the cell population at the macroscopic level depend on an intricate interplay between diffusive forces and chemotactic aggregation. A classical model for bacterial motion is the Patlak-Keller-Segel model and variants where the chemotactic fluxes are derived analytically from a mesoscopic description of the dynamics at the individual level and possibly internal pathways [12, 9, 8] . Here, we focus on a minimal model with discontinuous advection speed and, for the sake of simplicity, in one dimension of space (see [20] and discussion below). The cell density is denoted by ρ(t, x), and the chemoattractant concentration by S(t, x). Cluster formation of cells in quasi-equilibrium can be modeled by a version of the classical Patlak-Keller-Segel model with linear diffusion and a discontinuous drift term,

where χ > 0 is the chemoattractant sensitivity, and α ≥ 0 is the natural decay of the chemical. The system (1.1) is equipped with an initial condition ρ(0, x) = ρ 0 (x), whose regularity and decay at infinity will be discussed below.

This model differs from the standard Patlak-Keller-Segel model for which the advection speed is essentially linear with respect to the chemical gradient, χ∂ x S , see [2, 1, 18] for recent reviews. The discontinuous nonlinearity χsign(∂ x S) is the signature of strong amplification of signal variations at the individual level. This is the extreme point of a family of chemotaxis models with non-linear dependency upon the chemical gradient that can be applied to shallow and steep gradients, see [19] for the original derivation and [10] for a biological validation. We also refer to [13] for a review on mathematical modeling of chemotaxis (see, in particular Model (M7)), and [23] for a review dedicated to bacterial collective motion. The sign function was also used to reproduce traveling bands of bacteria with good agreement in [20] , which is the early motivation for the present work.

This strong discrepancy -χsign(∂ x S) versus χ∂ x S -requires a specific approach to handle the stability analysis of stationary states of (1.1). There exists a family of stationary states to model (1.1) given by (1.2) ρ ∞ (x) = M χ 2 e −χ|x−a| , (M, a) ∈ R * + × R .

The first paramater M is linked to the conservation property of (1.1a), and is fully determined by the mass of the initial condition, M = ρ 0 (x) dx since the equation for the cell density is conservative. The second parameter, a is linked to the invariance by translation, and is also determined by the initial condition, though in a non-trivial way. It is interesting to notice that the natural length scale L for the typical size of a cell cluster associated to the state (1.2) is L = χ −1 , independent of the total mass of cells M . This is in agreement with the observations reported in [16] , where the typical size of the cell cluster varies little with the number of bacteria. This is in opposition with the standard Patlak-Keller-Segel model, for which the length scale is L = (χM ) −1 due to the homogeneity of the problem (twice the number of cells results in twice the quantity of chemoattractant, and this in turn increases twofold the advection speed).

In the present work, we assume M = 2/χ without loss of generality (in order to cancel the prefactor in (1.2)). Further, we assume α > 0 as the case α = 0 can be treated in a more straighforward way with different methods, see Remark (5) following Theorem 1.3, or Appendix B for more details.

We stress that we are not concerned here with existence, uniqueness and regularity issues for solutions to (1.1). We assume that solutions are sufficiently regular for our calculations to hold, that is, the cell density is continuous at all times t ≥ 0, ρ(t, ·) ∈ C 0 (R), while its derivative ∂ x ρ develops jump discontinuities at points where ∂ x S changes sign (see below for more details).

Our goal is to prove local stability for the nonlinear problem (1.1) around the family of stationary states (1.2). The difficulty here arises from the discontinuity of the advection speed χsign(∂ x S). On the one hand, the choice of the sign function rules out a direct linearization of the non-linear term. On the other hand, (1.1a) is a piecewise linear advection-diffusion equation, up to the knowledge of those points where ∂ x S changes sign. We base our strategy on the latter observation. More precisely, we will crucially use the following preliminary result:

The dynamics of x(t) inherit from the condition ∂ x S(t, x(t)) = 0, that is

x(t))ẋ(t) = 0 . The formula above can be combined with the following representation of S = K α * ρ for α > 0, in order to derive a suitable expression forẋ. Here, K α is the fundamental solution of (1.1b):

This formulation enables to derive the following dynamics ofẋ in the moving frame, as a function of the half derivatives of the cell density at the interface:

.

(see Lemma 2.2 for details). It is then natural to reformulate (1.1) in the moving frame y = x − x(t), writing forρ(t, y) = ρ(t, x) andS(t, y) = S(t, x):

In doing so, we focus on the stability analysis of the shape of cell densityρ, and separate this question from the dynamics of the maximum point x of S.

The stationary state for (1.6) is then simply given by

We are now in a position to state our main result. Let H 1 χ be the weighted space of relative energy equipped with the following norm

. Theorem 1.3. Let M = 2/χ and α > 0. The family of stationary states {e −χ|x−a| } a∈R to model (1.1) is locally nonlinearly stable in the following sense: There exists ε 0 > 0 depending on χ and α such that, for all initial data satisfying

, the constant C depends onρ 0 , ε 0 , χ and α, and the rate γ can be chosen arbitrarily below the following upper bound (at the expense of increasing the prefactor C for larger choices of γ):

We make the following observations:

(1) The limit x ∞ exists but has no explicit value, up to our knowledge. Its precise dependence on the initial configuration and parameters of the model is not known. Our analysis is not meant to derive a rate of convergence of x(t) → x ∞ .

(2) The smallness condition ρ 0 −ρ ∞ H 1 χ ≤ ε 0 does not control the initial value ofẋ(0) which involves pointwise values ∂ yρ0 (0 ± ). Therefore, it is possible that x(t) has large variations, meaning that x(0) and x ∞ are far apart. In fact, the convergence of x(t) is obtained by means of the dissipative structure of the parabolic equation (1.6a).

(3) The upper bound of the convergence rate γ 0 in (1.10) is between χ 2 8 (α → ∞) and χ 2 4 (α → 0). However, our estimates on the convergence rate are not optimal, as resulting from successive inequalities. To obtain an optimal bound via functional inequalities, one could combine the improved Poincaré inequality (3.1) with the interpolation inequality (3.5) into a single inequality and seek optimizers. We do not follow this approach here as the two separate inequalities contain meaningful structure. (4) The limit α → ∞ is quite singular, as it is formally equivalent to the following problem (after scaling of S → S/α and the identification in the vanishing viscosity limit S = ρ):

However, we have no insight about the above problem. Regarding (1.9), it should be noted that γ can be chosen independently of α, but the prefactor C becomes singular as α → ∞ in our methodology, due to the control of non-linear contributions. Method. We emphasize one key methodological contribution. The problem (1.6) is equipped with two conservation laws, corresponding to the homogeneity of the problem (1.1a), and its invariance by translation (see Section 2.3 for more details):

Notably, there is a discrepancy between the weights, here 1 and e − √ α|y| , turning into e −χ|y| and e −(χ+ √ α)|y| in relative energy, see (1.8). As a result, it is not obvious which is the correct choice of functional space to work in. The less restrictive option e −(χ+ √ α)|y| might be considered.

However, this results in a deterioration of the lower bound on the convergence rate in (1.10). It is not even clear that all values of α can be encompassed. The alternative choice e −χ|y| turns out to be much more satisfactory. At the core of our method is an improved version of the standard Poincaré inequality with exponential weight, that allows to transfer the information given by the second conservation law to the appropriate weighted space, see Proposition 3.1.

Motivation and Perspectives. Our initial motivation comes from the mathematical modeling and analysis of concentration waves of chemotactic bacteria. We refer to [23] for a comprehensive review of modeling of bacteria colonies. The present work is rooted in [20] where the following model was proposed for the propagation of density bands under the conjunct effect of two chemotactic signals (a communication signal S secreted by the cells, and a nutrient signal N depleted by the cells).

x N − γρN . The connection to (1.1) is the following. Firstly, in the absence of a food source, the nutrient signal N is omitted, resulting in a zero speed traveling band, i.e. a stationary cluster, as reported in [16] . Secondly, the communication signal S is assumed to be in quasi-stationary equilibrium, that is (1.1b). We also consider D ρ = D S = 1 without loss of generality (after appropriate non-dimensionalization), and we denote χ S = χ the chemosensitivity.

The agreement with experimental data was shown to be very satisfactory [20] . Moreover, the model (1.12) was derived from a kinetic transport model suitable to describe the individual run-and-tumble motion of bacteria. The kinetic model agreed very well with another set of experimental data [21] .

It is one of the key perspectives of this work to prove stability of the stationary cluster solution of the kinetic-transport equation (whose existence is a particular case of [5] ). Note that exponential relaxation to equilibrium was established in [7] for the linear problem, where the communication signal S is prescribed, using hypocoercive techniques. This was extended in [15] to any dimension of space.

A second key perspective consists in proving non-linear stability of the traveling concentration waves, solution of (1.12), see [20] and [5] for existence of such waves, resp. for the parabolic problem and for the kinetic transport problem.

From now on, we work in the frame centered at x(t): y = x − x(t). Also, we drop the symbol ∼ for notational convenience.

2.1. Uniqueness of chemoattractant peak. In this section, we will show that the unique maximum property for the chemoattractant concentration S is preserved by the flow, and provide a proof for the corresponding main result Proposition 1.1.

To show this, we "decouple" the dynamics to the right of the maximum from the dynamics to the left. This can be done since in the moving frame, ∂ y S(t, y) solves the following equation on the half-space:

where the source term ∂ y ρ(t, y) is given via the solution to the following PDE:

, y > 0 . We denote by A the fundamental solution of (2.1), such that

Namely, we have (although the exact expression can be omitted)

The following reformulation is of key importance to control the sign of ∂ y S on R + . Lemma 2.1. Let S be a solution to (2.1). Then ∂ y S satisfies the following boundary value problem:

By definition of A, the first contribution

On the other hand, we have

y S satisfies the following boundary value problem: 

Similarly, we find that the second contribution S 2 :

Thus,

The value of c(t) depends on the dynamics ofẋ, see (1.6a). We develop its expression explicitly, focusing on the right hand side by using the continuity of the flux at the interface:

where the second line follows from (1.6b) and (2.1). Substituting this expression into the one above, we obtain the result.

Proof of Proposition 1.1. For ε > 0, define

In particular, we have ∂ y S ε (t, 0) = −ε < 0 for all t > 0, by definition, and ∂ y S ε (0, y) ≤ −ε < 0 for all y > 0, by assumption. For any fixed ε > 0, we aim to prove that ∂ y S ε (t, y) < 0 for all t, y > 0. Let us suppose by contradiction that there exists a smallest time t 0 > 0 and closest location y 0 > 0 such that ∂ y S ε (t 0 , y 0 ) = 0. It is necessarily an interior maximum value of ∂ y S ε with respect to y, so that the following conditions are verified:

By letting ε → 0, we conclude that ∂ y S(t, y) ≤ 0 for all t, y > 0. Then, we can use the PDE (2.4) again to show that the latter inequality is strict: in fact, it is a drift-diffusion equation on ∂ y S(t, y) with negative source term and non-positive initial data. We just proved that the solution remains non-positive. By application of the strong maximum principle, it cannot have an interior maximum point, so ∂ y S(t, y) is negative for t, y > 0.

A symmetric argument to the one presented above shows that ∂ y S(t, y) > 0 for t > 0 and y < 0, which concludes the proof of Proposition 1.1.

Dynamics of the chemoattractant peak. The dynamics of the point x(t) are given by the following lemma.

This is actually equivalent to the second conservation law (1.11b) after integration by parts. Differentiating with respect to time and substituting (1.6a), we obtain

where the flux

is continuous at y = 0 (the jump in the first derivative compensates exactly the change of sign). Simplifying the last term by integration by parts on the first component of j, we get

The two middle contributions vanish because of the conservation law (2.6). In addition, by continuity we have

which is equivalent to (2.5) using (1.6b).

Reformulation of the problem and conservation laws. In fact, using the continuity of the flux j defined in (2.7) at y = 0, we can reformulate the dynamics of x(t) as follows,

Substituting this expression into (1.6), we find that the dynamics of the cell density ρ are governed by

This formulation strongly suggest to work with relative densities:

Then, u satisfies the system ∂ t u(t, y) = e χ|y| ∂ y e −χ|y| ∂ y u(t, y) + ∂ y u(t, 0) ∂ 2 y S(t, 0)

Note that ∂ y u is continuous at y = 0, since the C 1 discontinuity of ρ has been exactly compensated in (2.10) .

From now on, we shall work with (2.11). This system admits the stationary state u ∞ (y) ≡ 1, and two conservation laws:

• conservation of mass:

(recall that the mass is fixed to M = 2/χ), and • centering frame: Setting λ := χ + √ α and using the notation · r for the weighted average,

the two conservation laws can be written equivalently as

In this section, we derive H 1 energy estimates to measure dissipation in (2.11). Interestingly, we find that dissipation always occurs, as the non-local interaction contribution is overwhelmed by heat dissipation. Our analysis relies on the following Poincaré-type inequality which is designed to handle the discrepancy in the exponential rates of the pair of conservation laws (resp. χ and λ = χ + √ α).

Proposition 3.1 (Poincaré inequality with unusual normalization). Assume λ ≥ χ, then for any w ∈ L 1 λ such that w ′ ∈ L 2 χ ,

Moreover, the constant 4 χ 2 is optimal.

Inequality (3.1) is the standard Poincaré inequality with exponential weight when λ = χ. In that case, the optimal constant 4 χ 2 can be deduced from spectral analysis of the linear operator −w ′′ + χ(signy)w ′ in the Hilbert space L 2 χ = L 2 (e −χ|y| ) which admits a spectral gap (0, χ 2 4 ), with an isolated value 0, and essential spectrum beyond χ 2 4 . Inequality (3.1) is an improvement of the standard Poincaré inequality (λ = χ), simply because

It may appear surprising at first glance that the optimal constant does not depend on the exponent λ ≥ χ. This is a consequence of the fact that the optimal constant is never reached in (3.1) because the exponential weight has critical decay, and the putative optimal functions are not in L 2 χ . This yields some room for improving (3.2) as in (3.1). We are not aware of the occurrence of (3.1) in the literature. There are some examples of weighted Poincaré inequalities with non standard averaging, see for instance [4, Appendix A.2] . We can also report the extremal case λ = +∞ which coincides with the following Hardy inequality:

The link between the latter and the classical Poincaré inequality with exponential weights is pointed out in [14, Corollary 1.3 ], see also [3] . Inequality (3.1) can be viewed as a continuous deformation connecting (3.2) and (3.3) with a constant optimal constant. Remark 3.2. We recall below the equivalence between (3.3) and the classical Hardy inequality. We restrict to the half-line {y > 0} and to χ = 1 for the sake of clarity. We make the change of variable x = ke y on both sides:

Letting k → 0 we recover the classical Hardy inequality on the half-line.

It is noticeable that the optimal constant in (3.1) does not depend on λ, provided it is strictly greater than χ. For the proof of Proposition 3.1 we present two arguments: first some insights based on spectral analysis (not complete enough to provide a rigorous proof), followed by a complete (technical) proof of Proposition 3.1 based on the reformulation of (3.1) as a quadratic form. Note that we can always restrict our arguments to χ = 1 without loss of generality by rescaling y to χy.

We begin with an analysis of the spectrum for the operator associated to inequality (3.1). Studying carefully the critical functions w of (3.1) (such that (i) w λ = 0 and (ii) 1 2 w 2 e −|y| dy = 1), we find that they must satisfy the following problem:

Here, the scalars µ and ν are Lagrange multipliers corresponding to the two constraints (i) and (ii). Actually, this problem is equivalent to the next one after the change of unknown w(y) = W (y)e |y|/2 :

This problem can be explicitly solved after tedious computations by means of the fundamental solution (details not shown). On the one hand, the constant ν must be adjusted to ensure w L 2 1 = 1. On the other hand, the condition w λ = 0 cannot be satisfied if µ < 1 4 . When µ = 1 4 , we have the following expression (general solution of the homogeneous problem augmented by a particular solution):

where the constant C can be adjusted to ensure w λ = 0 when λ > 1. However, one immediately notices that w / ∈ L 2 1 , simply because 1 4 belongs to the essential spectrum of −w ′′ + (signy)w ′ . Nonetheless, considering the quotient (3.7) R −R |w ′ (y)| 2 e −|y| dy R −R |w(y)| 2 e −|y| dy for arbitrary large values of R, we see that the additional contribution of the particular solution O e (1−λ)|y| becomes negligible, so that the quotient converges to the constant 1 4 as R → +∞, independently of λ. This is the reason why the optimal constant in (3.1) does not depend on λ in the case λ > 1.

Next, we provide a complete proof of Proposition 3.1 by direct computations. It relies on the following technical lemma, whose proof is postponed to Appendix A.

where M λ denotes the cumulative density function,

Then Ω λ,χ is non-negative and symmetric,

and we can rewrite the left-hand side of the Poincaré inequality (3.1) as

Proof of Proposition 3.1. We assume again χ = 1, denote Ω λ := Ω λ,1 , and introduce the new function W (x) = w ′ (x)e −|x|/2 . It then follows from Lemma 3.3 that the Poincaré inequality (3.1) is equivalent to

The quadratic form on the left hand side can be reformulated as follows:

] Ω λ (x 1 , x 2 )e |x 1 |/2 e |x 2 |/2 dx 1 dx 2 + |W (x)| 2 Ω λ (x, y)e |x|/2 e |y|/2 dx dy .

Therefore, in regard to (3.9) and thanks to non-negativity of Ω λ provided by Lemma 3.3, it is sufficient to prove the pointwise inequality (3.10) (∀x) e |x|/2 Ω λ (x, y)e |y|/2 dy ≤ 2 .

In fact, the left-hand-side above is independent of λ. To see this, note that M λ (−y) = 1−M λ (y), and so for any R > 0,

Since this expression is independent of λ for all R > 0, we obtain (M λ (y) − M 1 (y)) e |y|/2 dy = 0 .

Hence, it is enough to check the inequality (3.10) for λ = 1, i.e. Hence,

Taking R → ∞, and multiplying by e x/2 yields the desired bound (3.12) . Note that equality is achieved only in the limit x → ∞. Similarly, simplifying (3.11) for x < 0 and replacing x with −x, one obtains again (3.12) . This concludes the proof of the Poincaré inequality (3.1).

To analyse the stability of solutions to (2.11), we define v as the perturbation around the equilibrium,

Denoting w(t, y) := ∂ y v(t, y), the perturbation v satisfies

and the two conservation laws rewrite as

Next, we introduce the weigthed energy functionals

|∂ y w(t, y)| 2 e −χ|y| dy .

In short, we denote

The weighted H 2 norm G is introduced in order to control the pointwise term w(t, 0) in (3.13b) that cannot be controlled in H 1 regularity. Our main estimate is the following. Proposition 3.4 (Entropy dissipation). The dissipation of F along solutions v(t, y) to equation (3.13) is given by

where J is a higher-order (genuinely non-linear) contribution defined in (4.3).

Proof. Differentiating F (t) along solutions of (3.13), substituting the evolution of ∂ t w, and integrating by parts, we havė

We evaluate the terms I 1 and I 2 separately. Expanding I 1 , we find

Dealing with I 2 , we separate the quadratic terms from the higher order contributions,

We can rearrange the quadratic term I ′ 2 as follows

This concludes the proof, with J(t) := J 1 (t) + J 2 (t).

In order to turn the entropy dissipation result Proposition 3.4 into exponential relaxation to equilibrium (see Section 4), we will make use of the improved Poincaré inequality (3.1) together with the following lemma:

Proof. We first perform some singular integration by parts in order to relate f (0) and f ′ (y).

By Young's inequality, we find

Raising it to the square, we obtain the result.

This section is devoted to the proof of Theorem 1.3, with the help of the previous energy estimates. We proceed in two steps. Firstly, we prove that the shape converges to equilibrium (meaning that u converges to a constant unit value). Secondly, we prove that the center x(t) converges to a finite value (not determined).

Proof of Theorem 1.3 -Step 1: convergence of the shape. Applying Lemma 3.5 to w ∈ L 1 λ with w ′ ∈ L 2 χ , and noting that w λ = 0 by the second conservation law, we have

Combining with Proposition 3.4, we find

where the higher-order contributions in J(t) reads as follows,

The first integral involves ∂ y w∂ y v L 1 χ , ∂ y wv L 1 χ , w∂ y v L 1 χ and wv L 1 χ , which are all controlled quadratically by E, F and G. In addition, the prefactor w(t, 0) is controlled by G 1/2 following equation (4.1). Similarly, the second contribution in (4.3) follows the same pattern, as both v(t) λ and v(t, 0) are controlled by interpolation using Lemma 3.5:

and the conservation law v(t) χ = 0. It follows that

for some constant C > 0 depending on χ and α, provided that F is small enough. Therefore,

for some constant C depending on χ and α, provided that F is small enough. Now, by the improved Poincaré inequality Proposition 3.1 and the second conservation law, we obtain

. Moreover, the classical Poincaré inequality with exponential weight together with the first conservation law yield

. The important point is to isolate the terms with higher derivatives, that is, the terms that are controlled by G(t). In fact, using the two inequalities above, the bound (4.4) can be simplified to the following estimate:

As a consequence, we can control the relaxation of F (t) using the estimate (4.2): assuming that F (0) is small enough, we get

We conclude that F (t), hence E(t), relaxes exponentially fast to zero, provided that F (0) is initially small enough. Moreover, the rate of convergence is asymptotically close to γ 0 as given in (1.10) . To conclude the proof of the exponential decay (1.9), we simply note that

Remark 4.1. It follows from (4.6) that the asymptotic rate of convergence γ can be chosen arbitrarily close to the upper bound γ 0 , however at the expense of increasing the prefactor as mentioned in Theorem 1.3. Additionally, it also means that in Theorem 1.3 we can choose any ε 0 that satisfies 0 < ε 0 < 4γ 2 0 /C 2 where C > 0 as defined in (4.6). More precisely, if initially ρ 0 −ρ ∞ 2 H 1 χ is close to the upper bound, i.e. F (0) = 4 C 2 (γ 0 − δ) 2 with δ > 0 arbitrarily small, then the above argument still provides decay to equilibrium, with a potentially very slow rate of convergence 0 < γ ′ ≤ δ. It follows that for any arbitrarily small 0 < ε < ε 0 , there exists a large enough time T ε > 0 such that F (T ε ) ≤ ε. From (4.6) we conclude

that is, the prefactor becomes large as T ε → ∞ together with γ → γ 0 .

Proof of Theorem 1.3 -Step 2: convergence of the center. We recall the dynamics (2.8)-(2.10) of the center x(t):ẋ

.

On the one hand, we have (2.11b):

By the same argument as above, −∂ 2 y S(t, 0) is bounded below uniformly for t > 0 provided that F (0) is small enough, simply because the term λv(t, 0) − √ α v(t) λ is controlled by F (t) which is exponentially decaying in that case.

On the other hand, the value w(t, 0) is estimated by the following basic inequality: It follows that

Therefore,

where the last estimate follows from Young's inequality. This concludes the proof of the bound (4.8).

We deduce immediately thatẋ is integrable, since F is exponentially decaying, whereas G is part of the dissipation ofḞ by (4.2) and (4.5). 

We performed a numerical investigation of system (1.1a)-(1.1b). We used a Lagrangian formulation, based on the inverse of the cumulative distribution function of ρ, that is, we define Π and its inverse X such that:

By using the formulation of S as a convolution with the fundamental solution of the elliptic problem (1.1b), we obtain the following closed equation on X(t, η):

sign (X(t, η) − X(t,η)) e − √ α|X(t,η)−X(t,η)| dη .

This formulation is well suited for the numerical simulation of interacting particles in one dimension of space, see e.g. [11, 6] . We discretize the cumulative mass η i = i∆η, for some step mass ∆η > 0. The discretization in mass rather than space provides more accuracy in regions of higher density. We use an implicit Euler scheme for the diffusion part, and an explicit scheme for the interaction part of the equation. Being given the solution X t (i) at time t and at cumulative mass η = i∆η, we compute X t+∆t (i) by solving

using the ©Octave function fsolve. Then the associated discretized density ρ is recovered from X by doing the reverse transformation using ρ(t, X t (i)) = 2∆η/ (X t (i + 1) − X t (i − 1)).

Typical results are shown in Figure 1 for an asymmetrical initial density ρ 0 which is a relatively large perturbation of the equilibrium profile, associated with a initial signal S 0 that does not satisfy the condition of having a unique critical point (as in Proposition 1.1), see Figure 1 (d). However, after some transient period, it does admit a unique critical point. This suggests that the convergence holds beyond a perturbative regime as analyzed in the present work. We can also see from the numerical results that the theoretical rate of convergence is slightly underestimated.

We can draw a couple of perspectives from this numerical study. (i) It would be interesting to prove that no more than a critical point of S can persist in the long time asymptotics, so that any initial configuration falls into the scope of Proposition 1.1 after some time, as in Figure  1 . (ii) The rate of convergence (1.10) can certainly be improved with a more cautious analysis, and the identification of a suitable functional inequality. and we can rewrite the left-hand side of the Poincaré inequality (3.1) as

Proof. We assume χ = 1 without loss of generality and denote Ω λ = Ω λ,1 . We can reformulate the left-hand-side of (3.1) as follows:

where we changed the order of integration to obtain the last line. Denoting

and expanding the last square, we find

wherẽ Ω λ (y, x 1 , x 2 ) = M λ (x 1 )M λ (x 2 )1 x 1 <y 1 x 2 <y − M λ (x 1 )N λ (x 2 )1 x 1 <y 1 x 2 >y − N λ (x 1 )M λ (x 2 )1 x 1 >y 1 x 2 <y + N λ (x 1 )N λ (x 2 )1 x 1 >y 1 x 2 >y .

Let us consider the case x 1 < x 2 (the argument for x 1 > x 2 is similar), and using the identity N λ = 1 − M λ , we can simplify the expression above:

Ω λ (y, x 1 , x 2 ) = M λ (x 1 )M λ (x 2 )1 x 2 <y − M λ (x 1 )(1 − M λ (x 2 ))1 x 1 <y<x 2

By exchanging the order of integration, we find

where Ω λ (x 1 , x 2 ) = 1 2 Ω λ (y, Next, we show Ω λ ≥ 0. It is enough to show non-negativity for x < y by symmetry. By direct calculation,

Suppose on the one hand that y > 0, then Ω λ (x, y) ≥ M λ (x)(M λ (y) − M 1 (y)) ≥ 0, because λ → M λ (y) is increasing for y > 0. Suppose on the other hand that y < 0, hence x < 0 as well (recall x < y by assumption). Therefore,

which is non-negative because λ → M λ (x) is decreasing for x < 0. This concludes the proof.

Appendix B. L 1 -stability in the case α = 0

In the case α = 0 dynamical arguments together with known stability results for shock waves [22] provide global L 1 stability of solutions. We sketch the argument here for a more general signal response function φ : R → R with the biologically reasonable assumptions that φ is odd, bounded, non-increasing and regular enough. For u[∂ x S] denoting the chemotactic flux, we consider the slightly more general bacterial chemotaxis model

Given a symmetric velocity set V ⊂ R, the macroscopic flux u can be related to the microscopic behavior of individual cells via the signal response function, u[∂ x S](t, x) = 1 |V | v∈V vφ(v∂ x S(t, x))dv , an expression that was derived rigorously from a mesoscopic model in [20] by means of parabolic scaling techniques. Choosing the stiff response function φ(x) = −sign(x) yields model (1.1) considered in this work (with χ := v∈V |v| dv/|V |). What allows us to handle the case α = 0 with an alternative method is a reformulation of system (B.1) as a viscous scalar conservation law (SCL). Substituting ρ = −∂ 2

x S into (B.1a), denoting z(t, x) := ∂ x S(t, x) , f (r) := − 1 |V | v∈V Φ(vr) dv with Φ the antiderivative of the signal response function φ, and integrating (B.1a) once w.r.t.

x, we obtain the viscous SCL (B.2) ∂ t z + ∂ x f (z) = ∂ 2 x z . The fundamental solution of (B.1b) is given by K 0 (x) = −|x|/2, and since we fixed the cell mass at ρ dx = 2/χ, we obtain the far-field conditions for any t > 0, where the integration constant is determined by the far-field condition and using the fact that Φ is even. Whilst no explicit stationary state can be found, existence follows from an implicit argument.

Toward a mathematical theory of Keller-Segel models of pattern formation in biological tissues

On the parabolic-elliptic Patlak-Keller-Segel system in dimension 2 and higher

Exponential Integrability and Transportation Cost Related to Logarithmic Sobolev Inequalities

Hypocoercivity and sub-exponential local equilibria

Chemotactic waves of bacteria at the mesoscale

Particle approximation of the one dimensional Keller-Segel equation, stability and rigidity of the blow-up

Confinement by biased velocity jumps: Aggregation of Escherichia coli. Kinetic and Related Models

Model hierarchies for cell aggregation by chemotaxis

Derivation of hyperbolic models for chemosensitive movement

Analysis of chemotactic bacterial distributions in population migration assays using a mathematical model applicable to steep or shallow attractant gradients

Lagrangian numerical approximations to one-dimensional convolution-diffusion equations

The diffusion limit of transport equations derived from velocity-jump processes

A users guide to PDE models for chemotaxis

Quand est-ce que des bornes de Hardy permettent de calculer une constante de Poincar exacte sur la droite ? Annales de la facult des sciences de

On a linear runs and tumbles equation

Motility of Escherichia coli cells in clusters formed by chemotactic aggregation

L 1 stability of travelling waves with applications to convective porous media flow

Mathematical models for chemotaxis and their applications in self-organisation phenomena

Transport models for chemotactic cell populations based on individual cell behavior

Mathematical Description of Bacterial Traveling Pulses

Directional persistence of chemotactic bacteria in a traveling concentration wave

L 1 -stability of nonlinear waves in scalar conservation laws

Overview of Mathematical Approaches Used to Model Bacterial Chemotaxis II: Bacterial Populations

Let V ⊂ R be the velocity set (which is a symmetric open interval), and let φ : R → R be odd, non-increasing and in C −1 (R). Then system (B.1) admits at least one stationary state (ρ ∞ , S ∞ )

This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No 639638). FH was partially supported by the von Karman postdoctoral instructorship at California Institute of Technology, and through the Engineering and Physical Sciences Research Council (UK) grant number EP/H023348/1 for the University of Cambridge Centre for Doctoral Training, the Cambridge Centre for Analysis. The authors benefited from fruitful discussion with Jean Dolbeault and Ivan Gentil about Proposition 3.1. FH is grateful to Camille, Marine and Constance Bichet, and Joachim Schmitz-Justen and Rita Zimmermann for their hospitality during the SARS-CoV-2 outbreak that allowed to finish this project.

Appendix A. Reformulation of the improved Poincaré inequalitywhere M λ denotes the cumulative density function,Then Ω λ,χ is non-negative and symmetric,Applying the mean value theorem to the continuous function g :. W.l.o.g. set Φ(0) = 0, and so Φ ≤ 0 and v * = 0. F is a convex function of Z ∞ (since φ is assumed non-increasing) with the two zeros z + and z − , and it follows that the above system has two equilibria, z + = −1/χ being an attractor, and z − = 1/χ a repeller. For −1/χ < Z ∞ < 1/χ, the derivative Z ′ ∞ is negative and so has the desired behavior ensuring that the far-field condition is satisfied. As the two stationary points lie on the same trajectory in phase space, a qualitative analysis of the phase portrait provides existence of a solution Z ∞ satisfying (B.4).Using the reformulation as a viscous SCL above, we can show asymptotic stability with respect to the L 1 -distance.This result yields more than just stability in the sense that a solution z(t, x) initially close to the stationary state Z ∞ remains close for all times. In fact, Theorem B.2 holds for any initial perturbation, providing global L 1 -stability.The proof of Theorem B.2 makes use of the linear C 0 -semigroup T t : L ∞ (R) → L ∞ (R) corresponding to the Cauchy problem for (B.2), sending an initial datum z 0 ∈ L ∞ (R) to its corresponding solution z(t, x) with z(0, x) = z 0 (x). If f in (B.2) is at least C 2 , the semigroup (T t ) t≥0 enjoys the four Co-Properties: comparison (a ≤ b a.e. implies T t a ≤ T t b a.e), contraction ( T t a − T t b 1 ≤ a − b 1 ), conservation ( T t a dx = a dx) and constants (if a is constant, then T t a = a). These properties follow from the fact that the flux term ∂ x f (z) in (B.2) may be interpreated as a lower order perturbation of the heat equation ∂ t z = ∂ 2x z, see [22] . The contraction property immediately yields the decayand it remains to show that the limit is indeed zero as claimed.Proof. The key approach is to first focus on the case where the initial datum lies between two shifts of the profile Z ∞ as this case can be treated by means of dynamical systems theory using compactness of trajectories, ω-limits and Lasalle's invariance principle (Theorem 3 in [22, Chapter 7, Section 3.2], which relies on Lemma 4 in the same section). The general idea for this result is due to Osher and Ralston, see [17] . It then remains to show that the set of initial conditions considered in Theorem B.2 is included in the L 1 -closure of the set of L ∞ functions sandwiched between two arbitrary shifts of Z ∞ . Indeed, for z 0 ∈ L 1 + Z ∞ such that z + ≤ z 0 ≤ z − , the integral ∞ y |z 0 (x) − Z ∞ (x)| dx vanishes as y → ∞, and similarly for the distance towards −∞. So for any ε > 0 there exist s, t ∈ R such thatFor η > 0 small enough such that η < ε/(3χ|s − t| −1 ), we definez(x) equal to (1 − χη)z 0 (x) on (t, s) and equal to Z ∞ elsewhere. ThenAnd for big enough k ∈ R and small enough j ∈ R, we have