Content area
Bacterial colonies can form a wide variety of shapes and structures based on ambient and internal conditions. To help understand the mechanisms that determine the structure of and the diversity within these colonies, various numerical modeling techniques have been applied. The most commonly used ones are continuum models, agent-based models, and lattice models. Continuum models are usually computationally fast, but disregard information at the level of the individual, which can be crucial to understanding diversity in a colony. Agent-based models resolve local details to a greater level, but are computationally costly. Lattice-based approaches strike a balance between these two limiting cases. However, this is known to come at the price of introducing undesirable artifacts into the structure of the colonies. For instance, square lattices tend to produce square colonies even where an isotropic shape is expected. Here, we aim to overcome these limitations and we therefore study lattice-induced orientational symmetry in a class of hybrid numerical methods that combine aspects of lattice-based and continuum descriptions. We characterize these artifacts and show that they can be circumvented through the use of a disordered lattice which derives from an unstructured fluid. The main advantage of this approach is that the lattice itself does not imbue the colony with a preferential directionality. We demonstrate that our implementation enables the study of colony growth involving millions of individuals within hours of computation time on an ordinary desktop computer, while retaining many of the desirable features of agent-based models. Furthermore, our method can be readily adapted for a wide range of applications, opening up new avenues for studying the formation of colonies with diverse shapes and complex internal interactions.
Introduction
In nature, bacteria often form highly structured colonies or biofilms. Colonies grown in laboratory experiments reveal a wide variety of morphologies that depend on environmental conditions, such as substrate hardness and nutrient availability. For example, Bacillus subtilis forms disks, branches, concentric rings, and even a chiral branched pattern, where each branch is twisted with the same handedness [1–4]. In some cases, such structures may play a functional role, e.g., they may affect susceptibility to antimicrobial agents [5–9], prevent invasion by other bacteria [10] or increase resistance to predation [11]. In addition, the development and morphology of colonies has been shown to be important for evolutionary processes. Initially well-mixed populations of Escherichia coli form monoclonal patches as they grow out [12], and the colony morphology affects the genetic diversity that is maintained [13]. The position of individuals in the colony and thus the spatial structure of the bacterial population are also important for competition between bacterial strains [14,15] and their “social” interactions, including between antibiotic susceptible and resistant strains [9]. For these reasons, among others, the mechanisms of colony and biofilm growth remain topics of broad interest.
Experiments have provided many valuable insights into these spatially structured processes, but they can be technically challenging to design and control. As a result, numerical modeling has often been used to complement such studies and gain deeper understanding. Three broad approaches to this modeling have emerged: (i) Early studies into colony morphology formation employed continuous models, where bacteria are represented by a density field [16–20]. These models inherently only consider locally averaged and population-level properties of bacteria and specifically discard their individual nature. This is a severe limitation in studying evolutionary processes, where the genotypes and phenotypes of individual organisms play a crucial role. (ii) Individual bacteria are modeled directly in agent-based models. This has the advantage of having full control over internal states and tracking of individuals [13,21–24]. In addition, features such as bacterial shape and mechanical interaction can be described in detail by prescribing the equations of motion [25–31]. However, agent-based modeling is typically computationally costly [13,32]. (iii) An intermediate solution is to constrain the bacteria to move on a lattice [33–46]. This sacrifices some details of the local interactions in favor of a substantial speed advantage, as it avoids resolving the full dynamics of the individual bacteria. However, constraining bacteria to a lattice may introduce undesirable features in the solution [35,47], henceforth broadly referred to as ‘lattice artifacts’.
While these artifacts are often ignored, some earlier research has recognized them, in models with two [35,48] and three spatial dimensions [48]. Earlier reserach has also addressed these artifacts by using various disordered lattices [19,20,35,49,50]. However, these disordered lattices have in common that they are generated by perturbing a square reference lattice. In addition, potential lattice artifacts have not been quantified in the various models. In this work, we therefore revisit the problem of lattice-based modeling of bacterial colony growth, with the specific goal to simulate a large number of individuals while overcoming lattice-induced symmetry artifacts. Specifically, we explore a class of hybrid lattice-based models [51–54], develop a method to detect and quantify lattice artifacts in the shape of the simulated bacterial colonies using a convex hull, and propose a new way of generating disordered lattices based on the structure of a fluid.
In brief, in this paper we will show that implementing a hybrid lattice-based model on the commonly used square lattice does result in significant lattice artifacts. These cannot be easily eliminated by tuning the neighborhoods of lattice sites. We then compare previously proposed disordered lattices and demonstrate that while these lattices are effective at removing lattice artifacts in models that only take the binary neighborhood relationship between lattice sites into account, they do not fully eliminate them in models that include additional geometric properties of the lattice, such as distances between sites. Additionally, we highlight the challenge of detecting lattice artifacts through visual inspection of the resulting colony shapes alone. Next, we introduce an off-lattice fluid-derived approach to generate disordered lattices and show that, using this method, no lattice-induced symmetries are detectable across the range of models that we have considered. Finally, we demonstrate that hybrid lattice-based models can achieve colony sizes that are difficult to simulate with off-lattice agent-based models. We conclude by discussing how our approach facilitates the simulation of diversification and competition in large-scale colonies.
Results
We are interested in simulating a class of hybrid lattice-based models, where one or more continuous fields are coupled to discrete bacteria. The bacteria can have a variety of properties and interactions. An overview of this class of hybrid models is given in Fig 1. In these hybrid models, a lattice is used to constrain the positions of the bacteria (green disks in Fig 1) and to numerically integrate the continuous fields (orange, yellow, blue, and green slices in Fig 1). Below, we explore multiple approaches to define lattices and which lattice sites are neighbors. We evaluate these lattices with a minimal model of bacterial colony development, the full details of which are given in the Methods section. This minimal hybrid lattice-based model consists of bacteria coupled to a single continuous nutrient field, governed by a reaction-diffusion equation. The bacteria consume nutrients, i.e., they acts as sinks to this field. Once a bacterium has consumed an amount of nutrients above a threshold ng, it can divide and produce one daughter cell. This daughter cell is randomly placed on any empty neighboring lattice site. The bacterium is unable to divide when there are no empty neighboring lattice sites available for the daughter cell. During division, the nutrients of the parent cell are split stochastically between parent and daughter.
[Figure omitted. See PDF.]
The tilted white field illustrates how bacterial agents (green disks) are constrained to a lattice (black dots). The behavior of the bacteria is determined by various properties and internal processes; the green box lists potential examples. The red lines indicate the Voronoi cells associated with the lattice sites, which define the neighborhood and are used to compute the dynamics of the continuous fields. The bacteria can be coupled to multiple fields, as illustrated by the orange, yellow, blue, and green slices. For example, the orange slice represents the concentration field of a nutrient that is consumed by the bacteria. Its dynamics are governed by a diffusion equation, for which the bacteria serve as sinks, here indicated as gray disks.
The remainder of this section is structured as follows. First, we demonstrate that square lattices induce artifacts on resulting colony shapes. We also discuss our method to quantify these lattice artifacts. Next, we show that disordered lattices can diminish the observed lattice artifacts. We investigate two previously proposed methods for generating such a disordered lattice and characterize the quality of the solutions. We notice that there are shortcomings in certain cases, which inspires us to propose a new, fluid-derived method, which we compare in detail with the established methods. We conclude by demonstrating the computational performance of the minimal hybrid model by comparing it directly to an off-lattice (agent-based) approach.
Artifacts caused by using square lattices
We started by performing simulations of the hybrid model on a square lattice, where we only include nearest neighbors. This choice is also known as the Von Neumann neighborhood and is shown in the inset of Fig 2. Two distinct environmental conditions were considered by changing the initial (reduced) nutrient concentration c0: a nutrient-poor environment with c0 = 0.7 and a nutrient-rich environment with c0 = 3.0. The simulation units and remaining parameters can be found in Table 1.
[Figure omitted. See PDF.]
The nearest neighbors on a square grid are as indicated by the blue arrows in the top-left inset. (a) Simulation with low initial nutrient concentration () where is the minimum amount of nutrient that a bacterium needs to consume in order to divide. (b) Simulation with high initial nutrient concentration (). The remaining parameter choices are specified in Table 1. The green shapes in the center of both figures show a single representative colony, where the color of each bacterium indicates the number of divisions in the lineage of the bacterium (central color bar). The maximum division number is 131 for colony (a) and 106 for colony (b). The outer ring shows a histogram for the distribution of normal vectors to the convex hulls fitted to colonies. The error bars show the standard error of the mean resulting from 500 independent realizations. The gray scale serves to guide the eye to regions of high incidence. The red circle provides the ideal, isotropic, distribution with value . In both figures, the histograms indicate that the colonies tend to have a diamond shape.
[Figure omitted. See PDF.]
The results of these simulations are given in Fig 2. We can distinguish two colony morphologies: branched in a nutrient-poor environment (Fig 2a) and dense in a nutrient-rich environment (Fig 2b). These morphologies look qualitatively similar to the morphologies observed for B. subtilis colonies [1,56], and in other numerical studies [16,20,21,55,57]. However, it is visually clear that both colonies possess an overall diamond shape. This diamond shape is clearly an artifact from the underlying lattice, in line with previous work on this topic [35,48,58].
To quantify the discrete orientational symmetries in Fig 2, we repeated these simulations 500 times with random seeds and fitted a convex hull to each of the resulting colonies. We then determined the normal vector to each segment of the convex hulls. Subsequently, we created a histogram of the angles of these vectors, where each vector was weighted by the length of the associated convex hull segment. This histogram is shown in the outer ring of the panels of Fig 2. A more detailed description of this procedure is given in the Methods section along with an overview in S2 Figure. These distributions show that the shapes of the colonies in both the nutrient-poor and nutrient-rich environments are significantly different from the isotropic situation, which is indicated by the red circle. Furthermore, the previously observed four-fold symmetry of the diamond shape is clearly exposed.
We attempted to diminish these lattice artifacts by changing the definition of the neighborhood of a lattice site. Instead of the Von Neumann neighborhood, we now use the Moore neighborhood, which includes both the four nearest neighbors and the four next-nearest (diagonal) neighbors. To tune the relative frequency with which daughter cells are placed on a next-nearest neighbor site, we assign a weight to the next-nearest neighbor. Whenever division occurs, an empty site is then selected from the ones neighboring the parent cell, with the probability of the next-nearest neighbor sites weighted by . The value of affects the shape of the resulting colony and can be tuned to reduce the lattice artifacts.
In Fig 3, a selection of dense colonies with high initial nutrient concentration (c0 = 3) and of branched colonies with low initial nutrient concentration (c0 = 0.7) simulated with various choices of are given. It should be noted that we only apply the Moore neighborhood to the growth algorithm. The reaction-diffusion equations are solved on the Von Neumann neighborhood as before. For , we obtain square colonies instead of the diamond-shaped colonies of Fig 2. As decreases, colonies with high initial nutrient concentrations become less square and more circular. In contrast, the diamond shape reemerges in colonies with low initial nutrient concentrations. We do not recover the diamond shape in the dense colonies. The reason is that, even if is infinitesimally small, growth in diagonal directions still occurs when none of the lateral sites are empty. In the less dense branched colonies, the lateral sites are more likely to remain empty, making diagonal expansion rarer at infinitesimal . Consequently, at low , the growth of the branched colony becomes similar to the growth observed in colonies where only nearest-neighbor sites are considered, as seen in Fig 2.
[Figure omitted. See PDF.]
Simulations in the top row were all performed at high initial nutrient concentration, as in Fig 2b. The bottom row shows simulations performed at low initial nutrient concentration, as in Fig 2a. The neighborhood of a lattice site includes nearest and next-nearest neighbors as indicated in the bottom-right inset. The probability that a daughter cell is placed on a next nearest-neighbor site is weighted by , where gives equal probability between placement on next-nearest and nearest-neighbor sites. The representation and parameter choices are otherwise the same as in Fig 2. Visual inspection suggests that there is no single choice for wd that can resolve discrete orientational symmetry in both growth regimes.
To assess whether lattice artifacts can be resolved with this method we applied the convex hull procedure to the case, as this seemed to provide the most circular colonies at high initial nutrient concentrations. For completeness, we considered both the nutrient-poor and nutrient-rich environments. The results are shown in Fig 4. In the nutrient-poor environment of Fig 4a, the diamond-shaped lattice artifacts we observed in Fig 2 remain, as expected. In contrast, the diamond pattern disappears in the nutrient-rich environment of Fig 4b. However, the histogram around the dense morphology in Fig 4b still shows clear discrete orientational symmetry. Even when we focus exclusively on the high initial nutrient concentration case, we were unable to tune the model to remove the orientational symmetry imposed by the lattice.
[Figure omitted. See PDF.]
The neighborhood of a lattice site is indicated by the blue arrows in the top-left inset. The weight of growth into the diagonal next nearest-neighbor sites was set to . The selection of parameters and representation is otherwise the same as in Fig 2. (a) For the low initial nutrient concentration the maximum division number is 129. (b) The maximum division number is 86 for the high initial nutrient concentration. In both figures, the histograms demonstrate that the shapes of the colonies are affected by the symmetries of the lattice.
Taken together, these findings suggest that modifying the growth algorithm alone is not sufficient to remove symmetry-imposing lattice artifacts. Instead of changing the model, a better approach might be to change the underlying lattice. While one could consider other regular lattices, the symmetries of the triangular lattice used by Kozlovsky et al. [17] in their continuous model for colony growth remain evident in some of their solutions. Therefore, we expect that other regular lattices should produce similar artifacts.
Disordered lattices remove discrete orientational symmetries in colony growth
The artifacts that we observed are caused by the the underlying symmetries of the lattice. We used a square lattice, which has four-fold orientational order, which is expressed in the shape of the resulting colonies. A possible solution to overcome these artifacts is to use a disordered lattice. Ideally, such a lattice should obey certain requirements. In particular, the lattice should have a minimum distance between sites. This ensures that bacteria have a minimum size, and also prevents the integration of the reaction-diffusion equations from becoming unstable [59,60]. Other than that, the lattice should not possess additional symmetries or long-range order. Creating a disordered lattice with the desired properties requires some care. We start by detailing the steps involved in generating the lattice before turning to the main result of this work.
Several methods have been proposed for generating a set of pseudo-random coordinates to construct a disordered lattice (Fig 5). The Poisson random lattice proposed by Christ et al. [61] consists of a set of randomly chosen points with uniform distribution. A limitation of this method is that there is no minimum distance between neighbors. Without such a minimum distance, the effective size of bacteria determined by the lattice varies significantly. In addition, the solution to the continuous part of the hybrid model can become numerically unstable. Moukarzel and Herrmann [49] improved on this with their vectorizable random lattice (VRL). A schematic overview of this lattice is given in Fig 5a. The pseudo-random coordinates of the lattice are generated by placing one site of the disordered lattice randomly in each of the cells of a square reference lattice. A minimum distance between sites l0 can be imposed by restricting the placement domains to a subsquare centered in each cell. We will refer to this lattice as the restricted VRL. Tucker [50] adapted the VRL by changing the way l0 is determined. Instead of restricting where a site can be placed within each reference cell, the positions of two sites are redrawn when they are too close together. This lattice will be referred to as the redrawn VRL and is shown in Fig 5b.
[Figure omitted. See PDF.]
The red dots indicate the randomly drawn lattice sites and the orange lines indicate the associated Voronoi tessellation in all figures. For the (a) restricted and (b) redrawn Vectorizable Random Lattice (VRL) proposed by Moukarzel and Herrmann [49], and Tucker [50] respectively, the gray squares indicate the reference cells. Each reference cell contains exactly one lattice site. The blue double arrow indicates the tuneable minimum separation between lattice sites l0. In (a) the minimal separation l0 is ensured by randomly placing the lattice sites inside the subdomains of the reference cell indicated as green squares. Our new method (c) uses snapshots of an off-lattice fluid of bidisperse soft disks to generate the disordered grid.
Another system that is known to be disordered and have a well behaved particle-particle separation is a dense fluid of disks [62]. Therefore, we turned to a standard model fluid to generate the coordinates to our grid. We simulated a dense two-dimensional (2D) system of bidisperse Weeks-Chandler-Anderson (WCA) particles [63]. We chose a bidisperse system to prevent crystallization at high densities [64]. After equilibration, we took a snapshot of the configuration and used the particle centers of mass as sites in our disordered lattice. An overview of this fluid-derived lattice is given in Fig 5c and a more detailed description of the fluid simulations can be found in the Methods section.
For all three disordered lattices, we generated a regular Voronoi tessellation from the generated pseudo-random coordinates. This partitioned our 2D domain into cells. Two lattice sites were considered neighbors when their associated Voronoi cells share a wall.
We computed the probability density functions (PDF) of distances between neighboring lattice sites and of Voronoi cell volumes for all three types of disordered lattices to check whether they meet our requirements. Here, we define our unit of distance as the average distance between lattice sites. In addition, we investigated isotropy of the lattices by computing the pair distribution function
(1)
where is the position of lattice site i and N the total number of lattice sites. Expressed in polar coordinates, is a measure of the probability there is a site at a distance r and angle θ from another site. For isotropic lattices , where with A the total area of the lattice and g(r) the pair correlation function commonly used to describe fluids.
Fig 6 shows PDFs and pair distribution functions calculated for representative instances of the three disordered lattice types. In these examples, l0 = 0.5 was used for both versions of the VRL. We note that the redrawn VRL gives broader distributions for both the separation between neighbor sites and the volume of the Voronoi cells compared to the restricted VRL and the fluid-derived lattice (top row). The redrawn VRL also has a sharp cutoff in neighbor separations at l0. The pair distribution function (bottom row) of the restricted VRL shows clear anisotropy and long-range order. Although less pronounced, anisotropy and long-range order are also present in the redrawn VRL. The pair distribution function of the fluid-derived lattice shows no clear θ dependence, indicating that the lattice is isotropic. In addition, the fluid-derived lattice shows several oscillations at low r, indicating some short-range structure. These oscillations are typical for dense fluids and correspond to coordination layers around particles. These layers are related to the geometry of densely packed disks and would therefore also be expected in a 2D colony of disk-like bacteria. We conclude that only the fluid-derived lattice satisfies all requirements.
[Figure omitted. See PDF.]
The top row shows the probability density functions (PDFs) of the distances between neighbors (blue) and the volumes of the Voronoi cells associated to the lattice sites (orange). The bottom row gives the pair distribution functions of the lattice sites in polar coordinates on a symmetrical log-scale. For the two vectorizable random lattices (VRLs), we chose l0 = 0.5.
Next, we investigated whether the lattices can help resolve the artifacts we observed above. Since lattice artifacts can be present in both the agent-based and continuous components of the model, we first looked at the continuous part separately. We solved the bacterial growth model proposed by Kitsunezaki [19], as this is a well-established model that can serve as a test case. The Kitsunezaki model consists of three density fields: the density of active motile bacteria b, of inactive bacteria s, and of a nutrient n. These fields are governed by three reaction-diffusion equations given by
(2)(3)(4)
where μ is the rate of bacterial differentiation to the inactive state and D0 the bacterial diffusion coefficient. The non-linear diffusion term in Eq (2) is intended to model variable bacterial motility due to production of lubricant by the bacteria themselves. We set D0 = 0.1, , and the initial nutrient concentration to n0 = 1 everywhere. Our initial condition for the b field was a disk of radius 5 with b = 1 placed in the center of the simulation domain. Everywhere else we set b = 0. Finally, we initialized s = 0 everywhere. We used zero-flux boundary conditions for all fields.
We solved the Kitsunezaki model 500 times on each of the lattices to quantify any potential lattice artifacts with the convex-hull method described in the Methods section. Each simulation was performed on a new realization of the VRLs with l0 = 0.5. We also attempted to solve the Kitsunezaki model for some l0<0.5. However, this resulted in numerical instability for with timestep , originating from the spatial derivatives in the diffusion terms. This numerical instability can be mitigated at the cost of increased computation time by choosing a smaller timestep . We randomly shifted the fluid-derived lattice such that the initial disk of the b field was placed at a different position on the fluid-derived lattice each time.
A representative final colony and the normal vector distributions are plotted in Fig 7. We see a pronounced diamond shape in the solutions on the restricted VRL in Fig 7a. This diamond shape is much less pronounced in the colonies resulting on the redrawn VRL (Fig 7b) and hard to recognize in individual instances. However, the histogram does reveal the same four-fold symmetry that we recognize from our results for the agent-based colonies with the diamond-shaped artifact. These results show that some of the structure of the square reference lattice used to generate the VRLs is present in the solutions to the Kitsunezaki model. On the other hand, the solutions on the fluid-derived lattice in Fig 7c show no discrete orientational symmetry.
[Figure omitted. See PDF.]
The Kitsunezaki model [19] was solved on the three disordered lattices with parameters D0 = 0.1 and . The initial condition for the density of active motile bacteria b was a disk of radius 5 with b = 1 and b = 0 everywhere else. The nutrient density field n was initialized with n = 1 everywhere. The minimum distance between lattice sites in the VRLs l0 was set to 0.5. The central figures show a single representative result of where s is the density of inactive bacteria. The outer rings show the distribution of normal vectors to the fitted convex hulls as in Fig 2 based on 500 (independent) realizations.
Next, we investigated how well the minimal hybrid model works on the various disordered lattices. We performed simulations with both high and low c0 and applied the convex-hull procedure for each c0 on all disordered lattices. Representative final colony shapes and the histograms of normal vectors to the fitted convex hulls are given in Fig 8. We see that all three disordered lattices are effective in removing the discrete orientational symmetry observed in Figs 2, 4, and 7.
[Figure omitted. See PDF.]
The top row shows simulations with high initial nutrient concentration (c0 = 3) and the bottom row shows simulations with low initial nutrient concentrations (c0 = 0.7). The minimum separation between lattice sites in the VRLs l0 was set to 0.5. The selection of parameters and representation is otherwise the same as in Fig 2.
A key distinction between the continuous and agent components of the hybrid model lies in how they interact with the underlying lattice. The agent component of the model has a binary neighbor relation: a site is a neighbor or it is not. The continuous component of the model accounts for distances between sites, sizes of Voronoi cell walls, and Voronoi cell areas. This difference between the continuous and agent components of the hybrid model suggests that an agent-based model that includes more aspects of the lattice than the agent component of our hybrid model would be affected by lattice artifacts caused by the VRLs.
To test this hypothesis, we formulated a model of motile bacteria with cell-to-cell adhesion. We assume that the length of the Voronoi cell wall between lattice sites serves as a reasonable proxy for the contact area between bacterial cells, an assumption that is more applicable to the contact area between cells in confluent tissues than to bacterial colonies. However, this assumption allows us to incorporate more of the lattice geometry into the agent component of our hybrid model. The simulation begins with a fixed number of bacteria arranged in a colony, and no cell division is included in the model. The energy H of a configuration of bacteria is determined by the adhesion between bacteria and is given by:
(5)
where Nb is the number of bacteria, J determines the strength of adhesion per unit contact length between bacteria, and Lij is the length of the Voronoi cell wall between the sites occupied by bacteria i and j. Note that Lij = 0 when bacteria i and j are not on neighboring lattice sites. Eq (5) says that the energetically most favorable position for a bacterium is to be surrounded by other bacteria on all sides. We then use a Monte Carlo method described in the Methods section to find the equilibrium configuration of the bacterial colony.
We performed simulations of this model using the three disordered lattices and characterized the resulting equilibrium configurations by fitting convex hulls to the main cluster as detailed in the Methods section. We performed simulations starting with a circular colony of radius r = 30 and starting with a square colony with edge length for each of the lattices and ran the simulations long enough such that the convex hull normal vector histograms resulting from the square and the circular initial conditions became indistinguishable. The combined results of the 1000 simulations on each lattice are given in Fig 9. The lattice-induced four-fold symmetry has reappeared on the restricted VRL in Fig 9a as expected. We do not detect this lattice artifact on either the redrawn VRL or the fluid-derived lattice of Fig 9b and 9c, respectively. This shows that the geometric properties of the lattice are indeed important for the formation of the discrete orientational symmetry.
[Figure omitted. See PDF.]
The adhesion strength per unit contact area between cells was set to J = 4. The central figures (green) show a representative final colony shape. The representation of the outer rings constructed from 1000 repetitions is as in Fig 2.
Numerical performance of our implementation
We have seen that our fluid-based method can overcome lattice-artifact issues. Let us now turn to the performance of our implementation. We performed a series of simulations with different system sizes to understand the scaling of our method. We ran the simulation three times for each system size and recorded the mean system time when the mean nutrient concentration was lower than 0.01 or at least 70% of the lattice sites were occupied. All simulations were run on a workstation equipped with an Intel Core i7-7700 CPU running at a maximum frequency of 4.2 GHz, and 16 GB of system memory. The results are plotted in Fig 10.
[Figure omitted. See PDF.]
The scaling was determined both for simulations with high initial nutrient concentration (blue dots) and for simulations with low initial nutrient concentration (orange dots). The connecting lines are to guide the eye. The exponent of the power law fit is (dashed line). The longest running simulation consisted of 3,034,566 lattice sites.
The hybrid model exhibits power-law scaling with system size, with a scaling exponent of . There is some deviation from power-law scaling at smaller system sizes ( lattice sites), primarily due to the Voronoi tessellation performed at the start of each simulation. The largest system tested, with low initial nutrient concentration, contained 3,034,566 lattice sites with a final colony consisting of 1,639,299 bacteria. This simulation took approximately 12.5 hours to compute and corresponds to a real-world colony about 2 mm in diameter (assuming bacteria that are 1 μm in diameter). Computation time can be further reduced by parallelization, especially since the continuous component of the model is well-suited for this.
In order to see how the performance of the hybrid lattice-based model with fluid-derived lattice compares to that of a more detailed off-lattice agent-based method, we performed benchmark simulations using the iDynoMiCS software package [22]. This software package was developed to simulate bacterial biofilms in a flow cell. For the iDynoMiCS simulations, we used the simulation parameters as in Ref [13]. In the off-lattice agent-based simulations, we set the bulk nutrient concentration of the limiting nutrient to 10−2 g/L to model dense, non-branching colonies. We initialized the simulations by seeding 200 bacteria at the bottom of the rectangular simulation domain. The lateral boundaries of the simulation domain are periodic and the lower boundary is a wall. The upper boundary is connected to the bulk nutrient reservoir. Half the initial bacteria were colored blue and the other half red. Bacteria passed their color on to their daughter cells.
We adapted our hybrid lattice-based model to feature a rectangular simulation domain to better match the environment simulated by iDynoMiCS. The left and right boundaries are periodic and the top and bottom are walls with zero-flux boundary conditions as before. We initialized the hybrid model by seeding 200 bacteria at the bottom of the simulation domain and we set the initial nutrient concentration c0 to 3.0 to model dense, non-branching colony morphologies. The red and blue colors are the same as in the iDynoMiCS simulations. The remaining parameters are as in Table 1. We tracked the computation time as the colonies increased in size. These results of the agent-based and lattice-based simulations can be found in Fig 11.
[Figure omitted. See PDF.]
Computation time of the hybrid lattice-based model (orange dots) is compared with an agent-based off-lattice simulation of a similar system (blue dots). The connecting lines are to guide the eye. Colony size is plotted against computation time during a single simulation. A power law was fitted to both curves (dashed lines). The exponent fitted to the off-lattice curve was and the exponent fitted to the on-lattice curve was . The off-lattice simulation was performed using the iDynoMiCS software package [22] with parameters chosen as in Ref [13] and with bulk limiting nutrient concentration g/L. The parameters used in the hybrid lattice-based simulation are initial nutrient concentration c0 = 3.0 and the remaining parameters as in Table 1. The insets show the final colonies where half the initial bacteria were colored red and the other half blue. Color was inherited from parent to daughter cells. The final configuration of the off-lattice simulation contains 41,697 agents and took approximately 11 hours to compute. The final configuration of our hybrid on-lattice simulation contains 42,015 agents and took approximately 19 seconds to compute.
Our hybrid lattice-based mode for bacterial growth reaches similar colony sizes in significantly (often orders-of-magnitude) less time than the off-lattice agent-based growth model implemented by iDynoMiCS. Our hybrid model grew to 42,015 bacteria in approximately 19 seconds, while the iDynoMiCS simulation reached 41,697 bacteria in about 11 hours. We also compared simulations in the branching regime. However, the two simulation methods resulted in different colony morphologies. We discuss this discrepancy in further detail in S1 Appendix.
There are some caveats to the above comparison of computation time. For example, iDynoMiCS is written in Java, whereas our hybrid lattice-based model is implemented in C++. Furthermore, some of the physics included in iDynoMiCS, like the shoving interaction between bacteria, is absent in the minimal hybrid model. However, the significant difference in computation time suggests it is worthwhile to consider whether a problem of interest can be studied with lattice-based methods instead of a more detailed agent-based model.
Discussion and conclusion
To characterize lattice artifacts in simulations of growing bacterial colonies we have focused on a broad class of hybrid lattice-based models. Previous studies utilizing similar hybrid models used either a square lattice [42,52] or a hexagonal lattice [51]. In our investigation of the minimal hybrid model using a square lattice, we observed a distinct four-fold orientational symmetry in the shape of resulting colonies, resembling patterns seen in the Eden model [33,48]. This prompted us to explore whether disordered lattices could eliminate such artifacts.
We found that all three disordered lattices that we considered successfully resolved the lattice-induced orientational symmetry in the hybrid agent-based model. We had anticipated that some degree of orientational symmetry would remain in the restricted vectorizable random lattice (VRL), based on the following two observations. (i) We saw significant orientational symmetry in solutions to a continuous model on this lattice. (ii) The restricted VRL itself maintains four-fold orientational symmetry. A key distinction between the continuous and agent components of the minimal hybrid model lies in their use of the underlying lattice: the continuous component includes the distance between sites, the area of Voronoi cells, and the length of Voronoi cell walls, while the agent component only requires a binary neighborhood relationship. It is clear that the disordered lattices can effectively establish the appropriate binary neighborhood relationships for the agent component, where the weighted Moore neighborhood on the square lattice failed to do so. Further investigation is needed to determine the exact properties that these neighborhood relations should exhibit to suppress artifacts.
The four-fold orientational symmetry reappeared when we introduced a geometric property of the lattice in the agent component of the model of motile bacteria with cell-to-cell adhesion. From this, we conclude that the impact of symmetries in the underlying lattice strongly depends on the specific model used. Therefore, one should not assume that a specific choice of lattice will suffice, and it is essential to explicitly check for relevant lattice artifacts.
It is important to note that any lattice will impose some structure on numerical modeling results, which may not be suitable for the system of interest. The fluid-derived lattice can be adapted to better suit the characteristics of the model system. For instance, constructing the lattice using simulations of spherocylinders with viscoelastic interactions [25–27] might be more appropriate to modeling a system of rod-shaped bacteria. Another possible approach would be to base the lattice on an experimental realization of the system of interest. In this case, the lattice would be reconstructed from measurements of the structure of the experimental result.
We acknowledge that using a lattice-based model entails simplifying the system and omitting direct mechanical interactions between bacteria. Nevertheless, it may be possible to incorporate these interactions in an effective, approximate manner. For instance, as bacteria grow and increase in size, they exert mechanical pressure on neighboring cells, effectively pushing them apart. One approach to capturing this ‘shifting’ behavior in a lattice-based framework is to implement a rearrangement algorithm. This algorithm begins by moving a bacterium from the interior of the colony to an already occupied neighboring lattice site, creating a vacant lattice site within the colony. The rearrangement algorithm then resolves the resulting excess bacterium by sequentially displacing other bacteria across the lattice, ultimately relocating the excess to the colony’s boundary. To guide this process, the cumulative nutrient consumption of a cell, used as a proxy for local biomass, can inform the redistribution of cells. The resulting biomass density would direct the shifting algorithm, enabling space to be created in the colony’s interior to facilitate continued cell replication. Alternatively, the shifting algorithm could rely on metrics such as the shortest path length to the colony boundary to guide cell movement on the lattice as proposed by Drasdo and Block [35,65]. Similarly, the behavior of anisotropic bacteria could be approximated by assigning each bacterium a director, representing its orientation. These directors can be incorporated into the model’s rules to enable features such as directed motion and alignment with neighboring cells. Devising lattice-based approximations of mechanical interactions between bacteria could be more challenging than implementation of these dynamics in an off-lattice agent-based method. Despite this challenge, the approximate lattice-based method has the potential to substantially enhance computational performance. A detailed implementation and validation of the proposed dynamics are left for future work.
Let us now turn to the efficiency of our approach. One motivation for our study was to simulate large systems, a task traditionally achievable through continuous models. For instance, the continuous model of Schwarcz et al. [20] successfully simulated colonies up to 5 cm in radius (equivalent to approximately 1010 cells [66]), demonstrating the feasibility of modeling large-scale systems. However, the increased availability of observations of individual microbes has significantly advanced our understanding of bacterial behavior at the individual level [67,68]. Incorporating these individual-level properties into continuous models presents a challenge. Furthermore, when simulations are applied to study evolutionary processes, the discrete nature of individual bacteria becomes crucial.
Agent-based methods are more appropriate than continuous models to tackle questions involving inheritance and diversity, as they keep track of individuals. However, simulating the large system sizes needed for such studies with such methods presents significant challenges. This is especially the case when using off-lattice approaches, which have to also propagate the equations of motion. Typically, off-lattice agent-based systems consist of 10,000 to 100,000 agents [25,55,69]. Some studies have managed to simulate larger systems by optimizing the model and not actively simulating every individual agent. For example, Young et al. [24] were interested in the effect of active layer dynamics on the morphology of growing biofilms. They accelerated their simulations by removing inactive cells far behind the growing front. This allowed them to reach system sizes of 105 bacteria. The NUFEB code, developed by Li et al., is notable for enabling the simulation of exceptionally large bacterial populations, with up to 107 individuals [32]. This implementation builds on the parallelization offered by the popular simulation packages LAMMPS and OpenFOAM, but can only simulate such large systems effectively by making use of considerable computational resources.
A method by which large numbers of bacteria can be simulated with genetic inheritance, but without the need for high-performance computing, is lattice-based modeling. Such approaches have been considered for system sizes in excess of 100,000 bacteria [38]. Mukhamadiarov et al. [43] even considered a colony consisting of bacteria. Our approach stands out from these, as we use disordered lattices as opposed to regular lattices, which removes discrete orientational symmetries. We also considered the efficiency of the lattice-based approach compared to a full off-lattice simulation. Our results showed that the fluid-lattice-based hybrid model is orders of magnitude faster than a similar off-lattice simulation implemented by iDynoMiCS. This makes it clear that simulations involving millions of bacteria are accessible in hours on a regular desktop machine.
Our method can be naturally extended to three dimensions by providing the code with a three-dimensional lattice. The rules governing bacterial behavior, as well as the partial differential equations for the continuous fields, can be specified in the same way as in the two-dimensional case. However, this extension comes with an increase in computational cost. The primary reason for this is the growth in the number of lattice sites: when a two-dimensional lattice is extended to three dimensions with height h, the total number of sites increases approximately by a factor of h. In addition, the scaling of computation time with system size becomes slightly less favorable in three dimensions, due to the larger number of neighboring sites per lattice point, which increases the number of calculations required per site. Despite this added complexity, we still expect the three-dimensional lattice-based method to significantly outperform equivalent three-dimensional off-lattice agent-based simulations. This is because the bacterial dynamics in our model remain considerably simplified compared to those in off-lattice approaches. Alternatively, since pattern formation in expanding colonies primarily occurs in the plane, a pseudo-three-dimensional version of our model can be employed. In this approach, the two-dimensional lattice is modified to allow multiple bacteria to occupy a single site, with the number of bacteria per site representing the local colony height. The computational performance of this pseudo-three-dimensional model is expected to be comparable to that of the original two-dimensional version.
In summary, we have studied a minimal hybrid model for bacterial populations that constrains individual bacteria to a lattice and couples them to a continuous field. We were able to simulate millions of individual bacteria in a matter of hours on a regular desktop computer. This is far beyond what off-lattice codes can do, as we have argued above. In addition, we demonstrated that fluid-derived disordered lattices can be effective in mitigating lattice-induced symmetries. In future work, the hybrid model can be further extended to include features like active motility, lubricant production, and anisotropic bacteria. These extensions could help shed light on the formation of chiral colonies [3,4,70] and colonies consisting of concentric rings [1,2]. The larger accessible population sizes and the individual nature of the bacteria in our model will allow us to apply our model to questions relating to evolutionary processes in spatially structured populations. One example of this is selection in spatially structured populations due to a spatially inhomogeneous selection pressure, which could be relevant in the evolution of antimicrobial resistance [71–73]. We hope that the hybrid approach on a disordered lattice will prove to be a useful tool for the community interested in modeling spatially structured populations.
Methods
Definition of the hybrid bacterial-growth model
Here, we will provide details of the hybrid bacterial-growth model used in Figs 2, 3, 4, and 8. The bacteria are modeled as individual agents at positions . The nutrients are modeled as a continuous field c governed by a reaction-diffusion equation:
(6)
The first term describes the diffusion of nutrient with diffusion constant D. The second term describes the consumption of nutrient by all bacteria. The bacteria appear as point sinks, implemented by Dirac delta peaks at the center of each bacterium. The function f(c) describes the concentration-dependent nutrient uptake rate, which we chose to be specified by the Monod equation:
(7)
with the maximum nutrient consumption rate and K the half-saturation constant. The Monod Eq (7) states that, at low nutrient concentrations, the nutrient consumption rate increases with the nutrient concentration. When the nutrient concentration is high, bacteria consume the nutrients at a fixed maximum rate .
The positions of the bacteria are constrained to a two-dimensional (2D) lattice. Each lattice site can be occupied by at most one bacterium. Each bacterium is assigned an internal state n which describes the amount of nutrients consumed. Therefore, the rate of change of n is also given by the Monod Eq (7). When n exceeds a threshold , the bacterium has consumed enough nutrient to produce one daughter cell. This daughter cell is randomly placed on an empty neighbor lattice site. If there are no unoccupied neighbor sites, the bacterium is unable to produce the daughter cell. The daughter cell receives nd nutrients where nd is chosen from the interval with uniform probability. We set for all simulations. The parent cell retains nutrients. At every time step, the order in which these division events are executed is randomized.
Generating a disordered lattice
In the Results section, we evaluated three types of disordered lattice: reduced vectorizable random lattices (VRLs) as proposed by Moukarzel and Herrmann [49], redrawn VRLs as proposed by Tucker [50], and fluid-derived random lattices. We generated the VRLs as described in the Results section (Fig 5). The fluid-derived random lattices were obtained by molecular-dynamics simulations of a dense (2D) fluid. We simulated 40,000 soft disks (bidisperse sample with stoichiometry 1:1 and size ratio 1.25) in a square 2D domain with periodic boundary conditions at area fraction . The interaction between the disks is governed by the WCA potential [63]:
(8)
where ε is the interaction strength, rij the distance between particles i and j, and with di the diameter of particle i. We ran the simulations with where kb is the Boltzmann constant and T the temperature. We used the center of mass coordinates of the disks as pseudo-random coordinates for our lattice. When we needed more than 40,000 lattice sites, we extended the lattice by copying the original 40,000 coordinates until our lattice had the desired size.
From our set of pseudo-random coordinates, we obtained the desired disordered lattice by connecting neighboring sites. Here, neighbors of a lattice site are defined using a Delaunay triangulation. This is constructed from the regular Voronoi tessellation corresponding to our pseudo-random coordinates. The Voronoi tessellation partitions the domain into cells, where cell i is the set of all points that are closer to lattice site i than to any other site. Here, we use the Voro++ library developed by Rycroft [74] to construct the tessellation. By connecting all the sites that share a Voronoi cell wall, we obtain the associated Delaunay triangulation, which will serve as our disordered lattice. Two sites are neighbors if they are directly connected in the Delaunay triangulation. We rescaled the coordinates in the final configurations such that the average distance between sites equals the unit of distance .
Numerical integration of reaction-diffusion equations
The reaction-diffusion Eq (6) is discretized using the same lattice that constrains the bacteria. We always use the forward Euler method for time integration and use circular zero-flux boundary conditions. We use different approaches to calculate the spatial derivatives of Eq (6), that depend on the specifics of the lattice.
For square lattices, we calculate the gradient of the nutrient field using a finite differences scheme. We only consider the von Neumann neighborhood of each lattice site, resulting in a five-point central differences scheme. The discretized version of Eq (6) thus becomes
(9)
where is the area of a grid cell with lattice spacing , the time-step size, D the nutrient diffusion coefficient, the nutrient concentration at at time , and the indicator function that equals 1 when a bacterium is present and 0 otherwise. The boundary conditions are implemented with the method proposed in Ref [75].
For the disordered lattices, we use a finite-volumes approach [76]. The lattice needs to be subdivided into control volumes each containing one lattice site. We use the Voronoi cells we calculated when constructing the lattice for this purpose. Then Eq (6) discretized on the disordered lattice becomes
(10)
with the nutrient concentration on the lattice site labeled with index i at time , Ai the area of the Voronoi cell associated with lattice site i, N number of lattice sites, D the nutrient diffusion coefficient, Lij the length of the Voronoi cell wall separating sites i and j, which equals 0 when i and j are not neighboring sites, dij the distance between sites i and j, and the time-step size. The function is the same as in Eq (9).
Monte Carlo simulations of motile bacteria with cell-to-cell adhesion
The equilibrium configurations presented in Fig 9 were obtained as follows. The model assumes a fixed number of motile bacteria. The interaction between bacteria is given by cell-to-cell adhesion and is captured by the Hamiltonian of Eq (5). We then use the Metropolis-Hastings Monte Carlo method [77,78] to find equilibrium configurations. Trial moves consisted of attempts to move a randomly selected bacterium to a random neighbor site. If the neighbor site is occupied, the trial move is rejected. If not the trial move is accepted with probability
(11)
where o and n indicate the old and new configuration respectively, No and Nn the number of neighbor sites of the selected bacterium in the old and new configurations, and H the Hamiltonian governing the interaction between bacteria. The ratio of No and Nn is necessary to ensure detailed balance in disordered lattices where the number of neighbors is not the same for each site. Each Monte Carlo simulation consisted of trial moves.
Quantifying lattice-induced orientational symmetry
We quantify the observed lattice artifacts by fitting a convex hull to each individual final colony shape. We use 500 samples of each simulation to gather statistics. Each simulation was performed on a different realization of the VRLs and initialized on a different position on the fluid-derived lattice. We use the SciPy package to fit the convex hulls [79]. Next, we determine the normal vectors to each segment of the fitted convex hull. These vectors are then used to construct a histogram of their angles, with each vector weighted by the length of its corresponding convex hull segment. The height of the k-th bin in the histogram, denoted hk, is then given by
(12)
where Wk is the set of all normal vector weights in bin k. The standard deviation of the k-th bin is calculated as
(13)
Finally, the histogram is normalized to yield a probability density. A schematic overview of this procedure is given in S2 Figure.
Supporting information
S1 Text. Differences in colony morphology in simulations with iDynoMiCS and the hybrid lattice-based method.
https://doi.org/10.1371/journal.pone.0330491.s001
S2 Fig. Overview of the convex hull procedure for detecting orientational symmetries in the shape of simulated bacterial colonies.
https://doi.org/10.1371/journal.pone.0330491.s002
Acknowledgments
We are grateful to Hossein Nemati and Kim William Torre for useful discussions.
References
1. 1. Ràfols I. Formation of concentric rings in bacterial colonies [MSc thesis]. Japan: Chuo University; 1998. http://users.sussex.ac.uk/ir28/docs/natsci/Rafols-MSc-Thesis-1998.pdf
2. 2. Wakita J, Itoh H, Matsuyama T, Matsushita M. Self-affinity for the growing interface of bacterial colonies. J Phys Soc Jpn. 1997;66(1):67–72.
* View Article
* Google Scholar
3. 3. Ben-Jacob E, Tenenbaum A, Shochet O, Avidan O. Holotransformations of bacterial colonies and genome cybernetics. Physica A: Statistical Mechanics and its Applications. 1994;202(1–2):1–47.
* View Article
* Google Scholar
4. 4. Ben-Jacob E, Cohen I I, Shochet O, Tenenbaum A, Czirók A, Vicsek T. Cooperative formation of chiral patterns during growth of bacterial colonies. Phys Rev Lett. 1995;75(15):2899–902. pmid:10059433
* View Article
* PubMed/NCBI
* Google Scholar
5. 5. Brown MR, Allison DG, Gilbert P. Resistance of bacterial biofilms to antibiotics: a growth-rate related effect?. J Antimicrob Chemother. 1988;22(6):777–80. pmid:3072331
* View Article
* PubMed/NCBI
* Google Scholar
6. 6. Ashby MJ, Neale JE, Knott SJ, Critchley IA. Effect of antibiotics on non-growing planktonic cells and biofilms of Escherichia coli. J Antimicrob Chemother. 1994;33(3):443–52. pmid:8040110
* View Article
* PubMed/NCBI
* Google Scholar
7. 7. Xu KD, Stewart PS, Xia F, Huang CT, McFeters GA. Spatial physiological heterogeneity in Pseudomonas aeruginosa biofilm is determined by oxygen availability. Appl Environ Microbiol. 1998;64(10):4035–9. pmid:9758837
* View Article
* PubMed/NCBI
* Google Scholar
8. 8. Ito A, Taniuchi A, May T, Kawata K, Okabe S. Increased antibiotic resistance of Escherichia coli in mature biofilms. Appl Environ Microbiol. 2009;75(12):4093–100. pmid:19376922
* View Article
* PubMed/NCBI
* Google Scholar
9. 9. Frost I, Smith WPJ, Mitri S, Millan AS, Davit Y, Osborne JM, et al. Cooperation, competition and antibiotic resistance in bacterial colonies. ISME J. 2018;12(6):1582–93. pmid:29563570
* View Article
* PubMed/NCBI
* Google Scholar
10. 10. Nadell CD, Drescher K, Wingreen NS, Bassler BL. Extracellular matrix structure governs invasion resistance in bacterial biofilms. ISME J. 2015;9(8):1700–9. pmid:25603396
* View Article
* PubMed/NCBI
* Google Scholar
11. 11. Matz C, McDougald D, Moreno AM, Yung PY, Yildiz FH, Kjelleberg S. Biofilm formation and phenotypic variation enhance predation-driven persistence of Vibrio cholerae. Proc Natl Acad Sci U S A. 2005;102(46):16819–24. pmid:16267135
* View Article
* PubMed/NCBI
* Google Scholar
12. 12. Hallatschek O, Hersen P, Ramanathan S, Nelson DR. Genetic drift at expanding frontiers promotes gene segregation. Proc Natl Acad Sci U S A. 2007;104(50):19926–30. pmid:18056799
* View Article
* PubMed/NCBI
* Google Scholar
13. 13. Young E, Allen RJ. Lineage dynamics in growing biofilms: spatial patterns of standing vs. de novo diversity. Front Microbiol. 2022;13:915095. pmid:35966660
* View Article
* PubMed/NCBI
* Google Scholar
14. 14. Deforet M, Carmona-Fontaine C, Korolev KS, Xavier JB. Evolution at the edge of expanding populations. Am Nat. 2019;194(3):291–305. pmid:31553215
* View Article
* PubMed/NCBI
* Google Scholar
15. 15. Zöllner R, Oldewurtel ER, Kouzel N, Maier B. Phase and antigenic variation govern competition dynamics through positioning in bacterial colonies. Sci Rep. 2017;7(1):12151. pmid:28939833
* View Article
* PubMed/NCBI
* Google Scholar
16. 16. Golding I, Kozlovsky Y, Cohen I, Ben-Jacob E. Studies of bacterial branching growth using reaction–diffusion models for colonial development. Physica A: Statistical Mechanics and its Applications. 1998;260(3–4):510–54.
* View Article
* Google Scholar
17. 17. Kozlovsky Y, Cohen I, Golding I, Ben-Jacob E. Lubricating bacteria model for branching growth of bacterial colonies. Phys Rev E. 1999;59(6):7025–35.
* View Article
* Google Scholar
18. 18. Matsushita M, Hiramatsu F, Kobayashi N, Ozawa T, Yamazaki Y, Matsuyama T. Colony formation in bacteria: experiments and modeling. Biofilms. 2004;1(4):305–17.
* View Article
* Google Scholar
19. 19. Kitsunezaki S. Interface dynamics for bacterial colony formation. J Phys Soc Jpn. 1997;66(5):1544–50.
* View Article
* Google Scholar
20. 20. Schwarcz D, Levine H, Ben-Jacob E, Ariel G. Uniform modeling of bacterial colony patterns with varying nutrient and substrate. Physica D: Nonlinear Phenomena. 2016;318–319:91–9.
* View Article
* Google Scholar
21. 21. Nadell CD, Foster KR, Xavier JB. Emergence of spatial structure in cell groups and the evolution of cooperation. PLoS Comput Biol. 2010;6(3):e1000716. pmid:20333237
* View Article
* PubMed/NCBI
* Google Scholar
22. 22. Lardon LA, Merkey BV, Martins S, Dötsch A, Picioreanu C, Kreft J-U, et al. iDynoMiCS: next-generation individual-based modelling of biofilms. Environ Microbiol. 2011;13(9):2416–34. pmid:21410622
* View Article
* PubMed/NCBI
* Google Scholar
23. 23. Bonachela JA, Nadell CD, Xavier JB, Levin SA. Universality in Bacterial Colonies. J Stat Phys. 2011;144(2):303–15.
* View Article
* Google Scholar
24. 24. Young E, Melaugh G, Allen RJ. Active layer dynamics drives a transition to biofilm fingering. NPJ Biofilms Microbiomes. 2023;9(1):17. pmid:37024470
* View Article
* PubMed/NCBI
* Google Scholar
25. 25. Farrell FDC, Hallatschek O, Marenduzzo D, Waclaw B. Mechanically driven growth of quasi-two-dimensional microbial colonies. Phys Rev Lett. 2013;111(16):168101. pmid:24182305
* View Article
* PubMed/NCBI
* Google Scholar
26. 26. Warren MR, Sun H, Yan Y, Cremer J, Li B, Hwa T. Spatiotemporal establishment of dense bacterial colonies growing on hard agar. Elife. 2019;8:e41093. pmid:30855227
* View Article
* PubMed/NCBI
* Google Scholar
27. 27. Volfson D, Cookson S, Hasty J, Tsimring LS. Biomechanical ordering of dense cell populations. Proc Natl Acad Sci U S A. 2008;105(40):15346–51. pmid:18832176
* View Article
* PubMed/NCBI
* Google Scholar
28. 28. Langeslay B, Juarez G. Microdomains and stress distributions in bacterial monolayers on curved interfaces. Soft Matter. 2023;19(20):3605–13. pmid:37161525
* View Article
* PubMed/NCBI
* Google Scholar
29. 29. Los R, Echten D v H t, Nordemann G, Wehrens M, Tans SJ, Idema T. Defect dynamics in growing bacterial colonies. arXiv preprint 2022. http://arxiv.org/abs/2003.10509
* View Article
* Google Scholar
30. 30. Ghosh P, Mondal J, Ben-Jacob E, Levine H. Mechanically-driven phase separation in a growing bacterial colony. Proc Natl Acad Sci U S A. 2015;112(17):E2166-73. pmid:25870260
* View Article
* PubMed/NCBI
* Google Scholar
31. 31. Grant MAA, Wacław B, Allen RJ, Cicuta P. The role of mechanical forces in the planar-to-bulk transition in growing Escherichia coli microcolonies. J R Soc Interface. 2014;11(97):20140400. pmid:24920113
* View Article
* PubMed/NCBI
* Google Scholar
32. 32. Li B, Taniguchi D, Gedara JP, Gogulancea V, Gonzalez-Cabaleiro R, Chen J, et al. NUFEB: a massively parallel simulator for individual-based modelling of microbial communities. PLoS Comput Biol. 2019;15(12):e1007125. pmid:31830032
* View Article
* PubMed/NCBI
* Google Scholar
33. 33. Eden M. A two-dimensional growth process. Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, vol. 4. 1961. p. 223–39.
34. 34. Ermentrout GB, Edelstein-Keshet L. Cellular automata approaches to biological modeling. J Theor Biol. 1993;160(1):97–133. pmid:8474249
* View Article
* PubMed/NCBI
* Google Scholar
35. 35. Drasdo D. Coarse graining in simulated cell populations. Advs Complex Syst. 2005;08(02n03):319–63.
* View Article
* Google Scholar
36. 36. Kusch I, Markus M. Mollusc shell pigmentation: cellular automaton simulations and evidence for undecidability. Journal of Theoretical Biology. 1996;178(3):333–40.
* View Article
* Google Scholar
37. 37. Edmonds CA, Lillie AS, Cavalli-Sforza LL. Mutations arising in the wave front of an expanding population. Proc Natl Acad Sci U S A. 2004;101(4):975–9. pmid:14732681
* View Article
* PubMed/NCBI
* Google Scholar
38. 38. Möbius W, Murray AW, Nelson DR. How obstacles perturb population fronts and alter their genetic structure. PLoS Comput Biol. 2015;11(12):e1004615. pmid:26696601
* View Article
* PubMed/NCBI
* Google Scholar
39. 39. Li LH, Kardar M. Specialization at an expanding front. Phys Rev E. 2023;108(3):L032402. pmid:37849149
* View Article
* PubMed/NCBI
* Google Scholar
40. 40. Korolev KS, Avlund M, Hallatschek O, Nelson DR. Genetic demixing and evolution in linear stepping stone models. Rev Mod Phys. 2010;82(2):1691–718. pmid:21072144
* View Article
* PubMed/NCBI
* Google Scholar
41. 41. Kuhr J-T, Leisner M, Frey E. Range expansion with mutation and selection: dynamical phase transition in a two-species Eden model. New J Phys. 2011;13(11):113013.
* View Article
* Google Scholar
42. 42. van Dijk B, Hogeweg P. In silico gene-level evolution explains microbial population diversity through differential gene mobility. Genome Biol Evol. 2015;8(1):176–88. pmid:26710854
* View Article
* PubMed/NCBI
* Google Scholar
43. 43. Mukhamadiarov RI, Ciarchi M, Olmeda F, Rulands S. Clonal dynamics of surface-driven growing tissues. Phys Rev E. 2024;109(6–1):064407. pmid:39021023
* View Article
* PubMed/NCBI
* Google Scholar
44. 44. Abrudan MI, Smakman F, Grimbergen AJ, Westhoff S, Miller EL, van Wezel GP, et al. Socially mediated induction and suppression of antibiosis during bacterial coexistence. Proc Natl Acad Sci U S A. 2015;112(35):11054–9. pmid:26216986
* View Article
* PubMed/NCBI
* Google Scholar
45. 45. Swat MH, Thomas GL, Belmonte JM, Shirinifard A, Hmeljak D, Glazier JA. Multi-scale modeling of tissues using CompuCell3D. Methods Cell Biol. 2012;110:325–66. pmid:22482955
* View Article
* PubMed/NCBI
* Google Scholar
46. 46. Dukovski I, Bajić D, Chacón JM, Quintin M, Vila JCC, Sulheim S, et al. A metabolic modeling platform for the computation of microbial ecosystems in time and space (COMETS). Nat Protoc. 2021;16(11):5030–82. pmid:34635859
* View Article
* PubMed/NCBI
* Google Scholar
47. 47. Nemati H, de Graaf J. The cellular Potts model on disordered lattices. Soft Matter. 2024;20(42):8337–52. pmid:39283268
* View Article
* PubMed/NCBI
* Google Scholar
48. 48. Batchelor MT, Henry BI. Limits to Eden growth in two and three dimensions. Physics Letters A. 1991;157(4–5):229–36.
* View Article
* Google Scholar
49. 49. Moukarzel C, Herrmann HJ. A vectorizable random lattice. J Stat Phys. 1992;68(5–6):911–23.
* View Article
* Google Scholar
50. 50. Tucker LJ. A new computational approach to simulate pattern formation in Paenibacillus dendritiformis bacterial colonies. UC San Diego; 2010. https://escholarship.org/uc/item/6x91c78r
51. 51. Borer B, Ciccarese D, Johnson D, Or D. Spatial organization in microbial range expansion emerges from trophic dependencies and successful lineages. Commun Biol. 2020;3(1):685. pmid:33208809
* View Article
* PubMed/NCBI
* Google Scholar
52. 52. Gerlee P, Anderson ARA. A hybrid cellular automaton model of clonal evolution in cancer: the emergence of the glycolytic phenotype. J Theor Biol. 2008;250(4):705–22. pmid:18068192
* View Article
* PubMed/NCBI
* Google Scholar
53. 53. van Hoek MJA, Hogeweg P. In silico evolved lac operons exhibit bistability for artificial inducers, but not for lactose. Biophys J. 2006;91(8):2833–43. pmid:16877514
* View Article
* PubMed/NCBI
* Google Scholar
54. 54. Colizzi ES, van Dijk B, Merks RMH, Rozen DE, Vroomans RMA. Evolution of genome fragility enables microbial division of labor. Cold Spring Harbor Laboratory; 2021. https://doi.org/10.1101/2021.06.04.447040
55. 55. Rana N, Ghosh P, Perlekar P. Spreading of nonmotile bacteria on a hard agar plate: Comparison between agent-based and stochastic simulations. Phys Rev E. 2017;96(5–1):052403. pmid:29347735
* View Article
* PubMed/NCBI
* Google Scholar
56. 56. Kawasaki K, Mochizuki A, Matsushita M, Umeda T, Shigesada N. Modeling spatio-temporal patterns generated by Bacillus subtilis. J Theor Biol. 1997;188(2):177–85. pmid:9379672
* View Article
* PubMed/NCBI
* Google Scholar
57. 57. Matsushita M, Wakita J, Itoh H, Watanabe K, Arai T, Matsuyama T, et al. Formation of colony patterns by a bacterial cell population. Physica A: Statistical Mechanics and its Applications. 1999;274(1–2):190–9.
* View Article
* Google Scholar
58. 58. Freche P, Stauffer D, Stanley HE. Surface structure and anisotropy of Eden clusters. J Phys A: Math Gen. 1985;18(18):L1163–8.
* View Article
* Google Scholar
59. 59. Chan TF. Stability analysis of finite difference schemes for the advection-diffusion equation. SIAM J Numer Anal. 1984;21(2):272–84.
* View Article
* Google Scholar
60. 60. Duchemin L, Eggers J. The Explicit–Implicit–Null method: Removing the numerical instability of PDEs. Journal of Computational Physics. 2014;263:37–52.
* View Article
* Google Scholar
61. 61. Christ NH, Friedberg R, Lee TD. Random lattice field theory: general formulation. Nuclear Physics B. 1982;202(1):89–125.
* View Article
* Google Scholar
62. 62. Tanaka H, Tong H, Shi R, Russo J. Revealing key structural features hidden in liquids and glasses. Nat Rev Phys. 2019;1(5):333–48.
* View Article
* Google Scholar
63. 63. Weeks JD, Chandler D, Andersen HC. Role of repulsive forces in determining the equilibrium structure of simple liquids. The Journal of Chemical Physics. 1971;54(12):5237–47.
* View Article
* Google Scholar
64. 64. Watanabe H, Yukawa S, Ito N. Size-dispersity effects in two-dimensional melting. Phys Rev E Stat Nonlin Soft Matter Phys. 2005;71(1 Pt 2):016702. pmid:15697759
* View Article
* PubMed/NCBI
* Google Scholar
65. 65. Block M, Schöll E, Drasdo D. Classifying the expansion kinetics and critical surface dynamics of growing cell populations. Phys Rev Lett. 2007;99(24):248101. pmid:18233492
* View Article
* PubMed/NCBI
* Google Scholar
66. 66. Nelson DE, Young KD. Penicillin binding protein 5 affects cell diameter, contour, and morphology of Escherichia coli. J Bacteriol. 2000;182(6):1714–21. pmid:10692378
* View Article
* PubMed/NCBI
* Google Scholar
67. 67. Stewart PS, Franklin MJ. Physiological heterogeneity in biofilms. Nat Rev Microbiol. 2008;6(3):199–210. pmid:18264116
* View Article
* PubMed/NCBI
* Google Scholar
68. 68. Hellweger FL, Bucci V. A bunch of tiny individuals—individual-based modeling for microbes. Ecological Modelling. 2009;220(1):8–22.
* View Article
* Google Scholar
69. 69. Los R, Fecker T, Touw PAM, Tatenhove-Pel RJV, Idema T. Time of first contact determines cooperator success in a cross-feeding consortium. 2024. https://www.biorxiv.org/content/10.1101/2024.05.13.593921v2
70. 70. Cohen I, Ben-Jacob E. Orientation field model for chiral branching growth of bacterial colonies. 2000. http://arxiv.org/abs/cond-mat/0008446
71. 71. Greulich P, Waclaw B, Allen RJ. Mutational pathway determines whether drug gradients accelerate evolution of drug-resistant cells. Phys Rev Lett. 2012;109(8):088101. pmid:23002776
* View Article
* PubMed/NCBI
* Google Scholar
72. 72. Baym M, Lieberman TD, Kelsic ED, Chait R, Gross R, Yelin I, et al. Spatiotemporal microbial evolution on antibiotic landscapes. Science. 2016;353(6304):1147–51. pmid:27609891
* View Article
* PubMed/NCBI
* Google Scholar
73. 73. Hermsen R. The adaptation rate of a quantitative trait in an environmental gradient. Phys Biol. 2016;13(6):065003. pmid:27902491
* View Article
* PubMed/NCBI
* Google Scholar
74. 74. Rycroft CH. VORO++: a three-dimensional voronoi cell library in C++. Chaos. 2009;19(4):041111. pmid:20059195
* View Article
* PubMed/NCBI
* Google Scholar
75. 75. Hunt B. Finite difference approximation of boundary conditions along irregular boundaries. Numerical Meth Engineering. 1978;12(2):229–35.
* View Article
* Google Scholar
76. 76. Barth T, Herbin R, Ohlberger M. Finite volume methods: foundation and analysis. encyclopedia of computational mechanics second edition. Wiley; 2017. p. 1–60. https://doi.org/10.1002/9781119176817.ecm2010
77. 77. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. The Journal of Chemical Physics. 1953;21(6):1087–92.
* View Article
* Google Scholar
78. 78. Hastings WK. Monte Carlo sampling methods using Markov chains and their applications. Biometrika. 1970;57(1):97–109.
* View Article
* Google Scholar
79. 79. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods. 2020;17(3):261–72. pmid:32015543
* View Article
* PubMed/NCBI
* Google Scholar
Citation: Verhoef B, Hermsen R, de Graaf J (2025) Fluid-derived lattices for unbiased modeling of bacterial colony growth. PLoS One 20(8): e0330491. https://doi.org/10.1371/journal.pone.0330491
About the Authors:
Bryan Verhoef
Roles: Formal analysis, Methodology, Software, Validation, Visualization, Writing – original draft
E-mail: [email protected]
Affiliations: Institute for Theoretical Physics, Utrecht University, Princetonplein, Utrecht, The Netherlands, Theoretical Biology Group, Department of Biology, Utrecht University, Padualaan, Utrecht, The Netherlands
ORICD: https://orcid.org/0009-0000-3625-4596
Rutger Hermsen
Roles: Methodology, Supervision, Writing – review & editing
Affiliations: Theoretical Biology Group, Department of Biology, Utrecht University, Padualaan, Utrecht, The Netherlands, Centre for Complex Systems Studies, Utrecht University, Leuvenlaan, Utrecht, The Netherlands
Joost de Graaf
Roles: Conceptualization, Methodology, Supervision, Writing – review & editing
Affiliations: Institute for Theoretical Physics, Utrecht University, Princetonplein, Utrecht, The Netherlands, Centre for Complex Systems Studies, Utrecht University, Leuvenlaan, Utrecht, The Netherlands
[/RAW_REF_TEXT]
[/RAW_REF_TEXT]
[/RAW_REF_TEXT]
[/RAW_REF_TEXT]
[/RAW_REF_TEXT]
[/RAW_REF_TEXT]
[/RAW_REF_TEXT]
1. Ràfols I. Formation of concentric rings in bacterial colonies [MSc thesis]. Japan: Chuo University; 1998. http://users.sussex.ac.uk/ir28/docs/natsci/Rafols-MSc-Thesis-1998.pdf
2. Wakita J, Itoh H, Matsuyama T, Matsushita M. Self-affinity for the growing interface of bacterial colonies. J Phys Soc Jpn. 1997;66(1):67–72.
3. Ben-Jacob E, Tenenbaum A, Shochet O, Avidan O. Holotransformations of bacterial colonies and genome cybernetics. Physica A: Statistical Mechanics and its Applications. 1994;202(1–2):1–47.
4. Ben-Jacob E, Cohen I I, Shochet O, Tenenbaum A, Czirók A, Vicsek T. Cooperative formation of chiral patterns during growth of bacterial colonies. Phys Rev Lett. 1995;75(15):2899–902. pmid:10059433
5. Brown MR, Allison DG, Gilbert P. Resistance of bacterial biofilms to antibiotics: a growth-rate related effect?. J Antimicrob Chemother. 1988;22(6):777–80. pmid:3072331
6. Ashby MJ, Neale JE, Knott SJ, Critchley IA. Effect of antibiotics on non-growing planktonic cells and biofilms of Escherichia coli. J Antimicrob Chemother. 1994;33(3):443–52. pmid:8040110
7. Xu KD, Stewart PS, Xia F, Huang CT, McFeters GA. Spatial physiological heterogeneity in Pseudomonas aeruginosa biofilm is determined by oxygen availability. Appl Environ Microbiol. 1998;64(10):4035–9. pmid:9758837
8. Ito A, Taniuchi A, May T, Kawata K, Okabe S. Increased antibiotic resistance of Escherichia coli in mature biofilms. Appl Environ Microbiol. 2009;75(12):4093–100. pmid:19376922
9. Frost I, Smith WPJ, Mitri S, Millan AS, Davit Y, Osborne JM, et al. Cooperation, competition and antibiotic resistance in bacterial colonies. ISME J. 2018;12(6):1582–93. pmid:29563570
10. Nadell CD, Drescher K, Wingreen NS, Bassler BL. Extracellular matrix structure governs invasion resistance in bacterial biofilms. ISME J. 2015;9(8):1700–9. pmid:25603396
11. Matz C, McDougald D, Moreno AM, Yung PY, Yildiz FH, Kjelleberg S. Biofilm formation and phenotypic variation enhance predation-driven persistence of Vibrio cholerae. Proc Natl Acad Sci U S A. 2005;102(46):16819–24. pmid:16267135
12. Hallatschek O, Hersen P, Ramanathan S, Nelson DR. Genetic drift at expanding frontiers promotes gene segregation. Proc Natl Acad Sci U S A. 2007;104(50):19926–30. pmid:18056799
13. Young E, Allen RJ. Lineage dynamics in growing biofilms: spatial patterns of standing vs. de novo diversity. Front Microbiol. 2022;13:915095. pmid:35966660
14. Deforet M, Carmona-Fontaine C, Korolev KS, Xavier JB. Evolution at the edge of expanding populations. Am Nat. 2019;194(3):291–305. pmid:31553215
15. Zöllner R, Oldewurtel ER, Kouzel N, Maier B. Phase and antigenic variation govern competition dynamics through positioning in bacterial colonies. Sci Rep. 2017;7(1):12151. pmid:28939833
16. Golding I, Kozlovsky Y, Cohen I, Ben-Jacob E. Studies of bacterial branching growth using reaction–diffusion models for colonial development. Physica A: Statistical Mechanics and its Applications. 1998;260(3–4):510–54.
17. Kozlovsky Y, Cohen I, Golding I, Ben-Jacob E. Lubricating bacteria model for branching growth of bacterial colonies. Phys Rev E. 1999;59(6):7025–35.
18. Matsushita M, Hiramatsu F, Kobayashi N, Ozawa T, Yamazaki Y, Matsuyama T. Colony formation in bacteria: experiments and modeling. Biofilms. 2004;1(4):305–17.
19. Kitsunezaki S. Interface dynamics for bacterial colony formation. J Phys Soc Jpn. 1997;66(5):1544–50.
20. Schwarcz D, Levine H, Ben-Jacob E, Ariel G. Uniform modeling of bacterial colony patterns with varying nutrient and substrate. Physica D: Nonlinear Phenomena. 2016;318–319:91–9.
21. Nadell CD, Foster KR, Xavier JB. Emergence of spatial structure in cell groups and the evolution of cooperation. PLoS Comput Biol. 2010;6(3):e1000716. pmid:20333237
22. Lardon LA, Merkey BV, Martins S, Dötsch A, Picioreanu C, Kreft J-U, et al. iDynoMiCS: next-generation individual-based modelling of biofilms. Environ Microbiol. 2011;13(9):2416–34. pmid:21410622
23. Bonachela JA, Nadell CD, Xavier JB, Levin SA. Universality in Bacterial Colonies. J Stat Phys. 2011;144(2):303–15.
24. Young E, Melaugh G, Allen RJ. Active layer dynamics drives a transition to biofilm fingering. NPJ Biofilms Microbiomes. 2023;9(1):17. pmid:37024470
25. Farrell FDC, Hallatschek O, Marenduzzo D, Waclaw B. Mechanically driven growth of quasi-two-dimensional microbial colonies. Phys Rev Lett. 2013;111(16):168101. pmid:24182305
26. Warren MR, Sun H, Yan Y, Cremer J, Li B, Hwa T. Spatiotemporal establishment of dense bacterial colonies growing on hard agar. Elife. 2019;8:e41093. pmid:30855227
27. Volfson D, Cookson S, Hasty J, Tsimring LS. Biomechanical ordering of dense cell populations. Proc Natl Acad Sci U S A. 2008;105(40):15346–51. pmid:18832176
28. Langeslay B, Juarez G. Microdomains and stress distributions in bacterial monolayers on curved interfaces. Soft Matter. 2023;19(20):3605–13. pmid:37161525
29. Los R, Echten D v H t, Nordemann G, Wehrens M, Tans SJ, Idema T. Defect dynamics in growing bacterial colonies. arXiv preprint 2022. http://arxiv.org/abs/2003.10509
30. Ghosh P, Mondal J, Ben-Jacob E, Levine H. Mechanically-driven phase separation in a growing bacterial colony. Proc Natl Acad Sci U S A. 2015;112(17):E2166-73. pmid:25870260
31. Grant MAA, Wacław B, Allen RJ, Cicuta P. The role of mechanical forces in the planar-to-bulk transition in growing Escherichia coli microcolonies. J R Soc Interface. 2014;11(97):20140400. pmid:24920113
32. Li B, Taniguchi D, Gedara JP, Gogulancea V, Gonzalez-Cabaleiro R, Chen J, et al. NUFEB: a massively parallel simulator for individual-based modelling of microbial communities. PLoS Comput Biol. 2019;15(12):e1007125. pmid:31830032
33. Eden M. A two-dimensional growth process. Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, vol. 4. 1961. p. 223–39.
34. Ermentrout GB, Edelstein-Keshet L. Cellular automata approaches to biological modeling. J Theor Biol. 1993;160(1):97–133. pmid:8474249
35. Drasdo D. Coarse graining in simulated cell populations. Advs Complex Syst. 2005;08(02n03):319–63.
36. Kusch I, Markus M. Mollusc shell pigmentation: cellular automaton simulations and evidence for undecidability. Journal of Theoretical Biology. 1996;178(3):333–40.
37. Edmonds CA, Lillie AS, Cavalli-Sforza LL. Mutations arising in the wave front of an expanding population. Proc Natl Acad Sci U S A. 2004;101(4):975–9. pmid:14732681
38. Möbius W, Murray AW, Nelson DR. How obstacles perturb population fronts and alter their genetic structure. PLoS Comput Biol. 2015;11(12):e1004615. pmid:26696601
39. Li LH, Kardar M. Specialization at an expanding front. Phys Rev E. 2023;108(3):L032402. pmid:37849149
40. Korolev KS, Avlund M, Hallatschek O, Nelson DR. Genetic demixing and evolution in linear stepping stone models. Rev Mod Phys. 2010;82(2):1691–718. pmid:21072144
41. Kuhr J-T, Leisner M, Frey E. Range expansion with mutation and selection: dynamical phase transition in a two-species Eden model. New J Phys. 2011;13(11):113013.
42. van Dijk B, Hogeweg P. In silico gene-level evolution explains microbial population diversity through differential gene mobility. Genome Biol Evol. 2015;8(1):176–88. pmid:26710854
43. Mukhamadiarov RI, Ciarchi M, Olmeda F, Rulands S. Clonal dynamics of surface-driven growing tissues. Phys Rev E. 2024;109(6–1):064407. pmid:39021023
44. Abrudan MI, Smakman F, Grimbergen AJ, Westhoff S, Miller EL, van Wezel GP, et al. Socially mediated induction and suppression of antibiosis during bacterial coexistence. Proc Natl Acad Sci U S A. 2015;112(35):11054–9. pmid:26216986
45. Swat MH, Thomas GL, Belmonte JM, Shirinifard A, Hmeljak D, Glazier JA. Multi-scale modeling of tissues using CompuCell3D. Methods Cell Biol. 2012;110:325–66. pmid:22482955
46. Dukovski I, Bajić D, Chacón JM, Quintin M, Vila JCC, Sulheim S, et al. A metabolic modeling platform for the computation of microbial ecosystems in time and space (COMETS). Nat Protoc. 2021;16(11):5030–82. pmid:34635859
47. Nemati H, de Graaf J. The cellular Potts model on disordered lattices. Soft Matter. 2024;20(42):8337–52. pmid:39283268
48. Batchelor MT, Henry BI. Limits to Eden growth in two and three dimensions. Physics Letters A. 1991;157(4–5):229–36.
49. Moukarzel C, Herrmann HJ. A vectorizable random lattice. J Stat Phys. 1992;68(5–6):911–23.
50. Tucker LJ. A new computational approach to simulate pattern formation in Paenibacillus dendritiformis bacterial colonies. UC San Diego; 2010. https://escholarship.org/uc/item/6x91c78r
51. Borer B, Ciccarese D, Johnson D, Or D. Spatial organization in microbial range expansion emerges from trophic dependencies and successful lineages. Commun Biol. 2020;3(1):685. pmid:33208809
52. Gerlee P, Anderson ARA. A hybrid cellular automaton model of clonal evolution in cancer: the emergence of the glycolytic phenotype. J Theor Biol. 2008;250(4):705–22. pmid:18068192
53. van Hoek MJA, Hogeweg P. In silico evolved lac operons exhibit bistability for artificial inducers, but not for lactose. Biophys J. 2006;91(8):2833–43. pmid:16877514
54. Colizzi ES, van Dijk B, Merks RMH, Rozen DE, Vroomans RMA. Evolution of genome fragility enables microbial division of labor. Cold Spring Harbor Laboratory; 2021. https://doi.org/10.1101/2021.06.04.447040
55. Rana N, Ghosh P, Perlekar P. Spreading of nonmotile bacteria on a hard agar plate: Comparison between agent-based and stochastic simulations. Phys Rev E. 2017;96(5–1):052403. pmid:29347735
56. Kawasaki K, Mochizuki A, Matsushita M, Umeda T, Shigesada N. Modeling spatio-temporal patterns generated by Bacillus subtilis. J Theor Biol. 1997;188(2):177–85. pmid:9379672
57. Matsushita M, Wakita J, Itoh H, Watanabe K, Arai T, Matsuyama T, et al. Formation of colony patterns by a bacterial cell population. Physica A: Statistical Mechanics and its Applications. 1999;274(1–2):190–9.
58. Freche P, Stauffer D, Stanley HE. Surface structure and anisotropy of Eden clusters. J Phys A: Math Gen. 1985;18(18):L1163–8.
59. Chan TF. Stability analysis of finite difference schemes for the advection-diffusion equation. SIAM J Numer Anal. 1984;21(2):272–84.
60. Duchemin L, Eggers J. The Explicit–Implicit–Null method: Removing the numerical instability of PDEs. Journal of Computational Physics. 2014;263:37–52.
61. Christ NH, Friedberg R, Lee TD. Random lattice field theory: general formulation. Nuclear Physics B. 1982;202(1):89–125.
62. Tanaka H, Tong H, Shi R, Russo J. Revealing key structural features hidden in liquids and glasses. Nat Rev Phys. 2019;1(5):333–48.
63. Weeks JD, Chandler D, Andersen HC. Role of repulsive forces in determining the equilibrium structure of simple liquids. The Journal of Chemical Physics. 1971;54(12):5237–47.
64. Watanabe H, Yukawa S, Ito N. Size-dispersity effects in two-dimensional melting. Phys Rev E Stat Nonlin Soft Matter Phys. 2005;71(1 Pt 2):016702. pmid:15697759
65. Block M, Schöll E, Drasdo D. Classifying the expansion kinetics and critical surface dynamics of growing cell populations. Phys Rev Lett. 2007;99(24):248101. pmid:18233492
66. Nelson DE, Young KD. Penicillin binding protein 5 affects cell diameter, contour, and morphology of Escherichia coli. J Bacteriol. 2000;182(6):1714–21. pmid:10692378
67. Stewart PS, Franklin MJ. Physiological heterogeneity in biofilms. Nat Rev Microbiol. 2008;6(3):199–210. pmid:18264116
68. Hellweger FL, Bucci V. A bunch of tiny individuals—individual-based modeling for microbes. Ecological Modelling. 2009;220(1):8–22.
69. Los R, Fecker T, Touw PAM, Tatenhove-Pel RJV, Idema T. Time of first contact determines cooperator success in a cross-feeding consortium. 2024. https://www.biorxiv.org/content/10.1101/2024.05.13.593921v2
70. Cohen I, Ben-Jacob E. Orientation field model for chiral branching growth of bacterial colonies. 2000. http://arxiv.org/abs/cond-mat/0008446
71. Greulich P, Waclaw B, Allen RJ. Mutational pathway determines whether drug gradients accelerate evolution of drug-resistant cells. Phys Rev Lett. 2012;109(8):088101. pmid:23002776
72. Baym M, Lieberman TD, Kelsic ED, Chait R, Gross R, Yelin I, et al. Spatiotemporal microbial evolution on antibiotic landscapes. Science. 2016;353(6304):1147–51. pmid:27609891
73. Hermsen R. The adaptation rate of a quantitative trait in an environmental gradient. Phys Biol. 2016;13(6):065003. pmid:27902491
74. Rycroft CH. VORO++: a three-dimensional voronoi cell library in C++. Chaos. 2009;19(4):041111. pmid:20059195
75. Hunt B. Finite difference approximation of boundary conditions along irregular boundaries. Numerical Meth Engineering. 1978;12(2):229–35.
76. Barth T, Herbin R, Ohlberger M. Finite volume methods: foundation and analysis. encyclopedia of computational mechanics second edition. Wiley; 2017. p. 1–60. https://doi.org/10.1002/9781119176817.ecm2010
77. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. The Journal of Chemical Physics. 1953;21(6):1087–92.
78. Hastings WK. Monte Carlo sampling methods using Markov chains and their applications. Biometrika. 1970;57(1):97–109.
79. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods. 2020;17(3):261–72. pmid:32015543
© 2025 Verhoef et al. This is an open access article distributed under the terms of the Creative Commons Attribution License: http://creativecommons.org/licenses/by/4.0/ (the “License”), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.