1. Computational and Systems Biology
Download icon

Spatiotemporal establishment of dense bacterial colonies growing on hard agar

  1. Mya R Warren
  2. Hui Sun
  3. Yue Yan
  4. Jonas Cremer
  5. Bo Li  Is a corresponding author
  6. Terence Hwa  Is a corresponding author
  1. University of California, San Diego, United States
  2. California State University, Long Beach, United States
  3. Fudan University, China
Research Article
  • Cited 1
  • Views 745
  • Annotations
Cite this article as: eLife 2019;8:e41093 doi: 10.7554/eLife.41093

Abstract

The physical interactions of growing bacterial cells with each other and with their surroundings significantly affect the structure and dynamics of biofilms. Here a 3D agent-based model is formulated to describe the establishment of simple bacterial colonies expanding by the physical force of their growth. With a single set of parameters, the model captures key dynamical features of colony growth by non-motile, non EPS-producing E. coli cells on hard agar. The model, supported by experiment on colony growth in different types and concentrations of nutrients, suggests that radial colony expansion is not limited by nutrients as commonly believed, but by mechanical forces. Nutrient penetration instead governs vertical colony growth, through thin layers of vertically oriented cells lifting up their ancestors from the bottom. Overall, the model provides a versatile platform to investigate the influences of metabolic and environmental factors on the growth and morphology of bacterial colonies.

https://doi.org/10.7554/eLife.41093.001

Introduction

