Confinement by biased velocity jumps: aggregation of Escherichia coli

We investigate a linear kinetic equation derived from a velocity jump process modelling bacterial chemotaxis in the presence of an external chemical signal centered at the origin. We prove the existence of a positive equilibrium distribution with an exponential decay at infinity. We deduce a hypocoercivity result, namely: the solution of the Cauchy problem converges exponentially fast towards the stationary state. The strategy follows [J. Dolbeault, C. Mouhot, and C. Schmeiser, Hypocoercivity for linear kinetic equations conserving mass, Trans. AMS 2014]. The novelty here is that the equilibrium does not belong to the null spaces of the collision operator and of the transport operator. From a modelling viewpoint it is related to the observation that exponential confinement is generated by a spatially inhomogeneous bias in the velocity jump process.


Introduction
Unbiased velocity randomization by a jump process or by Brownian motion combined with acceleration by the force field produced by a confining potential can lead to convergence to an invariant probability measure, if the confinement is strong enough to balance the dispersive effect of velocity randomization. For kinetic transport models of this kind, convergence to confined equilibria has been studied extensively leading to strong convergence results with algebraic [10,14] and later also exponential [17,30] convergence rates. This is strongly related to the corresponding macroscopic description by Fokker-Planck equations of drift-diffusion type [4].
In this work a related type of particle dynamics is considered, where confinement is achieved by a biased velocity jump process, where the bias replaces the confining acceleration field. The motivation comes from kinetic transport models for the chemotactic motility of microorganisms driven by gradients of chemo-attractors. The prototypical example is the bacterium Escheria coli, whose swimming pattern has been described as run-and-tumble [7,9], meaning that periods of straight running alternate with periods of reorientation (tumbling). Since typically tumble-periods are short compared to run-periods, models with instantaneous velocity jumps seem reasonable, but see [33,13]. In the presence of a spatial chemo-attractant gradient, this stochastic process is biased upwards the gradient, although E. coli is too small to reliably measure the gradient along its length. An explanation for this phenomenon is that E. coli is able to measure gradients in time along its path and increases its tumbling frequency, if it experiences decreasing chemo-attractant concentrations. This produces the desired drift, even if the outward velocity after tumbling events is unbiased [8].
This and similar motility behavior types have been incorporated into kinetic transport models of chemotaxis [2,41,16,18,48], and a connection to versions of the classical Patlak-Keller-Segel (PKS-) model [35,42] has been made by the macroscopic diffusion limit [3,31]. If chemotaxis acts as a means of signalling between cells, it can be responsible for various types of pattern formation with aggregation as the most important and basic outcome. Two examples of pattern formation are stable clusters of bacteria [38] and travelling pulses of bacterial colonies [1,43,44]. For corresponding kinetic transport models with nonlinear coupling to a reaction-diffusion model for the chemo-attractant, the macroscopic diffusion limit produces nonlinear versions of the PKSmodel [12,28,43]. We emphasize that hydrodynamic limits can also be derived from the kinetic model [16,19,32].
Aggregation patterns of signalling E. coli have been observed and simulated in [38]. The stochastic simulations have considered a prescribed peaked stationary chemo-attractant concentration. A kinetic model corresponding to these simulations will be considered here: The velocity set V is supposed to be bounded and rotationally symmetric. The right hand side of (1.1) is the turning operator. It describes the velocity changes due to tumbling. The turning kernel K(x, v) is the rate of changing from velocity v to a different velocity v at position x. The x-dependence contains the influence of the chemo-attractant. The fact that the turning kernel only depends on the incoming (pre-tumbling) velocity means a complete randomization of velocity at tumbling events. This is actually not describing the experimental evidence precisely. Whereas independence of the outgoing velocity distribution from the chemoattractant gradient seems to be a valid assumption, some directional persistence of E. coli has been observed, i.e. an outgoing velocity distribution biased by the incoming velocity [7]. This effect might have important quantitative consequences on the efficiency of chemotactic foraging strategies [45,37,40], but apparently does not change the qualitative picture. The main deficiency of (1.1) as a model for the experiments of [38] is the lack of a nonlinear coupling with an equation for the chemo-attractant, describing production by the cells, diffusion, and decay. The restriction to the linear problem has purely mathematical reasons. The following sections will show that proving the existence and stability of aggregated stationary solutions already poses significant difficulties for a simple version of (1.1). Extensions to nonlinear models are the subject of ongoing investigations. Numerical studies [46,11,20] show that convergence to a steady state can be expected under suitable assumptions. We highlight the well-balanced numerical scheme proposed by Gosse in [24,Section 9.4]. This focuses on the local kinetic equilibria of (1.1) in a spatial finite volume with inflow boundary conditions, very much in the spirit of the present work.
The turning rate takes the larger value 1 + χ for cells moving away from x = 0, and the smaller value 1 − χ for cells moving towards x = 0. Equilibrium distributions of the turning operator, i.e. , multiples of 1/K, have a jump at x = 0 and are therefore not solutions of (1.1). Stationary solutions have to balance the turning operator with the transport operator v∂ x . Existence and uniqueness (up to a multiplicative constant) of a stationary solution will be proven in the following section. As an illustration we have a short look at the two-velocity model, also known as the Cattaneo model for chemotaxis [26,31,27,29,39,23,32], Here, steady states are easily found explicitly as multiples of g + (x) = g − (x) = e −2χ|x| . The exponential decay with respect to position carries over to the model with V = [−1/2, 1/2], which is also proven in the following section. We also refer to [34] for the analysis of a discrete velocity jump process and its stationary distribution based on the WKB expansion. Section 3 is concerned with the decay to equilibrium as t → ∞, employing the modified entropy approach of [17] for abstract hypocoercive operators [47]. This approach is based on a decomposition of the generator of the dynamics into a symmetric negative semidefinite operator and a skewsymmetric operator, where the latter needs to satisfy a certain coercivity condition on the null space of the former (called 'macroscopic coercivity' in [17], compare also to 'instability of the hydrodynamic description' in [15]), and the former needs to be coercive on the complement of its null space ('microscopic coercivity' in [17]). The set of equilibria has to be the intersection of the null spaces of both operators. In the present case, this does not permit the splitting into turning and transport operators, as usual for kinetic equations. The main preparatory step is therefore the definition of an appropriate splitting, before the modified entropy method is applied.
The terminology 'microscopic' and 'macroscopic' refers to an asymptotic limit, where the characteristic time scale of the symmetric operator is assumed to be much smaller than that of the skew-symmetric operator. This limit is carried out in the first part of Section 4. Typically the separation of time scales can be achieved by a macroscopic re-scaling of length in kinetic equations, which is not the case in the present situation by the redefinition of the splitting between operators. The asymptotic limit has therefore not much biological relevance. A second, more realistic macroscopic limit is carried out in the second part of Section 4, where smallness of the parameter χ is assumed, and length and time is rescaled diffusively. Both macroscopic limits produce drift-diffusion equations, whose diffusivities and convection velocities are different, but with the same qualitative behavior. This brings us back to the beginning of the introduction, since it shows that the macroscopic behavior created by biased velocity jumps is the same as for unbiased jumps combined with a confining potential.

