Aggregated Steady States of a Kinetic Model for Chemotaxis

A kinetic chemotaxis model with attractive interaction by quasistationary chemical signalling is considered. The special choice of the turning operator, with velocity jumps biased towards the chemical concentration gradient, permits closed ODE systems for moments of the distribution function of arbitrary order. The system for second order moments exhibits a critical mass phenomeneon. The main result is existence of an aggregated steady state for supercritical mass.


1.
Introduction. Chemotaxis, the movement of biological agents influenced by gradients of chemical concentrations, is a ubiquitous process in biological systems. On the other hand, the production or degradation of chemicals is at the basis of standard signalling mechanisms between individuals. This produces a nonlinear feedback which, together with chemotactic motility, may drive self-organization processes in groups of agents.
A typical example, observed in many bacterial and amoeboid species, is aggregation driven by the production of a diffusible chemical, and chemotactic movement biased towards the direction of the gradient of the chemical concentration (see the large literature on Dictyostelium discoideum or, for bacteria, [9]). Since motility usually has a random component, it is a standard question in this situation, if the attractive mechanism is strong enough to overcome the dispersion caused by the random motility component.
The type of mathematical models mostly depends on the nature of the motility process. The standard assumption of Brownian motion with drift, the latter determined by chemotaxis, leads to a version of the classical Patlak-Keller-Segel (PKS) model [11,8], where a convection-diffusion equation for the agent density is coupled with a reaction-diffusion equation for the chemical concentration. For certain bacterial species a description by a velocity jump process is more appropriate, whence the convection-diffusion equation of the PKS model is replaced by a kinetic transport equation [10]. The PKS model can typically be recovered as a macroscopic limit [4,7]. However, some observed phenomena are only explainable by kinetic models [12].
Three types of long time behavior can be observed in mathematical models. If the random motion of agents dominates, this leads to dispersion, i.e. the same qualitative behavior as for the heat equation. For dominating attractive effects, the agent density either has a nontrivial aggregated long-time limit, or it blows up in finite time, typically in a concentration event. The two-dimensional parabolic-elliptic PKS model (i.e. with a quasistationary equation for the chemical concentration) has been thoroughly analyzed with respect to these questions. It shows a critical mass phenomenon: Among the initial data with finite variance those with the total mass below a critical value lead to dispersion and those with supercritical mass to finite time blow-up [1]. At the blow-up time strong solutions cease to exist, but a continuation by measure solutions is possible as limiting case of regularized models [6,13]. The corresponding dichotomy has been shown to exist also in kinetic transport models [2]. The situation for other versions of the PKS model and, in particular, for kinetic models is less clear.
Motivated by experimental results for E. coli [9], a linear kinetic model with given aggregated chemical concentration has been analyzed in [3]. The existence of a nontrivial steady state and its dynamic stability have been proven (the latter by employing the methodology of [5]). The present work can be seen as a continuation, where the nonlinear coupling with a quasistationary model for the chemical is added. The main result is a critical mass phenomenon, but with a dichotomy between dispersion and the existence of an aggregated steady state. Consider the system a one-dimensional kinetic model for chemotaxis, where the cells with phase space density f (x, v, t) and macroscopic density and flux, produce the chemoattractant with density S(x, t). The dynamics of the chemoattractant (diffusion, production, and decay) is assumed to be fast (and therefore modelled as quasistationary) compared to the dynamics of the cells. We consider two choices for the turning kernel: For both models, we assume α, β, γ, κ, D > 0. In Model A, cells decide about reorientation by scanning the chemoattractant density in the directions of possible post-turning velocities. In Model B, they scan in the direction of possible velocity changes. Whereas in Model A the outcome of a turning event is independent from the pre-turning velocity, Model B can be motivated by directional persistence of the agents. This obviously makes it harder to reach a target, and it is one of the results of this work that Model B does not lead to cell aggregation. Both models have not been derived systematically from the microscopic behavior of a particular cell type. However, they are reasonable from a qualitative point of view, and they have the remarkable mathematical property that the evolution of moments can be computed by solving linear constant coefficient ODEs without solving the full equations (see Section 3). We observe that the rescaling eliminates all parameters, i.e., (1), (2) becomes with We consider the Cauchy problem with initial conditions The solution S of (2) is defined as the convolution product of the decaying fundamental solution of −∂ 2 x + id with ρ: The initial datum is assumed to possess moments of up to second order: We choose a reference frame such that the first order moments vanish initially: We shall show that for Model A the (x-and v-) moments of f up to any fixed order satisfy closed systems of linear, constant coefficient ODEs. The system of second order moments exhibits a critical mass phenomenon. If the total mass M is below a critical value, the second order moments grow indefinitely with time, whereas for large enough mass they converge to finite values. The corresponding system for Model B always produces growing second order moments. Therefore we shall concentrate on model A after this observation. It turns out that also the higher order moment systems exhibit a critical mass phenomenon, however with the critical mass increasing with the moment order. Since stationary solutions may be the limits when time tends to infinity of the solutions to the Cauchy problem associated to (1), (2), this suggests a mass dependent decay behavior of the steady state. However, a precise characterization is still open.
The rest of this work is structured as follows: Existence and uniqueness of global solutions of the Cauchy problem is proven in Section 2 for Models A and B. In Section 3 it is shown for Model A that in principle all moments of the solution can be obtained by recursively solving systems of linear ODEs, where all moments of fixed order are coupled. The second order system exhibits a critical mass phenomenon. Whereas for small mass the second order moments tend to infinity, they converge to a steady state, if the mass is above a critical value, indicating the existence of an aggregated stationary solution in the latter case and dispersion in the former. For higher order systems, roughly speaking the critical mass increases. It is shown that at each order the system is stable for mass large enough. On the other hand, for each fixed mass, stability only holds up to a finite order. This indicates algebraic decay of possible stationary solutions with mass dependent decay. However, the details of the dynamics are unclear. Higher order systems might exhibit increasing and decaying modes as well as oscillations. For this reason, we cannot deduce convergence to a stationary solution from the moment systems. For Model B the situation is simpler and less interesting. The second order moments again solve a closed ODE system, whose solutions always tend to infinity. Therefore we conjecture that in Model B dispersion always dominates and restrict the further analysis to Model A. A formal asymptotics for large mass, equivalent to a macroscopic limit, of Model A is performed in Section 4. In the large time limit of the resulting macroscopic convection equation, concentration with repect to position and exponential decay with respect to velocity is obtained, indicating anisotropic decay behavior in phase space also for finite mass. The existence of a smooth steady state for Model A with supercritical mass (of the second order moment system) is proven in Section 5. The natural question of a precise characterization of the decay behavior of the stationary solution remains open. So far it resisted both our attempts of a direct asymptotic analysis and of a reconstruction from results on the moments. Another question, which remains the subject of future work, is a characterization of the long time behavior when dispersion dominates, i.e. for subcritical mass in Model A and for all cases in Model B. Finally, it should be possible to extend some of the results to corresponding multi-dimensional models.
Proof. For every T > 0, let It can be computed explicitly as For ρ ∈ X T nonnegativity of R A (ρ) is obvious, and the mass conservation property follows by integration of (13) with respect to x, implying R : is Lipschitz with Lipschitz constant 1. This implies for ρ 1 , ρ 2 ∈ X T , after a change of variables, Thus, for T < ln 2 M the map R A is a contraction on X T . This proves local solvability. As a consequence of the uniform bound in L 1 (R 2 ) the solution can be extended indefinitely in time steps of length T .
) to the Cauchy problem associated to Model B, i.e., with A solution f to the Cauchy problem associated to Model B is directly obtained as the limit of the increasing sequence (f j ) defined by f 0 = 0 and f j+1 given from f j as the solution to where ρ fj = f j dv. f j+1 is explicitly given from f j by , v , s)dv ds.
Consequently it can be proven by induction that (f j ) is nonnegative, and non de- and f j ≤ f j+1 by (17). Moreover, denoting by m j (t) = f j (x, v, t)dxdv and integrating (16) with respect to (x, v) ∈ R 2 , it holds m j+1 = m 2 j − M m j+1 , so that it can be proven by induction that And so, by the monotone convergence theorem, (f j ) converges in L 1 to a nonnegative function f . It implies that the limit in And so, f is a solution of (14). Conservation of mass from (14) implies that f (x, v, t)dxdv = M . The limit of (17) implies that f ∈ C([0, ∞), L 1 + (R 2 )). f is the unique nonnegative solution of (14) since its construction makes it minimal among the nonnegative solutions of (14). Then also S has finite moments up to order N , and with the following relations hold: Proof. The result can be shown by straightforward computations. We multiply the differential equation for S by |x| k and x k , and use to show the boundedness of the moments of S and a). Then b) is a consequence of the substitution y = x + v and of the binomial theorem. If we concentrate on moments of order N , then the lemma implies where LOT (Lower Order Terms) stands for terms only depending on moments of order lower than N . Now we introduce moments of solutions f of (3), (4) with respect to x and v: With the help of Lemma 1 and with we can derive differential equations for the moments. The first and obvious one is mass conservation:Ȧ 0,0 = 0 =⇒ A 0,0 = M . As a consequence, the turning operator of Model A after elimination of the unknown S by (6), can now be written as For the first order moments, we obtaiṅ because of (8). For the moments of order two it gets more interesting: Application of the Routh-Hurwitz criterion to the characteristic polynomial of the coefficient matrix shows that for M < 2 at least one positive eigenvalue exists, whereas for M > 2 all eigenvalues have negative real parts. Thus, in the latter case all solutions converge to the steady state For the higher order moments we only concentrate on the highest order terms on the right hand sides: for 0 ≤ n ≤ N . This is a linear ODE system with constant coefficients and an inhomogeneity only depending on lower order moments. This shows that all moments can be computed recursively.
If the coefficient matrices of all systems up to order N only have eigenvalues with negative real parts, then all moments of order up to N have finite limits as t → ∞. We have already shown above that for N = 2 this property holds, iff M > 2.
The characteristic polynomial of the order N coefficient matrix can be written as Denoting the roots of r N by µ 1 , . . . , µ N ∈ C (multiple roots allowed), we found N more roots of p N : Obviously, all the N + 1 roots we found have negative real parts for large enough M . Proof. For fixed λ and M , Therefore, for N large enough there exists λ > 0 such that sign(p N (λ)) = (−1) N . On the other hand, completing the proof. Combination of the existence theorem 1 with the previous results leads to the propagation of moments: Corollary 1. Let the assumptions of Theorem 1 hold and let (1 + |x| N + |v| N )f I ∈ L 1 (R 2 ) for an N ≥ 1. Then the solution f of (9), (10) satisfies (1 + |x| N + |v| N )f ∈ L ∞ loc (R + ; L 1 (R 2 )). If N ≥ 2 and M > 2, then (1 + |x| 2 + |v| 2 )f ∈ L ∞ (R + ; L 1 (R 2 )).
For Model B, the computations are similar but a little more involved. As for Model A, A 0,0 = M and A 1,0 = A 0,1 = 0 hold. The 2 nd order moments satisfy the closed ODE systemȦ By their definition and by the Cauchy-Schwarz inequality, for nonvanishing initial data f I their initial values satisfy A 2,0 (0)A 0,2 (0) > A 1,1 (0) 2 and A 2,0 (0), A 0,2 > 0. A straightforward computation gives This guarantees that A 2,0 and A 0,2 remain positive for all times. As a consequence, d dt (M A 2,0 + A 0,2 ) = 2M (M + A 2,0 ) ≥ 2M 2 , so M A 2,0 + A 0,2 tends to infinity. Now assume A 2,0 remains bounded. Then A 0,2 → ∞, implying by the second differential equation A 1,1 → ∞ and, thus, by the first equation the contradiction A 2,0 → ∞. Assume on the other hand that A 0,2 remains bounded. Then A 2,0 → ∞, implying by the second differential equation A 1,1 → −∞ and, thus, by the third equation the contradiction A 0,2 → ∞. So both A 2,0 and A 0,2 become unbounded as t → ∞, meaning that the chemotactic effect is not strong enough to prevent dispersion. This is not too surprising, since Model B only supports velocity changes in the direction of the target, whereas Model A supports post-turning velocities in this direction. For this reason we concentrate on Model A for the rest of this work.

