Multiple asymptotics of Kinetic Equations with Internal States

TL;DR

This paper investigates the derivation of macroscopic models from microscopic assumptions in bacterial movement, revealing how different asymptotics lead to classical or flux-limited Keller-Segel equations and introducing new equilibrium equations with singular solutions.

Contribution

It provides a mathematical framework connecting microscopic internal states to macroscopic chemotaxis models, highlighting the impact of response stiffness and adaptation time on resulting equations.

Findings

Derivation of macroscopic equations from microscopic internal states.

Identification of conditions leading to Keller-Segel and flux-limited Keller-Segel models.

Discovery of new equilibrium equations with singular solutions.

Abstract

The run and tumble process is well established in order to describe the movement of bacteria in response to a chemical stimulus. However the relation between the tumbling rate and the internal state of bacteria is poorly understood. The present study aims at deriving models at the macroscopic scale from assumptions on the microscopic scales. In particular we are interested in comparisons between the stiffness of the response and the adaptation time. Depending on the asymptotics chosen both the standard Keller-Segel equation and the flux-limited Keller-Segel (FLKS) equation can appear. An interesting mathematical issue arises with a new type of equilibrium equation leading to solution with singularities.

Figures

Tables

Table 1. Table 1. The numerical accuracy for c 0 , 2 subscript 𝑐 0 2 c_{0,2} calculated by equation ( 4.4 ) and < Q > i subscript expectation 𝑄 𝑖 <Q>_{i} for different values of Λ ¯ 0 = 2.0 subscript ¯ Λ 0 2.0 \overline{\Lambda}_{0}=2.0 and 0.5, where the tumbling rate Λ 0 subscript Λ 0 \Lambda_{0} is considered, i.e., χ = 0 𝜒 0 \chi=0 in equation ( 4.18 ). The number of mesh intervals Δ ​ v Δ 𝑣 \Delta v and the time step size Δ ​ t Δ 𝑡 \Delta t are set as J = I 𝐽 𝐼 J=I and Δ ​ t = Δ ​ y / 2 ​ G Δ 𝑡 Δ 𝑦 2 𝐺 \Delta t=\Delta y/2G , respectively. The parameter of the external chemical gradient G = 1.0 𝐺 1.0 G=1.0 is fixed. The reference quantities c 0 , 2 ∗ superscript subscript 𝑐 0 2 c_{0,2}^{*} and < Q > i ∗ superscript subscript expectation 𝑄 𝑖 <Q>_{i}^{*} are the results obtained with I = 3200 𝐼 3200 I=3200 .

	${\overline{\Lambda }}_0$=2.0			${\overline{\Lambda }}_0$=0.5
$I$	$\frac{\left|c_{0,2}^I-c_{0,2}^{\ast }\right|}{c_{0,2}^{\ast }}$			$\frac{\left|c_{0,2}^I-c_{0,2}^{\ast }\right|}{c_{0,2}^{\ast }}$
200	9.1$\times {10}^{-5}$	1.2$\times {10}^{-4}$		8.1$\times {10}^{-4}$	4.2$\times {10}^{-2}$
400	2.2$\times {10}^{-5}$	2.9$\times {10}^{-5}$		5.8$\times {10}^{-4}$	2.5$\times {10}^{-2}$
800	5.4$\times {10}^{-6}$	6.9$\times {10}^{-6}$		3.4$\times {10}^{-4}$	1.3$\times {10}^{-2}$
1600	1.1$\times {10}^{-6}$	1.4$\times {10}^{-6}$		1.5$\times {10}^{-4}$	5.5$\times {10}^{-3}$

Equations

$$\frac{\mathrm{d}m}{\mathrm{d}t}=F(m,M(t,x)).$$

$${\partial }_{t}p+v\cdot {\nabla }_{x}p+{\partial }_m[F(m,M)p]$$

$$Q[m,M](p)=\frac{1}{||\mathbb{V}||}{\int }_{\mathbb{V}}\left[\lambda (m,M,v,v^{\prime })p(t,x,v^{\prime },m)-\lambda (m,M,v^{\prime },v)p(t,x,v,m)\right]\mathrm{d}{v}^{\prime }.$$

$${\partial }_t\rho -{\nabla }_x\cdot \left(D{\nabla }_x\rho \right)+{\nabla }_x\left(\rho \varphi (|{\nabla }_{x}S|)\nabla S\right)=0,$$

$${\partial }_{t}S-D_S\Delta S+\alpha S=\rho ,$$

$$\frac{{\rm d}m}{\,{\rm d}t}=F(m,M(t,x))=\frac{1}{\tau}\big{(}M(t,x)-m\big{)},$$

$$y=M(t,x)-m,$$

$${\partial }_{t}f+v\cdot {\nabla }_{x}f+({\partial }_{t}M+v\cdot {\nabla }_{x}M){\partial }_{y}f-\frac{1}{\tau }{\partial }_y(yf)$$

$$G={\nabla }_{x}M.$$

$$⟨f⟩=\frac{1}{||\mathbb{V}||}{\int }_{\mathbb{V}}f(t,x,v^{\prime },y)d{v}^{\prime },\lambda (y)={\lambda }_0\Lambda \left(\frac{y}{\delta }\right).$$

$${\partial }_{t}f+v\cdot {\nabla }_{x}f+{\partial }_y\left(\left(v\cdot G-\frac{y}{\tau }\right)f\right)$$

$$\tilde{x}=x/L_0,\tilde{v}=v/V_0,\tilde{t}=t/t_0,\tilde{y}=y/\delta ,$$

$$\sigma {\partial }_{\tilde{t}}f+\tilde{v}\cdot {\nabla }_{\tilde{x}}f+{\partial }_{\tilde{y}}\left(\left(\tilde{v}\cdot \frac{\tilde{G}}{\delta }-\frac{\tilde{y}}{\tilde{\tau }}\right)f\right)$$

$${\tilde{\lambda }}_0={\lambda }_0{L}_0/V_0,\tilde{\tau }=\tau /(L_0/V_0),\tilde{G}=G{L}_0,\sigma =L_0/(t_0{V}_0).$$

$$\sigma {\partial }_{t}f+v\cdot {\nabla }_{x}f+{\partial }_y\left(\left(v\cdot \frac{G}{\delta }-\frac{y}{\tau }\right)f\right)$$

$${\lambda }_0=\frac{1}{ϵ}\text{with}ϵ\ll 1.$$

$$\delta =\tau =ϵ={\lambda }_0^{-1},\sigma ={ϵ}^{\alpha -1},\alpha =1\text{or}2.$$

$$\Lambda (y)={\Lambda }_0(y)+ϵ{\Lambda }_1(y)>0,{\Lambda }_0,{\Lambda }_1\in C_b(\mathbb{R}),$$

$${ϵ}^{\alpha }{\partial }_t{q}_{ϵ}+ϵv\cdot {\nabla }_x{q}_{ϵ}+{\partial }_y\left((v\cdot G-y)q_{ϵ}\right)=\Lambda (y)(⟨q_{ϵ}⟩-q_{ϵ}),(y,v)\in \mathbb{R}\times \mathbb{V}.$$

$${\partial }_y\left((v\cdot G-y)Q_0\right)={\Lambda }_0(y)(⟨Q_0⟩-Q_0),Q_0\ge 0,{\iint }_{\mathbb{R}\times \mathbb{V}}{Q}_0\mathrm{d}y\mathrm{d}v=1.$$

$$v_0={\int }_{\mathbb{R}}{\int }_{\mathbb{V}}v{Q}_0\mathrm{d}v\mathrm{d}y.$$

$${\partial }_t{\rho }_0+{\nabla }_x\cdot \left(v_0{\rho }_0\right)=0.$$

$${\partial }_y\left((v\cdot G-y)q_0\right)={\Lambda }_0(y)(⟨q_0⟩-q_0).$$