Bacteria often form dense biofilms with complex spatiotemporal structures (Costerton et al., 1995; Nadell et al., 2016; O'Toole et al., 2000; Stoodley et al., 2002). Mechanical and biochemical interactions, together with cell growth, motility, and signaling, are some of the common elements underlying the rich variety of patterns and behaviors observed. Biofilms often play important roles in diverse settings ranging from environment to human health (Costerton et al., 1999; Jayaraman and Wood, 2008; Potera, 1999). But they are notoriously difficult to study experimentally because of their opaqueness, high heterogeneity and complex organization, involving multiple spatial and temporal scales (Roberts et al., 2015; Stewart and Franklin, 2008). In addition, biofilm-bound bacteria alter their micro-environment by secreting various polysaccharides, forming heterogeneous matrices of filaments that bind cells together within biofilms (Branda et al., 2005; Flemming and Wingender, 2010).

Over the years, various computational models have been constructed to capture different aspects of biofilm development (Alpkvist et al., 2006; Espeso et al., 2015; Ginovart et al., 2002; Klapper and Dockery, 2002; Kreft et al., 2001; Kreft et al., 1998; Picioreanu et al., 2004; Seminara et al., 2012; Tierra et al., 2015). However, most of these models are ‘descriptive’ in nature – the complexity of the biofilms makes it difficult to make quantitative comparison between experimental data and model predictions. In recent years, an increasing body of literature has been devoted to simpler, stripped down versions of the biofilm which can be more readily compared to experimental studies. The simplest among these is the growth of a simple bacterial colony on hard agar surface, with cells pushing against each other by the force of their own physical growth, without motility and without extracellular polysaccharides (Boyer et al., 2011; Cole et al., 2015; Farrell et al., 2013; Ghosh et al., 2015; Grant et al., 2014; Jayathilake et al., 2017; Rudge et al., 2013; Rudge et al., 2012; Volfson et al., 2008) In addition to serving as simpler models of biofilms, the growth of such colonies has been increasingly used in recent years as a model of microbial range expansion in studies of population genetics and ecology (Hallatschek et al., 2007; Hallatschek and Nelson, 2010; Korolev et al., 2012). Although the growth of such simple colonies has been investigated experimentally many decades ago (Cooper et al., 1968; Lewis and Wimpenny, 1981; Mitchell and Wimpenny, 1997; Palumbo et al., 1971; Pirt, 1967; Reyrolle and Letellier, 1979; Wimpenny, 1979), surprisingly, there has not yet been a common quantitative understanding of the basic elements controlling their growth, for example what factors determine the radial and vertical expansion speeds.

In this study, we develop a conceptually simple, yet physically realistic three-dimensional computational model, incorporating the elements of nutrient diffusion, cell-cell and cell-agar mechanical interactions, and introducing a unique cell-level model of surface tension. Our model is efficiently implemented with a parallel algorithm, enabling the simulation of a colony comprising a few million cells within 24 hr. The model is able to capture many observed features of the growing colonies, including the conic shape, the linear growth of the colony radius and height, and their dependence on the cell growth rate. Extensive analysis of the results reveals key driving forces underlying these observations, especially on the role of surface tension and the dynamic form of cell-agar friction, allowing us to make distinct predictions on how various biochemical and mechanical effects alter physiological features of the colony and generate macroscopic spatiotemporal patterns of the growing colony. To guide the construction of our model and validate our simulations, we conducted a series of experiments on the growth of colonies on agar using non-motile E. coli. A set of minimum media with various carbon sources was used to vary the cell growth rate.

Results

Experimental results

Experiments were performed using E. coli K12 strain EQ59, which is non-motile and harbors constitutive GFP expression; see 'Experimental Methods'. Each colony was inoculated as a single cell from batch culture growing in mid-log phase on 1.5% (w/v) agar with glucose minimal media, and incubated, covered, at 37°C for up to 1.5 days. The colony height profile was periodically monitored using a confocal microscope (see 'Experimental Methods'), and the result was highly repeatable; see Figure 1—figure supplement 1. Starting with a single cell, the colony remained a single layer through the first 13 hours (Figure 1AB), buckling into a second layer at around t=14 h at a radius of ~75 μm (Figure 1C–E and F). It then developed into a 3D colony over time, maintaining an approximate conic shape through the ensuing 10-15 hours after buckling (Figure 1G). During this period which we refer to as the ‘establishment phase’, the colony radius increased linearly in time with a constant radial speed VR45.2 μm/h and the colony height increased also linearly at a vertical speed VH12.4 μm/h (Figure 1H), reaching a radius of 500 μm and a height of 150 μm by t=24 h. As the colony grew further, the gain in height slowed down while radial expansion continued at the same speed (Figure 1H and Figure 1—figure supplement 1), leading to a significant flattening of colony morphology. In this study, we focus on the relatively simple establishment phase defined by 14t24 h, where both the radial and vertical growth are linear.

Figure 1 with 1 supplement see all
Experimental observations of the growth and morphology of a bacterial colony.

(A–E) Confocal images of an E. coli colony harboring GFP expression growing on 1.5% agar (glucose minimal medium) taken at various time after seeding (t=0). The center of the colony is indicated by the red dot. Single- and multi-layer regions are distinguished by red circles based on fluorescence intensity; see 'Experimental Methods'. (F) The radius of the first (red) and second layer (orange) of the colony, as well as their difference (green), versus time. (G) The cross-sectional profile of the growing bacterial colony at indicated time after single-cell inoculation. (H) After the buckling at around t=13 h, the colony radius (red symbols) increased at a constant speed VR=45.2 μm/h (red line), while the colony height (blue symbols) increased linearly with speed VH=12.4 μm/h (blue line). The latter slowed down some time after t=24 h. (I) The dependence of the radial speed VR (red symbols) and the vertical speed VH (blue symbols) on cell growth rate (x-axis), for colonies grown in minimal medium with 8 different carbon sources (Supplementary file 1-Table S1): glucose (O); arabinose (□); mannitol (); maltose (◇); fructose (); melibiose (); sorbitol (); mannose (☆). The lines are best linear fit of the data.

https://doi.org/10.7554/eLife.41093.002

We further probed the growth of colony using saturating amounts of different carbon sources, each supporting a different batch culture growth rate, spanning the range 0.5 h-1 to 1 h-1; see Supplementary file 1-Table S1. The radial and vertical expansion speeds obtained in the linear growth regime are plotted in Figure 1I against the batch culture growth rate in the respective medium. Our findings of vertical linear growth disagree with earlier finding by Pirt (Pirt, 1967) which was first questioned by Wimpenny (Lewis and Wimpenny, 1981Wimpenny, 1979). However, the latter reported much larger radial expansion speeds than ours, suggesting that their study might be in a very different regime dominated by swarming motility (Wu et al., 2011).

Simulation results and analysis

To describe the morphology and dynamics of these growing colonies in the linear regime (the establishment phase), we focus on several main elements in the process: the supply of nutrient and interaction driven by the physical growth of cells. We construct a minimal, multiscale, three-dimensional model consisting of the diffusion of nutrient through the agar and the colony; the growth, division, and movement of individual cells; and the cell-cell, cell-agar, cell-surface mechanical interactions that generate forces driving cell movement; see Figure 2. A salient summary of the model is provided in Materials and methods. As will be descripted, a unique aspect of this model is the implementation of the surface tension, which enables us to capture bulk as well as single layer effects. We use the data from our experiments and literature to estimate the range of key parameters in the model, and implement our model using various numerical techniques. Details of the model and numerical methods are given in Appendix 1. Through the bulk of the study described below, a standard set of parameters were used (Supplementary file 1-Tables S2-S4); effects due to variation of parameter values are discussed towards the end.

Schematics of cell-cell, cell-agar, cell-fluid, and surface tension forces investigated in this study.

Green area indicates the colony (with cells in yellow). Blue area indicates the agar.

https://doi.org/10.7554/eLife.41093.006

Radial and vertical growth of the colony

We start by examining how fast the colony expands radially and vertically. We run a simulation with the batch culture growth rate λS=1.0 h-1, which corresponds approximately to the growth of E. coli in glucose minimal medium (Supplementary file 1-Table S1). We use a substrate concentration Cs=0.5 mM here and will vary this parameter later. From Figure 3A, we see that the number of cells in a colony increases exponentially for approximately 10 hours before it slows down. From Figure 3B, we see that the cross-sectional profiles of the colony preserve their shapes and are evenly separated at equal time intervals for t12 h, suggesting a constant expansion of the colony in the radial and vertical directions by t=12 h, similar to the experimental profiles in Figure 1G. (The spatial cell density inside the colony is constant, ~0.68 ρcell, throughout the interior of the colony; see Figure 3—figure supplement 1.) Detail of the profile at the colony periphery appears to be different. This is due to an approximate height assignment based simply on thresholding the fluorescence intensity to obtain the global height profile. This thresholding procedure does not capture height at the periphery where it is one to a few layers in thickness. Figure 3C provides a quantitative picture of the colony radius (R, defined as the average radius of the bottom layer of the colony) and colony height (H, defined as the height at the center of the colony). At early time, t6 h, the colony expands radially, while the height remains close to zero, indicating that the colony is comprised of a thin layer (see discussion in 'Radial expansion – quantitative analysis'). At around t=7 h (indicated by the blue arrow in Figure 3C), the height starts to increase, indicating the occurrence of ‘buckling’. Details of this transition is shown in Figure 3D and E–I; they correspond well to the experimental patterns observed in Figure 1F and A–E. In particular, the model generates a constant width for the single-layer annulus region at the periphery, recapitulating report of a constant monolayer region by earlier mechanical study (Su et al., 2012). Moreover, the model captures the dynamical details around the buckling transition (compare Figures 3D and 1F), which exhibits an initial fast increase of the annulus width resulting from the initial non-compact nature of the cells forming the second layer; see Figure 3J–M. After that point, both the colony radius and height increase linearly with time, with radial expansion speed VR18 μm/h and vertical ascending speed VH6 μm/h; see Figure 3C. Thus, our model captures the linear increase of both the colony radius and height observed experimentally (Figure 1H). To understand the origin of these behaviors, we will analyze below the model output, first pictorially and then quantitatively. The lower numerical values of the speeds obtained from simulations are due to parameter settings chosen to limit computational time; this will be discussed in 'Parameter dependence'.

Figure 3 with 1 supplement see all
The simulated growth and morphology of a bacterial colony.

(A) A semi-log plot of number of cells vs. time showing the exponential growth of the population starting from a single cell at t=0. Red line shows exponential growth rate of 0.96 h-1. (B) Cross sections of the growing colony at various times after seeding of a single cell. (C) Plots of the radius (red) and height (blue) vs. time up to t=20 h, showing that, after an initial transient period of ~10 h, the growing colony increases linearly at the radial speed VR18 μm/h (red line) and vertical speed VH6.0 μm/h. The blue arrow at t=6 h indicates the time when the colony height starts to increase. (D) The radius of the first (red) and second (orange) layer of the colony as well as their difference (green) vs. time. (E–M) Top view of the colony at various time. Cells in the bottom layer (blue) and upper layers (yellow) are fitted into red circles. The time evolution of buckling phenomenon is captured in detail in (J–M).

https://doi.org/10.7554/eLife.41093.007

Vertical rise – a pictorial view

We first focus on factors driving the linear vertical rise of the colony. We start with a pictorial view of the cell configuration and motion inside the colony. Figure 4A shows a snapshot of cell configuration in a vertical slice through the center of the colony, taken at time t=20 h which is well in the steady linear growth regime. The colors distinguish the gross orientations of the cells. The model shows that cells near the top surface are oriented parallel to the colony surface (shown in cyan), while cells away from the top surface are mostly oriented vertically (shown in yellow). A detailed view of the top surface of the colony generated from the simulation is shown in Figure 4—figure supplement 1A. This prediction is validated by confocal scan of the colony in experiment as shown in Figure 4—figure supplement 1B.

Figure 4 with 2 supplements see all
The cross-sectional anatomy of a simulated colony.

(A) Snapshot of cross-sectional view of the colony at t=20 h. Cyan represents horizontally oriented cells (45o with z-axis); Yellow represents vertically oriented cells (45o with z-axis). (B) Fraction of vertically oriented cells averaged over z vs radius. (C) A side view of the azimuthally averaged director field, indicating the orientation of the rod-like cells. (D) A side view of the azimuthally averaged velocity field. (E) Vertical component of velocity, Vz, at various values of z along the center of the colony. Increase in vertical speed is seen only for the bottom 10 µm (F) A cross-sectional view of the colony, color representing the time since last division. Purple and blue represent cells that have not divided for the past 10 h, and red represents the actively dividing cells. (G) A cross-sectional view of the local growth rate in the colony, with the color bar showing the values of local growth rate. A disc-shaped 'growth zone' is revealed by the red color at the bottom of the colony.

https://doi.org/10.7554/eLife.41093.009

The model shows a thin region at the periphery of the colony in which all cells are oriented in plane. This region governs radial growth and will be discussed more in the next section. Away from the periphery into the colony interior, more and more cells stand up vertically. The azimuthally averaged angle from the agar surface is plotted against the radial position in Figure 4B. However, the internal verticalization took some time to develop (Figure 4—figure supplement 2); appreciable fraction of cells (50%) picked up vertical orientation only when the radius reached 250 µm.

To characterize the spatial variation in cell orientation more quantitatively, we coarse-grain the local director fields n(r,t) (as described in Appendix A1.5) for the snapshot of Figure 4A. In Figure 4C, we plot the orientation of the azimuthally averaged director field, coarse-grained over boxes of size 4 µm × 4 µm over the rz-plane. We see that the orientation is vertical in the colony interior, but changes to be parallel to the colony surface in a transition zone of ~50 μm into the surface along the radial direction.

Next, we examine the coarse-grained velocity field vr=(vxr, vyr, vzr) whose azimuthal average is shown as arrows in Figure 4D. The velocity field points in the vertical direction throughout most of the colony, even at the top surface where cells are oriented parallel to the colony surface according to Figure 4C. Very close to the periphery in the bottom layer, the velocity field turns sideway; it is oriented planarly there and will be discussed below in the context of radial growth. As indicated by the length of the arrows, the vertically oriented velocity increases in magnitude away from the agar. This is illustrated by the plot of vertical velocity at different height z at the center of the colony, that is Vzz=vz(0,0,z), in Figure 4E. We see that Vz increases through a thin region of height HS10 μm. The vertical ascension speed is saturated for z>HS, meaning that above this thickness HS, cells move up steadily.

Another way to visualize the vertical growth of the colony is to show the ‘age’ of cells in a cross-sectional view (Figure 4F). In this plot, the age of a cell is defined as the time duration since the last division of the cell, with red being the youngest and purple being the oldest. We see that cells at the bottom of the colony are all young (red), indicating that the bottom layer is constantly dividing. In contrast, the oldest cells (purple) occupy the top/center region of the colony, and the next oldest age groups (blue, green, etc.) are located in different layers below the purple top region.

Together, the above results suggest a simple picture of the vertical colony growth: The cells are oriented vertically (except those close to the surface) and are pushed up by growing cells located within a 10 µm thick growth zone at the bottom; they stop dividing once pushed out of the growth zone. This picture is verified in Figure 4G, where the cross-sectional plot of the local growth rate shows a clear growth zone of ~10 μm (red region) confined to the bottom of the colony.

Vertical rise – quantitative analysis

This disc-shaped growth zone at the bottom of the colony may be intuitive, since cells at the bottom of the colony are in direct contact with the agar and hence have the best access to the nutrients. A planar growth zone is in fact required to support the observed linear increase of colony radius and height (during the period t = 12-24 hours in Figure 1): As the colony has the shape of a cone (Figures 1G and 3B), its volume is given by VcolonyR2HR3. Assuming that the increase of the colony size is due to a portion of cells growing at the maximal possible rate (λS) in a growth zone of volume Vgrowth(t), then ddtVcolonyVgrowth(t) leads to VgrowthR2, that is a disc. The thickness of this growth zone is of interest because it controls the vertical ascension speed. As the local growth rate is merely a 'readout' of the nutrient concentration according to Equation 3 in Materials and methods, we look into the penetration of nutrients into the colony, which determines the thickness of the growth zone. In Figure 5A, we plot the vertical nutrient concentration profile at the center of the colony, CctrzC(0, 0, z), at various times t during colony growth. In the linear growth regime (for t12 h), the profile Cctrz is essentially stationary. As shown in Figure 5—figure supplement 1 and Appendix A2.3, this stationary profile drops quadratically at small heights (i.e. close to the agar surface), and exponentially at larger heights (top of the colony), with the crossover between these two dependences occurring at the height scale HS such that CctrHS=KS, the Monod constant appearing in Equation 3; see Appendix A2.3. Since the local growth rate drops substantially where the nutrient concentration is below KS, we can take the value HS as the thickness of the vertical growth zone, leading to the vertical ascending speed: VHHSλS.

Figure 5 with 2 supplements see all
Vertical penetration of nutrients.

(A) The profiles of nutrient concentration Cctr(z)=C(0,0,z) along the z-axis at different times. (B) The profile Cctr(z) in the uniform scale vs. that in the rescaled z-axis, zλS. (C) The numerical results for the height velocity VH vs. the batch culture growth rate λS with a fixed cell division length ldiv (open circles) and variable ldiv (asterisk), respectively. The square root fit for the open circles (solid line) is given by the expression VH=5.5λS+0.6; the linear fit for circles with λS0.5 h-1 (dashed line) is given by the expression VH=3.2λS+2.9.

https://doi.org/10.7554/eLife.41093.012

Detailed analysis of Equations 1 and 3 in Materials and methods shows that the scale of the stationary profile is set by D+/λS; cf. Appendix A2.3. This is verified in Figure 5B where the stationary profile Cctr(z) is computed by repeating the simulation for different growth rate λS: the profile collapses into the same curve for different values of λS when plotted against zλS; see Figure 5—figure supplement 2 for the same profiles without rescaling in z. Given this scaling, we expect the thickness of growth zone to decrease as HS1/λS for faster growth (due to stronger nutrient depletion), leading to a sublinear dependence of the vertical ascending speed, VHλS. Our numerical result on the growth of vertical height is shown as open blue symbols in Figure 5C. The results are well fitted by the square-root dependence on λS; see the solid line. In Figure 1I, we attempted to test the dependence of the vertical ascension speed on growth rate experimentally, by growing colony in different carbon sources supporting different growth rates. Unfortunately, the most distinguishing regime of the predicted square-root relation, for λS<0.4 h1, is difficult to realize by changing carbon sources. However, if we just fit the data in Figure 5C for λS>0.5 h1, we obtain a weak linear dependence (dashed line) that resembles the experimental data in Figure 1I obtained. Note that the overall scale of the vertical ascending speed is 2-fold smaller in the simulation. This is attributed to the smaller nutrient concentrations used in the model compared to experiment, as will be discussed further below in the section of parameter dependence.

Radial expansion – a pictorial view

We first study the case mimicking glucose medium, corresponding to the simulation result shown in Figures 3 and 4. Since cells at the bottom grow substantially (Figure 4G), we plot the cell configuration for the bottom layer of the colony at t=20 h in Figure 6A; the same color code as Figure 4A is used, with vertically oriented cells shown in yellow and horizontally oriented cells in cyan. The periphery is seen to be largely cyan while the interior is more yellowish, suggesting that cells at the interior of the bottom layer are already oriented vertically, consistent with the cross-sectional view shown in Figure 4A. We again coarse-grain the local director field n(r,t) for the snapshot of Figure 6A. Figure 6B shows the planar projection of this director field in the bottom layer, where each bar indicates the average cellular orientation of cells in a region. We observe an annular region of ~100 μm in width near the periphery, where the director field has a significant in-plane component, directed mostly along the radial direction, except at the outermost boundary, where the director field has a great azimuthal component. Towards the inner boundary of the annulus, the in-plane component becomes smaller in magnitude. Interior to the annulus, the in-plane projection of the director field vanishes, confirming that they are largely oriented vertically.

Coarse-grained view of director, velocity and pressure fields in the bottom layer of colony.

(A) The bottom view of a simulated colony. Color scheme is the same as in Figure 4A. (B) Bars show the planar component of the coarse-grained director field at the bottom layer. (C) Arrows show the planar component of the coarse-grained velocity field at the bottom layer. Colors in (B) and (C) indicate the local pressure; see scale bar on the far right. The pressure is expressed in unit of P0=γsurf/w0 where γsurf is the surface tension, the main force underlying pressure build-up in our model.

https://doi.org/10.7554/eLife.41093.015

Next, we examine the coarse-grained velocity field v(r,t) for the bottom layer of cells shown in Figure 6A, with the x-y projection of v(r,t) shown as arrows in Figure 6C. We observe a narrow annulus of non-vanishing velocity field (arrows with finite length) at the outermost edge pointing radially outward; see also the side view provided in Figure 4D. Since the in-plane component of velocity vanishes inside the annulus (turning to vertical speed as shown already in Figure 4D), the driving force for radial colony expansion reside solely in the narrow annulus where cells are oriented planarly (Figure 6B). Just as the thickness of the growth zone determines the vertical ascension speed, here the width of the annulus determines the colony’s radial expansion speed.

So, what controls the annulus width? Or, equivalently, what determines the transition of velocity to the vertical orientation in the interior? Qualitatively, the difference between the peripheral and interior regions can be appreciated by looking at the coarse-grained pressure field P(r,t) experienced by the bottom layer, indicated by the color in Figure 6BC. This pressure is zero at the colony outer edge, and gradually builds up in the interior due to the physical growth of cells inside the closely packed colony. Where pressure is high, cells are oriented vertically and move vertically. This analysis thus suggests that increased pressure due to the physical growth of cells, which itself results from friction exerted by the substrate on the expanding cells, eventually forces cells to buckle and flow upward, manifested by the reorientation of cell directors in the vertical direction. Once the flow turns upward, pressure does not build up further due to the lack of friction with the agar surface. Since the upward flow is resisted by the surface tension, we conclude that pressure maxes out in this case at a level that is mostly determined by the surface tension. Below, we investigate quantitatively this buckling phenomenon.

Radial expansion – quantitative analysis

First, we examine the nutrient profile at the colony agar interface for growth on glucose. As can be seen from Figure 7A, the nutrient concentration is reduced underneath the colony. However, the actual concentration (Figure 7BC) is still much larger than KS of glucose uptake (dashed line in Figure 7BC), so that cells at the bottom do not experience nutrient depletion. In fact, at the colony periphery, nutrient concentration is close to the bulk level (Figure 7D).

Spatiotemporal nutrient profiles.

(A) The xz cross sectional view of the nutrient concentration inside the colony and in the agar, at time t=20 h. (B) The nutrient profile at the agar surface for different times. The nutrient concentration vs. time at the center (C) and at the periphery (D) of the colony at the agar surface.

https://doi.org/10.7554/eLife.41093.016

To elucidate the determinants of buckling, we plot in Figure 8A the azimuthal-averaged radial velocity profile Vr(Δr) for the bottom layer of cells, for a range of (signed) distances Δr into the edge of the colony; see Appendix 1 Equations (A1.5.1 and A1.5.2) for the definitions of Vr and Δr. This radial velocity profile, which is stationary for t12 h, is nearly zero in the colony interior, but increases almost linearly within a ~20 μm region at the outermost periphery of the colony. Since the radial expansion speed of the colony VR is simply Vr at Δr=0, we see that the width of this annulus together with the slope of Vr(Δr) set the radial expansion speed of the entire colony.

Figure 8 with 2 supplements see all
Physical characteristics near the outer periphery of the colony.

(A) The azimuthally averaged local speed of radial expansion Vr vs. the signed distance Δr from the edge of the colony (cf. Equations (A1.5.1) in Appendix 1) at various times. (B) The azimuthally averaged local height H vs. Δr at various times. (C) The azimuthally averaged local speed of vertical expansion Vz vs.Δr, averaging over t= 12, 16, 20 hr. (D) The azimuthally averaged local pressure P vs. Δr at time t= 20 h. As in Figure 6, pressure is expressed in unit of P0=γsurf/w0. (E) The azimuthally averaged height H vs. Δr at growth rate λS= 0.5 h-1 and 1.0 h-1. (F) The simulated colony horizontal expansion speed VR vs. the batch culture growth rate λS with a fixed ldiv(red open circles) and a growth-rate dependent ldiv (red closed circles) using dynamic friction. The dashed line that fits the open circles is given by VR=16.9λS+0.8; the solid line that fits the closed circles for λS0.5 h-1 is VR=5.7λSldiv-0.2; the dash dotted line fits the closed circles with VR= 22.1λS- 5.2. In these expressions, the speeds are in unit of μm/h and growth rate in  h-1. For comparison, we also include simulated VR vs. λS with a growth-rate dependent ldiv, for a model with static friction alone (see Equations (A1.4.6) of Appendix 1) between cell and agar (blue triangles).

https://doi.org/10.7554/eLife.41093.017

To understand what goes on in this peripheral region, we examine in Figure 8B the azimuthal-averaged height profile of the colony, h(Δr), which is also stationary for t12 h, with hΔr1 μm in the ~20 μm periphery region. This indicates that this periphery region is comprised of a single layer of cells lying horizontally on agar. In this monolayer region, the increase of Vr can be understood analytically, as we explain now. By mass conservation, the rate of local cell volume increase is balanced by the divergence of the velocity field, that is V=λ, where λ is the local mass growth rate (Klapper and Dockery, 2002). Through most of the monolayer region (except close to the inner edge), the vertical velocity Vz is negligible (Figure 8C). Hence, Vr satisfies

1rrVr'=λ.

Throughout the periphery region, the growth rate is essentially the maximal growth rate, that is λλS, since the nutrient concentration in this region is set by the boundary value Cs which well exceeds the Monod constant KS; cf. Figure 8—figure supplement 1. Solving the above equation in the annulus in the limit |Δr|R yields a linear dependence,

VrΔrVR+λSΔr

where we used the definition of the radial expansion speed VR= VrΔr=0. This solution is indicated by the red line in Figure 8A, which is in agreement with the numerical data, with a small discrepancy for small Vr attributed to the neglected vertical velocity at the inner periphery.

Given the linear radial velocity profile (cf. the previous equation) in the peripheral monolayer region, the width of this region Wb, defined as the largest value of |Δr| where VrΔr=0, sets the radial expansion speed since VRλSWb. We call this width the 'buckling width' since in the outer most ring region of the colony of size being this buckling width, cells form a monolayer, expanding with the speed VR, while the interior cells that are away from the colony edge by this buckling width flow up vertically; cf. Figure 8—figure supplement 2. The magnitude of the buckling width is set by the radial forces exerted on the monolayer of cells. As these cells grow outward, they experience frictions from the agar substrate as well as surface tension that holds them down flat. These two forces lead to the accumulation of pressure, which is built up from the periphery. Figure 8D shows the azimuthally averaged pressure P(Δr) for the bottom layer of cells. At the inner edge of the monolayer region, pressure reaches a critical value that surface tension can no longer hold cells in a single layer. There, some cells buckle into the vertical dimension, leading to vertical flow of cells and forming multiple layers (Figure 8—figure supplement 2), alleviating the further build-up of pressure. It is interesting to observe that this buckling phenomenon is already evident early on during transition from monolayer growth to multiple layers, as shown in Figure 3D. The 20 µm annulus of monolayer at the periphery is set soon after the initial buckling transition, at around t=8 h (Figure 3D), setting the pace of radial expansion.

Quantitative details of the buckling transition depend on the form of the cell-agar friction. Two types of friction have been used in the cell-modeling literature, one which depends linearly on the cell-agar velocity, known as viscous or static friction (Farrell et al., 2013; Ghosh et al., 2015), and the other which saturates to a constant set by the magnitude of the normal force (in this case, resulting from the surface tension). The latter is referred to as dynamic friction; see Appendix 1.4. The two forms can be distinguished by comparing the buckling width Wb at different radial expansion speed VR: Static friction would increase for increased VR, leading to decreased buckling width, whereas dynamic friction would not be affected by the radial expansion speed. Experimentally, we characterized VR in sugars supporting different growth rates λS, and VR is seen to increase linearly with λS (Figure 1I), suggesting a constant Wb, and hence the applicability of dynamic friction. This dependence is tested by running simulations with the dynamic friction form (Equations 7a and 7b in 'Computational Model') for different growth rate λS. The buckling width Wb is indeed not dependent on λS (Figure 8E), and the radial expansion speed VR is indeed linear in λS (open red symbols and dashed red line, Figure 8F). In contract, static friction leads to a much weaker dependence of VR on λS (blue triangles in Figure 8F).

The linear dependence on λS seen in the experimental data in Figure 1I (red symbols) however exhibits a noticeable horizontal offset. This offset likely results from an additional effect we have not included into the model so far: The size of the cells is dependent on their growth rate, with faster growth rate being longer and wider (Jun and Taheri-Araghi, 2015; Nanninga and Woldringh, 1985; Taheri-Araghi et al., 2015). By repeating the established dependence of cell size on growth rate (see Equation (A2.2.3) in Appendix 2) for different values of λS, we recover a nonlinear dependence of VR on λS (Figure 8F, filled red circles and solid red line). Note that a similar horizontal offset is obtained as the experimental data in Figure 1I if we do a linear fit using the data with λS>0.5 h1 (dotted red line). On the other hand, the growth-rate dependence of cell sizes has no noticeable effect on the vertical ascension speed (filled blue symbols, Figure 5C) since the growth zone thickness HS1/λS does not depend on ldiv (Appendix 2.3).

Parameter dependance

The preceding analysis shows that the vertical expansion speed of the colony depends on the thickness of the vertical growth zone which is set by the nutrient penetration depth, while the radial expanding speed depends on the width of the monolayer annulus which is set by the onset of the buckling transition but not the nutrient. The sizes of these growth zones are therefore dependent on the magnitudes of the physical parameters in different ways: We expect changing the cell-agar friction to affect the onset of the buckling transition and hence the radial expansion speed VR, but not the vertical ascension speed VH. Conversely, we expect changing the nutrient concentration Cs to change the vertical nutrient penetration length and hence VH, but not VR. These expectations are indeed reproduced by the full colony simulation using different parameter values from the ones discussed so far, with the nutrient concentraiton Cs doubled in Figure 9A (only VH increased), and with the cell-agar friction quartered in Figure 9B (only VR increased). These predictions are further tested experimentally, by varying the glucose concentraton used in the agar (Figure 9C), and by repeating experiments in reduced agar densities (Figure 9D) which we expect to reduce the cell-agar friction. The observed changes are very much in line with the expectations of the model shown in Figure 9AB. These results serve to validate the very important qualitative results of our study, that radial grow of the colony is not limited by nutrients while the vertial growth is limited by nutrients.

Parameter dependence of colony growth characteristics.

Simulation results using the full model with 2x increase in glucose concentration (panel A) and 4x decrease in all frictional parameters (panel B) for VR (red bars) and VH (blue bars). Specifically, in panel (A) we fix the friction at a high level (with μca=0.8, μcc=0.1, γcc,t=10000 μm1h1, and hran=0.1 μm), and use Cs=0.5 mM as the lower glucose concentration, Cs=1.0 mM as the higher glucose concentration. In panel (B), we fix the glucose concentration at the lower level (Cs=0.5 mM), and vary the friction, from the higher value of μca=0.8, μcc=0.1, γcc,t=10000 μm-1h-1hran=0.1 μm used in (A) to the lower value of μca=0.2, μcc=0.025, γcc,t=2500 μm-1h-1hran=0.025 μm. The corresponding experimental results are shown in panels C and D: In (C), glucose concentration was varied with agar density fixed at 1.5%. In (D), agar density was varied with glucose fixed at 0.2% (w/v). The data for VH in panel C is consistent with a square root dependence on nutrient concentration (blue line) expected from the basic analysis in Figure 5.

https://doi.org/10.7554/eLife.41093.020

We note that the actual values of radial and vertical expansion speeds obtained (VR=17.2 μm/h and VH=5.8 μm/h), for the standard parameter set used (Supplementary file 1-Tables S2–S4) through the bulk of the study, are approximately 2x lower than the range of values obtained in experiments. The results of Figure 9AB show that the experimental range could be obtained simply by adjusting the combinations of parameters. We did not do that – the parameter set giving smaller VR and VH was chosen due to computational constraints: Higher nutrient concentrations requires longer computational time to reach the linear steady state due to the larger nutrient penetration depth. Similarly, lower cell-agar friction would lead to colony spreading too rapidly in the radial diretion, thus requiring larger simulation sizes and hence again longer computational time. Their combination becomes difficult to investigate at the level of details done in this study. The particular values of frictional coefficients in the standard parameter set have been chosen so that the colony retains similar aspect ratio as observed in experiments, but with both VR and VH being about half of the experimentally observed values for growth on glucose medium. As computing power continues to increase, these models should soon be able to reach sizes comparable to realistic colonies with realistic parameters.

Discussion

In this work, we presented a detailed quantitative study of the growth of a bacterial colony on hard agar surface starting from a single cell. For non-motile bacteria incapable of producing extracellular polysaccharides, the colony is driven primarily by the force of their own growth. Key factors involved are nutrient diffusion, mechanical interactions between cells, friction between cell and agar, and the surface tension holding the cells to the agar. We developed a continuum model for nutrient diffusion and implemented it with a multi-resolution numerical technique. With a discrete agent-based model, we captured mechanical interactions, including elasticity and dynamic friction. Most importantly, the surface tension of the liquid in the colony is implemented by introducing a restoring force on cells protruding from a smoothened colony surface.

Our model is able to capture quantitatively some of the characteristic features observed for bacterial colony growth, including the conic shape of the colony, the linear expansion of colony radius and height, and both the linear and sublinear dependence of the speed of radial expansion and that of vertical expansion, respectively, on the cell growth rate. The model makes a number of important predictions on the expanding colony as summarized in Figure 10: The growth zone is predicted to be disc-like and extended throughout the bottom of the colony, contrary to common belief (see below). Radial growth is driven by cells at the outer perimeter of the growth zone; these cells are predicted to form a thin layer, oriented parallel to the agar due to the downward pull of surface tension, with the width of the region determined by the onset of the buckling transition (which occurs when radial compression due to cell-agar friction overwhelms the surface tension). In the colony interior, cells are predicted to orient vertically and are mainly pushed upward by elongating cells in the bottom growth zone.

A schematic summary of key mechanisms in the growth of an E.coli colony.

After an initial, exponential monolayer growth, buckling occurs at the center of the colony. Cells then grow actively only in the bottom layers (red vertical arrows) whose thickness (HS) is determined by the nutrient penetration level (dashed blue line). Cells lying above them are passively pushed up. Throughout this yellow triangular region, cells are oriented vertically. Near the colony edge (cyan region), the cells are oriented planarly and grow outward (horizontal red arrow) continuously in a spread mode to expand the colony in the radial direction. The width of this annulus (Wb) is determined by mechanical effects arising from the surface tension which pulls the thin layer of cells into the agar, and cell-agar friction which builds up the pressure from the outer edge of the layer, eventually causing buckling at an inner radius where cells transition to the vertical orientation (the green region). These two characteristic parameters, HS and Wb, set the speeds of radial and vertical expansions, VR and VH, respectively, as shown in red. The growth rate dependence of these parameters is shown in blue.

https://doi.org/10.7554/eLife.41093.022

Capturing all these behaviors within a single model and with a fixed set of parameters is a non-trivial task despite the seeming simplicity of this problem. Many aspects of our model are taken from what are commonly adopted in the extensive literature devoted to this class of problems over the past decade (Boyer et al., 2011; Cole et al., 2015; Farrell et al., 2013; Ghosh et al., 2015; Grant et al., 2014; Jayathilake et al., 2017; Rudge et al., 2013; Rudge et al., 2012; Volfson et al., 2008). These include the basic modeling of metabolism and cell growth (Cole et al., 2015; Farrell et al., 2013; Rudge et al., 2012), and the use of Hertzian elasticity to describe cell-cell elastic interaction (Boyer et al., 2011; Farrell et al., 2013; Ghosh et al., 2015; Grant et al., 2014; Volfson et al., 2008), all incorporated as computational power increases to reach ever increasing colony sizes (Cole et al., 2015; Rudge et al., 2013; Rudge et al., 2012). Unique to our study is the treatment of mechanical interactions, specifically friction and cell-level surface tension, which we believe are at the root of all behaviors described above, including the forms of radial and vertical colony growth. A key result of our study is that the linear radial growth is driven by the growth of a thin layer of radially oriented cells located at the colony periphery, whose width is determined by mechanical buckling. Although the linear radial expansion of bacterial colonies has been known for about 50 years (Pirt, 1967), for a long time this was attributed to a ring-shaped growth zone at the outer colony periphery due to nutrient diffusion (Lewis and Wimpenny, 1981; Pirt, 1967; Wimpenny, 1979). Only quite recently has the notion been made that mechanical effects might also lead to linear radial growth (Farrell et al., 2013; Su et al., 2012). (Su et al., 2012) showed experimental results that implicated the interplay of forces in the radial expansion of colonies. (Farrell et al., 2013) proposed mechanical effects as a colloquial rationalization of numerical results generated by toy models with unrealistic details, for example a ‘gravity-like’ adhesion force acting on all cells in the colony. In our study, the adhesion of cells to the agar surface is provided by the surface tension of the liquid surrounding cells in the colony. We introduce a novel cell-level model of surface tension which acts only on cells at the colony surface, distinct from common models of surface tension which depends on the macroscopic curvature of the colony surface and cannot describe thin layers. It is this unique surface tension model that enables us to capture the dynamics from the initial single-layer cell growth, through buckling, to the growth of a macroscopic colony. This cell-level surface tension, responsible for pressing cells into the agar thereby generating friction that eventually causes buckling, cell reorientation and vertical colony growth, is thus the source of all mechanical interactions in the colony. A strong, uniform force such as the ones used in (Farrell et al., 2013) would lead to artificially flattened colonies, especially at the colony center where the height is the highest, since the force is proportional to the height in that model.

We regard the characterization of colony growth for different nutrients (which give rise to different cell growth rates) as a unique contribution by our study. The knowledge of the dependence of colony growth on cell growth allows us to discriminate different models of colony growth. As an example, an important component of our model that makes a quantitative difference to the outcome is the form of the friction used. Viscous drag (i.e., friction proportional to the velocity difference) is the form adopted in most models of cell dynamics (Farrell et al., 2013; Ghosh et al., 2015; Rudge et al., 2012). We instead adopt a form commonly used in modeling granular solids (Brilliantov et al., 1996; Cundall and Strack, 1979; Kuwabara and Kono, 1987; Shäfer et al., 1996). It involves a static friction depending on relative velocity, capped by a dynamic friction which is independent of the velocity. This form, introduced in one of the first models of 2D colony growth (Volfson et al., 2008), exerts a pressure which is independent of the speed of radial expansion, leading to a growth rate-independent buckling width and hence a radial expansion speed that is proportional to cell growth rate, in agreement with our experiments. In contrast, a model based on static friction would have the buckling width reducing with increasing cell growth rate, giving a sublinear dependence of radial expansion speed on cell growth rate which is not compatible with the data in Figure 1I. Indeed, in a model with static friction alone, a much weaker growth-rate dependence of radial expansion speed was obtained (Figure 8F blue triangles). Along a different line, Fisher-Kolmogorov (FK) dynamics has been used as a phenomenological model to describe radial colony expansion, and has been successful in describing certain spatial patterns formed in growing colonies (Cao et al., 2016). However, FK dynamics would predict a square-root dependence of the radial expansion speed on the cell growth rate (Fisher, 1937; Kolmogorov et al., 1937), which will need to be reformulated to conform to the observed dependences.

In addition to the well-known linear radial growth, the linear vertical growth of the colony is dissected for the first time qualitatively here since it was first reported (Lewis and Wimpenny, 1981; Wimpenny, 1979). Our analysis shows that the vertical expansion speed is limited by the depth that nutrient can penetrate upward into the colony from the agar. Accompanying our result of vertical growth is the predicted vertical orientation of cells in the colony interior, which transitions from the radial orientation at the outer periphery (i.e., the monolayer zone colored in cyan in Figure 10).

Cell verticalization has been observed experimentally for Vibrio parahaemolyticus (Enos-Berlage and McCarter, 2000) and for Vibrio Cholerae (Beroz et al., 2018; Yan et al., 2016). In both cases, vertical orientations could be seen already for very small bacterial colonies, possibly due to their production of extracellular polysaccharide substance (EPS). In this work, verticalization is predicted to occur for plain bacterial colonies as well, without the need of any EPS, but at much larger colony sizes. We have not been able to observe verticalization directly for our colonies due to multiple scattering associated with very dense colonies we are studying. This is left as a challenge for future studies.

In our model, verticalization results from an interplay among colony surface tension, cell-agar friction and the physical force of expansion due to cell growth. (Beroz et al., 2018) also introduced a discrete model to describe cell verticalization. In their model, verticalization resulted from a similar mechanical instability due to the interplay between in-plane compression force and cell-agar adhesion. Due to the different energy barriers against verticalization, the length scales of verticalization between our model and that of (Beroz et al., 2018) are very different: The colonies in Beroz et al. spread very slowly radially (~3μm/h), and verticalization occurs at a colony radius of ~5-10 μm. Colonies in our model spread much faster (~14 μm/h), and substantial verticalization occurs at a radius of ~250 μm; see Figure 4—figure supplement 2.

Although we have restricted our study to colonies growing in rather simple conditions, insight from our model can be readily used to make qualitative predictions in a variety of other conditions. Generally, we expect the radial expansion speed to be controlled by the buckling width and vertical expansion speed be controlled by the thickness of the growth zone. Thus, if agar hardness or ambient humidity is changed, the effect on air-liquid surface tension is expected to affect the buckling width and the ratio of the radial and vertical expansion speeds, hence changing the colony aspect ratio. Also, during later stages of colony growth when oxygen becomes limiting in the colony interior, the obligatory excretion of large amounts of fermentation product associated with anaerobiosis is predicted to lower the pH in colony interior and thereby slow down vertical colony growth while not affecting the radial growth. Our observations shown in Figure 1H and Figure 1—figure supplement 1 are in qualitative agreement with the expectation. A quantitative study of this late regime (t>24 h for the growth condition used in Figure 1H) requires a much more detailed model of anaerobic metabolism, pH effect, and growth transition kinetics, well beyond the scope of the current study, and will be reported elsewhere. Note that recent series of colony-based microbial range expansion studies (Hallatschek et al., 2007Hallatschek and Nelson, 2010; Korolev et al., 2012), which involve much larger colony sizes and longer periods of colony growth, are likely in this late regime where vertical growth has ceased. Nevertheless, the radial expansion of these large colonies may still be governed by the same factors discussed in this work.

While our work is exclusively on bacterial colonies without EPS, key results we learned from this study shed light on the more complex dynamics of heterogeneous biofilms. First, we establish that the radial growth of our colonies is not limited by nutrient as commonly believed, but by the interplay of surface tension and cell-agar friction (Figure 9). Given that biofilms have typically much lower bacterial densities, nutrient limitation will be even less of a problem. Also, EPS secreted by the bacteria could modify both the surface tension and cell-agar friction to control the radial expansion speed. Second, nutrient supply is limiting for the vertical growth of our colonies (Figure 9AC). This becomes less of a problem for the loosely packed biofilms. Moreover, biofilms are said to form channels in their interior (Wilking et al., 2013), which would further alleviate the supply of nutrient, thereby allowing for faster vertical expansion. Finally, verticalization of cells in the interior, which is important for vertical growth but occurs at rather large colony sizes according to our model (Figure 4—figure supplement 1), also occurs in biofilms but at much smaller colony sizes (Beroz et al., 2018; Enos-Berlage and McCarter, 2000; Yan et al., 2016). While the precise nature of the forces driving verticalization may be different in the two cases, the underlying origins may be similar — mechanical instability due to in-plane compression resulting from colony expansion and cell-agar friction. In light of these comparisons, we see that the additional ingredients provided by biofilms enable the colonies to expand faster both horizontally and vertically.

The model presented here, with results quantitatively comparable to experimental data, can be used to interpret large-scale data being generated by high-throughput colony growth assays to track the growth of different strains in different conditions (Takeuchi et al., 2014). Our model can be used as a launching pad, not only to include the more complex effects of metabolism and cell growth mentioned here, but also other factors such as extracellular matrix to allow the simulation of biofilms, and multiple interacting species to explore microbial ecology in compact space. Finally, it will be an interesting challenge to develop coarse-grained hydrodynamic models that incorporate the unique features of surface tension and dynamic friction discussed here, and capture the radial and vertical colony growth characteristics, both their temporal behaviors and their dependences on cell growth rates and other environmental factors.

Materials and methods

Experimental methods

Bacterial strain

The strain of E. coli K12 used in all the experiments reported in this work, EQ59, was derived from NCM3722 (Lyons et al., 2011), with deletion of the motA gene to remove bacterial motility and harboring constitutive GFP expression. We note that biofilm formation is highly suppressed in NCM3722, as acquired nonsense mutations within both the bsg and csg operons prevent the synthesis of extracellular cellulose and curli needed to support biofilm (Lyons et al., 2011; Serra et al., 2013).

To make strain EQ59, we cloned the gfp gene from pZE12G (Levine et al., 2007) into the KpnI/BamHI sites of the plasmid pKD13-rrnBT:Ptet (Klumpp et al., 2009), yielding the plasmid pKDT_Ptet-gfp. The fragment 'kmr:rrnBT:Ptet-gfp' present in pKDT_Ptet-gfp was PCR amplified, gel purified and then electroporated into EQ42 cells (Klumpp et al., 2009), expressing the λ-Red recombinase. The cells were incubated with shaking at 37°C for 1 hour and then applied onto LB+Km agar plates. The Kmr colonies were verified for the 'kmr:rrnBT:Ptet-gfp' substitution for the 67 bp intS/yfdG intergenic region between 117th and 51st nucleotides relative to the start codon of yfdG by colony PCR and subsequently by sequencing. The chromosomal region carrying 'kmr:rrnBT:Ptet-gfp' in EQ42 was then transferred to EQ54 (that is NCM3722ΔmotA) (Kim et al., 2012) by P1 transduction, yielding strain EQ59, in which the gfp gene is constitutively expressed in the absence of TetR.

Growth medium

Phosphate-buffered media (N- C-) was used for both batch culture and colony growth as described in Csonka et al. (1994). Various carbon sources were used as specified in Supplementary file 1-Table S1. The concentration of all carbon sources used was 0.2% (w/v) unless otherwise specified. 10 mM of NH4Cl was added as the sole nitrogen source. The agar concentration used was 1.5% (w/v) unless otherwise specified. 20 mL of molten agar gel was poured into 60 mm diameter dishes to a final thickness of approximately 7 mm, and allowed to cool at room temperature. Agar plates were sealed in plastic and stored at 4°C until use.

Cell growth

Batch culture growth was performed in a 37°C water-bath shaker (220 rpm). Cells from a fresh colony in a LB plate were inoculated into LB broth and grown for several hours at 37°C as seed cultures. Seed cultures were then transferred into the desired minimal medium and grown overnight at 37°C as pre-cultures. For batch culture growth rate measurements, overnight pre-cultures were diluted to OD6000.01 in the same minimal medium and grown at 37°C as experimental cultures. After two doublings, OD measurements were taken at various time over a 10-fold increase (i.e., from 0.04 to 0.4), and the growth rate was determined from a linear fit of ln(OD) vs. time.

Colonies were seeded on the agar gel as single cells. The pre-culture (prepared as above) was diluted to OD60010-6. 10 μL of culture (containing approximately 10 cells) was spread over pre-warmed plates and transferred immediately to a 37°C incubator for growth. Petri dishes remained covered at all times, except during periodic measurements with a confocal microscope, in order to minimize moisture loss.

Microscopy

Colonies were imaged with a Leica TCS SP8 inverted confocal microscope placed within an incubated box at 37°C. Samples were grown in covered petri dishes stacked on one side of the box. Each was moved to the microscope objective for periodic measurements. They were immediately covered once measurement was done. For the measurements, the dishes were uncovered and measurements were taken from the top (air) side to obtain a complete 3D image of the colony. GFP was excited with a 488 nm diode laser, and fluorescence was detected with a 10×/0.3 objective and a high sensitivity HyD SP GaAsP detector. For a large colony, an xy-montage was created and stitched together to form a single 3D image using the ImageJ Grid/Collection Stitching plugin.

Image analysis

The colony shape was obtained from the 3D confocal image using custom Matlab software. Under aerobic conditions, the bacterial fluorescence was spatially uniform near the top surface of the colony, and the surface height, h(x,y), could be reconstructed by simply thresholding the intensities: for each (x,y) position, the height was defined by the top pixel whose intensity was greater than the threshold. To account for the fact that fluorescence varied somewhat with growth conditions (sugar, agar concentration, etc.), this threshold was rescaled by the maximum fluorescence of the colony for each condition.

Furthermore, to capture the radius of the single- and multi-layers at early time of colony development (Figure 1A–E), we analyze the image intensity of the colony as the follows: for each stencil of 5×5 pixels centered at pixel (i,j), we count the number of pixels whose intensity is above a threshold, and call it ni,j. Pixel (i,j) is assigned as type 1 if 16>ni,j3, and as type 2 if ni,j16, indicating the pixel belonging to single- or multi-layer region, respectively. We then estimated the inner radius rinner and outer radius router of the colony by the formulas rinner=rμm/px Npx2/π and router=rμm/px(Npx1+Npx2)/π, where Npx1 and Npx2 are the total numbers of pixels of type 1 and 2, respectively, and rμm/px0.84 is the ratio of µm per pixel in our confocal image.

Colony growth curves

Colony growth was monitored by measuring an individual colony at intervals of 1—4 hours. The radial growth curve, R(t), was extremely reproducible from colony to colony on the same agar plate and from day to day on different plates, up to a small offset in time, tl, reflecting a variable lag time, of up to two hours before colony growth began. To monitor the colony growth over long periods of time, we started identical colonies at seed times ts separated by several hours. Growth curves extending over a period of multiple days could be obtained by stitching together R(t-tl-ts) at times where they overlapped. This stitching procedure is illustrated in Figure 1—figure supplement 1. For example, in Figure 1—figure supplement 1A, there are three different symbols: triangles, squares, and circles. Each symbol represents data from one colony. They are seeded several hours apart and are plotted together with respect to their respective starting time. The data thus shows that the colony development is highly repeatable and can be put together to reconstruct the overall dynamics which spans a long period. In most cases, at least three separate colonies are measured concurrently for each (short) time span, and three separate time spans were stitched together in a series.

Computational model

Continuum model of nutrient dynamics

We assume that the growth of cells in the colony is limited by a single type of nutrient (the carbon source), and use a continuum scalar field C=C(r,t) to represent the nutrient concentration at a spatial location r=(x,y,z) and time t. Agar, which contains the nutrient and which cannot be penetrated by cells (at the dense concentrations used in out experiments), is confined to the region z<0, while cells grow on top of the agar in the region z>0, and bounded by the colony surface Γ01 to be defined below; see Figure 11. Nutrient diffuses in the two compartments, agar and colony, according to the diffusion equations

Schematic of the computational box and different regions in the model of simulation.

The computational box is Ω=(L, L)×(L, L)×(a, b), where all L, a, and b are positive numbers in the units of length. This box is divided into the air region Ω0, colony region Ω1, and agar region Ω2=(L, L)×(L, L)×(a, 0), respectively. See Supplementary file 1-Table S5 for typical values of L, a, and b used in our simulations. The colony surface or colony-air interface Γ01 separates the colony from air. The plane z = 0 in the computational box is divided into two parts. One is the interface that separates the colony from agar, and is denoted by Γ12. The remaining part, denoted Γ02, separates the air from agar. Note that, since the bacterial colony grows with time t, all the air region Ω0, the colony region Ω1, the colony-air interface Γ01, and the colony-agar interface Γ02 depend on time t.

https://doi.org/10.7554/eLife.41093.023
(1) tC=D+ΔC ρλ/Y for z>0,
(2) tC=DΔC for z<0,

with the distinct diffusion coefficients D+ in the interstitial space between cells in the colony above the agar, and D- inside the agar. The second term on the right-hand side of Equation 1 describes the rate of nutrient consumption by growing cells. Here, ρ= ρ(r,t) is the local cell mass density (total mass of cells in a unit volume of space) and Y is the yield factor. For simplicity, we shall approximate the spatially and temporally varying cell mass density ρ= ρ(r,t) by a constant value ρ0. Above the spatial scale of a few cell lengths, the spatial variation in ρ is <5% within the colony; see Figure 3—figure supplement 1. The upper boundary of the colony (Γ01 in Figure 11) is defined by thresholding the density; see Appendix 1.2. The local mass growth rate λ=λ(r,t) is given by Monod kinetics

(3) λr,t=λSCr,tCr,t+KS ,

where λS is the batch culture growth rate for cells in a medium saturated with some sugar S, and KS is the Monod constant for the sugar S. At the (mean) interface (z=0) between the colony and agar substrate, we have the continuity of the nutrient concentration and its flux:

(4) C=C+andDzC=D+zC+attheagarcolonyinterfacez=0,

where the symbols C- and C+ indicate the nutrient concentration on the agar side (z<0) and colony side (z>0), respectively. Equations 1–4 are supplemented by boundary conditions imposed on the boundaries of a computational region comprising of both the colony and agar regions. We impose the flux-free boundary condition on the parts Γ01, Γ02, and Γb, and the Dirichlet boundary condition C=Cs on the lateral wall of agar region Γs; cf. Figure 11. The parameter Cs mimics the nutrient concentration far away from the colony. It is one of the key parameters in our study.

Discrete model for cell growth, division and movement

In addition to modeling the nutrient as a continuum, we use a discrete, agent-based model to describe the growth, division, and movement of cells, as well as the interactions of cells with each other and with the environment. In this agent-based model, each E. coli cell is represented by a sphero-cylinder, comprised of a cylinder with hemispherical caps of diameter (also called cell width) w0 on its two ends; see Figure 12A. For a given cell i at a given time t, we use a position vector ri(t) to denote the center-of-mass of the cell, a unit vector ni(t) to denote its orientation (the direction along the cylindrical axis), and li(t) to denote the length of the cylinder between the two centers of hemispherical caps. Each cell starts off with the same cylindrical length l0. We assume that during cell growth, the width w0 is fixed, and the cylinder length of a cell increases at a rate λ~t. We call this the cell elongation rate. It is proportional to the mass growth rate λt of the cell with a geometrical proportionality factor σλ~= σλ; see Appendix 1.3.

Schematic view of cell growth and division.

(A) A sphero-cylinder model of an E.coli cell. Here, w0 is the diameter of each of the hemispheres, p and q are the centers of these hemispheres, l=|pq| is the cylindrical length of the cell, n=(pq)/l is the unit vector along the cylindrical axis of the cell, and rc is the center of the cell. (B) Cell division. Once the cylindrical length of a cell reaches a critical value ldiv, the cell divides into two daughter cells. The two centers of hemispheres of the mother cell become the centers of hemispheres of the daughter cells. Each of these two daughter cells has the cylindrical length l0 with fluctuations, where l0 is a constant cylindrical length for any new born cell and any initial cell in the simulation. Fluctuations of angular velocities are also introduced for the daughter cells; cf. Appendix 1.3.

https://doi.org/10.7554/eLife.41093.024

The mass growth rate is calculated based on the nutrient concentration at the center ri(t) of the cell i at time t, that is λit=λ(rit,t), according to Equation 3. The growth of cylindrical length li(t) is then given by the growth equation

(5) ddtli=σλrit,tli(t).

Once the cylindrical length li(t) reaches a critical value ldiv, the cell divides into two daughter cells with cylindrical lengths being l0 with small fluctuation; see Figure 12B and Appendix 1.3. For different growth media supporting different growth rates λS, the value of ldiv is in general growth-rate dependent (Jun and Taheri-Araghi, 2015; Taheri-Araghi et al., 2015), the consequences of which are discussed above in 'Radial expansion – quantitative analysis'.

The position and orientation of cell i change according to its velocity vi and angular velocity ωi, which follow Newton’s second law

(6) Mitvi=FinetandIitωi=Tinet

where Mi and Ii are the mass and moment of inertia of the cell, and Finet and Tinet are the net force and net torque, respectively, exerted on that cell. As cells grow, divide and move, the colony region (defined by the part of boundary Γ01 in Figure 11) expands. The nutrient concentration in the new domain requires an update by solving the boundary-value problem of Equations 1–4 again.

Discrete model for interaction forces

The net force Finet and torque Tinet exerted on cell i arise from (a) cell-cell mechanical interaction, (b) cell-agar interaction (if the cell touches the agar surface), (c) cell-fluid interaction, and (d) surface tension (if the cell is on top of the colony); cf. Figure 2. Below we briefly describe each component used in our model. Details are provided in Appendix 1.4.

(a) Cell-cell interaction

In the interior of a colony, two cells interact only if they are in direct physical contact, characterized here by the overlap δcc in their sphero-cylinder cell boundaries; see Figure 13A. At the point of contact, the cell-cell interaction force Fcc is decomposed into the normal and tangential components, of magnitudes Fcc,n and Fcc,t as defined in and Appendix 1.4a. The normal force includes the Hertzian elasticity force with magnitude Fcc,elasw0δcc3/2 (Hertz, 1882; Johnson, 1985). Additionally, the normal and tangential force each has a dissipation component, of magnitude Fcc,disp, n and Fcc,disp,t, respectively, describing the effect of friction against cell movement. In the cell modeling literature (Farrell et al., 2013; Ghosh et al., 2015; Rudge et al., 2013; Rudge et al., 2012), these dissipation forces are often taken to be viscous. (We shall include such viscous force in (c) below.) In our model, we found it necessary to further include static and dynamic friction as described below.

Schematic view of cell-cell and cell-agar interactions.

(A) Cell-cell interaction. Two cells, centered at rc and rc', respectively, are in contact with each other. The shortest distance between the central cylindrical lines of the two cells is d=aa' with a and a' two points on the cylindrical central lines of these two cells, respectively. The amount of the overlap of these two cells is δcc=w0-d. We shall denote by rcc the center of the line segment connecting a and a'. The total interaction force, exerted at center rcc, is the sum of the normal force Fcc,n in the normal direction ncc=(a-a')/d and the tangential force Fcc,t in a direction orthogonal to ncc that is determined by the relative velocities of these two cells; see the details in Appendix 1.4 on force calculations. (B) Cell-agar interaction. A cell, centered at rc, touches the agar surface that has the mean position at z=0 and the roughness hran (the maximum fluctuation around the mean), with δca the amount of the overlap of the cell and agar. If the center of the hemispherical cap of the cell corresponding to the end that dips into the agar is p=xa, ya, za, then δca=w0/2 za. Denote rca=(xa, ya, za-w0/2), which is the midpoint of the line segment along the vertical line passing through the point p between z=0 and z=-δca. The total cell-agar interaction force Fca, exerted at the center rca, is the superposition of a normal force Fca,n in the vertical direction and the tangential force Fca,t in a direction along the xy plane that is determined by the velocity of the cell; cf. Appendix 1.4.

https://doi.org/10.7554/eLife.41093.025

We follow standard models of granular solids (Brilliantov et al., 1996; Cundall and Strack, 1979; Kuwabara and Kono, 1987; Shäfer et al., 1996), first introduced to cell modeling by Tsimring and his collaborators (Volfson et al., 2008). To model static friction, we adopt a fictitious drag force whose normal and tangential components, Fcc,disp, n and Fcc,disp,t, respectively, are taken to be proportional to the normal and tangential components of the relative cell-cell velocity, vcc,n and vcc,t. We use Fcc,disp,nδcc1/2vcc,n and Fcc,disp,tδccvcc,t, where the additional dependences on the overlap δcc captures the dependence on contact area; see Figure 13A. To implement dynamic friction, we cap the tangential dissipation by the static yield criterion, that is Fcc,disp,tmax=μccFcc,elas, where μcc is the dynamic frictional coefficient; see Appendix 1.4a. Thus, the full cell-cell interaction force is given by

(7a) Fcc,n=Fcc,elas+Fcc,disp,n ,
(7b) Fcc,t=minFcc,disp,t , μccFcc,elas.

As we see in ‘Radial expansion – quantitative analysis’, a dynamic friction form imposed by Equation 7b provides a natural explanation of the experimental observation that the radial velocity of the colony is independent of the cell growth rate.

(b) Cell-agar interaction

The force exerted on a cell in contact with the agar, Fca, can be similarly calculated as sketched in Figure 13B. The same forms of the elastic and frictional forces are used as in Equations 7a and b. Note that to break the planar symmetry and facilitate buckling of cell layer into the vertical direction, we introduced certain roughness to the agar surface, characterized by the roughness parameter hran which is the maximum fluctuation of the agar surface around its mean (z=0).

(c) Cell-fluid interaction

A cell also interacts with the surrounding fluid and experiences viscous drag. This is given by the Stokes drag force Fvisc, which is proportional to the velocity of the cell v. We note that in high-density colonies such as those studied here, dissipation due to viscous drag is significantly less than the cellular friction force.

(d) Surface tension

Surface tension is a critical factor determining the dynamics of an expanding colony. It is frequently treated as a property of a composite fluid comprising of cells plus the surrounding fluid (Grimson and Barker, 1993; Zhang et al., 2008). Alternatively, the liquid phase is ignored altogether, and surface tension is assumed to arise from attractive interactions between the cells themselves (Ben-Jacob et al., 2000). In both cases, surface tension reflects the curvature of the macroscopic colony profile. However, such coarse-grained treatments of surface tension cannot describe the initial layer-by-layer growth of the colony arising from buckling (Figure 1A–E), nor can they capture thin layers surrounding the periphery of large colonies (Figure 14). In order to capture these effects, we endeavor here to model the surface tension experienced by cells in a colony as a boundary force, that is as force experienced by discrete cells at the colony boundary due to increased surface tension of the continuum liquid these cells are immersed in.

Schematic view of surface tension acting on cells at the colony boundary.

(A) a snapshot of part of a growing colony from a typical computer simulation. The red curve defines the macroscopic colony-air interface which is used in updating the nutrient profile; cf. Γ01 in Figure 11 and Appendix 1.1. Cells on the top of the colony are held down by the surface tension force. (B) A sequence of configurations of colony during the early stage of growth. Top: Initial cells grow exponentially and form a monolayer. All cells in this layer are held down by the surface tension force. Middle: As more cells are born in the monolayer, the frictions between these cells and the rough agar surface increase. The competition between such frictions and the surface tension force that pulls down cells leads to an accumulation of the lateral pressure in the cells. The monolayer buckles up once that the surface tension can no longer hold down all the cells in the monolayer. Such buckling occurs at certain radial distance from the center of colony; cf. Figure 8D. Bottom: Once the monolayer buckles, the colony starts to grow vertically. In the meantime, the buckling region moves outward as the colony expands radially. (C) The parameters used in the definition of surface tension force. The parameter hcell is the height of a cell that sticks out of the macroscopic colony surface; it is measured from the mean height of the agar surface z=0 to the 'tallest point' in the cell. The blue and red lines describe the macroscopic water level (indicated by hw) and colony height (indicated by h), respectively. The parameter δh is used to control how tightly the surface tension holds back those cells on the top of the colony. See more details in Appendix 1.4c.

https://doi.org/10.7554/eLife.41093.026

For E. coli growing on hard agar, the cells themselves have no appreciable attraction to one another. They are instead held together at high densities in a colony above the agar through the surrounding liquid they share (blue color in Figure 14A): Liquid is pulled into and retained in the colony through the osmotic effects of hydrophilic molecules on cell surface (Seminara et al., 2012), wetting the surface of cells in the colony including those at the boundary. One can think of the cellular density within the colony as determined by the osmotic balance between the colony and the agar. This can be readily observed, as colonies grown on lower density agar are more liquid-like.

In the same way, the cohesion of a three-dimensional colony is maintained by the interaction of cells with the surrounding liquid at the air-liquid boundary. As shown in Figure 14A, the red curve indicates a smooth air-liquid boundary preferred by minimization of the liquid surface tension. Wherever a cell protrudes sharply out of the smooth surface, it drags out the liquid surrounding the cell, resulting in increased liquid surface tension, and hence a restoring force Fsurf acting on that cell. A detailed treatment of these physical effects, requiring both the cell configuration and the air-liquid boundary, is computationally untenable. Here, we do not model the liquid explicitly, but retain its effect on cells at the colony boundary via the restoring force. In Appendix 1.4c, we describe a toy model calculation which yields a saturating restoring force,

(8) Fsurf,0=πγsurfw0 ,

whose magnitude is proportional to the width of the cell w0 rather than the (much smaller) macroscopic curvature of the colony boundary. As shown in 'Radial expansion – quantitative analysis' and illustrated in Figure 14B, this large cell-level surface tension is able to hold a large group of cells in a monolayer above the agar surface, until the pressure inside the expanding monolayer (due to friction against motion on agar surface) exceeds a critical level to overcome the liquid surface tension resisting vertical protrusion, resulting in the 'buckling' of the monolayer into multiple layers.

To implement the surface tension force at the single-cell level in our model, we first compute the coarse-grained colony height hx,y,t (red curves in Figure 14AC) from the cell configurations. Then we compute the height of the coated liquid (blue curve in Figure 14C) hwx,y,t by adding δh to the colony height. This thickness δh depends primarily on the agar hardness, being larger for softer agar where cells are less tightly bound by the liquid. For each cell whose maximum height hcell exceeds the liquid height hw, we impose a restoring force normal to the liquid surface; see Figure 14C. As the magnitude of the saturating restoring force Fsurf,0 (Equation 8) is independent of the height difference δzhcell-hw for δz>0, the restoring force can be mathematically written as Fsurf=Fsurf,0u(δz), where u is the Heaviside step function. To avoid numerical instability, we make a linear extrapolation between Fsurf=0 and Fsurf=Fsurf,0 over a narrow transition region 0<δz<w0/10, which is 1/10 of the cell width.

Pressure calculation

Once all the individual forces exerted on a cell i described above are calculated, the net force Finet and the corresponding torque Tinet are calculated. Moreover, the pressure on the cell i can be also calculated as

(9) Pi=Vi-1jFjirji

where Vi is the volume of cell i, the index j runs through all the different forces experienced by cell i, and rji are the corresponding displacement vectors from the points where the forces are exerted to the cell center.

Coarse-grained variables

We define coarse-grained fields of cell spatial mass density ρ(r,t), velocity v(r,t), directors n(r,t), and pressure P(r,t) by averaging the corresponding individual quantities over small regions in the colony (e.g., a few finite-difference grid boxes); see Appendix 1.5 for details.

Model parameters

We fix the coefficient of nutrient diffusion in the agar to be D-=600 μm2/s, which is typical for the diffusion of small molecules in solution (Beuling et al., 1996; Cole et al., 2015). The diffusion coefficient in the colony is much smaller due to the fact that bacterial cells are not permeable to most sugars. We take D+=90 μm2/s with the influence of volume fraction and tortuosity; see Supplementary file 1-Table S2 for details. We take the value of the yield factor for different sugars used to be that for glucose, which is Y=0.5 gCDW/gglucose (Payne, 1970). As shown above in 'Radial and vertical growth of the colony', the local cell mass density is found to vary only mildly around an average of 0.68 ρcell, with  ρcell=0.137×10-12 gCDW/μm3 being the cell dry weight density (Basan et al., 2015). The Monod constant is taken to be that for glucose, KS=20 μM (Monod, 1949). The cell dividing length ldiv=3 μm and the diameter w0=1 μm are fixed in all the simulations unless otherwise indicated. The constant value Cs of nutrient concentration in the boundary conditions for the diffusion equations and the batch culture growth rate λs are used as control parameters in our simulations. Other parameters that are crucial to the colony patterns and growth dynamics include various friction coefficients. See Supplementary file 1-Tables S2–S4 for the definitions and estimated values of all the parameters.

Numerical implementation

We use an iteration algorithm for our simulations. It has two loops. The main loop, the 'outer loop', consists of the following 3 steps: (i) define the colony region using the local spatial cell density ρ=ρ(r,t); (ii) update the nutrient concentration by solving the diffusion equations in steady state; and (iii) simulate cell growth, division, and movement over a small-time increment. The last step has its own loop, the 'inner loop', consisting of the following steps: update the local cell growth rate by Equation 3; simulate cell growth and division; and compute the forces and torques on cells, update the cell velocities and angular velocities, and update all the cell positions. We use the velocity Verlet algorithm, a commonly used molecular dynamics simulations of macromolecules, to update the cell velocities and positions (Frenkel and Smit, 2002). The inner loop is determined with a time step Δt. Usually, we run through one main loop per 100—1,000 inner loops. In updating the nutrient concentration, we use the finite difference to discretize the equations and the Jacobi or Gauss-Seidel relaxation method to solve the resulting systems of linear equations. We use multi-resolution adaptive grids for a large computational domain, and use the OpenMP for parallelizing our code. See Appendix 1 for details. On a multi-processor (14-16 processors) computer, the simulation can reach a colony of a few million cells in 24 hours. We have placed the major and basic parts of our C++ codes in the repository GitHub (Warren et al., 2019; copy archived at https://github.com/elifesciences-publications/CellsMD3D).

Appendix 1

Simulation Model and Methods

We model the colony expansion through the coupling of the growth, division, and movement of individual cells with the diffusion and reaction of nutrients and wastes. The growth of an individual cell within a short time period is described by a linear growth equation. The local growth rate varies spatially and temporally, and is determined locally by the cell density and nutrient supply. Once a cell grows into a critical size, it divides itself into two daughter cells with some randomness in their cylindrical lengths and orientations. In the meantime, growing and dividing cells push each other, generating mechanical forces. These forces, together with those arising from the cell-agar interactions, cell-liquid interactions, and surface tension, determine the movement of individual cells described by Newton’s law of motion. At any instant of time, the coarse-graining of all cells through their spatial positions determines the cellular colony region. Nutrients and wastes diffuse in both the agar and colony regions, while their reactions only occur in the colony region. In our current study, we only include one species of nutrient and we do not consider any wastes.

A1.1 Set Up and Main Algorithm

Our computational box is Ω=(-L,L)×(-L,L)×(-a,b), where all L, a, and b are positive numbers in the units of length; cf. Figure 11 in the main text. It is divided into the air region Ω0, colony region Ω1, and agar region Ω2=(-L,L)×(-L,L)×(-a,0), respectively. The colony surface or colony-air interface Γ01 separates the colony from air. The plane z=0 in the computational box is divided into two parts. One is the interface that separates the colony from agar, and is denoted by Γ12. The remaining part, denoted Γ02, separates the air from agar. Note that, since the bacterial colony grows with time t, all the air region Ω0, the colony region Ω1, the colony-air interface Γ01, and the colony-agar interface Γ02 depend on time t. All the simulations of cell growth, division, and movement, and the force calculations are done in the colony region which expands with time. The reaction-diffusion equation for nutrient is solved in both the colony region Ω1 and the agar region Ω2 (where there is no reaction). The growth rate is also defined everywhere in the colony region.

We cover the computational box with a finite difference grid with a grid size hgrid. We generate random and small heights at the grid points on the mean agar surface (z=0) to construct a rough agar surface. These random heights are in the range [0,hran] with hran an input number representing the possible maximum height. Such a rough surface will be used to calculate the interaction force between a cell and the agar. Initially, we set the nutrient concentration to be a nonzero constant in the agar region but zero elsewhere. We also randomly distribute a set of initial cells on the agar surface, and define their initial velocities and angular velocities to be all zero.

Our simulation continues through time iteration that consists of two loops. The main loop, or outer loop, consists of the following steps:

(1) Generate the boundary of colony;

(2) Update the nutrient concentration and cell growth rate;

(3) Simulate the cell growth, division, and movement.

The last step is carried out through an inner loop, a time iteration with time step Δt, that consists of the following steps:

(3.1) Simulate the cell growth and division;

(3.2) Compute the forces and torques on cells, and update the cell velocities and angular velocities with a half time step;

(3.3) Update all the cell positions;

(3.4) Compute again the cell forces and torques, and update the cell velocities and angular velocities with the other half time step;

(3.5) Set t:=t+Δt and continue with Step (3.1).

Note that Steps (3.2)–(3.4) are the velocity Verlet algorithm (cf. Frenkel and Smit, 2002) used for the simulation of cell movement. Practically, we update the nutrient concentration once every 1001000 time steps of calculations of cell growth, division, and movement.

A1.2 Nutrient and Growth Rate Update

Cell density and colony boundary

Given the positions of all the cells at time t, we define the (local) volume fraction ϕ=ϕ(r,t)(r=(x,y,z)) of the cells in each finite difference grid box above z=0 by

ϕ= sum of volumes of the cells inside the grid box volume of grid box.

A cell is inside the grid box if the center of this cell is in this box. The volume of a cell is given by the formula in Equation (A1.3.1) in section A1.3 below. By averaging over nearest grid boxes, we obtain the volume fraction of cells at each grid point. We then define the cell mass density (at grid points) to be

(A1.2.1) ρ(r,t)=ϕ(r,t)ρcell,

where ρcell is the constant mass density of a typical mature cell, and can be estimated from experiment; cf. Supplementary file 1-Table S3. We also define the colony region Ω1=Ω1(t) to be that with ϕ(r,t)>0; cf. Figure 11 in the main text.

Nutrient update: Reaction-diffusion equations and boundary conditions

The concentration field C=C(r,t) of the nutrient is defined spatially on the colony and agar regions Ω1=Ω1(t) and Ω2, respectively; cf. Figure 11 in the main text. It is governed by the following system of reaction-diffusion equations, interface conditions, and boundary conditions:

(A1.2.2) Ct=D+ΔCλSρYCC+KSinΩ1(t),
(A1.2.3) Ct=DΔCinΩ2,
(A1.2.4) C=C+onΓ12(t),
(A1.2.5) DCz=D+C+zonΓ12(t),
(A1.2.6) Cn=0onΓ01(t)Γ02(t)Γb,
(A1.2.7) C=CsonΓs.

The first two equations, Equation (A1.2.2) and Equation (A1.2.3), are the reaction-diffusion equation for the concentration in the colony region Ω1(t) and the diffusion equation for the concentration in the agar region Ω2, respectively, where D- and D+ are the corresponding diffusion coefficients, Y is the yield factor, λS is the batch culture growth rate, KS is the capacity constant (the Monod constant) for sugar, and ρ=ρ(r,t) is the local cell mass density (cf. Equation (A1.2.1)). As the density of cells is rather uniform (cf. discussions in Appendix A2.1), we use a constant density value ρ0 to approximate ρ=ρ(r,t). We take this constant ρ0 to be the cell dry weight (CDW) per unit volume of the colony; cf. Supplementary file 1-Table S3. Equation (A1.2.4) and Equation (A1.2.5) are the interface conditions for the concentration across the colony-agar interface Γ12(t) (i.e., z=0), where + and - denote the colony side and agar side, respectively. The last two equations, Equation (A1.2.6) and Equation (A1.2.7), are the boundary conditions. On the agar-air interface Γ01, the colony-air interface Γ02, and the bottom of agar Γb, we impose the flux-free boundary condition, with /n denoting the derivative in the normal direction along the corresponding part of the boundary, pointing from the colony or agar to the air region or pointing downward from the bottom of agar. On Γs, the lateral faces of the agar region, we prescribe a constant value Cs of nutrient concentration that represents the maximum nutrient that the system supplies.

In each time step, we update the nutrient by solving the steady-state reaction-diffusion equations with the corresponding boundary conditions, that is, Equation (A1.2.2)–(A1.2.7) with C/t set to be 0. We use an iterative scheme to solve the equations with the previous nutrient concentration as the initial solution. This iterative scheme is constructed based on solving the corresponding time-dependent equations with the fixed colony region, and the time here is only a numerical parameter. We use the forward Euler’s method to discretize this numerical time (cf. Gustafsson et al., 2013Morton and Mayers, 1995). The spatial derivatives of concentration are discretized with central differencing schemes. For a grid point that is near the colony-air interface but is outside the colony region, we assign a value of nutrient concentration by interpolating the values of such concentration at nearby grid points inside the colony. To treat a large computational region Ω2 that represents the agar region, we use a nested, multi-level, finite difference grid as shown in Appendix 1—figure 1. Techniques of interpolation are employed to discretize the diffusion equation on grid points at the interface of grids with different levels. In each numerical time step, we sweep the grid points from top down to those at the interface z=0, and further down to the bottom Γb, and then from grids on the bottom Γb up to those at z=0 and further up to the top. The numerical time iteration stops if the difference between the concentration fields of the current and previous numerical time steps is smaller than a given tolerance errconv; cf. Supplementary file 1-Table S5. To ensure the numerical stability, we chose a numerical time step that approximately satisfies the CFL condition (Gustafsson et al., 2013Morton and Mayers, 1995). We use OpenMP parallelization for both nutrient update and cell activities. Our simulations speed up more than 15 times in one 2.5 GHz Intel Xeon cluster node with 12 cores and 24 CPU processors.

Appendix 1—figure 1
Schematic view of a nested finite difference grid for the agar region Ω2.
https://doi.org/10.7554/eLife.41093.029
Growth rate update

Once the nutrient concentration field C(r,t) is updated, we can define the spatially and temporally varying local mass growth rate λ=λ(r,t) by

(A1.2.8) λ(r,t)=λSC(r,t)C(r,t)+KS.

This is first defined at each grid point r and then at any other point in the colony region by interpolation.

A1.3 Cell Growth, Division, and Movement

We model an underlying E. coli bacterial cell as a sphero-cylinder; cf. Figure 12A in the main text. We denote by p and q the centers of the two hemispheres. We call =p-q the cell cylindrical length which can vary with time t. We also denote by n=(q-p)/ the unit vector pointing from one center of hemisphere p to the other q, and call it the direction of the cell. Note that the center of mass of the cell is rc=(p+q)/2. We denote by w0 the diameter of each of the hemispheres. The volume and mass of the cell are given by

(A1.3.1) Vcell=14πw02+16πw03 andMcell=ρcellVcell,

respectively, where ρcell is the constant cell mass density introduced in Equation (A1.2.1). We shall assume that all cells have the same diameter w0 of hemispheres, independent of time.

Cell growth

With our assumption, a cell grows only in its cylindrical length =(t) but not in its diameter of hemispheres. The growth of the cell is governed by the growth equation for the cell cylindrical length

(A1.3.2) d(t)dt=λ~(rc,t)(t).

Here, λ~(rc,t) is the cell elongation rate evaluated at the cell center rc. The elongation rate λ~ is proportional to the mass growth rate λ defined in Equation (A1.2.8): λ~(r,t)=σλ(r,t). This effective relation between the two growth rates will be discussed after we describe the cell division. Numerically, we use the first-order approximation to obtain the cell cylindrical length at time t+Δt by

(t+Δt)=(t)+λ~(rc,t)Δt,

where Δt is the simulation time step. Initially at t=0, all the cells start with a constant cylindrical length 0.

Cell division

When a cell of centers of hemispheres p and q grows long enough, with its cylindrical length div for some critical value div>w0, it divides into two daughter cells of cylindrical lengths 1 and 2, respectively; cf. Figure 12B in the main text. The two centers of hemispheres of the mother cell become the centers of hemispheres of the daughter cells. The lengths 1 and 2 of these daughter cells are given by

1=12divw02+ηranand2=12divw02ηran,

where η(-0.5,0.5) is a random variable and ran is an input positive number representing the maximum fluctuation of cell cylindrical length during the cell division; cf. Supplementary file 1-Table S3. In average, all new born cells have the same cylindrical length 0; and we take it to be the same for all the initial cells. This implies that div=20+w0; cf. Figure 12B in the main text. Note that before division, the mother cell has volume πw02div/4+πw03/6; after division, the total volume of two daughter cells is πw02div/4+πw03/12. There is πw03/12 volume loss every time due to division. For a fixed cell aspect ratio div:w0=3:1, which is what we have in our simulations, the volume loss is about 9%. Note also that we have used a constant dividing length div independent of the batch culture growth rate λS. The effect of this simplification is discussed in Appendix A2.2.

The centers of hemispheres of the two daughter cells are given by (cf. Figure 12B in the main text)

p1=p,q1=p-1(p-q),p2=q+2(p-q),q2=q.

Moreover, these daughter cells inherit the velocity from their mother cell. But the angular velocities of these two daughter cells are set to be

ω1=(0,0,0)andω2=ωran(0,0,ξ),

respectively, where ξ(-0.5,0.5) is a random variable and ωran is an input positive number that is the maximum fluctuation of the angular velocity; cf. Supplementary file 1-Table S3.

In the cell division, each of the two daughter cells also experiences the rotational fluctuation. Let us fix the center of mass rc of a daughter cell. We rotate the vector p-rc with p the center of a hemisphere of this daughter cell. (For notational simplicity, we use rc and p for this daughter cell; and they are different from those of the mother cell.) To do so, we construct a local Cartesian coordinate system with the origin at rc and the three unit coordinate vectors e1, e2, and e3 by e3=(prc)/||prc,

e1=(p-rc)×argmini{ex,ey,ez}(p-rc)i(p-rc)×argmini{ex,ey,ez}(p-rc)i,

and e2=e3×e1, where ex,ey,ez are the unit coordinate vectors in the original coordinate system. Then we rotate the vector p-rc by an angle of φ around the axis e2, and θ around e3 successively. Here, φ is a random variable in the range [0,φranπ], where φran is a constant that sets the magnitude of fluctuations of φ (cf. Supplementary file 1-Table S3) and θ[0,2π) is a random number. If we denote by pnew the new center of hemisphere after the rotation, then

pnew-rc=p-rc(sinφcosθe1+sinφsinθe2+cosφe3).

This equation allows us to find the coordinates of pnew in the original coordinate system. The coordinates of the other center of hemisphere qnew are given by qnew=2rc-pnew.

The two growth rates

We assume that, in the life span of a cell, the cell mass growth rate λ and the cell elongation rate λ~ stay constant. The growth of mass satisfies the equation M(t)=M(0)eλt. If the mass doubling time is td, then M(td)=2M(0) and hence λtd=ln2. On the other hand, the doubling time is the time that a new born cell, which in average has the cylindrical length 0, grows as its cylindrical length reaches the dividing length div. Note that div=20+w0. From the growth equation Equation (A1.3.2), we have then div=l0eλ~td. Therefore,

λ~td=lndiv0=ln2divdiv-w0.

Finally, λ~=σλ, where

σ=λ~λ=ln(div/0)ln2=ln(2div/(div-w0))ln2.

For the fixed dividing cell aspect ratio div:w0=3:1, we have σ=ln3:ln2.

Cell movement

Consider a cell at some time instant. Let us denote its centers of hemispheres by pold and qold, its cylindrical length by old=pold-qold, its direction by nold=(qold-pold)/old, and its mass by Mold. Let us also denote by vold and ωold the velocity and angular velocity, respectively, at the center of the cell. We apply the velocity-Verlet algorithm to update the cell position, velocity, and angular velocity with the simulation time step Δt.

We first calculate the force Fhalf and torque Thalf of the cell. Details of such calculations are given below in section A1.4. We then calculate the velocity vhalf and angular velocity ωhalf for a half time step:

vhalf=vold+Δt2FhalfMold,
Thalf,n=(Thalfnold)nold,
Thalf,t=Thalf-Thalf,n,
ωhalf=ωold+Δt2(Thalf,nIold,n+Thalf,tIold,t),

where Iold,n and Iold,t are the moments of inertia of the cell along the directions nold and its orthogonal, respectively. By direct calculations and using the constant density ρcell in the place of the mass density of an underlying cell, we have

Iold,n=ρcellπw04(132old+160w0),
Iold,t=1480ρcellπw02(4w03+15oldw02+20old2w0+10old3).

We now update the cell positions for the entire time step Δt by updating the centers of hemispheres

pnew=pold+Δt(vhalf+ωhalf×pold-qold2),
qnew=qold+Δt(vhalf+ωhalf×qold-pold2).

We update the force and torque of the new cell with the centers of hemispheres pnew and qnew by the procedure of force calculations described below in section A1.4 to get Fnew and Tnew. Finally, we update the velocity and angular velocity for the second half time step to get

vnew=vhalf+Δt2FnewMnew,
Tnew,n=(Tnewnnew)nnew,
Tnew,t=Tnew-Tnew,n,
ωnew=ωhalf+Δt2(Tnew,nInew,n+Tnew,tInew,t),

where Mnew, nnew=(qnew-pnew)/new, and new=pnew-qnew are the mass, direction, and cylindrical length, respectively, of the new cell with the centers of hemispheres pnew and qnew, and the moments of inertia Inew,n and Inew,t can be calculated similarly using the cylindrical length new.

A1.4 Force Calculations

Mechanical forces exerted on a cell include the elastic and dissipative forces from the cell-cell mechanical interactions, the elastic and dissipative forces from the cell-agar interaction if the cell is in contact with the agar, the surface tension force if the cell is on the top of the colony, and the viscous force from the interaction between the cell and the surrounding liquid. For a given cell, we shall denote these forces by Fcc, Fca, Fsurf, and Fvisc, respectively. So, the total force acting on the cell is

F=Fcc+Fca+Fsurf+Fvisc.

Note that most cells are in the interior of the colony; and they only experience the forces from the cell-cell and cell-liquid interactions.

(a) Cell-cell interaction force

When two cells are in direct contact, they generate the cell-cell interaction force. As in Volfson et al. (2008), we describe such forces using a standard model in granular solids (Brilliantov et al., 1996; Cundall and Strack, 1979; Kuwabara and Kono, 1987; Makse et al., 2004; Shäfer et al., 1996), with some adjustment based on experimental considerations on bacterial cells. An important part of our force scheme is a detailed treatment of the cell-cell frictional force, which together with the cell-agar friction are responsible for crucial mechanical behaviors such as buckling of the bacterial colony. While a similar grain-grain friction is commonly included in models of granular solids, however, friction in cellular models is often described to only include a purely viscous force (Farrell et al., 2013; Ghosh et al., 2015).

Let us fix a cell in the colony centered at rc, and call it a primary cell. Let us also fix a neighboring cell centered at rc. We denote by d the minimal distance between the central cylindrical line segments of the two cells, and by a and a the points on these line segments, respectively, that reach this minimal distance:

d=a-a;

cf.Figure 13A in the main text. We will describe an algorithm of finding the minimum distance d and these two points a and a at the end of this part. We denote by ncc the unit vector along the direction from a to a, by rcc the center of these two points, and by δcc the indentation size (i.e., the amount of overlap) between the two cells, respectively:

ncc=1d(aa);rcc=12(a+a);δcc={w0dif w0d>0,0if  w0d0.

Let v, v and ω, ω be the velocities and angular velocities of the primary and neighboring cells, respectively. The velocities V and V of the two cells at the midpoint rcc are then given by

V=v+ω×(rccrc)andV=v+ω×(rccrc),

respectively. Denote the difference of these velocities by

vcc=V-V.

We call the unit vector ncc the normal direction in the interaction of these two cells. We specify the tangential direction to be the unit vector

τcc=vcc-(vccncc)nccvcc-(vccncc)ncc,

if the denominator is nonzero, that is the relative velocity vcc has a nonzero tangential component. (Otherwise, τcc can be any unit vector orthogonal to ncc.) Note that the tangential direction τcc depends on the direction of the relative velocity vcc.

Let us assume now that these two cells are in direct contact with each other, that is δcc>0. The total interaction force Fcc between these two cells exerted at the center rcc and the corresponding torque Tcc about the axis rcc-rc, are given, respectively, by

Fcc=Fcc,n+Fcc,t,
Tcc=(rcc-rc)×Fcc.

Here, Fcc,n is the normal force in the direction ncc and Fcc,t the tangential force in the direction τcc. They are defined by

(A1.4.1) Fcc,n=[23kccw0δcc3/2-γcc,nMeffδcc(vccncc)]ncc,
(A1.4.2) Fcc,t=-min{γcc,tMeffδcc1/2|vccτcc|,23μcckccw0δcc3/2}τcc.

Note that the final form of the force Fcc is the same as that in the main text; cf. Equations 7a and b there.

The first part in the normal force Fcc,n is the Hertz contact force resulting from the elastic collision of the two cells and pointing from the neighboring cell to the primary cell. Here, we approximate the cells as spheres of radius w0/2. In such a case, the classical Hertz contact force is (4/3)kccw0/4δcc3/2, which is exactly what we have, where kcc is the reduced (or effective) elastic constant and w0/2 is the reduced radius. If E and ν are Young’s module and Poisson’s ratio of cells, then kcc=E/(2(1ν2)) (cf. Johnson, 1985Landau and Lifshitz, 1986Popov, 2010). The second part in the normal force Fcc,n is the friction force in the normal direction ncc due to the relative motion of the two cells, where γcc,n is the (normal) static friction coefficient and

Meff=ρcellVcellVcellVcell+Vcell

is the reduced mass.

This form of the normal force Equation (A1.4.1) has been used in Volfson et al. (2008). Similar forms of such a normal force can be found in the literature of granular solids (cf. Brilliantov et al., 1996; Herrmann and Luding, 1998; Kuwabara and Kono, 1987; Makse et al., 2004; Shäfer et al., 1996; Silbert et al., 2001). In particular, Kuwabara and Kono (1987) and Brilliantov et al. (1996) derived the normal force between two spherical granular grains, assuming such grains are viscoelastic. The frictional part in their derived normal force is proportional to Rredδv, where Rred is the reduced radius and equals R/2 if both spheres have the same radius R,δ is the amount of overlap of the grains (same as our δcc), and v is the relative speed of the head-on collision of the grains (same as our |vccncc|). Since we have always used a fixed radius of the spherical caps of a cell in our simulations, with or without the factor Rred does not change noticeably our simulation results. However, unlike in granular material simulations where the grain size is fixed, a cell in the bacterial colony can have different mass and size at different stages. So, including the dependence on the effective mass is reasonable.

The magnitude of the tangential force Fcc,t in Equation (A1.4.2), where γcc,t and μcc are the (tangential) static and dynamic friction constants, respectively, is the minimum of a tangential friction force that depends on the relative velocity vcc and the dynamic friction force that is proportional to the Hertz contact force which is part of the normal force Fcc,n in Equation (A1.4.1). This is different from that used in Volfson et al. (2008): we only include the Hertz contact force here, not the entire normal force, as the elastic collision force can be dominated in cell-cell interactions.

To understand the form of the tangential component of the friction force, let us first consider the friction for two solid objects in contact. Coulomb’s law for such friction states that the static tangential friction force Fstatic required to move the objects relative to each other and the dynamic tangental friction force Fdynamic between the objects to maintain such motion once it is initiated are both proportional to the normal force Fnormal pressing the objects together (Popov, 2010):

Fstatic=μsFnormalandFdynamics=μdFnormal,

where μs and μd are the static and dynamics friction coefficients, respectively. In general, μd<μs and μdμs; cf. Appendix 1—figure 2A. Here we assume for simplicity that μd=μs. In addition, we assume that the friction force is proportional to the relative velocity of the two cells in the tangential direction, with the proportion constant γ, when the tangential relative velocity is small; cf. Appendix 1—figure 2B. This allows the inclusion in the tangential friction of the dependence of the tangential velocity when it is small in magnitude and can be damped away quickly in the dynamics. Such an assumption is supported by our experimental observation that the radial velocity of colony expansion is linear in the batch culture growth rate but is independent of the individual cell velocity, as the ‘buckling length’ is independent of such local cell velocity. Our specific form of the tangential friction for small value of relative tangential velocity, that is the first part in the minimum of the tangential force Fcc,t in Equation (A1.4.2), is the same as that in Volfson et al. (2008). It is different from a form, popularly used in simulations of granular solids, that only includes the linear dependence on the relative tangential velocity but not the effective mass and the overlap distance δcc (cf. Shäfer et al., 1996) and the references therein. Moreover, for simplicity of simulations, we have also neglected the history dependence in the tangential velocity (Cundall and Strack, 1979; Shäfer et al., 1996).

Appendix 1—figure 2
The standard (A) and modified (B) model for friction between two objects, as a function of tangential relative velocity between the two objects, where μs and μd are the static and dynamic friction constants, respectively.
https://doi.org/10.7554/eLife.41093.030

We end this part with a method of computing the minimum distance d between the central cylindrical line segments γ1 and γ2 of the two cells and the corresponding points on these segments that reach this minimum distance. For notational convenience, let us denote by pi and qi the position vectors of the centers of hemispherical caps of the cell i with i=1 and 2, respectively. Let ui=qi-pi(i=0,1). We parameterize the central cylindrical line segments γi by ri(t)=pi+tui for 0t1(i=1,2). Denote p0=p1-p2 and define the distance-square function

f(s,t)=r1(s)-r2(t)2=su1-tu2+p02.

Clearly, f is minimized in [0,1]×[0,1] by some point (s0,t0)[0,1]×[0,1]. The minimum distance d>0, and the two points a1 and a2 on the two line segments reaching this distance are then given by

(A1.4.3) d2=f(s0,t0)=mins,t[0,1]f(s,t)anda1=r1(s0),a2=r2(t0).

We first assume that the two lines are not parallel, that is u1×u20. In this case, f(s,t) is a strictly convex and quadratic function with a constant symmetric positive definite Hessian matrix. By setting sf(s,t)=0 and tf(s,t)=0, we find that f is minimized in 2 at (s0,t0) with

s0=(p0×u2)(u1×u2)u1×u22andt0=(p0×u1)(u1×u2)u1×u22.

If (s0,t0)[0,1]×[0,1], then we obtain d and a1,a2 by Equation (A1.4.3). Otherwise, we compute the minimum value of f on each of the four sides of the square [0,1]×[0,1] and compare these values to find the global minimum points (s0,t0)[0,1]×[0,1] and the corresponding minimum value of f on this square. Consider, for instance, the side s=0 and 0t1. The function f(0,t) for all t is found to be minimized at t1=(a0u2)/u22. If t1[0,1], then the minimum value of f on this side is f(0,t1). Otherwise, this value will be the minimum of f(0,0) and f(0,1).

We now assume that the two line segments γ1 and γ2 are parallel. In this case, the minimum distance d is achieved by the minimum of the distance from p1 to the second line segment γ2 and that from q1 to γ2:

d=min{dist(p1,γ2),dist(q1,γ2)}.

Each distance from a point to a line segment can be found by minimizing a convex quadratic function on [0,1], similar to and simpler than the previous case.

(b) Cell-agar interaction force

When a cell is in contact with the agar surface (including the case of the cell overlaps partially with the agar region), a cell-agar interaction force Fca is generated and can be modeled similar to the cell-cell interaction force. Assume one end of the cell dips into the agar region; cf. Figure 13B in the main text. Let p=(xa,ya,za) be the center of the spherical cap corresponding to that end of the cell. We denote by δca the indentation depth: δca=w0/2-za. We also denote rca=(xa,yz,za-w0/2), which is the midpoint of the line segment along the vertical line passing through the point p between z=0 and z=-δca. As before, we denote by vca the relative velocity at the center rca. It is given by

vca=v+ω×(rca-rc),

where v is the velocity of the cell at its center rc and ω is the angular velocity of the cell. The normal direction is now the positive z-direction; we denote the unit vector along this direction by nca=(0,0,1). The tangential direction is defined through the relative velocity vca by

τca=vca-(vcanca)ncavca-(vcanca)nca,

if the denominator is nonzero. (Otherwise, τca can be any unit vector orthogonal to nca.)

Similar to the cell-cell interaction, the total cell-agar interaction force exerted on the midpoint rca and the corresponding torque are given by

Fca=Fca,n+Fca,t,
Tca=(rca-rc)×Fca.

Here Fca,n and Fca,t are the forces normal and tangential to the mean colony-agar interface z=0. They are given by

(A1.4.4) Fca,n=[223kcaw0δca3/2-γca,nMeffδca(vcanca)]nca,
(A1.4.5) Fca,t=-min{γca,tMcellδca1/2|vcaτca|,223μcakcaw0δca3/2}τca.

where kca is an elastic constant in the Hertzian stress, γca,n and γca,t are static friction constants, μca is the dynamic friction constant, and Mcell=ρcellVcell is the cell mass with Vcell the cell volume. Note that the factor 2 in the Hertz contact force part is different from those for the cell-cell interaction case. Here, the mean agar surface can be treated as a sphere of radius which leads to the reduced radius of the cell-agar system to be w0/2.

In our numerical implementation, we do not decide which end of the cell dips into the agar region. Instead, we compute the corresponding forces at both ends and add them together.

(c) Surface tension

Bacterial cells are hydrophilic. They are coated with water molecules. Once a bacterial cell sticks out of the colony surface, the tension between the air and water surface generates the surface tension force that brings down the cell. Such surface tension has long been recognized as a critical component in colony growth. Existing models of surface tension for colony, however, rarely treat the cells and the surrounding liquid as distinct media. Rather, the surface tension is frequently treated as a property of a composite fluid of cells plus liquid (Grimson and Barker, 1993; Zhang et al., 2008). Alternatively, the liquid phase is ignored and surface tension is assumed to arise from attractive interactions between the cells themselves (Farrell et al., 2013). In both cases, the surface tension scales with the macroscopic curvature of the colony. Here, we endeavor to model the surface tension force as a boundary force between the discrete cells and the continuum liquid.

The surface tension force can be calculated using the virtual work principle, dW=γsurfdA, where W is the work done by the water surface, γsurf is the water-vapor surface tension, and dA is the change in water area A as the cell is raised with an additional dh; cf. Appendix 1—figure 3. Here, we approximate a cell by a sphere of diameter w0. (The notation for the radius is R here.) The surface area is A=2πRh and hence dA=2πRdh. (In the case of a disk in the two-dimensional setting, the change in area is dA=2πrRdθ=2πR2sinθdθ. Note that h=Rsinθ and dh=Rsinθdθ. Hence dA=2πRdh.) As a result, the surface tension force is

F=dWdh=2πγsurfR=πw0γsurf.
Appendix 1—figure 3
Schematic of the derivation of the surface tension force.
https://doi.org/10.7554/eLife.41093.031

Therefore, the surface tension is a constant force when a cell sticks out of the water. However, when a cell is below the water level, it should not experience any surface tension force. As a result, the surface tension force of a cell is discontinuous across the water surface. To reconcile this discontinuity, we make an approximation that the surface tension force increases linearly with the height until the height reaches some critical value hc0.1w0, when it saturates to its maximum value of πw0γsurf.

Specifically, let us consider a cell on top of the colony. Let rc denote the center, and p and q the two centers of hemispheres of such a cell; cf. Figure 14C in the main text. We define the total surface tension force Fsurf and torque Tsurf on this cell to be

Fsurf=Fsurf,p+Fsurf,q,
Tsurf=(p-rc)×Fsurf,p+(q-rc)×Fsurf,q,

where

Fsurf,i=πw0γsurfmin{max{0,(hi-zi)nsnz+δh}0.1w0,1}ns,i=p or q.

Here, hi is the water level at i(i=p or q), nz=(0,0,1) is the unit vector along the z direction, zi=inz is the z-component of i, ns is the unit vector of the colony surface, and δh is a constant determining how tightly the surface tension holds the cells. The water level is coarse-grained. First, on each horizontal grid square i,j=(iΔx,(i+1)Δx)×(jΔy,(j+1)Δy), we set the water level, hi,j, at the center of this grid square to be the maximum of the z-coordinate of the two centers of hemispherical caps of cells whose centers are in the grid column i,j×[0,b]. (Recall that z=b is the top of our computational box; cf. Figure 11 in the main text.) Then we construct the coarse-grained water level everywhere by a continuous and piecewise linear interpolation at all hi,j.

(d) Viscous force

The cells in the colony experience a drag force—the Stokes drag force—due to their interactions with the surrounding liquid. Such viscous force Fvisc exerted at a cell by the surrounding liquid and the corresponding torque Tvisc are given by (note that w0 is the diameter)

Fvisc=3πμliqw0vandTvisc=πμliqw03ω,

respectively, where v and ω are the velocity and angular velocity, respectively, at the center of mass of the cell, and μliq is the liquid viscosity.

We finish our description of forces with a remark on the static and dynamic frictions. A static friction is proportional to the cell speed, while a dynamic friction is the same as the static friction for small speed but saturates after the speed is large; cf. Appendix 1—figure 2. In our tangential friction forces that arise from the cell-cell and cell-agar interactions, the saturation is controlled by capping through the elastic force; cf. Equation (A1.4.2) and Equation (A1.4.5). To compare our dynamic friction forces with static friction forces alone, we shall consider a static friction model where the tangential friction forces in Equation (A1.4.2) and Equation (A1.4.5) are replaced by the following static frictions:

(A1.4.6) Fcc,t=-γcc,tMeffδcc1/2|vccτcc|τccandFca,t=-γca,tMcellδca1/2|vcaτca|τca,

respectively. The difference between the static and dynamic friction models is shown in Figure 8F, and the related discussions are given in Discussion in the main text.

A1.5 Coarse-Grained Variables

To better present our simulation results, we need to coarse-grain the cell volume fraction ϕ, pressure field P, velocity field v, and director field n over a subregion 𝒢 of the colony region. Examples of a subregion 𝒢 include:

  • The union of a few grid boxes for coarse-graining in the entire colony;

  • A small box at the agar-colony interface (iΔx,(i+1)Δx)×(jΔy,(j+1)Δy)×(0,Δz) for some i, j for coarse-graining around such interface;

  • A small box in a vertical layer (-Δx,Δx)×(jΔy,(j+1)Δy)×(kΔz,(k+1)Δz) or (iΔx,(i+1)Δx)×(-Δy,Δy)×(kΔz,(k+1)Δz) for some i, j, and k for coarse-graining around the cross-section x=0 or y=0, respectively; and

  • A cylindrical ‘cube’ (iδr,(i+1)δr)×(jδθ,(j+1)δθ)×(kΔz,(k+1)Δz) in the cylindrical coordinates (r,θ,z) for some small δr>0, δθ>0, and Δz>0, and for some i, j, and k, for azimuthal coarse-graining.

Now let us fix a subregion 𝒢 in the colony. Let f denote the pressure P or a component of the velocity field v, and denote by fi such a quantity at the center of cell i. We define the coarse-grained average of f over 𝒢 to be

f¯(𝒢)=rci𝒢firci𝒢1,

where rci is the center of mass for the cell i. Similarly, we define the coarse-grained volume fraction over the subregion 𝒢 to be

ϕ¯(𝒢)=rci𝒢ViVol(𝒢),

where Vi is the volume of cell i and Vol(𝒢) is the volume of the subregion 𝒢. Note that, if 𝒢 is a grid box, then ϕ¯(𝒢) is the same as the volume fraction ϕ on that box; cf. section A1.2. Note also that the cell density ρ can be coarse-grained following its relation with the volume fraction ϕ; cf. Equation (A1.2.1).

The director field n cannot be coarse-grained simply by summing over the directors in a subregion, since the director of a given cell can have two possible directions and the sum can lead to artificial cancellations. Here, we define the coarse-grained director field n(𝒢) over a given subregion 𝒢 to be a maximizer of the maximization problem

maxnrci𝒢|nni|2subjectton=1.

By the Lagrange multiplier method, this leads to an eigenvalue problem for a 3-by-3 matrix, with the maximizer n(𝒢) being a unit eigenvector that corresponds to the largest eigenvalue.

Let δr>0 be small. Let Nθ2 be an integer and define δθ=2π/Nθ. We denote 𝒢i,j,k=(iδr,(i+1)δr)×(jδθ,(j+1)δθ)×(kΔz,(k+1)Δz) in the cylindrical coordinates (r,θ,z), We define the azimuthal average of a scalar field f by

f~(iδr,kΔz)=1Nθj=0Nθ1f¯(𝒢i,j,k),

where f¯(𝒢i,j,k) is the the coarse-grained average of f over 𝒢i,j,k. The azimuthal average of a vector field can be defined componentwise.

Given any point in the agar-colony interface with the polar coordinates (r,θ), we define

(A1.5.1) Δr=rR(θ),where R(θ)=maxρ(r,θ)>0r,

where ρ(r,θ) is the local cell density projected onto the colony-agar interface. This is the negative distance between this point and the colony edge along the ray of angle θ. For each integer j with 0jNθ, we denote by ij the largest integer not exceeding (R(jδθ)+Δr)/δr. We then define the azimuthal average of the radial component Vr of the velocity field v by

(A1.5.2) Vr(Δr)=1Nθj=0Nθ1v¯(𝒢ij,j,0)rj,

where rj=(cos(jδθ),sin(jδθ)). Other azimuthal-averaged quantities in terms of Δr can also be defined similarly.

https://doi.org/10.7554/eLife.41093.028

Appendix 2

Additional Results of Simulation and Analysis

Unless otherwise stated, all the notations and terms are the same as defined in the main text and Appendix 1.

A2.1 Constant Density and Volume Fraction

For each of the batch culture growth rates λS=1h-1 and 0.5h-1, we simulated the growth of bacterial colony. Figure 3—figure supplement 1 in the main text shows our simulation results on the volume fractions and their course-grained values. We observe that the volume fraction of cells in the colony has the mean value around 0.68 with the standard deviation around 0.03 for both of the batch culture growth rates λS=1h-1 and 0.5h-1. This suggests that we can treat the cell volume fraction as a constant inside the colony. Hence, we can also approximate well the density ρ by a constant ρ0; cf. Equation (A1.2.1) in Appendix 1.

A2.2 Effect of λS-Dependence on Cell Dividing Length

It has been known that cell dividing length div may vary with the batch culture growth rate λS according to the relation (cf. Donachie, 1968; Jun and Taheri-Araghi, 2015Wallden et al., 2016)

(A2.2.1) Vdiv=V0eλS/λS0,

where Vdiv is the corresponding dividing volume of the cell, V0 is a constant, and λS0=1.0h-1. In our simplified computational model, we have not included such variations. Here, we study the effect of the λS-dependence on the cell dividing length div by computer simulations.

Based on recent experimental observations on the size on E. coli cell (Jun and Taheri-Araghi, 2015; Si et al., 2017), we have the ratio

(A2.2.2) div:w0=3:1

for all the batch culture growth rate λS. This and Equation (A2.2.1), together with the formula for cell volume (cf. Equation (A1.3.1) in Appendix 1), then imply that the cell dividing length should be given by

11324πdiv3=V0eλS/λS0.

Setting λS=λS0=1h-1 in Equation (A2.2.1) and using the assumption Equation (A2.2.2), we obtain also that

11324πdiv,03=V0e,

where div,0 is the cell dividing length for λS=λS0. We now assume that this cylindrical length is div,0=3μm (Jun and Taheri-Araghi, 2015; Si et al., 2017). The above two equations then provide us with the cell dividing length

(A2.2.3) div=27eλS/λS0-1.

We chose λS=0.1h-1,0.2h-1,,1.0h-1. For each of these values of λS, we ran simulations with the variable dividing length div determined by Equation (A2.2.3). We also ran simulations with a fixed diving length div=3μm; cf. Supplementary file 1-Table S3. In Figure 5C in the main text, we plot the (constant) vertical ascending speed VH vs. λS for both fixed (open circles) and variable (filled circles) dividing lengths. We observe that the results from a fixed dividing length are consistent with those from a variable dividing length. In Figure 8F in the main text, we see that the (constant) radial velocity VR obtained with a fixed dividing length is close to that with a variable dividing length, but the discrepancy is more significant than the case of VH. In Appendix 2—figure 1 below, we plot the velocities VR and VH for variable dividing length, and fit the simulation data with λS0.5h-1. We observe that the straight line that fits simulated VR intersects the x-axis at 0.2h-1. This agrees well with the experimental result plotted in Figure 1H in the main text. Hence, the inclusion of λS-dependence on the cell dividing length can better describe experiment.

Appendix 2—figure 1
Experimental measurement on VR and VH with various batch culture growth rates.

Data are fitted with the stright line VR=51.27λS7.96 and the squre root curve VH=12.7λS for the redial and vertical speeds of expansion VR and VH respectively. (B) Simulation results on VR and VH with various batch culture growth rates and with a variable diviting length ldiv. data are fitted for λS0.5h1 using the straight line VR=22.1λS5.24 the squre-root curve VH=6.12λS, respectively.

https://doi.org/10.7554/eLife.41093.033

We now show that using a fixed cell dividing length div independent of the batch culture growth rate λS is a reasonable simplification for our simulations. We fixed the batch culture growth rate λS=1.0h-1 and the constant concentration in the boundary condition Cs=2.0mM. We then distributed randomly 625 cells at time t=0h. At t=9.0h, there are around 0.16 million cells in the colony. We picked randomly 2,000 of them from the bottom layer, and then tracked the local mass growth rate of each of these cells from t=9.0h to t=15h. If a cell divides during this time period, we kept track one of its two daughter cells. We found that the local growth rates of about 80% of these cells change from high values, larger than 90% of λS, to low values, smaller than 10% of λS, during this time period of colony growth. This indicates that the majority of the cells go through a complete transition in local growth rates.

We now consider those 80% cells that experience the transition from high to low growth rates. Appendix 2—figure 2A is the histogram of the number of cell generations (i.e., the number of divisions) of these cells. It is clear that during the transitioning time period, most of these cells did not divide or only divided once, signaling the sharpness of the high-to-low growth rate transition. To better understand such sharpness, we selected randomly 100 cells which complete the high-to-low transition, and tracked the growth rate of each of these cells during the time period from t=9.0h to t=16h. In Appendix 2—figure 2B, we plot the local growth rate vs. shifted time for each of these 100 cells. For a given cell, we shifted the time so that the growth rate of λS/2 was reached at the shifted time 0. All these indicate that cells pass the transitioning region in a short time, and it is therefore reasonable to assume a fixed ldiv and w0 for the entire simulation.

Appendix 2—figure 2
Statistics of cells undergoing growth transition.

(A) The histogram of the number of generations of the cells that experienced high-to-low growth rate transition.(B) Curves of growth rates vs. shifted time for 100 cells randomly selected from those 80% cells.

https://doi.org/10.7554/eLife.41093.034

A2.3 A One-Dimensional Model for Nutrient Penetration

The vertical ascending velocity and the pattern formation of the colony depend largely on how deep the nutrient from agar can penetrate into the interior of colony. Near the agar surface, cells grow fast with abundant nutrient supply. Away from the agar surface, the growth of cells in the colony is much limited due to the lack of nutrient. In this section, we construct a simplified 1D model of nutrient diffusion and reaction, and show that the nutrient concentration decays quadratically in the region where it is above the Monod constant KS but decays exponentially in the region where it is below KS. We introduce the nutrient penetration level HS to be the height (i.e., the z coordinate) at which the nutrient concentration takes the value KS, and analyze how HS affects the colony growth.

Since the colony is rather thin, the variations of the nutrient concentration in the x and y directions are rather small compared with that in the z direction. Therefore, we can assume that the steady-state concentration C=C(z) depends only on the z variable, and approximate the full Laplacian ΔC in the colony region by the 1D Laplacian zzC in the z variable. Consequently, we consider the following 1D model for the nutrient concentration C=C(z):

D+C=ρ0λSYCC+KSfor z>0,
C(0)=C0andC()=0,

where a prime denotes the spatial derivative and C0>0 is a controlling parameter that describes the amount of nutrient concentration from the agar substrate. Note that we have replaced the local cell density ρ by the constant density ρ0; cf. the description of nutrient update in Appendix A1.2. Setting

C~=CKSandz~=zβwithβ=(D+KSYρ0λS)1/2,

we obtain the non-dimensionalized form of the model

(A2.3.1) C~=C~C~+1for z~>0,
(A2.3.2) C~(0)=C~0andC~()=0,

where C~0=C0/KS.

We observe that there is a unique solution C~=C~(z~) to this boundary-value problem that is nonnegative for all z~0. This solution is the unique minimizer of the convex functional

I[u]=0[12u2+u-ln(1+u)]𝑑z~

of all nonnegative functions u such that both u and u are square-integrable on (0,), and u(0)=C~(0) and u()=0. Note that the solution C~=C~(z~) is a monotonically decreasing function of z~>0. For otherwise, C~ would reach a local maximum at some z~m>0 with C~(z~m)<0 but C~(z~m)0. This is impossible by Equation (A2.3.1). We also observe that C~()=0, for otherwise C~ would be negative and stay strictly away from 0 as C′′~0, leading eventually to C~()=- which would contradict the fact that C~()=0.

Now, multiplying both sides of Equation (A2.3.1) by C~, we get

12ddz~(C~)2=ddz~[C~-ln(C~+1)].

Integrating both sides of this equation from 0 to z~>0, we obtain

12(C~(z~))2-12(C~(0))2=[C~(z~)-ln(C~(z~)+1)]-[C~0-ln(C~0+1)].

Sending z~, we get

-12(C~(0))2=-[C~0-ln(C~0+1)].

The combination of the above two equations leads to

(A2.3.3) 12(C~(z~))2=C~(z~)ln(C~(z~)+1)z~>0.

We now study how fast C~=C~(z~) decays. Assume C~(z1)1 for some z~10. Then C~(z)1 for all z~z~1. Noting that

xln(x+1)x2/3if 0x1,

and that C~0, we have by Equation (A2.3.3) that -C~2/3C~ for all z~z~1. This leads to the exponential decay

(A2.3.4) C~(z~)C~(z1)e2/3(z~z~1)z~z~1.

If C~0=C~(0)1, then we can have z~1=0 in Equation (A2.3.4) and the concentration C~=C~(z~) decays exponentially.

Suppose now C~0>1. Since C~=C~(z~) decreases monotonically and C~()=0, there exists a unique z~0=z0/β>0 such that C~(z~0)=1 (i.e., C(z0/β)=KS) and C~(z~)1 for all z~z~0. Thus the inequality (Equation (A2.3.4)) holds with z~1=z~0 and C~(z~1)=C~(z~0)=1. We show now that the concentration C~=C~(z~) decays quadratically in [0,z~0]. Precisely, we shall prove that

(A2.3.5) [max(C~0-12z~,0)]2C~(z~)(C~0-ln(e/2)2z~)2 z~[0,z~0].

By Equation (A2.3.3) and the fact that C~0, we have -C~2C~ in [0,z~0]. This leads to (C~)-1/2. Integrating both sides of this inequality from 0 to z~, we obtain the first inequality in Equation (A2.3.5). Note that

(ln2)x-ln(x+1)0 x1.

Thus, by Equation (A2.3.3) and the fact that C~0, we obtain that -C~2ln(e/2)C~, and further that (C~)-(ln(e/2))/2 in [0,z~0]. An integration of both sides of this inequality from 0 to z~ then leads to the second inequality in Equation (A2.3.5).

In Figure 5—figure supplement 1 in the main text, we see that the nutrient concentration Cctr along the z-axis is described well by a quadratic function for CctrKS (i.e., C~>1) and that the nutrient concentration decays exponentially for CctrKS (i.e., C~1).

The position z~0>0, defined by C~(z~0)=1, is the (rescaled) vertical level across which the nutrient concentration transitions from the quadratic decay described in Equation (A2.3.5) to the exponential decay described by Equation (A2.3.4) with z~1=z~0. The nutrient penetration level we have defined is exactly HS=βz~0. Setting z~=z~0 in Equation (A2.3.5), we see that z~0=O(C~0) if C~01. Therefore,

HS=βz~0=O(D+C0Yρ0λS)=O(1λS).

Since the nutrient decays exponentially above HS, only those cells below the height HS grow with the maximum growth rate λS. Therefore, the speed of vertical expansion VH is given by VHHSλSλS; see the discussion in the main text: Vertical rise—quantitative analysis, in Simulation Results and Analysis.

https://doi.org/10.7554/eLife.41093.032

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
    Determination of biofilm diffusion coefficients using micro-electrodes
    1. EE Beuling
    2. JC Van den Heuvel
    3. SPP Ottengraf
    (1996)
    Progress in Biotechnology 11:31–38.
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
    Factors affecting the growth of bacterial colonies on agar plates
    1. AL Cooper
    2. AC Dean
    3. C Hinshelwood
    (1968)
    Proceedings of the Royal Society of London. Series B, Biological sciences 171:175–199.
    https://doi.org/10.1098/rspb.1968.0063
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
    Understanding molecular simulation (2nd Ed)
    1. D Frenkel
    2. B Smit
    (2002)
    New York: Academic Press.
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
    Ueber die Berührung fester elastischer Körper
    1. H Hertz
    (1882)
    J fur die Reine und Angew Math 1882:156–171.
  32. 32
  33. 33
  34. 34
    Contact Mechanics, Contact Mechanics
    1. KL Johnson
    (1985)
    Cambridge University Press.
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
    Investigation of the equation of diffusion combined with increasing of the substance and its application to a biology problem
    1. AN Kolmogorov
    2. IG Petrovsky
    3. NS Piskunov
    (1937)
    Bulletin of Moscow State University Series A: Mathematics and Mechanics 1:1–25.
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
    Theory of Elasticity
    1. LD Landau
    2. EM Lifshitz
    (1986)
    Oxford: Pergamon.
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
    Numerical Solution of Partial Differential Equations
    1. KW Morton
    2. DF Mayers
    (1995)
    Cambridge University Press.
  52. 52
  53. 53
    Cell growth, genome duplication and cell division
    1. N Nanninga
    2. CL Woldringh
    (1985)
    In: N Nanninga, editors. Molecular Cytology of Escherichia Coli. London, New York: . Academic Press. pp. 259–318.
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
    Localization of active microorganisms in cheese by autoradiography
    1. J Reyrolle
    2. F Letellier
    (1979)
    Applied and environmental microbiology 38:1162–1165.
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
    Granular flow down an inclined plane: bagnold scaling and rheology
    1. LE Silbert
    2. D Ertaş
    3. GS Grest
    4. TC Halsey
    5. D Levine
    6. SJ Plimpton
    (2001)
    Physical Review E, 64, 10.1103/PhysRevE.64.051302.
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
  82. 82
  83. 83

Decision letter

  1. Richard A Neher
    Reviewing Editor; University of Basel, Switzerland
  2. Naama Barkai
    Senior Editor; Weizmann Institute of Science, Israel
  3. Richard A Neher
    Reviewer; University of Basel, Switzerland

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Spatiotemporal establishment of growing bacterial colonies" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Richard Neher as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Naama Barkai as the Senior Editor.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

All reviewers agreed that your analysis of bacterial colony growth is an important step towards a quantitative understanding of this important aspect of microbiology. The manuscript clearly shows which mechanisms and forces are important in different locations and stages. Friction causes cells in the center of the colony to verticalize leading to vertical growth of the colony. This growth is restricted to a thin layer close to agar due to nutrient depletion. Radial growth occurs in a concentric ring at the periphery of the colony where cells are oriented horizontally. Your work nicely integrates observations, simulations, and analysis. However, there are a number of points that we would like to see addressed.

Essential revisions:

1) Verticalization of cells in the center of the colony is plausible and compelling in the simulation and the model, but strong evidence for this process in experiments is lacking. More direct evidence for cell orientation is particularly important as others (Su et al., 2012) have observed horizontal layers of cells. These authors used a constrained 2d geometry, but it is conceivable that surface tension has a similar effect. Enos-Berlage et al., 2000 provide some evidence in this direction. On a related semantic note, "buckling" is often more narrowly understood as the failure of a continuous material under compression – "verticalization" might be a better term.

2) A recent paper (Beroz et al., 2018) studies a very similar problem with interesting parallels and differences. The results you report should be compared to those by Beroz et al.

3) You stress the importance of incorporating both static and dynamic friction. To support this conclusion, we would like to see more explicitly that other/simpler models of friction lead to qualitatively different conclusions at odds with the experimental observations. An extended discussion of the implications of different forms of friction for colony or biofilm growth would be welcome.

4) The differences between dense homogeneous colonies and biofilms should be articulated more clearly. Many biofilms contain widely separated bacteria connected by polymeric substances. What concrete lessons can be learned from colony growth about biofilm growth?