Existence and exponential decay of stationary solutions
We seek a nonnegative stationary state g ∈ L 1 + (R × V ) satisfying For simplicity, we assume V = [−1/2, 1/2], such that |V | = 1. More generally we may assume that we are given a probability measure dν(v) which is compactly supported. Then we replace dv with dν(v) in the following. We denote We make the following assumptions on the turning kernel K(x, v): , and piecewise continuous, with a possible jump only at v = 0. These conditions are satisfied by our main example (1.2). Theorem 1. Under the assumptions (H1-H4) there exists a nontrivial stationary state g(x, v) solution to (2.1). It is positive, bounded and symmetric: g(x, v) = g(−x, −v). There exists α > 0, a positive velocity profile G, and a constant C > 0 such that In addition, we can state precisely the asymptotic behaviour of the stationary state g as x → +∞, namely variable separates in the limit: This is a specific feature of the so-called Milne problem in radiative transfer theory [6,21]. The analogy between the existence of a stationary state for equation (2.1) and the Milne problem is the cornerstone of the proof of Theorem 1.
Proposition 2. The stationary state constructed in Theorem 1 converges exponentially fast towards a multiple of the asymptotic profile e −αx G(v) as x → +∞. More precisely, there exists a constant C 0 such that the following estimate holds true, where β is the positive root of 1 V v 2 G 2 dv , and the constant H(g) is given by the formula: Remark 3. The higher-dimensional case is left open. We refer to [20] for accurate numerical simulations of (1.1) in the two-dimensional case, where the velocity set is the unit sphere S 1 . In this work the authors clearly observe convergence towards a spherically symmetric stationary state.
Proof of Theorem 1. The proof is divided into several steps. The first step consists in deriving the correct asymptotic behavior of g as |x| → +∞. In the second step we make the link between our stationary problem and the so-called conservative Milne problem in the half-space [6]. The core of the proof is a regularization process in the velocity variable at x = 0 (step 3), which enables us to apply the Krein-Rutman Theorem (step 4). The proof shares some similarity with analogous problems in homogenization theory (see for instance [5,25]). However the connection with the Milne problem is new up to our knowledge.
Step#1. Exponential decay as x → +∞. We make the following ansatz Substituting this ansatz in (2.1) yields Clearly the profile G is characterized up to a constant factor. We opt w.l.o.g. for the renormalization Therefore the exponent α is characterized by the dispersion relation: We now prove that this defines a unique α > 0 such that (∀v ) −1 dv < 0 (the average speed is negative on the far right: this is the confinement effect). (iii) The function J is convex on the admissible range of α.
As a consequence there exists a unique α ∈ 0, inf such that J(α) = 1.
By symmetry, a similar ansatz can be made on the left side: Step#2. Connection with the conservative Milne problem. For x ≥ 0 we define u by g(x, v) = e −αx G(v)u(x, v). We deduce from (2.1) that it satisfies the following equation, The profile g(0, v) = G(v)u(0, v) is unknown, of course. On the other hand, its knowledge is sufficient to reconstruct the entire function u. This is the purpose of the Milne problem. It states that for a given profile u(0, v) defined for v > 0 only, there exists a unique bounded function u, defined over R + × V , solution of (2.7).
satisfying the pointwise estimate, Proof of Lemma 4. This result is classical (see [6] and references therein). We recall the main lines of the proof for the sake of completeness.
For ε > 0 we associate the perturbed problem which possesses a unique solution in L ∞ (R + × V ), as can be proven using a fixed point argument as explained below. The Duhamel formulation for (2.10) reads where λ ε (v) = G(v) −1 +ε, and the macroscopic quantity is defined by For a given inflow data ϕ, we associate the map T ε from L ∞ (R + ) into itself, where u ε is defined by (2.11). Then T ε is a contraction with rate The fixed point of this map is a solution of (2.10). The maximum principle applied to (2.10) implies that u ε satisfies (2.9). Therefore, as ε → 0, we can extract a subsequence that converges in L ∞ (R + )-weak * to a solution u of (2.8), which satisfies Next we show that any bounded solution of (2.7) satisfies the estimate (2.12), implying uniqueness. We define U (x, v) = G(v) K(v) e αx , which also satisfies the differential equation in (2.8). It corresponds in fact to the function g = 1 K , which is a trivial solution of (2.1) for x > 0. For ε > 0 we consider w = u−εU . It satisfies (2.7) with inflow data ϕ ε = ϕ−εU (0, ·). Since u is bounded and −εU → −∞ as x → +∞, it is clear that w has a maximum in R + ×V . Notice that w is not necessarily continuous at v = 0, but piecewise continuity is sufficient. The maximum cannot be attained in (0, ∞) × V by the maximum principle. Therefore it is attained at x = 0 and we get A similar estimate holds for −u − εU . Letting ε → 0 we deduce that u satisfies (2.9), and therefore also (2.12).
Step#3. Compactness. The expected symmetry property g(0, v) = g(0, −v) of the equilibrium distribution motivates the definition of the fixed point operator B : with the Albedo operator (Aϕ)(v) = u(0, v), where u is the unique bounded solution of (2.8).
Lemma 5. The operator B is compact and positive.
Proof. To prove compactness, we first define the macroscopic quantity as above, We have u ∈ L ∞ (R + × V ) and from (2.8) also v∂ x u ∈ L ∞ (R + × V ). The one-dimensional averaging lemma [22] implies A ∈ C 0,θ (R + ) for all θ ∈ (0, 1), with the corresponding Hölder semi-norm bounded in terms of ϕ L ∞ . We deduce from the Duhamel representation formula that u(0, v) is uniformly continuous for v < 0, with a modulus of continuity which depends only on ϕ L ∞ and the modulus of continuity of vG(v) on V − . Positivity is immediate from the Duhamel formula since A(x) > 0 for x > 0, as soon as ϕ ≥ 0 and ϕ = 0.
Step#4. Conclusion. In order to apply the Krein-Rutman Theorem [36], we consider the restriction of B to continuous functions on V + , since the interior of L ∞ (V + ) is empty. Lemma 5 also holds for The Krein-Rutman theorem states that the operator B| C 0 (V + ) possesses a simple dominant eigenvalue λ ∈ R together with a positive eigenfunction ϕ: Bϕ = λϕ. The conservation property for (2.1) yields λ = 1 by the following argument: We denote by u ∈ L ∞ (R + × V ) the solution of the Milne problem (2.8) with inflow data u(0, v) = ϕ(v). We define accordingly g( implying λ = 1 by the positivity of g. It is straightforward to check that g, symmetrically extended to R × V , is a solution of (2.1), satisfying (2.2) as a consequence of (2.9).
Proof of Proposition 2. We adapt the method of [6]. We shall use a quantitative energy/energy dissipation approach with respect to the space variable. First, we integrate (2.7) against K + G 2 , in order to derive a non trivial conservation, This comes in addition to the zero-flux relation which is a straightforward consequence of equation (2.1) after integration with respect to the velocity variable. We observe that the relation (2.14) combined with (2.15) and (2.5) yields Secondly, we multiply (2.7) by 2K + G 2 (u − H(u)), where the constant H(u) is defined such that We obtain H(u)), and we have used the conservation (2.16). We define two auxiliary quantities, The dissipation inequality (2.17) reads ∂ x J(x) + 2κE(x) ≤ 0. On the other hand, multiplying (2.7) by 2vG 2 (u − H(u)), we get, where we have used (2.5), (2.15) and V vG dv = 0. This reads also ∂ x E(x) = 2αE(x) − 2J(x). All together this yields the damped second order inequality Complemented with the information that E is bounded, this is sufficient to prove that E decays exponentially fast. In fact, let β > 0 be the positive root of 1 2 β 2 + αβ − 2κ = 0, the functioñ We deduce that for all 0 ≤ x < y, we havẽ Hence we conclude that, Since E(y) is uniformly bounded, letting y → +∞ this yields that E(x) decays exponentially fast, i.e. the estimate (2.3) holds true up to a modification of the constant C 0 , and the abuse of notation H(g) = H(u).