$${\partial }_t{\rho }_{ϵ}+{\nabla }_x\cdot {\int }_{\mathbb{R}}{\int }_{\mathbb{V}}v{q}_{ϵ}\mathrm{d}v\mathrm{d}y=0,$$

$$f_{ϵ}(t,x,y,v)=q_{ϵ}(t,x+v_{0}t,y,v),$$

$${ϵ}^2{\partial }_t{f}_{ϵ}+ϵ(v-v_0)\cdot {\nabla }_x{f}_{ϵ}+{\partial }_y\left((v\cdot G-y)f_{ϵ}\right)=\Lambda (y)(⟨f_{ϵ}⟩-f_{ϵ}).$$

$$T_{0}g={\partial }_y\left((v\cdot G-y)g\right)-{\Lambda }_0(y)(⟨g⟩-g).$$

$$\underset{y\to \pm |G|}{\lim }\frac{g^2}{Q_0}(v,y)=0\text{for all v∈V.}$$

$${\int }_{\mathbb{R}}{\int }_{\mathbb{V}}\frac{g}{Q_0}{T}_{0}g\mathrm{d}v\mathrm{d}y\ge \frac{1}{2}{\int }_{\mathbb{R}}{\int }_{\mathbb{V}}{\Lambda }_0(y){\left(\frac{g}{Q_0}-⟨\frac{g}{Q_0}⟩\right)}^2{Q}_0(v)\mathrm{d}v\mathrm{d}y\ge 0.$$

$${\int }_{\mathbb{R}}{\int }_{\mathbb{V}}\frac{g}{Q_0}{T}_{0}g\mathrm{d}v\mathrm{d}y$$

$$=\frac{1}{2}{\int }_{\mathbb{R}}{\int }_{\mathbb{V}}{\left(\frac{g}{Q_0}\right)}^2{\partial }_y\left((v\cdot G-y)Q_0\right)\mathrm{d}v\mathrm{d}y-{\int }_{\mathbb{R}}{\int }_{\mathbb{V}}\frac{g}{Q_0}{\Lambda }_0(y)(⟨g⟩-g)\mathrm{d}v\mathrm{d}y,$$

Full text

1. Introduction

Chemotaxis is the movement of bacteria in response to a chemical stimulus. Individual bacteria moves by alternating forward-moving runs and reorienting tumbles. The velocity jump model has been proposed to describe the switching between these two states, in which the two processes, run-to-tumble and tumble-to-run, are modelled by two Poisson processes. The frequencies of these two poisson processes are determined by intracellular molecular biochemical pathway. Therefore, to study the chemotaxis behaviour quantitatively, it is crucial to understand the response of bacteria to signal changes and relate this information to the switching frequency. The mechanism has been well understood for Escherichia coli (E. coli) chemotaxis and other bacteria using similar strategies to move have also been observed ([Duffy97, Hazelbauer12]). In what follows we will focus on the behavior of E. Coli.

It is known that E. coli responds to signal changes in two steps: excitation and adaptation. Excitation is when E. coli rapidly changes the tumbling frequency as it detects external signal changes, while the slow adaptation allows the cell to relax back to the basal tumbling frequency (the frequency when the intracellular chemical reactions are at equilibrium) [JOT]. In the simplest description of the biochemical pathways, a single variable $m$ is used to represent the intracellular methylation level. The methylation has an equilibrium level $M(t,x)$ which is a function of extra-cellular chemical concentration. Using $F(m,M)$ as the adaptation rate, the equation for the adaptation process has the form

[TABLE]

The switching frequencies of run-to-tumble and tumble-to-run are determined by both $m$ and $M(t,x)$. The tumbling time can be ignored since it is usually much shorter than the running time, thus one can combine the two successive processes run-to-tumble and tumble-to-run together and assume that once the bacterium stops run it will immediately choose randomly some direction and start running again. Tumbling frequency is used to describe this Poisson process and it depends on the methylation level $m$ and its equilibrium $M(t,x)$.

In order to understand the relation between individual bacteria movement and their population level behaviour, pathway-based kinetic-transport model has been proposed in [ErbanOthmer04, SWOT]. This model governs the evolution of the probability density function $p(t,x,v,m)$ of the bacteria at time $t$, position $x\in{\mathbb{R}}^{d}$, velocity $v\in\mathbb{V}$, and methylation level $m$. The velocity space $\mathbb{V}$ is $[-V_{0},V_{0}]$ for $d=1$ and $V_{0}\mathbb{S}^{d}$ for $d\geq 2$ (the sphere with radius $V_{0}$). The general form of the kinetic equation is

[TABLE]

where the tumbling term $Q[m,M](p)$ satisfies

[TABLE]

Here $\lambda(m,M,v,v^{\prime})$ is the methylation-dependent tumbling frequency from $v^{\prime}$ to $v$ and $||\mathbb{V}||=\int_{\mathbb{V}}dv$.

The pathway-based kinetic equation in (1.1) is at the mesoscopic level. It can bridge the microscopic (individual) and macroscopic (population) level models by using moment closure or asymptotic analysis [ErbanOth2, SWOT, ST, STY, Xue2]. In this paper, we are interested in building such connections of some macroscopic models with (1.1). Particular attention will be paid to the flux-limited Keller-Segel (FLKS) model. The most standard and popular macroscopic population level model describing the dynamics of the bacteria density is the classical Keller-Segel (KS) model [Patlak, KS1, KS2]. There has been extensive mathematical studies of the KS equation [BBW2015]. Analytically it has been discovered that depending on the cell number, solutions of KS model can either undergo smooth dispersion or blow up in finite time [CW08, BCL09]. The reason that blow-up could happen is because the drift velocity is proportional to the gradient of the external chemicals $|\nabla S|$. Therefore it is not bounded when $|\nabla S|\to\infty$. This is in discrepancy with reality since the population level drift velocity is expected to saturate when the chemical gradient is large. Hence in [HP09] a more physically relevant FLKS model is proposed, which has the form

[TABLE]

where $\rho(t,x)$ is the cell density at time $t$ and position $x$, $\phi(|\nabla_{x}S|)|\nabla_{x}S|$ is a bounded function in $|\nabla_{x}S|$. Unlike the classical KS model, the solution to (1.3) exists globally in time [Chertock, BW2017]. Interesting features of the FLKS are specific stiff response induced unstability and existence of traveling waves (stable or unstable), [PeYa2018, CPY2018].

The classical KS model has been recovered from the kinetic equation (1.1) as diffusion limits, as well as the FLKS [PVW2018]. In [ErbanOthmer04, XO], the authors derived the KS equation by incorporating the linear adaptation. Recently, intrinsically nonlinear signaling pathway are considered [STY, Xue2]. The derivations in [ErbanOthmer04, STY, SWOT, XO, Xue2] are formal and are based on moment closure techniques. When the internal state is not far from its expectation, the moment closure method provides the correct behavior of E. coli on the population level. However, the closeness assumption is only valid when the chemical gradient is small so that the internal states of different bacteria can concentrate and yield the KS equation [SWOT, ST, Xue1]. In general this assumption does not always apply for other scales of the chemical gradient and adaptation time. Our goal is to derive different macroscopic models for different chemical gradient and adaptation time by asymptotic analysis.

We assume that the tumbling frequency $\lambda$ depends only on $M(t,x)-m$. This assumption is valid for E. coli chemotaxis [SWOT]. Moreover, we use the following linear model for the adaptation

[TABLE]

where $\tau$ gives the characteristic time scale of the adaptation process [ErbanOthmer04]. Introducing a new variable

[TABLE]

and letting $f(t,x,v,y)=p(t,x,v,m)$, one can rewrite the original model (1.1) as

[TABLE]