5) The source code should be made available. Ideally on a repository like github.

https://doi.org/10.7554/eLife.41093.039

Author response

Essential revisions:

1) Verticalization of cells in the center of the colony is plausible and compelling in the simulation and the model, but strong evidence for this process in experiments is lacking. More direct evidence for cell orientation is particularly important as others (Su et al., 2012) have observed horizontal layers of cells. These authors used a constrained 2d geometry, but it is conceivable that surface tension has a similar effect. Enos-Berlage et al., 2000 provide some evidence in this direction. On a related semantic note, "buckling" is often more narrowly understood as the failure of a continuous material under compression – "verticalization" might be a better term.

At a first glance, the formation of horizontally oriented cell layers observed by Su et al., 2012, is in apparent contradiction to vertically oriented cells in the interior of a colony predicted by our model. However, Su et al.’s images were taken for very thin colonies (no more than a few layers in thickness). For such thin colonies, our model actually predicts cells to lie mostly horizontally; see Author response image 1A and B. According to our model, appreciable verticalization of cells occurs only after the colony thickness reaches certain height of several layers of cells; see the inset of Author response image 1C.

Author response image 1
Lack of verticalization for thin colonies.

(A) Reproduction of Figure 2A from (Su et al., 2012). (B) Bottom view of a simulated colony at t=7h. (C) The fraction of vertically oriented cells (defined by the angle with the z-axis less than 45 degrees) vs. height from the simulation.

