Table Of Content,
Arrested phase separation in reproducing bacteria: a generic route to pattern
formation?
M. E. Cates,1 D. Marenduzzo,1 I. Pagonabarraga,2 and J. Tailleur1
1SUPA, School of Physics and Astronomy, University of Edinburgh, Mayfield Road, Edinburgh EH9 3JZ, UK
2Departament de F´ısica Fonamental, Universitat de Barcelona - Carrer Mart´ı Franqu´es 1, 08028-Barcelona, Spain
(Dated: January 4, 2010)
We present a generic mechanism by which reproducing microorganisms, with a diffusivity that
depends on the local population density, can form stable patterns. It is known that a decrease of
swimming speed with density can promote separation into bulk phases of two coexisting densities;
thisisopposedbythelogisticlawforbirthanddeathwhichallowsonlyasingleuniformdensityto
be stable. The result of this contest is an arrested nonequilibrium phase separation in which dense
droplets or rings become separated by less dense regions, with a characteristic steady-state length
0
scale. Cell division mainly occurs in the dilute regions and cell death in the dense ones, with a
1
0 continuous flux between these sustained by the diffusivity gradient. We formulate a mathematical
2 modelofthisinacaseinvolvingrun-and-tumblebacteria,andmakeconnectionswithawiderclassof
mechanismsfordensity-dependentmotility. Nochemotaxisisassumedinthemodel,yetitpredicts
n
the formation of patterns strikingly similar to those believed to result from chemotactic behavior.
a
J
4 Microbialandcellularcoloniesareamongthesimplest is diffusivity-independent.) Other non-equilibrium ef-
examples of self-assembly in living organisms. In na- fects, such as ratchet physics, have also been observed
]
E ture, bacteria are often found in concentrated biofilms, and used either to rectify the density of bacteria [16–18]
P matsorothercolonytypes,whichcangrowintospectac- or to extract work from bacterial assemblies [19].
. ular patterns visible under the microscope [1, 2]. Also Some aspects of bacterial patterning show features
o
i in the laboratory, bacteria such as E. coli and S. ty- common to other nonequilibrium systems, and a cru-
b phimirium form regular geometric patterns when they cial task is to identify the key ingredients that control
q- reproduceandgrowonaPetridishcontainingagelsuch their development. In many equilibrium and nonequi-
[ as agar. These patterns range from simple concentric librium phase transitions an initial instability creates
rings to elaborate ordered or amorphous arrangements density inhomogeneities; these coarsen, leading eventu-
1
of dots [3–6]. Their formation results from collective be- ally to macroscopic phase separation [20]. The situa-
v
7 haviourdrivenbyinteractionsbetweenthebacteria,such tion observed in bacterial assemblies often differs from
5 as chemotactic aggregation [6], competition for food [8] this; long-lived patterns emerge with fixed characteristic
0 or changes in phenotypes according to density [11]. The length scales, suggesting that any underlying phase sep-
0 question as to whether general mechanisms lie behind arationissomehowarrested. Thestrongdiversityofbio-
1. this diversity of microscopic pathways to patterning re- logicalfunctionsmetinexperimentshasledtoanequally
0 mains open. diverse range of proposed phenomenological models [5–
0 Unlike the self-assembly of colloidal particles, pattern 7, 9–11] to account for such effects. Most of them rely
1 formationinmotilemicroorganismsandotherlivingmat- on the coupling of bacteria with external fields (food,
:
v ter is typically driven by non-equilibrium rather than chemoattractant, stimulant, etc.), and many involve a
i thermodynamic forces. Indeed, the dynamics of both di- large number of parameters due to the complexity of the
X
lute and concentrated bacterial fluids is already known specific situation of interest. The most common mech-
ar tobevastlydifferentfromthatofasuspensionsofBrow- anism used to explain the bacterial patterns is chemo-
nian particles. For instance, suspensions of active, self- taxis [6]: the propensity of bacteria to swim up/down
propelled particles, have been predicted to exhibit gi- gradients of chemoattractants/repellants. This explana-
ant density fluctuations [12, 13], which have been ob- tionissowellestablishedintheliteratureforatleasttwo
served experimentally [14]. Similarly, an initially uni- organisms (E. coli and S. typhimuridium [6]) that obser-
form suspension of self-propelled particles performing a vation of similar patterns in other species might defensi-
“run-and-tumble” motion like E. coli has recently been bly be taken as evidence for a chemotactic phenotype.
showntheoreticallytoseparateintoabacteria-richanda Here we identify a very general mechanism that can
bacteria-poor phase, provided that the swimming speed leadtopatternformationinbacterialcoloniesandwhich
decreases sufficiently rapidly with density [15]. This is may encompass a large class of experimental situations.
akin to what happens in the spinodal decomposition of This mechanism involves a density-dependent motility,
binaryimmisciblefluids,buthasnocounterpartinasys- giving rise to a phase separation which is then arrested,
tem of Brownian particles interacting solely by density- on a well defined characteristic length scale, by the birth
dependent diffusivity. (The latter obey the fluctuation- and death dynamics of bacteria. For definiteness we
dissipation theorem, ensuring that the equilibrium state will work this through for a particular model of bacte-
2
rialrun-and-tumblemotion, involvingaswimspeedthat ever, a non-uniform swimming speed v(r) also results in
depends (via unspecified interactions) on local bacterial a mean drift velocity V = −vτ∇v [21] which here gives
density. This gives pattern similar to those observed V = −D(cid:48)(ρ)∇ρ/2 [15]. This contribution is crucial to
in experiments [5, 6]. However, the basic mechanism phase separation [15] and will again play a major role
—density-dependent motility coupled to logistic popu- here. However this term has no counterpart in ordinary
lationgrowth—isnotlimitedtothisexample. Ourwork Brownianmotion(evenifparticleshavevariablediffusiv-
demonstratesthatchemotaxisper seisnotaprerequisite ity) and was accordingly overlooked in previous studies
forobservingwhataresometimescolloquiallyreferredto which relied on phenomenological equations involving a
as ‘chemotactic patterns’. density-dependent diffusivity and no drift [8].
It is indeed remarkable that density-dependent swim Coupling the diffusion-drift equation for run-and-
speed and logistic growth alone are sufficient to create tumble bacteria, as derived in [15], with the logistic
some of the specific pattern types previously identified growth term, the full dynamics is then given by:
with specific chemotaxis mechanisms. In mechanistic
∂ρ(r,t)
terms, we find that the logistic growth dynamics effec-
= ∇·[D (ρ)∇ρ(r,t)] (1)
tively arrests a spinodal phase separation that is known ∂t e
(cid:18) (cid:19)
tofollowfromadensity-dependentswimspeed [15]. Put ρ(r,t)
+αρ(r,t) 1− −κ∇4ρ(r,t)
differently,aninitiallyuniformbacterialpopulationwith ρ
0
small fluctuations will aggregate into droplets, but these
willnotcoarsenfurtheronceacharacteristiclengthscale where the ‘effective diffusivity’ is
isachieved,atwhichaggregationandbirth/deatheffects
1
come into balance. Starting instead from a small inocu- D (ρ)=D(ρ)+ ρD(cid:48)(ρ) (2)
lum, we predict formation of concentric rings which, un- e 2
dersomeconditions,atleastpartiallybreakupintospots
This results from the summed effects of the true diffu-
at late times [5].
sive flux −D(ρ)∇ρ and the non-linear drift flux ρV. In
Toexemplifyourgenericmechanismwewillstartfrom
Eq.(1) we have also introduced a phenomenological sur-
a minimal model of run-and-tumble bacteria, that can
facetensionκ>0,whichcontrolsgradientsinthebacte-
run in straight lines with a swim speed v and randomly
rialdensity. Suchacontributionhasbeenshowntoarise
changedirectionataconstanttumblingrateτ−1 [21,22].
when the speed of a bacterium depends on the average
To this we add our two key ingredients: a local density-
density in a small local region around it, rather than a
dependent motility, and the birth/death of bacteria, the
strictlyinfinitesimalone[15]. Eq.(1)neglectsnoise,both
latter accounted for through a logistic growth model.
in the run-and-tumble dynamics and in the birth/death
Of course, bacteria can interact locally in various ways,
process. The former noise source conserves density and
ranging from steric collisions [15] to chemical quorum-
should become irrelevant at the experimental time scale
sensing [6]. (Indeed a nonspecific dependence of motility
ofdays. Ontheotherhand,thenon-conservativenoisein
on bacterial density was previously argued to be central
the birth and death dynamics may be more important,
to bacterial patterning by Kawasaki [8].) Here we fo-
and we have verified that our results are robust to its
cus on the net effect of all such interactions on the swim
introductionatsmalltomoderatelevels. Numericalsim-
speed v(ρ), which we assume to decrease with density
ulations of Eq. (1) have been performed with standard
ρ. This dependence might include the local effect of a
finite difference methods (although noise does require
secreted chemoattractant (such as aspartate [3–5] which
careful treatment, as in [23, 24]), with periodic bound-
causesaggregation,effectivelydecreasingv)butdoesnot
ary conditions used throughout. For definiteness, all our
assume one. simulations have been carried out with v(ρ)=v e−λρ/2,
0
In addition to their run-and-tumble motion, real bac- where v >0 is the swim speed of an isolated bacterium
0
teria continuously reproduce, at a medium-dependent and λ > 0 controls the decay of velocity with density.
growthratewhichrangesfromaboutonereciprocalhour The precise form of v(ρ) is however not crucial for the
in favourable environments such as Luria broth to sev- phenomenology presented here, and the instability anal-
eral orders of magnitude lower for ‘minimal’ media such ysis offered below does not assume it.
as M9. In bacterial colonies patterns may evolve on The logistic population dynamics alone would cause
timescalesofdays[6],overwhichsuchpopulationgrowth the bacterial density to evolve toward a uniform den-
dynamics can be important. sity, ρ(r) = ρ , which constitutes a fixed point for the
0
We now derive continuum equations for the local den- proposed model. Although this homogeneous configu-
sity ρ(r,t) in a population of run-and-tumble bacteria, ration is stable in the absence of bacterial interactions,
with swim speed v(ρ), growing at a rate α(1 − ρ/ρ ). it has been shown [15] that, without logistic growth, a
0
The latter represents a sum of birth and death terms, in density-dependent swim speed v(ρ) leads to phase sepa-
balance only at ρ=ρ . At large scales in a uniform sys- rationviaaspinodalinstabilitywheneverdv/dρ<−v/ρ.
0
tem,themotionofindividualbacteriaischaracterizedby By Eq. (2) this equates to the condition D < 0, and it
e
adiffusivityD(ρ)=v(ρ)2τ/d,whereτ−1 isthetumbling is indeed obvious that the diffusive part of Eq. (1) is
rate and d the dimensionality [21, 22]. Crucially how- unstable for negative D . It is important, clearly, that
e
3
D can be negative although D is not. This holds for a
e
much wider class of nonequilibrium models than the one
studied here; we return to this point at the end of the
paper.
For the choice of v(ρ) made in our simulations, we
have D = D(ρ)[1−ρλ/2] and the flat profile will thus
e
become unstable for ρ above 2/λ. We have confirmed
0
this numerically, and find that upon increasing ρ , the
0
uniform state becomes (linearly) unstable, evolving in a
1D geometry into a series of “bands” of high bacterial
density separated by low density regions. Depending on
theparameters,thistransitioncanbecontinuous(super- √
Figure2: ThreeplotsofΛ (q)for|D (ρ )|/ ακ=1;2;3. At
critical), with the onset of a harmonic profile whose am- q e 0
the transition, only one critical mode q=q is unstable.
plitude grows smoothly with ρ , or discontinuous (sub- c
0
critical) with strongly anharmonic profiles (see Fig. 1).
The flat profile ρ = ρ is thus stable if Λ ≤ 0 for all
0 q
q and is unstable otherwise. From the expression for
D (ρ ), Eq. (2), one sees that instability occurs if
e 0
ρ D(cid:48)(ρ ) D (ρ )
Φ≡− 0 0 ≥1 and − √e 0 ≥2 (4)
2D(ρ ) ακ
0
At the onset of the instability only one mode is unsta-
(cid:112)
ble, with wavevector q = 2α/|D (ρ )|, as can be seen
c e 0
in Figure 2. The first condition in Eq. (4), Φ ≥ 1, is
Figure 1: Growth of the instability in the supercritical (left)
equivalent to the requirement that D < 0 given pre-
and subcritical cases (right). The three lines correspond to e
viously. From the dispersion relation, Eq. (3), we see
three successive times. A small perturbation around ρ (red
0 that the resulting destabilization is balanced by the sta-
line) growth toward harmonic or anharmonic patterns in the
bilizing actions of bacterial reproduction and the surface
supercritical or subcritical case, respectively. Left: Super-
criticalcase(α=κ=0.01,λ=0.02,ρ =15,D =1,times: tension at large and small wavelength, respectively. The
0 0
102, 103, 104). Right: Subcritical case (α = κ = 0.005, unstablemodesthusliewithinabandq1 <q <q2 where
(cid:112) (cid:112)
λ=0.02, ρ =11, D =1, times: 3.102, 3.103, 105) q (cid:39) q ≡ α/|D (ρ )| and q (cid:39) q ≡ |D (ρ )|/κ set
0 0 1 α e 0 2 κ e 0
the wavelengths below and above which the stabilizing
The transition to pattern formation arising from effects of bacteria reproduction and the surface tension
Eq.(1)isafullynonequilibriumone: itisnotpossibleto can compete with the destabilizing effect of the negative
writedownaneffectivethermodynamicfreeenergywhich diffusivity,respectively. Forunstablemodestoexist,one
would lead to this equation of motion. Nonetheless, it is needs q1 ≤ q2; restoring prefactors, this yields 2qα ≤ qκ
possibletounderstandwhythebirth/deathprocesseffec- which is the second criterion in Eq. (4). This analysis is
tively arrests the spinodal decomposition induced by the consistentwiththeviewthatphaseseparationisarrested
density-dependent swim speed. The latter tends to sep- by the birth/death dynamics, which stabilizes the long
aratethesystemintohighandlowdensitydomainswith wavelength modes (Λ0 = −α), while the phenomenolog-
densities either side of ρ . (Without the logistic term, ical tension parameter κ primarily fixes the interfacial
0
these would coarsen with time.) Bacteria thus tend to structure of the domains, not their separation.
be born in the low density regions and to die in the high We now consider more closely the parameters control-
densityregions. Tomaintainasteadystate,theyhaveto ling the transition to pattern formation. For definite-
travel from one to the other: balancing the birth/death ness,weaddressthespecificcaseusedforoursimulations,
terms by the diffusion-drift transport flux between the D(ρ)=D(0)exp(−λρ). To put Eq. (1) in dimensionless
domains then sets a typical scale beyond which domain form, we define rescaled time, space and density as
coarsening can no longer progress. Were any domain to (cid:16)α(cid:17)1/4 ρ
becomemuchlarger,thedensityatitscentrewouldsoon t˜=αt; ˜r= r; u= (5)
κ ρ
regress towards ρ , re-triggering the spinodal instability 0
0
locally. (Thisiscloselyreminiscentofwhathappensina The equation of motion now reads
thermodynamic phase separation when the supersatura- u˙ =∇·[Re−2Φu(1−Φu)∇u]+u(1−u)−∇4u (6)
tion is continuously ramped [25].) √
To better understand the onset of the instability, let where R ≡ D0/ ακ and Φ = λρ0/2 are the two re-
uslinearizeEq.(1)aroundρ(r)=ρ andworkinFourier maining dimensionless control parameters. Meanwhile
0
space. Defining ρ(r)=ρ +(cid:80) δρ exp(iq·r) yields: the conditions (4) for pattern formation become
0 q q
exp(2Φ)
δ˙ρq =Λqδq; Λq =−α−q2De(ρ0)−κq4 (3) Φ≥1; R≥Rc =2 Φ−1 (7)
4
These relations, combined with the preceding linear sta- available become less reliable [26]. We emphasize again
bilityanalysis,defineaphasediagraminthe(R,Φ)plane thatthebasicmechanismforpatterningpresentedabove
(Figure3)thatagreesremarkablywellwithnumericalre- does not depend on the precise form chosen for v(ρ).
sults for systems prepared in a (slightly noisy) uniform Quantitativelyhowever,Eq.(8),andthefrontierbetween
initial state. subcriticality and supercriticality, do depend on the de-
tails of the interplay between the nonlinearity in v(ρ)
and the logistic growth term. We leave further analysis
of such model-specific features to future work.
While the amplitude equation is more easily devel-
oped in 1D, the stability analysis offered above is valid
in higher dimensions and it is natural to ask what hap-
pens in 2D, which is the relevant geometry for Petri dish
studieswithgrowingbacterialcolonies. Fig. 4showsthe
simulated time evolution of ρ(r,t) for a system started
with small random fluctuations around the equilibrium
density ρ , with other parameters as in Fig. 1. Per-
0
haps not surprisingly, bands are replaced by droplets of
the high density phase dispersed in a low density back-
ground at large times. This is the typical steady state
obtained with a near-uniform starting condition. How-
ever,thestructureandorganizationofthebacterialdrops
inthesteadystatedependsonthepoint(R,Φ)chosenin
thephasediagram. Generally,thecloserthesystemisto
the supercritical instability curve, the more ordered the
patterns. For instance we have observed an essentially
crystalline distribution of bacterial drops, which devel-
ops defects and eventually becomes amorphous on mov-
ingfurtherawayfromthephaseboundary(Fig. 3,insets
tomainpanel). Forparticularchoicesofparameters,our
model can also admit other steady state patterns. Close
Figure 3: Top: Phase diagram in the (R,Φ) plane. The to the supercritical line, where the phase transition is
outer region corresponds to stable behavior whereas within continuous, we can obtain long-lived stripes, whereas for
thecurve, patterningoccurs. Thesolidlineisthetheoretical fixed large values of R and Φ close to the (right) sub-
phaseboundary–Eq(7)–whichaccuratelyfitsthenumerics critical phase boundary, we have also observed ‘inverted
(black squares). The blue and red sections correspond to droplets’ with a high density lawn punctuated by low
continuous and discontinuous transitions respectively. The density ‘holes’.
two magenta dots correspond to two 2D simulations which
In these 2D geometries initialized from a near-uniform
show ordered harmonic patterns close to supercriticality and
state, droplets can coalesce in the early stages, while
amorphous patterns otherwise. Bottom-left: Transition in
the supercritical regime. The blues lines correspond to the at late times the dynamics is governed by evaporation-
theory–Eq(8)–whereasthesquarescomesfromsimulations condensation events (see Figs. 4 and 5). However it is
(Φ = 1.5; 1.35; 1.2 from top to bottom). Bottom-right: already apparent from Fig. 4 that coarsening eventually
TransitioninthesubcriticalregimeforΦ=1.06andΦ=1.7 stops and the droplets reach rather well-defined steady
(bottom to top). state sizes and centre-to-centre distances. This can be
quantified by looking at the time evolution of the char-
Closetothetransition, theemergentsteady-statepat- acteristic domain size, L(t), which we have computed as
tern can be studied using an amplitude equation (see the inverse of first moment (times 2π) of the structure
Supporting Information). Introducing ε=(R−R )/R , factor [20]. Fig. 5 suggests that L(t) at late times even-
c c
one gets in 1D that for 1.08≤Φ≤1.58, the transition is tually stops increasing and reaches a steady state value.
supercritical (continuous) and the steady state is given (The visible steps in domain size mark discrete evapo-
by ration events involving smaller bacterial droplets; pre-
sumably L(t) would become smooth for a large enough
u(cid:39)1+A(ε)cos(x); system.)
18(1−Φ)2 (8) These droplet patterns in steady state are very simi-
A2(ε)=ε lar to those observed for E. coli in a liquid medium or
34Φ4−56Φ3−24Φ2+31Φ+19
S. typhimurium in semi-solid agar (0.24% water-agar in
which agrees with simulations (Fig 3, bottom left). Out- Ref. [5]) when starting from a uniform distribution [6].
side this range, the transition becomes subcritical (dis- For the E. coli case, interactions are believed to come
continuous, Fig 3, bottom right) and the analytical tools from chemoattractant, emitted by the bacteria them-
5
Figure 4: Numerical results for a 2D simulation with size equal to 16×16, λ=0.3, D =1, α=0.01, κ=0.001 and ρ =10.
0 0
Times corresponding to the snapshots are (in simulation units, from left to right): 5, 12, 19.6 and 50.6.
our simulations with a single small droplet of high den-
sityρ,wefindthatasimilarlypatternedbacterialcolony
structure develops. First, the bacteria spread radially
(through a Fisher-like wave), forming an unstructured
lawn with the highest density at the center. This back-
ground density increases logistically until the onset of
instability via our generic phase-separation mechanism;
with circular symmetry, the instability causes concentric
ringsofhighbacterialdensitytosuccessivelydevelopthat
are very stable in time (Figs. 6c and 6d). The patterns
observed at later times again depend on the position of
the parameters in the (R,Φ) plane. If we fix a value of
R, e.g. 100, larger values of Φ in the unstable region
Figure 5: Plot of the characteristic domain size, L(t) as a lead to rings being very stable. For smaller values of Φ,
function of time for a system in the inhomogeneous phase, on the other hand, effectively corresponding to weaker
withinitiallyrandomdensityfluctuationsaroundρ0. Param- interactions between the bacteria, we observe that rings
eters were: α = 0.01, λ = 0.27, κ = 0.001, while the system
initiallyformbutrapidlydestabilizethroughasecondary
size was 20×20. The solid line corresponds to a single run,
modulation of the bacterial density along them. This
while the dashed line is an average over 6 runs. The steps in
eventually breaks the rings into a series of drops. The
the single run curve correspond to evaporation-condensation
innerringsdestabilizefirst, andthesystemevolveseven-
events, highlighted by black squares in the snapshots shown
tually to the same steady state as found starting from
in the figure (before and after one of the steps respectively,
arrowsindicatepositionsontheplotcorrespondingtothetwo a uniform density, composed of drops with well defined
snapshots). characteristic size and separation. All this phenomenol-
ogy is strikingly reminiscent of the dynamics observed
by Woodward et al [5] for S. typhimirium, where rings
are stable at large concentrations of potassium succinate
selves, that is not degraded over time [3, 4]. The
(a ‘stimulant’ which promotes pattern formation), but
chemoattractant distribution should become more and
break up into drops at smaller ones. Our model shows a
more uniform so that these interactions decay to zero
similarmorphologicalchangewhendecreasingΦ,i.e. the
as time proceeds. In our framework this is analagous
strength of the interactions.
to decreasing Φ, which will turn any initially unsta-
ble state into a homogeneous one, and can thus explain Different views are possible concerning the ability of
thatthepatternsobservedexperimentallyfadewithtime our generic model to reproduce the observed chemotac-
(whereas in our simulations Φ remains constant and the tic patterns of E. coli and S. typhimurium [6]. One pos-
pattern are stable indefinitely). E. coli in a semi-solid sibility is that Eq. (1), with the interpretation we have
mediumalsoexhibitsdropletpatternsofhighsymmetry. given for it, actually does embody the important physics
Inourframework,suchpatternsresultfromacontinuous of pattern formation in these organisms. Indeed it is
transition, close to the supercritical line. well accepted that bacteria in the high density concen-
The growth of bacterial colonies of S. typhimurium tric rings are essentially non-motile [27]. The precise
starting from a small inoculum of bacterial cells in semi- mechanism leading to this observation is not well under-
solidagarleadstoquitespecific(transientbutlong-lived) stood [6], but it is possible that the chemotactic mecha-
patterns, with the bacteria accumulating in concentric nismmainlyactstoswitchoffmotilityathighdensity. If
rings that can subsequently fragment into a pattern of so,byfocussingsolelyonthisaspect(withacorrespond-
dots [5, 6]. Once again, although these patterns are be- ingly vast reduction in the parameter space from that
lieved to stem from a chemotactic mechanism [6], we of explicit chemotactic models [5, 6]) our model might
find they can arise in principle without one, so long as capture the physics of these chemotactic patterns in a
our two basic ingredients of density-suppressed motility highly economical way. Interestingly, our model is es-
and logistic growth are both present. Indeed, initializing sentially local, whereas chemotaxis in principle mediates
6
Figure 6: Dynamics of formation of ”chemotactic patterns” in 2D, starting from a single small bacterial droplet in the middle
of the simulation sample. Top row: Formation of ”chemotactic rings” in a system with α=0.1, λ=0.33, and κ=0.001. The
simulation box has size 40×40. The snapshots correspond to times equal to (in simulation units, from left to right) 10, 50,
100 and 270. Bottom row: Breakage of ”rings” into ”dots”. The four snapshots correspond to the time evolution of a system
with α=0.1, λ=0.26, and κ=0.001. We show a 40×40 fraction of the simulation box, with the boundaries far away and
not affecting the pattern. The snapshots correspond to times equal to (in simulation units, from left to right) 10, 70, 290 and
1220.
interactions between bacteria that are nonlocal in both with D a many-body diffusivity and µ the correspond-
space and time. It is not clear whether such nonlocal- ing mobility. This alone ensures that no drift velocity
ity is essential for the chemotactic models in [5, 6] or can arise purely from gradients in D. In nonequilbrium
if fast-variables approximations and gradient expansions systems, one should expect generically to find such drift
would reduce these models (which invole between 9 and velocities, and the run-and-tumble model is merely one
12 parameters) into Eq.[1]. In this case, we would still instance of this. Accordingly one can expect in princi-
have in Eq.[1] a highly economical model for chemotac- ple to find cases of negative D in other microorganisms
e
ticpattern-formationorganisms,possiblywithadifferent showing distinctly different forms of density-dependent
interpretation of D and κ. self-propulsion.
e
Alternatively, the success of our local model for these To summarize, we have studied the dynamics of a sys-
chemotactic organisms might be largely coincidental. tem of reproducing and interacting run-and-tumble bac-
But in that case, such a sparse model should be easily teria, in the case where interactions lead to a decreasing
falsified, for instance by using the linear stability anal- local swim speed with increasing local density. We have
ysis to relate the typical length scale of the patterns thereby identified a potentially generic mechanism for
to microbial parameters. This length scale is of order pattern formation in which an instability towards phase
(cid:112)
2π/q = 2π |D |/α, with |D | (cid:39) D, a typical bacterial separation, caused by the tendency for bacteria to move
c e e
diffusioncoefficient(D ∼O(100µm2s−1)forE.coli[22]). slowly where they are numerous, is arrested by the birth
Using the previously quoted growth rate α(cid:39)1 hr−1, we and death dynamics of bacterial populations. We have
get a ring separation of ∼ 1 mm, in order-of-magnitude shown that these two ingredients alone are enough to
agreement with the experimental value [5]. This test capturemanyofthepatternsobservedexperimentallyin
couldperhapsbesharpenedusefullybyalteringthegrow- bacterial colonies – including some that have only pre-
ing medium so as to change α. viously been explained using far more complex models
involving specific chemotactic mechanisms. Indeed, if
More generally, our analysis of Eq. (1) shows that the
motility decreases steeply enough with density, then a
main prerequisite for pattern formation, assuming the
spatially homogeneous bacterial population becomes un-
presence the logistic growth term, is negativity of the
stable to density fluctuations leading to the formation of
effective diffusion constant D . For run-and-tumble dy-
e bands (1D) or droplets (2D). The length scale of the re-
namics, D < 0 was shown to arise for a sufficiently
e sultingpatternissetbyabalancebetweendiffusion-drift
strong decay of swim-speed with density; it does so be-
fluxes and the logistic relaxation of the population den-
cause spatial variations in the true diffusivity D(ρ) cre-
sitytowardsitsfixed-pointvalue. Startinginsteadfroma
ate a drift flux ρV = −ρD(cid:48)(ρ)∇ρ/2 which can overcom-
smallinitialdropletofbacteria,wepredicttheformation
pensate the true diffusive flux −D∇ρ [15]. Negative D
e ofconcentric rings, eachofwhich mayeventuallyfurther
could, however, equally arise for any density dependent
separate into droplets.
nonequilibrium diffusion process. Indeed, the principle
ofdetailedbalance,whichholdsonlyforequilibriumsys- Inseveralwellstudiedsystems,suchcharacteristicpat-
tems, leads to the Einstein relation, that D = k Tµ terns are (with good reason) believed to be the direct
B
7
result of chemotactic behavior [5, 6]. It is therefore re- for ε>0 and small. The dynamics now reads
markable that they can also arise purely from the inter-
playofdensity-dependentdiffusivityandlogisticgrowth, w˙ =Lw−2ε∂2w+g(w) (A5)
x
without explicit reference to the dynamics (or even the
presence) of a chemoattractant. This suggests that sim- where L=−(1+∂2)2 is the linear part of the evolution
x
ilar patterns might arise in organisms having no true operator at the transition, 2ε∂2w gives an extra linear
x
chemotactic behavior at all. Such patterns could then part due to the perturbation (ε > 0) and g(w) is the
be the result of local chemical signalling without gradi- non-linear part:
ent detection (quorum-sensing, not chemotaxis) or even
purely physical interactions (steric hindrance), either of g(w)=−w2 −∂ (cid:104)2(1+ε)(cid:0)(1+ w )e−2w−1(cid:1)∂ w(cid:105)
whichcouldinprincipleproducetherequireddependence Φ x Φ−1 x
of motility on density. Last, a motility decreasing with (A6)
density is only one of the many mechanism that could
lead to D(cid:48)(ρ) < 0 and our analysis would apply equally
e
to all such cases. 1. Amplitude equation
The simplest version of our model allows identifica-
tion of just two dimensionless parameters that control As usual with the amplitude equation approach, we
theentirepattern-formingprocess. Inbothhomogeneous expandwinpowerseriesoftheperturbationεandstudy
and centrosymmetric geometries, this gives predictions Eq. (A5) order by order. As shown below (Eqs (21-23)),
forhowthepatterntypedependsoninteractionstrength the correct expansion is
which are broadly confirmed by experimental data. This
suggests that some of the diverse patterns formed by w =U ε1/2+U ε+U ε3/2+... (A7)
0 1 2
colonies of motile bacteria could have a relatively uni-
versal origin. Expanding (A6) to the order ε3/2 and substituting
in (A5) yields order by order:
Acknowledgments −LU = 0 (A8)
0
U2 3−2Φ
We thank Otti Croze for discussions. We acknowledge −LU1 = − Φ0 − Φ−1 ∂x2U02 (A9)
funding from EPSRC EP/E030173. MEC holds a Royal
2U U 3−2Φ
Society Research Professorship. IP acknowledges the −LU = −2∂2U − 0 1 − 2∂2U U(A10)
2 x 0 Φ Φ−1 x 0 1
Spanish MICINN for financial support (FIS2008-04386).
4Φ−2
− ∂2U3
3Φ−1 x 0
Appendix A
Equation (A8) can be easily solved and yields
We show here how to derive the amplitude equa- U =Aeix+A∗e−ix (A11)
0
tion (8). Let us start from the dimensionless equation
of motion (6)
Theamplitudeoftheperturbationwearetryingtoderive
is thus 2|A|. Equation (A9) can also be solved directly:
u˙ =∇[Re−2Φu(1−Φu)∇u]+u(1−u)−∇4u (A1)
U =Beix+B∗e−ix+C+De2ix+D∗e−2ix (A12)
and recall the two conditions for patterning Eq. (7): 1
where B can de determined from higher order equations
Φ>1; Rexp(−2Φ)(Φ−1)>2 (A2)
(but does not interest us here), and C and D are given
by
To analyze precisely the transition, we derive below the
steady-state limit of the amplitude equation in 1D. By
2|A|2 A2 3−2Φ 1
inspection one sees that the unperturbed steady-state C =− ; D = (4 − ) (A13)
Φ 9 Φ−1 Φ
of (A1) is given by u = 1. To characterize the ampli-
tude of the perturbation around u = 1, we introduce
as can be checked by direct substitution in Eq. (A9).
u=1+w/Φ so that w evolves with
Equation (A10) does not always have a solution. In-
w˙ =−∂x[Re−2Φ(Φ−1)(1+Φw−1)e−2w∂xw]−w(1+wΦ)−∂x4w dyieeeldd,athmeulatpipplleicoaftieoinx,o(fsiLncetoLeainxy=fu0ncatniodnLUi2s lcinaenanro)t.
(A3) The r.h.s. however does contain a multiple of eix whose
We are interested by the vicinity of the transition where prefactor must thus vanish. This gives a condition for
the expansion to provide a proper steady-state solution
Re−2Φ(Φ−1)=2(1+ε) (A4) of the problem. Let us summarize the contributions of
8
√
the different terms to the prefactor of eix in the r.h.s of ε. One could look for a more general expansion:
Eq. (A10)
−2∂2U yields 2A (A14) w =U εα+U ε2α+U ε3α (A20)
x 0 0 1 2
3−2Φ 2 (cid:16)3−2Φ 1(cid:17)
( 2∂ − )U U yields 2A|A|2 − (A15)
Φ−1 x Φ 0 1 Φ−1 Φ
In this case, the expansion of equation (A5) yields two
×(cid:16)1(43−2Φ − 1)− 2(cid:17) power series: (cid:80)R εαk and (cid:80)W εαk+1. For the two se-
k k
9 Φ−1 Φ Φ
riestogivetermsthatcanbalanceeach-other, oneneeds
4Φ−2 Φ−2
− ∂2U3 yields 4 A|A|2 (A16) α+1=kα for k ≥2 and thus
3Φ−1 x 0 Φ−1
The sum of these terms vanishes only if 1
α= (A21)
(cid:16) (cid:17) k−1
A 9Φ2(Φ−1)2+2|A|2(34Φ4−56Φ3−24Φ2+31Φ+19) =0
(A17)
Thecandidatesforαarethus1; 1/2; 1/3; .... Notethat
and thus either
α ≤ 1 implies 2α+1 ≥ 3α. We can therefore stop the
9Φ2(1−Φ)2 expansion at 3α and 2α+1 to get the first three terms
A=0; or |A|2 =
2(34Φ4−56Φ3−24Φ2+31Φ+19)in the expansion of equation (A5)
(A18)
Let us first try α = 1. The order by order the expan-
Finally, the first order in the amplitude equation yields
sion yields
√
w(x)=2|A| εcos(x−x ) (A19)
0
L2U = 0, O(ε) (A22)
where x is a constant. Note that by construction 0
|A|2 > 00 and a non-zero solution only exists for Φ ∈ L2U = −U02 − 3−2Φ∂2U2−2∂2U , O(2ε()A23)
[1.08439,1.59237]. For these values of Φ, Eq. (A18) 1 Φ Φ−1 x 0 x 0
and (A19) work very well, as can be checked in figure
7. Outside this range the transition becomes subcritical
Equation(A22)yieldsU =Aeikx+A∗e−ikxbutequation
and the standard approach does not work anymore. Al- 0
(A23)cannotbesolvedsincethereisanon-zeromultiple
ternative treatments have been proposed but are not as
of eikx on the r.h.s. (−2∂2U ) which cannot result from
reliable (see ref [20] for more details). Interestingly, we x 0
the application of L2 to any function. Thus α=1 is not
see that the order of the transition and the amplitude of
an option.
theperturbationdependonhownon-lineartermsing(w)
balance the linear growth term −2ε∂2w in (A5). Since For α ≤ 1/3, then α +1 > 1 ≥ 3α. There is thus
x
the former depends on the non-linear relation v(ρ), we no contribution of −2ε∂2w to the first three orders in
x
do not expect equation (A18) to be generic, as opposed the expansion of (A5). In particular, the two first order
to the stability analysis which can be expressed solely in are still given by (A8) and (A9), whereas the third order
terms of De(ρ0) and its derivative. is given by (A10) without the term linear in U0. This
means that the contribution (A14) is not present and
the prefactor of eikx in the r.h.s. of (A10) only contains
2. What is the correct expansion? multiples of |A|2A. The resolvability condition (A17) is
thus of the form A|A|2f(Φ) = 0 which implies |A| = 0.
√
In (A7), we expanded w in power series of ε, thus The only expansion which yields a result is thus for α=
assuming that the amplitude is an analytic function of 1/2.
[1] Shapiro J. A. (1995) The significance of bacterial colony 376:49-53.
patterns. BioEssays 17: 597-607. [5] Woodward D.E. et al. (1995) Spatiotemporal patterns
[2] Harshey R. M. (2003) Bacterial motility on a surface: generated by Salmonella Typhimurium. Biophys. J.
Many ways to a common goal. Ann. Rev. Microbiol. 68:2181-2189.
57:249-273. [6] Murray J.D., Mathematical Biology, Vol. 2, Springer-
[3] Budrene E. O., Berg H. C. (1991) Complex patterns Verlag, New York (2003).
formed by motile cells of Escherichia-coli. Nature 349: [7] Ben-Jacob E., Cohen, I., Levine H. (2000) Cooperative
630-633. self-organization of microorganisms Adv. in Phys. 49,
[4] BudreneE.O.,BergH.C.(1995)Dynamicsofformation 395-554.
ofsymmetricalpatternsbychemotacticbacteria.Nature [8] Kawasaki K. et al (1997) Modeling Spatio-Temporal
9
Figure 7: We simulated equation (A1) for systems of size L = 400 with periodic boundary conditions, using several values
of Φ (From left to right, top to bottom, Φ = 1.06; 1.07; 1.12; 1.2; 1.35; 1.5; 1.7; 1.95; 2.5). The steady-state w(x) was then
decomposed in Fourier series w(x)=a +(cid:80)N/2a cos(2πnx/L)+b sin(2πnx/L), where N is the number of data points and
0 n n n
the cut-off when n→∞ is given by the Nyquist frequency. The blue points correspond to the amplitude of the largest mode:
√
w =max a2 +b2. When the transition is continuous, we compare these points with the results of the amplitude equation
q √ n n n
w =2 ε|A|,where|A|issolutionof (A18)(redlines)andtheagreementisexcellent. ForΦ>1.58orΦ<1.08,thetransition
q
is clearly discontinuous.
Patterns Generated by Bacillus subtilis. J. theor. Biol. and rectification of dilute bacteria. Europhys. Lett.
188:177-185. 86:60002.
[9] Tyson R., Lubkin, S.R., Murray, J.D. (1999) A minimal [19] Angelani L, Di Leonardo R, Ruocco G (2009) Self-
mechanism for bacterial pattern formation. Proc. Roy. StartingMicromotorsinaBacterialBath.Phs.Rev.Lett.
Soc. Lon. B 266:299-304. 102:048104.
[10] Brenner,M.P.,Levitov,L.S.,Budrene,E.O.(1998)Phys- [20] Chaikin P.M., Lubenski T.C, Principles of condensed
ical mechanisms for chemotactic pattern formation by matter physics, Cambridge University Press, Cambridge
bacteria. Biophys. J. 74:1677-1693. (1995).
[11] Espiov S.E., Shapiro J.A. (1998) Kinetic model of Pro- [21] Schnitzer M.J. (1993) Theory of continuum random
teus mirabilis swarm colony. J. Math. Biol. 36:249-268. walks and application to chemotaxis. Phys. Rev. E
[12] Toner J., Tu Y.H., Ramaswamy S. (2005) Hydrodynam- 48:2553-2568.
icsandphasesofflocks.Ann.Phys.(N.Y.)318: 170-244. [22] Berg H.C., E. coli in Motion, Springer (2003).
[13] Ramaswamy S, Simha R.A., Toner J. (2003) Active ne- [23] DickmanR.(1994)Numericalstudyofafieldtheoryfor
matics on a substrate: giant number fluctuations and directed percolation. Phys. Rev. E 50:4404-4409.
long-time tails. Europhys. Lett. 62: 196-202. [24] Dornic I., Chate H., Munoz M.A. (2005) Integration of
[14] Narayan V., Ramaswamy S., Menon N. (2007) Long- Langevinequationswithmultiplicativenoiseandthevi-
lived giant number fluctuations in a swarming granular ability of field theories for absorbing phase transitions.
nematic. Science 317: 105-108. Phys. Rev. Lett. 94: 100601.
[15] TailleurJ,CatesM.E.(2008)Statisticalmechanicsofin- [25] Cates M.E., Vollmer J., Wagner A., Vollmer D. (2002)
teracting run-and-tumble bacteria. Phys. Rev. Lett. 100: Phase separation in binary fluid mixtures with contin-
218103. uously ramped temperature. Phil. Trans. Roy. Soc. A
[16] P. Galajda, J. Keymer, P. Chaikin, R. Austin (2008) J. 361:793-804.
Bacteriol. 189:8704-8707. [26] P. Becherer, A.N. Morozov, W. van Saarloos, (2009)
[17] P. Galajda et al. (2008) Funnel ratchets in biology at Probing a subcritical instability with an amplitude ex-
low Reynolds number: choanotaxis. J. Modern Optics pansion: An exploration of how far one can get. Physica
55:3413-3422. (2008) D 238:1827-1840.
[18] Tailleur J, Cates M-E (2009) Sedimentation, trapping, [27] MittalN.,BudreneE.O.,BrennerM.P.,vanOudernaar-
10
den, A. (2003) Motility of Escherichia coli in clusters Sci. USA 100:13259-13263.
formed by chemotactic aggregation. Proc. Natl. Acad.