In this paper, we only consider the case when $M(t,x)$ is independent of time and $\nabla_{x}M$ is uniform in space. This special case is physically interesting: for E. coli chemotaxis, since $M(t,x)$ is related to the extra-cellular attractant profile by a logarithmic dependency, a uniform $\nabla_{x}M$ corresponds to the exponential environment as in the experiment in [Kalinin]. Let $G$ be the constant vector given by

[TABLE]

Let $\lambda_{0}$ be the characteristic rate of tumbling and $\delta$ the stiffness of the response. We denote

[TABLE]

Then equation (1.6) becomes

[TABLE]

To perform the asymptotic analysis, we introduce the non-dimensional variables

[TABLE]

where $L_{0}$, $V_{0}$, and $t_{0}$ are the characteristic length, speed, and time respectively. Then the equation 1.8 is reformulated as

[TABLE]

with the parameters given by

[TABLE]

Since there will be no confusion, in the remainder of the paper we drop the tilde sign for simplicity. The final form of the kinetic equation under consideration is

[TABLE]

Our goal in this paper is to start from (1.10) and systematically derive macroscopic limits from it: assuming $\lambda_{0}$ is large and $G=O(1)$, we consider the seamless reorientation by runs and tumbles with different scalings of $\tau$ and $\delta$. More specifically, we investigate the following three cases:

• I.

Fast adaptation and stiff response: both $\tau$ and $\delta$ are as small as $\lambda_{0}^{-1}$;

• II.

Very fast adaptation and very stiff response: both $\tau$ and $\delta$ are much smaller than $\lambda_{0}^{-1}$;

• III.

Moderate adaptation and moderate response, i.e. both $\tau$ and $\delta$ are of $O(1)$.

The common theme in these limits is to decide when the classical KS and the FLKS equations will occur. A brief summary of our result is

1. (1)

In Case I, both the hyperbolic model and a FLKS type model can be found. In particular, when the leading order behaviour of the bacteria population satisfies a hyperbolic model, we further consider the motion of its front profile and find a FLKS model in the moving frame. 2. (2)

In Case II, a FLKS type model can be derived. 3. (3)

In Case III the solution will tend to the solution of a classical KS model.

The rest of the paper is organized as follows. In Section 2, we use asymptotic analysis to formally derive the leading-order macroscopic equations from (1.1). Some properties of the leading order distribution are rigorously shown in Section 3. In Section 4, we show numerical results that are consistent with the analytical properties in Section 3. Finally we conclude in Section 5.

2. Formal Asymptotic Limits

In this section we show that both classical and flux-limited Keller-Segel equations can arise in different physical regimes characterized by the stiffness of the chemotactic response and the rate of adaptation. In all the cases considered below we fix $\lambda_{0}$ as

[TABLE]

2.1. Fast adaptation and stiff response

In this case we assume that

[TABLE]

We also assume that the tumbling frequency $\Lambda$ has the specific structure

[TABLE]

where $\Lambda_{0}$ is a strictly positive continuous function and $\Lambda_{1}$ is independent of $\epsilon$. The space $C_{b}({\mathbb{R}})$ is the collection of continuous and bounded functions on ${\mathbb{R}}$.

Consider the scaled equation

[TABLE]

The leading-order term of $q_{\epsilon}$ is of the form $\rho(t,x)Q_{0}(y,v)$ where $Q_{0}$ satisfies

[TABLE]

The existence and uniqueness of such $Q_{0}$ is shown in Theorem 3.1.

Remark 2.1*.*

In Theorem 3.1 we show that $Q_{0}$ is compactly supported on $[-|G|,|G|]\times\mathbb{V}$ and is strictly positive in $(-|G|,|G|)\times\mathbb{V}$. In the rest of the current section, all the integration involving $Q_{0}$ are performed over $(-|G|,|G|)\times\mathbb{V}$, although we often write the integration domain as ${\mathbb{R}}\times\mathbb{V}$ for the sake of simplicity in notation.

Let $v_{0}$ be the leading-order average velocity given by

[TABLE]

In the case where $v_{0}\neq 0$, by letting $\alpha=1$ one can derive that the leading-order approximate of $q_{\epsilon}$ satisfies a pure transport or hyperbolic equation. More specifically, we have the formal limit

Theorem 2.1**.**

Suppose $\alpha=1$ and $\Lambda$ satisfies (2.1) together with the assumptions in Theorem 3.1. Suppose the average velocity $v_{0}$ defined in (2.4) is nonzero. Let $q_{\epsilon}$ be the non-negative solution to (2.2). Suppose $q_{\epsilon}\to q_{0}$ in $L^{1}(\!\,{\rm d}v\,{\rm d}y\,{\rm d}x)$. Then $q_{0}=\rho_{0}(t,x)Q_{0}(y,v)$ and $\rho_{0}$ satisfies the hyperbolic equation

[TABLE]

Proof.

Since (2.2) is a linear equation, we can pass $\epsilon\to 0$ and obtain $q_{0}$ as a non-negative $L^{1}$ solution to the steady-state equation such that

[TABLE]

By the uniqueness shown in Theorem 3.1, $q_{0}$ is a multiple of $Q_{0}$ with a coefficient $\rho_{0}\in L^{1}(\!\,{\rm d}x)$. The conservation law associated with (2.2) and $\alpha=1$ has the form

[TABLE]

The limiting equation (2.5) is obtained by passing $\epsilon\to 0$ and using that $q_{0}=\rho_{0}Q_{0}$. ∎

To observe nontrivial diffusive behaviour from equation (2.2), we make a change of variable such that

[TABLE]

where again $v_{0}$ is the average speed defined in (2.4). Physically speaking, we consider the motion of the front profile of $q_{\epsilon}$. Now we let $\alpha=2$. Then the equation for $f_{\epsilon}$ reads

[TABLE]

Before showing the formal limit for $f_{\epsilon}$, we prove a lemma of the classical entropy-estimate type:

Lemma 2.1**.**

Let $T_{0}$ be the operator defined as

[TABLE]

Let $Q_{0}$ be the unique solution in (2.3). Then for any $g$ that makes sense of all the integrals and satisfies that

[TABLE]

we have

[TABLE]

Note that a particular case for (2.8) to hold is when $g/Q_{0}\in L^{\infty}({\mathbb{R}}\times\mathbb{V})$.

Proof.

The proof follows from a direct calculation. Indeed, by integration by parts,

[TABLE]

where condition (2.8) is used to guarantee that the boundary terms vanish in the integration by parts. Using the equation for $Q_{0}$, we have

[TABLE]

The equal sign holds only when $g/Q_{0}$ is independent of $v$. ∎

Equipped with Lemma 2.1, we show the formal asymptotic limit of (2.6) in the following theorem:

Theorem 2.2**.**

Let $f_{\epsilon}$ be the non-negative solution to (2.6) and $\rho_{\epsilon}=\int_{\mathbb{R}}\int_{\mathbb{V}}f_{\epsilon}(t,x,y,v)\,{\rm d}v\,{\rm d}y$ such that

[TABLE]

for some $f_{0}(t)\in L^{1}(\,{\rm d}v\,{\rm d}y\,{\rm d}x)$. Then $f_{0}=\rho_{0}(t,x)Q_{0}(y,v)$ and $\rho_{0}$ satisfies the drift-diffusion equation

[TABLE]

where the diffusion coefficient and drift velocity are given by

[TABLE]

*Moreover, $D_{0,2}$ is strictly positive. *

We conjecture that $c_{0,2}$ is bounded in $G$, thus rising again a FLKS type equation. A possible route toward a proof is as follows. We write

[TABLE]

and it remains to establish that $T^{-1\ast}_{0}(v-v_{0})$ is uniformly bounded in $G$ in $L^{\infty}$. Numerics sustain this boundedness, see Fig. 5.

Proof.

First, we use the same argument as in Theorem 2.1 to deduce that there exists $\rho_{0}\in L^{1}(\!\,{\rm d}x)$ such that $f_{0}=\rho_{0}Q_{0}$ and $\rho_{\epsilon}\to\rho_{0}$ in $L^{1}(\!\,{\rm d}x)$. Next, the conservation law for (2.6) reads

