Aggregated Steady States of a Kinetic Model for Chemotaxis Anne Nouri∗ and Christian Schmeiser† 6 1 0 Abstract 2 Akineticchemotaxismodelwithattractiveinteractionbyquasistationarychemicalsignalling n isconsidered. Thespecialchoiceoftheturningoperator,withvelocityjumpsbiasedtowardsthe a chemical concentration gradient, permits closed ODE systems for moments of the distribution J function of arbitrary order. The system for second order moments exhibits a critical mass 8 2 phenomeneon. The mainresultis existence ofanaggregatedsteady state for supercriticalmass. ] P 1 Introduction A . h Chemotaxis, the movement of biological agents influenced by gradients of chemical concentrations, t is a ubiquitous process in biological systems. On the other hand, the production or degradation of a m 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. 1 v A typical example, observed in many bacterial and amoeboid species, is aggregation driven by 2 the production of a diffusible chemical, and chemotactic movement biased towards the direction of 1 the gradient of the chemical concentration (see the large literature on Dictyostelium discoideum or, 7 7 for bacteria, [9]). Since motility usually has a random component, it is a standard question in this 0 situation, if the attractive mechanism is strong enough to overcome the dispersion caused by the . 1 random motility component. 0 The type of mathematical models mostly depends on the nature of the motility process. The 6 1 standard assumption of Brownian motion with drift, the latter determined by chemotaxis, leads v: to a version of the classical Patlak-Keller-Segel (PKS) model [11, 8], where a convection-diffusion i equation for theagent density is coupled with a reaction-diffusion equation for the chemical concen- X tration. For certain bacterial species a description by a velocity jump process is more appropriate, r a whence the convection-diffusion equation of the PKS model is replaced by a kinetic transport equa- tion [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 heatequation. Fordominatingattractiveeffects,theagentdensityeitherhasanontrivialaggregated long-timelimit,oritblowsupinfinitetime,typically inaconcentration event. Thetwo-dimensional parabolic-elliptic PKS model (i.e. with a quasistationary equation for the chemical concentration) hasbeenthoroughly analyzed withrespecttothesequestions. Itshows acritical massphenomenon: Among the initial data with finite variance those with the total mass below a critical value lead ∗Aix-Marseille University, CNRS, Centrale Marseille, I2M UMR 7373, 13453 Marseille, France, [email protected] †Universityof Vienna, Oskar-Morgenstern-Platz 1, 1090 Wien, Austria, [email protected] 1 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 ∂ f +v∂ f = T[S](v′ → v,x,t)f′−T[S](v → v′,x,t)f dv′, (1) t x R Z −D∂2S = βρ −(cid:0)γS, (cid:1) (2) x f a one-dimensional kinetic model for chemotaxis, where the cells with phase space density f(x,v,t) and macroscopic density and flux, ρ (x,t) = f(x,v,t)dv and, respectively, j (x,t) = vf(x,v,t)dv, f f R R Z Z producethe 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: Model A: T[S](v → v′,x,t) = κS(x+αv′,t), Model B: T[S](v → v′,x,t) = κS(x+α(v′ −v),t). 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. Both models have not been derived systematically fromthemicroscopic behaviorof aparticular cell type. However, they arereasonable 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 v x αγ2 γ/D t → αt, v → , x → , f → f , S → S , α γ/D γ/D κβD κ p eliminates all parameters, ip.e., (1), (2) becopmes ∂ f +v∂ f = T[S](v′ → v,x,t)f′−T[S](v → v′,x,t)f dv′, (3) t x R Z −∂2S = ρ −S,(cid:0) (cid:1) (4) x f with Model A: T[S](v → v′,x,t) = S(x+v′,t), Model B: T[S](v → v′,x,t) = S(x+v′−v,t). 2 We consider the Cauchy problem with initial conditions f(x,v,0) = f (x,v) ≥ 0 for x,v ∈ R. (5) I The solution S of (2) is defined as the convolution product of the decaying fundamental solution of −∂2+id with ρ: x 1 S[ρ](x,t) = e−|x−y|ρ(y,t)dy. (6) 2 R Z The initial datum is assumed to possess moments of up to second order: ρ dx = f dvdx = M < ∞, I I R R R Z Z Z |x|2ρ dx = |x|2f dvdx < ∞, I I R R R Z Z Z |v|2f dvdx < ∞. (7) I R R Z Z We choose a reference frame such that the first order moments vanish initially: j dx = vf dvdx = 0, I I R R R Z Z Z xρ dx = xf dvdx = 0. (8) I I R R R Z Z Z 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 plan of the paper is the following. Global in time solutions to the Cauchy problem are determined in Section 2 for models A and B. The long term behavior of their moments is studied in Section 3 and proven to depend on the total mass of the solution. A formal asymptotics of the solution to Model A is performed for large mass in Section 4. The existence of a smooth steady state for Model A with supercritical mass (of the second order moment system) is proven in Section 5. 2 The Cauchy problem Theorem 1 Given f ∈ L1(R2) and M = f (x,v)dxdv, there is a unique solution I + R2 I f ∈ C([0,∞),L1 (R2)) to the Cauchy problem associated to Model A, i.e., + R ∂ f +v∂ f = Q (f), f(t= 0) = f , (9) t x A I with 1 Q (f)(v,x) = S[ρ ](x+v)ρ (x)−Mf(x,v), S[ρ](x) = ρ(y)e−|x−y|dy, (10) A f f 2 R Z where ρ (x,t) = f(x,v,t)dv. f R R 3 Proof. For every T > 0, let X := ρ ∈ C [0,T];L1(R) : ρ(x,t)dx = M, ∀t ∈ [0,T] , T + R (cid:26) Z (cid:27) (cid:0) (cid:1) equipped with the natural norm k·k , and let R (ρ) = f dv, where f is the solution of XT A R ∂ f(x,v,t)+v∂ f(x,v,t) = ρ(x,t)S[ρ](xR+v,t)−Mf(x,v,t), (11) t x f(x,v,0) = f (x,v). (12) I It can be computed explicitly as t R (ρ)(x,t) = e−Mt f (x−vt,v)dv+ e−Ms S[ρ](x+v(1−s),t−s)ρ(x−vs,t−s)dvds. (13) A I ZR Z0 ZR For ρ ∈ X nonnegativity of R (ρ) is obvious, and the mass conservation property follows by T A integration of (13) with respect to x, implying R : X → X . T T The idea is to show that R is a contraction with respect to k · k . First we observe that A XT S[ρ](x,t)dx = M for ρ∈ X and that ρ7→ S[ρ] as a map from L1(R) to L1(R) is Lipschitz with R T Lipschitz constant 1. This implies for ρ ,ρ ∈ X , after a change of variables, 1 2 T R T kR (ρ )−R (ρ )k ≤ e−Ms ρ (ξ,t−s)|S[ρ ]−S[ρ ]|(η,t−s) A 1 A 2 XT 1 1 2 Z0 ZR2 (cid:0)+S[ρ ](η,t−s)|ρ −ρ |(ξ,t−s) dξdηds 2 1 2 ≤ 2(1−e−MT)kρ −ρ k . 1 2 XT (cid:1) Thus, for T < ln2 the map R is a contraction on X . This proves local solvability. As a M A T consequence of the uniform bound in L1(R2) the solution can be extended indefinitely in time steps of length T. Theorem 2 Given f ∈ L1(R2) and M = f (x,v)dxdv, there is a unique solution I + R2 I f ∈ C([0,∞),L1 (R2)) to the Cauchy problem associated to Model B, i.e., + R ∂ f +v∂ f = Q (f), f(t = 0) = f , (14) t x B I with 1 Q (f)(v,x) = S[ρ ](x+v−v′)f(x,v′)dv′ −Mf(x,v), S[ρ](x) = ρ(y)e−|x−y|dy, (15) B f 2 R Z Z where ρ (x,t) = f(x,v,t)dv. f R Proof. A solutioRn f to the Cauchy problem associated to Model B is directly obtained as the limit of the increasing sequence (f ) defined by f = 0 and f given from f as the solution to j 0 j+1 j ∂ f (x,v,t)+v∂ f (x,v,t) = S[ρ ](x+v−v′,t)f (x,v′,t)dv′−Mf (x,v,t), (16) t j+1 x j+1 fj j j+1 Z f (x,v,0) = f (x,v), j+1 I where ρ = f dv. f is explicitly given from f by fj j j+1 j R t f (x,v,t) = e−Mtf (x−vt,v)+ eM(s−t) S[ρ ](x+v(s−t+1)−v′,s)f (x+v(s−t),v′,s)dv′ds. j+1 I fj j Z0 ZR (17) 4 Consequently it can be proven by induction that (f ) is nonnegative, and non decreasing since j f = 0 ≤ f and f ≤ f imply 0 ≤ S[ρ ] ≤ S[ρ ] and f ≤ f by (17). Moreover, denoting 0 1 j−1 j fj−1 fj j j+1 by m (t) = f (x,v,t)dxdv and integrating (16) with respect to (x,v) ∈ R2, it holds j j R m′ = m2−Mm , j+1 j j+1 so that it can be proven by induction that m (t) ≤ M, j ∈ N. (18) j And so, by the monotone convergence theorem, (f ) converges in L1 to a nonnegative function f. It j impliesthatthelimitinL1 of S[ρ ](x+v−v′,t)f (x,v′,t)dv′ is S[ρ ](x+v−v′,t)f(x,v′,t)dv′. fj j f And so, f is a solution of (14). Conservation of mass from (14) implies that f(x,v,t)dxdv = M. (cid:0)R (cid:1) R The limit of (17) implies that f ∈ C([0,∞),L1 (R2)). f is the unique nonnegative solution of (14) + R since its construction makes it minimal among the nonnegative solutions of (14). Indeed, if there wereanother nonnegative solution f˜, thenf ≤ f˜and f(x,v,t)dxdv = f˜(x,v,t)dxdv = M would imply that f˜= f. R R 3 Evolution of moments The closedness of the equations for the moments for Model A relies on the following result Lemma 1 Let 0 ≤ n≤ N, let ρ ∈ L1(R) have finite moments up to order N, i.e., + |x|kρ(x)dx < ∞, k = 0,...,N , R Z and let S be the bounded solution of ∂2S = S −ρ, i.e., x 1 S(x) = e−|x−y|ρ(y)dy. 2 R Z Then also S has finite moments up to order N, and with R := xkρ(x)dx, S := xkS(x)dx, k = 0,...,N , k k R R Z Z the following relations hold: a) S = R +k(k−1)S , k = 0,...,N , k k k−2 n n b) xN−nvnS(x+v)ρ(x)dxdv = (−1)n−kS R 0≤ n ≤ N. k N−k R R k Z Z k=0(cid:18) (cid:19) X Proof. The result can be shown by straightforward computations. We multiply the differential equation for S by |x|k and xk, and use |x|k∂2Sdx = k(k−1) |x|k−2Sdx, and xk∂2Sdx = k(k−1)S , x x k−2 R R R Z Z Z 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. 5 If we concentrate on moments of order N, then the lemma implies S = R +LOT , N N and xN−nvnS(x+v)ρ(x)dxdv = ((−1)n +δ )R R +LOT , n,N 0 N R R Z Z 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: A (t):= xmvnf(x,v,t)dvdx. m,n R R Z Z With the help of Lemma 1 and with − xN−nvn(v∂ f)dvdx = (N −n)A , x N−n−1,n+1 R R Z Z wecanderivedifferentialequationsforthemoments. Thefirstandobviousoneismassconservation: A˙ = 0 =⇒ A =M . 0,0 0,0 As a consequence, the turning operator of Model A after elimination of the unknown S by (6), can now be written as Q (f)(x,v) = S[ρ ](x+v)ρ(x)−Mf(x,v). A f For the first order moments, we obtain A˙ = A , A˙ = −MA =⇒ A = A = 0, 1,0 0,1 0,1 0,1 1,0 0,1 because of (8). For the moments of order two it gets more interesting: A˙ = 2A , 2,0 1,1 A˙ = A −MA −MA , 1,1 0,2 1,1 2,0 A˙ = 2MA −MA +2M2. 0,2 2,0 0,2 Application oftheRouth-Hurwitzcriterion tothecharacteristic polynomialof thecoefficient 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 2M 2M2 (A ,A ,A ) = ,0, . 2,0 1,1 0,2 M −2 M −2 (cid:18) (cid:19) For the higher order moments we only concentrate on the highest order terms on the right hand sides: A˙ = (N −n)A +((−1)n +δ )MA −MA +LOT , N−n,n N−n−1,n+1 n,N N,0 N−n,n 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. 6 If the coefficient matrices of all systems up to order N only have eigenvalues with negative real parts, then all moments of order upto N have finitelimits 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 N−1(−M −λ)n p (λ) = −λ(−M −λ)N +MN! +(−1)NMN!, N n! n=0 X and the determinant of the coefficient matrix is thus given by N−1 Mn p (0) =(−1)NMN!q (M), q (M) = 1+ (−1)N−n . N N N n! n=0 X As we know already and can also be seen from q (M) = 2−M, the 2nd-order coefficient matrix has 2 a zero eigenvalue for M = 2. The same is true for 3rd order (q (M) = M −M2/2 = (2−M)M/2) 3 but, surprisingly, not for 4th order. The function q (M) = 2−M+M2/2−M3/6 has a unique zero 4 M ∈ (2,3). 4 Conjecture: The functions q have unique positive zeroes M , building an increasing sequence, N N which tends to infinity. The Nth-order linear system above is stable, iff M >M . N If the conjecture is true then, for every fixed M > 0, only a finite number of moments tends to a bounded value as t → ∞. This would indicate an M-dependent decay of the equilibrium distribution with stronger decay for larger values of M. The essential parts of the conjecture can be proved: Lemma 2 For fixed N, and M large enough, all roots of p have negative real parts. N Proof. First we look for eigenvalues, which remain bounded as M → ∞. For fixed λ, p (λ) N = (−1)N−1(λ+N)+O(M−1), MN which provides a first root λ = −N +O(M−1). 0 Next we look for roots λ = −M −µ with µ bounded as M → ∞. It is straightforward to show p (−M −µ) N µn N = r (µ)+O(M−1), with r (µ) = (−1)N + . N N N!M n! n=0 X Denoting the roots of r by µ ,...,µ ∈C (multiple roots allowed), we found N more roots of p : N 1 N N λ = −M −µ +O(M−1), j = 1,...,N . j j Obviously, all the N +1 roots we found have negative real parts for large enough M. Lemma 3 For fixed M, and N large enough, there exists at least one positive root of p . N 7 Proof. For fixed λ and M, p (λ) N ≈ e−M−λ+(−1)N , as N → ∞. N!M Therefore, for N large enough there exists λ > 0 such that sign(p (λ)) = (−1)N. On the other N hand, lim p (λ)(−1)N−1 = ∞, N λ→∞ 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 ∈ L1(R2) for an N ≥ 1. I Then the solution f of (9), (10) satisfies (1+|x|N +|v|N)f ∈ L∞(R+; L1(R2)). loc If N ≥ 2 and M > 2, then (1+|x|2+|v|2)f ∈ L∞(R+; L1(R2)). For Model B, the computations are similar but a little more involved. As for Model A, A = M 0,0 and A = A = 0 hold. The 2nd order moments satisfy the closed ODE system 1,0 0,1 A˙ = 2A , 2,0 1,1 A˙ = A −MA , 1,1 0,2 2,0 A˙ = 2M (M +A −A ) . (19) 0,2 2,0 1,1 By their definition and by the Cauchy-Schwarz inequality, for nonvanishing initial data f their I initialvaluessatisfyA (0)A (0) > A (0)2 andA (0),A > 0. Astraightforwardcomputation 2,0 0,2 1,1 2,0 0,2 gives d (A A −A2 ) = 2MA (M +A ). dt 2,0 0,2 1,1 2,0 2,0 This guarantees that A and A remain positive for all times. 2,0 0,2 It is also easily seen that the Jacobian of the right hand side of (19) has at least one positive eigenvalue. TheonlysteadystatehasnegativeA -andA -componentsandcanthereforeneverbe 2,0 0,2 reached. Thus,allsolutioncomponentstendtoinfinityexponentially, meaningthatthechemotactic effect is not strong enough to prevent dispersion. For this reason we concentrate on Model A for the rest of this work. 4 Formal asymptotics for large mass With the rescaling f → Mf, S → MS, Model A takes the form ∂ f +v∂ f = M(S[ρ ](x+v)ρ−f) , (20) t x f with M now taking the role of an inverse Knudsen number. The rescaled version of the steady states for the moments are 2 2M A = , A = 0, A = , (21) 2,0,∞ 1,1,∞ 0,2,∞ M −2 M −2 8 which suggests an equilibrium state concentrating with respect to x as M → ∞. As M → ∞, formally f(x,v,t) → f (x,v,t) = ρ (x,t)S[ρ ](x+v,t). Mass conservation gives 0 0 0 ∂ ρ −∂ (xρ ) = 0. t 0 x 0 Obviously, we have ρ (x,t) → δ(x) as t → ∞ and, thus, 0 1 lim f (x,v,t) = e−|v|δ(x). 0 t→∞ 2 This is in agreement with the limit as M → ∞ in (21). 5 Stationary solutions In this section we first prove in Theorem 3 the existence of even nonnegative L1solutions to the stationary problem, then their C∞ regularity in Theorem 4. Theorem 3 For any M >2 there is an even nonnegative L1(R2) solution f of v∂ f(x,v) = ρ (x)S[ρ ](x+v)−Mf(x,v), f(x,v)dxdv = M. (22) x f f R2 Z Proof. As a consequence of Lemma 1, it holds that S[ρ](x+v)dv = ρ(y)dy, vS[ρ](x+v)dv = ρ(y)(y−x)dy, (23) R R R R Z Z Z Z v2S[ρ](x+v)dv = ρ(y)((y−x)2+2)dy. (24) R R Z Z Let j ∈ N∗ and M = M(1−e−2j −1/j). Our first goal is to prove the existence and uniqueness of j an even function f ∈ L1(R2), such that j + 1 f (x,v) = 0, for |x| > j, or |v| < , or |v| > 4j, (25) j j 1 v∂ f (x,v) = ρ (x)S[ρ ](x+v)−Mf (x,v), |x| <j, < |v| < 4j, (26) x j fj fj j j f (−j,v) =f (−j,−v), f (j,v) = f (j,−v), (27) j j j j f (x,v)dxdv ∈[M ,M]. (28) j j R2 Z Let K be the convex set K := ρ ∈ L1(R): ρ even, ρ(x)= 0 if |x|> j, ρ(x)dx ∈ [M ,M] . + j R (cid:26) Z (cid:27) Let the map T be defined on K by T(ρ) = F dv, where the restriction of F to R+ ×R is the R solution of R 1 F(x,v) = 0, for x > j, or |v| < , or |v| > 4j, (29) j 1 v∂ F(x,v) = ρ(x)S[ρ](x+v)−MF(x,v), 0 < x < j, < |v| < 4j, (30) x j 9 0 F(0,v) = eMj/v F(j,−v)− eMτρ(j +τv)S[ρ](j +τv+v)dτ , v > 0, (31) Z−j/v (cid:0) (cid:1) and F is extended by parity w.r.t. (x,v) to R−×R. Equations (30), (31) imply that F(j,v) = F(j,−v). (32) Denote by X the space of nonnegative L1 −4j,−1 functions with the weight |v|. j j The solution F ∈ L1(R+ ×R) of (29)–(31)(cid:16)e(cid:16)xists and(cid:17)i(cid:17)s unique, because the map + 1 γ ∈X → F(j,−v), −j < v < − , j j where F is the solution of (29), (30) and 0 F(j,v) = γ(v), v < 0, F(0,v) = eMj/v γ(−v)− eMτρ(j +τv)S[ρ|(j +τv+v)dτ , v > 0, Z−j/v ! is a contraction. Indeed, for any (γ1,γ2)∈ Xj2 with images (F1(j,−v),F2(j,−v))−4j<v<−1, we have j −1/j −1/j |v||F (j,−v)−F (j,−v)|dv = |v||γ (v)−γ (v)|e2Mj/vdv 1 2 1 2 Z−j Z−4j −1/j ≤e−M/2 |v||γ (v)−γ (v)|dv. 1 2 Z−j Due to the exponential form of (30), F is nonnegative. Hence, T(ρ) is nonnegative. Moreover, T(ρ) is even since F is even. Equation (30) holds on (−j,j)× (−4j,−1/j)∪(1/j,4j) since F, ρ and S are even functions. Integrating it on (−j,j)× (−4j,−1/j)∪(1/j,4j) and using (29) and (23), ρ (cid:0) (cid:1) implies that (cid:0) (cid:1) M F dxdv = M2− ρ(x) S[ρ](x+v)dvdx. ZR2 ZR Z|v|∈(0,1/j)∪(4j,∞) Moreover, j S[ρ](x+v)dv ≤ e−2j +1/j ρ(y)dy, |x|< j. Z|v|∈(0,1/j)∪(4j,∞) Z−j (cid:0) (cid:1) And so, T maps K into K. We claim that T is compact with respect to the L1 topology. Indeed, let (ρ ) be a sequence in K, i.e., such that n ρ (x)dx ∈ [M ,M], n ∈ N. n j R Z By definition of K, the sequence (T(ρ )) = F dv satisfies n R n (cid:0)R (cid:1) T(ρ )(x)dx ∈ [M ,M]. n j R Z Moreover, 1 |∂ F (x,v)| < jM ρ (x)+F (x,v) , |x| < j, < |v| < 4j, x n n n j (cid:0) (cid:1) 10