Spatiotemporal establishment of dense bacterial colonies growing on hard agar
Abstract
The physical interactions of growing bacterial cells with each other and with their surroundings significantly affect the structure and dynamics of biofilms. Here a 3D agentbased 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 nonmotile, non EPSproducing E. coli cells on hard agar. The model, supported by experiment on colony growth in different types and concentrations of nutrients, suggests that radial colony expansion is not limited by nutrients as commonly believed, but by mechanical forces. Nutrient penetration instead governs vertical colony growth, through thin layers of vertically oriented cells lifting up their ancestors from the bottom. Overall, the model provides a versatile platform to investigate the influences of metabolic and environmental factors on the growth and morphology of bacterial colonies.
https://doi.org/10.7554/eLife.41093.001Introduction
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, biofilmbound bacteria alter their microenvironment 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 threedimensional computational model, incorporating the elements of nutrient diffusion, cellcell and cellagar mechanical interactions, and introducing a unique celllevel 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 cellagar 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 nonmotile E. coli. A set of minimum media with various carbon sources was used to vary the cell growth rate.
Results
Experimental results
Experiments were performed using E. coli K12 strain EQ59, which is nonmotile and harbors constitutive GFP expression; see 'Experimental Methods'. Each colony was inoculated as a single cell from batch culture growing in midlog 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\mathrm{h}$ at a radius of ~$75\text{}\mu \mathrm{m}$ (Figure 1C–E and F). It then developed into a 3D colony over time, maintaining an approximate conic shape through the ensuing 1015 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}_{\mathrm{R}}\approx 45.2\mathrm{\mu}\mathrm{m}/\mathrm{h}$ and the colony height increased also linearly at a vertical speed ${V}_{\mathrm{H}}\approx 12.4\mathrm{\mu}\mathrm{m}/\mathrm{h}$ (Figure 1H), reaching a radius of $\sim 500\mathrm{\mu}\mathrm{m}$ and a height of $\sim 150\mathrm{\mu}\mathrm{m}$ by $t=24\mathrm{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\le t\le 24\mathrm{h}$, where both the radial and vertical growth are linear.

Figure 1—source data 1
 https://doi.org/10.7554/eLife.41093.005
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{\mathrm{h}}^{1}$ to $1{\mathrm{h}}^{1}$; see Supplementary file 1Table 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, threedimensional model consisting of the diffusion of nutrient through the agar and the colony; the growth, division, and movement of individual cells; and the cellcell, cellagar, cellsurface 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 1Tables S2S4); 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 ${\lambda}_{\mathrm{S}}=1.0{\mathrm{h}}^{1}$, which corresponds approximately to the growth of E. coli in glucose minimal medium (Supplementary file 1Table 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 crosssectional 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 singlelayer 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 noncompact 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'.
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\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.
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 coarsegrain 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, coarsegrained 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 coarsegrained 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 $z\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{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 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 µ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 crosssectional 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.
Vertical rise – quantitative analysis
This discshaped 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$ = 1224 hours in Figure 1): As the colony has the shape of a cone (Figures 1G and 3B), its volume is given by $V}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{y}}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}{R}^{2}H\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}{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 $\frac{d}{dt}{V}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{y}}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}{V}_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}\left(t\right)$ leads to $V}_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}{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: $V}_{\mathrm{H}}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}{H}_{\mathrm{S}}{\lambda}_{\mathrm{S}$.
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 $H}_{\mathrm{S}}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}1/\sqrt{{\lambda}_{\mathrm{S}}$ for faster growth (due to stronger nutrient depletion), leading to a sublinear dependence of the vertical ascending speed, $V}_{\mathrm{H}}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}\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 squareroot 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 squareroot relation, for $\lambda}_{\mathrm{S}}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0.4\text{}{\mathrm{h}}^{1$, is difficult to realize by changing carbon sources. However, if we just fit the data in Figure 5C for $\lambda}_{\mathrm{S}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0.5\text{}{\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 2fold 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\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 crosssectional view shown in Figure 4A. We again coarsegrain 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 inplane 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 inplane component becomes smaller in magnitude. Interior to the annulus, the inplane projection of the director field vanishes, confirming that they are largely oriented vertically.
Next, we examine the coarsegrained velocity field $\overrightarrow{v}(\overrightarrow{r},t)$ for the bottom layer of cells shown in Figure 6A, with the xy projection of $\overrightarrow{v}(\overrightarrow{r},t)$ shown as arrows in Figure 6C. We observe a narrow annulus of nonvanishing velocity field (arrows with finite length) at the outermost edge pointing radially outward; see also the side view provided in Figure 4D. Since the inplane 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 coarsegrained 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.
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}_{\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).
To elucidate the determinants of buckling, we plot in Figure 8A the azimuthalaveraged 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.
To understand what goes on in this peripheral region, we examine in Figure 8B the azimuthalaveraged 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 buildup 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 cellagar friction. Two types of friction have been used in the cellmodeling literature, one which depends linearly on the cellagar 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 TaheriAraghi, 2015; Nanninga and Woldringh, 1985; TaheriAraghi 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 $\lambda}_{\mathrm{S}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0.5\text{}{\mathrm{h}}^{1$ (dotted red line). On the other hand, the growthrate dependence of cell sizes has no noticeable effect on the vertical ascension speed (filled blue symbols, Figure 5C) since the growth zone thickness $H}_{\mathrm{S}}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}1/\sqrt{{\lambda}_{\mathrm{S}}$ does not depend on ${\mathcal{l}}_{\mathrm{d}\mathrm{i}\mathrm{v}}$ (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 cellagar 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 cellagar 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 cellagar 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—source data 1
 https://doi.org/10.7554/eLife.41093.021
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 1Tables 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 cellagar 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.
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 nonmotile 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 multiresolution numerical technique. With a discrete agentbased model, we captured mechanical interactions, including elasticity and dynamic friction. Most importantly, the surface tension of the liquid in the colony is implemented by introducing a restoring force on cells protruding from a smoothened colony surface.
Our model is able to capture quantitatively some of the characteristic features observed for bacterial colony growth, including the conic shape of the colony, the linear expansion of colony radius and height, and both the linear and sublinear dependence of the speed of radial expansion and that of vertical expansion, respectively, on the cell growth rate. The model makes a number of important predictions on the expanding colony as summarized in Figure 10: The growth zone is predicted to be disclike 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 cellagar 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 celllevel surface tension, which we believe are at the root of all behaviors described above, including the forms of radial and vertical colony growth. A key result of our study is that the linear radial growth is driven by the growth of a thin layer of radially oriented cells located at the colony periphery, whose width is determined by mechanical buckling. Although the linear radial expansion of bacterial colonies has been known for about 50 years (Pirt, 1967), for a long time this was attributed to a ringshaped 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 ‘gravitylike’ 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 celllevel 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 singlelayer cell growth, through buckling, to the growth of a macroscopic colony. This celllevel 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 rateindependent 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 growthrate dependence of radial expansion speed was obtained (Figure 8F blue triangles). Along a different line, FisherKolmogorov (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 squareroot 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 wellknown 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 (EnosBerlage 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, cellagar 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 inplane compression force and cellagar adhesion. Due to the different energy barriers against verticalization, the length scales of verticalization between our model and that of (Beroz et al., 2018) are very different: The colonies in Beroz et al. spread very slowly radially ($~3\mu \mathrm{m}/\mathrm{h}$), and verticalization occurs at a colony radius of $~510\mu \mathrm{m}$. Colonies in our model spread much faster ($~14\mu \mathrm{m}/\mathrm{h}$), and substantial verticalization occurs at a radius of $~250\mu \mathrm{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 airliquid 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\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}24\text{}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 colonybased 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 cellagar 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 cellagar 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; EnosBerlage and McCarter, 2000; Yan et al., 2016). While the precise nature of the forces driving verticalization may be different in the two cases, the underlying origins may be similar — mechanical instability due to inplane compression resulting from colony expansion and cellagar 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 largescale data being generated by highthroughput 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 coarsegrained hydrodynamic models that incorporate the unique features of surface tension and dynamic friction discussed here, and capture the radial and vertical colony growth characteristics, both their temporal behaviors and their dependences on cell growth rates and other environmental factors.
Materials and methods
Experimental methods
Bacterial strain
Request a detailed protocolThe 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 pKD13rrnBT:Ptet (Klumpp et al., 2009), yielding the plasmid pKDT_Ptetgfp. The fragment 'km^{r}:rrnBT:Ptetgfp' present in pKDT_Ptetgfp 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:Ptetgfp' 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:Ptetgfp' 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.
Growth medium
Request a detailed protocolPhosphatebuffered 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 1Table 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
Request a detailed protocolBatch culture growth was performed in a 37°C waterbath 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 precultures. For batch culture growth rate measurements, overnight precultures 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 10fold 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 preculture (prepared as above) was diluted to $O{D}_{600}\approx {10}^{6}$. $10\mu 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
Request a detailed protocolColonies 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 xymontage was created and stitched together to form a single 3D image using the ImageJ Grid/Collection Stitching plugin.
Image analysis
Request a detailed protocolThe 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 multilayers 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 $16\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{n}_{i,j}\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}3$, and as type 2 if ${n}_{i,j}\ge 16$, indicating the pixel belonging to single or multilayer 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 curves
Request a detailed protocolColony 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
Request a detailed protocolWe 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 $z\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0$, while cells grow on top of the agar in the region $z\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}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
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 righthand 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 $\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}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 $\left(z\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0\right)$ and colony side $\left(z\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}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 fluxfree 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.
Discrete model for cell growth, division and movement
Request a detailed protocolIn addition to modeling the nutrient as a continuum, we use a discrete, agentbased 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 agentbased model, each E. coli cell is represented by a spherocylinder, 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 centerofmass 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 $\stackrel{~}{\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 $: $\stackrel{~}{\lambda}=\sigma \lambda ;$ see Appendix 1.3.
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 growthrate dependent (Jun and TaheriAraghi, 2015; TaheriAraghi 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 boundaryvalue problem of Equations 1–4 again.
Discrete model for interaction forces
Request a detailed protocolThe 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) cellcell mechanical interaction, (b) cellagar interaction (if the cell touches the agar surface), (c) cellfluid 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) Cellcell interaction
Request a detailed protocolIn 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 spherocylinder cell boundaries; see Figure 13A. At the point of contact, the cellcell 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 $F}_{\mathrm{c}\mathrm{c},\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{s}}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}\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.
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 cellcell velocity, ${v}_{\text{cc},n}$ and ${v}_{\text{cc},t}$. We use $F}_{\mathrm{c}\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},n}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}{\delta}_{\mathrm{c}\mathrm{c}}^{1/2}{v}_{\text{cc},n$ and $F}_{c\mathrm{c},\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p},t}\phantom{\rule{thinmathspace}{0ex}}\propto \phantom{\rule{thinmathspace}{0ex}}{\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 cellcell interaction force is given by
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) Cellagar interaction
Request a detailed protocolThe force exerted on a cell in contact with the agar, ${\overrightarrow{F}}_{\mathrm{c}\mathrm{a}}$, 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}_{\mathrm{r}\mathrm{a}\mathrm{n}}$ which is the maximum fluctuation of the agar surface around its mean $(z=0)$.
(c) Cellfluid interaction
Request a detailed protocolA cell also interacts with the surrounding fluid and experiences viscous drag. This is given by the Stokes drag force ${\overrightarrow{F}}_{\text{visc}}$, which is proportional to the velocity of the cell $\overrightarrow{v}$. We note that in highdensity colonies such as those studied here, dissipation due to viscous drag is significantly less than the cellular friction force.
(d) Surface tension
Request a detailed protocolSurface 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 (BenJacob et al., 2000). In both cases, surface tension reflects the curvature of the macroscopic colony profile. However, such coarsegrained treatments of surface tension cannot describe the initial layerbylayer 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 liquidlike.
In the same way, the cohesion of a threedimensional colony is maintained by the interaction of cells with the surrounding liquid at the airliquid boundary. As shown in Figure 14A, the red curve indicates a smooth airliquid 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 airliquid 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 celllevel 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 singlecell level in our model, we first compute the coarsegrained 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 $\delta z\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}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 $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\delta z\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{w}_{0}/10$, which is $1/10$ of the cell width.
Pressure calculation
Request a detailed protocolOnce 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.
Coarsegrained variables
Request a detailed protocolWe define coarsegrained 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 finitedifference grid boxes); see Appendix 1.5 for details.
Model parameters
Request a detailed protocolWe 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 1Table 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}_{\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 1Tables S2–S4 for the definitions and estimated values of all the parameters.
Numerical implementation
Request a detailed protocolWe 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 smalltime 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 GaussSeidel relaxation method to solve the resulting systems of linear equations. We use multiresolution adaptive grids for a large computational domain, and use the OpenMP for parallelizing our code. See Appendix 1 for details. On a multiprocessor (1416 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/elifesciencespublications/CellsMD3D).
Appendix 1
Simulation Model and Methods
We model the colony expansion through the coupling of the growth, division, and movement of individual cells with the diffusion and reaction of nutrients and wastes. The growth of an individual cell within a short time period is described by a linear growth equation. The local growth rate varies spatially and temporally, and is determined locally by the cell density and nutrient supply. Once a cell grows into a critical size, it divides itself into two daughter cells with some randomness in their cylindrical lengths and orientations. In the meantime, growing and dividing cells push each other, generating mechanical forces. These forces, together with those arising from the cellagar interactions, cellliquid interactions, and surface tension, determine the movement of individual cells described by Newton’s law of motion. At any instant of time, the coarsegraining 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 $\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 colonyair 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 colonyair interface ${\mathrm{\Gamma}}_{01}$, and the colonyagar 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 reactiondiffusion 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.
A1.2 Nutrient and Growth Rate Update
Cell density and colony boundary
Given the positions of all the cells at time $t$, we define the (local) volume fraction $\varphi =\varphi (\overrightarrow{r},t)(\overrightarrow{r}=(x,y,z))$ of the cells in each finite difference grid box above $z=0$ by
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
where ${\rho}_{\mathrm{cell}}$ is the constant mass density of a typical mature cell, and can be estimated from experiment; cf. Supplementary file 1Table S3. We also define the colony region ${\mathrm{\Omega}}_{1}={\mathrm{\Omega}}_{1}(t)$ to be that with $\varphi (\overrightarrow{r},t)\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$; cf. Figure 11 in the main text.
Nutrient update: Reactiondiffusion equations and boundary conditions
The concentration field $C=C(\overrightarrow{r},t)$ of the nutrient is defined spatially on the colony and agar regions ${\mathrm{\Omega}}_{1}={\mathrm{\Omega}}_{1}(t)$ and ${\mathrm{\Omega}}_{2}$, respectively; cf. Figure 11 in the main text. It is governed by the following system of reactiondiffusion 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 ${\mathrm{\Omega}}_{1}(t)$ and the diffusion equation for the concentration in the agar region ${\mathrm{\Omega}}_{2}$, respectively, where ${D}_{}$ and ${D}_{+}$ are the corresponding diffusion coefficients, $Y$ is the yield factor, ${\lambda}_{\mathrm{S}}$ is the batch culture growth rate, ${K}_{\mathrm{S}}$ is the capacity constant (the Monod constant) for sugar, and $\rho =\rho (\overrightarrow{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 ${\rho}_{0}$ to approximate $\rho =\rho (\overrightarrow{r},t)$. We take this constant ${\rho}_{0}$ to be the cell dry weight (CDW) per unit volume of the colony; cf. Supplementary file 1Table S3. Equation (A1.2.4) and Equation (A1.2.5) are the interface conditions for the concentration across the colonyagar interface ${\mathrm{\Gamma}}_{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 agarair interface ${\mathrm{\Gamma}}_{01}$, the colonyair interface ${\mathrm{\Gamma}}_{02}$, and the bottom of agar ${\mathrm{\Gamma}}_{\mathrm{b}}$, we impose the fluxfree boundary condition, with $\partial /\partial n$ denoting the derivative in the normal direction along the corresponding part of the boundary, pointing from the colony or agar to the air region or pointing downward from the bottom of agar. On $\mathrm{\Gamma}}_{\mathrm{s}$, the lateral faces of the agar region, we prescribe a constant value ${C}_{\mathrm{s}}$ of nutrient concentration that represents the maximum nutrient that the system supplies.
In each time step, we update the nutrient by solving the steadystate reactiondiffusion equations with the corresponding boundary conditions, that is, Equation (A1.2.2)–(A1.2.7) with $\partial C/\partial t$ set to be 0. We use an iterative scheme to solve the equations with the previous nutrient concentration as the initial solution. This iterative scheme is constructed based on solving the corresponding timedependent 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 colonyair 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 ${\mathrm{\Omega}}_{2}$ that represents the agar region, we use a nested, multilevel, 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 ${\mathrm{\Gamma}}_{\mathrm{b}}$, and then from grids on the bottom ${\mathrm{\Gamma}}_{\mathrm{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 $er{r}_{\mathrm{conv}}$; cf. Supplementary file 1Table 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(\overrightarrow{r},t)$ is updated, we can define the spatially and temporally varying local mass growth rate $\lambda =\lambda (\overrightarrow{r},t)$ by
This is first defined at each grid point $\overrightarrow{r}$ and then at any other point in the colony region by interpolation.
A1.3 Cell Growth, Division, and Movement
We model an underlying E. coli bacterial cell as a spherocylinder; 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.
Cell growth
With our assumption, a cell grows only in its cylindrical length $\mathrm{\ell}=\mathrm{\ell}(t)$ but not in its diameter of hemispheres. The growth of the cell is governed by the growth equation for the cell cylindrical length
Here, $\stackrel{~}{\lambda}({\overrightarrow{r}}_{\mathrm{c}},t)$ is the cell elongation rate evaluated at the cell center ${\overrightarrow{r}}_{c}$. The elongation rate $\stackrel{~}{\lambda}$ is proportional to the mass growth rate $\lambda $ defined in Equation (A1.2.8): $\stackrel{~}{\lambda}(\overrightarrow{r},t)=\sigma \lambda (\overrightarrow{r},t)$. This effective relation between the two growth rates will be discussed after we describe the cell division. Numerically, we use the firstorder approximation to obtain the cell cylindrical length at time $t+\mathrm{\Delta}t$ by
where $\mathrm{\Delta}t$ is the simulation time step. Initially at $t=0$, all the cells start with a constant cylindrical length ${\mathrm{\ell}}_{0}$.
Cell division
When a cell of centers of hemispheres $\overrightarrow{p}$ and $\overrightarrow{q}$ grows long enough, with its cylindrical length $\mathrm{\ell}\ge {\mathrm{\ell}}_{\mathrm{div}}$ for some critical value ${\mathrm{\ell}}_{\mathrm{div}}>{w}_{0}$, it divides into two daughter cells of cylindrical lengths ${\mathrm{\ell}}_{1}$ and ${\mathrm{\ell}}_{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 ${\mathrm{\ell}}_{1}$ and ${\mathrm{\ell}}_{2}$ of these daughter cells are given by
where $\eta \in (0.5,0.5)$ is a random variable and ${\mathrm{\ell}}_{\text{ran}}$ is an input positive number representing the maximum fluctuation of cell cylindrical length during the cell division; cf. Supplementary file 1Table S3. In average, all new born cells have the same cylindrical length ${\mathrm{\ell}}_{0}$; and we take it to be the same for all the initial cells. This implies that ${\mathrm{\ell}}_{\mathrm{div}}=2{\mathrm{\ell}}_{0}+{w}_{0}$; cf. Figure 12B in the main text. Note that before division, the mother cell has volume $\pi {w}_{0}^{2}{\mathrm{\ell}}_{\mathrm{div}}/4+\pi {w}_{0}^{3}/6$; after division, the total volume of two daughter cells is $\pi {w}_{0}^{2}{\mathrm{\ell}}_{\mathrm{div}}/4+\pi {w}_{0}^{3}/12$. There is $\pi {w}_{0}^{3}/12$ volume loss every time due to division. For a fixed cell aspect ratio ${\mathrm{\ell}}_{\mathrm{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 ${\mathrm{\ell}}_{\mathrm{div}}$ independent of the batch culture growth rate ${\lambda}_{\mathrm{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)
Moreover, these daughter cells inherit the velocity from their mother cell. But the angular velocities of these two daughter cells are set to be
respectively, where $\xi \in (0.5,0.5)$ is a random variable and ${\omega}_{\mathrm{ran}}$ is an input positive number that is the maximum fluctuation of the angular velocity; cf. Supplementary file 1Table S3.
In the cell division, each of the two daughter cells also experiences the rotational fluctuation. Let us fix the center of mass ${\overrightarrow{r}}_{\mathrm{c}}$ of a daughter cell. We rotate the vector $\overrightarrow{p}{\overrightarrow{r}}_{\mathrm{c}}$ with $\overrightarrow{p}$ the center of a hemisphere of this daughter cell. (For notational simplicity, we use ${\overrightarrow{r}}_{\mathrm{c}}$ and $\overrightarrow{p}$ for this daughter cell; and they are different from those of the mother cell.) To do so, we construct a local Cartesian coordinate system with the origin at ${\overrightarrow{r}}_{\mathrm{c}}$ and the three unit coordinate vectors ${\overrightarrow{e}}_{1}$, ${\overrightarrow{e}}_{2}$, and ${\overrightarrow{e}}_{3}$ by $\overrightarrow{e}}_{3}=\left(\overrightarrow{p}{\overrightarrow{r}}_{c}\right)/\overrightarrow{p}{\overrightarrow{r}}_{c$,
and ${\overrightarrow{e}}_{2}={\overrightarrow{e}}_{3}\times {\overrightarrow{e}}_{1}$, where ${\overrightarrow{e}}_{x},{\overrightarrow{e}}_{y},{\overrightarrow{e}}_{z}$ are the unit coordinate vectors in the original coordinate system. Then we rotate the vector $\overrightarrow{p}{\overrightarrow{r}}_{\mathrm{c}}$ by an angle of $\phi $ around the axis ${\overrightarrow{e}}_{2}$, and $\theta $ around ${\overrightarrow{e}}_{3}$ successively. Here, $\phi $ is a random variable in the range $[0,{\phi}_{\mathrm{ran}}\pi ]$, where ${\phi}_{\mathrm{ran}}$ is a constant that sets the magnitude of fluctuations of $\phi $ (cf. Supplementary file 1Table S3) and $\theta \in [0,2\pi )$ is a random number. If we denote by ${\overrightarrow{p}}_{\mathrm{new}}$ the new center of hemisphere after the rotation, then
This equation allows us to find the coordinates of ${\overrightarrow{p}}_{\mathrm{new}}$ in the original coordinate system. The coordinates of the other center of hemisphere ${\overrightarrow{q}}_{\mathrm{new}}$ are given by ${\overrightarrow{q}}_{\mathrm{new}}=2{\overrightarrow{r}}_{\mathrm{c}}{\overrightarrow{p}}_{\mathrm{new}}.$
The two growth rates
We assume that, in the life span of a cell, the cell mass growth rate $\lambda $ and the cell elongation rate $\stackrel{~}{\lambda}$ stay constant. The growth of mass satisfies the equation $M(t)=M(0){e}^{\lambda t}$. If the mass doubling time is ${t}_{\mathrm{d}}$, then $M({t}_{\mathrm{d}})=2M(0)$ and hence $\lambda {t}_{\mathrm{d}}=\mathrm{ln}2.$ On the other hand, the doubling time is the time that a new born cell, which in average has the cylindrical length ${\mathrm{\ell}}_{0}$, grows as its cylindrical length reaches the dividing length ${\mathrm{\ell}}_{\mathrm{div}}$. Note that ${\mathrm{\ell}}_{\mathrm{div}}=2{\mathrm{\ell}}_{0}+{w}_{0}$. From the growth equation Equation (A1.3.2), we have then ${\mathrm{\ell}}_{\mathrm{div}}={l}_{0}{e}^{\stackrel{~}{\lambda}{t}_{\mathrm{d}}}.$ Therefore,
Finally, $\stackrel{~}{\lambda}=\sigma \lambda $, where
For the fixed dividing cell aspect ratio ${\mathrm{\ell}}_{\mathrm{div}}:{w}_{0}=3:1$, we have $\sigma =\mathrm{ln}3:\mathrm{ln}2$.
Cell movement
Consider a cell at some time instant. Let us denote its centers of hemispheres by ${\overrightarrow{p}}_{\mathrm{old}}$ and ${\overrightarrow{q}}_{\mathrm{old}}$, its cylindrical length by ${\mathrm{\ell}}_{\mathrm{old}}=\parallel {\overrightarrow{p}}_{\mathrm{old}}{\overrightarrow{q}}_{\mathrm{old}}\parallel $, its direction by ${\overrightarrow{n}}_{\mathrm{old}}=({\overrightarrow{q}}_{\mathrm{old}}{\overrightarrow{p}}_{\mathrm{old}})/{\mathrm{\ell}}_{\mathrm{old}},$ and its mass by ${M}_{\mathrm{old}}.$ Let us also denote by ${\overrightarrow{v}}_{\mathrm{old}}$ and ${\overrightarrow{\omega}}_{\mathrm{old}}$ the velocity and angular velocity, respectively, at the center of the cell. We apply the velocityVerlet algorithm to update the cell position, velocity, and angular velocity with the simulation time step $\mathrm{\Delta}t.$
We first calculate the force ${\overrightarrow{F}}_{\mathrm{half}}$ and torque ${\overrightarrow{T}}_{\mathrm{half}}$ of the cell. Details of such calculations are given below in section A1.4. We then calculate the velocity ${\overrightarrow{v}}_{\mathrm{half}}$ and angular velocity ${\overrightarrow{\omega}}_{\mathrm{half}}$ for a half time step:
where ${I}_{\mathrm{old},\mathrm{n}}$ and ${I}_{\mathrm{old},\mathrm{t}}$ are the moments of inertia of the cell along the directions ${\overrightarrow{n}}_{\mathrm{old}}$ and its orthogonal, respectively. By direct calculations and using the constant density ${\rho}_{\mathrm{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 $\mathrm{\Delta}t$ by updating the centers of hemispheres
We update the force and torque of the new cell with the centers of hemispheres ${\overrightarrow{p}}_{\mathrm{new}}$ and ${\overrightarrow{q}}_{\mathrm{new}}$ by the procedure of force calculations described below in section A1.4 to get ${\overrightarrow{F}}_{\mathrm{new}}$ and ${\overrightarrow{T}}_{\mathrm{new}}$. Finally, we update the velocity and angular velocity for the second half time step to get
where ${M}_{\mathrm{new}}$, ${\overrightarrow{n}}_{\mathrm{new}}=({\overrightarrow{q}}_{\mathrm{new}}{\overrightarrow{p}}_{\mathrm{new}})/{\mathrm{\ell}}_{\mathrm{new}}$, and ${\mathrm{\ell}}_{\mathrm{new}}=\parallel {\overrightarrow{p}}_{\mathrm{new}}{\overrightarrow{q}}_{\mathrm{new}}\parallel $ are the mass, direction, and cylindrical length, respectively, of the new cell with the centers of hemispheres ${\overrightarrow{p}}_{\mathrm{new}}$ and ${\overrightarrow{q}}_{\mathrm{new}}$, and the moments of inertia ${I}_{\mathrm{new},n}$ and ${I}_{\mathrm{new},t}$ can be calculated similarly using the cylindrical length ${\mathrm{\ell}}_{\mathrm{new}}.$
A1.4 Force Calculations
Mechanical forces exerted on a cell include the elastic and dissipative forces from the cellcell mechanical interactions, the elastic and dissipative forces from the cellagar 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 cellcell and cellliquid interactions.
(a) Cellcell interaction force
When two cells are in direct contact, they generate the cellcell interaction force. As in Volfson et al. (2008), we describe such forces using a standard model in granular solids (Brilliantov et al., 1996; Cundall and Strack, 1979; Kuwabara and Kono, 1987; Makse et al., 2004; Shäfer et al., 1996), with some adjustment based on experimental considerations on bacterial cells. An important part of our force scheme is a detailed treatment of the cellcell frictional force, which together with the cellagar friction are responsible for crucial mechanical behaviors such as buckling of the bacterial colony. While a similar graingrain friction is commonly included in models of granular solids, however, friction in cellular models is often described to only include a purely viscous force (Farrell et al., 2013; Ghosh et al., 2015).
Let us fix a cell in the colony centered at ${\overrightarrow{r}}_{\mathrm{c}}$, and call it a primary cell. Let us also fix a neighboring cell centered at ${\overrightarrow{r}}_{{\mathrm{c}}^{\prime}}.$ We denote by $d$ the minimal distance between the central cylindrical line segments of the two cells, and by $\overrightarrow{a}$ and ${\overrightarrow{a}}^{\prime}$ the points on these line segments, respectively, that reach this minimal distance:
cf.Figure 13A in the main text. We will describe an algorithm of finding the minimum distance $d$ and these two points $\overrightarrow{a}$ and ${\overrightarrow{a}}^{\prime}$ at the end of this part. We denote by ${\overrightarrow{n}}_{\mathrm{cc}}$ the unit vector along the direction from ${\overrightarrow{a}}^{\prime}$ to $\overrightarrow{a},$ by ${\overrightarrow{r}}_{\mathrm{cc}}$ the center of these two points, and by ${\delta}_{\mathrm{cc}}$ the indentation size (i.e., the amount of overlap) between the two cells, respectively:
Let $\overrightarrow{v}$, ${\overrightarrow{v}}^{\prime}$ and $\overrightarrow{\omega}$, ${\overrightarrow{\omega}}^{\prime}$ be the velocities and angular velocities of the primary and neighboring cells, respectively. The velocities $\overrightarrow{V}$ and ${\overrightarrow{V}}^{\prime}$ of the two cells at the midpoint ${\overrightarrow{r}}_{\mathrm{cc}}$ are then given by
respectively. Denote the difference of these velocities by
We call the unit vector ${\overrightarrow{n}}_{\mathrm{cc}}$ the normal direction in the interaction of these two cells. We specify the tangential direction to be the unit vector
if the denominator is nonzero, that is the relative velocity ${\overrightarrow{v}}_{\mathrm{cc}}$ has a nonzero tangential component. (Otherwise, ${\overrightarrow{\tau}}_{\mathrm{cc}}$ can be any unit vector orthogonal to ${\overrightarrow{n}}_{\mathrm{cc}}$.) Note that the tangential direction ${\overrightarrow{\tau}}_{\mathrm{cc}}$ depends on the direction of the relative velocity ${\overrightarrow{v}}_{\mathrm{cc}}.$
Let us assume now that these two cells are in direct contact with each other, that is ${\delta}_{\mathrm{c}\mathrm{c}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0.$ The total interaction force ${\overrightarrow{F}}_{\mathrm{cc}}$ between these two cells exerted at the center ${\overrightarrow{r}}_{\mathrm{cc}}$ and the corresponding torque ${\overrightarrow{T}}_{\mathrm{cc}}$ about the axis ${\overrightarrow{r}}_{\mathrm{cc}}{\overrightarrow{r}}_{\mathrm{c}}$, are given, respectively, by
Here, ${\overrightarrow{F}}_{\mathrm{cc},\mathrm{n}}$ is the normal force in the direction ${\overrightarrow{n}}_{\mathrm{cc}}$ and ${\overrightarrow{F}}_{\mathrm{cc},\mathrm{t}}$ the tangential force in the direction ${\overrightarrow{\tau}}_{\mathrm{cc}}.$ They are defined by
Note that the final form of the force ${\overrightarrow{F}}_{\mathrm{cc}}$ is the same as that in the main text; cf. Equations 7a and b there.
The first part in the normal force ${\overrightarrow{F}}_{\mathrm{cc},n}$ is the Hertz contact force resulting from the elastic collision of the two cells and pointing from the neighboring cell to the primary cell. Here, we approximate the cells as spheres of radius ${w}_{0}/2$. In such a case, the classical Hertz contact force is $(4/3){k}_{\mathrm{cc}}\sqrt{{w}_{0}/4}{\delta}_{\mathrm{cc}}^{3/2},$ which is exactly what we have, where ${k}_{\mathrm{cc}}$ is the reduced (or effective) elastic constant and ${w}_{0}/2$ is the reduced radius. If $E$ and $\nu $ are Young’s module and Poisson’s ratio of cells, then ${k}_{\mathrm{c}\mathrm{c}}=E/(2(1{\nu}^{2}))$ (cf. Johnson, 1985; Landau and Lifshitz, 1986; Popov, 2010). The second part in the normal force ${\overrightarrow{F}}_{\mathrm{cc},n}$ is the friction force in the normal direction ${\overrightarrow{n}}_{\mathrm{cc}}$ due to the relative motion of the two cells, where ${\gamma}_{\mathrm{cc},n}$ is the (normal) static friction coefficient and
is the reduced mass.
This form of the normal force Equation (A1.4.1) has been used in Volfson et al. (2008). Similar forms of such a normal force can be found in the literature of granular solids (cf. Brilliantov et al., 1996; Herrmann and Luding, 1998; Kuwabara and Kono, 1987; Makse et al., 2004; Shäfer et al., 1996; Silbert et al., 2001). In particular, Kuwabara and Kono (1987) and Brilliantov et al. (1996) derived the normal force between two spherical granular grains, assuming such grains are viscoelastic. The frictional part in their derived normal force is proportional to $\sqrt{{R}_{\mathrm{red}}}\sqrt{\delta}v$, where ${R}_{\mathrm{red}}$ is the reduced radius and equals $R/2$ if both spheres have the same radius $R,\delta $ is the amount of overlap of the grains (same as our ${\delta}_{\mathrm{cc}}$), and $v$ is the relative speed of the headon collision of the grains (same as our ${\overrightarrow{v}}_{\mathrm{cc}}\cdot {\overrightarrow{n}}_{\mathrm{cc}}$). Since we have always used a fixed radius of the spherical caps of a cell in our simulations, with or without the factor $\sqrt{{R}_{\mathrm{red}}}$ does not change noticeably our simulation results. However, unlike in granular material simulations where the grain size is fixed, a cell in the bacterial colony can have different mass and size at different stages. So, including the dependence on the effective mass is reasonable.
The magnitude of the tangential force ${\overrightarrow{F}}_{\mathrm{cc},t}$ in Equation (A1.4.2), where ${\gamma}_{\mathrm{cc},t}$ and ${\mu}_{\mathrm{cc}}$ are the (tangential) static and dynamic friction constants, respectively, is the minimum of a tangential friction force that depends on the relative velocity ${\overrightarrow{v}}_{\mathrm{cc}}$ and the dynamic friction force that is proportional to the Hertz contact force which is part of the normal force ${\overrightarrow{F}}_{\mathrm{cc},n}$ in Equation (A1.4.1). This is different from that used in Volfson et al. (2008): we only include the Hertz contact force here, not the entire normal force, as the elastic collision force can be dominated in cellcell interactions.
To understand the form of the tangential component of the friction force, let us first consider the friction for two solid objects in contact. Coulomb’s law for such friction states that the static tangential friction force ${\overrightarrow{F}}_{\mathrm{static}}$ required to move the objects relative to each other and the dynamic tangental friction force ${\overrightarrow{F}}_{\mathrm{dynamic}}$ between the objects to maintain such motion once it is initiated are both proportional to the normal force ${F}_{\mathrm{normal}}$ pressing the objects together (Popov, 2010):
where ${\mu}_{\mathrm{s}}$ and ${\mu}_{\mathrm{d}}$ are the static and dynamics friction coefficients, respectively. In general, $\mu}_{\mathrm{d}}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{\mu}_{\mathrm{s}$ and ${\mu}_{\mathrm{d}}\approx {\mu}_{\mathrm{s}};$ cf. Appendix 1—figure 2A. Here we assume for simplicity that ${\mu}_{\mathrm{d}}={\mu}_{\mathrm{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 $\gamma $, when the tangential relative velocity is small; cf. Appendix 1—figure 2B. This allows the inclusion in the tangential friction of the dependence of the tangential velocity when it is small in magnitude and can be damped away quickly in the dynamics. Such an assumption is supported by our experimental observation that the radial velocity of colony expansion is linear in the batch culture growth rate but is independent of the individual cell velocity, as the ‘buckling length’ is independent of such local cell velocity. Our specific form of the tangential friction for small value of relative tangential velocity, that is the first part in the minimum of the tangential force ${\overrightarrow{F}}_{\mathrm{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 ${\delta}_{\mathrm{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 ${\gamma}_{1}$ and ${\gamma}_{2}$ of the two cells and the corresponding points on these segments that reach this minimum distance. For notational convenience, let us denote by ${\overrightarrow{p}}_{i}$ and ${\overrightarrow{q}}_{i}$ the position vectors of the centers of hemispherical caps of the cell $i$ with $i=1$ and $2$, respectively. Let ${\overrightarrow{u}}_{i}={\overrightarrow{q}}_{i}{\overrightarrow{p}}_{i}(i=0,1).$ We parameterize the central cylindrical line segments ${\gamma}_{i}$ by ${\overrightarrow{r}}_{i}(t)={\overrightarrow{p}}_{i}+t{\overrightarrow{u}}_{i}$ for $0\le t\le 1(i=1,2).$ Denote ${\overrightarrow{p}}_{0}={\overrightarrow{p}}_{1}{\overrightarrow{p}}_{2}$ and define the distancesquare function
Clearly, $f$ is minimized in $[0,1]\times [0,1]$ by some point $({s}_{0},{t}_{0})\in [0,1]\times [0,1].$ The minimum distance $d\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$, and the two points ${\overrightarrow{a}}_{1}$ and ${\overrightarrow{a}}_{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 ${\overrightarrow{u}}_{1}\times {\overrightarrow{u}}_{2}\ne \overrightarrow{0}.$ In this case, $f(s,t)$ is a strictly convex and quadratic function with a constant symmetric positive definite Hessian matrix. By setting ${\partial}_{s}f(s,t)=0$ and ${\partial}_{t}f(s,t)=0$, we find that $f$ is minimized in ${\mathbb{R}}^{2}$ at $({s}_{0},{t}_{0})$ with
If $({s}_{0},{t}_{0})\in [0,1]\times [0,1]$, then we obtain $d$ and ${\overrightarrow{a}}_{1},{\overrightarrow{a}}_{2}$ by Equation (A1.4.3). Otherwise, we compute the minimum value of $f$ on each of the four sides of the square $[0,1]\times [0,1]$ and compare these values to find the global minimum points $({s}_{0},{t}_{0})\in [0,1]\times [0,1]$ and the corresponding minimum value of $f$ on this square. Consider, for instance, the side $s=0$ and $0\le t\le 1.$ The function $f(0,t)$ for all $t\in \mathbb{R}$ is found to be minimized at ${t}_{1}=({\overrightarrow{a}}_{0}\cdot {\overrightarrow{u}}_{2})/{\parallel {\overrightarrow{u}}_{2}\parallel}^{2}$. If ${t}_{1}\in [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 ${\gamma}_{1}$ and ${\gamma}_{2}$ are parallel. In this case, the minimum distance $d$ is achieved by the minimum of the distance from ${\overrightarrow{p}}_{1}$ to the second line segment ${\gamma}_{2}$ and that from ${\overrightarrow{q}}_{1}$ to ${\gamma}_{2}$:
Each distance from a point to a line segment can be found by minimizing a convex quadratic function on $[0,1]$, similar to and simpler than the previous case.
(b) Cellagar 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 cellagar interaction force ${\overrightarrow{F}}_{\mathrm{ca}}$ is generated and can be modeled similar to the cellcell interaction force. Assume one end of the cell dips into the agar region; cf. Figure 13B in the main text. Let $\overrightarrow{p}=({x}_{a},{y}_{a},{z}_{a})$ be the center of the spherical cap corresponding to that end of the cell. We denote by ${\delta}_{\mathrm{ca}}$ the indentation depth: ${\delta}_{\mathrm{ca}}={w}_{0}/2{z}_{a}.$ We also denote ${\overrightarrow{r}}_{\mathrm{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 point $\overrightarrow{p}$ between $z=0$ and $z={\delta}_{\mathrm{ca}}.$ As before, we denote by ${\overrightarrow{v}}_{\mathrm{ca}}$ the relative velocity at the center ${\overrightarrow{r}}_{\mathrm{ca}}$. It is given by
where $\overrightarrow{v}$ is the velocity of the cell at its center ${\overrightarrow{r}}_{\mathrm{c}}$ and $\overrightarrow{\omega}$ 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 ${\overrightarrow{n}}_{\mathrm{ca}}=(0,0,1).$ The tangential direction is defined through the relative velocity ${\overrightarrow{v}}_{\mathrm{ca}}$ by
if the denominator is nonzero. (Otherwise, ${\overrightarrow{\tau}}_{\mathrm{ca}}$ can be any unit vector orthogonal to ${\overrightarrow{n}}_{\mathrm{ca}}$.)
Similar to the cellcell interaction, the total cellagar interaction force exerted on the midpoint ${\overrightarrow{r}}_{\mathrm{ca}}$ and the corresponding torque are given by
Here ${\overrightarrow{F}}_{\mathrm{ca},n}$ and ${\overrightarrow{F}}_{\mathrm{ca},t}$ are the forces normal and tangential to the mean colonyagar interface $z=0.$ They are given by
where ${k}_{\mathrm{ca}}$ is an elastic constant in the Hertzian stress, ${\gamma}_{\mathrm{ca},n}$ and ${\gamma}_{\mathrm{ca},t}$ are static friction constants, ${\mu}_{\mathrm{ca}}$ is the dynamic friction constant, and ${M}_{\mathrm{cell}}={\rho}_{\mathrm{cell}}{V}_{\mathrm{cell}}$ is the cell mass with ${V}_{\mathrm{cell}}$ the cell volume. Note that the factor $\sqrt{2}$ in the Hertz contact force part is different from those for the cellcell interaction case. Here, the mean agar surface can be treated as a sphere of radius $\mathrm{\infty}$ which leads to the reduced radius of the cellagar 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, $\mathrm{d}W={\gamma}_{\mathrm{surf}}\mathrm{d}A$, where $W$ is the work done by the water surface, ${\gamma}_{\mathrm{surf}}$ is the watervapor surface tension, and $\mathrm{d}A$ is the change in water area $A$ as the cell is raised with an additional $\mathrm{d}h$; cf. Appendix 1—figure 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=2\pi Rh$ and hence $dA=2\pi Rdh$. (In the case of a disk in the twodimensional setting, the change in area is $\mathrm{d}A=2\pi rR\mathrm{d}\theta =2\pi {R}^{2}\mathrm{sin}\theta \mathrm{d}\theta $. Note that $h=R\mathrm{sin}\theta $ and $\mathrm{d}h=R\mathrm{sin}\theta \mathrm{d}\theta $. Hence $\mathrm{d}A=2\pi R\mathrm{d}h.$) As a result, the surface tension force is
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}_{\mathrm{c}}\approx 0.1{w}_{0}$, when it saturates to its maximum value of $\pi {w}_{0}{\gamma}_{\mathrm{surf}}$.
Specifically, let us consider a cell on top of the colony. Let ${\overrightarrow{r}}_{c}$ denote the center, and $\overrightarrow{p}$ and $\overrightarrow{q}$ the two centers of hemispheres of such a cell; cf. Figure 14C in the main text. We define the total surface tension force ${\overrightarrow{F}}_{\mathrm{surf}}$ and torque ${\overrightarrow{T}}_{\mathrm{surf}}$ on this cell to be
where
Here, ${h}_{\overrightarrow{i}}$ is the water level at $\overrightarrow{i}(\overrightarrow{i}=\overrightarrow{p}$ or $\overrightarrow{q}$), ${\overrightarrow{n}}_{z}=(0,0,1)$ is the unit vector along the $z$ direction, ${z}_{\overrightarrow{i}}=\overrightarrow{i}\cdot {\overrightarrow{n}}_{z}$ is the $z$component of $\overrightarrow{i}$, ${\overrightarrow{n}}_{\mathrm{s}}$ is the unit vector of the colony surface, and $\delta h$ is a constant determining how tightly the surface tension holds the cells. The water level is coarsegrained. First, on each horizontal grid square ${\mathcal{B}}_{i,j}=(i\mathrm{\Delta}x,(i+1)\mathrm{\Delta}x)\times (j\mathrm{\Delta}y,(j+1)\mathrm{\Delta}y)$, 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 ${\mathcal{B}}_{i,j}\times [0,b]$. (Recall that $z=b$ is the top of our computational box; cf. Figure 11 in the main text.) Then we construct the coarsegrained 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 force ${\overrightarrow{F}}_{\mathrm{visc}}$ exerted at a cell by the surrounding liquid and the corresponding torque ${T}_{\mathrm{visc}}$ are given by (note that ${w}_{0}$ is the diameter)
respectively, where $\overrightarrow{v}$ and $\overrightarrow{\omega}$ are the velocity and angular velocity, respectively, at the center of mass of the cell, and ${\mu}_{\text{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 cellcell and cellagar 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:
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 CoarseGrained Variables
To better present our simulation results, we need to coarsegrain the cell volume fraction $\varphi $, pressure field $P$, velocity field $\overrightarrow{v}$, and director field $\overrightarrow{n}$ over a subregion $\mathcal{G}$ of the colony region. Examples of a subregion $\mathcal{G}$ include:
The union of a few grid boxes for coarsegraining in the entire colony;
A small box at the agarcolony interface $(i\mathrm{\Delta}x,(i+1)\mathrm{\Delta}x)\times (j\mathrm{\Delta}y,(j+1)\mathrm{\Delta}y)\times (0,\mathrm{\Delta}z)$ for some $i$, $j$ for coarsegraining around such interface;
A small box in a vertical layer $(\mathrm{\Delta}x,\mathrm{\Delta}x)\times (j\mathrm{\Delta}y,(j+1)\mathrm{\Delta}y)\times (k\mathrm{\Delta}z,(k+1)\mathrm{\Delta}z)$ or $(i\mathrm{\Delta}x,(i+1)\mathrm{\Delta}x)\times (\mathrm{\Delta}y,\mathrm{\Delta}y)\times (k\mathrm{\Delta}z,(k+1)\mathrm{\Delta}z)$ for some $i$, $j$, and $k$ for coarsegraining around the crosssection $x=0$ or $y=0$, respectively; and
A cylindrical ‘cube’ $(i\delta r,(i+1)\delta r)\times (j\delta \theta ,(j+1)\delta \theta )\times (k\mathrm{\Delta}z,(k+1)\mathrm{\Delta}z)$ in the cylindrical coordinates $(r,\theta ,z)$ for some small $\delta r\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$, $\delta \theta \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$, and $\mathrm{\Delta}z\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$, and for some $i$, $j$, and $k$, for azimuthal coarsegraining.
Now let us fix a subregion $\mathcal{G}$ in the colony. Let $f$ denote the pressure $P$ or a component of the velocity field $\overrightarrow{v},$ and denote by ${f}_{i}$ such a quantity at the center of cell $i$. We define the coarsegrained average of $f$ over $\mathcal{G}$ to be
where ${\overrightarrow{r}}_{{\mathrm{c}}_{i}}$ is the center of mass for the cell $i$. Similarly, we define the coarsegrained volume fraction over the subregion $\mathcal{G}$ to be
where ${V}_{i}$ is the volume of cell $i$ and $\mathrm{Vol}(\mathcal{G})$ is the volume of the subregion $\mathcal{G}$. Note that, if $\mathcal{G}$ is a grid box, then $\overline{\varphi}(\mathcal{G})$ is the same as the volume fraction $\varphi $ on that box; cf. section A1.2. Note also that the cell density $\rho $ can be coarsegrained following its relation with the volume fraction $\varphi $; cf. Equation (A1.2.1).
The director field $\overrightarrow{n}$ cannot be coarsegrained 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 coarsegrained director field $\overrightarrow{n}(\mathcal{G})$ over a given subregion $\mathcal{G}$ to be a maximizer of the maximization problem
By the Lagrange multiplier method, this leads to an eigenvalue problem for a 3by3 matrix, with the maximizer $\overrightarrow{n}(\mathcal{G})$ being a unit eigenvector that corresponds to the largest eigenvalue.
Let $\delta r\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ be small. Let ${N}_{\theta}\ge 2$ be an integer and define $\delta \theta =2\pi /{N}_{\theta}$. We denote ${\mathcal{\mathcal{G}}}_{i,j,k}=(i\delta r,(i+1)\delta r)\times (j\delta \theta ,(j+1)\delta \theta )\times (k\mathrm{\Delta}z,(k+1)\mathrm{\Delta}z)$ in the cylindrical coordinates $(r,\theta ,z),$ We define the azimuthal average of a scalar field $f$ by
where $\overline{f}({\mathcal{\mathcal{G}}}_{i,j,k})$ is the the coarsegrained average of $f$ over $\mathcal{\mathcal{G}}}_{i,j,k$. The azimuthal average of a vector field can be defined componentwise.
Given any point in the agarcolony interface with the polar coordinates $(r,\theta )$, we define
where $\rho ({r}^{\mathrm{\prime}},\theta )$ is the local cell density projected onto the colonyagar interface. This is the negative distance between this point and the colony edge along the ray of angle $\theta .$ For each integer $j$ with $0\le j\le {N}_{\theta}$, we denote by $i}_{j$ the largest integer not exceeding $(R(j\delta \theta )+\mathrm{\Delta}r)/\delta r$. We then define the azimuthal average of the radial component $V}_{r$ of the velocity field $\overrightarrow{v}$ by
where ${\overrightarrow{r}}_{j}=(\mathrm{cos}(j\delta \theta ),\mathrm{sin}(j\delta \theta )).$ Other azimuthalaveraged quantities in terms of $\mathrm{\Delta}r$ can also be defined similarly.
Appendix 2
Additional Results of Simulation and Analysis
Unless otherwise stated, all the notations and terms are the same as defined in the main text and Appendix 1.
A2.1 Constant Density and Volume Fraction
For each of the batch culture growth rates ${\lambda}_{\mathrm{S}}=1{\mathrm{h}}^{1}$ and $0.5{\mathrm{h}}^{1}$, we simulated the growth of bacterial colony. Figure 3—figure supplement 1 in the main text shows our simulation results on the volume fractions and their coursegrained values. We observe that the volume fraction of cells in the colony has the mean value around $0.68$ with the standard deviation around $0.03$ for both of the batch culture growth rates ${\lambda}_{\mathrm{S}}=1{\mathrm{h}}^{1}$ and $0.5{\mathrm{h}}^{1}$. This suggests that we can treat the cell volume fraction as a constant inside the colony. Hence, we can also approximate well the density $\rho $ by a constant ${\rho}_{0}$; cf. Equation (A1.2.1) in Appendix 1.
A2.2 Effect of ${\lambda}_{\mathrm{S}}$Dependence on Cell Dividing Length
It has been known that cell dividing length ${\mathrm{\ell}}_{\mathrm{div}}$ may vary with the batch culture growth rate ${\lambda}_{\mathrm{S}}$ according to the relation (cf. Donachie, 1968; Jun and TaheriAraghi, 2015; Wallden et al., 2016)
where ${V}_{\mathrm{div}}$ is the corresponding dividing volume of the cell, ${V}_{0}$ is a constant, and ${\lambda}_{\mathrm{S0}}=1.0{h}^{1}$. In our simplified computational model, we have not included such variations. Here, we study the effect of the ${\lambda}_{\mathrm{S}}$dependence on the cell dividing length ${\mathrm{\ell}}_{\mathrm{div}}$ by computer simulations.
Based on recent experimental observations on the size on E. coli cell (Jun and TaheriAraghi, 2015; Si et al., 2017), we have the ratio
for all the batch culture growth rate ${\lambda}_{\mathrm{S}}.$ This and Equation (A2.2.1), together with the formula for cell volume (cf. Equation (A1.3.1) in Appendix 1), then imply that the cell dividing length should be given by
Setting ${\lambda}_{\mathrm{S}}={\lambda}_{\mathrm{S0}}=1{h}^{1}$ in Equation (A2.2.1) and using the assumption Equation (A2.2.2), we obtain also that
where ${\mathrm{\ell}}_{\mathrm{div},0}$ is the cell dividing length for ${\lambda}_{\mathrm{S}}={\lambda}_{\mathrm{S0}}$. We now assume that this cylindrical length is ${\mathrm{\ell}}_{\mathrm{div},0}=3\mu \mathrm{m}$ (Jun and TaheriAraghi, 2015; Si et al., 2017). The above two equations then provide us with the cell dividing length
We chose ${\lambda}_{\mathrm{S}}=0.1{\mathrm{h}}^{1},0.2{\mathrm{h}}^{1},\mathrm{\dots},1.0{\mathrm{h}}^{1}.$ For each of these values of ${\lambda}_{\mathrm{S}}$, we ran simulations with the variable dividing length ${\mathrm{\ell}}_{\mathrm{div}}$ determined by Equation (A2.2.3). We also ran simulations with a fixed diving length ${\mathrm{\ell}}_{\mathrm{div}}=3\mu $m; cf. Supplementary file 1Table S3. In Figure 5C in the main text, we plot the (constant) vertical ascending speed ${V}_{\mathrm{H}}$ vs. ${\lambda}_{\mathrm{S}}$ for both fixed (open circles) and variable (filled circles) dividing lengths. We observe that the results from a fixed dividing length are consistent with those from a variable dividing length. In Figure 8F in the main text, we see that the (constant) radial velocity ${V}_{\mathrm{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}_{\mathrm{H}}.$ In Appendix 2—figure 1 below, we plot the velocities ${V}_{\mathrm{R}}$ and ${V}_{\mathrm{H}}$ for variable dividing length, and fit the simulation data with ${\lambda}_{\mathrm{S}}\ge 0.5{\mathrm{h}}^{1}$. We observe that the straight line that fits simulated ${V}_{\mathrm{R}}$ intersects the $x$axis at $\sim 0.2{\mathrm{h}}^{1}$. This agrees well with the experimental result plotted in Figure 1H in the main text. Hence, the inclusion of ${\lambda}_{\mathrm{S}}$dependence on the cell dividing length can better describe experiment.
We now show that using a fixed cell dividing length ${\mathrm{\ell}}_{\mathrm{div}}$ independent of the batch culture growth rate ${\lambda}_{\mathrm{S}}$ is a reasonable simplification for our simulations. We fixed the batch culture growth rate ${\lambda}_{\mathrm{S}}=1.0{\mathrm{h}}^{1}$ and the constant concentration in the boundary condition ${C}_{\mathrm{s}}=2.0\mathrm{mM}.$ We then distributed randomly $625$ cells at time $t=0\mathrm{h}.$ At $t=9.0\mathrm{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\mathrm{h}$ to $t=15\mathrm{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 ${\lambda}_{\mathrm{S}}$, to low values, smaller than $10\%$ of ${\lambda}_{\mathrm{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 hightolow growth rate transition. To better understand such sharpness, we selected randomly $100$ cells which complete the hightolow transition, and tracked the growth rate of each of these cells during the time period from $t=9.0\mathrm{h}$ to $t=16\mathrm{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 ${\lambda}_{\mathrm{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}_{\mathrm{div}}$ and ${w}_{0}$ for the entire simulation.
A2.3 A OneDimensional Model for Nutrient Penetration
The vertical ascending velocity and the pattern formation of the colony depend largely on how deep the nutrient from agar can penetrate into the interior of colony. Near the agar surface, cells grow fast with abundant nutrient supply. Away from the agar surface, the growth of cells in the colony is much limited due to the lack of nutrient. In this section, we construct a simplified 1D model of nutrient diffusion and reaction, and show that the nutrient concentration decays quadratically in the region where it is above the Monod constant ${K}_{\mathrm{S}}$ but decays exponentially in the region where it is below ${K}_{\mathrm{S}}$. We introduce the nutrient penetration level ${H}_{\mathrm{S}}$ to be the height (i.e., the $z$ coordinate) at which the nutrient concentration takes the value ${K}_{\mathrm{S}}$, and analyze how ${H}_{\mathrm{S}}$ affects the colony growth.
Since the colony is rather thin, the variations of the nutrient concentration in the $x$ and $y$ directions are rather small compared with that in the $z$ direction. Therefore, we can assume that the steadystate concentration $C=C(z)$ depends only on the $z$ variable, and approximate the full Laplacian $\mathrm{\Delta}C$ in the colony region by the 1D Laplacian ${\partial}_{zz}C$ in the $z$ variable. Consequently, we consider the following 1D model for the nutrient concentration $C=C(z):$
where a prime denotes the spatial derivative and ${C}_{0}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ is a controlling parameter that describes the amount of nutrient concentration from the agar substrate. Note that we have replaced the local cell density $\rho $ by the constant density ${\rho}_{0}$; cf. the description of nutrient update in Appendix A1.2. Setting
we obtain the nondimensionalized form of the model
where ${\stackrel{~}{C}}_{0}={C}_{0}/{K}_{\mathrm{S}}.$
We observe that there is a unique solution $\stackrel{~}{C}=\stackrel{~}{C}(\stackrel{~}{z})$ to this boundaryvalue problem that is nonnegative for all $\stackrel{~}{z}\ge 0.$ This solution is the unique minimizer of the convex functional
of all nonnegative functions $u$ such that both $u$ and ${u}^{\prime}$ are squareintegrable on $(0,\mathrm{\infty})$, and $u(0)=\stackrel{~}{C}(0)$ and $u(\mathrm{\infty})=0.$ Note that the solution $\stackrel{~}{C}=\stackrel{~}{C}(\stackrel{~}{z})$ is a monotonically decreasing function of $\stackrel{~}{z}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$. For otherwise, $\stackrel{~}{C}$ would reach a local maximum at some ${\stackrel{~}{z}}_{\mathrm{m}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ with ${\stackrel{~}{C}}^{\mathrm{\prime}\mathrm{\prime}}({\stackrel{~}{z}}_{\mathrm{m}})\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0$ but $\stackrel{~}{C}({\stackrel{~}{z}}_{\mathrm{m}})\ge 0$. This is impossible by Equation (A2.3.1). We also observe that ${\stackrel{~}{C}}^{\prime}(\mathrm{\infty})=0$, for otherwise ${\stackrel{~}{C}}^{\prime}$ would be negative and stay strictly away from 0 as $\stackrel{~}{{C}^{\prime \prime}}\ge 0$, leading eventually to $\stackrel{~}{C}(\mathrm{\infty})=\mathrm{\infty}$ which would contradict the fact that $\stackrel{~}{C}(\mathrm{\infty})=0.$
Now, multiplying both sides of Equation (A2.3.1) by ${\stackrel{~}{C}}^{\prime}$, we get
Integrating both sides of this equation from 0 to $\stackrel{~}{z}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$, we obtain
Sending $\stackrel{~}{z}\to \mathrm{\infty}$, we get
The combination of the above two equations leads to
We now study how fast $\stackrel{~}{C}=\stackrel{~}{C}(\stackrel{~}{z})$ decays. Assume $\stackrel{~}{C}({z}_{1})\le 1$ for some ${\stackrel{~}{z}}_{1}\ge 0$. Then $\stackrel{~}{C}(z)\le 1$ for all $\stackrel{~}{z}\ge {\stackrel{~}{z}}_{1}.$ Noting that
and that ${\stackrel{~}{C}}^{\prime}\le 0$, we have by Equation (A2.3.3) that ${\stackrel{~}{C}}^{\prime}\ge \sqrt{2/3}\stackrel{~}{C}$ for all $\stackrel{~}{z}\ge {\stackrel{~}{z}}_{1}$. This leads to the exponential decay
If ${\stackrel{~}{C}}_{0}=\stackrel{~}{C}(0)\le 1$, then we can have ${\stackrel{~}{z}}_{1}=0$ in Equation (A2.3.4) and the concentration $\stackrel{~}{C}=\stackrel{~}{C}(\stackrel{~}{z})$ decays exponentially.
Suppose now ${\stackrel{~}{C}}_{0}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$. Since $\stackrel{~}{C}=\stackrel{~}{C}(\stackrel{~}{z})$ decreases monotonically and $\stackrel{~}{C}(\mathrm{\infty})=0$, there exists a unique ${\stackrel{~}{z}}_{0}={z}_{0}/\beta \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ such that $\stackrel{~}{C}({\stackrel{~}{z}}_{0})=1$ (i.e., $C({z}_{0}/\beta )={K}_{\mathrm{S}}$) and $\stackrel{~}{C}(\stackrel{~}{z})\le 1$ for all $\stackrel{~}{z}\ge {\stackrel{~}{z}}_{0}.$ Thus the inequality (Equation (A2.3.4)) holds with ${\stackrel{~}{z}}_{1}={\stackrel{~}{z}}_{0}$ and $\stackrel{~}{C}({\stackrel{~}{z}}_{1})=\stackrel{~}{C}({\stackrel{~}{z}}_{0})=1.$ We show now that the concentration $\stackrel{~}{C}=\stackrel{~}{C}(\stackrel{~}{z})$ decays quadratically in $[0,{\stackrel{~}{z}}_{0}].$ Precisely, we shall prove that
By Equation (A2.3.3) and the fact that ${\stackrel{~}{C}}^{\prime}\le 0$, we have ${\stackrel{~}{C}}^{\prime}\le \sqrt{2\stackrel{~}{C}}$ in $[0,{\stackrel{~}{z}}_{0}]$. This leads to ${(\sqrt{\stackrel{~}{C}})}^{\prime}\ge 1/\sqrt{2}$. Integrating both sides of this inequality from 0 to $\stackrel{~}{z}$, we obtain the first inequality in Equation (A2.3.5). Note that
Thus, by Equation (A2.3.3) and the fact that ${\stackrel{~}{C}}^{\prime}\le 0$, we obtain that ${\stackrel{~}{C}}^{\prime}\ge \sqrt{2\mathrm{ln}(e/2)\stackrel{~}{C}}$, and further that ${(\sqrt{\stackrel{~}{C}})}^{\prime}\le \sqrt{(\mathrm{ln}(e/2))/2}$ in $[0,{\stackrel{~}{z}}_{0}]$. An integration of both sides of this inequality from 0 to $\stackrel{~}{z}$ then leads to the second inequality in Equation (A2.3.5).
In Figure 5—figure supplement 1 in the main text, we see that the nutrient concentration ${C}_{\mathrm{ctr}}$ along the $z$axis is described well by a quadratic function for ${C}_{\mathrm{ctr}}\ge {K}_{\mathrm{S}}$ (i.e., $\stackrel{~}{C}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$) and that the nutrient concentration decays exponentially for ${C}_{\mathrm{ctr}}\le {K}_{\mathrm{S}}$ (i.e., $\stackrel{~}{C}\le 1$).
The position ${\stackrel{~}{z}}_{0}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$, defined by $\stackrel{~}{C}({\stackrel{~}{z}}_{0})=1$, is the (rescaled) vertical level across which the nutrient concentration transitions from the quadratic decay described in Equation (A2.3.5) to the exponential decay described by Equation (A2.3.4) with ${\stackrel{~}{z}}_{1}={\stackrel{~}{z}}_{0}.$ The nutrient penetration level we have defined is exactly ${H}_{\mathrm{S}}=\beta {\stackrel{~}{z}}_{0}$. Setting $\stackrel{~}{z}={\stackrel{~}{z}}_{0}$ in Equation (A2.3.5), we see that ${\stackrel{~}{z}}_{0}=O(\sqrt{{\stackrel{~}{C}}_{0}})$ if ${\stackrel{~}{C}}_{0}\gg 1.$ Therefore,
Since the nutrient decays exponentially above ${H}_{\mathrm{S}}$, only those cells below the height ${H}_{\mathrm{S}}$ grow with the maximum growth rate ${\lambda}_{\mathrm{S}}$. Therefore, the speed of vertical expansion ${V}_{\mathrm{H}}$ is given by ${V}_{\mathrm{H}}\propto {H}_{\mathrm{S}}{\lambda}_{\mathrm{S}}\propto \sqrt{{\lambda}_{\mathrm{S}}};$ see the discussion in the main text: Vertical rise—quantitative analysis, in Simulation Results and Analysis.
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/elifesciencespublications/CellsMD3D). Experimental source data files are provided for Figure1, Figure 1figure supplement 1, and Figure 9.
References

Threedimensional biofilm model with individual cells and continuum EPS matrixBiotechnology and Bioengineering 94:961–979.https://doi.org/10.1002/bit.20917

Inflating bacterial cells by increased protein synthesisMolecular Systems Biology 11:836.https://doi.org/10.15252/msb.20156178

Cooperative selforganization of microorganismsAdvances in Physics 49:395–554.https://doi.org/10.1080/000187300405228

Verticalization of bacterial biofilmsNature Physics 14:954–960.https://doi.org/10.1038/s4156701801704

Determination of biofilm diffusion coefficients using microelectrodesProgress in Biotechnology 11:31–38.

Buckling instability in ordered bacterial coloniesPhysical Biology 8:26008.https://doi.org/10.1088/14783975/8/2/026008

Biofilms: the matrix revisitedTrends in Microbiology 13:20–26.https://doi.org/10.1016/j.tim.2004.11.006

Model for collisions in granular gasesPhysical Review E 53:5382–5392.https://doi.org/10.1103/PhysRevE.53.5382

Factors affecting the growth of bacterial colonies on agar platesProceedings of the Royal Society of London. Series B, Biological sciences 171:175–199.https://doi.org/10.1098/rspb.1968.0063

Microbial biofilmsAnnual Review of Microbiology 49:711–745.https://doi.org/10.1146/annurev.mi.49.100195.003431

Differential growth of wrinkled biofilmsPhysical Review E 91:22710.https://doi.org/10.1103/PhysRevE.91.022710

Mechanically driven growth of quasitwodimensional microbial coloniesPhysical Review Letters 111:168101.https://doi.org/10.1103/PhysRevLett.111.168101

The wave of advance of advantageous genesAnnals of Eugenics 7:355–369.https://doi.org/10.1111/j.14691809.1937.tb02153.x

INDISIM, an individualbased discrete simulation model to study bacterial culturesJournal of Theoretical Biology 214:305–319.https://doi.org/10.1006/jtbi.2001.2466

The role of mechanical forces in the planartobulk transition in growing Escherichia coli microcoloniesJournal of The Royal Society Interface 11:20140400.https://doi.org/10.1098/rsif.2014.0400

A continuum model for the growth of bacterial colonies on a surfaceJournal of Physics A: Mathematical and General 26:5645–5654.https://doi.org/10.1088/03054470/26/21/006

Modeling granular media on the computerContinuum Mechanics and Thermodynamics 10:189–231.https://doi.org/10.1007/s001610050089

Ueber die Berührung fester elastischer KörperJ fur die Reine und Angew Math 1882:156–171.

Bacterial quorum sensing: signals, circuits, and implications for biofilms and diseaseAnnual Review of Biomedical Engineering 10:145–167.https://doi.org/10.1146/annurev.bioeng.10.061807.160536

Cellsize maintenance: universal strategy revealedTrends in Microbiology 23:4–6.https://doi.org/10.1016/j.tim.2014.12.001

Needbased activation of ammonium uptake in Escherichia coliMolecular Systems Biology 8:616.https://doi.org/10.1038/msb.2012.46

Finger Formation in Biofilm LayersSIAM Journal on Applied Mathematics 62:853–869.https://doi.org/10.1137/S0036139900371709

Investigation of the equation of diffusion combined with increasing of the substance and its application to a biology problemBulletin of Moscow State University Series A: Mathematics and Mechanics 1:1–25.

Selective sweeps in growing microbial coloniesPhysical Biology 9:26008.https://doi.org/10.1088/14783975/9/2/026008

Individualbased modelling of biofilmsMicrobiology 147:2897–2912.https://doi.org/10.1099/00221287147112897

Restitution coefficient in a collision between two spheresJapanese Journal of Applied Physics 26:1230–1233.https://doi.org/10.1143/JJAP.26.1230

The influence of nutrition and temperature on the growth of colonies of Escherichia coli K12Canadian Journal of Microbiology 27:679–684.https://doi.org/10.1139/m81105

The Growth of Bacterial CulturesAnnual Review of Microbiology 3:371–394.https://doi.org/10.1146/annurev.mi.03.100149.002103

BookNumerical Solution of Partial Differential EquationsCambridge University Press.

Spatial structure, cooperation and competition in biofilmsNature Reviews Microbiology 14:589–600.https://doi.org/10.1038/nrmicro.2016.84

BookCell growth, genome duplication and cell divisionIn: Nanninga N, editors. Molecular Cytology of Escherichia Coli. London, New York: . Academic Press. pp. 259–318.

Biofilm formation as microbial developmentAnnual Review of Microbiology 54:49–79.https://doi.org/10.1146/annurev.micro.54.1.49

Growth measurements on surface colonies of bacteriaJournal of general microbiology 66:137–143.https://doi.org/10.1099/00221287662137

Energy yields and growth of heterotrophsAnnual Review of Microbiology 24:17–52.https://doi.org/10.1146/annurev.mi.24.100170.000313

Particlebased multidimensional multispecies biofilm modelApplied and Environmental Microbiology 70:3024–3040.https://doi.org/10.1128/AEM.70.5.30243040.2004

A kinetic study of the mode of growth of surface colonies of bacteria and fungiJournal of General Microbiology 47:181–197.https://doi.org/10.1099/00221287472181

Localization of active microorganisms in cheese by autoradiographyApplied and environmental microbiology 38:1162–1165.

The Limitations of In Vitro Experimentation in Understanding Biofilms and Chronic InfectionJournal of Molecular Biology 427:3646–3661.https://doi.org/10.1016/j.jmb.2015.09.002

Computational modeling of synthetic microbial biofilmsACS Synthetic Biology 1:345–352.https://doi.org/10.1021/sb300031n

Cellulose as an architectural element in spatially structured Escherichia coli biofilmsJournal of Bacteriology 195:5540–5554.https://doi.org/10.1128/JB.0094613

Force Schemes in Simulations of Granular MaterialsJournal de Physique I 6:5–20.https://doi.org/10.1051/jp1:1996129

Physiological heterogeneity in biofilmsNature Reviews Microbiology 6:199–210.https://doi.org/10.1038/nrmicro1838

Biofilms as complex differentiated communitiesAnnual Review of Microbiology 56:187–209.https://doi.org/10.1146/annurev.micro.56.012302.160705

Cellsize control and homeostasis in bacteriaCurrent Biology 25:385–391.https://doi.org/10.1016/j.cub.2014.12.009

Multicomponent model of deformation and detachment of a biofilm under fluid flowJournal of The Royal Society Interface 12:20150045.https://doi.org/10.1098/rsif.2015.0045

The growth and form of bacterial coloniesJournal of General Microbiology 114:483–486.https://doi.org/10.1099/002212871142483

Phase Field Models for Biofilms. I. Theory and OneDimensional SimulationsSIAM Journal on Applied Mathematics 69:641–669.https://doi.org/10.1137/070691966
Article and author information
Author details
Funding
National Science Foundation (DMS1319731)
 Bo Li
National Science Foundation (DMS1620487)
 Bo Li
Natural Sciences and Engineering Research Council of Canada
 Mya R Warren
Simons Foundation (542387)
 Terence Hwa
Simons Foundation (522790)
 Hui Sun
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work was partially supported by the NSF through grants DMS1319731 and DMS1620487 to BL, by Simons Foundation through grants #542387 to TH and #522790 to HS. MRW was supported by a Canadian Natural Sciences and Engineering Research Council postdoctoral fellowship. This work used the NSF Extreme Science and Engineering Discovery Environment (XSEDE) through grant ACI1548562. We thank Agnese Seminara for useful discussions. We also thank the reviewers to bring to our attention the references Beroz et al. (Nat. Phys. 2018) and EnosBerlage and McCarter (J. Bacteriol. 2000).
Version history
 Received: August 15, 2018
 Accepted: February 20, 2019
 Version of Record published: March 11, 2019 (version 1)
Copyright
© 2019, Warren et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 7,060
 Page views

 1,000
 Downloads

 60
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cell Biology
 Computational and Systems Biology
Computer models of the human ventricular cardiomyocyte action potential (AP) have reached a level of detail and maturity that has led to an increasing number of applications in the pharmaceutical sector. However, interfacing the models with experimental data can become a significant computational burden. To mitigate the computational burden, the present study introduces a neural network (NN) that emulates the AP for given maximum conductances of selected ion channels, pumps, and exchangers. Its applicability in pharmacological studies was tested on synthetic and experimental data. The NN emulator potentially enables massive speedups compared to regular simulations and the forward problem (find drugged AP for pharmacological parameters defined as scaling factors of control maximum conductances) on synthetic data could be solved with average rootmeansquare errors (RMSE) of 0.47 mV in normal APs and of 14.5 mV in abnormal APs exhibiting early afterdepolarizations (72.5% of the emulated APs were alining with the abnormality, and the substantial majority of the remaining APs demonstrated pronounced proximity). This demonstrates not only very fast and mostly very accurate AP emulations but also the capability of accounting for discontinuities, a major advantage over existing emulation strategies. Furthermore, the inverse problem (find pharmacological parameters for control and drugged APs through optimization) on synthetic data could be solved with high accuracy shown by a maximum RMSE of 0.22 in the estimated pharmacological parameters. However, notable mismatches were observed between pharmacological parameters estimated from experimental data and distributions obtained from the Comprehensive in vitro Proarrhythmia Assay initiative. This reveals larger inaccuracies which can be attributed particularly to the fact that small tissue preparations were studied while the emulator was trained on single cardiomyocyte data. Overall, our study highlights the potential of NN emulators as powerful tool for an increased efficiency in future quantitative systems pharmacology studies.

 Computational and Systems Biology
 Neuroscience
Closedloop neuronal stimulation has a strong therapeutic potential for neurological disorders such as Parkinson’s disease. However, at the moment, standard stimulation protocols rely on continuous openloop stimulation and the design of adaptive controllers is an active field of research. Delayed feedback control (DFC), a popular method used to control chaotic systems, has been proposed as a closedloop technique for desynchronisation of neuronal populations but, so far, was only tested in computational studies. We implement DFC for the first time in neuronal populations and access its efficacy in disrupting unwanted neuronal oscillations. To analyse in detail the performance of this activity control algorithm, we used specialised in vitro platforms with high spatiotemporal monitoring/stimulating capabilities. We show that the conventional DFC in fact worsens the neuronal population oscillatory behaviour, which was never reported before. Conversely, we present an improved control algorithm, adaptive DFC (aDFC), which monitors the ongoing oscillation periodicity and selftunes accordingly. aDFC effectively disrupts collective neuronal oscillations restoring a more physiological state. Overall, these results support aDFC as a better candidate for therapeutic closedloop brain stimulation.