[TABLE]

By the definition of $v_{0}$ and $Q_{0}$, we have

[TABLE]

Hence  (2.12) can also be written as

[TABLE]

Then we can re-write equation (2.6) as

[TABLE]

By Theorem A.1, we apply the pseudo-inverse of $T_{0}$ and obtain that

[TABLE]

Note that every term on the right-hand side of the above equality is well-defined since the terms inside $T_{0}^{-1}$ all satisfy the orthogonality condition given in (A.2). By the limit $f_{\epsilon}-\rho_{\epsilon}Q_{0}\to 0$ we have

[TABLE]

Hence, it formally holds that

[TABLE]

Inserting such limit in the limiting form of the conservation law in  (2.12) gives (1.3) with $D_{0,2}$ and $c_{0,2}$ satisfying (2.10) and (2.11) respectively.

Finally, we show that $D_{0,2}$ is a positive definite matrix. To this end, denote

[TABLE]

or equivalently,

[TABLE]

Then for any arbitrary $\alpha\in{\mathbb{R}}^{d}$,

[TABLE]

Hence,

[TABLE]

By Lemma 2.1, if we multiply $(h\cdot\alpha)/Q_{0}$ to $T_{0}(h\cdot\alpha)$ and integrate in $(y,v)$, then

[TABLE]

This shows $\alpha^{T}D_{0,2}\alpha$ is non-negative and the equal sign of the above inequality holds only when $(h\cdot\alpha)/Q_{0}$ is independent of $v$. Now we show that for any $\alpha$ independent of $(y,v)$, the quantity $(h\cdot\alpha)/Q_{0}$ must depend on $v$. Hence the strict positivity of $\alpha^{T}D_{0,2}\alpha$ must hold. To this end, suppose on the contrary that there exists $\alpha_{0}\in{\mathbb{R}}^{d}$ independent of $(y,v)$ and $\beta(y)$ such that $h\cdot\alpha_{0}=\beta(y)Q_{0}$. Then

[TABLE]

which, by the definition of $Q_{0}$ simplifies to

[TABLE]

Since $\alpha_{0}$ is independent of $y,v$, the only solution to such equation is $\beta^{\prime}=0$ and $\alpha_{0}=0$. Hence $\alpha^{T}D_{0,2}\alpha>0$ for any $\alpha\neq 0$. This proves the strict positivity of the diffusive coefficient matrix $D_{0,2}$. ∎

In order to justify that (1.3) is indeed a flux-limited Keller-Segel equation, ideally we want to show that the drift velocity $c_{0,2}$ is generally nonzero, it is uniformly bounded in the chemical gradient $G$ and vanishes when $G$ approaches zero. In what follows we show that these properties are satisfied in the special case where $\Lambda_{0}$ is a constant. In Section 4 we give some numerical evidence that $c_{0,2}$ is uniformly bounded in $G$ even when $\Lambda_{0}$ depends on $y$.

Lemma 2.2**.**

Suppose $\Lambda_{0}$ is a constant. Then the drift velocity $c_{0,2}$ satisfy that

(a) $c_{0,2}$ is nonzero.

(a) $c_{0,2}$ is bounded for all $G$. Moreover, $c_{0,2}\to 0$ as $G\to 0$.

Proof.

(a) Since $\Lambda_{0}$ is a constant, $Q_{0}$ must satisfy the symmetry condition (3.4) in Theorem 3.1. Hence $v_{0}=0$. Moreover, we have the relation $\left(T^{-1}_{0}\right)^{\ast}v=\frac{1}{\Lambda_{0}}v$ since $v$ satisfies the equation

[TABLE]

The drift velocity $c_{0,2}$ is now reduced to

[TABLE]

We want to show that $\int_{\mathbb{V}}vQ_{0}(y,v)\,{\rm d}v$ as a function in $y$ is not zero. To this end, we integrate (2.3) in $v$. This gives

[TABLE]

Since $Q_{0}$ is compactly supported on $[-G,G]\times\mathbb{V}$, one must have

[TABLE]

Since $\int_{\mathbb{V}}Q_{0}\,{\rm d}v$ is not a zero function in $y$, we deduce that $\int_{\mathbb{V}}vQ_{0}(y,v)\,{\rm d}v$ is not a zero function in $y$ either. Since $\int_{\mathbb{V}}vQ_{0}\,{\rm d}v\in L^{1}(\!\,{\rm d}y)$, there exists $\Lambda_{1}\in L^{\infty}(\!\,{\rm d}y)$ such that

[TABLE]

(b) With the assumption that $\Lambda_{0}$ is a constant, we have derived the simplified form of $c_{0,2}$ in (2.15). Then

[TABLE]

Hence $c_{0,2}$ is uniformly bounded in $G$. Since $\Lambda_{1}$ is continuous and $supp\,Q_{0}=[-|G|,|G|]\times\mathbb{V}$, we have

[TABLE]

Hence $c_{0,2}$ vanishes as $G$ approaches zero. ∎

2.2. Very fast adaptation and very stiff response

In this subsection we show a second scaling where flux-limited Keller-Segel equations can also arise. We consider the regime that combines the scalings in [PTV2016] and [PVW2018]. In particular, let

[TABLE]

Then the scaled equation becomes

[TABLE]

Assume that $\Lambda$ satisfies the same condition as in [PVW2018] such that

[TABLE]

where $\overline{\Lambda}_{0}$ is a positive constant and $\Lambda_{1}\in C_{b}({\mathbb{R}})$. Denote ${\mathcal{P}}({\mathbb{R}}^{d}\times{\mathbb{R}}\times\mathbb{V})$ as the space of probability measures in $(x,y,v)$. The formal asymptotic limit is

Theorem 2.3**.**

Let $f_{\epsilon}$ be the solution to (2.16) with $\Lambda$ satisfying (2.17). Suppose there exists $f_{0}(t,\cdot,\cdot,\cdot)\in{\mathcal{P}}({\mathbb{R}}^{d}\times{\mathbb{R}}\times\mathbb{V})$ such that $f_{\epsilon}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}f_{0}$ as probability measures. Then

[TABLE]

Moreover, $\rho_{0}$ satisfies the flux-limited Keller-Segel equation

[TABLE]

with

[TABLE]

Proof.

Passing $\epsilon\to 0$ in (2.16) gives the equation for $f_{0}$ as

[TABLE]

By the non-negativity of $f_{0}$ (and a similar argument as in Theorem 3.1 showing its support), we have

[TABLE]

Hence, there exists $\rho_{0}(t,x)\geq 0$ such that

[TABLE]

The conservation law associated with (2.16) is

[TABLE]

To derive the limit of the $\frac{1}{\epsilon}$-term, we multiply (2.16) by $v$ and integrate in $(y,v)$. This gives

[TABLE]

Passing $\epsilon\to 0$ and applying (2.20), we have

[TABLE]

where $D_{0,3},c_{0,3}$ satisfy (2.19). ∎

Note that in general the drift velocity $c_{0,3}$ is non-zero. If $\mathbb{V}={\mathbb{S}}^{d}$ with $d\geq 1$, then by the rotational invariance, $c_{0,3}$ is along the direction of the chemical gradient $G$. Moreover, it is bounded for all $G$ as long as $\Lambda_{1}$ is a bounded function. Similar as in the proof of Lemma 2.2, for $\Lambda_{1}$ being continuous, we have $c_{0,3}\to 0$ as $G\to 0$. These observations justify that (2.18) is of the flux-limited Keller-Segel type.

2.3. Moderate Adaptation and Moderate Response

The classical Keller-Segel equation can also be derived from kinetic equations with the internal state. The particular scaling is for a slow adaption and a moderate gradient where

[TABLE]

The scaled equation is

[TABLE]

