Run-and-tumble bacteria slowly approaching the diffusive regime
Andrea Villa-Torrealba, Cristóbal Chávez Raby, Pablo de Castro, Rodrigo Soto
aa r X i v : . [ c ond - m a t . s o f t ] J un Run-and-tumble bacteria slowly approaching the diffusive regime
Andrea Villa-Torrealba, ∗ Crist´obal Ch´avez Raby, Pablo de Castro, † and Rodrigo Soto ‡ Departamento de F´ısica, Facultad de Ciencias F´ısicas y Matem´aticas,Universidad de Chile, Avenida Blanco Encalada 2008, Santiago, Chile (Dated: June 3, 2020)The run-and-tumble (RT) dynamics followed by bacterial swimmers gives rise first to a ballisticmotion due to their persistence, and later, through consecutive tumbles, to a diffusive process.Here we investigate how long it takes for a dilute swimmer suspension to reach the diffusive regimeas well as what is the amplitude of the deviations from the diffusive dynamics. A linear timedependence of the mean-squared displacement (MSD) is insufficient to characterize diffusion andthus we also focus on the excess kurtosis of the displacement distribution. Four swimming strategiesare considered: (i) the conventional RT model with complete reorientation after tumbling, (ii) thecase of partial reorientation, characterized by a distribution of tumbling angles, (iii) a run-and-reverse model with rotational diffusion, and (iv) a RT particle where the tumbling rate depends onthe stochastic concentration of an internal protein. By analyzing the associated kinetic equationsfor the probability density function and simulating the models, we find that for models (ii), (iii), and(iv) the relaxation to diffusion can take much longer than the mean time between tumble events,evidencing the existence of large tails in the particle displacements. Moreover, the excess kurtosiscan assume large positive values. In model (ii) it is possible for some distributions of tumbling anglesthat the MSD reaches a linear time dependence but, still, the dynamics remains non-Gaussian forlong times. This is also the case in model (iii) for small rotational diffusivity. For all models,the long-time diffusion coefficients are also obtained. The theoretical approach, which relies oneigenvalue and angular Fourier expansions of the van Hove function, is in excellent agreement withthe simulations.
I. INTRODUCTION
There are billions of different species of bacteria onEarth [1]. Because of adaption, their life and swimmingstyles vary across a multitude of distinct environmentsand conditions [2–6]. The vast majority have never beenresearched, and are thus dubbed
Microbial Dark Mat-ter [7]. On the other hand, the
E. coli bacteria continueto be extensively studied. Their motion is usually mod-eled as a run-and-tumble (RT) dynamics. In fact, theirflagella can rotate and propel the cell body in a “run”mode which can suddenly terminate whenever some ofthem reverse direction [8]. This leads to a quick reori-entation mode called “tumble”—which is then followedby another run—with an average tumbling angle of ap-proximately 70 ◦ [9]. In the case of marine bacteria, up to70% of them are thought to have a distribution of tum-bling angles peaked around 180 ◦ instead [10]. Examplesinclude S. putrefaciens and
P. haloplanktis [11], and thusin this case we can speak of a run-and-reverse motion.In his seminal work [12, 13], Berg showed that bacte-ria and other microswimmers performing run-and-tumblemotion develop, in the long term, a diffusive motion.If V is the characteristic run velocity and ν the tum-ble rate (or rotational diffusion coefficient, in the caseof mutant swimmers that tumble only very rarely), thediffusion coefficient scales as D ∼ V /ν , with a pref- ∗ [email protected] † [email protected] ‡ rsoto@dfi.uchile.cl actor that depends on the tumble properties. For ex-ample, in the case of three-dimensional Markovian swim-mers, i.e., each tumble is uncorrelated from previous onesand tumble events are distributed as a Poisson process, D = V / [3 ν (1 − h cos θ s i )], where θ s is the tumbling (or“scattering”) angle between the pre- and post-tumble di-rectors [12, 13]. In the case of E. coli , the data in Ref. [12]gives h cos θ s i ≃ . ν ≃ . − , V ≃ . µ m / s, whichresults in D ≃ µ m / s [14]. The diffusive descriptionof bacterial spreading is extensively used because of itssimplicity, which allows, for example, to couple this ran-dom dynamics with hydrodynamic flows and with thediffusion of nutrients and other chemicals, or to considercomplex geometrical restrictions (for recent applications,see [15–18]). Also, it is possible to include cell divisionand death by employing reaction–diffusion equations, asit is common in chemical and environmental engineeringto describe the spatiotemporal spreading of bacteria [19–21]. Finally, nonlinear effects as a density-dependentdiffusion coefficient are key to describe motility-inducedphase separation [22, 23]. At short times, on the otherhand, the swimmers’ persistent motion gives rise to aballistic motion. Na¨ıvely, the crossover time T cross be-tween the ballistic and diffusive regimes is expected tobe relatively small and to scale as ν − . In this articlewe thoroughly show that, depending on tumbling strate-gies and parameters, the prefactor of this scaling can bequite large and thus the non-diffusive regime can per-sist for long times. This can happen even if the mean-squared displacement (MSD) reaches a linear time de-pendence relatively quickly since having MSD ∼ t is anecessary but not sufficient condition for being in the dif-fusive regime. Note that, associated with T cross , there isa spatial scale L cross = p MSD( T cross ) where the diffusivedescription is not valid. Simulations of active Brownianparticles (ABPs) with large values of L cross show that forlengths of this order a non-diffusive regime indeed arises[16].Our motivation is to quantitatively study the disper-sal process of bacteria. With that purpose in mind, weconsider several run-and-tumble models which are dis-tinct in swimming strategy and compare how slowly thesemicroswimmers approach the diffusive regime. We alsoprovide the spreading dynamics for temporal and spa-tial scales smaller than T cross and L cross , respectively.The swimming strategies considered here are different notonly in terms of the distribution of tumbling angles butalso in whether or not the tumbling rate remains constantover time. In particular we consider the Tu–Grinsteinmodel [24], where the concentration of a phosphorylatedinternal protein named CheY-P changes stochasticallywith time [25], affecting the tumbling rate exponentially.Previous studies have discussed departures from diffusionby using the MSD for run-and-tumble swimmers [26] andthrough the excess kurtosis of the displacement distribu-tion for ABPs [27–29]. More recently, Ref. [30] has stud-ied the non-Gaussian behavior of interacting run-and-tumble particles in the context of active polymer chainsand lattice models, where the authors considered sim-pler tumbling processes and employed analytical meth-ods which are based on solving the associated Langevinequation. In the case of the present work, our analy-sis is done by performing simulations and derivations ofboth the MSD and the excess kurtosis for the differentRT models, aiming to appropriately determine how longthe system takes to reach the diffusive regime. Further-more, the analytical part is carried out from associatedkinetic equations, with Fokker-Planck terms to describerotational diffusion and the evolution of the protein con-centration [31] coupled with a Lorentz term to accountfor the tumbling [32–34]. The simulations are essentiallynumerical implementations of the stochastic rules of mo-tion, i.e., Langevin dynamics. In all cases, we will con-sider two spatial dimensions.The paper is organized as follows. Section II brings ourreview and further development of general theoretical as-pects that will be used throughout the paper. In SectionIII we consider three distinct swimming strategies withconstant tumbling rate. Section IV brings a thoroughanalysis of the case with stochastic tumbling rate. Ourconclusions and a discussion are presented in Section V.Finally, the appendix gives technical details about thesimulations. II. GENERAL THEORETICAL ASPECTS
We start by presenting commonly used model-independent expressions which will be essential in thefollowing sections. From these results we will then derivea general framework to more clearly extract how slowly the diffusive regime is approached. Consider a single bac-terium, initially located at the origin with random orien-tation and internal state. The object of study is ρ ( r , t ),the bacterial density at vector position r at time t ob-tained by averaging over different realizations and initialstates. For this initial condition [ ρ ( r ,
0) = δ ( r )] the bac-terial density is called the van Hove function [35]. TheMSD is h r ( t ) i = Z d r r ρ ( r , t ) . (1)When at long times the diffusive regime is achieved, thedensity obeys ∂ρ∂t = D ∇ ρ, (2)where D is the diffusion coefficient, with solution in twospatial dimensions ρ ( r , t ) = 14 πDt e − r / Dt . (3)Equation (3) implies that h r ( t ) i ∼ t and the diffusioncoefficient is obtained with Einstein’s relation [36], D = lim t →∞ h r ( t ) i t . (4)Calculations become easier to perform through the def-inition of˜ ρ ( k , s ) ≡ Z ∞ d t e − st Z d r e − i k · r ρ ( r , t ) (5)as the Laplace–Fourier transform of ρ ( r , t ), where k is theFourier wave vector and s is the Laplace complex vari-able. Similarly to what is derived in Ref. [35], the secondspatial moment (MSD) and the fourth spatial moment in2D can be calculated, respectively, from h r ( t ) i = L − ( − ∂ ∂k ˜ ρ ( k , s ) (cid:12)(cid:12)(cid:12)(cid:12) k =0 ) (6)and h r ( t ) i = L − ( ∂ ∂k ˜ ρ ( k , s ) (cid:12)(cid:12)(cid:12)(cid:12) k =0 ) , (7)where L − denotes the inverse Laplace transform opera-tor used to bring the result back to the time t domain.The corresponding long-time diffusion coefficient D canbe expressed as [35] D = lim ω → lim k → ω k Re [˜ ρ ( k , iω )] , (8)where ω is real, k ≡ | k | , and Re(˜ ρ ) denotes the real partof ˜ ρ .In the diffusive regime, not only the MSD must growlinearly, but also the displacement distribution must beGaussian. In order to measure the non-Gaussianity ofa particle’s displacement distribution, i.e., a departurefrom the diffusive regime, we will be interested in theexcess kurtosis, defined in 2D by γ ( t ) ≡ h r ih r i − . (9)The excess kurtosis is dimensionless and vanishes for aGaussian distribution of displacements r . For isotropicdistributions, negative values of γ indicate that the distri-bution decays faster than a Gaussian for large displace-ments, while positive values implies that the distribu-tion presents heavy tails. Notice that for one and threedimensions, one would need to subtract 3 and 5/3, re-spectively, instead of 2 in Eq. (9). At short times, whenthe motion is ballistic, h r i equals h r i and, therefore, γ ( t →
0) = −
1, as it is verified for all models presentedin the next sections.At instances where we give explicit expressions for thesecond and fourth moments, we will omit similar expres-sions for the excess kurtosis since they are lengthy andprovide no further information. Nevertheless, expressionsfor their limits as well as their plots will be given and dis-cussed.
A. Extracting the excess kurtosis tail
In general, we will see in the next sections that h r i and h r i approach their asymptotic regimes with exponentialand subdominant polynomial corrections. As a result,the excess kurtosis (9) approaches zero as γ ( t ) ∼ X n a n t − β n e − µ n t , (10)with particular sets of coefficients a n , exponents β n ≥ µ n ≥ µ n and β n are small or zero. From the definition (9), exponentialfactors can come from either the second or fourth mo-ment. For example, for the second moment, we will seein the next sections that h r ( t ) i ∼ Dt + ∞ X n =0 c n t λ n e − µ n t , (11)which in Laplace space gives for small s h e r ( s ) i ∼ Ds + ∞ X n =0 c n λ n !( s + µ n ) λ n (12)and similarly for h e r ( s ) i . Hence, the exponents µ n are recognized as minus the poles of ∂ k ˜ ρ ( k , s ) | k =0 and ∂ k ˜ ρ ( k , s ) | k =0 , and the power exponents λ n are associatedwith pole multiplicity. The slowest decaying mode will be identified as the smallest µ n . For the majority of themodels considered in this article, the long-time behaviorof excess kurtosis can be explicitly obtained in real time.However, for the last model, we will need to extract itfrom Laplace space, as there is no closed expression for γ ( t ). III. CONSTANT TUMBLING RATE
We will now examine three separate limiting cases ofthe well known Markovian run-and-tumble model. Con-sider a particle moving in two spatial dimensions, forwhich tumbling occurs at a constant rate ν . That is,the random walker moves with a constant speed V alonga body-axis ˆn = (cos θ, sin θ ) that can change abruptly ata tumble event, suddenly decorrelating its orientation—in the case of E. coli the duration of the tumble is aboutten times smaller than the duration of the runs [13] andso it is taken as zero here. The new random orientationis chosen with a kernel W ( θ, θ ′ ) that sets the probabilitythat the swimmer changes between two specified orien-tation angles θ and θ ′ at a tumble. We will assume thatthe space is isotropic, hence, the kernel only depends onthe angle difference, W ( θ, θ ′ ) = w ( θ s ), where w is aneven periodic function and θ s ≡ θ ′ − θ is the tumblingangle. In addition to that, the model’s particle is sub-ject to thermal rotational diffusion with coefficient D r .Thus, in the meantime between two consecutive tumblesthe orientation will change slowly and diffusively. The ki-netic equation for the distribution function f = f ( r , θ, t )is [32–34, 37] ∂f∂t + V ˆn ·∇ f = ν Z π w ( θ − θ ′ ) f ( r , θ ′ , t )d θ ′ − ν f + D r ∇ ˆn f, (13)where the distribution function is normalized such that ρ ( r , t ) = R π f ( r , θ, t )d θ . The kernel satisfies R w ( θ )d θ =1, which guarantees that the density ρ is conserved. Wenotice that some of the MSD results in this section arealready present in some form in Refs. [8, 26, 32–34, 37],but they will be developed here either as calibration ofour methodology or to facilitate comparisons against newexpressions such as for the excess kurtosis and with sim-ulations. An entirely new discussion in which the MSDplays only a limited role is provided. A. Conventional run-and-tumble model withcomplete reorientation
We start with the limiting case where D r = 0 and thereis complete reorientation after tumbling, that is, w ( θ s ) =1 / π , the simplest version of the run-and-tumble model.Although no known microswimmer reorients completelyafter a tumble event, this model will serve to calibrateour methodology, as mentioned. In this case the kineticequation for the probability density function f ( r , θ, t ) isjust ∂f∂t + V ˆn · ∇ f = ν π Z π f ( r , θ ′ , t )d θ ′ − ν f. (14)The initial condition is f ( r , θ,
0) = δ ( r ) / π , meaning thatthe initial orientation is random. With a view to ob-taining ˜ ρ ( k , s ) satisfying this kinetic equation, we moveto Laplace–Fourier space. This leads to Eq. (14) beingrewritten as( s + iV k · ˆn + ν ) ˜ f = 12 π (1 + ν ˜ ρ ) , (15)where we used that˜ ρ ( k , s ) = Z π ˜ f ( k , θ, s )d θ (16)for the Laplace–Fourier transform ˜ f ( k , θ, s ) of the distri-bution function. Therefore by isolating ˜ f ( k , θ, s ) and in-tegrating over θ we obtain a closed equation for ˜ ρ , whichgives˜ ρ ( k , s ) = 1 p ( s + ν ) + V k − ν , (17)= 1 s − V k s ( s + ν ) + (3 s + 2 ν ) V k s ( s + ν ) + O ( k ) , (18)where in the second line we made a Taylor expansion in k to easily identify the poles associated to the secondand fourth moments. We can now use the equations inSection II to obtain our desired quantities. The MSD is h r i = 2 V ν (cid:0) ν t + e − ν t − (cid:1) , (19)from which one can either use Einstein’s relation (4) ordirectly apply Eq. (8) to obtain D = V ν , (20)which is a widely known result [38–40]. The fourth spa-tial moment reads h r i = 4 V ν (cid:2) (cid:0) ν t − ν t + 3 (cid:1) + e − ν t ( ν t − (cid:3) , (21)allowing one to compute the excess kurtosis directlythrough (9). The kurtosis longest-standing exponentialgoes as exp( − ν t ), which does not present any singularbehavior.In Fig. 1 the above expressions for the MSD, the diffu-sion coefficient, and the excess kurtosis are tested againstour simulations, which have been performed by directlysolving the associated run-and-tumble motion equations(see the appendix for details on the simulation method).The agreement is excellent as expected since no approx-imations were made. FIG. 1. Conventional run-and-tumble model with completereorientation (to calibrate our methodology): theory (solidblack line) and simulation (circles) for the time evolution ofthe excess kurtosis γ and, in the inset, of the MSD (log-logscale). The dashed line is 4 Dt where the diffusion coefficient D is given by Eq. (20). Units are chosen such that V = ν = 1. B. Partial reorientation
We now generalize the previous analysis to the case ofpartial reorientation while keeping D r = 0. In this casethe kernel is no longer uniformly distributed between 0and 2 π and it is fully characterized by its cosine Fouriercomponents σ n ≡ h cos ( nθ s ) i = Z π − π d θ s w ( θ s ) cos ( nθ s ) , n ≥ , (22)which, in the previous case, vanish completely. Thismodel accounts for many flagellated bacteria and uni-cellular algae [41]. For the case of E. coli , the kernelhas been measured [12], giving σ ≃ .
33 [14]. It canbe shown that the mean-squared displacement dependson σ only [8]. Thus, for the purpose of computing thisquantity, only the average value σ matters and so we donot need to worry about the whole shape of w . However,we show below that the excess kurtosis and the crossovertime to reach the diffusive regime depend also on σ .The Laplace–Fourier transform of the kinetic equa-tion (13) for this case is( s + iV k · ˆn + ν ) ˜ f = 12 π + ν Z π w ( θ − θ ′ ) ˜ f ( k , θ ′ , s )d θ ′ , (23)where we used the same initial condition as in Sec. III A.To solve it, we expand the distribution function in Fouriermodes˜ f ( k , θ, s ) = ∞ X n =0 [ h n cos( nθ ) + g n sin( nθ )] , (24)where the coefficients h n and g n depend on k and s , andare to be determined by plugging the solution into thekinetic equation. Taking k = k ˆx , it is clear that the sinemodes will vanish identically, and so we can set g n = 0from now on. The convolution integral can be expressedas Z π w ( θ − θ ′ ) ˜ f ( k , θ ′ , s )d θ ′ = ∞ X n =0 σ n h n cos( nθ ) . (25)By keeping terms up to n = 2, we truncate the Fourierseries, which allows us to obtain a closed expression for˜ ρ ( k , s ). For the sake of presentation the long result isexpressed as an expansion up to fourth order in k . Thishas no implications as no higher-order derivative in k willbe required. We have˜ ρ ( k , s ) = 1 s − V k s ( ν (1 − σ ) + s )+ (3 s − ν ( σ − V k s ( ν (1 − σ ) + s ) ( ν (1 − σ ) + s ) + O ( k ) . (26)Using the expressions of Sec. II, the MSD is h r i = 2 V ν (1 − σ ) (cid:20) ν t + e − ν (1 − σ ) t − − σ ) (cid:21) , (27)with diffusion coefficient D = V ν (1 − σ ) , (28)which is a well known result [34, 37]. The fourth momentis h r i = 8 V ν (cid:20) ν t (1 − σ ) + e − ν (1 − σ ) t ( σ − σ ) (1 − σ ) − σ + 2 ( σ − σ − σ + 10 σ − − σ ) (1 − σ ) − ν (3 σ − σ − te − ν (1 − σ ) t (1 − σ ) ( σ − σ ) + ν ( σ − σ + 3) t (1 − σ ) ( σ − − (cid:0) σ − σ + 2) σ + 6 σ + 2 σ + 1 (cid:1) e − ν (1 − σ ) t (1 − σ ) ( σ − σ ) (cid:21) . (29)While for the complete-reorientation kernel ofSec. III A there is a single relaxation time, for a gen-eral kernel two relaxation rates appear: ν = ν (1 − σ )and ν = ν (1 − σ ). In this regard, the complete-reorientation case is singular since the two relaxationtimes merge, increasing the multiplicity of the corre-sponding pole in (26). This implies that, while for thecomplete-reorientation case the excess kurtosis decayspurely exponentially as γ ∼ exp( − ν t ), here γ is the sumof two leading terms, exp( − ν t ) /t and exp( − ν t ) /t , ex-cept for the singular case where both rates are equal, inwhich γ ∼ exp( − ν , t ).The approach to a linear time dependence in the MSDis controlled by the relaxation time T = 1 /ν , whichdiverges when the average tumbling angle is small. Nat-urally, in this case, when swimmers deviate little in each FIG. 2. Run-and-tumble model with partial reorientation.(a) Theory (solid black lines) and simulation (circles) for (thenegative of) the excess kurtosis as a function of time in log-log scale. The MSDs are shown in the inset in linear scale tohighlight the departure between the parameters. The valuesof ∆ are indicated in degrees: 164 ◦ , 262 ◦ , and, the complete-reorientation limit, 360 ◦ (green, blue, and red, respectively).(b) MSD in log-log scale for a kernel uniformly distributedaround both 0 ◦ and 180 ◦ with ∆ = 20 ◦ as defined in the maintext, giving σ = 0 and σ ≃ .
99. Theory (solid black line),simulation (circles), and the linear part of the MSD (dashedblue line). Insets: probability distribution function of the x -displacement at different time instants as from simulations(solid lines are normalized Gaussian distributions with thesame mean and variance as the corresponding data). Theexcess kurtosis (not shown) changes from negative to positiveat t ≃ .
76. Units are chosen such that V = ν = 1. tumble event the persistence is enhanced, implying alarge diffusion coefficient. Importantly, also the ampli-tude of the non-diffusive term diverges when σ ≈ − ∆ / , ∆ / T = 1 /ν , appears in thefourth moment given by Eq. (29). Both the relaxationtime and the associated amplitude diverge when σ ≈ γ with amplitudes that can be quite large asthey read (2 T − T ) /t . It can be seen that the smallerthe ∆ the slower is the excess kurtosis approach to zero.Note that, for this kernel, the two relaxation times scaleas 1 / ( ν ∆ ) and are of a similar order, implying thatboth conditions for the diffusive regime to be valid—thelinear increase of the MSD and a small excess kurtosis—are attained in the same timescale.It is possible, however, that T and T decouple if σ ≈ σ is far from 1. Thenalthough the MSD reaches the linear regime rapidly,the excess kurtosis remains finite and positive for longtimes, implying that the diffusion equation is not validin this period. This situation occurs, for example, if thetumbling angles distribution is sharply centered aroundboth 0 and 180 ◦ . Figure 2b shows the case where θ s is uniformly distributed in the ranges [ − ∆ / , ∆ /
4] and[180 ◦ − ∆ / , ◦ + ∆ / θ s is randomly drawn out of these two ranges. Inthis case, a small ∆ leads indeed to a sharp separationof the time scales T and T . As a result, we can see inFig. 2b that the MSD becomes linear in time even if thedisplacement distribution is still strongly non-Gaussian,as revealed by the insets. For this class of kernels, tum-bling gives rise for a single swimmer to a one-dimensionalrandom walk along ˆn and only slowly, with a rate propor-tional to the dispersion of tumbling angles around 0 and180 ◦ , i.e., σ , the process evolves to a two-dimensionaldiffusion. For a collection of swimmers initially seededat r = 0, the intermediate dynamics for T < t < T willtherefore be diffusive only in the radial direction.Also, our analytical results for the MSD in Eq. (27)and for the excess kurtosis [from Eq. (9) using (27) and(29)] are compared against the simulations in Fig. 2. We highlight that they agree well with simulations (even forsmall values of ∆, not shown) despite the approximationmade in truncating the Fourier series up to n = 2. C. Run-and-reverse with thermal rotationaldiffusion
As already mentioned, up to 70% of marine bacteriaare believed to have a distribution of tumbling anglespeaked around 180 ◦ [10]. The soil bacteria Bradyrhizo-bium diazoefficiens has also been shown to perform thiskind of tumbling [42]. In the limiting case known asrun-and-reverse dynamics, which we consider now, theparticle’s tumble can only lead to the exactly oppositemotion direction. In this limit it becomes physically un-reasonable to neglect thermal diffusion and so we willtake D r >
0; otherwise the swimmer will indefinitelyperform a one-dimensional random walk. In this model, w ( θ s ) = δ ( θ s − π ), and hence the kinetic equation reads ∂f∂t + V ˆn · ∇ f = ν f ( r , θ + π, t ) − ν f + D r ∂ f∂θ , (30)where we notice that the indicated instance of f is eval-uated at θ + π , while the other ones are evaluated at θ as per usual. After the Laplace–Fourier transform is ap-plied and using the same initial condition as before, i.e., f ( r , θ,
0) = δ ( r ) / π , we obtain( s + iV k · ˆn + ν ) ˜ f − ν ˜ f ( k , θ + π, s ) − D r ∂ ˜ f∂θ = 12 π . (31)As in the previous case, we expand ˜ f in a Fourier series[Eq. (24)], where again g n = 0 by symmetry. We truncatethe series keeping only the terms n ≤ f ( k , θ, s ) over θ , we obtain˜ ρ ( k , s ) = 4 (4 D r + s ) ( D r + 2 ν + s ) + V k V k D r + 20 D r s + 8 ν s (4 D r + s ) + 16 sD + 3 sV k + 4 s . (32)Upon using the formulae in Section II, we find that the MSD is given by h r i = 2 V ( D r + 2 ν ) h ( D r + 2 ν ) t + e − ( D r +2 ν ) t − i , (33)with diffusion coefficient D = V D r + 2 ν ) , (34)while the fourth moment is h r i = V (cid:20) D − ν − ν D r D ( D r + 2 ν ) + 16 t ( D r + 2 ν ) + 8 ν t − D r tD r ( D r + 2 ν ) + e − D r t D (3 D r − ν ) − e − ( D r +2 ν ) t (cid:0) D (2 ν t + 49) − ν D r (11 ν t + 19) + 15 D t + 12 ν (2 ν t + 3) (cid:1) (3 D r − ν ) ( D r + 2 ν ) . (35)The MSD can rapidly reach a regime where it grows lin-early with time, but for small rotational diffusion, theprocess remains non-Gaussian, with large positive valuesof γ . Similarly to the previous case, when the scatter-ing angle is narrowly distributed around 0 ◦ and 180 ◦ ,reversions at rate ν induce a one-dimensional diffusivemotion along the director axis, but an authentic two-dimensional diffusion is only achieved at a typical time1 / (4 D r ) as the axis changes direction. For a perfect one-dimensional random walk, i.e. for D r = 0, the excesskurtosis equals 1 as γ = h x i / h x i − −
2. Forsmall enough rotational diffusion, the excess kurtosis be-comes positive, having a peak that can be quite large(see Fig. 3), reflecting this quasi one-dimensional motion.Positive excess kurtosis with a very slow decay also ap-pear in the similar case of partial reorientation withoutrotational diffusion (Sec. III B) for σ ≈
1, as can be seendirectly from Eq. (29).
FIG. 3. Run-and-reverse model with rotational diffusion: the-ory (solid black lines) and simulation (circles) for the excesskurtosis γ as a function of time (and for the MSD in the inset).The values of D r are indicated: 0 .
01, 0 .
1, and 1 (green, blue,and red, respectively). Units are chosen such that V = ν = 1.The MSDs are shown in linear scale to highlight the departurebetween the parameters. IV. STOCHASTIC TUMBLING RATE
In bacteria like
E. coli , the tumbling process is trig-gered by a reversion in the sense of rotation (fromcounter-clockwise, CCW, to clockwise, CW) of one orseveral flagella. As a result, the flagella bundle dissem-bles and the propulsion thrust is lost [43]. By analyz-ing the biochemistry of the molecular motor, Tu andGrinstein proposed that the tumbling process can be de-scribed as a two state activated system, where the freeenergy barrier to transit from the CCW to the CW statedepends sensibly on the concentration inside the bacte-rial body of the so-called CheY-P protein, denoted by[ Y ] [24]. In the Tu–Grinstein model the tumble rate is ν = ¯ ν exp( − G ([ Y ]) /k B T ), where G is the free energy bar-rier and ¯ ν a constant. Expanding G around the average value [ Y ], they propose ν ( X ) = ν e αX , (36)where X ( t ) = ([ Y ]( t ) − [ Y ]) /σ Y corresponds to the fluc-tuations in concentration normalized to σ Y , the standarddeviation of [ Y ]. Finally, ν absorbs all the prefactors.Note that ν has been used in the previous sections to de-note the tumbling rate of models without stochasticity,that is, where the tumbling rate is constant over time.Here we use it with exactly the same meaning: in thelimit where α → ν ( X ) → ν . Theparameter α is positive [44] and quantifies the sensitiv-ity of the system to changes in the protein concentration.This phosphorylated protein has a small production rate,with a long memory time T , and consequently X is welldescribed by the Ornstein-Uhlenbeck processd X d t = − XT + r T ξ ( t ) , (37)where ξ is an additive zero-mean Gaussian white noisewith correlation h ξ ( t ) ξ ( t ′ ) i = δ ( t − t ′ ). By tracking sev-eral individual E. coli bacteria it has been possible to fitthe model parameters to T = 19 . ν = 0 .
65 s − , and α = 1 .
62 [45]. The same experiments gave for the ro-tational diffusivity D r = 0 .
025 s − and for the tumbling σ = 0 . D r ≪ ν and that σ ≈ X .With X as a new variable of the distribution function,the kinetic equation for f = f ( r , θ, X, t ) reads ∂f∂t + V ˆn · ∇ f = 1 T (cid:20) ∂ f∂X + ∂ ( Xf ) ∂X (cid:21) + ν ( X )2 π Z π f ( r , θ ′ , X, t )d θ ′ − ν ( X ) f. (38)where the distribution function is normalized such that ρ ( r , t ) = R f ( r , θ, X, t )d θ d X .Once again, we change to the Laplace–Fourier spaceand so Eq. (38) becomes s ˜ f − π ) / e − X / + iV k · ˆn ˜ f =1 T " ∂ ˜ f∂X + ∂ ( X ˜ f ) ∂X + ν ( X ) (cid:20) ˜ g ( k , X, s )2 π − ˜ f (cid:21) , (39)where ˜ f stands for ˜ f ( k , θ, X, s ), and we have made useof the definition˜ g ( k , X, s ) ≡ Z π ˜ f ( k , θ ′ , X, s )d θ ′ (40)and the initial condition f ( r , θ, X, t = 0) = 1(2 π ) / e − X / δ ( r ) , (41)which indicates that the internal variable X is in equilib-rium. We propose the solution˜ f ( k , θ, X, s ) = ∞ X n =0 G n ( X ) ˜ f n ( k , θ, s ) , (42)where the coefficients ˜ f n ( k , θ, s ) do not depend on X and G n ( X ) ≡ e − X / H n ( X/ √ , (43)in which H n is the Hermite polynomial of order n [suchthat H ( x ) = 1, H ( x ) = 2 x , . . . ] [46]. Using the eigen-value equation for the Hermite polynomials allows us towrite1 T " ∂ ˜ f ( k , θ, X, s ) ∂X + ∂ ( X ˜ f ( k , θ, X, s )) ∂X = − T ∞ X n =0 ne − X / H n ( X/ √
2) ˜ f n ( k , θ, s ) . (44)Since our goal is to find the Laplace–Fourier transformof ρ ( r , t ), i.e., ˜ ρ ( k , s ), which does not depend on θ , it ishelpful to define ˜ g n ( k , s ) = R π ˜ f n ( k , θ, s )d θ .At this point we proceed by plugging the above equa-tions into Eq. (39), then multiplying by H m ( X/ √ X . One obtains ∞ X n =0 A mn ( k , θ, s ) ˜ f n ( k , θ, s ) = c m + ∞ X n =0 B mn ˜ g n ( k , s ) , (45)where A mn ( k , θ, s ) ≡ n n ! √ πδ mn (cid:16) s + iV k · ˆn + nT (cid:17) + ν J mn , (46) B mn ≡ ν π J mn , c m ≡ π √ δ m , (47)where the δ ij are Kronecker deltas and J mn ≡ Z ∞−∞ e − y e √ αy H n ( y ) H m ( y )d y. (48)The linear Eqs. (45) can be solved for ˜ f n in terms of˜ g n . Integrating over θ gives now a closed linear set ofequations for ˜ g n , which can be directly solved. Notingthat ˜ ρ = ˜ g (which can be seen through the orthogonalitybetween H and H n ), one obtains˜ ρ ( k , s ) = Z π (cid:20) √ π ( A − ) ( k , θ, s )+ √ π ∞ X m,n =0 ( A − ) m ( k , θ, s ) B mn ˜ g n ( k , s ) (cid:21) d θ. (49)To obtain explicit expressions, Eq. (49) is truncatedat a certain order n = m = N max . The greater the N max the higher is the order of a polynomial in α thatappears in J mn . Hence, increasing N max one increasesthe range in α over which the theory is valid. However,the greater the N max the more complicated are the ele-ments of the inverse of A , which eventually need to beintegrated in θ . Therefore N max also affects how compli-cated it is the ˜ ρ ( k , s ) over which one needs to apply theinverse Laplace transform as well as to compute limits.As it turns out, those complications grow rapidly with N max , with the case N max = 0 being the only one thatwe have treated fully analytically. The ˜ ρ ( k , s ) obtainedby expanding up to this order is identical to the con-ventional RT case (17), provided that one considers thetumbling rate to be ν exp ( α / X . See Section IV A for therelated analysis of the limits T → T → ∞ .For N max = 1 new physics is found. Although in-volved, it is possible to obtain an explicit expression for˜ ρ ( k , s ) from where the diffusion coefficient is obtainedusing Eq. (8), D = V ν e α / (cid:18) α ν T ν T (cid:19) . (50)It is not possible, however, to analytically perform theinverse transforms of the second and fourth moments. In-stead, they are calculated by applying a semi-numericalinverse Laplace transform method for comparison withsimulations. As one can see in Fig. 4, the analytical re-sults in this case agree very well with simulations up toa significant value of α . The higher the α , the higher thepeak in the excess kurtosis.Despite the aforementioned complications in obtainingan expression for ˜ ρ ( k , s ), we can still use the method inSection II A to extract the late-time exponential decayof the excess kurtosis. The second and fourth momentsshare poles, and upon reversing their sign we obtain ν = ν [1 + α (1 / − T ν ) + O ( α )] , (51) ν = ν [1 + 1 / ( T ν ) + α (3 / T ν ) + O ( α )] . (52)While the second rate remains finite for all values of theparameters, ν decays linearly with T . Therefore thegreater the protein memory the longer it will take fordiffusion to be achieved. But the behavior for varying α depends on where in the memory T range we are: if T ν > / ν also decays with α , meaning a slowerapproach to diffusion, and if T ν is smaller than thatthen increasing α speeds up the approach. To evaluatethe importance of this eventual slow approach to diffu-sion, we compute the multiplicity and amplitude of theassociated pole, obtaining h e r ( s ) i ∼ V [1 + ν T (1 + 2 ν T ) α / O ( α )] ν ( s + ν ) (53)This implies that at long times the excess kurtosis expo-nential decay is γ ∼ exp( − ν t ), with an amplitude thatgrows with α , in agreement with the results shown inFig. 4. (a) (b) FIG. 4. Run-and-tumble model with stochastic tumbling rate.(a) Excess kurtosis as a function of time for ν T = 12 . α . The case α = 1 .
62, which corresponds tothe experimentally fitted value for
E. coli is shown in the insetas it falls out of scale. The circles are the simulation resultswhile the solid lines are the small- α theoretical prediction(only for α = 0 and α = 0 . D scaled by the interpolating expression (56) as a function of τ = ( ν T − / ( ν T + 1) for different values of α . The pointsat τ = ± V = ν = 1. A. Zero- and infinite-memory limits
In the limiting case of very small memory time T , X fluctuates rapidly and the tumble rate is effectively anaverage of (36) over all possible values of X , that is, h ν i = ν exp ( α / f forsmall T as ˜ f = ˜ f + T ˜ f + O ( T ) and ˜ g = ˜ g + T ˜ g + O ( T ),and replacing these into the Laplace–Fourier-transformedkinetic equation (39). For O (1 /T ) we obtain a simple dif-ferential equation in X for ˜ f whose solution can be castas ˜ f = e − X / a ( k , θ, s ) where a ( k , θ, s ) is some coeffi- cient function independent of X . At O ( T ) the equationreads iV k · ˆn e − X / a + s e − X / a + ν e − X / αX a − ν e − X / αX b √ π − e − X / (2 π ) / = ∂ ˜ f ∂X + ∂ ( X ˜ f ) ∂X , (54)where b ( k , s ) ≡ R π a ( k , θ, s )d θ . The RHS can be viewedas a differential operator D acting on ˜ f , where the ker-nel of the adjoint operator D † is 1. Thus, upon us-ing the Fredholm Alternative theorem, setting the X -integral of the LHS to zero, one gets the conventionalRT equation (14) with tumbling rate h ν i . Therefore, D T → = V / [2 ν exp ( α / T → ∞ limit is also interesting and roughly corre-sponds to the experimentally fitted values for the E. coli ,for which ν T ≃ .
3. In this case a particle starts with acertain protein concentration (and hence a certain tum-bling rate) as determined by X , which is then kept fixedat all times. The system is therefore equivalent to consid-ering a “polydisperse dilute fluid”, that is, a set of non-interacting particles, where each one has a fixed tumblingrate ν i drawn from a continuous distribution. Thus theaveraged diffusion coefficient is D T →∞ = * V ν i + = V ν Z ∞−∞ e − ( αX + X / d X = V e α / ν , (55)where we notice the opposite sign in the exponential ar-gument in comparison to the T → T and the small- α expansion (50)can be interpolated in a compact expression D i = V ν exp (cid:20) α ( ν T − ν T + 1) (cid:21) . (56)By changing T between its two limits we change τ ≡ ( ν T − / ( ν T + 1) in such a way that τ ∈ [ − ,
1] and,hence we have the bounds D T → ≤ D ≤ D T →∞ . Sim-ulations with different values of α and T show that thisinterpolating expression is good for small α across dis-tinct orders of T (see Fig. 4b). V. CONCLUSIONS
Here we reviewed and extended general theoreticalmethods as well as performed simulations to investi-gate the approach to diffusion of run-and-tumble bacteriawithin four models: conventional run-and-tumble, partialreorientation, run-and-reverse with rotational diffusion,and stochastic tumbling rate. By focusing on the mean-squared displacement and on the excess kurtosis bothanalytically and computationally, we have extracted theeffects of basic model parameters on how slowly diffusionis reached. The methods have been presented in a waythat makes them easy to be translated into other mod-els of particle dispersal. Although we have worked in 2D0for the sake of simplicity, 3D generalizations should bestraightforward to perform [8]. Furthermore, since manytracking experiments are performed in quasi-2D geome-tries [40], our results are directly applicable.For the conventional RT model with complete reorien-tation we obtained that the excess kurtosis approacheszero exponentially with a rate equal to the tumblingrate ν . However, for the other models, new time scalesappear, which can make the approach to the diffusiveregime much slower. For the case of partial reorientation,the new time scales depend on the averages σ = h cos θ s i and σ = h cos 2 θ s i of the tumbling angle θ s , and di-verge when either of them approaches one. This happenswhen the θ s distribution is sharply peaked around both0 and 180 ◦ . For the run-and-reverse model the new timescale is given by the inverse of the rotational diffusiv-ity, D r . When D r ≪ ν , swimmers remain performing aone-dimensional random walk for a long time and transitslowly to the full diffusive motion. Finally, the stochas-tic tumbling rate model, which describes the dynamicsof E. coli , is characterized by two parameters: the sen-sitivity α of the tumbling rate to the concentration fluc-tuations of a relevant protein and the memory time T of this concentration fluctuations. Analytical results areobtained as an expansion for small α , in which case longrelaxation times, eventually diverging, are obtained forlong memory times. Simulations are in excellent agree-ment. In this model we also compute the long-time dif-fusion coefficient, finding an expression valid for small α and any value of T .Concomitantly, when the relaxation times grow, thesame happens with the amplitude of the excess kurtosis,implying that the swimmer dispersion remains largelynon-Gaussian for long times, even though the MSD canalready increase linearly with time. The emergence oflarge relaxation times to reach the vanishing of the excesskurtosis implies that diffusion or reaction–diffusion equa-tions cannot be used to describe bacterial dispersion atintermediate times and distances. Instead, kinetic theoryor discrete element method simulations could be used.This becomes relevant in the design of microrobots forbioengineering applications [47] which include, for exam-ple, killing pathogenous bacteria [48] or removing toxicheavy metals from contaminated water [49].By simulating with the experimentally obtained E. coli values for the partial reorientation model, ν = 1 . − and h cos θ s i ≈ .
33 [9], we estimate that the time toreach an excess kurtosis γ ( t ) such that | γ ( t ) | = 0 .
05 is t ≈ . V , as expected. A similar analysis can be done for themodel with stochastic tumbling rate by using the previ-ously mentioned values T = 19 . ν = 0 .
65 s − , and α = 1 .
62, and by setting D r = 0 and σ = 0 [45]. Thisgives | γ ( t ) | = 0 .
05 at t ≈ . VI. ACKNOWLEDGMENTS
This research is supported by Fondecyt Grant No.1180791 (R.S.) and by the Millennium Nucleus Physicsof Active Mater of ANID (Chile).
APPENDIX: SIMULATIONS
The four models considered in this paper have beensimulated in 2D in the following direct way. First, unitsare chosen such that V = ν = 1. Each particle starts atthe origin r = at t = 0. In the time step ∆ t = 0 . V and varying orientation θ , whose initialvalue is uniformly distributed between 0 and 2 π . Aftereach time step a new orientation is chosen if, and onlyif, a newly drawn random number between 0 and 1 issmaller than ν ∆ t where ν is the tumbling rate. The waythe new θ is chosen as well as the values of ν depend oneach model as follows.For the first three models (Sec. III), a constant ν = ν is taken. In the case of complete re-orientation, any new θ is drawn randomly between 0 and 2 π , just like the ini-tial orientation. In the partial re-orientation model, anew θ is determined from a tumbling angle θ s that isdrawn from one of two distributions: in the first casethe tumbling angle is uniformly distributed in the range[ − ∆ / , ∆ / θ s is uniformly distributed in the ranges [ − ∆ / , ∆ / ◦ − ∆ / , ◦ +∆ / θ by an amount of π at each tumble. Be-sides, as in any simple implementation of rotational dif-fusion, θ changes further at each time step (that is, notonly at each tumble) by the amount √ D r ∆ tη where η is a normally distributed stochastic variable of mean 0and variance 1.In the last model (Sec. IV) we remind that, althougheach new orientation is chosen between 0 and 2 π as inthe complete reorientation model, the tumbling rate ν isno longer constant but rather it follows ν ( X ) = ν e αX .1The initial value of the stochastic variable X is also nor-mally distributed with mean 0 and variance 1. It thenevolves at each time step following a simple Euler-likescheme appropriate to its stochastic differential equa-tion (37), i.e., at each step X changes by the amount − X ∆ t/T + p t/T ξ , with ξ normally distributed withmean 0 and variance 1, again. A large number of timesteps is then performed for each bacterium until a chosen total physical time t is reached. This whole time series isrepeated 2 × times in order to provide the averagedbehavior of the bacteria. In the case of the log-log scaledata presented in Fig. 2a the simulations were repeated2 × , instead. For the numerical values of the diffusioncoefficients, the total simulated time was t = 1000, butvirtually identical results can be obtained with t = 200,which confirms that the MSD decays much more quicklythan the excess kurtosis. [1] Brendan B Larsen, Elizabeth C Miller, Matthew KRhodes, and John J Wiens. Inordinate fondness multi-plied and redistributed: the number of species on Earthand the new pie of life. The Quarterly Review of Biology ,92(3):229–265, 2017.[2] William B Whitman, David C Coleman, and William JWiebe. Prokaryotes: the unseen majority.
Proceedingsof the National Academy of Sciences , 95(12):6578–6583,1998.[3] Stilianos Louca, Florent Mazel, Michael Doebeli, andLaura Wegener Parfrey. A census-based estimate ofEarth’s bacterial and archaeal diversity.
PLoS Biol ,17(2):e3000106, 2019.[4] Oliver Pohl, Marius Hintsche, Zahra Alirezaeizanjani,Maximilian Seyrich, Carsten Beta, and Holger Stark. In-ferring the chemotactic strategy of p. putida and E. coliusing modified Kramers–Moyal coefficients.
PLoS Com-putational Biology , 13(1):e1005329, 2017.[5] Fran¸cois Detcheverry. Generalized run-and-turn mo-tions: From bacteria to L´evy walks.
Physical ReviewE , 96(1):012415, 2017.[6] Maximilian Seyrich, Zahra Alirezaeizanjani, CarstenBeta, and Holger Stark. Statistical parameter inference ofbacterial swimming strategies.
New Journal of Physics ,20(10):103033, 2018.[7] Graham F Hatfull. Dark matter of the biosphere: theamazing world of bacteriophage diversity.
Journal of Vi-rology , 89(16):8107–8110, 2015.[8] Johannes Taktikos, Holger Stark, and Vasily Zaburdaev.How the motility pattern of bacteria affects their disper-sal and chemotaxis.
PLoS ONE , 8(12):e81936, 2013.[9] Howard C Berg.
Random walks in biology . PrincetonUniversity Press, 1993.[10] Jens Efsen Johansen, Jarone Pinhassi, Nicholas Black-burn, Ulla Li Zweifel, and ˚Ake Hagstr¨om. Variability inmotility characteristics among marine bacteria.
AquaticMicrobial Ecology , 28(3):229–237, 2002.[11] Greg M Barbara and James G Mitchell. Bacterialtracking of motile algae.
FEMS Microbiology Ecology ,44(1):79–87, 2003.[12] Howard C Berg and Douglas A Brown. Chemotaxis inEscherichia coli analysed by three-dimensional tracking.
Nature , 239(5374):500–504, 1972.[13] Howard C Berg.
E. coli in Motion . Springer Science &Business Media, 2008.[14] Peter S Lovely and FW Dahlquist. Statistical measures ofbacterial motility and chemotaxis.
Journal of TheoreticalBiology , 50(2):477–496, 1975.[15] L Angelani, R Di Leonardo, and M Paoluzzi. First- passage time of run-and-tumble particles.
The EuropeanPhysical Journal E , 37(7):59, 2014.[16] Lorenzo Caprini, Fabio Cecconi, and Umberto MariniBettolo Marconi. Transport of active particles in anopen-wedge channel.
The Journal of Chemical Physics ,150(14):144903, 2019.[17] Caleb G Wagner, Michael F Hagan, and AparnaBaskaran. Steady-state distributions of ideal activeBrownian particles under confinement and forcing.
Jour-nal of Statistical Mechanics: Theory and Experiment ,2017(4):043203, 2017.[18] Maximilian Seyrich, Andrzej Palugniok, and HolgerStark. Traveling concentration pulses of bacteria in ageneralized Keller–Segel model.
New Journal of Physics ,21(10):103001, 2019.[19] Stephen A Gourley and Yang Kuang. A delay reaction–diffusion model of the spread of bacteriophage infection.
SIAM Journal on Applied Mathematics , 65(2):550?–566,2004.[20] Nick F Britton.
Reaction–Diffusion Equations and TheirApplications to Biology . Academic Press, 1986.[21] Messoud Efendiev.
Evolution Equations Arising in theModelling of Life Sciences . Springer Basel, 2013.[22] J Tailleur and ME Cates. Statistical mechanics of inter-acting run-and-tumble bacteria.
Physical review letters ,100(21):218103, 2008.[23] Michael E Cates and Julien Tailleur. Motility-inducedphase separation.
Annu. Rev. Condens. Matter Phys. ,6(1):219–244, 2015.[24] Yuhai Tu and G Grinstein. How white noise generatespower-law switching in bacterial flagellar motors.
Physi-cal Review Letters , 94(20):208101, 2005.[25] Subrata Dev and Sakuntala Chatterjee. Run-and-tumblemotion with steplike responses to a stochastic input.
Physical Review E , 99(1):012402, 2019.[26] M Reza Shaebani and Heiko Rieger. Transient anoma-lous diffusion in run-and-tumble dynamics.
Frontiers inPhysics , 7:120, 2019.[27] Borge ten Hagen, Sven van Teeffelen, and HartmutL¨owen. Brownian motion of a self-propelled particle.
Journal of Physics: Condensed Matter , 23(19):194119,2011.[28] Xu Zheng, Borge ten Hagen, Andreas Kaiser, Meil-ing Wu, Haihang Cui, Zhanhua Silber-Li, and Hart-mut L¨owen. Non-gaussian statistics for the motion ofself-propelled Janus particles: Experiment versus theory.
Physical Review E , 88(3):032304, 2013.[29] Urna Basu, Satya N Majumdar, Alberto Rosso, andGr´egory Schehr. Active Brownian motion in two dimen- sions. Physical Review E , 98(6):062121, 2018.[30] Stefanie Put, Jonas Berx, and Carlo Vanderzande. Non-gaussian anomalous dynamics in systems of interactingrun-and-tumble particles.
Journal of Statistical Mechan-ics: Theory and Experiment , 2019(12):123205, 2019.[31] K Martens, L Angelani, R Di Leonardo, and L Bocquet.Probability distributions for the run-and-tumble bacte-rial dynamics: An analogy to the Lorentz model.
TheEuropean Physical Journal E , 35(9):84, 2012.[32] David Saintillan. The dilute rheology of swimming sus-pensions: A simple kinetic model.
Experimental Mechan-ics , 50(9):1275–1281, 2010.[33] David Saintillan. Rheology of active fluids.
Annu. Rev.Fluid Mech. , 50:563–592, 2018.[34] Rodrigo Soto.
Kinetic Theory and Transport Phenom-ena . Oxford Master Series in Physics. Oxford UniversityPress, 2016.[35] Jean Pierre Boon and Sidney Yip.
Molecular hydrody-namics . Courier Corporation, 1991.[36] Robert Mazo.
Brownian Motion, Fluctuations, Dynamicsand Applications . Oxford Univ. Press, 2002.[37] Jonathan Saragosti, Pascal Silberzan, and Axel Buguin.Modeling E. coli tumbles by rotational diffusion. Impli-cations for chemotaxis.
PloS ONE , 7(4), 2012.[38] Pawel Romanczuk, Markus B¨ar, Werner Ebeling, Ben-jamin Lindner, and Lutz Schimansky-Geier. ActiveBrownian particles.
The European Physical Journal Spe-cial Topics , 202(1):1–162, 2012.[39] Johannes Taktikos, Vasily Zaburdaev, and Holger Stark.Modeling a self-propelled autochemotactic walker.
Phys-ical Review E , 84(4):041924, 2011.[40] Clemens Bechinger, Roberto Di Leonardo, HartmutL¨owen, Charles Reichhardt, Giorgio Volpe, and GiovanniVolpe. Active particles in complex and crowded environ-ments.
Reviews of Modern Physics , 88(4):045006, 2016.[41] Roman Stocker and William M Durham. Tumbling forstealth?
Science , 325(5939):400–402, 2009.[42] J Ignacio Quelas, M Julia Althabegoiti, Celia Jimenez-Sanchez, Augusto A Melgarejo, Ver´onica I Marconi,El´ıas J Mongiardini, Sebasti´an A Trejo, Florencia Men-gucci, Jos´e-Julio Ortega-Calvo, and An´ıbal R Lodeiro.Swimming performance of Bradyrhizobium diazoefficiensis an emergent property of its two flagellar systems.
Sci-entific Reports , 6:23841, 2016.[43] Xiaobing Chen and Howard C Berg. Torque-speed rela-tionship of the flagellar rotary motor of Escherichia coli.
Biophysical Journal , 78(2):1036–1041, 2000.[44] Philippe Cluzel, Michael Surette, and Stanislas Leibler.An ultrasensitive bacterial motor revealed by mon-itoring signaling proteins in single cells.
Science , 287(5458):1652–1655, 2000.[45] Nuris Figueroa-Morales, Rodrigo Soto, Gaspard Junot,Thierry Darnige, Carine Douarche, Vincent A Martinez,Anke Lindner, and ´Eric Cl´ement. 3d spatial explorationby E. coli echoes motor temporal variability.
PhysicalReview X , 10(2):021004, 2020.[46] George B Arfken and Hans J Weber. Mathematical meth-ods for physicists, 1999.[47] Hakan Ceylan, Joshua Giltinan, Kristen Kozielski, andMetin Sitti. Mobile microrobots for bioengineering ap-plications.
Lab on a Chip , 17(10):1705–1724, 2017.[48] Diana Vilela, Morgan M Stanton, Jemish Parmar,and Samuel S´anchez. Microbots decorated with silvernanoparticles kill bacteria in aqueous media.
ACS Ap-plied Materials & Interfaces , 9(27):22093–22100, 2017.[49] Diana Vilela, Jemish Parmar, Yongfei Zeng, Yanli Zhao,and Samuel S´anchez. Graphene-based microbots for toxicheavy metal removal and recovery from water.
Nano Let-ters , 16(4):2860–2866, 2016.[50] Igor S Aranson and Lev S Tsimring. Pattern formationof microtubules and motors: Inelastic interaction of polarrods.
Physical Review E , 71(5):050901, 2005.[51] David Saintillan and Michael J Shelley. Instabilities andpattern formation in active particle suspensions: kinetictheory and continuum simulations.
Physical Review Let-ters , 100(17):178103, 2008.[52] Aparna Baskaran and M Cristina Marchetti. Hydrody-namics of self-propelled hard rods.
Physical Review E ,77(1):011920, 2008.[53] Rodrigo Soto and Ramin Golestanian. Run-and-tumbledynamics in a crowded environment: Persistent exclusionprocess for swimmers.
Physical Review E , 89(1):012706,2014.[54] Pablo de Castro and Peter Sollich. Phase separation dy-namics of polydisperse colloids: a mean-field lattice-gastheory.
Phys. Chem. Chem. Phys. , 19:22509–22527, 2017.[55] Pablo de Castro and Peter Sollich. Critical phase be-havior in multi-component fluid mixtures: Completescaling analysis.
The Journal of Chemical Physics ,149(20):204902, 2018.[56] Pablo de Castro and Peter Sollich. Phase separationof mixtures after a second quench: composition hetero-geneities.
Soft Matter , 15(45):9287–9299, 2019.[57] Pablo Souza de Castro Melo.
Phase separation of poly-disperse fluids . King’s College London, 2019.[58] Fernando Parisio. Reverse rotations in the circularlydriven motion of a rigid body.
Physical Review E ,78(5):055601, 2008.[59] Pablo de Castro and Fernando Parisio. Role of viscousfriction in the reverse rotation of a disk.