Mechanical interactions in bacterial colonies and the surfing probability of beneficial mutations

5 MB·

by  Fred F. Farrell

Mechanical interactions in bacterial colonies and the surfing probability of
beneficial mutations
Fred F. Farrell1, Matti Gralka2, Oskar Hallatschek2,3, and Bartlomiej Waclaw1,4
1SUPA School of Physics and Astronomy, The University of Edinburgh, Mayfield Road, Edinburgh
EH9 3JZ, UK
2Department of Physics, University of California, Berkeley, CA 94720, USA
7 3Department of Integrative Biology, University of California, Berkeley, California 94720, USA
1
0 4Centre for Synthetic and Systems Biology, The University of Edinburgh, Edinburgh, UK
2

n
a January 13, 2017
J

2
1
Abstract 1 Introduction

]
h
p Bacteria are the most numerous organisms on Earth
-
capable of autonomous reproduction. They have
o
bi Bacterial conglomerates such as biofilms and micro- colonised virtually all ecological niches and are able
. colonies are ubiquitous in nature and play an impor- to survive harsh conditions intolerable for other organ-
s
c tant role in industry and medicine. In contrast to isms such as high salinity, low pH, extreme tempera-
i
s well-mixed, diluted cultures routinely used in micro- tures, or the presence of toxic elements and compounds
y
bial research, bacteria in a microcolony interact me- [1]. Many bacteria are important animal or human
h
p chanically with one another and with the substrate to pathogens, but some bacteria find applications in the
[ which they are attached. Despite their ubiquity, little industry as waste degraders [2] or to produce fuels and

1 isknownabouttheroleofsuchmechanicalinteractions chemicals [3]. In all these roles, biological evolution of
v on growth and biological evolution of microbial popu- microbes is an undesired side effect because it can dis-
7
3 lations. Here we use a computer model of a microbial rupt industrial processes or lead to the emergence of
3 colony of rod-shaped cells to investigate how physical new pathogenic [4] or antibiotic-resistant strains [5].
3
interactionsbetweencellsdeterminetheirmotioninthe Experimental research on bacterial evolution has
0
. colony, this affects biological evolution. We show that been traditionally carried out in well-stirred cultures
1
0 the probability that a faster-growing mutant “surfs” at [6, 7]. However, in their natural environment bacte-
7 the colony’s frontier and creates a macroscopic sector ria often form aggregates such as microcolonies and
1 depends on physical properties of cells (shape, elastic- biofilms. Such aggregates can be found on food [8],
:
v ity, friction). Although all these factors contribute to teeth (plague), on catheters or surgical implants [9],
i
X thesurfingprobabilityinseeminglydifferentways,they inside water distribution pipes [10], or in the lungs of
r allultimatelyexhibittheireffectsbyalteringtherough- people affected by cystic fibrosis [11]. Bacteria in these
a
ness of the expanding frontier of the colony and the aggregates adhere to one another and the surface on
orientation of cells. Our predictions are confirmed by which they live, form layers of reduced permeability to
experiments in which we measure the surfing probabil- detergents and drugs, and sometimes switch to a dif-
ity for colonies of different front roughness. Our results ferent phenotype that is more resistant to treatment
show that physical interactions between bacterial cells [12, 13, 14]; this causes biofilms to be notoriously diffi-
play an important role in biological evolution of new cult to remove.
traits, and suggest that these interaction may be rele- An important aspect of bacteria living in dense con-
vanttoprocessessuchasde novoevolutionofantibiotic glomerates is that they do not only interact via chem-
resistance. ical signaling such as quorum sensing [15] but also
1
(a) (b) (c)
1um
e growth
variable m
length ti
division
F
-F
a small, random “kick”
Figure 1: (a) Illustration of the computer algorithm. Bacteria are modelled as rods of varying length and constant
diameter. When a growing rod exceeds a critical length, it splits into two smaller rods. (b) A small simulated colony. (c)
The same colony with nutrient concentration shown as different shades of gray (white = maximal concentration, black =
minimal); the cells are represented as thin green lines.
through mechanical forces such as when they push experimental scenario in which bacteria grow on the
away or drag other bacteria when sliding past them. surface of agarose gel infused with nutrients. We have
Computer simulations [16, 17, 18, 19] and experiments previously demonstrated [17] that this model predicts
[20, 21, 22, 23, 24] have indicated that such mechanical a non-equilibrium phase transition between a regular
interactions play an important role in determining how (circular) and irregular (branched) shape of a radially
microbial colonies grow and what shape they assume. expanding colony of microbes, and that it can be used
However, the impact of these interactions on biological to study biological evolution in microbial colonies [34].
evolution has not been explored. Here, weusethismodeltoshowthatthesurfingproba-
A particularly interesting scenario relevant to micro- bilityofabeneficialmutationdependsprimarilyonthe
bial evolution in microcolonies and biofilms is that of roughness of the expanding front of the colony, and to
a range expansion [25] in which a population of mi- alesserextendonthethicknessofthefrontandcellular
crobes invades a new territory. If a new genetic variant ordering at the front. We also investigate how mechan-
arises near the invasion front, it either “surfs” on the ical properties of cells such as elasticity, friction, and
front and spreads into the new territory, or (if unlucky) adhesion affect these three quantities. We corroborate
it lags behind the front and forms only a small “bub- some of our results by experiments in which we vary
ble” in the bulk of the population [26]. This stochas- the roughness of the growing front and show that it
tic process, called “gene surfing”, has been extensively influences the surfing probability as expected.
studied [27, 28, 29, 30, 31, 32, 33, 34] but these works
have not addressed the role of mechanical interactions
between cells. Many of the existing models do not 2 Computer model
consider individual cells [27], assume Eden-like growth
[31], or are only appropriate for diluted populations of WeuseacomputermodelsimilartothatfromRefs.[17,
motile cells described by reaction-diffusion equations 23, 34], with some modifications. Here we discuss only
similar to the Fisher-Kolmogorov equation [35]. On thegenericalgorithm; moredetailswillbegiveninsub-
the other hand, agent-based models of biofilm growth, sequent sections where we shall talk about the role of
which have been applied to study biological evolution each of the mechanical factors.
in growing biofilms [36, 37, 38], use very simple rules to We assume that bacteria form a monolayer as if the
mimic cell-cell repulsion which neglect important phys- colony was two-dimensional and bacteria always re-
ical aspects of cell-cell and cell-substrate interactions mained attached to the substrate. This is a good ap-
such as adhesion and friction. proximation to what occurs at the edge of the colony
In this work, we use a computer model of a growing and, as we shall see, is entirely justifiable because the
microbial colony to study how gene surfing is affected edge is the part of the colony most relevant for biolog-
by the mechanical properties of cells and their environ- ical evolution of new traits. We model cells as sphe-
ment. Inourmodel, non-motilebacteriagrowattached rocylinders of variable length and constant diameter
to a two-dimensional permeable surface which delivers d = 2r = 1µm (Fig. 1a). Cells repel each other with
0
nutrients to the colony. This corresponds to a common normal force determined by the Hertzian contact the-
2
ory: F = (4/3)Er1/2h3/2 where h is the overlap dis- Name Value Units
0
tance between the walls of the interacting cells, and E Nutrient diffusion constant D 50 µm2/h
plays the role of the elastic modulus of the cell. The Nutrient concentration c0 1 a.u.
dynamics is overdamped, i.e. the linear/angular veloc- Nutrient uptake rate k 1 – 3 a.u./h
ityisproportionaltothetotalforce/totaltorqueacting Young modulus E 100 kPa
on the cell: Elongation length vl 4 µm/h
Cell diameter 1 µm
d(cid:126)r
i = F(cid:126)/(ζm), (1) Average max. inter-cap distance l 4 µm
dt c
Damping coefficient ζ 500 Pa·h
dφ
i
= τ/(ζJ). (2)
dt
Table 1: Default values of the parameters of the model.
In the above equations (cid:126)r is the position of the centre This gives ≈ 30min doubling time and the average length
i
of bacterium ≈ 3µm. If not indicated otherwise, all results
of mass of cell i, φ is the angle it makes with the x
i
axis, F(cid:126) and τ are the total force and torque acting on presented have been obtained using these parameters.
the cell, m and J are its mass and the momentum of
inertia (perpendicular to the plane of growth), and ζ
(Fig. 2a), and a colony growing in a narrow (width L)
is the damping (friction) coefficient. We initially as-
but infinitely long vertical tube with periodic bound-
sume that friction is isotropic, and explore anisotropic
ary conditions in the direction lateral to the expanding
friction later in Sec. 4.3.
front (Fig. 2d). While the radial expansion case rep-
Bacteria grow by consuming nutrients that diffuse in
resents a typical experimental scenario, only relatively
the substrate. The limiting nutrient concentration dy-
small colonies (106 cells as opposed to > 108 cells in
namics is modelled by the diffusion equation with sinks
a real colony [34]) can be simulated in this way due
corresponding to the bacteria consuming the nutrient:
to the high computational cost. The second method
∂c (cid:18)∂2c ∂2c(cid:19) (cid:88) (growth in a tube) enables us to simulate growth for
= D + −k δ((cid:126)r −(cid:126)r). (3)
∂t ∂x2 ∂y2 i longer periods of time at the expense of confining the
i colony to a narrow strip and removing the curvature
Here(cid:126)r = (x,y), c = c((cid:126)r,t)isthenutrientconcentration of the growing front. This has however little effect on
at position (cid:126)r and time t, D is the diffusion coefficient the surfing probability of faster-growing mutants if the
of the nutrient, and k is the nutrient uptake rate. The width L of the tube is sufficiently large.
initial concentration c((cid:126)r,0) = c . Figure 1b, shows a snapshot of a small colony; the
0
Acellelongatesataconstantratev aslongasthelo- concentrationofthelimitingnutrientisalsoshown. Ta-
l
cal nutrient concentration is larger than a certain frac- ble 1 shows default values of all parameters used in
tion (>1%) of the initial concentration. When a grow- the simulation. Many of these parameters have been
ing cell reaches a pre-determined length, it divides into taken from literature data for the bacterium E. coli
two daughter cells whose lengths are half the length of [34], but some parameters such as the damping co-
the mother cell. The critical inter-cap distance l efficient must be estimated indirectly [17]. We note
cap−cap
at which this occurs is a random variable from a Gaus- that the assumed value of the diffusion constant D is
sian distribution with mean (cid:96) and standard deviation unrealistically small; the actual value for small nutri-
c
±0.15 (cid:96) . Varying (cid:96) allows us to extrapolate between ent molecules such as sugars and aminoacids would be
c c
quasi-sphericalcells(e.g. yeastsS. cerevisaeorthebac- ∼ 106µm2/h, i.e., four orders of magnitude larger. Our
terium S. aureus) and rod-shaped cells (e.g. E. coli or choice of D is a compromise between realism and com-
P. aeruginosa), whereas the randomness of l ac- putational cost; we have also showed in Ref. [17] that
cap−cap
counts for the loss of synchrony in replication that oc- the precise value of the diffusion coefficient is irrelevant
curs after a few generations (the coefficient of variation intheparameterregimeweareinterestedhere. Wealso
ofthetimetodivision∼ 0.1−0.2[39,40,41]). Thetwo note that in reality cessation of growth in the center
daughter cells have the same orientation as the parent of the colony and the emergence of the growing layer
cell, plus a small random perturbation to prevent the may be due to the accumulation of waste chemicals,
cells from growing in a straight line. pH change etc., rather than nutrient exhaustion. Here
We use two geometries in our simulations: a radially we focus on the mechanical aspects of growing colonies
expanding colony that starts from a single bacterium and do not aim at reproducing the exact biochemistry
3
(a) (d) (e)
27 µm 62 µm 140 µm
20
N=63 N=784 N=5450
m)15
ρ (µ10
5
0
0 5 10 15 20
(f) y (mm)
m)150 (b) m/h)8 (c) 20
adius (µ15000 peed (µ246 h (µm)11505
r s
0 0 0
0 5 10152025 0 5 10 15 20 25 0 1 2 3 4
time (h) time (h) y (mm)
Figure2: (a)Snapshotsofaradially-growingsimulatedcolonytakenatdifferenttimes(sizes),fork =2. Growingbacteria
are bright green, quiescent (non-growing) bacteria are dark green. (b) The radius of the colony increases approximately
linearly in time. (c) The expansion speed tends to a constant value for long times. (d) Example configuration of cells
from a simulation in a tube of width L = 80µm. The colony expands vertically. h is the thickness of the growing layer
(Eq. (4)), ρ is the roughness of the front (Eq. (5)). (e,f) Thickness and roughness as functions of the position y of the
front, for L=1280µm and k =2.5, and for 10 indepedent simulation runs (different colours).
of microbial cells, as long as the simulation leads to the Measuring surfing probability. For each pair of
formation of a well-defined growth layer (as observed mutantandwildtype, amixedstartingpopulationwas
experimentally). prepared that contained a low initial frequency P of
i
mutantshavingaselectiveadvantages. Colonygrowth
was initiated by placing 2µl of the mixtures onto plates
and incubated until the desired final population size
3 Experiments
was reached. The initial droplet radius was measured
to compute the number of cells at the droplet perime-
Experiments were performed as described in our previ- ter. The resulting colonies were imaged with a Zeiss
ous work [34]. Here we provide a brief description of AxioZoom v16. The number of sectors was determined
these methods. by eye. The surfing probability was calculated using
Strains and growth conditions. For the mixture Eq. (10).
experiments measuring surfing probability, we used Timelapse movies. For single cell-scale timelapse
pairs of microbial strains that differed in fluorescence movies, we used a Zeiss LSM700 confocal microscope
color and a selectable marker. The selective differ- with a stage-top incubator to image the first few layers
ence between the strains was adjusted as in [34] us- of most advanced cells in growing S. cerevisiae and E.
ing low doses of antibiotics. The background strains coli colonies between a coverslip and an agar pad for
and antibiotics used were E. coli DH5α with tetracy- about four hours, taking an image every minute.
cline, E. coli MG1655 with chloramphenicol, and S. Measuring roughness. Imagesofatleast10equal-
cerevisiae W303 with cycloheximide. Selective differ- sized colonies per condition were segmented and the
ences were measured using the colliding colony assay boundary detected. The squared radial distance δr2
[32]. E. coli strains were grown on LB agar (2%) between boundary curve and the best-fit circle to the
medium (10g/L tryptone, 5g/L yeast extract, 10g/L colony was measured as a function of the angle and
NaCl)ateither37◦Cor21◦C.S.cerevisiae experiments averaged over all possible windows of length l. The
were performed on either YPD (20g/L peptone, 10g/L resultingmeanδr2 wasaveragedoverdifferentcolonies.
yeast extract, 20g/L glucose) or CSM (0.79g/L CSM Images of moving fronts at the single-cell level from
(Sunrise media Inc.), 20 g/L glucose) at 30◦C. 20g/L the timelapse movies were first segmented using a local
agar was added to media before autoclaving. Antibi- adaptative threshold algorithm to identify cells. The
otics were added after autoclaving and cooling of the front was found by the outlines of cells directly at the
media to below 60◦C. front. For all possible windows of length l, a line was
4
fitted to the front line and the mean squared distance
fromthebest-fitlinewasmeasured,asinRef.[27]. The
resulting mean squared distance was averaged over all
(a)
windows of length l and all frames.
4 Simulation results
(b)
4.1 Growth and statistical properties of the
simulated colony
Wenowdiscussthepropertiesofoursimulatedcolonies.
When the colony is small, all bacteria grow and repli-
(c)
cate. As the colony expands, the nutrient becomes de-
pleted in the centre of the colony because diffusion of Figure 3: The frontier of the colony for three different
the nutrient cannot compensate its uptake by growing nutrient uptake rates k = 1.8 (a), k = 2.2 (b) and k =
cells. This causes cessation of growth in the centre. 2.6 (c). The thickness of the growing layer (bright green)
When this happens, growth becomes restricted to a decreases only moderately (1.64×) from h = 13.5±0.1µm
for k = 1.8 to h = 8.2±0.1µm for k = 2.6, but this has
narrow layer at the edge of the colony, see Fig. 2a,
a large impact on the front roughness which changes from
and Supplementary Video 1. The radius of the colony
ρ = 2.1±0.2µm to ρ = 9.3±0.4µm, correspondingly. For
increases approximately linearly in time (Fig. 2b,c).
k = 2.6 the growing layer begins to loose continuity and
The presence of a “growing layer” of cells and the linear
splits into separate branches.
growth of the colony’s radius agree with what has been
observed experimentally [42, 34].
Statisticalpropertiesofthegrowinglayercanbecon-
veniently studied using the “tube-like” geometry. Fig-
After a short transient the expansion velocity, the
ure 2d shows a typical configuration of cells at the
nutrient profile, and other properties of the growing
colony’sfrontier(seealsoSupplementaryVideo2). The
layer stabilize and vary little with time (Fig. 2e,f). It
growing layer can be characterized by its thickness h
is therefore convenient to choose a new reference frame
and roughness ρ which we calculate as follows. We first
co-moving with the leading edge of the colony. Since
rasterize the growing front of the colony using pixels
cells that lag behind the front do not replicate, we do
of size 1 × 1µm, and find the two edges of the front:
nothavetosimulatethesecellsexplicitly. Thisdramat-
the upper one (the colony edge) {y+} and the lower
i ically speeds up simulations and enables us to study
one (the boundary between the growing and quiescent
stripes of the colony of width L > 1mm and length
cells) {y−}. We then calculate the average thickness as
i > 10mm.
h = 1 (cid:88)L min (cid:113)(i−j)2+(y+−y−)2. (4) We have shown previously [17] that the thickness of
L j=1,...,L i j the growing layer of cells is controlled by the nutrient
i=1
concentration c , nutrient uptake rate k, growth rate b,
0
This method takes into account that the growing layer
and elasticity E of cells. This in turn affects the rough-
can be curved and does not have to run parallel to the
ness of the leading edge of the colony, see Fig. 3, where
x axis1. Similarly, we calculate the average roughness
we vary the uptake rate k while keeping the remain-
as
(cid:118) ing parameters constant. Figure 4 shows that front
(cid:117) L
(cid:117)1 (cid:88) thickness decreases and its roughness increases with in-
ρ = (cid:116) (y+−Y+)2 , (5)
L i creasing k; eventually, when a critical value k ≈ 2.5 is
c
i=1
crossed,thegrowingfrontsplitsintoseparatebranches.
where Y+ = (1/L)(cid:80) y+. Note that all quantities This transition has been investigated in details in Ref.
i i
(L,Y+,y+,y−) are in pixels and not µm. [17]. Although this scenario can be realized experimen-
i i
tally [43, 44], here we focus on the “smooth” regime
1Alternatively,hcanbedefinedastheareaofthecolonythat
in which colonies do not branch out and the frontier
containsreplicatingcellsdividedbytheinterfacelengthL. Both
methods produce similar results. remains continuous.
5
(a) (b) (c)
)ess ρµm (12235050 LLLL====1361624200080 ess h (µm) 11220505 )ess ρµm (12235050
n 10 n n 10
gh 5 ck 5 gh 5
rou 0 thi 0 rou 0
1.0 1.5 2.0 2.5 3.0 1.0 1.5 2.0 2.5 3.0 5 10 15 20
uptakerate k uptakerate k thickness h (µm)
Figure 4: Thickness and roughness of the growing layer for different front lengths (tube widths) L=160 (red), L=320
(green), L = 640 (blue), and L = 1280 µm (purple). (a) Thickness h decreases as the nutrient uptake rate k increases.
h does not depend on the length L of the front. (b) Roughness ρ increases with both k and L. (c) Roughness versus
thickness; different points correspond to different k from the left and middle figure.
4.2 Surfing probability of a beneficial mu-
tation
Whenamutationarisesatthecolony’sfrontier, itsfate
(a) (b) (c)
can be twofold [27, 34]. If cells carrying the new mu-
tation remain in the active layer, the mutation “surfs” Figure 5: The fate of mutants. Left and middle panels
on the moving edge of the colony and the progeny of showdifferentfatesofasectoroffitter(s=0.1)mutantcells
(red)inacolonyof“wild-type” cells(green). Thesectorcan
the mutant cell eventually forms a macroscopic “sec-
either expand (left panel) or collapse and become trapped
tor” (Fig. 5). On the other hand, if cells carrying the
in the bulk when random fluctuations cause mutant cells to
mutation leave the active layer, the mutation becomes
lag behind the front (middle panel). Right panel shows a
trapped as a bubble in the bulk of the colony [26]. Due
sector with larger (s = 0.5) growth advantage; significantly
to the random nature of replication and mixing at the fastergrowthofmutantcellsleadstoa“bump” atthefront.
front, surfing is a stochastic process; a mutation re- In all cases k =1.8,L=160µm.
mains in the active layer in the limit t → ∞ with some
probability P which we shall call here the surfing
surf
probability. with some small probability per division. The simula-
Surfing is a softer version of fixation - a notion from tion finishes when either fixation (all cells in the grow-
population genetics in which a mutant takes over the ing layers becoming mutants) or extinction (no mutant
population. The soft-sweep surfing probability has cells in the growing layer) is achieved. Before inserting
therefore a hard-selection-sweep counterpart, the fix- the mutant cell, the colony is simulated until the prop-
ation probability, which is the probability that the new erties of the growing layer stabilize and both thickness
mutation spreads in the population so that eventually and roughness reach steady-state values. The simu-
all cells have it. Both surfing and fixation probabili- lation is then repeated many times and the probabil-
ties depend on the balance between selection (how well ity of surfing is estimated from the proportion of runs
the mutant grows compared to the parent strain) and leading to fixation of the mutant in the growing layer.
genetic drift (fluctuations in the number of organisms Snapshots showing different fates (extinction, surfing)
due to randomness in reproduction events) [45]. In the of mutant sectors are shown in Fig 5.
previous work [34] we showed that P increased ap- Surfing probability depends on the position of
surf
proximately linearly with selective advantage s – the thecellinthegrowinglayer. InRef. [34]weshowed
difference between the growth rate of the mutant and that the surfing probability strongly depends on how
the parent strain. Here, we study how the properties deeply in the growing layer a mutant was born. Here
of the active layer affect P for a fixed s. we would like to emphasize this result as it will become
surf
Wefirstrunsimulationsintheplanar-frontgeometry important later. Let ∆ be the distance from the edge
inwhicharandomcellpickedupfromthegrowinglayer of the colony to the place the mutant first occurred.
of cells with probability proportional to its growth rate Figure 6 shows the probability density P(∆|surf) that
is replaced by a mutant cell with selective advantage a cell was born a distance ∆ behind the colony front,
s > 0. This can be thought of as mutations occurring given that it went on to surf on the edge of the expand-
6
(a) (a) 0.10
3 k=1.6
0.08
rf) k=2.0
su 2 k=2.4 0.06
Δ| urf
P( 1 Ps
0.04 L=160
L=320
0 0.02 L=640
0 1 2 3 4
L=1280
Δ (µm) 0.00
8 10 12 14 16
(b) 4
h (µm)
s=0.00 s=0.1
) 3
rf s=0.02 s=0.2 (b) 0.10
u
Δ|s 2 s=0.05 s=0.5 LL==136200
( 0.08
P L=640
1
0.06 L=1280
0 Psurf
0 1 2 3 4 0.04
Δ (µm)
0.02
Figure 6: (a) P(∆|surf) for L = 160µm, selective advan-
0.00
tage s = 0.02, and different k = 1.6,2.0,2.4. (b) P(∆|surf)
2 4 6 8 10 12
for L = 160µm, k = 2.0, and different selective advantages
s = 0,0.02,0.05,0.1,0.2,0.5. Only mutants from the first ρ (µm)
layer of cells have a significant chance of surfing.
Figure 7: (a)P fordifferentthicknesshofthegrowing
surf
layer, for s=0.02 and L=160,320,640,1280 µm (different
colours). (b) the same data as a function of front roughness
ing colony. It is evident that only cells born extremely
ρ. Between103 and104 simulationswereperformedforeach
close to the frontier have a chance to surf. Cells born
data point to estimate P .
surf
deeper must get past the cells in front of them. This
is unlikely to happen, even if the cell has a significant
growth advantage, as the cell’s growth will also tend to how the surfing probability P varies as a function of
surf
push forward the cells in front of it. This also justifies the thickness and the roughness of the front. P in-
surf
why we focus on 2d colonies; even though real colonies creases with increasing thickness h and decreases with
are three-dimensional, all interesting dynamics occurs increasing roughness ρ. We know from Fig. 4 that
at the edge of the colony, made of a single layer of cells. thickness and roughness are inversely correlated so this
Given that surfing is restricted to the first layer of reciprocal behaviour is not surprising. An interesting
cells, and the distribution P(∆|surf) is approximately question is which of the two quantities, roughness or
the same for all explored parameter sets (different k thickness, directly affects the probability of surfing?
and s), it may seem to be a waste of computer time From a statistics point of view, thickness h seems to
to study the fate of mutants that occurred deeply in be a better predictor of P because data points for
surf
the growing layer. To save the time, and to remove the the same h but for different L correlate better. How-
effect the front thickness has on P (thicker layer = ever, it could be that it is actually front roughness that
surf
lower overall probability), we changed the way of intro- directly (in the causal sense) affects the surfing proba-
ducingmutants. Insteadofinsertingmutantsanywhere bility and that P and h are anti-correlated because
surf
in the growing layer, we henceforth inserted them only of the relationship between h and ρ.
at the frontier. We performed two computer experiments to address
Roughness of the front is more predictive of the above question. First, we simulated a colony that
P than its thickness. Using the new method of had a very low and constant roughness ρ ≈ 1, inde-
surf
introducing mutants (only the first layer of cells), we pendently of front’s thickness. This was achieved by
run simulations for s = 0.02 and for different widths L introducing an external force F = −gy acting on the
y
andnutrientuptakerateskasinFig. 4. Figure7shows centre of mass of each cell, where g > 0 was a “flatten-
7
ingfactor” whosemagnitudedeterminedthestrengthof
suppressionofdeviationsfromaflatfront. P plotted
surf
in Figure 8, left, as a function of h for two cases: “nor-
mal”, rough front, and “flattened” front, demonstrates
that the surfing probability does not depend on h in
the case of flat front.
Second, we varied roughness while keeping thickness
constant. This was done by measuring front roughness
in each simulation step, and switching on the external
“flattening” force Fy = −gy if the roughness was larger (a)
than a desired value ρ . Figure 8, right, shows that
max 0.08
althoughthicknessremainsthesameforalldatapoints,
P decreases with increasing roughness. 0.06
surf urf
We can conclude from this that it is the increase Ps
0.04
in the roughness, and not decreasing thickness, that
lowers the surfing probability for thinner fronts (larger 0.02 normal
flattened
nutrient intake rate k). However, the data points in
0.00
Fig. 7, right, from different simulations do not collapse
8 10 12 14 16
onto a single curve as it would be expected if average,
h (µm)
large-scale front roughness was the only factor.
Local roughness predicts P . According to the (b)
surf
normal
theoryofRef. [29],thedynamicsofamutantsectorcan
0.08
flattened
be described by a random process similar to Brownian
motion in which the sector boundaries drift away from urfurf 0.06
eachotherwithconstantvelocity. Thevelocitydepends PPss
0.04
on the growth advantage s whereas the amplitude of
random fluctuations in the positions of boundary walls 0.02
is set by the microscopic dynamics at the front. We
0.00
reasoned that these fluctuations must depend on the
2 4 6 8 10
roughness ρ of the frontier, and that a mutant sector
ρ (µm)
should be affected by front roughness when the sec-
tor is small compared to the magnitude of fluctuations. Figure 8: (a) P as the function of front thickness h for
surf
This means that local roughness ρ(l), determined over the normal (black) and flattened front (red, g = 500), for
L=320µm. We vary the nutrient uptake rate k =1.6...2.8
the length l of the front, should be more important
to simulate fronts of different thickness. The flat front has
than the global roughness ρ(L). We calculated the lo-
roughness ρ between 0.84 and 1.0 for all k. (b) P for the
cal roughness as surf
normal (black) and flattened front (blue) as the function of
(cid:118) roughness ρ. The flattened front has approximaly the same
n (cid:117) i+l
1 (cid:88)(cid:117)1 (cid:88) thickness for all data points (h between 10.0 and 10.3µm).
ρ(l) = (cid:116) (y+−Y+)2. (6)
n l j Thepointscorrespondtomaximumroughnesssettoρ =
max
i=1 j=i
2,3.5,5, and 7, for k = 2.6; the actual (measured) ρ differs
HereY+ istheaverageheightoftheinterfaceand{y+} very little from these values.
i
are the vertical coordinates (interface height) of the
points at the leading edge, obtained as in Section 4.1.
Figure 9 shows that P for different L now collapse
surf
ontoasinglecurve,foralllengthsl ≈ 10...100µmover
which roughness has been calculated.
Orientation of cells affects P . So far we have
surf
focused only on the macroscopic properties of the lead-
ing edge of the colony, completely neglecting its gran-
ular nature due to the presence of individual cells. Re-
8
call that in our model each cell is rod-shaped, and the 4.3 Surfing probability and the mechanical
direction in which it grows is determined by the ori- properties of bacteria
entation of the rod. Figure 10a shows that cells at
Our results from the previous section demonstrate that
the leading edge assume orientations slightly more par-
surfing is affected by (i) the roughness of the growing
allel to the direction of growth (vertical) in the flat-
layer, (ii) the orientation of cells, (iii) the thickness of
tened front than in the normal simulation. A natural
the growing layer if mutations occur inside the growing
question is how does cellular alignment affects P ,
surf
layer and not only at its edge. To show this, we varied
independently of the roughness? To answer this ques-
thickness, roughness, and orientation of cells by using
tion, we simulated a modified model, in which external
ad hocexternalforcesflatteningoutthefrontorforcing
torque τ = −τ sin[(φ − φ ) mod π] was ap-
max preferred
the cells to order in a particular way. In this section
plied to the cells, forcing them to align preferentially
we will investigate what parameters of the model affect
in the direction φ . We investigated two forced
preferred
surfing in the absence of such artificial force fields.
alignments: φ = 0 corresponding to cells par-
preferred
Thickness of the growing layer. If cells are pro-
allel to the x axis and hence to the growing edge of
hibitedtoformmultiplelayers,asinour2dsimulations,
the colony, and φ = π/2 which corresponds to
preferred
thickness h can be determined from the parameters of
the vertical orientation of cells (perpendicular to the
the model by a simple dimensional analysis. Assum-
growing edge).
ing that h is proportional to the characteristic scale
Figure 10b compares these two different modes with over which the nutrient concentration and cell density
previous simulations with no external torque, for ap- reaches bulk values [17], we can approximate h by
proximately the same thickness and roughness of the
(cid:115)
growing layer. It is evident that the orientation of cells E
h ≈ (1/β−1)3/4, (7)
strongly affects the surfing probability: horizontally- (ζ/a)φ
forced cells have ∼ 3x smaller P compared to the
surf
where E is the elastic modulus of the bacterium (Pa),
normal case, which in turn has P ∼ 5x smaller than
surf
a is the average area per cell (µm2), ζ is the friction
vertically-forced cells.
coefficient (Pa·h), φ is the replication rate (h−1), and
Shorter cells have higher P than long cells.
surf β < 1 is a dimensionless ratio of the nutrient consump-
To check how the aspect ratio of cells affect P , we
surf tionratetobiomassproductionrate(i.e. newbacteria):
simulated cells whose maximal length was only 2µm
β = (kρ )/(φc ). Equation (7) shows that thickness h
0 0
and the minimal separation before the spherical caps
increases with increasing cell stiffness (larger E) and
was zero, i.e., the cells became circles immediately af-
replication rate φ, and decreases with increasing nutri-
ter division. As before we selected a set of k’s such
entuptakek andincreasingfrictionζ. Theaspectratio
that the thickness and roughness were approximately
of the cells does not affect h in our model. Equation
the same for all simulations. In order to make a fair
(7) suggests that the thickness of the growing layer can
comparison between “short rods” and “long rods” from
be conveniently controlled in an experiment by varying
previous simulations, thickness and roughness were ex-
temperature or growth medium (which both affect the
pressed in cell lengths rather than in µm. This was
growth rate), or by varying the nutrient concentration
done by dividing both h and ρ by the average length of
c . We shall use the first two methods when discussing
0
a cell measured for cells from the growing layer. Figure
the experimental verification of our theory.
10c show that short rods have a much higher surfing
Orientationofcells. Ausefulmeasureoftheglobal
probability than long rods.
alignment of cells in the colony is the order param-
In all previous simulations, even for short rods, cells eter S = (cid:10)cos2(φ−Φ)(cid:11). Here φ is the angle a cell
remembered their orientation from before division and makes with the x-axis and Φ is the angular coordinate
growth always initially occurred in that direction. To of the vector normal to the front; this is to remove a
see whether this has any impact on P , we consid- trivial contribution to S due to the curvature of the
surf
ered a scenario in which the new direction of growth front caused by roughness. According to this defini-
is selected randomly and does not correlate with the tion, S = 1 if all cells are perfectly vertically aligned
direction prior to division. Figure 10c shows that P (in the direction of growth), S = 0 if they are hori-
surf
almost does not change regardless whether a short cell zontal (parallel to the front), and S = 1/2 if their ori-
randomly changes its orientation after division or not. entations are random. It turns out that changing the
9
(a) (b) (c)
0.10 0.10 0.10
L=160
0.08 0.08 0.08 L=320
urf 0.06 urf 0.06 urf 0.06 L=640
Ps Ps Ps L=1280
0.04 0.04 0.04
0.02 0.02 0.02
ℓ=10 ℓ=35 ℓ=98
0.00 0.00 0.00
0.8 0.9 1.0 1.1 1.2 1.3 1.0 1.5 2.0 2.5 3.0 2 3 4 5 6
ρ(ℓ) (µm) ρ(ℓ) (µm) ρ(ℓ) (µm)
Figure 9: P as the function of local roughness ρ(l) of the growing layer, for different sizes L=160,320,640,1280 µm
surf
(as in Fig. 7) and s=0.02. Left: l=10, middle: l=35, right: l=98 µm. For each l, data points for different L collapse
onto a single curve.
(a) flatten, g=500
normal
vertical
horizontal
(b) (c)
P
surf P
0.35 surf
h µm) 0.14 h cells)
0.30 1 2( (6
8 0.12 4
0.25 4 2
0 0.10 0
0.20 ρ µm) ρ cells)
4( 0.08 4(
0.15 2 0.06 2
0 0
0.10
0.04
0.05 0.02
no forced horizontal vertical long rods short rods short rods
alignment alignment alignment randomized
Figure 10: (a) Orientation of cells (colours as in the circle in the upper-right corner) in the growing layer for different
models. (b, c) Comparison of fixation probabilities for different cellular alignments at the front, for approximately the
same thickness and roughness, both of which were controlled by varying k. To achieve this, different k needed to be used
inpanels(b,c)andhencethetwopanelscannotbedirectlycompared. InallcasesL=320µm,s=0.02. Forhorizontally-
and vertically-forced cells, τ =10000. Short cells have a maximum length of 2µm; upon division, they become circles
max
of diameter 1µm.
10