In this case we only need to require that $\Lambda\in C^{1}({\mathbb{R}})$, in particular, it does not have to satisfy the specific form in (2.1). The formal asymptotic limit is

Theorem 2.4**.**

Suppose $\Lambda\in C^{1}({\mathbb{R}})$ and $\Lambda(y)>0$ for all $y\in{\mathbb{R}}$. Let $f_{\epsilon}$ be the non-negative solution to (2.21). Suppose there exists $f_{0}(t,\cdot,\cdot,\cdot)\in{\mathcal{P}}({\mathbb{R}}^{d}\times{\mathbb{R}}\times\mathbb{V})$ such that $f_{\epsilon}\stackrel{{\scriptstyle\ast}}{{\rightharpoonup}}f_{0}$ as probability measures. Then $f_{0}=\rho_{0}(t,x)\delta_{0}(y)$ and $\rho_{0}$ satisfies the Keller-Segel equation

[TABLE]

with the diffusion coefficient and drift velocity given by

[TABLE]

Proof.

The conservation law associated with (2.21) is

[TABLE]

We apply the Hilbert expansion for equation (2.21). Write formally

[TABLE]

Then the leading order $f_{0}$ satisfies

[TABLE]

Hence $f_{0}$ is independent of $v$. The next order in equation (2.21) gives

[TABLE]

The solvability condition requires that

[TABLE]

Since $f_{0}$ is independent of $v$, this reduces to

[TABLE]

By the non-negativity of $f_{0}$, we must have

[TABLE]

Hence, $f_{0}$ concentrates at $y=0$, that is, there exists $\rho_{0}\geq 0$ such that

[TABLE]

Divide (2.21) by $\Lambda$, and then multiply by $v$ and integrate in $(y,v)$. We obtain

[TABLE]

Passing $\epsilon\to 0$, we formally obtain that

[TABLE]

Apply such limit in (2.22) gives the regular Keller-Segel where the drift velocity $c_{0,4}$ is linear in the chemical gradient $G$ as long as $\Lambda^{\prime}(0)\neq 0$. ∎

Remark 2.2*.*

One can consider a fourth case with slow adaptation and stiff response where

[TABLE]

The scaled equation is

[TABLE]

If we change the variable $y$ to $z=\epsilon y$, then the model becomes

[TABLE]

A special choice of $\Lambda=\overline{\Lambda}_{0}+\epsilon\Lambda_{1}$ with $\overline{\Lambda}_{0}$ will give rise to a pure diffusion equation, as can be seen by letting $\Lambda$ to be a constant in Case III. The case of actual interest for further studies is of course when the tumbling rate takes the general form $\Lambda(\frac{z}{\epsilon})$ as in section 2.3.

3. Existence and Properties of $Q_{0}$

In several scaling limits, the equilibrium $Q_{0}$ occurs as the solution of the eigenfunction problem

[TABLE]

Here, we prove existence and uniqueness of $Q_{0}$ assuming that $\Lambda_{0}$ may depend continuously on $y$ and satisfies, for two constants $0<\lambda_{1}\leq\lambda_{2}$, the bounds

[TABLE]

We also derive some basic properties of $Q_{0}$. The precise statement is

Theorem 3.1** (Existence and regularity of $Q_{0}$).**

Suppose $\Lambda$ is continuous and satisfies (3.3), then there exists a unique probability solution $Q_{0}$ of equation (3.1). Moreover, $Q_{0}$, is compactly supported on $[-|G|,\,|G|]\times\mathbb{V}$, strictly positive on $(-|G|,\,|G|)\times\mathbb{V}$ and satisfies

(a) $\left\langle Q_{0}\right\rangle\in L^{\infty}(-|G|,\,|G|)$ and $Q_{0}$ is Lipschitz continuous in compact subsets of ${\mathbb{R}}\times\mathbb{V}\backslash\{v\cdot G=y\}$;

(b) If $\lambda_{2}<1$, then $Q_{0}$ blows up at the diagonal $\{v\cdot G=y,\;|y|<|G|\}$, with a rate no less than $\left\lvert v\cdot G-y\right\rvert^{\lambda_{2}-1}$;

(c) If $\lambda_{1}>1$, then $Q_{0}\in L^{\infty}([-|G|,\,|G|]\times\mathbb{V})$;

(d) If $\lambda_{1}<1$, then $Q_{0}$ blows up at the diagonal at most as

[TABLE]

Furthermore, for $y$ fixed, $Q_{0}$ only depends on $v$ in the direction $G$. In the special case where $\Lambda$ is an even function in $y$, $Q_{0}$ is even in $(y,v)$,

[TABLE]

Remark 3.1*.*

Additionally, one can check that for $\lambda_{2}=1$ the blow-up is at least logarithmic and for $\lambda_{1}=1$, the blow-up is at most logarithmic.

Proof.

The proof is organized as follows. We build a probability solution by a time evolution method and denote integration by $Q_{0}\,{\rm d}v\,{\rm d}y$ rather than ${\rm d}Q_{0}(y,v)$. Then, we prove that such solution satisfies the announced properties. Finally we prove uniqueness.

Compact Support   First we show that any probability solution to (3.1) must be compactly supported on $[-|G|,\,|G|]\times\mathbb{V}$. Indeed, integrating (3.1) in $v$ gives

[TABLE]

Therefore, there exists a constant $\alpha_{0}$ such that

[TABLE]

We claim that $\alpha_{0}=0$ by the following observations:

[TABLE]

For any $y>|G|$, repeating the above argument, we have

[TABLE]

Therefore $\int_{\mathbb{V}}Q_{0}(y,v)\,{\rm d}v=0$ for any $y>|G|$ which implies $Q_{0}(y,v)=0$ if $y>|G|$. Similar argument applies when $y\leq-|G|$. Hence $Q_{0}$ is compactly supported on $[-|G|,\,|G|]\times\mathbb{V}$.

Existence   Second, we show the existence of a solution of (3.1) with a support contained in $[-|G|,\,|G|]\times\mathbb{V}$. To this end, let $\epsilon>0$ be arbitrary and consider the evolution equation

[TABLE]

where the initial data satisfies

[TABLE]

The global existence of a non-negative solution to (3.5) is a classical matter. Moreover,

[TABLE]

We can show that $h_{\epsilon}$ also has a support included in $[-|G|,\,|G|]\times\mathbb{V}$. To this end, let $\phi\in C^{1}({\mathbb{R}})$ be a non-negative function such that

[TABLE]

Because we always have $\phi^{\prime}(y)(v\cdot G-y)\leq 0$, multiplying $\phi(y)$ to (3.5) and integrate in $(y,v)$ gives

[TABLE]

Using the compact support property for $h^{in}(y,v)$, we conclude that

[TABLE]

The compact support property follows.

Now, consider the family of probability measures $\{h_{\epsilon}\}$ on $[0,1]\times{\mathbb{R}}\times\mathbb{V}$. Being compactly supported, it is tight. Therefore, there exists a probability measure $h_{0}(t,x,v)$ such that, after extraction of a subsequence,

[TABLE]

Take the limit $\epsilon\to 0$ in equation (3.5), we have

[TABLE]

Define the probability measure $Q_{0}$ by

[TABLE]

then it satisfies the equation (3.1) with support contained in $[-|G|,\,|G|]\times\mathbb{V}$.

The next steps are to prove that this measure is an $L^{1}$ function with the announced properties.

Some general bounds on $Q_{0}$**   In the sequel, we will make use of the following representation of any weak solution of (3.1). We first deal with the values of $v$ and $y$ such that $y<v\cdot G$. For $\lambda>0$ to be chosen later, we multiply equation (3.1) by $(v\cdot G-y)^{-\lambda}$ and use the chain rule to obtain**

[TABLE]

The choices $\lambda=\lambda_{2}$ and $\lambda=\lambda_{1}$ yield successively

[TABLE]

[TABLE]