https://doi.org/10.7554/eLife.41093.036

We have added Figure 4—figure supplement 2, which reproduces Author response image 1B and C, and also the following sentence in “Vertical rise – a pictorial view” in the main text: (See also our response to the next comment.)

“However, the internal verticalization took some time to develop (Figure 4—figure supplement 2); appreciable fraction of cells (50%) picked up vertical orientation only when the radius reached 250μm.”

We thank the reviewer for pointing out the very interesting study by Enos-Berlage et al. which we were not aware of previously. These authors showed very clear confocal images of cells standing up. However, it is very difficult for us to obtain images of similar qualities for our colonies. Before we go into our results, let us explain the experimental issues involved since some of the reviewers are theorists.

1) They put a slide on top of the colony, which enabled them to use a (probably) water objective to get better resolution. If we do that to our colony, our colony collapses. Why the difference? Their colonies have EPS, which tends to make cells stick together more. (They say it has good structural integrity, because of EPS).

2) Because our colonies produce no EPS, our colonies are much more densely packed and appear opaque. It is therefore more difficult to distinguish between different cells.

3) We were only able to take high resolution image of our colonies from the bottom, through a layer of agar. To do so, we had to use long working distance objectives, which compromises z-resolution. Also, we had to significantly reduce the agar thickness, with possible physiological consequences.

