Spatiotemporal establishment of dense bacterial colonies growing on hard agar

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.


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.

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 V R » 45:2 m=h and the colony height increased also linearly at a vertical speed V H » 12: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 14 t 24 h, where both the radial and vertical growth are linear.
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, 1981;Wimpenny, 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.

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 l 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 C s ¼ 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 t ! 12 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 following figure supplement is available for figure 3: Figure supplement 1. The spatially varying cell density (per unit volume of colony) is related to the spatially varying cell volume fraction f by ¼ f cell , where the volume fraction f is defined as the volume of all cells in a unit volume of the colony and cell is the constant mass density of a typical mature cell; cf. Appendix 1.2 on nutrient update. Figure 3 continued on next page 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, t 6 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 V R » 18 m=h and vertical ascending speed V H » 6 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'.

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 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 mm.
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 mm Â 4 mm 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 v ! r ; v z r ! 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 V z z ð Þ ¼ v z 0; 0; z ð Þ, in Figure 4E. We see that V z increases through a thin  region of height H S » 10 m. The vertical ascension speed is saturated for z > H S , meaning that above this thickness H S , cells move up steadily. Another way to visualize the vertical growth of the colony is to show the 'age' of cells in a crosssectional 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 mm 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 V colony / R 2 H / R 3 . Assuming that the increase of the colony size is due to a portion of cells growing at the maximal possible rate (l S ) in a growth zone of volume V growth t ð Þ, then d dt V colony / V growth t ð Þ leads to V growth / R 2 , 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, C ctr z ð Þ C 0; 0; z ð Þ, at various times t during colony growth. In the linear growth regime (for t ! 12 h), the profile C ctr z ð Þ 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 H S such that C ctr H S ð Þ ¼ K S , the Monod constant appearing in Equation 3; see Appendix A2.3. Since the local growth rate drops substantially where the nutrient concentration is below K S , we can take the value H S as the thickness of the vertical growth zone, leading to the vertical ascending speed: Detailed analysis of Equations 1 and 3 in Materials and methods shows that the scale of the stationary profile is set by ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi D þ =l S p ; cf. Appendix A2.3. This is verified in Figure 5B where the stationary profile C ctr z ð Þ is computed by repeating the simulation for different growth rate l S : the profile collapses into the same curve for different values of l S when plotted against z ffiffiffiffiffi l S p ; 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 H S / 1= ffiffiffiffiffi l S p for faster growth (due to stronger nutrient depletion), leading to a sublinear dependence of the vertical ascending speed, V H / ffiffiffiffiffi l S p . 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 l 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 l S < 0:4 h À1 , is difficult to realize by changing carbon sources. However, if we just fit the data in Figure 5C for l S > 0:5 h À1 , 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.
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 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 K S 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).
To elucidate the determinants of buckling, we plot in Figure 8A the azimuthal-averaged radial velocity profile V r Dr ð Þ for the bottom layer of cells, for a range of (signed) distances Dr into the edge of the colony; see Appendix 1 Equations (A1.5.1 and A1.5.2) for the definitions of V r and Dr. This radial velocity profile, which is stationary for t ! 12 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 V R is simply V r at Dr ¼ 0, we see that the width of this annulus together with the slope of V r Dr ð Þ set the radial expansion speed of the entire colony. To understand what goes on in this peripheral region, we examine in Figure 8B the azimuthalaveraged height profile of the colony, h Dr ð Þ, which is also stationary for t ! 12 h, with h Dr ð Þ » 1 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 V r 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 r ! Á V ! ¼ l, where l 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 V z is negligible ( Figure 8C). Hence, V r satisfies 1 r rV r ð Þ 0 ¼ l: Throughout the periphery region, the growth rate is essentially the maximal growth rate, that is l » l S , since the nutrient concentration in this region is set by the boundary value C s which well exceeds the Monod constant K S ; cf. Figure 8-figure supplement 1. Solving the above equation in the annulus in the limit jDrj ( R yields a linear dependence, where we used the definition of the radial expansion speed 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 V r 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 W b , defined as the largest value of jDrj where V r Dr ð Þ ¼ 0, sets the radial expansion speed since V R / l S Á W b . 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 V R ; while the interior cells that are away from the colony edge by this buckling width flow up vertically; cf. 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 Dr ð Þ 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 6, pressure is expressed in unit of P 0 ¼ g surf =w 0 . (E) The azimuthally averaged height H vs. Dr at growth rate l S ¼ 0:5 h À1 and 1:0 h À1 . (F) The simulated colony horizontal expansion speed V R vs. the batch culture growth rate l S with a fixed l div (red open circles) and a growth-rate dependent l div (red closed circles) using dynamic friction. The dashed line that fits the open circles is given by V R ¼ 16:9l S þ 0:8; the solid line that fits the closed circles for l S ! 0:5 h À1 is V R ¼ 5:7l S l div À 0:2; the dash dotted line fits the closed circles with V R ¼ 22:1l 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 V R vs. l S with a growth-rate dependent l div , for a model with static friction alone (see Equations (A1.4.6) of Appendix 1) between cell and agar (blue triangles). DOI: https://doi.org/10.7554/eLife.41093.017 The following figure supplements are available for figure 8:  Figure 4A, overlaid with coarse-grained velocity field (zoomed in view of the same periphery region in Figure 4D). DOI: https://doi.org/10.7554/eLife.41093.019 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 mm 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 W b at different radial expansion speed V R : Static friction would increase for increased V R , leading to decreased buckling width, whereas dynamic friction would not be affected by the radial expansion speed. Experimentally, we characterized V R in sugars supporting different growth rates l S , and V R is seen to increase linearly with l S (Figure 1I), suggesting a constant W b , 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 l S . The buckling width W b is indeed not dependent on l S ( Figure 8E), and the radial expansion speed V R is indeed linear in l S (open red symbols and dashed red line, Figure 8F). In contract, static friction leads to a much weaker dependence of V R on l S (blue triangles in Figure 8F).
The linear dependence on l 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 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 l S , we recover a nonlinear dependence of V R on l 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 l S > 0:5 h À1 (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 H S / 1= ffiffiffiffiffi l S p does not depend on l div (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 V R , but not the vertical ascension speed V H . Conversely, we expect changing the nutrient concentration C s to change the vertical nutrient penetration length and hence V H , but not V R . These expectations are indeed reproduced by the full colony simulation using different parameter values from the ones discussed so far, with the nutrient concentraiton C s doubled in Figure 9A (only V H increased), and with the cell-agar friction quartered in Figure 9B (only V R 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. We note that the actual values of radial and vertical expansion speeds obtained (V R ¼ 17:2 m=h and V H ¼ 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 V R and V H 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 V R and V H 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, level (with ca ¼ 0:8; cc ¼ 0:1; g cc;t ¼ 10000 m À1 h À1 , and h ran ¼ 0:1 m), and use C s ¼ 0:5 mM as the lower glucose concentration, C s ¼ 1:0 mM as the higher glucose concentration. In panel (B), we fix the glucose concentration at the lower level (C s ¼ 0:5 mM), and vary the friction, from the higher value of ca ¼ 0:8; cc ¼ 0:1; g cc;t ¼ 10000 m À1 h À1 , h ran ¼ 0:1 m used in (A) to the lower value of ca ¼ 0:2; cc ¼ 0:025; g cc;t ¼ 2500 m À1 h À1 , h ran ¼ 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 V H in panel C is consistent with a square root dependence on nutrient concentration (blue line) expected from the basic analysis in Figure 5. DOI: https://doi.org/10.7554/eLife.41093.020 The following source data is available for figure 9: Source data 1. Experimental data on the horizontal and vertical colony expansion speeds at various glucose and agar concentrations. DOI: https://doi.org/10.7554/eLife.41093.021 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.
Capturing all these behaviors within a single model and with a fixed set of parameters is a nontrivial 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 cellcell 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 vertically oriented cells planarly oriented cells cells in slanted orientation cell velocity Figure 10. 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 (H S ) 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 (W b ) 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, H S and W b , set the speeds of radial and vertical expansions, V R and V H , respectively, as shown in red. The growth rate dependence of these parameters is shown in blue. DOI: https://doi.org/10.7554/eLife.41093.022 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 (~3m=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., 2007;Hallatschek 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 similarmechanical 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.

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 'km r :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 l-Red recombinase. The cells were incubated with shaking at 37˚C for 1 hour and then applied onto LB +Km agar plates. The Km r colonies were verified for the 'km r :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 'km r :rrnBT: Ptet-gfp' in EQ42 was then transferred to EQ54 (that is NCM3722DmotA) (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 NH 4 Cl 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 OD 600 » 0: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 OD 600 » 10 À6 . 10 L of culture (containing approximately 10 cells) was spread over prewarmed 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 n i;j . Pixel i; j ð Þ is assigned as type 1 if 16 > n i;j ! 3, and as type 2 if n i;j ! 16, indicating the pixel belonging to single-or multi-layer region, respectively. We then estimated the inner radius r inner and outer radius r outer of the colony by the formulas r inner ¼ r m=px ffiffiffiffiffiffiffiffiffiffiffiffiffiffi N px2 =p p and r outer ¼ r m=px , where N px1 and N px2 are the total numbers of pixels of type 1 and 2, respectively, and r m=px » 0:84 is the ratio of mm 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, t l , 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 t s separated by several hours. Growth curves extending over a period of multiple days could be obtained by stitching together R t À t l À t s ð Þ 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 G 01 to be defined below; see Figure 11. Nutrient diffuses in the two compartments, agar and colony, according to the diffusion equations 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 (G 01 in Figure 11) is defined by thresholding the density; see Appendix 1.2. The local mass growth rate l ¼ l r ! ; t is given by where l S is the batch culture growth rate for cells in a medium saturated with some sugar S, and K S 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: 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 G 01 ; G 02 ; and G b ; and the Dirichlet boundary condition C ¼ C s on the lateral wall of agar region G s ; cf. Figure 11. The parameter C s 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) w 0 on its two ends; see Figure 12A. For a given cell i at a given time t, we use a position vector r  Table S5 for typical values of L, a, and b used in our simulations. The colony surface or colony-air interface G 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 G 12 . The remaining part, denoted G 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 G 01 , and the colony-agar interface G 02 depend on time t. DOI: https://doi.org/10.7554/eLife.41093.023 to denote the center-of-mass of the cell, a unit vector n ! i t ð Þ to denote its orientation (the direction along the cylindrical axis), and l i t ð Þ to denote the length of the cylinder between the two centers of hemispherical caps. Each cell starts off with the same cylindrical length l 0 . We assume that during cell growth, the width w 0 is fixed, and the cylinder length of a cell increases at a rate l t ð Þ: We call this the cell elongation rate. It is proportional to the mass growth rate l t ð Þ of the cell with a geometrical proportionality factor s: l ¼ sl; see Appendix 1.3.
The mass growth rate is calculated based on the nutrient concentration at the center r . The growth of cylindrical length l i t ð Þ is then given by the growth equation Once the cylindrical length l i t ð Þ reaches a critical value l div , the cell divides into two daughter cells with cylindrical lengths being l 0 with small fluctuation; see Figure 12B and Appendix 1.3. For different growth media supporting different growth rates l S , the value of l div is in general growth-rate dependent 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 v ! i and angular velocity ! ! i , which follow Newton's second law where M i and I i are the mass and moment of inertia of the cell, and F ! net i and T ! net i 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 G 01 in Figure 11

(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 d cc in their sphero-cylinder cell boundaries; see Figure 13A. At the point of contact, the cell-cell interaction force F ! cc is decomposed into the normal and tangential components, of magnitudes F cc;n and F cc;t as defined in and Appendix 1.4a. The normal force includes the Hertzian elasticity force with magnitude F cc;elas / ffiffiffiffiffi ffi w 0 p d 3=2 cc (Hertz, 1882;Johnson, 1985). Additionally, the normal and tangential force each has a dissipation component, of magnitude F cc;disp; n and F cc;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.
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, F cc;disp; n and F cc;disp;t , respectively, are taken to be proportional to the normal and tangential components of the relative cell-cell velocity, v cc;n and v cc;t . We use F cc;disp;n / d 1=2 cc v cc;n and F cc;disp;t / d cc v cc;t , where the additional dependences on the overlap d 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 F max cc;disp;t ¼ cc F cc;elas , where cc is the dynamic frictional coefficient; see Appendix 1.4a. Thus, the full cell-cell interaction force is given by F cc;n ¼ F cc;elas þ F cc;disp;n ; (7a) F cc;t ¼ min F cc;disp;t ; cc F cc;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, F ! ca , 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 h ran 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 F ! visc , 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. 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 F surf 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, whose magnitude is proportional to the width of the cell w 0 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 h x; 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) h w x; y; t ð Þ by adding dh to the colony height. This thickness dh 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 h cell exceeds the liquid height h w , we impose a restoring force normal to the liquid surface; see Figure 14C. As the magnitude of the saturating restoring force F surf;0 (Equation 8) is independent of the height Figure 14. 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. G 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 h cell 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 h w ) and colony height (indicated by h), respectively. The parameter d 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. DOI: https://doi.org/10.7554/eLife.41093.026 difference dz h cell À h w for dz > 0, the restoring force can be mathematically written as F surf ¼ F surf;0 Á u dz ð Þ, where u is the Heaviside step function. To avoid numerical instability, we make a linear extrapolation between F surf ¼ 0 and F surf ¼ F surf;0 over a narrow transition region 0 < dz < w 0 =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 F ! net i and the corresponding torque T ! net i are calculated. Moreover, the pressure on the cell i can be also calculated as where V i is the volume of cell i, the index j runs through all the different forces experienced by cell i, and r ! ji 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 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 m 2 =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 m 2 =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 g CDW =g glucose (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 g CDW =m 3 being the cell dry weight density (Basan et al., 2015). The Monod constant is taken to be that for glucose, K S ¼ 20 M (Monod, 1949). The cell dividing length l div ¼ 3 m and the diameter w 0 ¼ 1 m are fixed in all the simulations unless otherwise indicated. The constant value C s of nutrient concentration in the boundary conditions for the diffusion equations and the batch culture growth rate l 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 Dt. 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).
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Data availability
The simulation data files are large, hence we do not include it. All the simulation data can be generated from running the source code in GitHub (https://github.com/huiprobable/CellsMD3D; copy archived at https://github.com/elifesciences-publications/CellsMD3D). Experimental source data files are provided for Figure1, Figure 1-figure supplement 1, and Figure 9.
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 G 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 G 12 . The remaining part, denoted G 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 G 01 , and the colony-agar interface G 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 h grid : 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; h ran with h ran 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 Dt, 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 þ Dt 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 100-1000 time steps of calculations of cell growth, division, and movement.

Cell density and colony boundary
Given the positions of all the cells at time t, we define the (local) volume fraction f ¼ fðr; tÞ of the cells in each finite difference grid box above z ¼ 0 by f ¼ 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 ðr; tÞ ¼ fðr; tÞ cell ; (A1.2.1) 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 fð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: The first two equations, Equation (A1.2.2) and Equation (A1.2.3), are the reactiondiffusion 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, l S is the batch culture growth rate, K S 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 G 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 G 01 , the colony-air interface G 02 , and the bottom of agar G b , we impose the flux-free boundary condition, with q=qn 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 G s , the lateral faces of the agar region, we prescribe a constant value C s 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 qC=qt 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., 2013;Morton 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 G b , and then from grids on the bottom G 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 err conv ; 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., 2013;Morton 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. 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 l ¼ lðr; tÞ by lðr; tÞ ¼ l S Cðr; tÞ Cðr; tÞ þ K S : (A1.2.8) This is first defined at each grid pointr 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 byp andq the centers of the two hemispheres. We call ' ¼ kp Àqk the cell cylindrical length which can vary with time t. We also denote byñ ¼ ðq ÀpÞ=' the unit vector pointing from one center of hemispherep to the otherq, and call it the direction of the cell. Note that the center of mass of the cell isr c ¼ ðp þqÞ=2: We denote by w 0 the diameter of each of the hemispheres. The volume and mass of the cell are given by 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 w 0 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. Here,lðr c ; tÞ is the cell elongation rate evaluated at the cell centerr c . The elongation ratel is proportional to the mass growth rate l defined in Equation (A1.2.8):lðr; tÞ ¼ slð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 þ Dt by where Dt 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 hemispheresp andq grows long enough, with its cylindrical length ' ! ' div for some critical value ' div >w 0 , 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 where h 2 ðÀ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 ¼ 2' 0 þ w 0 ; cf. Figure 12B in the main text. Note that before division, the mother cell has volume pw 2 0 ' div =4 þ pw 3 0 =6; after division, the total volume of two daughter cells is pw 2 0 ' div =4 þ pw 3 0 =12. There is pw 3 0 =12 volume loss every time due to division. For a fixed cell aspect ratio ' div : w 0 ¼ 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 l 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)p 1 ¼p;q 1 ¼p À ' 1 'p Àq ð Þ;p 2 ¼q þ ' 2 'p Àq ð Þ;q 2 ¼q: Moreover, these daughter cells inherit the velocity from their mother cell. But the angular velocities of these two daughter cells are set to bẽ ! 1 ¼ ð0; 0; 0Þ and! 2 ¼ ! ran ð0; 0; Þ; respectively, where 2 ðÀ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 massr c of a daughter cell. We rotate the vectorp Àr c with p the center of a hemisphere of this daughter cell. (For notational simplicity, we user c andp 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 atr c and the three unit coordinate vectorsẽ 1 ,ẽ 2 , andẽ 3 byẽ 3 ¼p Àr c andẽ 2 ¼ẽ 3 Âẽ 1 , whereẽ x ;ẽ y ;ẽ z are the unit coordinate vectors in the original coordinate system. Then we rotate the vectorp Àr c by an angle of ' around the axisẽ 2 , and aroundẽ 3 successively. Here, ' is a random variable in the range ½0; ' ran p, where ' ran is a constant that sets the magnitude of fluctuations of ' (cf. Supplementary file 1- Table S3) and 2 ½0; 2pÞ is a random number. If we denote byp new the new center of hemisphere after the rotation, theñ p new Àr c ¼ kp Àr c k sin ' cos ẽ 1 þ sin ' sin ẽ 2 þ cos 'ẽ 3 ð Þ: This equation allows us to find the coordinates ofp new in the original coordinate system. The coordinates of the other center of hemisphereq new are given byq new ¼ 2r c Àp new : The two growth rates We assume that, in the life span of a cell, the cell mass growth rate l and the cell elongation ratel stay constant. The growth of mass satisfies the equation MðtÞ ¼ Mð0Þe lt . If the mass doubling time is t d , then Mðt d Þ ¼ 2Mð0Þ and hence lt d ¼ ln 2: 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 ¼ 2' 0 þ w 0 . From the growth equation Equation (A1.3.2), we have then ' div ¼ l 0 el td : Therefore, For the fixed dividing cell aspect ratio ' div : w 0 ¼ 3 : 1, we have s ¼ ln 3 : ln 2.

Cell movement
Consider a cell at some time instant. Let us denote its centers of hemispheres byp old andq old , its cylindrical length by ' old ¼ kp old Àq old k, its direction byñ old ¼ ðq old Àp old Þ=' old ; and its mass by M old : Let us also denote byṽ old 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 Dt: We first calculate the forceF half and torqueT half of the cell. Details of such calculations are given below in section A1.4. We then calculate the velocityṽ half and angular velocity! half for a half time step: where I old;n and I old;t are the moments of inertia of the cell along the directionsñ old 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 We now update the cell positions for the entire time step Dt by updating the centers of hemispheresp new ¼p old þ Dtṽ half þ! half Âp old Àq old 2 ; q new ¼q old þ Dtṽ half þ! half Âq old Àp old 2 : We update the force and torque of the new cell with the centers of hemispheresp new and q new by the procedure of force calculations described below in section A1.4 to getF new and

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 F static ¼ sFnormal andF dynamics ¼ 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 g, 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 forceF cc;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 d 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). We end this part with a method of computing the minimum distance d between the central cylindrical line segments g 1 and g 2 of the two cells and the corresponding points on these segments that reach this minimum distance. For notational convenience, let us denote byp i andq i the position vectors of the centers of hemispherical caps of the cell i with i ¼ 1 and 2, respectively. Letũ i ¼q i Àp i We parameterize the central cylindrical line segments g i byr i ðtÞ 1 p i þ tũ i for 0 t 1 Denotep 0 ¼p 1 Àp 2 and define the distance-square function f ðs; tÞ ¼ kr 1 ðsÞ Àr 2 ðtÞk 2 ¼ ksũ 1 À tũ 2 þp 0 k 2 :

Relative tangential velocity
Clearly, f is minimized in ½0; 1 Â ½0; 1 by some point ðs 0 ; t 0 Þ 2 ½0; 1 Â ½0; 1: The minimum distance d > 0, and the two pointsã 1 andã 2 on the two line segments reaching this distance are then given by We first assume that the two lines are not parallel, that isũ 1 Âũ 2 6 ¼0: In this case, f ðs; tÞ is a strictly convex and quadratic function with a constant symmetric positive definite Hessian matrix. By setting q s f ðs; tÞ ¼ 0 and q t f ðs; tÞ ¼ 0, we find that f is minimized in R 2 at ðs 0 ; t 0 Þ with s 0 ¼ ðp 0 Âũ 2 Þ Á ðũ 1 Âũ 2 Þ kũ 1 Âũ 2 k 2 and t 0 ¼ ðp 0 Âũ 1 Þ Á ðũ 1 Âũ 2 Þ kũ 1 Âũ 2 k 2 : If ðs 0 ; t 0 Þ 2 ½0; 1 Â ½0; 1, then we obtain d andã 1 ; 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 ðs 0 ; t 0 Þ 2 ½0; 1 Â ½0; 1 and the corresponding minimum value of f on this square. Consider, for instance, the side s ¼ 0 and 0 t 1: The function f ð0; tÞ for all t 2 R is found to be minimized at t 1 ¼ ðã 0 Áũ 2 Þ=kũ 2 k 2 . If t 1 2 ½0; 1, then the minimum value of f on this side is f ð0; t 1 Þ: Otherwise, this value will be the minimum of f ð0; 0Þ and f ð0; 1Þ. We now assume that the two line segments g 1 and g 2 are parallel. In this case, the minimum distance d is achieved by the minimum of the distance fromp 1 to the second line segment g 2 and that fromq 1 to g 2 : d ¼ minfdist ðp 1 ; g 2 Þ; dist ðq 1 ; g 2 Þg: 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 forceF ca 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. Letp ¼ ðx a ; y a ; z a Þ be the center of the spherical cap corresponding to that end of the cell. We denote by d ca the indentation depth: d ca ¼ w 0 =2 À z a : We also denoter ca ¼ ðx a ; y z ; z a À w 0 =2Þ, which is the midpoint of the line segment along the vertical line passing through the pointp between z ¼ 0 and z ¼ Àd ca : As before, we denote byṽ ca the relative velocity at the centerr ca . It is given bỹ v ca ¼ṽ þ! Â ðr ca Àr c Þ; whereṽ is the velocity of the cell at its centerr c 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ñ ca ¼ ð0; 0; 1Þ: The tangential direction is defined through the relative velocityṽ ca bỹ t ca ¼ṽ ca À ðṽ ca Áñ ca Þñ ca kṽ ca À ðṽ ca Áñ ca Þñ ca k ; if the denominator is nonzero. (Otherwise,t ca can be any unit vector orthogonal toñ ca .) Similar to the cell-cell interaction, the total cell-agar interaction force exerted on the midpointr ca and the corresponding torque are given bỹ F ca ¼F ca;n þF ca;t ; T ca ¼ ðr ca Àr c Þ ÂF ca : HereF ca;n andF ca;t are the forces normal and tangential to the mean colony-agar interface z ¼ 0: They are given byF where k ca is an elastic constant in the Hertzian stress, g ca;n and g ca;t are static friction constants, ca is the dynamic friction constant, and M cell ¼ cell V cell is the cell mass with V cell the cell volume. Note that the factor ffiffi ffi 2 p 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 w 0 =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 ¼ g surf dA, where W is the work done by the water surface, g 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 1figure 3. Here, we approximate a cell by a sphere of diameter w 0 : (The notation for the radius is R here.) The surface area is A ¼ 2pRh and hence dA ¼ 2pRdh. (In the case of a disk in the two-dimensional setting, the change in area is dA ¼ 2prRd ¼ 2pR 2 sin d. Note that h ¼ R sin and dh ¼ R sin d. Hence dA ¼ 2pRdh:) As a result, the surface tension force is Appendix 1-figure 3. Schematic of the derivation of the surface tension force. 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 h c » 0:1w 0 , when it saturates to its maximum value of pw 0 g surf .
Specifically, let us consider a cell on top of the colony. Letr c denote the center, andp andq the two centers of hemispheres of such a cell; cf. Figure 14C in the main text. We define the total surface tension forceF surf and torqueT surf on this cell to be F surf ¼F surf;p þF surf;q ; T surf ¼ ðp Àr c Þ ÂF surf;p þ ðq Àr c Þ ÂF surf;q ; whereF surf;ĩ ¼ pw 0 g surf min max f0; ðhĩ À zĩÞñ s Áñ z þ dhg 0:1w 0 ; 1 & 'ñ s ;ĩ ¼p orq: Here, hĩ is the water level atĩ orq),ñ z ¼ ð0; 0; 1Þ is the unit vector along the z direction, zĩ 1 i Áñ z is the z-component ofĩ,ñ s is the unit vector of the colony surface, and dh is a constant determining how tightly the surface tension holds the cells. The water level is coarse-grained. First, on each horizontal grid square B i;j ¼ ðiDx; ði þ 1ÞDxÞ Â ðjDy; ðj þ 1ÞDyÞ, we set the water level, h i;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 B 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 h i;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 forceF visc exerted at a cell by the surrounding liquid and the corresponding torque T visc are given by (note that w 0 is the diameter)F visc ¼ À3p liq w 0ṽ andT visc ¼ Àp liq w 3 0! ; respectively, whereṽ 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: F cc;t ¼ Àg cc;t M eff d 1=2 cc jṽ cc Át cc jt cc andF ca;t ¼ Àg ca;t M cell d 1=2 ca jṽ ca Át ca jt ca ; (A1.4.6) 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 f, pressure field P, velocity fieldṽ, and director fieldñ over a subregion G of the colony region. Examples of a subregion G include: . The union of a few grid boxes for coarse-graining in the entire colony; . A small box at the agar-colony interface ðiDx; ði þ 1ÞDxÞ Â ðjDy; ðj þ 1ÞDyÞ Â ð0; DzÞ for some i, j for coarse-graining around such interface; . A small box in a vertical layer ðÀDx; DxÞ Â ðjDy; ðj þ 1ÞDyÞ Â ðkDz; ðk þ 1ÞDzÞ or ðiDx; ði þ 1ÞDxÞ Â ðÀDy; DyÞ Â ðkDz; ðk þ 1ÞDzÞ for some i, j, and k for coarse-graining around the cross-section x ¼ 0 or y ¼ 0, respectively; and . A cylindrical 'cube' ðidr; ði þ 1ÞdrÞ Â ðjd; ðj þ 1ÞdÞ Â ðkDz; ðk þ 1ÞDzÞ in the cylindrical coordinates ðr; ; zÞ for some small dr > 0, d > 0, and Dz > 0, and for some i, j, and k, for azimuthal coarse-graining. Now let us fix a subregion G in the colony. Let f denote the pressure P or a component of the velocity fieldṽ; and denote by f i such a quantity at the center of cell i. We define the coarse-grained average of f over G to be f ðGÞ ¼ wherer ci is the center of mass for the cell i. Similarly, we define the coarse-grained volume fraction over the subregion G to be fðGÞ ¼ where V i is the volume of cell i and VolðGÞ is the volume of the subregion G. Note that, if G is a grid box, then fðGÞ is the same as the volume fraction f on that box; cf. section A1.2. Note also that the cell density can be coarse-grained following its relation with the volume fraction f; cf. Equation (A1.2.1).
The director fieldñ 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ñðGÞ over a given subregion G to be a maximizer of the maximization problem max n X rc i 2G jñ Áñ i j 2 subject to kñk ¼ 1: By the Lagrange multiplier method, this leads to an eigenvalue problem for a 3-by-3 matrix, with the maximizerñðGÞ being a unit eigenvector that corresponds to the largest eigenvalue.
Let dr > 0 be small. Let N ! 2 be an integer and define d ¼ 2p=N . We denote G i;j;k ¼ ðidr; ði þ 1ÞdrÞ Â ðjd; ðj þ 1ÞdÞ Â ðkDz; ðk þ 1ÞDzÞ in the cylindrical coordinates ðr; ; zÞ; We define the azimuthal average of a scalar field f bỹ where f ðG i;j;k Þ is the the coarse-grained average of f over G 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 Dr ¼ r À RðÞ; where RðÞ ¼ max where ðr 0 ; Þ 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 0 j N , we denote by i j the largest integer not exceeding ðRðjdÞ þ DrÞ=dr. We then define the azimuthal average of the radial component V r of the velocity fieldṽ by 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 V R 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 V H : In Appendix 2-figure 1 below, we plot the velocities V R and V H for variable dividing length, and fit the simulation data with l S ! 0:5 h À1 . We observe that the straight line that fits simulated V R intersects the x-axis at~0:2 h À1 . This agrees well with the experimental result plotted in Figure 1H in the main text. Hence, the inclusion of l Sdependence on the cell dividing length can better describe experiment. Appendix 2-figure 1. Experimental measurement on V R and V H with various batch culture growth rates. Data are fitted with the stright line V R ¼ 51:27l S À 7:96 and the squre root curve V H ¼ 12:7 ffiffiffiffiffi l S p for the redial and vertical speeds of expansion V R and V H respectively. (B) Simulation results on V R and V H with various batch culture growth rates and with a variable diviting length l div . data are fitted for l S ! 0:5 h À1 using the straight line V R ¼ 22:1l S À 5:24 the squre-root curve V H ¼ 6:12 ffiffiffiffiffi l S p ; respectively. We now show that using a fixed cell dividing length ' div independent of the batch culture growth rate l S is a reasonable simplification for our simulations. We fixed the batch culture growth rate l S ¼ 1:0 h À1 and the constant concentration in the boundary condition C s ¼ 2:0 mM: We then distributed randomly 625 cells at time t ¼ 0 h: At t ¼ 9:0 h; 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:0 h to t ¼ 15 h: 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 l S , to low values, smaller than 10% of l 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:0 h to t ¼ 16 h: 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 l 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 l div and w 0 for the entire simulation.