Table Of ContentCollective behavior of quorum-sensing run-and-tumble particles in confinement
Markus Rein, Nike Heinß, Friederike Schmid, and Thomas Speck
Institut fu¨r Physik, Johannes Gutenberg-Universita¨t Mainz, Staudingerweg 7-9, 55128 Mainz, Germany
We study a generic model for quorum-sensing bacteria in circular confinement. Every bacterium
produces signaling molecules, the local concentration of which triggers a response when a certain
threshold is reached. If this response lowers the motility then an aggregation of bacteria occurs,
whichdiffersfundamentallyfromstandardmotility-inducedphaseseparationduetothelong-ranged
nature of the concentration of signal molecules. We analyze this phenomenon analytically and
by numerical simulations employing two different protocols leading to stationary cluster and ring
morphologies, respectively.
PACSnumbers: 05.40.-a,87.17.Jj,87.18.Gh
Introduction. Motilityandlocomotionarebasicaswell bacteriathemotilityisnolongerafunctionofthedensity
as challenging tasks for microorganisms exploring com- but of the local concentration of the autoinducers. The
6 plexaqueousenvironments[1],andnaturehasdeveloped dependence of the local concentration on the sources (or
1 a range of diverse strategies for this purpose. For exam- sinks in the case of catalytic swimmers) is strongly non-
0
ple,thespermcellsofseaurchinsfindtheeggbymoving local and long-ranged, which precludes a mapping onto
2
along helical paths, the curvature of which is controlled equilibrium phase separation.
n by the concentration of a chemoattractant [2, 3]. On InthisLetterwestudyhowpatternscanemergebased
a
the other hand, bacteria might use quorum sensing to onmotilitychangeseveninpopulationswithaconserved
J
respond to changes in their environment [4]. The proba- number of members. To this end, we combine a simple
7
blymostfamousexampleisthemarinebacteriumV. fis- model for bacteria dynamics with quorum sensing. Bac-
] cheri,whichcontrolsbioluminescenceinaccordancewith teria(ormoregenerally, particles)solelyinteractviasig-
t
f its population density. To this end bacteria measure the naling molecules. Particles are confined and we observe
o
localconcentrationofcertainsignalmolecules,calledau- aggregationinthecenteroftheconfinementmediatedby
s
. toinducers, which are emitted by other bacteria. the autoinducers (see Refs. 20, 21 for experiments and
t
a Moving along a chemical gradient of chemoattractants numericalresultsinmorecomplexconfininggeometries).
m or repellents is called chemotaxis. The arguably most This aggregation is in contrast to other collective be-
- famous and best studied model for chemotaxis is the havior of confined active particles like the self-organized
d
Keller-Segel model [5, 6], which consists of two coupled pump in a harmonic trap [22, 23] and the aggregation at
n
partial differential equations, one for the density of dif- walls [24–26]. By combining numerical simulations and
o
c fusing bacteria and one for the concentration of signal analytical theory, we show that the aggregation is deter-
[ molecules. The Keller-Segel model has become a cor- mined by a set of universal parameters that depend on
nerstone to study pattern formation (such as rings and system size.
1
v spots in E. coli [7, 8] and S. typhimurium [9]) and self- Model. We model the bacteria as run-and-tumble par-
2 organization in general [10]. It is not restricted to bacte- ticles moving in two dimensions above a substrate. The
1 ria, e.g., chemotactic behavior has also been reported dynamics mimics straight “runs” due to, e.g., synchro-
5
for self-propelled colloidal particles [11–13], which are nizedflagellainterruptedbyrandom“tumble”events[27,
1
phoretically driven by the catalytic decomposition of, 28]. The equations of motion are
0
. e.g., hydrogen peroxide playing the role of the chemical
1 r˙ =ve +µ F , (1)
signal. k k 0 k
0
6 It has been argued that chemotaxis is not the only where Fk is the force and µ0 the bare mobility. Ev-
1 route to self-organization of motile cells and bacteria, ery particle has an orientation e = (cosϕ ,sinϕ )T
k k k
: and that similar patterns are observed in an arrested along which it is propelled with speed v. This orienta-
v
i motility-induced phase transition in combination with tion remains fixed for an exponentially distributed ran-
X
bacteria reproduction [14, 15]. Such a scenario relies on dom waiting time with mean τ , after which a tumble
r
r a positive feedback through a density-dependent motil- event occurs. We assume the tumbling to occur instan-
a
ity with “slow” bacteria in dense environments [16] so taneously and pick a new, uniformly distributed, orien-
thattheymaymoveagainstadensitygradient. Itallows tation ϕ [29].
k
an effective equilibrium description in terms of a coarse- Every particle produces autoinducers with rate γ.
grained population density as long as the motility is a These signal molecules with concentration cˆ(r,t) diffuse
local function, which seems to be a good assumption for with diffusion coefficient D . While the actual particles
c
short-rangedphysicalinteractions,e.g.,forself-propelled move in two dimensions and are confined by a circular
colloidal particles [17–19]. However, for quorum-sensing confinementwithradiusR(e.g.,duetoasemi-permeable
2
Mean-field theory. We first consider run-and-tumble
particles that only interact via sensing the autoinduc-
ers. Itisthensufficienttoconsidertheone-pointdensity
ψ(r,ϕ,t) of position and orientation, which obeys the
dynamical equation
1 1
∂ ψ =−∇·(veψ)− ψ+ ρ (4)
t τ 2πτ
r r
(cid:82)2π
with particle density ρ(r,t) = dϕ ψ(r,ϕ,t) corre-
0
sponding to the zeroth moment of ψ. The first moment
(cid:82)2π
p(r,t) = dϕ eψ(r,ϕ,t) describes the orientational
0
density of the run-and-tumble particles. From Eq. (4)
we obtain the adiabatic solution p = −1τ ∇(vρ) drop-
2 r
ping the time derivative and neglecting the dependence
on the second moment.
For the analytical treatment we further approximate
FIG. 1: Snapshot of a configuration with N =1000 quorum-
sensing run-and-tumble particles. The color code indicates v = v(c) and τr = τr(c), i.e., the response depends on
thelocalconcentrationcˆofsignalmoleculesfeltbyeachpar- the local average concentration c(r). It is instructive to
ticle. first consider the Keller-Segel model, which, assuming
constant τ , follows in the limit of a weak perturbation
r
of the velocity v(c) ≈ v¯+v(cid:48)(c¯)(c−c¯) around a uniform
membrane),theautoinducerscanpenetratethiswalland
concentration c¯. Eq. (4) implies ∂ ρ = −∇·(vp), which
t
permeate the semispace above the substrate. The time
leads to
evolution of the concentration is thus described by
∂ ρ=∇·(D∇ρ−χρ∇c) (5)
N t
(cid:88)
∂ cˆ=D ∇2cˆ+γ δ(r−r ). (2)
t c k afterinsertingtheadiabaticsolutionfortheorientational
k=1
density. This is the Keller-Segel model together with
Inthefollowingweassumethatthemoleculardiffusionis ∂tc=Dc∇2c+γρδ(z), where the source term in Eq. (2)
muchfasterthanthemotionofthemuchlargerparticles has been replaced by the density ρ3D(r,z) = ρ(r)δ(z).
sothatthereisaninstantaneousstationaryconcentration The two coefficients D = 21τrv¯2 and χ = −12τrv¯v(cid:48)(c¯) are
of autoinducers the effective diffusivity and chemotactic sensitivity, re-
spectively. This result demonstrates that chemotactic
γ (cid:88)N 1 behavior can be achieved simply by changing the magni-
cˆ(r)= . (3)
4πD |r−r | tudeofthespeeddependingonthedifferencebetweenthe
c k
k=1 localandafixedreferenceconcentrationwithoutsensing
the concentration gradient.
Due to the autoinducers eventually leaving the confine-
In the following, however, we are rather interested in
ment the concentration remains finite. The collective
collective effects that arise because of large (discontinu-
behavior of the particles is controlled through τ (r) =
r
ous) changes of the speed due to some particles reaching
τ (cˆ(r)) and v(r) = v(cˆ(r)), which both can vary spa-
r
athreshold,whichisthusbeyondthescopeoftheKeller-
tially through their dependence on the concentration of
Segel model. In the steady state, ∂ ψ =0 and we obtain
autoinducersc(r). Averagingoverparticlepositionsgives t
from Eq. (4)
the average concentration c(r)=(cid:104)cˆ(r)(cid:105).
For the numerical integration of Eq. (1) we employ
τ
a fixed time step δt. After propagating all particles p=− r∇(vρ), 0=∇·(vp), (6)
2
along their orientation, a new random orientation is as-
signed with probability δt/τ . Instead of an ideal hard where we have neglected the dependence on the second
r
wall, we employ the Weeks-Chandler-Andersen (WCA) moment. Forsimplicity, theconfinementisnowmodeled
potential [30] leading to a small but finite “thickness” through the no-flux boundary condition n·p| =0 with
R
d=5·10−3Rofthewall. Particleswithoutward-pointing wall normal n, which ignores the trapping of particles at
orientations remain trapped within the wall until the the wall due to their persistent motion. To compare the
orientations, due to their rotational diffusion, point in- theory with the numerical results, we define the effective
wards again. Fig. 1 shows a snapshot of the system bulk density ρ = N /(πR2) with N the average
0 bulk bulk
for N = 1000 particles after relaxation to the steady number of run-and-tumble particles inside r < R deter-
state [SM]. mined in the simulations.
3
*
FIG. 2: Clustering of slow particles. (a) The particle density ρ(r) for reduced speed α = 0.4 and persistence length (cid:96) = 0.01
obtained numerically for N = 1000 particles (symbols) and from the mean-field theory (dashed line). The solid line is a fit
to Eq. (11). A sketch of a cluster with radius r and density ρ is shown in the inset, where ρ is the density of the dilute
∗ 1 2
region surrounding the cluster. (b) Densities ρ (dense, right branch) and ρ (dilute, left branch) for threshold c¯= 1.45c .
1 2 0
The numerical results slightly depend on the persistence length (cid:96), still the overall agreement with the theoretical prediction is
excellent. (c) Numerical phase diagram of reduced speed α vs. threshold c¯for (cid:96) = 0.01. The color bar indicates the density
difference(ρ −ρ )/ρ ,whereasopencirclesindicatethatnoformationofaclusterhasoccurred. Alsoshownisthetheoretical
1 2 0
prediction (solid lines).
Exploiting the mentioned time-scale separation be- in the simulations are shown for α = 0.4 and N = 1000
tween molecular diffusion and particle motion, the con- particles. It clearly shows a higher inner density corre-
centration of autoinducers follows from Eq. (2) with sponding to the cluster and a lower outer density.
∂ c=0andthesourcetermagainreplacedbythedensity The theoretical density is a step function that can be
t
ρ (r,z) = ρ(r)δ(z), ∇2c(r) = −(γ/D )ρ(r)δ(z), which obtainedasfollows. Theradiusr oftheclusterisdeter-
3D c ∗
isPoisson’sequation. Sincetheconcentrationc(r)deter- mined by the condition c(r )= c¯. Since vρ =const, the
∗
minesthespeedviav(c)itiscoupledwithEq.(6),which density in both regions is constant with dilute density
can now be solved for the density profile ρ(r). In the re- ρ = αρ . Taking into account the conservation of the
2 1
mainder of this Letter we will discuss the two situations bulk density
of one and two thresholds.
Piecewiseconstantspeed. Wenowspecializetoapiece- 1 (cid:90) R
ρ = dr 2πrρ(r), (8)
wise constant speed v(c). Suppose that there are two re- 0 πR2
0
gions with different speeds. Within each region we find
the density of the cluster follows as
∇·p = 0 from Eq. (6) and hence the normal compo-
nents of p across the interface have to be equal. For ρ
ρ = 0 (9)
non-vanishing p this would imply a steady particle cur- 1 α+(1−α)x2
∗
rent, which is excluded by the no-flux boundary condi-
tion. Hence, we conclude that p = 0 and, within our with x = r /R. The radially symmetric solution of the
∗ ∗
theory neglecting fluctuations, vρ=const. Interestingly, Poisson equation reads
the reorientation time τ drops out and does not influ-
r
ence the steady state. This is quite in contrast to self- γ (cid:90) R
c(r)= dr(cid:48) K(r(cid:48)/r)ρ(r(cid:48)) (10)
propelled particles with volume exclusion, the collective πD
c 0
behavior of which is strongly influenced by orientation
relaxation [16, 19]. with kernel K(x) = xK(x2) for x < 1 and K(x) =
Specifically, we introduce a threshold concentration c¯ K(x−2) for x > 1, where K(x) is the complete elliptic
above which the particles slow down by a factor α(cid:54)1, integral of the first kind. Inserting the step profile ρ(r)
we determine self-consistently the radius r of the clus-
(cid:40) ∗
v (cˆ<c¯), ter and thus ρ and ρ for given c¯, α, and ρ . There is
v(cˆ)= 0 (7) 1 2 0
v0α (cˆ(cid:62)c¯) a lower bound c0 < c¯ to the threshold below which no
clustering is possible. It is obtained by considering a ho-
with reference speed v . In the following, the impor- mogenous density ρ = ρ with interface at r = R and
0 1 0 ∗
tance of the directed motion is captured by the persis- thus c =c(R)=γρ R/(πD ).
0 0 c
tence length (cid:96) = v τ /R divided by the radius of the In contrast to the theoretical step profiles, the numer-
0 r
confinement. Due to the slow decay of cˆ(r), we expect ical profiles show a finite interfacial width due to fluctu-
thatthesteadystatewillhavearadialsymmetrywithan ations. Despite the fact that these fluctuations are not
innerdenseclusterandadiluteouterregion. InFig.2(a) accounted for in the theory, the densities of the cluster
thedensityprofilespredictedbythetheoryandmeasured andthediluteouterregionarepredictedquiteaccurately
4
(at least for not too small persistence lengths). The full
profile is well fitted by the empirical expression
ρ +ρ ρ −ρ (cid:18)r−r(cid:48) (cid:19)
ρ(r)= 1 2 − 1 2 tanh ∗ , (11)
2 2 2w
from which the densities ρ and the width w of the
1,2
profile can be extracted. Fitted positions r(cid:48) < r of the
∗ ∗
interface are smaller than the prediction r . In Fig. 2(b)
∗
the densities ρ are plotted as a function of reduced
1,2
speed α for fixed threshold concentration c¯=1.45c and
0
severalvaluesofthepersistencelength,whichshowgood
agreement with the theoretical curve.
The phase diagram in the α–c¯ plane is shown in
Fig. 2(c). Besides the lower threshold c there is also an
0
upper threshold depending on α beyond which no clus-
tering is possible anymore. After a bit of algebra one
finds for the concentration at the interface
c αE(x2)+(1−α)x
∗ = ∗ ∗ (12)
c α+(1−α)x2
0 ∗
with E(x) the complete elliptic integral of the second
kind. This function has a maximum for 0 (cid:54) x (cid:54) 1,
∗
whichmeansthatathresholdhigherthanthismaximum
cannot be reached and, therefore, no clustering is possi-
ble. Again, the theoretical predictions for the parameter
space where clustering is possible agree very well with
the numerical observations as shown in Fig. 2(c).
Rings. More complicated morphologies can also be re- FIG.3: Formationofaringfortwothresholdsc¯2 >c¯1. Shown
alized. To this end we study numerically the effect a is(a)theparticledensitiesρ(r),(b)theconcentrationprofiles
c(r), and (c) the actual average speeds (cid:104)v(cˆ)(cid:105) as functions of
secondthresholdc¯ >c¯ hasontheclustering, wherefor
2 1 the radial distance r/R for three different values of the up-
concentrations cˆ(cid:62) c¯ the particles again move with the
2 per threshold c¯ with the lower threshold c¯ = 1.4c held
2 1 0
higher speed v . The resulting profiles for density ρ(r),
0 fixed. Symbols indicate simulation results for N =1000 par-
concentration c(r), and actual speed (cid:104)v(cˆ)(cid:105) are shown in ticles, solid lines correspond to the mean-field predictions.
Fig.3forafixedfirstthreshold c¯ =1.4c . Ifc¯ liesout- Parameters are: persistence length (cid:96) = 0.01 and reduced
1 0 2
side the region indicated in Fig. 2(c) where clustering is speed α = 0.3. The inset in (b) shows a snapshot of the
possible then basically no change is observed. If the sec- ring (c¯2 = 1.6c0, colors as in Fig. 1), in the inset of (c) the
protocol is sketched.
ond threshold lies within the clustering region the inner
part of the cluster is depleted, leading to the formation
ofastationarydensering[SM].Suchringshavebeenob-
served, e.g., for E. coli [7] but have been attributed to such that the product vρ remains constant. As shown
metabolizing the chemoattractant. in Fig. 3, the mean-field solution correctly captures the
Fig. 3(b) shows the corresponding average concentra- qualitativebehaviorwithaninnerregionwherev >v0α,
tions c(r) of the autoinducers. While these are roughly aring(offinitewidth)withinwhichv =v0α,andasharp
independent of c¯ outside the ring, the concentrations interface to the dilute outer region.
2
saturate at the second threshold c¯ inside the ring. Par- Conclusions. To conclude, we have presented a quan-
2
ticles that locally cross the threshold move faster so that titative theory for the collective behavior of quorum-
there is an effective “pressure” to reduce the inner den- sensing run-and-tumble particles. For one threshold we
sity. Indeed, Fig. 3(c) shows that the measured aver- havederivedspecificexpressionsfortheclustermorphol-
age speed (cid:104)v(cˆ)(cid:105) is higher than αv in the inner region, ogy in circular confinement and we have confirmed these
0
dropping with increasing r. The density maximum co- theoretical predictions in numerical simulations. For a
incides with the minimum of the speed before the speed second threshold we have found the formation of a ring,
again increases going towards the dilute region. While which is also correctly described by the mean-field the-
wedonothaveclosedanalyticalexpressions, wecanstill ory. Patterns typically ascribed to chemotaxis [8, 9] or
solve the mean-field equations numerically, whereby the motility-inducedphaseseparation[14]couldthusalsobe
speed varies between v α and v and the density follows the result of a quorum-sensing mechanism that changes
0 0
5
themotilityofsinglemicroorganismsinresponsetoenvi- St. Angelo, Y. Cao, T. E. Mallouk, P. E. Lammert, and
ronmentalchanges. Incontrastto(effective)equilibrium V. H. Crespi, J. Am. Chem. Soc. 126, 13424 (2004).
phaseseparation,thedensitiesρ ∝ρ areproportional [12] I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert,
1,2 0
and L. Bocquet, Phys. Rev. Lett. 108, 268303 (2012).
to the global density. Moreover, the lower threshold
[13] O. Pohl and H. Stark, Phys. Rev. Lett. 112, 238303
c ∝ R depends on the system size R due to the long-
0 (2014).
ranged concentration profile of the signaling molecules,
[14] M. E. Cates, D. Marenduzzo, I. Pagonabarraga, and
which implies that clustering in a sufficiently large sys- J. Tailleur, Proc. Natl. Acad. Sci. U.S.A. 107, 11715
temissuppressed. Asafirststepwehaveconsideredthe (2010).
mostbasiccombinationofquorumsensingwithasimpli- [15] M. E. Cates, Rep. Prog. Phys. 75, 042601 (2012).
fied model of directed motion. It will be interesting to [16] M.E.CatesandJ.Tailleur,Annu.Rev.Condens.Matter
Phys. 6, 219 (2015).
exploreothermorphologiesandtostudythebasicmech-
[17] I. Buttinoni, J. Bialk´e, F. Ku¨mmel, H. Lo¨wen,
anism for aggregation in more realistic models and test
C.Bechinger,andT.Speck,Phys.Rev.Lett.110,238301
the validity of the scenario we have found.
(2013).
TSacknowledgesfinancialsupportbytheDFGwithin [18] T. Speck, J. Bialk´e, A. M. Menzel, and H. Lo¨wen, Phys.
priorityprogramSPP1726(grantnumberSP1382/3-1). Rev. Lett. 112, 218304 (2014).
We thank ZDV Mainz for computing time on MOGON. [19] T. Speck, A. M. Menzel, J. Bialk´e, and H. Lo¨wen, J.
Chem. Phys. 142, 224109 (2015).
[20] S.Park,P.M.Wolanin,E.A.Yuzbashyan,P.Silberzan,
J. B. Stock, and R. H. Austin, Science 301, 188 (2003).
[21] E.J.Marsden,C.Valeriani,I.Sullivan,M.E.Cates,and
D. Marenduzzo, Soft Matter 10, 157 (2014).
[1] E.LaugaandT.R.Powers,Rep.Prog.Phys.72,096601
(2009). [22] R. W. Nash, R. Adhikari, J. Tailleur, and M. E. Cates,
[2] B. M. Friedrich and F. Ju¨licher, Proc. Natl. Acad. Sci. Phys. Rev. Lett. 104, 258101 (2010).
U.S.A. 104, 13256 (2007). [23] M.Hennes,K.Wolff,andH.Stark,Phys.Rev.Lett.112,
[3] J. F. Jikeli, L. Alvarez, B. M. Friedrich, L. G. Wilson, 238104 (2014).
R.Pascal,R.Colin,M.Pichlo,A.Rennhack,C.Brenker, [24] Y. Fily, A. Baskaran, and M. F. Hagan, Soft Matter 10,
and U. B. Kaupp, Nat. Commun. 6, 7985 (2015). 5609 (2014).
[4] M. B. Miller and B. L. Bassler, Annu. Rev. Microbiol. [25] X.Yang,M.L.Manning,andM.C.Marchetti,SoftMat-
55, 165 (2001). ter 10, 6477 (2014).
[5] E.F.KellerandL.A.Segel,J.Theo.Bio.26,399(1970). [26] F. Smallenburg and H. Lo¨wen, Phys. Rev. E 92, 032304
[6] M. P. Brenner, L. S. Levitov, and E. O. Budrene, Bio- (2015).
phys. J. 74, 1677 (1998). [27] M. Polin, I. Tuval, K. Drescher, J. P. Gollub, and R. E.
[7] J. Adler, Science 153, 708 (1966). Goldstein, Science 325, 487 (2009).
[8] E. O. Budrene and H. C. Berg, Nature 349, 630 (1991). [28] S.Y.Reigh,R.G.Winkler,andG.Gompper,SoftMat-
[9] D. Woodward, R. Tyson, M. Myerscough, J. Murray, ter 8, 4363 (2012).
E. Budrene, and H. Berg, Biophys. J. 68, 2181 (1995). [29] M. J. Schnitzer, Phys. Rev. E 48, 2553 (1993).
[10] T. Hillen and K. Painter, J. Math. Biol. 58, 183 (2009). [30] J.D.Weeks,D.Chandler,andH.C.Andersen,J.Chem.
[11] W.F.Paxton,K.C.Kistler,C.C.Olmeda,A.Sen,S.K. Phys. 54, 5237 (1971).