Despite the above difficulties, we attempted to directly visualize cellular orientations in the interior of our colonies. As seen in Author response image 2, cells near the periphery of the colony clearly have elongated shape while cells in the interior have round shape. This is consistent with the predicted vertical orientation in the interior. However, as mentioned above, in order to take confocal images at single-cell resolution, we had to use very thin agar plates (~0.3mm thickness). This significantly perturbs the colony growth characteristics, e.g., the radial and vertical growth, from the bulk of the results described in this study (which was done with 7mm thick agar). Many parameters can be affected for such different experimental settings, agar dryness, surface roughness, nutrient depletion, etc. Consequently, we do not feel comfortable to include these images in the manuscript, even though they support our prediction. Instead, we prefer to leave the result of verticalization in compact colony as a prediction, to be validated experimentally in details in a better designed future study.

Author response image 2
Experimental results on the confocal images of a bacterial colony at different positions, from the center to the edge.

The bacterial strain and the growth medium are the same as described in “Experimental Methods” in the main text, except that the thickness of agar dish is ~0.3mm. The picture was taken at roughly 24 hours after the initial inoculation.

https://doi.org/10.7554/eLife.41093.037

Finally, we thank the reviewer for suggesting the term “verticalization”. We now refer to it for the microscopic cell flow in the colony. However, we keep the term “buckling” to refer to the mechanical event itself, and still use “buckling width” to refer to the single layered, concentric ring region in the colony periphery.