Decay to equilibrium by hypocoercivity
Since (1.1) conserves the total mass, we expect lim t→∞ f (·, ·, t) = µ ∞ g (with g solving (2.1)), where the constant µ ∞ is chosen such that µ ∞ R×V g(x, v)dv dx = R×V f (x, v, 0)dv dx. Replacing f by f − µ ∞ g, we may assume µ ∞ = 0, and therefore in the following. The decay to equilibrium will be proved by employing the abstract procedure of [17]. It is based on a situation, where the equilibrium lies in the intersection of the null spaces of the collision and the transport operator. However, the equilibrium g is not in the null spaces of the collision operator Qf := V K f dv − Kf and of the transport operator v∂ x . Therefore, as a first step, collision and transport operators will be redefined.
Multiplying (1.1) by f /g, we get: where · is the norm on induced by the scalar product This motivates the definition of the symmetrized collision operator with the same dissipation: Hence we rewrite (1.1) as ∂ t f + Tf = Lf with the modified transport operator It is easily checked that Tg = Lg = 0 and that T is skew symmetric with respect to ·, · . Obviously, for fixed x, the null space of L is spanned by g(x, ·). The orthogonal projection to N (L) is given by We also observe implying ΠTΠ = 0 (by the property V vg dv = 0 of the equilibrium distribution), which is Assumption (H3) in the abstract setting of [17], Section 1.3. The so called 'microscopic coercivity' Assumption (H1) is the subject of the following result. Proof.
− Lf, f ≥ K min 2 The next step is 'macroscopic coercivity' (Assumption (H2) in [17]). It relies on the asymptotic behavior of g as |x| → ∞. By Theorem 1, the equilibrium distribution satisfies where α is the unique positive solution of the dispersion relation (2.6), and u min , u max are positive constants.