Integrating for $z\in(-\infty,y)$ (in fact $(-|G|,y)$ and we keep this bound $-|G|$ which is more convenient for later computations), we find

[TABLE]

[TABLE]

$Q_{0}$ does not carry mass on the diagonal.**   We conclude from (3.9) that $Q_{0}$ is an $L^{\infty}_{loc}$ function away from the diagonal where $v\cdot G=y$. However, this is not enough because our arguments below require that $Q_{0}$ does not carry mass on the diagonal. To do so, decompose $Q_{0}=H_{0}+\mu$ with $\mu$ a measure supported by $\{y=v\cdot G\}$ and $H_{0}\perp\mu$. Inserting this decomposition in the equation, we find**

[TABLE]

Since $\left\langle Q_{0}\right\rangle$ is independent of $v$, its restriction on $y=v\cdot G$ as a measure is zero. Hence $\mu=0$.

Bound $\left\langle Q_{0}\right\rangle\in L^{\infty}(-|G|,\,|G|)$**   We estimate $Q_{0}$ on the sets $v\cdot G-y>0$ and $v\cdot G-y<0$ separately since the diagonal does not carry mass. For $v\cdot G-y>0$, we integrate inequality (3.9) in $v$ and obtain**

[TABLE]

The inner integral in $v$ is bounded as

[TABLE]

and thus , for $z<y$, we conclude, using $c_{d}$ as a constant which may change from line to line and depend on $d,\;G$, $\lambda_{1}$, $\lambda_{2}$, that

[TABLE]

Therefore, we find

[TABLE]

As a first step, we now show that $\left\langle Q_{0}\right\rangle\in L^{2}(-|G|,\,|G|)$. To this end, the contribution to the $L^{2}$-norm of $Q_{0}$ on $\{v\cdot G-y>0\}$ satisfies, using Jensen’s inequality,

[TABLE]

where we have used the fact that $Q_{0}$ is a probability measure and the last inequality holds because

[TABLE]

The estimate for the integral where $v\cdot G-y<0$ is similar. Combining the two parts over $v\cdot G>y$ and $v\cdot G<y$, we conclude that $\left\langle Q_{0}\right\rangle\in L^{2}(-|G|,\,|G|)$.

Building on the $L^{2}$-bound of $\left\langle Q_{0}\right\rangle$ we can now show that $\left\langle Q_{0}\right\rangle\in L^{\infty}(-|G|,\,|G|)$. To this end, use (3.11) and the Cauchy-Schwarz inequality,

[TABLE]

since

[TABLE]

The $L^{\infty}$-bound of $\int_{v\cdot G<y}Q_{0}(y,v)\,{\rm d}v$ follows in a similar way and the details are omitted. Combining these two parts, we obtain that $\left\langle Q_{0}\right\rangle\in L^{\infty}(-|G|,\,|G|)$.

$Q_{0}$ is Lipschitz continuous away from the diagonal.** We now prove continuity of $Q_{0}$ and thus we also obtain that $Q_{0}=0$ for $|y|=|G|$ except for the diagonal points $v\cdot G=\pm|G|$. Such behaviour of $Q_{0}$ is depicted in the numerical results in Section 4. The proof for continuity uses the representation formula for the solution. We re-write the $Q_{0}$-equation as**

[TABLE]

Away from the diagonal $\{v\cdot G=y\}$, the exponential term is smooth in $v$ and continuous in $y$, and the right-hand side is in $L^{\infty}_{loc}(\!\,{\rm d}v\,{\rm d}y)$, which implies that $Q_{0}(y,v)$ is Lipschitz continuous on every compact set in ${\mathbb{R}}\times\mathbb{V}\setminus\{v\cdot G=y\}$.

Strict positivity of $\left\langle Q_{0}\right\rangle$ and $Q_{0}$**   Again we first deal with the values of $v$ and $y$ such that $y<v\cdot G$. Using inequality (3.8) one has**

[TABLE]

**As before we begin with the inner integral and estimate it, for $z\leq y$, as **

[TABLE]

where $c_{d}$ is again a constant depending on $d,\;G$, $\lambda_{1}$, $\lambda_{2}$ which changes from line to line. This yields

[TABLE]

The same calculation, for $y>v\cdot G$ gives,

[TABLE]

Finally we arrive at

[TABLE]

The strict positivity of $Q_{0}$ then follows from the lower bound on $\left\langle Q_{0}\right\rangle$ and (3.8).

Upper bound for $Q_{0}$** Using the upper bound on $\left\langle Q_{0}\right\rangle$ and (3.9), we conclude that for $y<v\cdot G$**

[TABLE]

and a similar inequality for $y>v\cdot G$ concludes the points (c) and (d).

Uniqueness, $v\cdot G$ dependency and symmetry**   Finally, we show that a probability solution to (3.1) must be unique. Suppose instead there are two probability solutions, denoted as $Q_{1}$ and $Q_{2}$. From the steps above, these are $L^{1}$ functions. To avoid the difficulty that $Q_{1}$, $Q_{2}$ vanish near the boundary $y=\pm G$, we introduce**

[TABLE]

The main advantage of $Q_{3}$ is the uniform boundedness given by

[TABLE]

Such bound is elusive a-priori for $Q_{1}/Q_{2}$, which renders integration by parts involving $Q_{1}/Q_{2}$ invalid. Note that by the linearity of (3.1), $Q_{3}$ is also a normalized non-negative solution which satisfies all the properties shown above. In particular, $Q_{3}$ is strictly positive in \big{(}(-|G|,\,|G|)\times\mathbb{V}\big{)}. The uniqueness is shown by a similar argument as proving the positivity of $D_{0,2}$ in Theorem 2.2. In particular, to make use of the entropy-type estimate, we multiply $Q_{1}/Q_{3}$ to the $Q_{1}$-equation and integrate in $(y,v)$. Since $Q_{1},Q_{3}$ are both solutions to (3.1), we apply the same estimate as in the proof of Lemma 2.1 and obtain that

[TABLE]

Therefore, $Q_{1}/Q_{3}$ is independent of $v$. Hence there exists a function $\gamma(y)$ such that

[TABLE]

Inserting such relation in the $Q_{1}$-equation gives

[TABLE]

We can now use the $Q_{3}$-equation to derive that

[TABLE]

This implies that $\gamma$ is a constant function in $y$. If we denote it as $\gamma_{0}$, then $Q_{1}=\gamma_{0}Q_{3}$. By the normalization conditions for both $Q_{1},Q_{3}$, we get $Q_{1}=Q_{3}$, which further implies that $Q_{1}=Q_{2}$, therefore the uniqueness.

From uniqueness, the symmetry in (3.4) follows immediately since $Q_{0}(-y,-v)$ is also a non-negative and normalized solution to (3.1). To show that $Q_{0}$ depends only on $v\cdot G$, we remark that the solution in $d=1$, with $v_{1}$ the direction of $G$ can be extended to a solution in $d$ dimension in $v$ (independent of the orthogonal directions to $v_{1}$) and this provides the unique $d$-dimensional solution. ∎

4. Numerical illustration on $Q_{0}$

We now numerically illustrate the properties established previously on the solution $Q_{0}$ of equation (3.1) with $d=1$ which means $\mathbb{V}=(-1,1)$. We also make the connection with the coefficients found for the continuum FLKS limits and compute the flux $c_{0,2}$ defined by (2.11).

In order to calculate $c_{0,2}$, we first solve the leading order equation  (3.1) to obtain $Q_{0}$, and then we solve the following equation of the next order by using $Q_{0}(y,v)$, i.e.,

[TABLE]

with the constraint

[TABLE]

The flux $c_{0,2}$ is finally obtained by the integration of the solution $h(y,v)$ such that

[TABLE]

We note that in the case when $\Lambda_{0}(y)$ is a constant, i.e., $\Lambda_{0}(y)=\overline{\Lambda}_{0}$, the flux $c_{0,2}$ is directly calculated from $Q_{0}$ by