2) A recent paper (Beroz et al., 2018) studies a very similar problem with interesting parallels and differences. The results you report should be compared to those by Beroz et al.

We believe the driving force for verticalization in our model is the same as that in the model of Beroz et al. It is a buckling instability driven by in-plane mechanical compression fueled by cell growth. The barrier for verticalization is however different. In Beroz et al. it is a cell-surface adhesion. In our case, it is the surface tension that leads to an effective force driving the colony into the agar. The two different mechanisms lead to quantitative differences in the process of verticalization. While Beroz et al. found cells to stand up even within a mono-layer, our model shows substantial verticalization only when the colony grows to a sufficient thickness, which is consistent with experimental observations by us and by Su et al.; see response to Essential revision 1 and Author response image 1.

With that said, however, we note that while our model applies directly to the compact colony cluster studied in our experiments, the experimental system that Beroz et al. modeled was actually a loosely packed colony. As far as we understand, cell-cell adhesion not modeled in Beroz et al. played a big role as described in Yan et al., 2016. It seems to us that additional assumptions are needed to justify the application of Beroz et al.’s model to their experimental systems.

We further note that, while the main focus of Beroz’s paper is on verticalization, the main focus of our study is to characterize and describe the temporal dynamical development of simple bacterial colonies in 3D. In particular, we focus on how radial expansion and vertical growth depend on various environmental parameters. We regard our key finding to be the elucidation of the origin of the linear radial expansion, which was long thought to be due to nutrient limitation. Instead, we show here that it is mechanical in origin, resulting from an interplay of colony surface tension and cell-agar friction.