Proof. A straightforward computation gives
By the boundedness of the velocity space, ρ g and m g satisfy estimates of the form (3.5), implying the existence of a positive constant c such that We claim that the measure ρ g dx permits a Poincaré inequality. This can be proved via bounds on the spectrum of Schrödinger operators. It is a consequence of the fact that ρ g can be bounded from above and below by multiples of a smooth function, whose logarithm is asymptotically linear as |x| → ∞ and Theorem A.1 of [47]. Therefore, there exists a positive constant λ M such that The approach of [17] relies on the 'modified entropy' and with a small positive constant ε. Its time derivative is given by d dt The first two terms on the right hand side are responsible for the decay of H[f ], noting microscopic coercivity (Lemma 6) and the fact that the macroscopic coercivity result of Lemma 7 implies Proof. The result follows from the estimate It remains to prove boundedness of AT(1−Π), which is equivalent to an elliptic regularity estimate. Proof. As in [17] we work on the adjoint. It is sufficient to prove boundedness of (AT) * = −T 2 Π(1 + (TΠ) * TΠ) −1 .
Introducing ϕ = (1 + (TΠ) * TΠ) −1 f , the scalar product of the equivalent equation with ϕ leads to the estimate ϕ , TΠϕ ≤ f . Application of Π to (3.6) gives We shall have to estimate the norm of satisfying |T 2 Πϕ| ≤ Ce −α|x| ∂ 2 x u ϕ + |∂ x u ϕ | , in terms of f , which is a weighted L 2 → H 2 regularisation result for (3.7). The first order derivative has already been taken care of by the bound for TΠϕ = gv∂ x u ϕ . We shall also need where we have used the equation (2.1), satisfied by g, and (3.5). Now (3.7) is rewritten as , and the proof is completed by taking L 2 -norms with the weight e α|x| , noting that m g ≥ ce −α|x| .
We have completed the verification of the assumptions of Theorem 2 of [17] and arrive at our main result: Theorem 10. Let a stationary positive solution g of (1.1) (unique up to a constant factor) satisfy (3.5), let f I ∈ L 2 (dv dx/g) (⊂ L 1 (dv dx)), and let Then the solution of (1.1) subject to f (t = 0) = f I satisfies with positive constants C and λ, only depending on χ ∈ (0, 1).

