Finite response time in stripe formation by bacteria with density-suppressed motility
aa r X i v : . [ q - b i o . CB ] M a y Finite response time in stripe formation by bacteria with density-suppressed motility
Xingyu Zhang ∗ and Namiko Mitarai † The Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100-DK, Denmark. (Dated: May 13, 2019)Genetically engineered bacteria to increase the tumbling frequency of the run-and-tumble motionfor the higher local bacterial density form visible stripe pattern composed of successive high and lowdensity regions on an agar plate. We propose a model that includes a simplified regulatory dynamicsof the tumbling frequency in individual cells to clarify the role of finite response time. We showthat the time-delay due to the response dynamics results in the instability in a homogeneous steadystate allowing a pattern formation. For further understanding, we propose a simplified two-statemodel that allows us to describe the response time dependence of the instability analytically. Weshow that the instability occurs at long wave length as long as the response time is comparable withthe tumbling timescale and the non-linearity of the response function to the change of the densityis high enough. The minimum system size to see the instability grows with the response time τ ,proportional to √ τ in the large delay limit. I. INTRODUCTION
Pattern formation in actively moving and growing liv-ing systems have been attracting interest in the physicscommunity in the last decades [1–3]. Especially, bacterialsystems provide unique opportunities to perform variouswell-controlled experiments thanks to their fast dynamicsand growth as well as the possibility of engineering themolecular network to control their action. One of theinteresting works has been done on the pattern forma-tion by bacterial motion and growth, where a theoreticalanalysis of growing bacteria population with density de-pendent motility [4] predicted various arrested patterns,including concentric ring of high and low bacterial den-sity regions. This was experimentally observed by Liu etal. [5], where
Escherichia coli was engineered to combinethe quorum sensing circuit of
Vibrio fisherei [6] with
E.coli ’s chemotaxis regulatory network [7, 8].
E. coli swims by the run-and-tumble motion, i.e., abacterium ”runs” roughly straight at a constant speed(about 10-20 µ m/sec[9]) between the ”tumbling” wherethe direction of the swimming is shuffled randomly. In awild type E. coli , the tumbling frequency decreases (in-creases) upon the increase (decrease) of attractant (re-pellant), but when a bacterium is left in a constant lig-and environment, the tumbling frequency shows adapta-tion to go back to the original tumbling frequency. Thischange of the tumbling frequency and adaptation enable
E. coli to move towards higher (lower) concentration ofthe attractant (repellant) [10, 11]. The regulation is trig-gered by the binding of the ligand to the receptor. Thebinding affects the phosphorylation of the protein CheY,which interacts with the flagella motor to change thetumbling frequency. Previous research has shown thatthe change of CheY concentration keeps the tumbling ∗ Present address: Department of Biology, Friedrich-Alexander-Universit¨at Erlangen-N¨urnberg, Erlangen, Germany. / Max-Planck-Zentrum f¨ur Physik und Medizin, Erlangen, Germany. † [email protected] duration to be constant and short ( ∼ ∼ Vibrio fisherei [6]. Theaccumulation of the quorum sensing molecules in the en-vironment due to high local density of the bacteria ac-tivates the production of a repressor protein CI, whichrepresses the production of the protein CheZ that is neg-atively correlated with the tumbling frequency. Thus theengineered bacteria tumble more frequently in the denseregion, reducing the diffusion constant of the cells. Whensuch bacteria is grown on a agar plate, the sequential es-tablishment of stripe pattern of the high- and low-celldensity regions were observed [5].This system has been previously modeled by using areaction-diffusion type equations that describes the bac-teria motility through the diffusion and the growth termis added to it [5, 13]. One of the important assumptionsin the model is that the bacterial density ρ obeys thefollowing type of the equations: ∂ρ∂t = ∇ µ ( h ) ρ + g ( ρ ) (1)where µ ( h ) is the quorum sensing molecule density ( h )dependent bacterial motility, and g ( ρ ) is the bacterialgrowth rate. The equation (1) assumes that the diffusionflux of the bacteria J is given by J = −∇ µ ( h ) ρ , with thequotum-sensing molecule dependent motility being insidethe gradient. If in contrast one assumes J = − µ ( h ) ∇ ρ ,the uniform density is the only solution for no flux steadystate J = 0, thus the model would not support the pat-tern. Clearly, in this case, the uniform density state isthe only no-flux steady state without growth ( g ( ρ ) = 0).It is not trivial when J should take the form J = −∇ µ ( h ) ρ . Indeed, the previous works [14, 15] derived thefollowing argument: For simplicity, let’s consider a onedimensional run-and-tumble model, where a bacteriumrun + or − direction with a constant speed v , and itchanges its direction with the frequency α ( h ) (i.e., thetumbling duration is ignored), where h is the signalingmolecule density that the tumbling frequency respondsto. Denoting the population that moves to the + ( − )direction as p + ( p − ), we have ∂p + ∂t = − v ∂p + ∂x − α ( h ) p + + α ( h ) p − , (2) ∂p − ∂t = v ∂p − ∂x + α ( h ) p + − α ( h ) p − . (3)Then the total density ρ = p + + p − obeys ∂ρ∂t = − ∂J∂x , (4) ∂J∂t = = − v ∂ρ∂x − α ( h ) J, (5)with J ≡ v ( p + − p − ). When we assume the relaxation offlux is fast hence ∂J/∂t = 0, we get J = − v / (2 α ( h )) ∇ ρ .Clearly, the no-flux steady state J = 0 supports onlythe uniform solution of a constant ρ , even though theeffective diffusion constant D = v / (2 α ( h )) depends ona parameter h that may vary over space.Of course, the run-and-tumbling model describedabove is a simplification of the reality, and the macro-scopic model of type eq. (1) may be derived from a morerealistic microscopic model. For example, it is well knownthat, if the running velocity v is depends on h , the abovecalculation will be modified to allow the non-uniformsteady state [14, 15]. Another way is to take into ac-count the finite duration of tumbling, which makes h de-pendent tumbling frequency α to produce non-uniformsteady states [16].One of the significant simplification in the describedmodel is to assume that the tumbling frequency wouldchange immediately upon the change of the quorum-sensing molecules [17]. In the engineered bacteria,quorum-sensing molecule indirectly regulate the tran-scription of cheZ gene. The time scale for the changeof the local density of the quorum-sensing molecule re-flected in the CheZ protein level is determined by thelifetime of CheZ protein and intermediate regulator pro-teins CI. When the protein is stable, as often assumedfor the repressor CI [18] unless there is a DNA damagethat activates the CI degradation [19], this response timescale can be as slow as the cell division time scale. Thisis in contrast to the wild type cell’s chemotaxis regula-tion, where the signal transfer is fast because it is doneby changing the protein activity by phospholyration andother modification, and not by changing the protein lev-els [8, 10]. For chemotaxis, many models have been pro-posed and analyzed under externally imposed chemicalgradient [20, 21], and even the self-aggregation via com-bination of quorum-sensing-like chemical and the chemo-taxis has been studied [22]. However, the response timein these situation is rather fast, and the adaptation of the tumbling frequency is the important part of thosemodeling.Recently, a more detailed model of the experimentalsystem [5] has been analyzed by Xue et al. [23]. Afterproposing model equations that take into account intra-cellular signaling including CheZ production and dilutionby growth for each cell, corresponding partial differen-tial equations were derived by assuming the separation oftime scales of intracellular signaling processes, run-and-tumble motions, and growth-induced pattern formation.With using realistic parameter values, the model success-fully reproduced the experimentally observed patterns.Interestingly, the authors analyzed the effect of CheZturnover timescale in the pattern, and reported that thestripe would not form when the time scale of CheZ degra-dation is 10-fold slower than the dilution by cell growth,while 10-fold faster response of CheZ results in sharperpatterns. This work clearly demonstrates the importanceof the response delay, but the relative complexity of themodel makes it difficult to obtain deeper understanding,such as the functional form of the delay-time dependenceof the instability threshold.Motivated by these observations, we propose a simplerun-and-tumble model with density dependent tumblingfrequency with time delay in the response. For simplic-ity, we set the growth term of the bacterium to be zerofor most of the analysis, and focus on the instability ofthe uniform solution, which is necessary for the stripeformation. This model considers the dynamics of the in-ternal state of a cell, and the internal state variable con-trols the tumbling frequency. We show that the modelcan support nonuniform density distribution as long asthe time delay is finite, and we analyze the role of thetime delay in supporting a non-uniform density profile.We then propose a simplified two-state model, where thebacteria can take only high or low tumbling frequency,but the ratio changes depending on the density with atime scale τ . The model allows us to analytically de-termine the condition for the instability of the uniformsolution, revealing that non-zero time delay and strongenough non-linearity can make uniform solution unsta-ble within a certain range of density, but the system sizerequired for the instability increases with the time de-lay. We study how this reflects in the stripe formationwith growth by adding the growth term in the model,and show that the coupling between the protein dilutionrate and the bacterial growth rate ensures that the pat-tern can be formed as long as the non-linearity of theresponse is strong enough. II. RUN-AND-TUMBLE MODEL WITH THEREGULATOR PROTEINA. Microscopic model
First, we consider a microscopic model, where individ-ual cell is moving either left or right at a speed v , and B a c t e r i a ρ / ρ c (b) 0 2 4 6 8 10 12Position x (mm)012 C he m i c a l c FIG. 1. A snapshot of the simulated microscopic model. (a)The bacterial density distribution over space. (b)The hori-zontal axis shows the position, and the vertical axis showsthe concentration of the regulator c , and the color representsthe probability density. changes the direction with the tumbling frequency α ( c i ).Here, c i is the concentration of the tumbling-frequencycontrolling protein in the i -th cell. We assume that thetumbling frequency is an increasing function of c i , with α ( c ) = α l + α h − α l c/c ) − m , (6)where α l and α h are the minimum and the maximumtumbling frequency, respectively, c characterizes thethreshold concentration of the regulatory protein for acell to be in high tumbling mode, and m > c i is transcriptionally regulated. The reg-ulation is through the quorum-sensing system in reality,but here we make a simplification that the active reg-ulator level is proportional to ρ ( t, x i ), where x i is theposition of the bacterium i . Under these assumptions,we model the activation of c i production by the bacterialdensity as dc i dt = k ρ ( t, x ) /ρ c ) − n − c i τ . (7)Here, k is the maximum production rate of c i , and n > ρ c for activation. The degradation timescale τ determines the timescale of the response to the change of the bacterial density. Combined with eq. (6),the bacteria in high density region will tumble more of-ten, resulting in the smaller motility.We have simulated this one-dimensional model. In thefollowing, we set the running velocity v = 15 µ m/s, thetumbling frequencies α h = 11 . α l = 0 . /s . Thisgives the effective diffusion constant D to be within therange from 10 µ m /s to 450 µ m /s. The response timescale τ was set to 400 s , which is reasonably long com-pare to the time scale tumbling. The non-linearlity of theresponses are set by choosing m = 10 and n = 1. Thresh-old protein level c is considered as a concentration unit.We set k = 2 c τ , which allows the maximum proteinlevel to reach twice of the regulation threshold c . Wediscretized the space into meshes of the size ∆ x = 20 µ m,and the local density ρ ( x ) was calculated by countingall the cells in the mesh and divided by mesh size. Thedensity threshold ρ c was set to 1 . /µ m. The model wassimulated by with a constant time step ∆ t = 0 . . × cells in the system size L = 1 . × µ m under the pe-riodic boundary condition. At time zero, the cells aredistributed uniformly in space with the protein level uni-formly distributed between 0 and c max . The snapshotwas taken after 10 s ≈
28 h. The figure clearly showsthe development of locally high density regions and lowdensity regions (Fig. 1a). Figure 1b shows that the cellsin the high density region has high level of the protein c i , resulting in the high tumbling frequency. This clearlydepicts that with the finite response time introduced bythe regulation (7) destabilizes the uniform density state. B. Continuum description of the model
Since the microscopic model is rather time consumingto simulate, we simulate the behavior in a wider range ofparameters by using a continuum version of the model.We consider the the master equation of the distributionfunction f ( t, x, c ), where f ( t, x, c ) dxdc is the probabilityfor bacteria to have regulator concentration between c and c + dc and to be at location between x and x + dx .It is normalized so that N = R ∞−∞ R ∞ f ( t, x, c ) dxdc givesthe total number of the bacteria cells in the system. Wefurther consider the noise in the protein concentration c by interpreting eq. (7) as a rate equation, in order tomake f ( t, x, c ) smooth in the c -space. Using the methodparallel to those described in eqs. (2)-(5) and performingthe system size expansion for the noise in c [24], we obtain ρ ( t, x ) ≡ Z ∞ f ( t, x, c ) dc, (8) ∂f∂t = − ∂J x ∂x − ∂∂c (cid:20) f (cid:18) k ρ/ρ c ) − n − cτ (cid:19)(cid:21) + 12Ω ∂ ∂c (cid:20) f (cid:18) k ρ/ρ c ) − n + cτ (cid:19)(cid:21) , (9) ∂J x ∂t = − v ∂f∂x − α ( c ) J x − ∂∂c (cid:20) J x (cid:18) k ρ/ρ c ) − n − cτ (cid:19)(cid:21) + 12Ω ∂ ∂c (cid:20) J x (cid:18) k ρ/ρ c ) − n + cτ (cid:19)(cid:21) . (10)This equation have a uniform steady state solution fora given average density ρ , which is given by f ss ( c ) = Ae − bc ( a + c ) ab − , (11) a = c max c
11 + ( ρ /ρ c ) − n , (12) b = 2Ω c . (13)Here, A is a normalization constant, defined by R ∞ f ss ( c ) dc = ρ (See eq. (8)). We analyze the sta-bility of this uniform solution by numerically simulatingthe model with an initial condition slightly perturbedfrom the uniform solution by small noise. The formed in-homogeneous pattern (Fig. 2ab)is comparable with thatof the microscopic model. We found that the instabilityrequires a large enough system size, and the critical sizeincrease with the response time scale (Fig. 2c). In thelarge τ limit, the critical system size L c appears to beproportional to √ τ . III. THE SIMPLIFIED TWO-STATE MODEL
The model proposed in the previous section was stilldifficult to tackle analytically. In order to have deeper un-derstanding of the role of the response time, we proposea simplified two-state model that capture the importantfeatures, and perform the stability analysis.Let us consider a run-and-tumble model where a cellcan take either the high tumbling state with the tumblingfrequency α h or the low tumbling state with the tumblingfrequency α l . In each state, a cell will ”run” at a velocity v between the tumbling, and switch the direction withthe given tumbling frequency. We introduce the density-controlled average tumbling frequency by assuming thatthe ratio r ( ρ ) between the low frequency state and thehigh frequency state in the steady state is given by afunction r ( ρ ) = 11 + ( ρ/ρ c ) n , (14) B a c t e r i a ρ / ρ c (b) 0 2 4 6 8 10 12Position x (mm)012 C he m i c a l c τ (s)0.30.60.91.2 L ( mm ) FIG. 2. Snapshot of the continuum model with the initialcondition of uniform density ρ = 1 with small randomness byadding uniformly distributed number in the range (0 , . × − ρ ) to the initial density. (a) The snapshot of bacterialdensity over space at time t s = 30000 s . (b)The horizontal axisshows the position, and vertical axis shows the concentrationof the regulator c , and the color represents the probabilitydensity f ( t s , x, c ). (c) The phase diagram of the instabilityas a function of the time delay τ and the system size L . Thesystem is simulated with the initial condition as in (a), and theinstability was detected by calculating the standard deviation σ of the density profile ρ f ( x ) after long time from the uniformsolution, where σ = L R L ( ρ ( x ) − ρ ) dx . The solid line isproportional to √ τ , which is shown as a guide for the eye. with ρ represents the local density of the bacteria, makingthe average tumbling frequency ¯ α ( ρ ) for a given density ρ to be ¯ α ( ρ ) = r ( ρ ) α l + (1 − r ( ρ )) α h . (15)To introduce the time scale of the state change, τ , we in-troduce the rate for an bacterium to switch from the high(low) tumbling state to the low (high) tumbling state tobe r ( ρ ) /τ [(1 − r ( ρ )) /τ ].By introducing the density of bacteria in the low (high)state ρ l ( ρ h ), we get ρ = ρ h + ρ l , (16) ∂ρ l ∂t = − ∂J l ∂x − τ [(1 − r ( ρ )) ρ l − r ( ρ ) ρ h ] , (17) ∂J l ∂t = − v ∂ρ l ∂x − τ [(1 − r ( ρ )) J l − r ( ρ ) J h ] − α l J l , (18) ∂ρ h ∂t = − ∂J h ∂x + 1 τ [(1 − r ( ρ )) ρ l − r ( ρ ) ρ h ] , (19) ∂J h ∂t = − v ∂ρ h ∂x + 1 τ [(1 − r ( ρ )) J l − r ( ρ ) J h ] − α h J h . (20)Note that, in the limit of infinitely fast response ( τ =0) where ρ l = r ( ρ ) ρ , ρ h = (1 − r ( ρ )) ρ , J l = r ( ρ ) J , and J h = (1 − r ( ρ )) J always hold, eqs. (16)-(20) reduces tothe no time delay case of eqs. (4) and (5) where α ( h )replaced with ¯ α ( ρ ). Clearly the model does not supportthe non-uniform steady state without the time delay.The two-state model (16)-(20) has a uniform steadystate solution ρ l = r ( ρ ) ρ , (21) ρ h = (1 − r ( ρ )) ρ , (22) J h = J l = 0 . (23)We now perform the linear stability analysis of theuniform steady state solution (21)-(23). By setting r = r ( ρ ), we assume ρ l = r ρ + ˜ ρ l , ρ h = (1 − r ) ρ + ˜ ρ h , J h = ˜ J h , J l = ˜ J l , linearize over the perturbation, andlook for the normal modes in the form ˜ ρ l ˜ J l ˜ ρ h ˜ J h = ˆ ρ l ˆ J l ˆ ρ h ˆ J h exp ( ωt − ikx ) . (24)The dispersion relations can be obtained by usingMathematica, as an example of the most unstable modeis shown in Fig. 3. The imaginary part of ω for the mostunstable mode is zero, consistent with the fact that thepattern does not travel. Fig. 3 clearly shows that theinstability happens at the long-wavelength (small k ).By expanding the dispersion relation for the small k ,we get an analytical expression ω ( k ) = − α l τ + 2( α h − α l ) τ ( r + r ′ ρ )2 α l r + 2 α h (1 − r ) + 4 α l α h τ v k + O ( k ) , (25)where r ′ = r ′ ( ρ ) with prime being the derivative by itsargument. The sign of the k term determines the sta-bility in the long wave length. Note that when τ = 0, ω ( k ) = − v α ( ρ ) k + O ( k ), i.e., the system is always lin-early stable in the long wavelength limit and the stabilityis characterized by the diffusion constant for the givendensity v / α ( ρ ). When τ > r ′ <
0, the signcan be positive, i.e., the system has effectively a negativediffusion constant that makes the uniform state unstable. -0.010-0.00500.005 0 1 2 3 4k c E i gen v a l ue ω wavenumber kRe ω Im ω FIG. 3. The eigenvalue ω for the most unstable mode as afunction of the wave number k . The response time τ is setto 400s, and the Hill coefficient n is set to 10. The meandensity ρ is set to be ρ c . Solid line shows the real part,which increase in proportional to k in the small- k limit. Thecritical wave number k c is given by eq. (27). Dotted line showsthe imaginary part of the mode, which stays zero. Especially, at ρ = ρ c , the expression simplifies, andthe condition for the instability can be written as n > α h − α l ) τ + 4 α h /α l − . (26)This shows that the non-linearity of ρ dependence neededfor the instability diverges if we take τ → +0 limit withkeeping the tumbling frequency difference ( α h − α l ) finite,consistent with the lack of instability without the delay.At the same time, with realistic parameters for E. coli , α h − α l is order of 10/sec and α h /α l is order 10, so aslong as the delay τ is order of a few seconds or larger, thehill coefficient n slightly larger than 2 will be sufficientfor instability.The critical wave-number k c > ω ( k c ) = 0, andgiven by k c ( ρ ) = p − [1 + 2 α l τ + 2( α h − α l ) τ ( r + r ′ ρ )] τ v , (27)which is a positive real number when the system is lin-early unstable. Eq. (27) shows that the smallest systemsize L c = 2 π/k c necessary to see the instability increaseswith τ , and scales with √ τ in the large τ limit. This isconsistent with the observation in the continuum modelas shown in Fig. 2c.At ρ = ρ c , we have L c ( ρ c ) = 2 πτ v p − [1 + 2 α l τ + ( α h − α l ) τ (2 − n ) / , (28)showing that the critical system size also dependent onthe non-linearity of the response n .We confirmed that the two state model simulated withthe parameters where the uniform solution is linearly un- B a c t e r i a den s i t y ρ / ρ c Position x (mm) 0 0.5 1 1.5 2 2.5 3 3.5 -40 -30 -20 -10 0 10 20 30 40(b) B a c t e r i a den s i t y ρ / ρ c Position x (mm)
FIG. 4. (a) Numerical simulation of two state model with α h = 11 . s − , α l = 0 . s − , v = 15 µms − , τ = 4 × s .Initial condition was ρ l = r ( ρ ) ρ , ρ h = (1 − r ( ρ )) ρ , J l = J h = 0. Small randomness was added to the ini-tial density by adding uniformly distributed number in therange (0 , . × − ρ ) to ρ l and ρ f . (b)Two state modelwith growth. α h = 11 . s − , α l = 0 . s − , v = 15 µms − , τ = 15min, and ρ max = 3 ρ c are used. The growth rate g is set to 1 /τ . Non-linearity of the response n is set to 10.For the initial condition, a small population of bacteria wasplaced in the center at t = 0 as follows: ρ ( x ) = d √ π e − x w with d = 0 . ρ c , w = 0 . ρ l ( x ) = r ( ρ ( x )) ρ ( x ), ρ h ( x ) =(1 − r ( ρ ( x )) ρ ( x ), J h ( x ) = J l ( x ) = 0. The snapshot wastaken at t = 1 × s ∼
28 hours. Note that the stripe patterndisappears in the middle due to the logistic term that allowsgrowth for low density region, as discussed in ref. [13]. stable shows a stripe pattern as shown in Fig. 4a. We ob-serve the formation of the high and low density regions,similar to the pattern observed in the original complexmodel (Fig. 2).Finally, we checked if the two-state model can give a stripe pattern in growing bacterial front, when thetime delay is comparable with the bacterial division time,which would be the situation where the time delay iscoming from the dilution of the stable proteins. For thispurpose, we simulate two-state model with the growth byadding a logistic growth term gρ l (1 − ρ/ρ max ) to eq. (17)and gρ h (1 − ρ/ρ max ) to eq. (19), respectively, with agrowth rate g (cf. [13]). We set g = 1 /τ and simulatedthe expansion of the bacterial population in one dimen-sion. The result shown in Fig. 4(b) confirmed that thesystem can obtain the stripe formation at the front witha high enough non-linearity. IV. SUMMARY AND DISCUSSION
We have shown that adding a response time delayto a simple run-and-tumble model with bacterial den-sity dependent tumbling frequency causes the instabilityin the uniform density profile. The instability is long-wave length, and the critical system size L c for instabil-ity grows with the time delay as L c ∝ √ τ in the large τ limit. These observations are reproduced in analyticallytractable two state model, where the explicit parameterdependence of L c was obtained. For the two state model,it was found that the non-linearity n required for the in-stability diverges in the τ → +0 limit. This is consistentwith the previous finding of lack of instability in run-and-tumble model with τ = 0. Finally, the two state modelwas extended to include the cellular growth, and it wasconfirmed that the stripe formation is possible even ifthe time delay τ is comparable with the generation time1 /g . The stripe formation was seen in the growth front,showing that instability is there for this relatively longtime delay.In biological systems, the typical timescales for dif-ferent types of regulations can vary significantly, andsometimes significant delay in response can happen. Thepresent work demonstrates the importance of taking intoaccount those delays in proper modeling of observed phe-nomena. ACKNOWLEDGEMENT
This work was funded by the Danish National ResearchFoundation (BASP: DNRF120). XZ and NM thank Lei-han Tang for fruitful discussions. [1] H. Fujikawa and M. Matsushita, Journal of the physicalsociety of japan , 3875 (1989).[2] E. Ben-Jacob, I. Cohen, and H. Levine, Advances inPhysics , 395 (2000).[3] J. A. Shapiro, Annual Reviews in Microbiology , 81(1998). [4] M. Cates, D. Marenduzzo, I. Pagonabarraga, andJ. Tailleur, Proceedings of the National Academy of Sci-ences , 11715 (2010).[5] C. Liu, X. Fu, L. Liu, X. Ren, C. K. Chau, S. Li, L. Xiang,H. Zeng, G. Chen, L.-H. Tang, et al. , Science , 238(2011). [6] C. M. Waters and B. L. Bassler, Annu. Rev. Cell Dev.Biol. , 319 (2005).[7] H. C. Berg, D. A. Brown, et al. , Nature , 500 (1972).[8] H. C. Berg, E. coli in Motion (Springer Science & Busi-ness Media, 2008).[9] U. Alon, L. Camarena, M. G. Surette, B. A. y Arcas,Y. Liu, S. Leibler, and J. B. Stock, The EMBO journal , 4238 (1998).[10] N. Barkai and S. Leibler, Nature , 913 (1997).[11] U. Alon, M. G. Surette, N. Barkai, and S. Leibler, Nature , 168 (1999).[12] A. Bren, M. Welch, Y. Blat, and M. Eisenbach, Pro-ceedings of the National Academy of Sciences , 10090(1996).[13] X. Fu, L.-H. Tang, C. Liu, J.-D. Huang, T. Hwa, andP. Lenz, Physical review letters , 198102 (2012).[14] M. J. Schnitzer, Physical Review E , 2553 (1993).[15] J. Tailleur and M. Cates, Physical review letters ,218103 (2008). [16] A. Curatolo, Collective behaviours in living systems :From bacteria to molecular motors , Ph.D. thesis, Uni-versite Paris Diderot (2017).[17] A. L. Andersen,