We have added the following sentences in the section Discussion in the main text to describe the similarities between our work and that by Beroz et al. on the cell verticalization:

“Cell verticalization has been observed experimentally for Vibrio parahaemolyticus (Enos-Berlage and McCarter, 2000) and for Vibrio Cholerae (Beroz et al., 2018; Yan et al., 2016). […] Due to the different energy barriers against verticalization, the length scales of verticalization between our model and that of (Beroz et al., 2018) are very different: The colonies in Beroz et al. spread very slowly radially (~3μm/h), and verticalization occurs at a colony radius of 5-10μm. Colonies in our model spread much faster (~14μm/h), and substantial verticalization occurs at a radius of ~250μm; see Figure 4—figure supplement 2.”

3) You stress the importance of incorporating both static and dynamic friction. To support this conclusion, we would like to see more explicitly that other/simpler models of friction lead to qualitatively different conclusions at odds with the experimental observations. An extended discussion of the implications of different forms of friction for colony or biofilm growth would be welcome.

Introducing the dynamic friction was motivated by the strong dependence of the radial expansion speed on the growth rate observed (Figure 1H of the main text, reproduced as Author response image 3A). As can be seen from our simulation results (filled circles in Figure 8F of the revised manuscript, reproduced as Author response image 3B), the model with dynamic friction captures the experimental observations well. If we adopt a model with the static friction alone, we would instead obtain results shown in Author response image 3C (filled circles), which is the same as blue triangles in the new Figure 8F in the main text, with a much weaker growth-rate dependence.