Macroscopic limits
The 'macroscopic limit' corresponding to the modified entropy method. The separation into microscopic and macroscopic contributions employed in the previous section can be motivated by a macroscopic limit, based on the assumption of a separation of time scales related to the collision and the transport operators. With an appropriate (parabolic) rescaling of time, this leads to studying the limit as ε → 0+ in Whereas for the standard transport operator v∂ x the scale separation can be achieved by a length rescaling, the introduction of the 'Knudsen number' ε is completely artificial in the present situation, since the modified collision and transport operators contain identical terms whose different weighting cannot be justified by scaling arguments. This is the reason for the quotation marks in the title of this subsection. The limit ε → 0 in the abstract equation (4.1) has been carried out formally in [17] under the assumptions ΠTΠ = 0, already used above, and that the restriction of the collision operator L to (1 − Π)H possesses an inverse J. The formal limits f 0 of solutions of (4.1) satisfy Note that the macroscopic coercivity estimate in Lemma 7 is related to the simplified version −(TΠ) * (TΠ) of the macroscopic evolution operator. Since f 0 (x, v, t) = ρ 0 (x, t)g(x, v)/ρ g (x), the above evolution equation is equivalent to an equation for ρ 0 . The additional requirement f ∈ (1 − Π)H determines µ uniquely.
Proof. With ν = V Kf dv, µ = V f /g dv, the equation Lf = h can be rewritten as Division by g and integration with respect to v gives an equation for ν − λµ leading to the claimed result. A second equation for µ and ν can be obtained by multiplication by K and integration. The coefficient matrix of the resulting system is non-invertible, and the solvability condition turns out to be V h dv = 0, i.e. h ∈ (1 − Π)H.
Proposition 12. Let L and T be defined by (3.2) and, respectively, (3.3). Then the equation (4.2) is equivalent to where D = D(x) is given by Remark 13. Note that the macroscopic evolution operator and the simplification −(TΠ) * TΠ only differ by the diffusivities D vs. m g . Under the assumption (3.5), both have the same behaviour as |x| → ∞.

Multiplication of the second equation by v and integration finally gives
The equilibrium distributions ρ g and e −3|x| of the two macroscopic limits share the exponential behavior for |x| → ∞ up to the decay rate (by (3.5)), which is a consequence of the shared asymptotic behavior of the diffusivities D/ρ g vs. 1/12 and drift velocities −D∂ x (1/ρ g ) vs. −sign (x)/4.