Formal asymptotics for large mass.
With the rescaling f → M f , S → M S, Model A takes the form with M now taking the role of an inverse Knudsen number. The rescaled version of the steady states for the moments are which suggests an equilibrium state concentrating with respect to x as M → ∞.
Obviously, we have ρ 0 (x, t) → δ(x) as t → ∞ and, thus, This is in agreement with the limit as M → ∞ in (22).

Stationary solutions.
In this section we first prove in Theorem 3 the existence of even nonnegative L 1 -solutions to the stationary problem, then their C ∞regularity in Theorem 4.
Theorem 3. For any M > 2 there exists an even nonnegative distributional solu- which is also a mild solution satisfying Proof. The strategy of the proof is to first solve a truncated problem and then to pass to the limit. Large (resp. small) values of the position (resp. velocity) variable will be truncated. At the resulting boundaries of the position domain, reflecting boundary conditions will be used. This perturbation also produces a perturbation of the total mass, which needs to be corrected. Let j ≥ 2 and Ω j := {(x, v) : |x| < j, |v| > 1/j}. Our first goal is to prove the existence and uniqueness of even functions f j ∈ L 1 Note that, with j → ∞, problem (23) is recovered at least formally, if the mass correction factor in (28) converges to 1, which will be shown below. Problem (25)-(28) will be solved by a fixed point iteration. Let K be the convex set K := ρ ∈ L 1 + (R) : ρ even, ρ(x) = 0 for |x| > j , and let the map T be defined on K by where F is the solution of Since the evenness of ρ ≥ 0 implies evenness of S[ρ] ≥ 0, this problem is invariant under (x, v) ↔ (−x, −v). Therefore it suffices to compute F for x > 0, subject to the symmetry condition and then to extend it by parity. With the explicit representations and the boundary conditions (32), (33), the inflow data at x = 0, v > 0, and at x = j, v < 0, can be computed explicitly: This determines F ≥ 0 completely with the consequence that T (ρ) ≥ 0 is well defined and even, implying that T maps K into itself. Integration of (31), using (30) and (32), implies that In the integral with respect to y, we use S[ρ](y) = 1 2 R e −|y−z| ρ(z)dz ≤ M/2. We obtain which shows that the correction factor in (29) converges to 1 as j → ∞. Before we use that, however, the problem for finite j still has to be solved.
Since ρ → S[ρ] is obviously continuous as a map from L 1 (R) to L ∞ (R), it is straightforward to show continuity of T with respect to the L 1 (R)-topology. We claim it is also compact. Since we have in Ω j , |∂ x ρ F | ≤ jM (ρ + ρ F ) follows, and the estimate (38) implies Therefore the map T : K → K is compact with respect to the L 1 -topology and has a fixed point by the Schauder theorem, i.e. there exists an even solution (f j , ρ j ) of (25)-(28). For the limit j → ∞, we note that the estimate (38) also holds for F = f j with the consequence Since the bound (39) for the derivative is not uniform in j, we shall use as a replacement bounds for the second order moments. From (26) we derive the moment system corresponding to the stationary version of (19): where the notation A j,m,n = R 2 x m v n f j (x, v)d(x, v) Consequently the only possible singular part of f is a Dirac measure at v = 0: with g ∈ L 1 + (R 2 ) and γ ∈ L 1 + (R). In the distributional formulation of (23), we choose a sequence of test functions ϕ n ∈ C ∞ 0 (R 2 ), such that ϕ n and ∂ x ϕ n are uniformly bounded, ϕ n (x, 0) → 1 for all x ∈ R, and |supp(ϕ n )| → 0 as n → ∞. Passing to the limit in the distributional formulation, only the last term remains, showing that γ = 0 and f = g ∈ L 1 (R 2 ).
We observe that (34)-(36) holds for ρ = ρ j and F = f j . The observations of above also justify passing to the limit j → ∞ in these formulas. A short computation then leads to the mild formulation (24), completing the proof. for arbitrary n and, therefore, ρ f ∈ C ∞ (R). Actually, we shall use the equivalent | ρ f (ξ)| ≤ c n (1 + ξ 2 ) n , ∀n ≥ 0 ,