Author response image 3
The dynamic and static friction models distinguished by the growth-rate dependence of radial expansion speed.

(A) From main text Figure 1F. Experimental results on the radial speed VR and the vertical speed VH with respect to various batch culture growth rates. (B) From Figure 8F in the main text (new version). Simulation results of radial speed VR of the colony, where both static and dynamic frictions are included in simulations. (C) New simulation results of radial speed VR of the colony, where only the static friction is included.

https://doi.org/10.7554/eLife.41093.038

The reason that dynamic friction is called for here is simple. A proportional relation between the radial expansion speed VR and growth-rate λs is obtained if the buckling width is the same in different growth conditions. (It is strictly a proportional relation for cells with a fixed division length, as shown by the open circles in Figure 8F in the revised version, because the relation becomes superlinear when the growth-rate dependence of the cell division length itself is taken into account, as shown by the filled circles.) The latter is produced by the dynamic friction, since the frictional force depends only on the normal force, rather than on the relative speed between cell and agar. In contrast with the static friction, the faster cells grow, the faster is the relative speed between cell and agar, hence stronger the frictional force. This leads to reduced dependence of VR on the growth-rate (compare Author response image 3B and Author response image 3C).

We have added the blue triangles in Figure 8F, same as Author response image 3C. In addition, we have added the following sentence in the section “Radial expansion – quantitative analysis” in “Simulation Results and Analysis”:

“In contract, static friction leads to a much weaker dependence of VR on λS (blue triangles in Figure 8F).”

We have also added the following sentence in Conclusion (in the paragraph where we discuss the fiction forces) in the main text:

“Indeed, in a model with static friction alone, a much weaker growth-rate dependence of radial expansion speed was obtained (Figure 8F blue triangles).”

4) The differences between dense homogeneous colonies and biofilms should be articulated more clearly. Many biofilms contain widely separated bacteria connected by polymeric substances. What concrete lessons can be learned from colony growth about biofilm growth?

Due to complex EPS, biofilm is notoriously difficult to model. In this work, we focus on bacterial colonies without EPS, in order to confront predictions based on model with key known ingredients (growth rate, nutrient diffusion, elastic cell-cell interactions, and friction forces) with quantitative experiments. Key results we learned from this study all shed light on the more complex biofilm dynamics:

1) One key finding is that radial growth of our colonies is not limited by nutrient as commonly believed, but by the interplay of surface tension and cell-agar friction. Given that biofilms have typically much lower bacterial densities, nutrient limitation will be even less of a problem. Also, the EPS secreted by the bacteria could modify both the surface tension and cell-agar friction to control the radial expansion speed.

2) Nutrient supply is limiting for the vertical growth of our colonies. This becomes less of a problem for the loosely packed biofilms. Moreover, biofilms are said to form channels in their interior, which will further alleviate the supply of nutrient, thereby allowing for faster vertical expansion.

3) Verticalization of cells in the interior, which is important for vertical growth but occurs at rather large colony sizes according to our model (Author response image 1C), also occurs in biofilms but at much smaller colony sizes as described by Yan et al., 2016, and Beroz et al., 2018. While the precise nature of the forces driving verticalization may be different in the two cases, the underlying origin appear to be similar – mechanical instability due to in-plane compression resulting from colony expansion and cell-agar friction.

In light of the above, we can see that the additional ingredients in biofilms provide the colonies with the means to expand even faster both horizontally and vertically.

We have added the following in the section Discussion (second to the last paragraph) of the main text:

“While our work is exclusively on bacterial colonies without EPS, key results we learned from this study shed light on the more complex dynamics of heterogeneous biofilms. […] In light of these comparisons, we see that the additional ingredients provided by biofilms enable the colonies to expand faster both horizontally and vertically.”

5) The source code should be made available. Ideally on a repository like github.

Our computer codes are developed using C++ with OpenMP parallelization. We have placed the major and basic parts of our codes in the repository GitHub. The address is: https://github.com/huiprobable/CellsMD3D

https://doi.org/10.7554/eLife.41093.040

Article and author information

Author details

  1. Mya R Warren

    Department of Physics, University of California, San Diego, La Jolla, United States
    Contribution
    Conceptualization, Investigation, Methodology
    Contributed equally with
    Hui Sun
    Competing interests
    No competing interests declared
  2. Hui Sun

    1. Department of Physics, University of California, San Diego, La Jolla, United States
    2. Department of Mathematics, University of California, San Diego, La Jolla, United States
    3. Department of Mathematics and Statistics, California State University, Long Beach, Long Beach, United States
    Contribution
    Investigation, Visualization, Methodology, Writing—original draft
    Contributed equally with
    Mya R Warren
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3137-3873
  3. Yue Yan

    1. Department of Mathematics, University of California, San Diego, La Jolla, United States
    2. School of Mathematical Sciences, Fudan University, Shanghai, China
    Present address
    School of Mathematics, Shanghai University of Finance and Economics, Shanghai, China
    Contribution
    Investigation, Visualization, Methodology
    Competing interests
    No competing interests declared
  4. Jonas Cremer

    Department of Physics, University of California, San Diego, La Jolla, United States
    Present address
    Groningen Biomolecular Sciences and Biotechnology Institute, University of Groningen, Groningen, Netherlands
    Contribution
    Investigation, Visualization, Methodology
    Competing interests
    No competing interests declared
  5. Bo Li

    Department of Mathematics, University of California, San Diego, La Jolla, United States
    Contribution
    Supervision, Writing—original draft, Writing—review and editing
    For correspondence
    bli@math.ucsd.edu
    Competing interests
    No competing interests declared
  6. Terence Hwa

    Department of Physics, University of California, San Diego, La Jolla, United States
    Contribution
    Supervision, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    hwa@ucsd.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1837-6842

Funding

National Science Foundation (DMS-1319731)

  • Bo Li

National Science Foundation (DMS-1620487)

  • Bo Li

Natural Sciences and Engineering Research Council of Canada

  • Mya R Warren

Simons Foundation (542387)

  • Terence Hwa

Simons Foundation (522790)

  • Hui Sun

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was partially supported by the NSF through grants DMS-1319731 and DMS-1620487 to BL, by Simons Foundation through grants #542387 to TH and #522790 to HS. MRW was supported by a Canadian Natural Sciences and Engineering Research Council post-doctoral fellowship. This work used the NSF Extreme Science and Engineering Discovery Environment (XSEDE) through grant ACI-1548562. We thank Agnese Seminara for useful discussions. We also thank the reviewers to bring to our attention the references Beroz et al. (Nat. Phys. 2018) and Enos-Berlage and McCarter (J. Bacteriol. 2000).

Senior Editor

  1. Naama Barkai, Weizmann Institute of Science, Israel

Reviewing Editor

  1. Richard A Neher, University of Basel, Switzerland

Reviewer

  1. Richard A Neher, University of Basel, Switzerland

Publication history

  1. Received: August 15, 2018
  2. Accepted: February 20, 2019
  3. Version of Record published: March 11, 2019 (version 1)

Copyright

© 2019, Warren et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 745
    Page views
  • 168
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    2. Microbiology and Infectious Disease
    Nir Drayman et al.
    Research Article
    1. Computational and Systems Biology
    2. Physics of Living Systems
    Aby Joseph et al.
    Tools and Resources