Bacterial chemotaxis considering memory effects: Derivation of the reaction-diffusion equations

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ “This is an Accepted Manuscript of an article published in Physical Review E on 15 May 2025, available at: https://doi.org/10.1103/PhysRevE.111.054409 .” Bacterial chemotaxis considering memory effects Manuel Mayo1, 2 and Rodrigo Soto2 1F´ısica Te´orica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain 2Departamento de F´ısica, Facultad de Ciencias F´ısicas y Matem´aticas, Universidad de Chile, Avenida Blanco Encalada 2008, Santiago, Chile (Dated: April 19, 2025) Bacterial chemotaxis for the case of E. coli is controlled by methylation of chemoreceptors, which in a biochemical pathway regulates the concentration of the CheY-P protein that finally controls the tumbling rate. As a consequence, the tumbling rate adjusts to changes in the concentration of relevant chemicals, such as to produce a biased random walk toward chemoattractants of against the repellers. Methylation is a slow process, implying that the internal concentration of CheY-P is not instantaneously adapted to the environment, and therefore the tumbling rate presents a memory. This implies that the Keller–Segel equations used to describe chemotaxis at the macroscopic scale, which assume a local relation between the bacterial flux and the chemical gradient, cannot be fully valid as memory and the associated nonlocal response are not taken into account. To derive the new equations that replace the Keller–Segel ones, we use a kinetic approach, in which a kinetic equation for the bacterial transport is written considering the dynamics of the protein concentration. When memory is large, the protein concentration field must be considered a relevant variable on equal foot as the bacterial density. Working out in detail the Chapman–Enskog method, the dynamical equations for these fields are obtained, which have the form of reaction-diffusion equations with flux and source terms depending on the gradients on the chemical signal. Also, the transport coefficients are obtained entirely in terms o the microscopic dynamics, showing important symmetry properties and giving their values of the case of E. coli. Solving the equations for an inhomogeneous signal it is shown that the response is nonlocal, with a smoothing length as large as 170 µm for E. coli. The homogeneous response and the relaxational dynamics are also studied in detail. For completeness, the case of small memory is also studied, in which case the Chapman–Enskog method reproduces the Keller–Segel equations, with explicit expressions for the transport coefficients. I. INTRODUCTION The pioneering works of Berg and coworkers gave the study of bacterial transport a quantitative character, with detailed mathematical models [1, 2]. The most relevant transport process is chemotaxis, where bacteria move along or against chemical gradients. The current understanding of chemotaxis of flagellated bacteria like E. coli is that they perform a biased random walk, with longer walks in the direction of the attractant and shorter in the opposite direction. For E. coli, this is achieved by modulating the tumble rate in response to changes in the chemical concentration as they move. As the tumbling process is stochastic [1, 3, 4], this description immediately implies that bacterial diffusion appears together with chemotaxis. At the macroscopic scale, bacterial chemotaxis is described by the Keller–Segel equations, which couple the dynamics of the bacterial density ρto the ligand concentration l(food, chemoattractant, or chemorepellent) [5, 6]. In presence of a ligand gradient, the bacterial flux has a chemotactic term, proportional to the ligand gradient, which is added to the diffusive flux. Hence, for the density one gets ∂ρ ∂t =−∇·J,(1a) with J=−D∇ρ+µρ∇l. (1b) In case cell division or death processes also take place, a source term can be added to the first equation. Here, D is the diffusion coefficient and µthe chemotactic mobility, which is positive for chemoattractants and negative for repellers. Coupled to this, the Keller–Segel system is complemented by an equation for the ligand concentration, which is a diffusion equation with a sink term representing the depletion of the ligand by bacteria. In this article, we will analyze the equation for the bacterial density and, therefore, we will assume that the ligand gradient is imposed. Eq. (1b) uses a linear coupling with the ligand, and more refined relations have been obtained. For example, the flux being a nonlinear function of ∇lor the mobility depending on the absolute value of the ligand [7–11]. The Keller–Segel equations can be derived from the run-and-tumble dynamics applying different techniques of nonequilibrium statistical mechanics [7, 8, 12–16]. For that it is assumed that gradients of the different fields are small compared to the scale of the microscopic agents (bacteria in this case). That is, the Keller–Segel equations result from a process of coarse graining to the socalled hydrodynamic scale. This process does not only generate the dynamical equations, but also the transport coefficients Dand µ. The tumbling rate in bacteria like E. coli is controlled by the concentration of the CheY-P protein inside the bacterial body and chemotaxis results from varying the equilibrium concentration of this protein as response to changes in the ligand concentration [4, 9, 17–19]. The 2 change of Chey-P concentration is a result of a biochemical pathway, where an important element in the methylation of the chemoreceptors. Methylation is a slow process and therefore the concentration of CheY-P presents a long memory time, introducing a new temporal scale of tens of seconds [10, 20–22] and an associated length scale of hundreds of micrometers [23, 24]. These scales are comparable to those appearing in microfluidic devices and in the complex natural environments where bacteria live, such as soils, organs, or pores (see, for example, [25–29]). Also, in sea water, microscale nutrient patches appear in the form of chemical pulses [30, 31]. As a consequence, the simple Keller–Segel equations are not expected to be valid on those scales and new macroscopic equations valid for those scales need to be derived. One possible approach could be to include terms of higher gradient order in the chemotactic flux or the use of nonlocal kernels in the flux. This has been the approach followed by some authors. In Ref. [32] the one-dimensional case is considered, which is extended into two dimensions at the expense of performing some uncontrolled approximations [33]. Similarly, in Ref. [34] the need of introducing additional variables make that the resulting macroscopic equations are not fully based on the microscopic model. In Refs. [35, 36], although starting from a kinetic model with memory the macroscopic equations are obtained in the long time limit, in which case the Keller-Segel equations are recovered. Finally, it must be noted that in those works, the fluctuations of the internal variables are completely neglected, which we show below is not correct, at least for the case of E. coli. Here, we adopt a different strategy to obtain the macroscopic equations. First, we show that the bacterial concentration of CheY-P, ρX, is a relevant field, as important to the density, to describe the chemotactic process. Hence, it is necessary to derive the coupled equations for ρand ρX, which is the main purpose of this article. As different bacteria strains can differ on their tumbling and methylation memory times, the latter not being necessarily large as in the case of E. coli, we will derive the hydrodynamic equations for both for cases of large and small memory. The article is organized as follows. In Sect. II we first present the mathematical model for chemotaxis of bacteria like E. coli at the individual scale, which is incorporated into a kinetic theory for an ensemble of bacteria in presence of chemical signals. We end this section presenting some relevant mathematical properties of the kinetic equation. In Sect. III we apply the ChapmanEnskog procedure, which is a systematic approach to derive coarse-grained equations for successive temporal scales, to the kinetic equation when the memory is not necessarily large. In this case, the usual Keller–Segel model is recovered. The Chapman-Enskog method for the case of large memory is worked out in Sect. IV, where new equations are obtained together with the associated transport coefficients. This results in Eqs. (91-93) that replace the Keller–Segel equations (1) and constitute the principal contribution of this article. Several relevant cases are studied in detail in Sect. V. The short section VI gives the values of the transport coefficients for the case of E. coli. Finally, conclusions are given in Sect. VII. II. KINETIC THEORY OF BACTERIAL CHEMOTAXIS A. Chemotactic model Bacteria like E. coli move in fluids following the socalled run-and-tumble dynamics [1]. In the run phase, bacteria move with roughly constant speed and direction, process that is interrupted by rapid reorientations called tumbles. The latter are initiated when there is a reversal in the direction of rotation (from counterclockwise, CCW, to clockwise, CW) of one or multiple flagella [2, 37]. This process of motor switch is governed by the concentration y(t) of the phosphorylated CheY protein (CheY-P) inside the bacterial body. In a simple description, based on the biochemistry of the molecular motor, Tu and Grinstein proposed a model to describe the tumbling process as a two-state activated system [4]. In this model, the activation free energy barrier ∆Gand, hence, the transition rate ν∼e−∆G/kT from the CCW to the CW state depend on the instantaneous concentration y(t). Considering a Taylor expansion of ∆Gand defining the normalized concentration deviation, X(t) = [y(t)− ⟨y⟩]/σy, it results ν=ν0eλX , where σy is the standard deviation of yand the dimensionless parameter λmeasures the sensitivity of the tumbling rate to changes in X. In the Tu and Grinstein model, Xis a Gaussian variable of null average and unit variance, with correlation time τ. In summary, Xis described therefore as an Ornstein–Uhlenbeck process. By tracking E. coli (RP437 bacteria in motility buffer supplemented with serine), this model was validated and allowed to extract the model parameters: ν0= 0.22 s−1,τ= 19 s, and λ= 1.62 [22]. As advanced, the memory time τ, which is governed by the methylation process, is long compared to the mean run time. Associated to it, there is a characteristic memory length L=V τ, where Vis the bacterial swim speed. Using V= 27 µm/s from Ref. [22], gives L= 500 µm. This relatively large correlation length, comparable to pore sizes in natural environments, generate nonlocal responses that we aim to consider in our description. Also, the large value of λindicates that the fluctuations in Xare important and their effects cannot be disregarded. Bacteria respond to chemotactic signals by modulating their tumbling rate. As these microorganisms are too small to measure gradients along their body, they integrate the ligand signal on their run, with a chemical pathway that can be modeled by three principal components: methylation level m(t), kinase activity a(t), and the aforementioned y(t). With these elements, the tumbling rate modulation can be summarized as follows. 3 Starting from a steady state with a=y= 0 and m adapted to the ligand concentration l, an abrupt change of ltriggers a rapid decrease of yand ato negative values. Later, on a larger methylation time scale τm(of the same order of τ)madapts to a new value, and aand ybecome null again. Hence, for a time τm,ybecomes negative as a response to changes on l. That is, yfollows the negative of the ligand time derivative. The dynamics of this process can be expressed mathematically using coupled Langevin equations for a(t), y(t), and m(t) [9, 17–19]. Theoretical analysis of these equation have studied the effects of memory on the response when time-dependent signals are imposed [19, 38, 39]. In Ref. [24], by adiabatically eliminating the fast modes of these equations, and keeping linear couplings, we showed that the chemotactic coupling can be described by a modification of the Tu and Grinstein dynamical equation for X, to incorporate the coupling with the ligand as ˙ X=−X+b˙ l τ+r2 τξ, (2) where in this Langevin equation ξis a white noise of correlation ⟨ξ(t)ξ(t′)⟩=δ(t−t′), bis the coupling constant to the ligand (positive for attractants and negative for repellers), and the square root prefactor guarantees that if ˙ lis constant, Xreaches a normal distribution of unit variance. This simple model has been improved in several aspects. Motor adaptation implies that the νdependence on Xis not exponential [20, 21, 40, 41]. Also, the fluctuation of CheY-P are not small, implying that the Langevin equation for Xneeds to include non-linear coupling terms [42, 43]. There are several models, still under study, and therefore we consider the most general description for the stochastic dynamics of X. In principle the noise intensity can also depend on X. With a change of variable it is always possible to make this intensity constant, resulting in the general equation ˙ X=−A(X, l) + B(X, l)˙ l τ+r2 τξ. (3) Here, due to the possible change of variable, Xis not exactly the normalized CheY-P concentration, but it is closely related to it. To avoid overloading the language, we will continue calling it the normalized CheY-P concentration. Aand Bare dimensionless functions of order one, and we keep τas the only relevant time scale. We assume that in absence of noise, for any value of l, there is a single stable fixed point. Also, to ensure the existence of a linear response regime, the scalings A(X)∼X and B(X)∼bshould be imposed for small X. And for the tumbling rate we take, ν=ν0C(X, l),(4) where Cis a monotonous increasing function of X, normalized such that C(0) = 1. In what follows, we will work with this general model, but for concrete results we will use the linear model, corresponding to Alin(X) = X, Blin(X) = b, Clin(X) = eλX.(5) Finally, Eq. (3) is coupled to the bacterial motion because ˙ lis the rate of change of lin the comoving frame of the swimmer. For a bacterium moving with speed Valong the director ˆ n, it is written as the Lagrangian derivative ˙ l=Vˆ n·∇l+∂l ∂t.(6) In this representation, chemotaxis is rationalized as follows. For simplicity, we consider the linear model, but the analysis is analogous for the general case. For a swimmer moving parallel to a chemoattractant gradient (ˆ n· ∇l > 0 and b > 0), the normalized protein X becomes negative on average (that is, ybecomes smaller than ⟨y⟩), reducing by Eq. (4) the tumbling rate, and exactly the opposite results for a swimmer moving against the gradient. The result is a biased run-and-tumble random walk, with longer runs in the direction of the gradient. Similarly, for a chemorepellent (b < 0), the motion is biased against the gradient. In absence of memory and fluctuations, Xadapts instantaneously the ligand rate of change, X=−b˙ l, resulting in the tumbling rate ν=ν0e−λb˙ l. This expression or similar ones have been used to describe chemotaxis at the kinetic level [16, 44– 49]. Here, we go beyond this approximation, considering the effects of fluctuations and memory in the chemotactic process. In Ref. [24] we proposed a kinetic equation that incorporated all the elements of the linear model, which we solved to study the stationary chemotactic mobility and the linear response to signals varying in space and time. Importantly, the obtained response is nonlocal in space and time, as an effect of the memory time. The method to build the solutions is not easy to adapt to different geometries and configurations, as it needs to fully describe the distribution function f(r,ˆ n, X, t) of having bacteria at position r, with director ˆ n, and normalized concentration of the CheY-P protein Xat time t. A simpler approach is to study the dynamics of slowly-varying fields analog to the Keller–Segel equations (1) for the density field. With that purpose, we will first derive the kinetic equation that describes the dynamics of bacterial suspensions coupled with chemotactic signals with the general model (3) and (4). Then, we will apply the Chapman– Enskog procedure to this kinetic equation, which is a systematic method to derive the macroscopic equations for the relevant set of slow fields. B. Kinetic equation The distribution function evolves in time by the motion of swimmers, the stochastic evolution of X, and tumbling. All these processes can be captured by the kinetic 4 equation ∂f ∂t +Vˆ n·∇f= F[f] + ˙ l τ ∂(Bf) ∂X !+T[f],(7) which describes the temporal evolution of the distribution function f(r,ˆ n, X, t) in d(2 or 3) spatial dimensions, In Eq. (7) and subsequent equations, unless necessary for clarity or to avoid ambiguity, we will omit the Xand largument of the functions to simplify the notation The left hand side (LHS) of the equation describes the persistent motion of bacteria with constant speed V=Vˆ n. The parenthesis on the right hand side (RHS) is a Fokker– Planck term corresponding to the Langevin equation (3) that describes the evolution of CheY-P. Here, F[f]≡1 τ∂2f ∂X2+∂(Af) ∂X =1 τb F[f] (8) is associated to the free evolution of X, while the second term in the parenthesis of Eq. (7) accounts for the coupling of CheY-P with the ligand. The last term in Eq. (7) is T[f]≡ν0C(X)Zdˆ n′w(ˆ n′·ˆ n)f(r,ˆ n′, X, t)−f(9) and describes the tumbling process as a Lorentz-type equation having a gain term and a loss term, with a tumble rate determined by Eq. (4). For simplicity, the differential for the integral over the directors will be denoted simply by dˆ ninstead of the more formal notation dd−1ˆ n. The tumbling kernel w, which gives the probability for a director ˆ n′to change into ˆ n, is assumed by isotropy to depend only on the relative angle between the two directors, and it is normalized such that Rdˆ nw(ˆ n′·ˆ n) = 1. In the next sections, we will show that main results will depend only on the first moment of w α1=Zdˆ n′w(ˆ n′·ˆ n)ˆ n′·ˆ n,(10) entering only as a numerical factor in the transport coefficients. A simple choice for wis to take it totally isotropic, meaning that the director after tumbling is chosen completely at random in the unit sphere. That is, wisotropic = 1/Ωd, where Ωdis a the area of the ddimensional unit sphere (Ω2= 2πand Ω3= 4π). In this case, α1,isotropic = 0. The tumbling of E. coli is not fully isotropic, with a larger preference toward persisting angles, resulting in α1,E.coli ≈0.33 [1]. Experiments show that the ligand concentration gradient also affects the tumbling process and the kernel wis also a functional of l[50–52]. Here, we will not consider these kind of models, which will be let for future work. Kinetic equations similar to Eq. (7) have been previously proposed that consider also the internal protein concentration as a relevant variable [53, 54]. In Ref. [53] the hydrodynamic equation are derived from the kinetic equation for times much longer than the memory time, resulting in the Keller–Segel equations with values of the transport coefficients that depend on the memory time. In Ref. [54] the kinetic equation is solved for stationary gradients, being able to obtain the chemotactic currents and their transient dynamics. The kinetic equation can be made dimensionless by taking ν0=V= 1 and b≡B(0) = 1, election that fixes the units of time, length, and ligand concentration. In this case, the most important dimensionless parameters of the model are the tumbling rate sensitivity λ≡dC(X)/dX|X=0, the dimensionless memory time ˆτ=ν0τ, and the dimensionless intensity of the ligand gradient c ∇l=bV ∇l. The tracking of E. coli in [22] allowed to fit its stochastic dynamic to the linear models, giving ˆτ= 4.2, meaning that bacteria keep the concentration of CheY-P and, therefore, the tumbling rate, constant in average for four tumbling events. As it has not been proven that ˆτis larger than one for all flagellated bacteria, we will analyze both the cases of large and small memory times. Finally, to help readability and the interpretation of the terms, in general we will keep dimensions in what follows. The eigenvalues of the Fokker–Planck operator, F[Un] = −γn τUn(11) satisfy 0 = γ0< γ1< γ2. . . and the eigenfunctions are Un(X) = un(X)ϕ(X), where ϕ(X)≡ϕ0 eΦ(X) Ωd ,(12) is the stationary solution of Eq. (8), with Φ(X) = −ZX dX′A(X′),(13) and ϕ0the normalization constant such that u0= 1 and Rdˆn RdXϕ(X)un(X)up(X) = δnp [55]. For the linear model ϕ=e−X2/2 Ωd√2π,un(X) = Hn(X/√2)/√2nn!, with Hn the Hermite polynomial of order n[such that H0(x) = 1, H1(x)=2x, . . . ][56], and γn=n[24]. When ˆτ≪1, the concentration of CheY-P responds rapidly to variations of the ligand and to fluctuations. Therefore, the relevant macroscopic field is the bacterial density ρ(r, t)≡Zdˆ nZdXf(r,ˆ n, X, t),(14) which is a conserved slow field. Associated to it, is bacterial current J(r, t)≡Zdˆ nZdXV ˆ nf(r,ˆ n, X, t).(15) On the other side, if the memory time is long, i.e. ˆτ≫1, the CheY-P concentration remains constant for 5 several tumbling events. The ordering of the eigenvalues of Fimply that the first CheY-P moment will be the next slowest mode. Therefore, in this case, besides ρand J, the relevant field for a macroscopic description is the CheY-P density ρX(r, t)≡Zdˆ nZdXu1(X)f(r,ˆ n, X, t),(16) which can be considered as a slowly evolving field. The associated current is JX(r, t)≡Zdˆ nZdXu1(X)f(r,ˆ n, X, t)Vˆ n.(17) As in this system there is no spontaneously broken symmetries or critical fields, there are no additional relevant fields [57, 58]. Starting from the kinetic equation (7), the objective of this article is to derive macroscopic equations for the relevant fields, as extensions to the Keller–Segel equations. For that, we will apply the Chapman–Enskog method, which is a systematic approach, where the macroscopic dynamical equations are obtained at different time scales [59, 60]. In Sect. III we will consider the case of short memory, where ρis the only conserved field and the application of the Chapman–Enskog method is standard. As it can be anticipated, the Keller–Segel model is recovered in this regime and the method provides expressions for the transport coefficients from the microscopic dynamics. In Sect. IV the case of long memory is worked out. Here, the slow fields are ρand ρXbut, as the CheY-P density is not strictly conserved, we need to apply a modified Chapman–Enskog procedure similar to that used in the study of granular gases, where energy is not conserved [61, 62]. The outcome of this analysis will be a coupled set of equations for ρand ρX, where memory manifests in the form of nonlocal responses. C. General mathematical properties Before obtaining the conservation equations, we study the properties of the Fokker–Planck and tumbling operators of the kinetic equation (7). First, we define the linear operators L0and L L0[f]≡F[f] + T[f],(18) L[f]≡F[f] + ˙ l τ ∂(Bf) ∂X +T[f],(19) where in L0we leave apart the coupling term because it has different symmetry properties. Note that the operators F,T,L0, and Lhave all units of inverse of time. Now we prove that L0is hermitian and semidefinite negative under the scalar product ⟨f|g⟩ ≡ Zdˆn ZdXϕ−1fg, (20) It is direct to show that ⟨g|F[f]⟩=1 τZdˆ nZdXϕ−1g∂2f ∂X2+∂(Af) ∂X  =−1 τZdˆ nZdXϕ−1∂g ∂X +Ag∂f ∂X +Af,(21) where we consider that the functions must decay in the limits (natural boundary conditions of Ref. [55]) and we use that Φ(X)′=−A(X), by using Eq. (13). For the tumbling part, we get ⟨g|T[f]⟩=−ν0 2Zdˆ nZdˆ n′ZdXϕ−1C ×w(ˆ n′·ˆ n) [g(ˆ n′, X)−g(ˆ n, X)] [f(ˆ n′, X)−f(ˆ n, X)] . (22) So, assuming that C(X)>0, the previous results imply that Fand T(and hence also L0) are hermitian and semidefinite-negative. It is direct to verify by simple substitution that L0[ϕ] = 0, implying that feq(ˆ n, X) = ρ0ϕ(X) is a stationary solution Eq. (7) in absence of the ligand. We now show this is the only stationary solution and is associated to the conservation of density. Let us assume that there is a different solution of L0[f] = 0, which we write as f(ˆ n, X) = ϕ(X)φ(ˆ n, X). It is direct that Zdˆ nZdXφ(ˆ n, X)L0[f] = ⟨f|F[f]⟩+⟨f|T[f]⟩.(23) As the LHS vanishes because fbelongs to the Ker of L0 and using that both Fand Lare semidefinite-negative, one obtains that ⟨f|F[f]⟩=⟨f|T[f]⟩= 0. As it is standard for the Fokker–Planck operator or from Eq. (21), the only solution of ⟨f|F[f]⟩= 0 is that f(ˆ n, X) = ϕ(X)φ(ˆ n). Also, from Eq. (22), ⟨f|T [f]⟩can only vanish if φis independent of ˆ n, that is a constant. With this, we have obtained that there is no other solution apart from ϕ. This result, together with the semidefinite-negative character of L0, implies that in absence of a ligand signal, the system relaxes toward the equilibrium distribution feq which is isotropic. The dynamical evolution of the different fields are obtained by computing the moments of the kinetic equation, that is, multiplying Eq. (7) by g(ˆ n, X) and integrating over ˆ nand X. Conserved fields are those such that the right hand side vanishes for all distribution functions f. For example, using g= 1 gives the density field, which is conserved. We now show that the bacterial density is the only conserved field for this system. For the tumbling operator, we note that Rdˆ ndX g(ˆ n, X)T[f] = ⟨gϕ|T [f]⟩, which using Eq. (22) gives that it vanishes for any fonly if gis independent of ˆ n. On the other hand, for the Fokker–Planck term with coupling, Zdˆ nZdXg(ˆ n, X) F[f] + ˙ l τ ∂Bf ∂X ! =−1 τZdˆ nZdX ∂g ∂X ∂f ∂X +fA +B˙ lf.(24) 6 Eq. (24) vanishes for all fonly if gis independent of X. These two results imply that when computing moments of the kinetic equation (7), the right hand side vanishes only if gis simultaneously independent of Xand ˆ n, that is g= 1. Therefore, density is the only conserved field. Importantly, Tand Fare isotropic (both commute with the rotation operator). This means that, for the angular dependence, solutions of L0can be expressed as a sum in Fourier series in two dimensions or spherical harmonics in three dimensions. But, Tand Fdo not commute, and we cannot use a common basis in Xfor L0that simplifies the analysis. We choose to use the eigenbasis (11) , which is a diagonal basis for the Fokker–Planck operator. Then, we will derive the hydrodynamic equations using solutions of the kinetic equation of the form f(ˆn, X) = X n,m |fnm⟩=X n,m fnmeinθUm(X).(25) Although the calculations will be performed both in two and three dimensions, for concreteness, we will present expressions in two dimensions, as in this case. In three dimensions the expansion is analogous, with the use of spherical harmonics instead of the Fourier angular modes. Applying L0over each element of the basis, |fmn⟩, we get, L0|fnm⟩=−γm τUm(X)einθfnm +ν0Um(X)C(X)Zd2θ′w(θ−θ′)einθ′−einθfnm. (26) Using the inner product defined in Eq.(20), we find the elements of operator matrix, Fn′m′ nm ≡ ⟨fn′m′|L0|fnm⟩=−cmm′n τδnn′,(27) where cmm′n=ν0τbmm′(1 −αn) + γmδmm′.(28) Here αn=Zdˆ n′w(ˆ n′·ˆ n)(ˆ n′·ˆ n)n,(29) where we note that, by normalization, α0= 1. The matrix bmm′= ΩdZ∞ −∞ dXϕ−1C(X)Um(X)Um′(X),(30) depends on the explicit chemotactic model and can be easily computed one the functions A,B, and Care given. For example, considering the linear model (5), this matrix is bmm′(λ) = 1 √2π2m+m′m!m′! ×Z∞ −∞ eλX−X2/2Hm(X/√2)Hm′(X/√2)dX, =eλ2/2   1λ··· λ1 + λ2··· . . .. . ....   ,(31) where we take into account that for linear model, the eigenbasis of the Fokker–Planck equation are the Hermite polinomials. This matrix is nondiagonal, but it becomes close to diagonal for small values of λ. Then, it is expected that the series (25) can be truncated with few number of Hermite polynomials for small λ. Consequently, this assumption will be considered to present some explicit expressions. In this presentation, we are not considering rotational diffusion as it is a process that does not respond to ligand concentration and is therefore not key to the chemotactic process. It can, however, modify the numerical values, hindering the chemotactic efficiency by introducing additional randomness into the bacterial motion. Anyhow, it can be included by simply adding a term Drτn2 in Eq. (28), where Dris the rotational diffusion coefficient. For the strain of E. coli studied in Ref. [22], the best fit gave that the dimensionless rotational diffusivity is ν0Dr≈0.12, and is therefore subdominant in the analysis. III. CHAPMAN–ENSKOG METHOD FOR SHORT-TIME MEMORY A. Chapman–Enskog expansion In this section, we consider short memory times, i.e., ˆτ≲1. This implies that the bacterial density ρis the only slow field. Integrating Eq. (7) over ˆ nand X, the RHS vanishes and we obtain the conservation equation associated with the density of bacteria, ∂ρ ∂t +∇·J= 0.(32) This equation is not closed for ρ, since it depends on the yet unknown Jfield (15), which in turns depends on f. The objective here is to find an explicit functional form f=f[ρ]. To do that, we assume that Eq. (7) admits normal solutions, i.e., that the normal solutions of the distribution function, f, vary over time and space through ρ, f(r,ˆ n, X, t) = f[ˆn, X|ρ(r, t)].(33) To find this functional dependence and, consequently, the macroscopic equation, we assume that spatial gradients 7 are small. For that, we introduce a small bookkeeping parameter ε, which characterizes the spatial gradients. In this regime, the conserved field is therefore a slowly evolving field. The distribution function can depend on ρand its gradients, which also vary slowly in space and time. The distribution function is represented as a series expansion in ε, f=f(0) +εf(1) +ε2f(2) +··· .(34) This expansion introduces some arbitrariness when fixing the constants and the number of solutions increases. In the Chapman–Enskog method [61, 62], this arbitrariness is solved by demanding that the ε0= 1 order exactly reproduces the hydrodynamic fields, i.e. Zd2ˆn ZdXf(k)=ρδ0k,∀k. (35) Next, we introduce a separation of time scales, t0= t, t1=εt, t2=ε2t, . . . , where we recall that we made time dimensionless by taking ν0= 1. Each temporal variable reflects the dynamics of different time scales (see Fig. 1). For example, on the tumbling time scale characterized by t=O(1), t1and t2will be small and the distribution function will evolve only with t0, with no dependence on t1or t2. For t=O(ε−1), on the hydrodynamic time scale linear with gradients, the dynamics will be described only by t1, which takes values of order unity, while t0has already saturated to large values and t2still has negligible values. Finally, the slowest scale t=O(ε−2), describes with t2the hydrodynamic regime that depends quadratically on the spatial gradients [60]. To change between times scales smoothly and that the expansion is well-defined mathematically, we impose that the distribution function is regular as ti→0 and ti→ ∞, for i= 0,1,2, . . . . We can thus write the distribution function as f=f[ˆn, X|ρ(r;t0, t1, . . . )]. Using the chain rule for temporal derivatives, we have ∂f ∂t →∂f ∂t0 +ε∂f ∂t1 +ε2∂f ∂t2 +. . . . (36) FIG. 1. Representation of the different time scales of the system. t0represents the time scale at which kinetic effects are relevant (in our case, tumbling time scale). In this regime, the distribution function evolves towards the local steady-state distribution. t1corresponds to the reversible hydrodynamic scale. Finally, t2is the slowest time scale, associated to diffusive processes. Consistent with the previous discussion, we rewrite the kinetic equation (7) as ∂f ∂t +εV ˆn ·∇f=L0[f] + ε˙ l τ ∂(Bf) ∂X ,(37) where εis also introduced in the Lagrangian derivative of l, meaning that the driving gradient evolves slowly in space and time, which is a plausible assumption since this field is externally controlled and molecular diffusion rapidly smoothes any initial large gradient. Substituting the series expansion (34) and the multiple time scales (36) into the kinetic equation (37), results in a hierarchy of equations for each order in ε. B. Zeroth order equation Considering the terms of order zero for εin Eq. (37) gives, ∂f(0) ∂ρ ∂ρ ∂t0 =L0[f(0)],(38) where we used that fdepends on time via the density field. Taking the first moment, that is, multiplying Eq. (38) by 1 and integrating over ˆ nand Xgives ∂ρ ∂t0 = 0,(39) where we use the properties shown in Sec. II C. This result means that ρdoes not change in the kinetic time scale, as expected, because neither tumbling or the evolution of Xchange the particle positions. Inserting back this result into Eq. (38), gives a closed equation for f(0), L0[f(0)]=0.(40) In Sec. II C we showed that the solution of this equation is f(0)(r,ˆn, t) = ρ(r, t)ϕ(X), where ρhas an arbitrary spatio-temporal dependence that will be determined in the next sections when analyzing the kinetic equation to higher orders in ε. C. First order equation Following the same steps, we obtain the equation to first order in ε ∂f(1) ∂t0 +∂f(0) ∂t1 +Vˆn ·∇f(0) =L0hf(1)i+˙ l τ ∂(Bf(0)) ∂X . (41) The first term in (41) vanishes trivially because fdepends on time via ρ, and this field does not depend on t0. Taking the first moment of the equation gives ∂ρ ∂t1 = 0,(42) 8 where the convective term gives a vanishing contribution by parity in ˆ nand the right hand side vanishes by the properties of L0and because f(0) vanishes for Xgoing to infinity. This result means that ρdoes not evolve neither on this time scale. Considering these results and calculating explicitly the terms that depend on f(0) in Eq. (41), we find a closed equation for f(1), L0[f(1)] =Vˆn ·h∇ρ+ρ τE1∇liϕ(X),(43) with E1(X)≡[A(X)B(X)−∂XB(X)] ,(44) where we used that Φ′(X) = −A(X). Here we are considering that the food is stationary, i.e., l≡l(r) to simplify the analysis of the solutions of f(1). Nevertheless, as we show below, the term associated with the explicit time derivative of the ligands does not contribute to the driftdiffusion dynamics, so it can be safely omitted. Using the isotropy of L0, we propose solutions of the form f(1) =V ν0 ˆn ·hO(X)∇ρ+Q(X)ρ τ∇li,(45) where the prefactor ν0has been included to make the unknowns Oand Qdimensionless. In principle, as in the kinetic theory of gases and granular gases, Oand Qare isotropic functions that could depend also on the magnitude of the vector. As ˆ nis unitary, this dependence is not present here. To obtain these functions explicitly, we introduce the expression of f(1) into Eq. (43), resulting in the two equations L0[Oˆn] = ν0ϕˆn,(46a) L0[Qˆn] = ν0E1ϕˆn.(46b) Notably, the isotropy of the operators imply that for any function ξ(X), L0[ξˆn] = ˆn τb L′ 0[ξ], where b L′ 0≡b F − ν0τ(1 −α1)C(X). Then, expanding in Fokker-Planck eigenfunctions, O(X) = PmOmUm(X), and similarly for Q, and applying the orthogonality conditions, a set of equations is obtained for the coefficients Omand Qm, for m= 0,1, . . . X m cmn1Om=−ν0τδn0,(47a) X m cmn1Qm=−ν0τΩdZdXunE1ϕ. (47b) The solution of these equations depends on the specific model, as specified by the model functions A,B, and C. For example, in the case of the linear model, truncating at m= 1, that is, keeping two polynomials, O0≈ − 1 eλ2/2"1 + λ2ˆτeλ2/2 (1 + eλ2/2ˆτ)#,(48a) O1≈λˆτ 1 + eλ2/2ˆτ,(48b) Q0≈bλˆτ 1 + eλ2 2ˆτ,(48c) Q1≈ − bˆτ 1 + eλ2 2ˆτ,(48d) where, for simplicity, in the above expressions, we have considered the case of totally isotropic tumbling, that is α1= 0. D. Hydrodynamic equation The kinetic equation to second order in εreads ∂f(2) ∂t0 +∂f(1) ∂t1 +∂f(0) ∂t2 +Vˆn·∇f(1) =L0hf(2)i+˙ l τ ∂(Bf(1)) ∂X . (49) As the density does not depend on t0and t1, the first and second terms vanish identically. Now, computing the first moment of the equation, gives again a vanishing right hand by the properties of L0and because f(1) vanishes for Xgoing to infinity. But, as f(1) is odd in ˆ n[Eq. (45)], the convective term now gives a finite contribution and the moment equation reads ∂ρ ∂t2 +∇·J= 0,(50) with the bacterial current J J(r, t) = Zdˆ nZdXf(1)(r,ˆ n, X, t)Vˆ n, =−D∇ρ+µρ∇l, (51) where Eq. (45) for f(1) has been used. Here, Dis the diffusion coefficient D=−V2Ωd dν0ZdX O(X) = −V2O0 dν0 (52) and µthe mobility, µ=bV 2Ωd dν0τZdX Q(X) = bV 2Q0 dˆτ.(53) The close expressions (52) and (53), depending on only one coefficient, result from the orthogonality of the base functions. Finally, if we had considered an explicit time dependence of the ligand in Eq. (43), f(1) would have had an additional term proportional to ∂tl, isotropic in ˆ n. But this term, would not contribute to the current when integrating over ˆ nin (51). 15 with ψ0=D12(D11µ21 −D12µ11) D11(D11D22 −D2 12),(99a) ψ1=D22µ11 −D12µ21 (D11D22 −D2 12),(99b) k0=sD11γ1 τ(D11D22 −D2 12).(99c) The response function Ψ(k, 0), has a constant term, associated to a local response, added to a Lorentzian form, indicating a nonlocal response with a smoothing length L0=k−1 0. In the memoryless limit (τ→0) the density response function becomes Ψ(k, 0) = ψ0+ψ1= µ11/D11 =µ/D, independent of k, indicating that the response is completely local and we recover the response of the Keller–Segel model [5]. In the linear model, the symmetry µ21 =µ11D22/D12 of the transport coefficients result in that ψ0=µ11/D11 and ψ1= 0, implying that the density response is purely nonlocal. Figure 4 shows the amplitude of the density response ψ0and the smoothing length L0for the linear model using the numerical solution of the equations (80) for the transport coefficients. Considering only up to order n= 1 in the Hermite series, we have the explicit expressions for the linear case ψ0=bν0(1 −α1)λeλ2/2 1 + (1 −α1)(1 + λ2)ˆτeλ2/2,(100) L0=V τ p1 + (1 −α1)(1 + λ2)ˆτeλ2/2.(101) Both the full numerical solution and the approximate analytical expressions indicate that the response amplitude decays with memory while that the smoothing length grows with memory. Indeed, memory makes that the agents miss the precise position of the chemotactic signal because, even if the reach an optimal position, they continue moving persistently with reduced tumbling rate in a wrong direction. As expected, ψ0grows with the sensitivity λ, but the dependence of the smoothing length with λis weaker. To exemplify the character of the response and compare with agent-based simulations, we consider the simple case of a stationary signal with a step profile in one direction, l(x) = l0+l1sgn(x), with sgn(x)≡x/|x|the sign function. Here, the stationary bacterial density and CheY-P concentration can be explicitly obtained by calculating the inverse Fourier transform in Eq. (98a). Note that although the response functions [Eq. (98a)] are valid for small k, the non-local effect decays with kand the Fourier transform converges for the imposed step function, giving ρ(x) = ρ0+ρ0l1sgn(x)hψ1+ψ0(1 −e−k0|x|)i,(102) ρX(x) = −ρ0 D11 D12 ψ0l1sgn(x)e−k0|x|.(103) 15 0 2 4 6 8 10 0.0 0.5 1.0 1.5 2.0 0 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 FIG. 4. Amplitude 0and smoothing length L0of the static response function [Eq. (98a)] for the linear model as a function of the dimensionless memory time ˆ⌧for =0.1 (blue), 1.0 (red), and 3.0 (green). The solid lines correspond to evaluating the transport coefficients using the numerical solution of the system of equations (80) truncating it up to n=10. In all plots, we consider a fully isotropic tumbling kernel, with ↵1=0. response 0and the smoothing length L0for the linear model using the numerical solution of the equations (80) for the transport coefficients. Considering only up to order n= 1 in the Hermite series, we have the explicit expressions for the linear case 0=b⌫0(1 ↵1)e2/2 1+(1↵1)(1 + 2)ˆ⌧e2/2,(100) L0=V⌧ p1+(1↵1)(1 + 2)ˆ⌧e2/2.(101) Both the full numerical solution and the approximate analytical expressions indicate that the response amplitude decays with memory while that the smoothing length grows with memory. Indeed, memory makes that the agents miss the precise position of the chemotactic signal because, even if the reach an optimal position, they continue moving persistently with reduced tumbling rate in a wrong direction. As expected, 0grows with the sensitivity , but the dependence of the smoothing length with is weaker. To exemplify the character of the response and compare with agent-based simulations, we consider the simple case of a stationary signal with a step profile in one direction, l(x)=l0+l1sgn(x), with sgn(x)⌘x/|x|the sign function. Here, the stationary bacterial density and CheY-P concentration can be explicitly obtained by calculating the inverse Fourier transform in Eq. (98a). Note that although the response functions [Eq. (98a)] are valid for small k, the non-local e↵ect decays with kand the Fourier transform converges for the imposed step function, giving ⇢(x)=⇢0+⇢0l1sgn(x)h 1+ 0(1 ek0|x|)i,(102) ⇢X(x)=⇢0 D11 D12 0l1sgn(x)ek0|x|.(103) Alternatively, these expressions can also be obtained by directly solving the linearized hydrodynamic equations in real space. We recall that the solution of the hydrodynamic equations in real or Fourier space, or even their numerical solution, is much simpler than to directly deal with the kinetic equation for the same geometry as in Ref. [24]. To compare, the Keller–Segel model gives for the density, the local response ⇢(x)=⇢0eµl(x)/D ⇡ ⇢0+⇢0 0l1sgn(x), where the linear response was given in the second expression. To verify the analytical predictions, agent-based simulations are performed in two dimensions for the linear model: bacteria move in two spatial dimensions and the directors are characterized by a single angle, which, at tumbles, is sorted uniformly in the interval [0,2⇡]. That is, the tumbling kernel is isotropic, implying ↵1= 0. These simulations involve the numerical solution of the stochastic equation (2) for Xand the position vector r⌘xˆ x+yˆ yusing the Euler–Heun scheme, following the method used in Ref. [24] to reproduce the dynamics. The system is initialized homogeneously and, after a transient to reach the steady state (longer than the memory time ⌧and the di↵usion time L2/D,whereLis the box size), ⇢(x), ⇢X(x), hXi(x), and h⌫i(x) are measured. The latter is directly measured counting the number of tumbles on x bins. Figure 5 presents the comparison between the simulation results and the theoretical predictions for the normalized bacterial excess density ⇢/⇢0⌘[⇢(x)⇢0]/⇢0, and the normalized protein density ⇢X(x)/⇢0,showinga very good agreement. Both in the theoretical and numerical results, the bacterial density displays a rounded profile, compared with the step response of the Keller–Segel model, and ⇢X presents a sharp discontinuity. At the hydrodynamic level, the steady state is characterized by J= 0. The Keller–Segel model [Eq. (1b)] implies that ⇢becomes flat as soon as lis uniform, giving rise to the discontinuous density profile. Considering memory, a vanishing mass flux gives D11r[(1 D12hXi/D11)⇢]=µ11⇢0rl, (104) where the linearized form of Eq. (93a) was used. Then, over the distance where hXi⌘⇢X/⇢has not relaxed, ⇢will not be uniform even if lis uniform. This is precisely the smoothing distance L0. This smoothing length can be well understood also in terms of the microscopic chemotactic model. For concreteness, consider a bacterium that crosses the boundary from the poor to the rich ligand region. In this case, lexperiences an abrupt increase, meaning that ˙ lis characterized by a Dirac delta function with positive weight. As a consequence, Xsuffers an instantaneous decrease as given by Eq. (3). This abrupt change in Xmanifests as a discontinuity in ⇢X. After that, as the ligand does not change, Xwill relax on a time ⌧to zero, which gives a length scale of the order of L0⇠V⌧. During this relaxation time, the tumbling rate is smaller than average, implying that bacteria moving from the poor to the rich ligand region tumble FIG. 4. Smoothing length L0(left) and amplitude ψ0(right) of the static response function [Eq. (98a)] for the linear model as a function of the dimensionless memory time ˆτfor λ= 0.1 (blue), 1.0 (red), and 3.0 (green). The solid lines correspond to evaluating the transport coefficients using the numerical solution of the system of equations (80) truncating it up to n= 10. In all plots, we consider a fully isotropic tumbling kernel, with α1= 0. Alternatively, these expressions can also be obtained by directly solving the linearized hydrodynamic equations in real space. We recall that the solution of the hydrodynamic equations in real or Fourier space, or even their numerical solution, is much simpler than to directly deal with the kinetic equation for the same geometry as in Ref. [24]. To compare, the Keller–Segel model gives for the density, the local response ρ(x) = ρ0eµl(x)/D ≈ ρ0+ρ0ψ0l1sgn(x), where the linear response was given in the second expression. To verify the analytical predictions, agent-based simulations are performed in two dimensions for the linear model: bacteria move in two spatial dimensions and the directors are characterized by a single angle, which, at tumbles, is sorted uniformly in the interval [0,2π]. That is, the tumbling kernel is isotropic, implying α1= 0. These simulations involve the numerical solution of the stochastic equation (2) for Xand the position vector r≡xˆ x+yˆ yusing the Euler–Heun scheme, following the method used in Ref. [24] to reproduce the dynamics. The system is initialized homogeneously and, after a transient to reach the steady state (longer than the memory time τand the diffusion time L2/D, where Lis the box size), ρ(x), ρX(x), ⟨X⟩(x), and ⟨ν⟩(x) are measured. The latter is directly measured counting the number of tumbles on x bins. Figure 5 presents the comparison between the simulation results and the theoretical predictions for the normalized bacterial excess density ∆ρ/ρ0≡[ρ(x)−ρ0]/ρ0, and the normalized protein density ρX(x)/ρ0, showing a very good agreement. Both in the theoretical and numerical results, the bacterial density displays a rounded profile, compared with the step response of the Keller–Segel model, and ρX presents a sharp discontinuity. At the hydrodynamic level, the steady state is characterized by J= 0. The Keller–Segel model [Eq. (1b)] implies that ρbecomes flat as soon as lis uniform, giving rise to the discontinuous 16 density profile. Considering memory, a vanishing mass flux gives D11∇[(1 −D12⟨X⟩/D11)ρ] = µ11ρ0∇l, (104) where the linearized form of Eq. (93a) was used. Then, over the distance where ⟨X⟩ ≡ ρX/ρ has not relaxed, ρwill not be uniform even if lis uniform. This is precisely the smoothing distance L0. This smoothing length can be well understood also in terms of the microscopic chemotactic model. For concreteness, consider a bacterium that crosses the boundary from the poor to the rich ligand region. In this case, lexperiences an abrupt increase, meaning that ˙ lis characterized by a Dirac delta function with positive weight. As a consequence, Xsuffers an instantaneous decrease as given by Eq. (3). This abrupt change in Xmanifests as a discontinuity in ρX. After that, as the ligand does not change, Xwill relax on a time τto zero, which gives a length scale of the order of L0∼V τ. During this relaxation time, the tumbling rate is smaller than average, implying that bacteria moving from the poor to the rich ligand region tumble less than those moving in the opposite direction. This imbalance, on the length scale L0, produces a net bacterial flux toward the ligand rich region that will stop once the resulting excess density generates a diffusive flux in the opposite direction, such as that the total bacterial flux vanishes. The final density profile will therefore be smoothed over the same length scale. Experimentally, the bacterial density is measurable, but ρXcan be difficult to obtain. A related observable, which is accessible to experiments [69, 70] is the positiondependent average tumbling rate ⟨ν⟩(r), which in the linear model equals ν0⟨eλX ⟩(r). To obtain it in the present formalism, we use the cumulant-generating function [71], which for the random variable Xwith parameter λcan be defined as K(λ) = ln ⟨eλX ⟩. Expanding this expression in Taylor series, we obtain K(λ) = ∞ X n=1 κn λn n!=λ⟨X⟩+λ2 2σ2 X+··· ,(105) where the omitted terms depend on higher cumulants κnof X. Taking into account that in the linear model Xis well described by a normalized Gaussian variable, i.e. the only nonzero cumulants are the first and the second, where the standard deviation is σX= 1, it results ln (⟨ν/ν0⟩(r)) = λ⟨X⟩(r) + λ2 2. Finally, we obtain for the average tumbling rate ⟨ν⟩(r) = ν0eλ2/2eλ⟨X⟩(r),(106a) =ν0eλ2/2eλρX(r)/ρ(r),(106b) where we used that ⟨X⟩=ρX/ρ. In the case of ligand with a step profile, ⟨X⟩=D11l1µ11sgn(x) D12D11 −D12l1µ11sgn(x)(1 −ek0|x|).(107) 14 which for the random variable Xwith parameter can be defined as K()=lnheXi.(104) Expanding this expression in Taylor series, we obtain K()= 1 X n=1 n n n!=hXi+2 22 X+···,(105) where the omitted terms depend on higher cumulants n of X. Taking into account that Xis well described by a normalized Gaussian variable, i.e. the only nonzero cumulants are the first and the second, where the standard deviation is X= 1, it results ln (h⌫/⌫0i(r)) = hXi(r)+2 2. Finally, we obtain for the average tumbling rate h⌫i(r)=⇣⌫0e2/2⌘ehXi(r),(106a) =⇣⌫0e2/2⌘e⇢X(r)/⇢(r),(106b) where we used that hXi=⇢X/⇢. In the case of ligand with a step profile, hXi=D11l1µ11sgn(x) D12D11 D12l1µ11sgn(x)(1 ek0|x|).(107) The average tumbling rate is displayed in Fig. 6, where we compare the values obtained in the simulation with the expression (106a) by using the explicit expression of hXigiven by Eq. (107), and with the expression (106b) using the measured values of ⇢(x) and ⇢X(x). The agreement between the simulations and the theoretical results is good, with better agreement using the simulation measures of hXi. D. Di↵usive dynamics As a final case of interest, we study the bacterial dynamics in absence of a chemotactic signal. In this case, Eqs. (52) and (53), with the transport relations Eqs. (71) and (78), reduce to @ @t✓˜⇢ ˜⇢X◆=M✓˜⇢ ˜⇢X◆.(108) for the Fourier modes ˜⇢⌘Zdr⇢eik·r,(109) ˜⇢X⌘Zdr⇢Xeik·r.(110) The dynamic matrix Mis M⌘✓D11k2D12k2 D12k2D22k2+1 ⌧◆,(111) 10.07.55.02.5 0.0 2.5 5.0 7.5 10.0 x 0.075 0.050 0.025 0.000 0.025 0.050 0.075 /0 10.07.55.02.5 0.0 2.5 5.0 7.5 10.0 x 0.4 0.2 0.0 0.2 0.4 X/0 FIG. 5. Stationary normalized density (top) and protein density (bottom) profiles generated for a step function signal. The (blue) circles represent the results of the simulations of particles with ˆ⌧=1.0, =1.0. In all cases, the particles are confined in a square box of size L= 20 with periodic boundary conditions, and with a chemotactic pulse of amplitude l1=0.4. The (red) solid lines are the theoretical predictions. The case of fully isotropic tumbling kernel, with ↵1=0,has been considered. Units have been fixed to V=⌫0=b=1. FIG. 6. Tumble rate average profile for the case described in Fig. 5. The blue circles represent the results of the simulations, the red solid line is the theoretical prediction Eq. (106a) by using hXigiven by Eq. (107), and the green dashed line is the theoretical prediction (106b) using the measured values of ⇢(x)and⇢X(x). where we used that D21 =D12. This matrix is symFIG. 5. Stationary normalized density (top) and protein density (bottom) profiles generated for a step function signal in the linear model. The blue circles represent the results of the simulations of particles with ˆτ= 1.0, λ= 1.0. In all cases, the particles are confined in a square box of size L= 20 with periodic boundary conditions, and with a chemotactic pulse of amplitude l1= 0.4. For these parameters L0= 0.43. The red solid lines are the theoretical predictions and the solid dashed line the prediction of the KS model. The case of fully isotropic tumbling kernel, with α1= 0, has been considered. Units have been fixed to V=ν0=b= 1. The average tumbling rate is displayed in Fig. 6, where we compare the values obtained in the simulation with the expression (106a) by using the explicit expression of ⟨X⟩given by Eq. (107), and with the expression (106b) using the measured values of ρ(x) and ρX(x). The agreement between the simulations and the theoretical results is good, with better agreement using the simulation measures of ⟨X⟩. D. Traveling wave For a more complex statio-temporal response, here we consider the case of a traveling chemotactic wave l(x, t) = l0eik(x−Vst). In this case, the response is obtained by substituting ω=V k in Ψ and ΨXgiven in Appendix C. Replacing in Eqs. (91) and (93a), gives for the linear response of the chemotactic current J=VsΨρ(k, Vsk)l0,(108) 17 14 which for the random variable Xwith parameter can be defined as K()=lnheXi.(104) Expanding this expression in Taylor series, we obtain K()= 1 X n=1 n n n!=hXi+2 22 X+···,(105) where the omitted terms depend on higher cumulants n of X. Taking into account that Xis well described by a normalized Gaussian variable, i.e. the only nonzero cumulants are the first and the second, where the standard deviation is X= 1, it results ln (h⌫/⌫0i(r)) = hXi(r)+2 2. Finally, we obtain for the average tumbling rate h⌫i(r)=⇣⌫0e2/2⌘ehXi(r),(106a) =⇣⌫0e2/2⌘e⇢X(r)/⇢(r),(106b) where we used that hXi=⇢X/⇢. In the case of ligand with a step profile, hXi=D11l1µ11sgn(x) D12D11 D12l1µ11sgn(x)(1 ek0|x|).(107) The average tumbling rate is displayed in Fig. 6, where we compare the values obtained in the simulation with the expression (106a) by using the explicit expression of hXigiven by Eq. (107), and with the expression (106b) using the measured values of ⇢(x) and ⇢X(x). The agreement between the simulations and the theoretical results is good, with better agreement using the simulation measures of hXi. D. Di↵usive dynamics As a final case of interest, we study the bacterial dynamics in absence of a chemotactic signal. In this case, Eqs. (52) and (53), with the transport relations Eqs. (71) and (78), reduce to @ @t✓˜⇢ ˜⇢X◆=M✓˜⇢ ˜⇢X◆.(108) for the Fourier modes ˜⇢⌘Zdr⇢eik·r,(109) ˜⇢X⌘Zdr⇢Xeik·r.(110) The dynamic matrix Mis M⌘✓D11k2D12k2 D12k2D22k2+1 ⌧◆,(111) 10.07.55.02.5 0.0 2.5 5.0 7.5 10.0 x 0.075 0.050 0.025 0.000 0.025 0.050 0.075 /0 10.07.55.02.5 0.0 2.5 5.0 7.5 10.0 x 0.4 0.2 0.0 0.2 0.4 X/0 FIG. 5. Stationary normalized density (top) and protein density (bottom) profiles generated for a step function signal. The (blue) circles represent the results of the simulations of particles with ˆ⌧=1.0, =1.0. In all cases, the particles are confined in a square box of size L= 20 with periodic boundary conditions, and with a chemotactic pulse of amplitude l1=0.4. The (red) solid lines are the theoretical predictions. The case of fully isotropic tumbling kernel, with ↵1=0,has been considered. Units have been fixed to V=⌫0=b=1. 10.07.55.02.5 0.0 2.5 5.0 7.5 10.0 x 1.4 1.5 1.6 1.7 1.8 1.9 2.0 hi FIG. 6. Tumble rate average profile for the case described in Fig. 5. The blue circles represent the results of the simulations, the red solid line is the theoretical prediction Eq. (106a) by using hXigiven by Eq. (107), and the green dashed line is the theoretical prediction (106b) using the measured values of ⇢(x)and⇢X(x). where we used that D21 =D12. This matrix is symFIG. 6. Tumble rate average profile for the case described in Fig. 5. The blue circles represent the results of the simulations, the red solid line is the theoretical prediction Eq. (106a) by using ⟨X⟩given by Eq. (107), and the green dashed line is the theoretical prediction (106b) using the measured values of ρ(x) and ρX(x). while in the Keller–Segel theory, the response is JKS =kµVs Dk2−ikVs l0.(109) Recent experiments performed with E. coli, showed that the resulting chemotactic current is not monotonic with the wave speed Vs, presenting a maximum at Vs≈ 8µm/s [72]. Figure 7 shows the predicted current using the linear model with the fitted values for E. coli (see Section VI), which is compared with the prediction of the Keller–Segel theory. Although the experiment is performed in a strong non-linear response regime, where the current saturates to large values, it is remarkable that the present model predicts the existence of the maximum for a wave velocity in the same order or magnitude, while the Keller–Segel theory gives a monotonically increasing current failing event to predict the existence of a maximum. These results are consistent with the analysis in Ref. [72], where it is shown that the origin of the maximum is the existence of a finite memory. E. Diffusive dynamics As a final case of interest, we study the bacterial dynamics in absence of a chemotactic signal. In this case, Eqs. (57) and (58), with the transport relations Eqs. (83) and (86), reduce to ∂ ∂t ˜ρ ˜ρX=−M˜ρ ˜ρX.(110) for the Fourier modes ˜ρ≡Zdrρeik·r,(111a) ˜ρX≡ZdrρXeik·r.(111b) Memory model Keller-Segel 0 10 20 30 40 50 0.00 0.02 0.04 0.06 0.08 0.10 0.12 FIG. 7. Normalized bacterial current as a response to a traveling chemotactic wave of speed Vs, computed using the values of the linear model for E. coli. The current is computed for a periodic box of length L= 800 µm, equal to the value in the experiment [72], meaning that Ψρis evaluated at k= 2π/L. In blue the prediction of this article and in red the prediction of the Keller–Segel theory (divided by 5 to help the comparison). The dynamic matrix Mis M≡D11k2−D12k2 −D12k2D22k2+γ1 τ,(112) where we used that D21 =D12. This matrix is symmetric, has real coefficients and is semidefinite positive, implying that it is diagonalizable with positive and real eigenvalues χi, and eigenvectors vi, with i= 1,2. Then, the general solution is of the form ˜ρ ˜ρX= 2 X j=1 civie−χit,(113) with the constants cidetermined by the initial conditions. Figure 8 shows the eigenvalues χ1,2as a function of kin logarithmic scale for the lineal model. The values of the sensitivity to fluctuations and the memory time used for the plot are λ= 1.62 and ˆτ= 4.2, corresponding to the experimental values of E. coli [22]. For small k, there is one eigenvalue that saturates to a constant and another one that goes as k2, reflecting a diffusive mode. For large k, now the two modes are diffusive, with new diffusion constants, one of those being smaller than the small wavevector diffusivity. To describe in more detail this behavior and determine the crossover wavelength k∗, we analyze separately the cases of small and large wavevectors. Expanding the eigenvalues in a Taylor series in k, and keeping terms up to k2, we get χ1=D11k2+O(k4),(114) χ2=γ1 τ+D22k2+O(k4),(115) 18 0.01 0.1 1 10 0.01 1 100 FIG. 8. Relaxation eigenvalues χi,i= 1,2 , as a function of kin log-log scale for the linear model. The blue and red solid lines correspond to χ1and χ2, respectively. The green dashed line represents k2, in order to identify the modes that presents diffusive behavior in the different regimes. The vertical dashed line corresponds to k∗≡(D22τ)−1/2, separating the different dynamical regimes. The values used are λ= 1.62, ˆτ= 4.2, and α1= 0.33. with their corresponding eigenvectors, truncated to leading order in k, v1=1 τD12k2/γ1,(116) v2=−τD12k2/γ1 −1.(117) Consistent with the figure, there is only one diffusive mode with diffusion coefficient D11, associated with the eigenvector v1, which is essentially the density mode. As discussed in Sec. V A, this diffusion coefficient equals the short memory value D. The other mode, which is associated to ρX, relaxes at a constant rate τin the limit of vanishing wavevectors. Note that in this regime, ρand ρXare practically decoupled. To study the regime of high k, we take the limit 1/τ → 0 in the general expression for χ1and χ2, obtaining lim 1/τ→0χ1,2= D11 +D22 ∓p(D11 −D22)2+ 4D2 12 2!k2. (118) The modes are indeed diffusive, with diffusion coefficients D∓given by the term inside the parenthesis. It is direct to verify that D−< D < D+, as was identified in the figure. Figure 9 displays these three diffusion coefficients for the linear model as a function of the model parameters, where it can be verified that D−can be notoriously smaller than D, implying long relaxation times. In this regime, the eigenvectors (not written by directly obtained from the dynamic matrix) couple ρand ρXat the same order, so both fields evolves diffusively. The crossover between the two regimes can be directly obtained by observation of the matrix [Eq. (112)] or from Eq. (115), resulting in k∗≡(D22τ/γ1)−1/2. In the case of large memory times, the associated crossover wavelength is also large. VI. NUMERICAL VALUES FOR E. COLI The parameters of the linear model for E. coli (RP437 bacteria in motility buffer supplemented with serine) were determined in Ref. [22]: V= 27 µm/s, ν0= 0.22 s−1,τ= 19 s, and λ= 1.62. Also, the first moment of the tumble kernel is known: α1≈0.33 [1]. The only remaining parameter is b, which depends on the specific ligand to be considered. With these results, it is possible to provide explicit values for the transport coefficients and related parameters discussed in previous sections. The diffusion coefficients are D11 = 1.3×103µm2/s, D12 =D21 = 0.81 ×103µm2/s, and D22 = 0.99 ×103µm2/s. The chemotactic mobilities are µ11/b = 0.42 ×102µm2/s2, µ12/b = 2.2×102µm2/s2,µ21/b = 0.52 ×102µm2/s2, and µ22/b = 2.6×102µm2/s2. The amplitude and smoothing length of the static response function are ψ0/b = 0.032 s−1and L0= 1.7×102µm, respectively. Finally, the large wavevector diffusivities are D+= 2.0× 103µm2/s and D−= 0.33 ×103µm2/s, which are valid for wavelengths smaller than L∗= 2π/k∗= 1.5×103µm. VII. CONCLUSIONS The bacterial chemotactic response can be rationalized in terms of a small number of variables associated to the internal concentration of relevant proteins inside the bacterial body, which control the tumble rate. In the case of E. coli, there is a strong scale separation between the relaxation rates of these proteins allowing for a reduction of the model to a single variable, the concentration of the CheY-P protein governed by a single memory time τ. Inspired by this model, but considering general nonlinear coupling and responses to the ligand, it is possible to write down a kinetic equation for an ensemble of bacteria in presence of a ligand field (food, chemoattractant, or chemorepellent), with arbitrary spatio-temporal dependence [Eq. (7)]. This kinetic equation considers the evolution of an effective internal variable Xthat govern the chemotactic process, which we assume is controlled by a single large memory time τ, while all other characteristic times are shorter. In the simpler model, Xis the CheY-P protein concentration inside the bacterial body, τis the methylation time of the chemorecetors, and coupling and responses to the ligand are assumed to be simple linear functions of X. This linear model is analyzed throughout the article as it allows to provide explicit results and because its parameters have been measured for E. coli. 19 19 0 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 0 2 4 6 8 10 0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 0 1 2 3 4 5 6 FIG. 9. Small kdimensionless di↵usion coefficient ˜ D=d⌫0D/V 2(blue), and large kdimensionless di↵usion coefficients ˜ D+=d⌫0D+/V 2(red) and ˜ D=d⌫0D/V 2(green) for the linear model as a function of the dimensionless memory ˆ⌧for selected values of the sensitivity . The solid lines correspond evaluating the transport coefficients using the numerical solution of the system of equations (80) truncating it up to n= 10. In all plots, we consider a fully isotropic tumbling kernel, with ↵1=0. derived. They have the form of reaction-di↵usion equations, which are coupled to the ligand field. Associated to these equations, the di↵usion and mobility transport coefficients are obtained in terms of the microscopic parameters. The derived equations are analyzed in simple regimes, aiming to highlight some relevant features. In the case of a uniform and stationary ligand gradient, a stationary chemotactic current is obtained. However, the coupling with ⇢Xgenerates for a transient ⌧, an e↵ective chemotactic mobility larger than the stationary one if the ensemble was previously placed in a region with a strong ligand gradient. Also, analyzing the linear response for signals with spatio-temporal dependence, a nonlocal response in time and space is obtained, which is absent in the Keller–Segel description. The associated smoothing length equals 170 µm for the case of E. coli,whichis comparable to many spatial features in natural and artificial microfluidic environments. Note that rotational di↵usion could reduce slightly the smoothing length. Experimentally, for the linear model, the nonlocal response can be determined by measuring the density profiles and the local tumbling rate when the ensemble is placed in presence of an inhomogeneous signal. Appropriate observables should be defined for the case of other chemotactic models. When analyzing the response to chemoattractant traveling waves, the response current presents a maximum for a particular wave velocity as has been observed in experiments [67], maximum that is absent in the Keller–Segel description. Finally, the di↵usive dynamics in absence of signal is analyzed, showing that the di↵usion coefficient is scale dependent, with a crossover at 1.5⇥103µm. Besides the cases considered in this article, these equations can be used is di↵erent relevant experimental or numerical configurations. The derivation of the hydrodynamic equations and the transport coefficients has been presented in detail for future use, in case other chemotactic and motility models need to be considered, for example, dealing with saturation of the ligand receptors, more complex chemotactic circuits, or bacterial variability (see for example, [9, 10]). Possible extensions to the model analyzed in this manuscript can include the modification of the tumbling kernel with ligand gradient [45–47], or chemokinesis, in which bacteria modify the propulsion speed as a response to ligand concentration [68, 69]. In some cases, only the transport coefficients will change but in others, for example, when dealing with bacterial interactions, also the hydrodynamic equations could be modified. The hydrodynamic equations derived in this manuscript need to be complemented with appropriate boundary conditions. A first approach would be simply to impose non-flux boundary conditions for both fields. However, the problem is far from trivial because bacteria and other microswimmers, either by simple persistence of motion or by hydrodynamic attraction, tend to accumulate on surfaces [70], where they modify their motion [71]. More research is needed to describe the interaction of bacteria with surfaces, considering memory e↵ects. In Refs. [72, 73] the case of bacteria responding to diffusing point sources taking place in natural environments is analyzed, showing that due to the spatio-temporal dependence of the signal, non-trivial responses can emerge. This and similar configurations are perfect examples where the equations derived here can be applied, either numerically or analytically. A proper analysis, considering di↵erent regimes goes beyond the scope of the manuscript and will be devoted to future research. FIG. 9. Small kdimensionless diffusion coefficient ˜ D=dν0D/V 2(blue), and large kdimensionless diffusion coefficients ˜ D+=dν0D+/V 2(red) and ˜ D−=dν0D−/V 2(green) for the linear model as a function of the dimensionless memory ˆτfor selected values of the sensitivity λ. The solid lines correspond evaluating the transport coefficients using the numerical solution of the system of equations (80) truncating it up to n= 10. In all plots, we consider a fully isotropic tumbling kernel, with α1= 0. The kinetic analysis shows that there is a emergent long length scale, which in the case of E. coli can be of several hundreds of micrometers, making it necessary to build a practical framework when the ligand field varies on these scales as it happens in natural an artificial microfluidic environments. Employing the Chapman– Enskog method, we derived the hydrodynamic equations that describe a bacterial ensemble as a formal expansion in spatial gradients. In the case the memory time is small, the methods yields to the standard Keller–Segel equation for the bacterial density ρ, with a chemotactic mobility and diffusion coefficient computed entirely in terms of the microscopic parameters [Eqs. (50-53)]. More relevant is the case of long memory time, as it has been determined experimentally for E. coli. In this case, besides the bacterial density, the density of the internal variable ρXemerges as a new relevant field, and the coupled equations for ρand ρXare derived. They have the form of reaction-diffusion equations, which are coupled to the ligand field [Eqs. (91-93)]. Associated to these equations, the diffusion and mobility transport coefficients are obtained in terms of the microscopic parameters. The derived equations are analyzed in simple regimes, aiming to highlight some relevant features. In the case of a uniform and stationary ligand gradient, a stationary chemotactic current is obtained. However, the coupling with ρXgenerates for a transient τ, an effective chemotactic mobility larger than the stationary one if the ensemble was previously placed in a region with a strong ligand gradient. Also, analyzing the linear response for signals with spatio-temporal dependence, a nonlocal response in time and space is obtained, which is absent in the Keller–Segel description. The associated smoothing length equals 170 µm for the case of E. coli, which is comparable to many spatial features in natural and artificial microfluidic environments. Note that rotational diffusion could reduce slightly the smoothing length. Experimentally, for the linear model, the nonlocal response can be determined by measuring the density profiles and the local tumbling rate when the ensemble is placed in presence of an inhomogeneous signal. Appropriate observables should be defined for the case of other chemotactic models. When analyzing the response to chemoattractant traveling waves, the response current presents a maximum for a particular wave velocity as has been observed in experiments [72], maximum that is absent in the Keller–Segel description. Finally, the diffusive dynamics in absence of signal is analyzed, showing that the diffusion coefficient is scale dependent, with a crossover at 1.5×103µm. Besides the cases considered in this article, these equations can be used is different relevant experimental or numerical configurations. The derivation of the hydrodynamic equations and the transport coefficients has been presented in detail for future use, in case other chemotactic and motility models need to be considered, for example, dealing with saturation of the ligand receptors, more complex chemotactic circuits, or bacterial variability (see for example, [9, 10]). Possible extensions to the model analyzed in this manuscript can include the modification of the tumbling kernel with ligand gradient [50–52], or chemokinesis, in which bacteria modify the propulsion speed as a response to ligand concentration [73, 74]. In some cases, only the transport coefficients will change but in others, for example, when dealing with bacterial interactions, also the hydrodynamic equations could be modified. The hydrodynamic equations derived in this manuscript need to be complemented with appropriate boundary conditions. A first approach would be simply to impose non-flux boundary conditions for both fields. However, the problem is far from trivial because bacteria and other microswimmers, either by simple persistence of motion or by hydrodynamic attraction, tend to accumulate on surfaces [75], where they modify 20 their motion [76]. More research is needed to describe the interaction of bacteria with surfaces, considering memory effects. In Refs. [77, 78] the case of bacteria responding to diffusing point sources taking place in natural environments is analyzed, showing that due to the spatio-temporal dependence of the signal, non-trivial responses can emerge. This and similar configurations are perfect examples where the equations derived here can be applied, either numerically or analytically. A proper analysis, considering different regimes goes beyond the scope of the manuscript and will be devoted to future research. ACKNOWLEDGMENTS This research is supported by Fondecyt Grants No. 1220536 (RS) and ANID Millennium Science Initiative Program NCN19 170, Chile. MM acknowledges financial support by grant ProyExcel 00505 funded by Junta de Andaluc´ıa and grant PID2021-126348NBI00 funded by MCIN/AEI/10.13039/501100011033/ and ERDF “A way of making Europe”. Appendix A: Analytical solution for the linear model in the long memory case Keeping up to n= 1, that is, considering two polynomials, the solution of Eqs. (80) is O0=−1 eλ2 2"1 1−α1 +λ2ˆτeλ2 2 1 + (1 −α1)eλ2 2ˆτ#,(A1a) O1=λˆτ 1 + (1 −α1)eλ2 2ˆτ,(A1b) P0=Q0/b =λˆτ 1 + (1 −α1)eλ2 2ˆτ,(A1c) P1=Q1/b =−ˆτ 1 + (1 −α1)eλ2 2ˆτ,(A1d) R0=R1= 0.(A1e) where the expression (31) for the bmm′matrix was used. If we truncate up to n= 2, that is, considering three polynomials, the results are O0=− e−λ2 2h(1 −α1)2eλ2λ4+ 2λ2+ 2ˆτ2+ (1 −α1)eλ2 2λ4+ 8λ2+ 6ˆτ+ 4i 2(1 −α1)h(1 −α1)2eλ2ˆτ2+ (1 −α1)eλ2 2(2λ2+ 3) ˆτ+ 2i,(A2a) O1= λˆτh(1 −α1)eλ2 2λ2+ 1ˆτ+ 2i h(1 −α1)2eλ2ˆτ2+ (1 −α1)eλ2 2(2λ2+ 3) ˆτ+ 2i,(A2b) O2=− λ2ˆτh(1 −α1)eλ2 2ˆτ−1i √2h(1 −α1)2eλ2ˆτ2+ (1 −α1)eλ2 2(2λ2+ 3) ˆτ+ 2i,(A2c) P0=Q0/b = λˆτh(1 −α1)eλ2 2λ2+ 1ˆτ+ 2i (1 −α1)2eλ2ˆτ2+ (1 −α1)eλ2 2(2λ2+ 3) ˆτ+ 2,(A3a) P1=Q1/b =− ˆτh(1 −α1)eλ2 22λ2+ 1ˆτ+ 2i h(1 −α1)2eλ2ˆτ2+ (1 −α1)eλ2 2(2λ2+ 3) ˆτ+ 2i,(A3b) P2=Q2/b =√2(1 −α1)eλ2 2λˆτ2 (1 −α1)22eλ2ˆτ2+ (1 −α1)2eλ2 2(2λ2+ 3) ˆτ+ 4,(A3c) 21 R0/b =2(1 −α1)λ2eλ2 2ˆτ2 2 (λ2+ 1) + (1 −α1)eλ2 2(λ2+ 2) λ2ˆτ−2(1 −α1)2eλ2ˆτ2,(A4a) R1/b =− 2λˆτhλ2+ 2 −2(1 −α1)eλ2 2ˆτi 2(1 −α1)2eλ2ˆτ2−(1 −α1)eλ2 2(λ2+ 2) λ2ˆτ−2 (λ2+ 1),(A4b) R2/b = 2√2ˆτhλ2+ 1 −(1 −α1)eλ2 2ˆτi 2(1 −α1)2eλ2ˆτ2−(1 −α1)eλ2 2(λ2+ 2) λ2ˆτ−2 (λ2+ 1).(A4c) Appendix B: Small memory time limit for the transport coefficients in J First, we recall that the operator b Fis not invertible, and its kernel is ϕ=U0. Then, the operator b L′ 0=b F −ν0τ(1 −α1)C(X) =∂2 ∂X2+∂ ∂X A(X)−ν0τ(1 −α1)C(X) (B1) becomes also non invertible and with the same kernel when ν0τ→0. However, in Eq. (79d) for P, the RHS is proportional to U1, which is orthogonal to the kernel, and therefore the solution of this equation is proportional to ν0τ. As a consequence D12 vanishes in this limit by Eq. (84b). In Eq. (79f), γ1is strictly positive making the total operator invertible. As a consequence, Ris proportional to ν0τand the associated transport coefficient µ12 vanishes in the limit of small memory by Eq. (84d). On the other hand, in Eq. (79c) for O, the RHS belongs to the kernel of b F. Therefore, in the limit ν0τ→0 both the LHS and the RHS go to zero simultaneously, resulting in a solution that is of order 1. By Eq. (84a), this implies that D11 remains finite in the limit. Finally, in Eq. (79e) for Q, the RHS is not orthogonal to the kernel of the singular operator, implying that the solution remains finite in the limit, and so it does µ11 by Eq. (84c). Appendix C: Linear response functions Ψρ(k, ω) = k2[γ1µ11 +k2(D22µ11 −D12µ21)τ−i(D12g1+µ11)τω] D11γ1k2+ (−D2 12 +D11D22)k4τ−i[γ1+ (D11 +D22)k2τ]ω−τω2,(C1a) ΨX(k, ω) = τ[−D12k4µ11 + (D11k2−iw)(k2µ21 +ig1ω)] −D11γ1k2+ (D2 12 −D11D22)k4τ+i[γ1+ (D11 +D22)k2τ]ω+τω2.(C1b) [1] H. C. Berg and D. A. Brown, Chemotaxis in Escherichia coli analysed by three-dimensional tracking, Nature 239, 500 (1972). [2] H. C. Berg, E. coli in Motion (Springer Science & Business Media, 2008). [3] E. Korobkova, T. Emonet, J. M. Vilar, T. S. Shimizu, and P. Cluzel, From molecular noise to behavioural variability in a single bacterium, Nature 428, 574 (2004). [4] Y. Tu and G. Grinstein, How white noise generates power-law switching in bacterial flagellar motors, Physical Review Letters 94, 208101 (2005). [5] E. F. Keller and L. A. Segel, Initiation of slime mold aggregation viewed as an instability, Journal of theoretical biology 26, 399 (1970). [6] E. F. Keller and L. A. Segel, Model for chemotaxis, Journal of theoretical biology 30, 225 (1971). [7] M. A. Rivero, R. T. Tranquillo, H. M. Buettner, and D. A. Lauffenburger, Transport models for chemotactic cell populations based on individual cell behavior, Chemical engineering science 44, 2881 (1989). [8] R. M. Ford and P. T. Cummings, On the relationship between cell balance equations for chemotactic cell populations, SIAM Journal on Applied Mathematics 52, 1426 (1992). [9] Y. Tu, T. S. Shimizu, and H. C. Berg, Modeling the chemotactic response of Escherichia coli to time-varying stimuli, Proceedings of the National Academy of Sciences 105, 14855 (2008). [10] Y. Tu, Quantitative modeling of bacterial chemotaxis: signal amplification and accurate adaptation, Annual re- 22 view of biophysics 42, 337 (2013). [11] J. Cremer, T. Honda, Y. Tang, J. Wong-Ng, M. Vergassola, and T. Hwa, Chemotaxis as a navigation strategy to boost range expansion, Nature 575, 658 (2019). [12] W. Alt, Biased random walk models for chemotaxis and related diffusion approximations, Journal of mathematical biology 9, 147 (1980). [13] M. J. Schnitzer, Theory of continuum random walks and application to chemotaxis, Physical Review E 48, 2553 (1993). [14] K. C. Chen, P. T. Cummings, and R. M. Ford, Perturbation expansion of alt’s cell balance equations reduces to segel’s one-dimensional equations for shallow chemoattractant gradients, SIAM Journal on Applied Mathematics 59, 35 (1998). [15] R. Bearon and T. Pedley, Modelling run-and-tumble chemotaxis in a shear flow, Bulletin of mathematical biology 62, 775 (2000). [16] T. Kasyap and D. L. Koch, Instability of an inhomogeneous bacterial suspension subjected to a chemoattractant gradient, Journal of fluid mechanics 741, 619 (2014). [17] F. Tostevin and P. Rein ten Wolde, Mutual information between input and output trajectories of biochemical networks, Physical Review Letters 102, 218101 (2009). [18] G. Lan, P. Sartori, S. Neumann, V. Sourjik, and Y. Tu, The energy–speed–accuracy trade-off in sensory adaptation, Nature Physics 8, 422 (2012). [19] S. Ito and T. Sagawa, Maxwell’s demon in biochemical signal transduction with feedback loop, Nature Communications 6, 7498 (2015). [20] J. Yuan, R. W. Branch, B. G. Hosu, and H. C. Berg, Adaptation at the output of the chemotaxis signalling pathway, Nature 484, 233 (2012). [21] C. Zhang, R. He, R. Zhang, and J. Yuan, Motor adaptive remodeling speeds up bacterial chemotactic adaptation, Biophysical journal 114, 1225 (2018). [22] N. Figueroa-Morales, R. Soto, G. Junot, T. Darnige, C. Douarche, V. A. Martinez, A. Lindner, and E. Cl´ement, 3D spatial exploration by E. coli echoes motor temporal variability, Physical Review X 10, 021004 (2020). [23] N. Figueroa-Morales, A. Rivera, R. Soto, A. Lindner, E. Altshuler, and ´ E. Cl´ement, E. coli “supercontaminates” narrow ducts fostered by broad run-time distribution, Science advances 6, eaay0155 (2020). [24] A. Villa-Torrealba, S. Navia, and R. Soto, Kinetic modeling of the chemotactic process in run-and-tumble bacteria, Physical Review E 107, 034605 (2023). [25] A. Sarkar, G. Georgiou, and M. M. Sharma, Transport of bacteria in porous media: I. an experimental investigation, Biotechnology and Bioengineering 44, 489 (1994). [26] H. Mao, P. S. Cremer, and M. D. Manson, A sensitive, versatile microfluidic assay for bacterial chemotaxis, Proceedings of the National Academy of Sciences 100, 5449 (2003). [27] R. M. Ford and R. W. Harvey, Role of chemotaxis in the transport of bacteria through saturated porous media, Advances in Water Resources 30, 1608 (2007). [28] N. A. Licata, B. Mohari, C. Fuqua, and S. Setayeshgar, Diffusion of bacterial cells in porous media, Biophysical journal 110, 247 (2016). [29] T. Bhattacharjee and S. S. Datta, Bacterial hopping and trapping in porous media, Nature communications 10, 2075 (2019). [30] N. Blackburn, T. Fenchel, and J. Mitchell, Microscale nutrient patches in planktonic habitats shown by chemotactic bacteria, Science 282, 2254 (1998). [31] R. Stocker, Marine microbes see a sea of gradients, science 338, 628 (2012). [32] G. Si, T. Wu, Q. Ouyang, and Y. Tu, Pathway-based mean-field model for escherichia coli chemotaxis, Physical review letters 109, 048101 (2012). [33] G. Si, M. Tang, and X. Yang, A pathway-based meanfield model for e. coli chemotaxis: Mathematical derivation and its hyperbolic and parabolic limits, Multiscale Modeling & Simulation 12, 907 (2014). [34] F. Chalub, Y. Dolak-Struss, P. Markowich, D. Oelz, C. Schmeiser, and A. Soreff, Model hierarchies for cell aggregation by chemotaxis, Mathematical Models and Methods in Applied Sciences 16, 1173 (2006). [35] C. Xue and H. G. Othmer, Multiscale models of taxisdriven patterning in bacterial populations, SIAM Journal on Applied Mathematics 70, 133 (2009). [36] C. Xue, Macroscopic equations for bacterial chemotaxis: integration of detailed biochemistry of cell signaling, Journal of mathematical biology 70, 1 (2015). [37] J. Taktikos, H. Stark, and V. Zaburdaev, How the motility pattern of bacteria affects their dispersal and chemotaxis, PLoS ONE 8, e81936 (2013). [38] Y. S. Dufour, X. Fu, L. Hernandez-Nunez, and T. Emonet, Limits of feedback control in bacterial chemotaxis, PLoS computational biology 10, e1003694 (2014). [39] J. Wong-Ng, A. Melbinger, A. Celani, and M. Vergassola, The role of adaptation in bacterial speed races, PLoS computational biology 12, e1004974 (2016). [40] P. Cluzel, M. Surette, and S. Leibler, An ultrasensitive bacterial motor revealed by monitoring signaling proteins in single cells, Science 287, 1652 (2000). [41] F. Bai, R. W. Branch, D. V. Nicolau Jr, T. Pilizota, B. C. Steel, P. K. Maini, and R. M. Berry, Conformational spread as a mechanism for cooperativity in the bacterial flagellar switch, science 327, 685 (2010). [42] R. Colin, C. Rosazza, A. Vaknin, and V. Sourjik, Multiple sources of slow activity fluctuations in a bacterial chemosensory network, Elife 6, e26796 (2017). [43] J. M. Keegstra, K. Kamino, F. Anquez, M. D. Lazova, T. Emonet, and T. S. Shimizu, Phenotypic diversity and temporal variability in a bacterial signaling network revealed by single-cell fret, Elife 6, e27455 (2017). [44] M. J. Schnitzer, Theory of continuum random walks and application to chemotaxis, Physical Review E 48, 2553 (1993). [45] K. C. Chen, R. M. Ford, and C. P. T., Cell balance equation for chemotactic bacteria with a biphasic tumbling frequency, Journal of Mathematical Biology 47, 518 (2003). [46] E. Lushi, R. E. Goldstein, and M. J. Shelley, Collective chemotactic dynamics in the presence of self-generated fluid flows, Physical Review E 86, 040902 (2012). [47] T. V. Kasyap and D. L. Koch, Chemotaxis driven instability of a confined bacterial suspension, Physical Review Letters 108, 038101 (2012). [48] R. N. Bearon and T. J. Pedley, Modelling run-and-tumble chemotaxis in a shear flow, Bulletin of Mathematical Biology 62, 775 (2000). [49] M. J. Tindall, P. K. Maini, S. L. Porter, and J. P. 23 Armitage, Overview of mathematical approaches used to model bacterial chemotaxis ii: bacterial populations, Bulletin of mathematical biology 70, 1570 (2008). [50] N. Vladimirov, D. Lebiedz, and V. Sourjik, Predicted auxiliary navigation mechanism of peritrichously flagellated chemotactic bacteria, PLoS computational biology 6, e1000717 (2010). [51] J. Saragosti, V. Calvez, N. Bournaveas, B. Perthame, A. Buguin, and P. Silberzan, Directional persistence of chemotactic bacteria in a traveling concentration wave, Proceedings of the National Academy of Sciences 108, 16235 (2011). [52] J. Saragosti, P. Silberzan, and A. Buguin, Modeling E. coli tumbles by rotational diffusion. Implications for chemotaxis, PloS ONE 7, 1 (2012). [53] A. Celani and M. Vergassola, Bacterial strategies for chemotaxis response, Proceedings of the National Academy of Sciences 107, 1391 (2010). [54] J. Long, S. W. Zucker, and T. Emonet, Feedback between motion and sensation provides nonlinear boost in runand-tumble navigation, PLoS computational biology 13, e1005429 (2017). [55] H. Risken, The Fokker-Planck equation: methods of solution and applications (Springer, 1997). [56] G. B. Arfken, H. J. Weber, and F. E. Harris, Mathematical methods for physicists: a comprehensive guide (Academic press, 2011). [57] D. Forster, Hydrodynamic fluctuations, broken symmetry, and correlation functions (CRC Press, 2018). [58] P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics (Cambridge University Press, 1995). [59] S. Chapman and T. G. Cowling, The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases (Cambridge university press, 1990). [60] R. Soto, Kinetic Theory and Transport Phenomena, Oxford Master Series in Physics (Oxford University Press, 2016). [61] N. Brilliantov and T. P¨oschel, Kinetic theory of granular gases (Oxford University Press, 2004). [62] V. Garz´o, Granular gaseous flows (Springer, 2019). [63] A. Villa-Torrealba, C. Ch´avez-Raby, P. de Castro, and R. Soto, Run-and-tumble bacteria slowly approaching the diffusive regime, Physical Review E 101, 062607 (2020). [64] N. Sela and I. Goldhirsch, Hydrodynamic equations for rapid flows of smooth inelastic spheres, to burnett order, Journal of Fluid Mechanics 361, 41 (1998). [65] N. Khalil, V. Garz´o, and A. Santos, Hydrodynamic burnett equations for inelastic maxwell models of granular gases, Physical Review E 89, 052201 (2014). [66] S. de Groot and P. Mazur, Non-equilibrium Thermodynamics, Dover Books on Physics (Dover Publications, 1984). [67] A. Santos, Transport coefficients of d-dimensional inelastic maxwell models, Physica A: Statistical Mechanics and its Applications 321, 442 (2003). [68] J. J. Brey, M. Garc´ıa de Soria, and P. Maynar, Breakdown of hydrodynamics in the inelastic maxwell model of granular gases, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 82, 021303 (2010). [69] S. Khan, S. Jain, G. P. Reid, and D. R. Trentham, The fast tumble signal in bacterial chemotaxis, Biophysical journal 86, 4049 (2004). [70] G. Junot, T. Darnige, A. Lindner, V. A. Martinez, J. Arlt, A. Dawson, W. C. Poon, H. Auradou, and E. Cl´ement, Run-to-tumble variability controls the surface residence times of e. coli bacteria, Physical Review Letters 128, 248101 (2022). [71] N. G. Van Kampen, Stochastic processes in physics and chemistry, Vol. 1 (Elsevier, 1992). [72] Z. Li, Q. Cai, X. Zhang, G. Si, Q. Ouyang, C. Luo, and Y. Tu, Barrier crossing in escherichia coli chemotaxis, Physical review letters 118, 098101 (2017). [73] J. P. Armitage and R. Schmitt, Bacterial chemotaxis: Rhodobacter sphaeroide and sinorhizobium melilotivariations on a theme?, Microbiology 143, 3671 (1997). [74] M. Garren, K. Son, J.-B. Raina, R. Rusconi, F. Menolascina, O. H. Shapiro, J. Tout, D. G. Bourne, J. R. Seymour, and R. Stocker, A bacterial pathogen uses dimethylsulfoniopropionate as a cue to target heatstressed corals, The ISME journal 8, 999 (2014). [75] A. P. Berke, L. Turner, H. C. Berg, and E. Lauga, Hydrodynamic attraction of swimming microorganisms by surfaces, Physical Review Letters 101, 038102 (2008). [76] E. Lauga, W. R. DiLuzio, G. M. Whitesides, and H. A. Stone, Swimming in circles: motion of bacteria near solid boundaries, Biophysical journal 90, 400 (2006). [77] A. M. Hein, D. R. Brumley, F. Carrara, R. Stocker, and S. A. Levin, Physical limits on bacterial navigation in dynamic environments, Journal of The Royal Society Interface 13, 20150844 (2016). [78] T. Jakuszeit, J. Lindsey-Jones, F. J. Peaudecerf, and O. A. Croze, Migration and accumulation of bacteria with chemotaxis and chemokinesis, The European Physical Journal E 44, 1 (2021).