[TABLE]

4.1. Numerical Scheme

We consider the lattice mesh system on the domain $[-G,G]\times[-1,1]\subset\mathrm{R}^{2}$ such that

[TABLE]

[TABLE]

with $\Delta y=G/I$ and $\Delta v=1/J$.

We choose the mesh system such that for each $y_{i}$, there exists a mesh point on the diagonal $y=vG$. More precisely, $J/I$ is an integer. Denote the mesh point on the diagonal by $(y_{i_{*}},v_{j_{*}})$. For each $v_{j}$, depending on if $y_{i}\leq v_{j}G$ or $y_{i}\geq v_{j}G$, the discretizations at the grid point $(y_{i},v_{j})$ are different. The details are as following.

For fixed $j$, when $y_{i}\leq v_{j}G$, equation (3.1) is discretized by using the first- and second-order backward difference scheme , i.e.,

[TABLE]

and for $y_{i}\geq v_{j}G$, using the first- and second-order forward difference scheme, i.e.,

[TABLE]

with the boundary condition

[TABLE]

Here, we write $Q_{i,j}=Q_{0}(y_{i},v_{j})$, $<Q>_{i}=<Q_{0}>(y_{i})$, $\Lambda_{0}^{i}=\Lambda_{0}(y_{i})$, and

[TABLE]

At each mesh point on the diagonal ($y_{i_{*}}$,$v_{j_{*}}$), we calculate two values for $Q_{i_{*},j_{*}}$ by using equation (4.7) and (4.8) with letting $W_{i_{*},j_{*}}Q_{i_{*},j_{*}}=0$ on each eqaution. The values calculated at the left side of the diagonal by (4.7), say $Q^{L}_{i_{*}}$, are used in the integration $\int_{v_{j_{*}}}^{1}Q_{0}(y_{i},v^{\prime})dv^{\prime}$ and the values calculated at the right side of the diagonal by (4.8), say $Q^{R}_{i_{*}}$, are used in the integration $\int_{-1}^{v_{j_{*}}}Q_{0}(y_{i},v^{\prime})dv^{\prime}$.

Thus, we calculate $<Q>_{i}$ ($i=1,\cdots,2I-1$) as

[TABLE]

and $<Q>_{0}=<Q>_{2I}=0$. Hereafter, the trapezoidal rule is used for the integration with respect to $v$ while Simpson’s rule is used for the integration with respect to $y$. Thus, for example, the integration $\overline{Q}=\int_{-G}^{G}<Q_{0}>(y)dy$ is calculated as

[TABLE]

To obtain the solution of Eqs. (4.7)–(4.9), we consider the following time-evolution semi-implicit scheme,

[TABLE]

for $y_{1}\leq y_{i}\leq v_{j}G$ with a uniform initial condition. This scheme uses a lower diagonal matrix which solves $Q^{n+1}_{i,j}$ very quickly. For $v_{j}G\leq y_{i}\leq y_{2I-1}$, we replace the backward difference in (4.13) with the forward difference (4.8), i.e.,

[TABLE]

which solves $Q^{n+1}_{i,j}$ by using the upper diagonal matrix.

At each time step, we calculate $\overline{Q}^{n}$ by equation (4.12) and normalize the solution of the time-evolution scheme (4.13) and (4.14), say $\widetilde{Q}^{n}_{i,j}$, by

[TABLE]

in order to satisfy the normalized condition $\int_{-G}^{G}\int_{-1}^{1}Q(y,v)dydv=1$.

We repeat the above process until the $Q_{i,j}^{n}$ satisfy the following convergence condition

[TABLE]

After we obtain the numerical solution $Q_{i,j}$, we solve equation (4.1) by using the same time-evolution scheme as (4.13) and (4.14) with the inhomogeneous term $\Lambda_{1}^{i}(<Q>_{i}-Q_{i,j})$. In order to satisfy the condition (4.2), we correct the solution of the time-evolution scheme at each time step, say $\widetilde{h}^{n}_{i,j}$, by

[TABLE]

where $r$ is calculated as $r=\int_{-G}^{G}\int_{-1}^{1}\widetilde{h}^{n}(y,v)dydv$.

4.2. Numerical Result

We carry out the numerical computation of (4.13)–(4.17) when the tumbling functions are written as,

[TABLE]

and

[TABLE]

Here, $\overline{\Lambda}_{0}$ is the mean tumbling rate and $\chi$ is the modulation amplitude of $\Lambda_{0}(y)$.

In the numerical computation, we set the time-step size as $\Delta t=\Delta y/2G$ and the numbers of mesh interval $\Delta y$ and $\Delta v$ as $I=J$. Table 1 shows the accuracy of the numerical results for various mesh systems in the case when $\Lambda_{0}$ is constant, i.e., $\Lambda_{0}=\overline{\Lambda}_{0}$ and $\chi=0$. It is seen that the flux $c_{0,2}$ and the moment $<Q_{0}>$ converge against the mesh interval for both $\overline{\Lambda}_{0}$=2.0 and 0.5. The accuracy of $c_{0,2}$ is approximately second order for $\overline{\Lambda}_{0}$=2.0 while it is approximately first order for $\overline{\Lambda}_{0}$=0.5. We also compare the results of $c_{0,2}$ obtained by (4.3) and by (4.4). The relative differences of the two values are 1.5$\times 10^{-6}$ for $\overline{\Lambda}_{0}=2$ and 2.2$\times 10^{-2}$ for $\overline{\Lambda}_{0}=0.5$ in the fine mesh system with $I=3200$. In the following of this section, we show the numerical results obtained in the mesh system with $I=1600$.

Figures 1 shows the numerical results of $Q_{0}$ for different values of $\overline{\Lambda}_{0}$ when the tumbling rate $\Lambda_{0}(y)$ is constant, i.e, $\chi=0$. It is seen that the distributions are origin symmetry, i.e., $Q(y,v)=Q(-y,-v)$ because the tumbling rate $\Lambda_{0}$ is independent of the internal state $y$. As the tumbling rate $\overline{\Lambda}_{0}$ decreases, $Q_{0}$ becomes more and more concentrated at the diagonal $y=vG$.

Figure 2 shows the one-dimensional distribution of $Q_{0}(y,v)$ against $(vG-y)^{-1}$ in the domain $vG-y>0$ when the velocity $v$ is fixed. It is clearly seen that $Q_{0}(y,v)$ blows up at the diagonal for $\overline{\Lambda}_{0}\leq 1$; the rates of the blowup are $(vG-y)^{\overline{\Lambda}_{0}-1}$ for $\overline{\Lambda}_{0}=0.5$ and logarithmic for $\overline{\Lambda}_{0}=1.0$ as is established in Theorem 3.1. On the other hand, for $\overline{L}ambda_{0}$=2.0, $Q_{0}$ is continuous at the diagonal $y=vG$. This is seen because the curves for $v$=0.75 and 0.5 (dashed lines) correspond to those for $v=-0.75$ and $-0.5$ in $vG-y<0$, respectively, due to the symmetry and the curves with the same values of $|v|$ converge to the same values at the diagonal $y=vG$.

Figure 3 shows the decay behavior of $Q_{0}(y,v)$ near $y=-G$ for a fixed velocity $v$. These plots illustrate the behavior $Q_{0}(y,v)=(|G|-|y|)^{\alpha}$ near $y=-G$ with $\alpha>1$. The power $\alpha$ is measured as $\alpha=1.6$ for $\overline{\Lambda}_{0}=0.5$, $\alpha=1.8$ for $\overline{\Lambda}_{0}=1.0$, and $\alpha=2.6$ for $\overline{\Lambda}_{0}=2.0$. Thus, the power increases with $\overline{\Lambda}_{0}$. All these numerical results illustrate the properties of $Q_{0}(y,v)$ established in Theorem 3.1.

