We start by examining how fast the colony expands radially and vertically. We run a simulation with the batch culture growth rate ${\lambda }_{\mathrm{S}}=1.0{\mathrm{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_{\mathrm{s}}=0.5\mathrm{m}\mathrm{M}$ 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\ge 12\mathrm{h}$, suggesting a constant expansion of the colony in the radial and vertical directions by $t=12\mathrm{h}$, similar to the experimental profiles in Figure 1G. (The spatial cell density inside the colony is constant, ~0.68 ${\rho }_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}$, throughout the interior of the colony; see Figure 3—figure supplement 1.) Detail of the profile at the colony periphery appears to be different. This is due to an approximate height assignment based simply on thresholding the fluorescence intensity to obtain the global height profile. This thresholding procedure does not capture height at the periphery where it is one to a few layers in thickness. Figure 3C provides a quantitative picture of the colony radius (R, defined as the average radius of the bottom layer of the colony) and colony height (H, defined as the height at the center of the colony). At early time, $t\le 6\mathrm{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\mathrm{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_{\mathrm{R}}\approx 18\mathrm{\mu }\mathrm{m}/\mathrm{h}$ and vertical ascending speed $V_{\mathrm{H}}\approx 6\mathrm{\mu }\mathrm{m}/\mathrm{h}$; see Figure 3C. Thus, our model captures the linear increase of both the colony radius and height observed experimentally (Figure 1H). To understand the origin of these behaviors, we will analyze below the model output, first pictorially and then quantitatively. The lower numerical values of the speeds obtained from simulations are due to parameter settings chosen to limit computational time; this will be discussed in 'Parameter dependence'.

Figure 3 with 1 supplement see all

Download asset Open asset

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\mathrm{h}$ which is well in the steady linear growth regime. The colors distinguish the gross orientations of the cells. The model shows that cells near the top surface are oriented parallel to the colony surface (shown in cyan), while cells away from the top surface are mostly oriented vertically (shown in yellow). A detailed view of the top surface of the colony generated from the simulation is shown in Figure 4—figure supplement 1A. This prediction is validated by confocal scan of the colony in experiment as shown in Figure 4—figure supplement 1B.

Figure 4 with 2 supplements see all

Download asset Open asset

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

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

Next, we examine the coarse-grained velocity field $\overrightarrow{v}\left(\overrightarrow{r}\right)={(v}_x\left(\overrightarrow{r}\right),v_y\left(\overrightarrow{r}\right),v_z\left(\overrightarrow{r}\right))$ 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\left(z\right)=v_z(\mathrm{0,0},z)$, in Figure 4E. We see that $V_z$ increases through a thin region of height $H_{\mathrm{S}}\approx 10\mathrm{\mu }\mathrm{m}$. The vertical ascension speed is saturated for ${\displaystyle z>H_{\mathrm{S}}}$, meaning that above this thickness $H_{\mathrm{S}}$, cells move up steadily.

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

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

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 ${\displaystyle 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 ${\displaystyle V_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{y}}\propto R^{2}H\propto R^3}$. Assuming that the increase of the colony size is due to a portion of cells growing at the maximal possible rate (${\lambda }_{\mathrm{S}}$) in a growth zone of volume $V_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}\left(t\right)$, then ${\displaystyle \frac{d}{dt}{V}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{y}}\propto V_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}\left(t\right)}$ leads to ${\displaystyle V_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}\propto 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_{\mathrm{c}\mathrm{t}\mathrm{r}}\left(z\right)\equiv C(0,0,z)$, at various times t during colony growth. In the linear growth regime (for $t\ge 12\mathrm{h}$), the profile $C_{\mathrm{c}\mathrm{t}\mathrm{r}}\left(z\right)$ 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_{\mathrm{S}}$ such that $C_{\mathrm{c}\mathrm{t}\mathrm{r}}\left(H_{\mathrm{S}}\right)=K_{\mathrm{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_{\mathrm{S}}$, we can take the value $H_{\mathrm{S}}$ as the thickness of the vertical growth zone, leading to the vertical ascending speed: ${\displaystyle V_{\mathrm{H}}\propto H_{\mathrm{S}}{\lambda }_{\mathrm{S}}}$.

Figure 5 with 2 supplements see all

Download asset Open asset

Detailed analysis of Equations 1 and 3 in Materials and methods shows that the scale of the stationary profile is set by $\sqrt{D_{+}/{\lambda }_S}$; cf. Appendix A2.3. This is verified in Figure 5B where the stationary profile $C_{\mathrm{c}\mathrm{t}\mathrm{r}}\left(z\right)$ is computed by repeating the simulation for different growth rate ${\lambda }_{\mathrm{S}}$: the profile collapses into the same curve for different values of ${\lambda }_{\mathrm{S}}$ when plotted against $z\sqrt{{\lambda }_{\mathrm{S}}}$; see Figure 5—figure supplement 2 for the same profiles without rescaling in z. Given this scaling, we expect the thickness of growth zone to decrease as ${\displaystyle H_{\mathrm{S}}\propto 1/\sqrt{{\lambda }_{\mathrm{S}}}}$ for faster growth (due to stronger nutrient depletion), leading to a sublinear dependence of the vertical ascending speed, ${\displaystyle V_{\mathrm{H}}\propto \sqrt{{\lambda }_{\mathrm{S}}}}$. Our numerical result on the growth of vertical height is shown as open blue symbols in Figure 5C. The results are well fitted by the square-root dependence on ${\lambda }_{\mathrm{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 ${\displaystyle {\lambda }_{\mathrm{S}}<0.4{\mathrm{h}}^{-1}}$, is difficult to realize by changing carbon sources. However, if we just fit the data in Figure 5C for ${\displaystyle {\lambda }_{\mathrm{S}}>0.5{\mathrm{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.

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\mathrm{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 $\overrightarrow{n}(\overrightarrow{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\mathrm{\mu }\mathrm{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.

Figure 6

Download asset Open asset

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

So, what controls the annulus width? Or, equivalently, what determines the transition of velocity to the vertical orientation in the interior? Qualitatively, the difference between the peripheral and interior regions can be appreciated by looking at the coarse-grained pressure field $P(\overrightarrow{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.

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_{\mathrm{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).

Figure 7

Download asset Open asset

To elucidate the determinants of buckling, we plot in Figure 8A the azimuthal-averaged radial velocity profile $V_r\left(\Delta r\right)$ for the bottom layer of cells, for a range of (signed) distances $\Delta r$ into the edge of the colony; see Appendix 1 Equations (A1.5.1 and A1.5.2) for the definitions of $V_r$ and $\Delta r$. This radial velocity profile, which is stationary for $t\ge 12\mathrm{h}$, is nearly zero in the colony interior, but increases almost linearly within a $~20\mathrm{\mu }\mathrm{m}$ region at the outermost periphery of the colony. Since the radial expansion speed of the colony $V_{\mathrm{R}}$ is simply $V_r$ at $\Delta r=0$, we see that the width of this annulus together with the slope of $V_r\left(\Delta r\right)$ set the radial expansion speed of the entire colony.

Figure 8 with 2 supplements see all

Download asset Open asset

To understand what goes on in this peripheral region, we examine in Figure 8B the azimuthal-averaged height profile of the colony, $h\left(\Delta r\right)$, which is also stationary for $t\ge 12\mathrm{h}$, with $h\left(\Delta r\right)\approx 1\mathrm{\mu }\mathrm{m}$ in the $~20\mathrm{\mu }\mathrm{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 $\overrightarrow{\nabla }\cdot \overrightarrow{V}=\lambda$, where $\lambda$ 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

Throughout the periphery region, the growth rate is essentially the maximal growth rate, that is $\lambda \approx {\lambda }_{\mathrm{S}}$, since the nutrient concentration in this region is set by the boundary value $C_{\mathrm{s}}$ which well exceeds the Monod constant $K_{\mathrm{S}}$; cf. Figure 8—figure supplement 1. Solving the above equation in the annulus in the limit $\left|\Delta r\right|\ll R$ yields a linear dependence,

where we used the definition of the radial expansion speed ${V_{\mathrm{R}}=V}_r\left(\Delta r=0\right)$. 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_{\mathrm{b}}$, defined as the largest value of $\left|\Delta r\right|$ where $V_r\left(\Delta r\right)=0$, sets the radial expansion speed since $V_{\mathrm{R}}\propto {\lambda }_S\cdot W_{\mathrm{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_{\mathrm{R}},$ while the interior cells that are away from the colony edge by this buckling width flow up vertically; cf. Figure 8—figure supplement 2. The magnitude of the buckling width is set by the radial forces exerted on the monolayer of cells. As these cells grow outward, they experience frictions from the agar substrate as well as surface tension that holds them down flat. These two forces lead to the accumulation of pressure, which is built up from the periphery. Figure 8D shows the azimuthally averaged pressure $P\left(\Delta r\right)$ for the bottom layer of cells. At the inner edge of the monolayer region, pressure reaches a critical value that surface tension can no longer hold cells in a single layer. There, some cells buckle into the vertical dimension, leading to vertical flow of cells and forming multiple layers (Figure 8—figure supplement 2), alleviating the further build-up of pressure. It is interesting to observe that this buckling phenomenon is already evident early on during transition from monolayer growth to multiple layers, as shown in Figure 3D. The 20 µm annulus of monolayer at the periphery is set soon after the initial buckling transition, at around $t=8\mathrm{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 ${\mathrm{W}}_{\mathrm{b}}$ at different radial expansion speed ${\mathrm{V}}_{\mathrm{R}}$: Static friction would increase for increased ${\mathrm{V}}_{\mathrm{R}}$, leading to decreased buckling width, whereas dynamic friction would not be affected by the radial expansion speed. Experimentally, we characterized ${\mathrm{V}}_{\mathrm{R}}$ in sugars supporting different growth rates ${\mathrm{\lambda }}_{\mathrm{S}}$, and ${\mathrm{V}}_{\mathrm{R}}$ is seen to increase linearly with ${\mathrm{\lambda }}_{\mathrm{S}}$ (Figure 1I), suggesting a constant ${\mathrm{W}}_{\mathrm{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 ${\mathrm{\lambda }}_{\mathrm{S}}$. The buckling width ${\mathrm{W}}_{\mathrm{b}}$ is indeed not dependent on ${\mathrm{\lambda }}_{\mathrm{S}}$ (Figure 8E), and the radial expansion speed ${\mathrm{V}}_{\mathrm{R}}$ is indeed linear in ${\mathrm{\lambda }}_{\mathrm{S}}$ (open red symbols and dashed red line, Figure 8F). In contract, static friction leads to a much weaker dependence of $V_{\mathrm{R}}$ on ${\lambda }_{\mathrm{S}}$ (blue triangles in Figure 8F).

The linear dependence on ${\lambda }_{\mathrm{S}}$ seen in the experimental data in Figure 1I (red symbols) however exhibits a noticeable horizontal offset. This offset likely results from an additional effect we have not included into the model so far: The size of the cells is dependent on their growth rate, with faster growth rate being longer and wider (Jun and Taheri-Araghi, 2015; Nanninga and Woldringh, 1985; Taheri-Araghi et al., 2015). By repeating the established dependence of cell size on growth rate (see Equation (A2.2.3) in Appendix 2) for different values of ${\lambda }_{\mathrm{S}}$, we recover a nonlinear dependence of $V_{\mathrm{R}}$ on ${\lambda }_{\mathrm{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 ${\displaystyle {\lambda }_{\mathrm{S}}>0.5{\mathrm{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 ${\displaystyle H_{\mathrm{S}}\propto 1/\sqrt{{\lambda }_{\mathrm{S}}}}$ does not depend on ${\mathcal{l}}_{\mathrm{d}\mathrm{i}\mathrm{v}}$ (Appendix 2.3).

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_{\mathrm{R}}$, but not the vertical ascension speed $V_{\mathrm{H}}$. Conversely, we expect changing the nutrient concentration $C_{\mathrm{s}}$ to change the vertical nutrient penetration length and hence $V_{\mathrm{H}}$, but not $V_{\mathrm{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_{\mathrm{s}}$ doubled in Figure 9A (only $V_{\mathrm{H}}$ increased), and with the cell-agar friction quartered in Figure 9B (only $V_{\mathrm{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.

Figure 9

Download asset Open asset

Figure 9—source data 1

Experimental data on the horizontal and vertical colony expansion speeds at various glucose and agar concentrations.

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

Download elife-41093-fig9-data1-v1.xlsx

We note that the actual values of radial and vertical expansion speeds obtained ($V_{\mathrm{R}}=17.2\mu \mathrm{m}/\mathrm{h}$ and $V_{\mathrm{H}}=5.8\mu \mathrm{m}/\mathrm{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_{\mathrm{R}}$ and $V_{\mathrm{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_{\mathrm{R}}$ and $V_{\mathrm{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.

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 $\lambda$-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 NCM3722$\Delta motA$) (Kim et al., 2012) by P1 transduction, yielding strain EQ59, in which the gfp gene is constitutively expressed in the absence of TetR.

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₄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.

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 $O{D}_{600}\approx 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 $O{D}_{600}\approx {10}^{-6}$. $10\mu L$ of culture (containing approximately 10 cells) was spread over pre-warmed plates and transferred immediately to a 37°C incubator for growth. Petri dishes remained covered at all times, except during periodic measurements with a confocal microscope, in order to minimize moisture loss.

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\times /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.

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\times 5$ pixels centered at pixel $(\mathrm{i},\mathrm{j})$, we count the number of pixels whose intensity is above a threshold, and call it $n_{i,j}$. Pixel $(\mathrm{i},\mathrm{j})$ is assigned as type 1 if ${\displaystyle 16>n_{i,j}\ge 3}$, and as type 2 if $n_{i,j}\ge 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_{\mu m/px}\sqrt{N_{px2}/\pi }$ and $r_{outer}=r_{\mu m/px}\sqrt{{(N}_{px1}+N_{px2})/\pi }$, where $N_{px1}$ and $N_{px2}$ are the total numbers of pixels of type 1 and 2, respectively, and $r_{\mu m/px}\approx 0.84$ is the ratio of µm per pixel in our confocal image.

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.

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(\overrightarrow{r},t)$ to represent the nutrient concentration at a spatial location $\overrightarrow{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 ${\displaystyle z<0}$, while cells grow on top of the agar in the region ${\displaystyle z>0}$, and bounded by the colony surface ${\Gamma }_{01}$ to be defined below; see Figure 11. Nutrient diffuses in the two compartments, agar and colony, according to the diffusion equations

Figure 11

Download asset Open asset

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, $\rho =\rho (\overrightarrow{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 $\rho =\rho (\overrightarrow{r},t)$ by a constant value ${\rho }_0$. Above the spatial scale of a few cell lengths, the spatial variation in $\rho$ is ${\displaystyle <5\mathrm{\%}}$ within the colony; see Figure 3—figure supplement 1. The upper boundary of the colony (${\Gamma }_{01}$ in Figure 11) is defined by thresholding the density; see Appendix 1.2. The local mass growth rate $\lambda =\lambda (\overrightarrow{r},t)$ is given by Monod kinetics

where ${\lambda }_{\mathrm{S}}$ is the batch culture growth rate for cells in a medium saturated with some sugar S, and $K_{\mathrm{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 ${\displaystyle \left(z<0\right)}$ and colony side ${\displaystyle \left(z>0\right)}$, 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 ${\Gamma }_{01},{\Gamma }_{02},{\mathrm{a}\mathrm{n}\mathrm{d}\Gamma }_{\mathrm{b}},$ and the Dirichlet boundary condition $C=C_{\mathrm{s}}$ on the lateral wall of agar region ${\Gamma }_{\mathrm{s}}$; cf. Figure 11. The parameter $C_{\mathrm{s}}$ mimics the nutrient concentration far away from the colony. It is one of the key parameters in our study.

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 ${\overrightarrow{r}}_i\left(t\right)$ to denote the center-of-mass of the cell, a unit vector ${\overrightarrow{n}}_i\left(t\right)$ to denote its orientation (the direction along the cylindrical axis), and ${\mathcal{l}}_i\left(t\right)$ to denote the length of the cylinder between the two centers of hemispherical caps. Each cell starts off with the same cylindrical length ${\mathcal{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 $\tilde{\lambda }\left(t\right).$ We call this the cell elongation rate. It is proportional to the mass growth rate $\lambda \left(t\right)$ of the cell with a geometrical proportionality factor $\sigma$: $\tilde{\lambda }=\sigma \lambda ;$ see Appendix 1.3.

Figure 12

Download asset Open asset

The mass growth rate is calculated based on the nutrient concentration at the center ${\overrightarrow{r}}_i\left(t\right)$ of the cell i at time $t$, that is ${\lambda }_i\left(t\right)=\lambda ({\overrightarrow{r}}_i\left(t\right),t)$, according to Equation 3. The growth of cylindrical length ${\mathcal{l}}_i\left(t\right)$ is then given by the growth equation

Once the cylindrical length ${\mathcal{l}}_i\left(t\right)$ reaches a critical value ${\mathcal{l}}_{\mathrm{d}\mathrm{i}\mathrm{v}}$, the cell divides into two daughter cells with cylindrical lengths being ${\mathcal{l}}_0$ with small fluctuation; see Figure 12B and Appendix 1.3. For different growth media supporting different growth rates ${\lambda }_{\mathrm{S}}$, the value of ${\mathcal{l}}_{\mathrm{d}\mathrm{i}\mathrm{v}}$ is in general growth-rate dependent (Jun and Taheri-Araghi, 2015; Taheri-Araghi et al., 2015), the consequences of which are discussed above in 'Radial expansion – quantitative analysis'.

The position and orientation of cell i change according to its velocity ${\overrightarrow{v}}_i$ and angular velocity ${\overrightarrow{\omega }}_i$, which follow Newton’s second law

where $M_i$ and $I_i$ are the mass and moment of inertia of the cell, and ${\overrightarrow{F}}_i^{\mathrm{n}\mathrm{e}\mathrm{t}}$ and ${\overrightarrow{T}}_i^{\mathrm{n}\mathrm{e}\mathrm{t}}$ 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 ${\Gamma }_{01}$ in Figure 11) expands. The nutrient concentration in the new domain requires an update by solving the boundary-value problem of Equations 1–4 again.

The net force ${\overrightarrow{F}}_i^{\mathrm{n}\mathrm{e}\mathrm{t}}$ and torque ${\overrightarrow{T}}_i^{\mathrm{n}\mathrm{e}\mathrm{t}}$ exerted on cell i arise from (a) cell-cell mechanical interaction, (b) cell-agar interaction (if the cell touches the agar surface), (c) cell-fluid interaction, and (d) surface tension (if the cell is on top of the colony); cf. Figure 2. Below we briefly describe each component used in our model. Details are provided in Appendix 1.4.

(a) Cell-cell interaction

Request a detailed protocol

In the interior of a colony, two cells interact only if they are in direct physical contact, characterized here by the overlap ${\delta }_{\mathrm{c}\mathrm{c}}$ in their sphero-cylinder cell boundaries; see Figure 13A. At the point of contact, the cell-cell interaction force ${\overrightarrow{F}}_{\mathrm{c}\mathrm{c}}$ is decomposed into the normal and tangential components, of magnitudes $F_{\mathrm{c}\mathrm{c},n}$ and $F_{\mathrm{c}\mathrm{c},t}$ as defined in and Appendix 1.4a. The normal force includes the Hertzian elasticity force with magnitude ${\displaystyle F_{\mathrm{c}\mathrm{c},\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{s}}\propto \sqrt{w_0}{\delta }_{\mathrm{c}\mathrm{c}}^{3/2}}$ (Hertz, 1882; Johnson, 1985). Additionally, the normal and tangential force each has a dissipation component, of magnitude $F_{\mathrm{c}\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},n}$ and $F_{\mathrm{c}\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},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.

Figure 13

Download asset Open asset

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_{\mathrm{c}\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},n}$ and $F_{\mathrm{c}\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},t}$, respectively, are taken to be proportional to the normal and tangential components of the relative cell-cell velocity, $v_{\text{cc},n}$ and $v_{\text{cc},t}$. We use ${\displaystyle F_{\mathrm{c}\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},n}\propto {\delta }_{\mathrm{c}\mathrm{c}}^{1/2}{v}_{\text{cc},n}}$ and ${\displaystyle F_{c\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},t}\propto {\delta }_{\mathrm{c}\mathrm{c}}{v}_{\text{cc},t}}$, where the additional dependences on the overlap ${\delta }_{\mathrm{c}\mathrm{c}}$ 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_{\mathrm{c}\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},t}^{\text{max}}={\mu }_{\mathrm{c}\mathrm{c}}{F}_{\mathrm{c}\mathrm{c},\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{s}}$, where ${\mu }_{\mathrm{c}\mathrm{c}}$ is the dynamic frictional coefficient; see Appendix 1.4a. Thus, the full cell-cell interaction force is given by

(7a) $F_{\mathrm{c}\mathrm{c},n}=F_{\mathrm{c}\mathrm{c},\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{s}}+F_{\mathrm{c}\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},n},$

(7b) $F_{\mathrm{c}\mathrm{c},t}=\min \left\{F_{c\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},t},{\mu }_{\mathrm{c}\mathrm{c}}{F}_{\mathrm{c}\mathrm{c},\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{s}}\right\}.$

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.

(d) Surface tension

Request a detailed protocol

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.

Figure 14

Download asset Open asset

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_{\text{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,

(8) $F_{\text{surf},0}=\pi {\gamma }_{\text{surf}}{w}_0,$

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\left(x,y,t\right)$ (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\left(x,y,t\right)$ by adding $\delta h$ to the colony height. This thickness $\delta h$ depends primarily on the agar hardness, being larger for softer agar where cells are less tightly bound by the liquid. For each cell whose maximum height $h_{\text{cell}}$ exceeds the liquid height $h_{\mathrm{w}}$, we impose a restoring force normal to the liquid surface; see Figure 14C. As the magnitude of the saturating restoring force $F_{\text{surf},0}$ (Equation 8) is independent of the height difference $\delta z\equiv h_{\text{cell}}-h_{\mathrm{w}}$ for ${\displaystyle \delta z>0}$, the restoring force can be mathematically written as $F_{\text{surf}}=F_{\text{surf},0}\cdot u\left(\delta z\right)$, where $u$ is the Heaviside step function. To avoid numerical instability, we make a linear extrapolation between $F_{\text{surf}}=0$ and $F_{\text{surf}}=F_{\text{surf},0}$ over a narrow transition region ${\displaystyle 0<\delta z<w_0/10}$, which is $1/10$ of the cell width.

Once all the individual forces exerted on a cell i described above are calculated, the net force ${\overrightarrow{F}}_i^{\mathrm{n}\mathrm{e}\mathrm{t}}$ and the corresponding torque ${\overrightarrow{T}}_i^{\mathrm{n}\mathrm{e}\mathrm{t}}$ 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 ${\overrightarrow{r}}_{ji}$ are the corresponding displacement vectors from the points where the forces are exerted to the cell center.

We define coarse-grained fields of cell spatial mass density $\rho (\overrightarrow{r},t)$, velocity $\overrightarrow{v}(\overrightarrow{r},t)$, directors $\overrightarrow{n}(\overrightarrow{r},t)$, and pressure $P(\overrightarrow{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.

We fix the coefficient of nutrient diffusion in the agar to be $D_{-}=600\mu {\mathrm{m}}^2/\mathrm{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\mu {\mathrm{m}}^2/\mathrm{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.5g_{\mathrm{C}\mathrm{D}\mathrm{W}}/g_{\mathrm{g}\mathrm{l}\mathrm{u}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{e}}$ (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{\rho }_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}$, with ${\rho }_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}=0.137\times {10}^{-12}{g}_{\mathrm{C}\mathrm{D}\mathrm{W}}/\mu {\mathrm{m}}^3$ being the cell dry weight density (Basan et al., 2015). The Monod constant is taken to be that for glucose, $K_{\mathrm{S}}=20\mathrm{\mu }\mathrm{M}$ (Monod, 1949). The cell dividing length ${\mathcal{l}}_{\mathrm{d}\mathrm{i}\mathrm{v}}=3\mathrm{\mu }\mathrm{m}$ and the diameter $w_0=1\mathrm{\mu }\mathrm{m}$ are fixed in all the simulations unless otherwise indicated. The constant value $C_{\mathrm{s}}$ of nutrient concentration in the boundary conditions for the diffusion equations and the batch culture growth rate ${\lambda }_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.

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 $\rho =\rho (\overrightarrow{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 $\Delta t$. Usually, we run through one main loop per 100—1,000 inner loops. In updating the nutrient concentration, we use the finite difference to discretize the equations and the Jacobi or Gauss-Seidel relaxation method to solve the resulting systems of linear equations. We use multi-resolution adaptive grids for a large computational domain, and use the OpenMP for parallelizing our code. See Appendix 1 for details. On a multi-processor (14-16 processors) computer, the simulation can reach a colony of a few million cells in 24 hours. We have placed the major and basic parts of our C++ codes in the repository GitHub (Warren et al., 2019; copy archived at https://github.com/elifesciences-publications/CellsMD3D).

Our computational box is $\mathrm{\Omega }=(-L,L)\times (-L,L)\times (-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 ${\mathrm{\Omega }}_0$, colony region ${\mathrm{\Omega }}_1$, and agar region ${\mathrm{\Omega }}_2=(-L,L)\times (-L,L)\times (-a,0)$, respectively. The colony surface or colony-air interface ${\mathrm{\Gamma }}_{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 ${\mathrm{\Gamma }}_{12}$. The remaining part, denoted ${\mathrm{\Gamma }}_{02}$, separates the air from agar. Note that, since the bacterial colony grows with time $t$, all the air region ${\mathrm{\Omega }}_0$, the colony region ${\mathrm{\Omega }}_1$, the colony-air interface ${\mathrm{\Gamma }}_{01}$, and the colony-agar interface ${\mathrm{\Gamma }}_{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 ${\mathrm{\Omega }}_1$ and the agar region ${\mathrm{\Omega }}_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_{\mathrm{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_{\mathrm{ran}}]$ with $h_{\mathrm{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 $\mathrm{\Delta }t$, that consists of the following steps:

(3.1) Simulate the cell growth and division;

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

(3.3) Update all the cell positions;

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

(3.5) Set $t:=t+\mathrm{\Delta }t$ and continue with Step (3.1).

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

We model an underlying E. coli bacterial cell as a sphero-cylinder; cf. Figure 12A in the main text. We denote by $\overrightarrow{p}$ and $\overrightarrow{q}$ the centers of the two hemispheres. We call $\mathrm{\ell }=\parallel \overrightarrow{p}-\overrightarrow{q}\parallel$ the cell cylindrical length which can vary with time $t$. We also denote by $\overrightarrow{n}=(\overrightarrow{q}-\overrightarrow{p})/\mathrm{\ell }$ the unit vector pointing from one center of hemisphere $\overrightarrow{p}$ to the other $\overrightarrow{q}$, and call it the direction of the cell. Note that the center of mass of the cell is ${\overrightarrow{r}}_{\mathrm{c}}=(\overrightarrow{p}+\overrightarrow{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 ${\rho }_{\mathrm{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.

Mechanical forces exerted on a cell include the elastic and dissipative forces from the cell-cell mechanical interactions, the elastic and dissipative forces from the cell-agar interaction if the cell is in contact with the agar, the surface tension force if the cell is on the top of the colony, and the viscous force from the interaction between the cell and the surrounding liquid. For a given cell, we shall denote these forces by ${\overrightarrow{F}}_{\mathrm{cc}}$, ${\overrightarrow{F}}_{\mathrm{ca}}$, ${\overrightarrow{F}}_{\mathrm{surf}}$, and ${\overrightarrow{F}}_{\mathrm{visc}}$, respectively. So, the total force acting on the cell is

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