Figure 4 shows the distributions of $Q_{0}(y,v)$ and $h(y,v)$ when the tumbling rate $\Lambda_{0}$ depends on the internal state $y$. It is clearly seen that the distributions are not anymore symmetric but rather concentrated along the diagonal $y=vG$ in $y>0$.

Finally, we show the drift velocity $c_{0,2}$ against the chemical gradient $G$ in Figure 5. Figure 5(a) shows the results for various values of $\overline{\Lambda}_{0}$ when the tumbling rate is constant, i.e., $\chi$=0, while figure 5(b) shows the results for various values of modulation $\chi$, where the tumbling rate $\Lambda_{0}$ is not constant. In both cases, when the gradient is small, say $G\lesssim 1$, the fluxes $c_{0,2}$ are almost linearly proportional to the gradient, $c_{0}\propto G$. However, they saturate for $G\gtrsim 10$ and approach to the constant values. These results illustrate the boundedness of the drift velocity $c_{0,2}$ when the tumbling rate $\Lambda_{0}$ is constant. Unexpectedly, we can also observe non-monotonic profiles of $c_{0,2}$ against the chemical gradient $G$ when the modulation is large, i.e., $\chi$=0.5 and 0.8.

5. Conclusion

We have systematically studied scaling limits of kinetic equations which describe run-and-tumble movement of bacteria with internal chemical pathways. The complexity of the phenomena stems form the different scales between space, time and velocity as well as the stiffness of the tumbling response and methylation adaptation. The question in these different asymptotic limits is to distinguish between the standard Keller-Segel equation and the flux-limitied Keller-Segel.

It appears that stiff-response, a physically appropriate regime, always leads to the FLKS equation. In particular this conclusion applies to fast adaptation and stiff response which corresponds to the measured paprameters for E. coli.

The dominating profile, given by the function $Q_{0}(y,v)$, is also remarkable with a possible blow-up along the diagonal $y=v|G|$. However, the physically relevant regime is when $\Lambda_{0}>1$ and then $Q_{0}(y,v)$ is smooth.

Several other possible scalings are still possible. Also we have ignored noise in the internal state, which can be a route to mathematically interesting analysis.

Appendix A Well-posedness of the corrector equation (2.14)

We show that the equation (2.14) for the corrector has a solution using the Fredholm theory. More generally, we consider the equation with a source term $R_{1}$ such that

[TABLE]

The well-posedness theorem states

Theorem A.1**.**

Let $Q_{0}$ be the unique non-negative solution given in Theorem 3.1. Suppose the function $R_{1}$ satisfied that $R_{1}\in L^{1}((-|G|,\,|G|)\times\mathbb{V})$ and

[TABLE]

and there exist constants $\mu_{0}>-1$ and $c_{0}>0$ such that

[TABLE]

Then equation (A.1) has a unique solution $h\in L^{1}((-|G|,\,|G|)\times\mathbb{V})$ such that

[TABLE]

Proof.

Re-write (A.1) as

[TABLE]

Then for $v\cdot G-y>0$, we have

[TABLE]

which gives

[TABLE]

Similarly, if $v\cdot G-y<0$, then

[TABLE]

which gives

[TABLE]

Using (A.5) and (A.6), we obtain the integral equation for $\left\langle h\right\rangle$ as

[TABLE]

Denote $K$ as the linear operator such that

[TABLE]

Denote the source term as

[TABLE]

Using the bound of $R_{1}$ in (A.3) and the estimate, for $v\cdot G-y>0$,

[TABLE]

we bound the first term in $R$ as

[TABLE]

Similar as the estimates in (3.10), we have

[TABLE]

where again the constant $c_{d}$ may change from line to line. Since $\mu_{0}>-1$, the first term in $R$ is bounded as

[TABLE]

Similar estimates hold for the second term in $R$. Overall we have

[TABLE]

Moreover, $\left\langle h\right\rangle$ satisfies

[TABLE]

We will show that the operator ${\rm{Id}}-K$ is Fredholm and $R\in({\rm Null}\,({\rm{Id}}-K^{\ast}))^{\perp}$ where the orthogonality is taken in $L^{2}(-|G|,|G|)$. To this end, we first find the kernel of $K$ as

[TABLE]

and the kernel of the adjoint operator $K^{\ast}$ is

[TABLE]

In order to show that $K$ is compact on $L^{2}(-|G|,\,|G|)$ we prove that $k\in L^{2}(\,{\rm d}y\,{\rm d}z)$ so that $K$ is a Hilbert-Schmidt operator. We only present the details for the $L^{2}$-bound of the first term in $k(y,v)$, for the bound of the second term follows in a similar way. Hence, by (A.8) and the upper/lower bounds of $\Lambda$, we have

[TABLE]

where the last step follows from (3.10). Hence $K$ is compact on $L^{2}(-|G|,\,|G|)$.

Next, we show that $R\in\left({\rm Null}\,({\rm{Id}}-K^{\ast})\right)^{\perp}$. The exponential term in $k^{\ast}$ can be simplified as

[TABLE]

and similarly,

[TABLE]

Hence, $k^{\ast}$ becomes

[TABLE]

The structure of the null space of ${\rm{Id}}-K^{\ast}$ can be seen by working with a modified kernel

[TABLE]

We claim that for any $y\in(-|G|,|G|)$, $\widetilde{k}^{\ast}$ satisfies the normalization such that

[TABLE]

Indeed, by using the spherical coordinates, the first term in $\int_{-|G|}^{|G|}\widetilde{k}^{\ast}(y,z)\,{\rm d}z$ satisfies

[TABLE]

where $J$ is the Jacobian in the spherical coordinates. The last inner integral is evaluated as

[TABLE]

Hence, the first term in $\int_{-|G|}^{|G|}\widetilde{k}^{\ast}(y,z)\,{\rm d}z$ is reduced to

[TABLE]

Similarly, the second integral in $\int_{-|G|}^{|G|}\widetilde{k}^{\ast}(y,z)\,{\rm d}z$ can be simplified as

[TABLE]

Summing up gives the normalization condition (A.10) for $\widetilde{k}^{\ast}$. Therefore, if we denote $\widetilde{K}^{\ast}$ as the operator with kernel $\widetilde{k}^{\ast}$, then the normalization condition guarantees that the only elements in ${\rm Null}\,({\rm{Id}}-\widetilde{K}^{\ast})$ are constants, which further implies that

[TABLE]

Hence, we have the equivalence such that

[TABLE]

Now we check such orthogonality condition indeed holds. This can be done through a direct computation using the definition of $R$ given in (A). We can also show it by observing that $R=\left\langle w\right\rangle$ where $w$ is the solution to

[TABLE]

Then (A.11) holds by integrating (A.12) in $y$ and using condition (A.2) for $R_{1}$.

Combining the orthogonality and the Fredholm property of ${\rm{Id}}-K$, we deduce that (A.9) has a solution in $L^{2}(-|G|,\,|G|)\cap\left({\rm Null}\,(I-K)\right)^{\perp}$. Denote this solution as $H_{1}$ and define

[TABLE]

Then

[TABLE]

which shows $h_{1}$ is the unique solution to (A.1) satisfying (A.4) with $\left\langle h_{1}\right\rangle\in L^{2}(-|G|,\,|G|)$. ∎

Acknowledgements: BP and WS would like to express their gratitude to the NSF research network grant RNMS11-07444 (KI-Net). WS would also like to thank Laboratoire Jacques-Louis Lions (LJLL) at Sorbonne University for their hospitality and Centre National de la Recherche Scientifique (CNRS) for providing generous support.

BP has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 740623). The research of WS is supported in part by NSERC Discovery Individual Grant R611626. MT is partially supported by NSFC 11871340. SY has received funding from the Japan Society for the Promotion of Science (JSPS) under the KAKENHI grant No 16K17554.