On the growth of non-motile bacteria colonies: an agent-based model for pattern formation

TL;DR

This paper introduces an off-lattice agent-based model for non-motile bacterial colony growth, capturing diverse pattern formations and analyzing their scaling and multifractal properties, advancing understanding beyond traditional lattice-based models.

Contribution

It presents a novel off-lattice simulation approach for bacterial colonies, revealing pattern formation mechanisms and scaling laws not addressed by previous models.

Findings

Simulated a wide range of colony morphologies matching experimental patterns.

Identified a scaling relationship between interface and total cell numbers.

Analyzed multifractal properties to characterize screening effects in growth patterns.

Abstract

In the growth of bacterial colonies, a great variety of complex patterns are observed in experiments, depending on external conditions and the bacterial species. Typically, existing models employ systems of reaction-diffusion equations or consist of growth processes based on rules, and are limited to a discrete lattice. In contrast, the two-dimensional model proposed here is an off-lattice simulation, where bacteria are modelled as rigid circles and nutrients are point-like, Brownian particles. Varying the nutrient diffusion and concentration, we simulate a wide range of morphologies compatible with experimental observations, from round and compact to extremely branched patterns. A scaling relationship is found between the number of cells in the interface and the total number of cells, with two characteristic regimes. These regimes correspond to the compact and branched patterns, which are exhibited for sufficiently small and large colonies, respectively. In addition, we characterise the screening effect observed in the structures by analysing the multifractal properties of the growth probability.

Tables

Table 1. Table 1: Numerical values used for the simulation.

Variable	Value [arbitrary units]
$kT$	$1$
$\Delta t$	$0.01$
$r_b$	$1$
$m$ (nutrient mass)	$0.0001$

Table 2. Table 2: Fractal dimension D f subscript 𝐷 𝑓 D_{f} . Values reported in experimental works, our model and DLA (the error is the standard deviation). Only the results for the most branched cases are included with C = 1 / 80 , D = 512 formulae-sequence 𝐶 1 80 𝐷 512 C=1/80,D=512 and C = 1 / 40 , D = 512 formulae-sequence 𝐶 1 40 𝐷 512 C=1/40,D=512 . Greater values than for the case of DLA are expected, where only one particle diffuses at a time, unlike our case where we have many ( C > 0 𝐶 0 C>0 ).

	$D_f$
Model [$C=1/80;D=512$]	$1.760\pm 0.004$
Model [$C=1/40;D=512$]	$1.778\pm 0.003$
Experiment Matsushita1990	$1.73\pm 0.02$
Experiment Ohgiwari1992	$1.70\pm 0.02$
DLA Hayakawa1987 ; Barabasi1995	$1.71\pm 0.01$

Table 3. Table 3: Some characteristic values of the generalised dimension D q subscript 𝐷 𝑞 D_{q} (the error is the standard deviation). DLA results found in Hayakawa1987 ; Matsushita1987 are included for comparison.

	$C=1/20;D=512$		$C=1/40;D=512$		$C=1/40;D=1024$
	C.M.	L.M.	C.M.	L.M.	C.M.	L.M.
$D_{q=0}$	$-$	$1.60\pm 0.01$	$-$	$1.63\pm 0.01$	$-$	$1.53\pm 0.01$
$D_{q=1}$	$1.08\pm 0.01$	$1.08\pm 0.05$	$1.09\pm 0.01$	$1.10\pm 0.09$	$1.11\pm 0.01$	$1.09\pm 0.08$
$D_{q=2}$	$1.00\pm 0.02$	$1.00\pm 0.09$	$1.00\pm 0.02$	$1.01\pm 0.12$	$1.04\pm 0.01$	$1.01\pm 0.13$
$D_{q\gg 1}$	$0.82\pm 0.04$	$0.75\pm 0.09$	$0.82\pm 0.04$	$0.79\pm 0.15$	$0.87\pm 0.04$	$0.78\pm 0.15$
$D_{q\ll -1}$	$-$	$19.0\pm 1.5$	$-$	$15.5\pm 1.5$	$-$	$18.7\pm 2.1$

Equations

$$D_q=\frac{1}{(q-1)}\underset{l\to 0}{\lim }\frac{\ln \sum _i{P}_i^q(l)}{\ln l}.$$

$$D_1=\underset{l\to 0}{\lim }\frac{\sum _i{P}_i(l)\log P_i(l)}{\ln l}.$$

$$S\sim \{\begin{array}{ll}{N}^{1/2},& N\ll N^{\ast }\\ N,& N\gg N^{\ast }\end{array},$$

$$S/S^{\ast }\sim \{\begin{array}{ll}(N/N^{\ast }{)}^{1/2},& N/N^{\ast }\ll 1\\ N/N^{\ast },& N/N^{\ast }\gg 1\end{array}.$$

$$f(x)\sim \{\begin{array}{ll}const,& x\ll 1\\ x,& x\gg 1\end{array},$$

$$S=N^{1/2}f[(\frac{N}{N^{\ast }}{)}^{1/2}].$$

Full text

1 Introduction

The concept of active matter is relatively new within soft matter physics; the fundamental units of this type of matter, called active agents, have the particularity of absorbing energy from their environment and dissipating it in order to move, grow or replicate, among other activities DeMagistris2015 . Most of the examples of active matter are biological, such as bacteria. Although they can be seen as the simplest living organisms, they present interesting behaviours, both individually and collectively.

Bacteria exhibit many different types of movement, depending on the species and the environment, which determines the macroscopic appearance of the colony. According to Henrichsen Henrichsen1972 , six types of motility can be identified: swimming, swarming, twitching, darting, gliding (newer studies subdivide this category MarkJ.McBride2001 ) and sliding. We focus on the last one, which is a mechanism produced by the expansive forces of the colony, in combination with special properties of the cell membranes characterised by low friction with the substrate on which they grow; the bacteria do not move by their own motors, but push each other by duplicating themselves and competing for the same spaces. Despite its simplicity, its importance has been pointed out in more complex bacterial processes, such as in the formation, dispersion and restructuring of biofilms Mazza2016 . Although we will not strictly enter the biofilms field in this work, it is important to mention they are the focus of numerous studies, due to the complexity of the processes involved, their resistance to hostile environments, and the challenge they present to medicine Costerton1995 ; Mazza2016 ; Hall-Stoodley2004 ; OToole2000 ; Flemming2010 .

At the end of the 80s, Matsuyama Matsuyama1989 , Fujikawa and Matsushita Matsushita1990 , showed that the patterns of bacterial colonies obtained in the laboratory could be fractal objects. The properties of their patterns depend on two main factors: the concentration of nutrients, which influences the growth rate of the colony, and the concentration of agar, which determines the hardness of the substrate, and therefore, the mobility of the bacteria. In the absence of special forms of motility, the patterns were classified in a two-dimensional phase diagram in which five characteristic patterns were identified: diffusion-limited aggregation-like (DLA-like), Eden-like, dense branching morphology (DBM), concentric ring and homogeneous disk-like. The experiments were performed mainly with the species Bacillus subtilis Matsushita1990 ; Matsushita1998 ; Ohgiwari1992 ; Fujikawa1992 . Without self-propulsion, only DLA and Eden-like patterns are expected.

At the theoretical level, continuous models are the most traditional and extended way of studying the patterns exhibited by bacteria colonies. In them, both bacteria and nutrients (or any other variable of interest) are represented by density functions per unit area, and the spatio-temporal evolution of the system is described by systems of reaction-diffusion equations Matsushita1998 ; Kawasaki1997 ; Mimura2000 ; Giverso2015 ; Golding1998 ; Matsushita1999 ; Kozlovsky1999 ; Matsushita2004 ; Lacasta1999 ; Tronnolone2018 ; Marrocco2010 . These models are successful in describing a wide range of patterns, although they are valid only at a mesoscopic scale. To represent growth at the microscopic level, microorganisms must be represented by discrete mobile entities (agents) Ben-Jacob1994 ; Farrell2013 ; Li1995 ; Melke2010 . In this scheme, on-lattice Li1995 and off-lattice Ben-Jacob1994 ; Farrell2013 ; Melke2010 ; Santalla approaches can be chosen.

The motivation of this work is to propose a microscopic model that can explain the experimental observations, based on the fact that sliding is dominated by the mechanical interaction between the bacterial cells. As has been said, continuous models work well only on a mesoscopic scale, whereas in agent-based models, if space is discretised as Euclidean networks, mechanical laws cannot be used. In spite of being computationally expensive, in this work, we choose an off-lattice model, in order to represent our agents as rigid bodies governed by laws of mechanics. Thus, we can analyse the growth of bacterial colonies on a microscopic level. The off-lattice approach also avoids anisotropies in the patterns exhibited by the colonies induced by the discretization of the space.

A typical way to characterise the complex structures that arise in surface growth is by means of the Hausdorff dimension, often referred to as the fractal dimension. However, the fractal dimension is not a unique descriptor, as it was shown that two structures may have the same fractal dimension but are fundamentally different Hayakawa1987 . In order to describe structures more deeply and unequivocally, the determination of the multifractal properties of an associated measure (e.g. growth probability) offers a suitable supplement to the sole measurement of their fractal dimension. This is an entropy-based approach Kinsner2005 , classified this way to differentiate it from other analyses that rely only on metric concepts. Here the scaling properties are analysed for variations in different parts of the pattern, which are overlooked by a simple measurement of the fractal dimension. As it is arduous to treat growth models with nonlocal rules analytically Barabasi1995 , it may be interesting to characterise an associated measure such as the growth probability to gain some insight about the process.

A way to describe the multifractal behaviour is through the generalised dimensions $D_{q}$ (also known as Rényi dimensions). If one covers the support of the measure (set of all points where the measure is positive) with a set of boxes of size $l$ and defines a probability $P_{i}(l)$ (integrated measure) in the $i$th box, the generalised dimensions $D_{q}$ correspond to the scaling exponents for the $q$th moments of $P_{i}$, defined by $\sum\limits_{i}P_{i}^{q}(l)\sim l^{(q-1)D_{q}}$. In this context, $q$ is typically referred to as the order $q$ of the generalised dimension $D_{q}$. Solving for $D_{q}$ and taking the limit of $l\rightarrow 0,$ the conventional expression for the generalised dimensions is given by

[TABLE]

For the case $q=1$, the L’Hôpital’s rule must be used; thus,

[TABLE]

The generalised dimensions are exponents that characterise the non-uniformity of the measure; the positive orders $q$ accentuate the regions with higher probabilities while the negative $q$’s the opposite.

In early works, $D_{q}$ was only defined for $q\geq 0$ Hentschel1983 , so the previous definition had the particularity that the first value of the set, i.e. $D_{q=0}$, corresponded to the Hausdorff dimension of the support (because all the boxes have the same weight). There are also specific names for other certain values. For example, $D_{q=1}$ is known as the information dimension, which is interesting in the case of diffusion-limited aggregations, since it can be physically interpreted as the fractal dimension of the active region, i.e., the unscreened region Hayakawa1987 ; $D_{q=2}$ is known as the correlation dimension; D_{q\rightarrow\text{\textpm}\text{\infty}} are known as the Chebyshev dimensions, which are calculated with the maximum and minimum probabilities, respectively; equivalences with other dimensions definitions can be made, even for fractional $q$ values Kinsner2005 .

This is not the only multifractal analysis possible. In other works, the singularity spectrum is computed Halsey1986 ; Chhabra1989a , which is closely related to the Rényi dimensions by a Legendre transform. Temporal fractals can also be studied, where the local scaling properties are now related to time behaviour Ivanov2009 . We will use the generalised dimensions to characterise quantitatively the patterns produced by our model.

2 Model

In this paper, we model the growth of non-motile bacterial colonies under different environmental conditions, specifically, nutrient concentration and nutrient diffusion. The growth rules are inspired by biology, as we capture the essential characteristics of bacteria without losing simplicity. We consider a two-dimensional and off-lattice space, which allows us to consider mechanical interactions between the agents, as we will explain below.

There are two kinds of particles in the model, nutrient particles and bacterial cells. Both of them have physical properties such as size, mass, position, velocity and might have applied forces. Nutrient particles are idealised as Brownian particles, so its initial velocities follow the Maxwell-Boltzmann distribution and evolve according to a Langevin equation of the form $m\dot{v}(t)=-\frac{\kappa T}{D}v(t)+f(t),$ where $\kappa$ is the Boltzmann’s constant, $T$ is the temperature and $D$ is the diffusion coefficient. The function $f(t)$ is a stochastic force whose components follow a Gaussian probability distribution with mean zero and standard deviation $\sigma=\kappa T\sqrt{2/D}$. The Langevin equation is numerically integrated using a small time step $\varDelta t$, following the explanations in The Fokker-Planck Equation by H. Risken Risken1989 . Nutrient particles are considered point-like, non-interacting with each other and with a small mass. The bacterial cells are modelled as rigid circles with radius $r_{b}$, so they can interact with each other through normal forces. The numerical values used in the simulations can be found in Table 1.

In addition to the physical properties, bacterial cells have two biological characteristics: they can be fed and reproduce. The first one is an interaction with the nutrients, which are absorbed by the cell when they are in contact. Reproduction is the process by which the bacterium duplicates: an identical copy of the cell is generated in the same position as the original, so they overlap. Then, they are disaggregated by opposing velocities in a random direction. As a consequence, these cells may collide elastically with neighbours, according to the mechanics of rigid bodies, as shown in Fig. 1. When the cells stop overlapping, they stop moving and become static because the medium is considered to be very viscous so that the momentum gained by the collisions is immediately dissipated. All the calculations to resolve collisions and overlaps are based on an iterative constraint solver introduced by Erin Catto Catto2005 .

In summary, reproduction causes both the movement of newborn cells and their neighbours, representing sliding motility.

The simulation begins with a single cell in the origin of coordinates in a 2D substrate and a given quantity of nutrient particles diffusing in space, according to the concentration and diffusion coefficient specified. When a nutrient particle touches the cell, it is absorbed and the bacteria duplicates. Now the colony is formed by two cells, which can absorb nutrients and reproduce. The process continues in this way and the colony grows progressively. All the bacteria have a time delay (20 integration time steps), during which they can’t duplicate; this rule ensures that no duplication occurs while newborn cells are still overlapping and it is consistent with biological observations, e.g., Bacillus subtilis species has a delay of $\sim 25$ minutes between duplications Golding1998 .

In order to keep the nutrient concentration constant, there is a ring that acts as a nutrient reservoir located at a given distance from the most external position of the bacteria. This distance increases progressively as the colony grows, so the separation is $60r_{b}$ at least. The nutrient concentration within the ring ($r>60r_{b}$) is kept constant at a specified value. But, closer to the colony ($r<60r_{b}$), the concentration drops due to nutrient absorption by the bacteria. Periodic boundaries conditions are considered for the outer side of the ring (Fig. 2).

The growth stops when the colony reaches a radius of $600r_{b}$, when characteristic patterns are fully developed. This implies that we have up to half a million cells forming the colony, depending on the parameters.

3 Results and discussions

3.1 First characterization of the structures

In order to see the variety of morphologies that the model can produce, we choose several different values of nutrient concentration and nutrient diffusion, and register the position of the bacteria along the perimeter of the colony over time and average over a hundred realizations for each set of parameters. A morphology diagram is shown in Fig. 3a, where it can be seen that it is possible to generate round and compact colonies, as well as ramified, going through a variety of intermediate patterns. Similar morphological crossover can be seen in Figs. 3b–e, which corresponds to the experiments carried out in Matsushita1998 ; Ohgiwari1992 . The fractal dimension of intermediate patterns in experiments was reported for the case of the most ramified one. In Table 2, we summarise the results found in the bibliography Matsushita1990 ; Ohgiwari1992 and ours (for the most ramified cases), which are in good agreement.

3.2 Scaling properties

Despite using different values for the parameters, it is observed that the curves of the number of bacteria at the interface S versus the total number N show two power-law regimes. The first regime corresponds to initial compact structures $S\sim N^{1/2}$, while the second regime corresponds to ramified structures with $S\sim N$, as shown in Fig. 4a. These two behaviours are characteristic of the Eden and DLA models, respectively.

To characterise the crossover between regimes, we compute the total number of bacteria $N^{*}$ at which the crossover takes place. To achieve this, we simply fit power-law functions in the tails of the $S$ vs $N$ curves and compute the intersection. More sofisticated methods for an automatic determination of crossovers can be found in Bashan2008 . After dividing $N$ by $N^{*}$ in each of the data sets, the $Y$-axis is divided by some value $S^{*}$ looking for a satisfactory collapse of the curves. We found that the best collapse occurs when $N^{*}=S^{*2}$, as shown in Fig. 4b.

It can be seen that $N^{*}$ depends on the diffusion $D$ and the concentration $C$, having an increasing relationship with both (Fig. 5).

Taking these observations into account, an attempt is made to establish a scaling law. We know that the behaviour of S is:

[TABLE]

where $N^{*}=N^{*}(C,D)$. The curves collapse dividing $N$ and $S$ by $N^{*}$ and $S^{*}$, respectively, so:

[TABLE]

Then, having validated the relationship $N^{*}=S^{*2}$ and proposing the scaling function

[TABLE]

the relation between $N$ and $S$ can be written as:

[TABLE]

This result suggests that if we allow $N$ to grow sufficiently, branches will always be generated, after a critical number of cells is reached, dependent on the parameters $C$ and $D$.

3.3 Multifractality of the growth probability

The growth probability of each region of the colony can give information about why a certain pattern displays. Every cell duplicates when a nutrient particle is captured, so the growth probability is associated with the probability that a diffusing particle reaches the site where the cell is. We use two methods to estimate this probability, focusing on the final stage of the colony. The first one consists on counting how many nutrient particles are absorbed by each cell without letting it duplicate, i.e., the colony is "frozen" and the growth probability of each cell is computed dividing this counting by the total of particles incorporated by the whole colony (we use approximately $10^{5}$ particles). We will refer to this method as C.M. The disadvantage with this method is that it does not estimate low probabilities well, because several million particles may be captured by the colony in total, but the internal regions may hardly incorporate any. Due to this, we also solve the Laplace equation $\nabla^{2}\phi=0$, where $\phi$ represents the nutrient concentration, by the relaxation method Li2003 , where $\phi=1$ at infinity and $\phi=0$ along the perimeter of the colony, as it can be seen that the growth probability is proportional to the gradient of the potential $\nabla\phi$ Hayakawa1987 . We use an iteration error of $10^{-5}$, after checking that the multifractal curves do not vary appreciably. In order to use this method, referred to as L.M. henceforth, properly, space has to be discretised, so some differences with the C.M. are expected.

In Fig. 6, it is shown how uneven is the number of nutrient particles consumed between the outer and inner regions of a ramified colony. This phenomenon is usually referred to as ‘shadowing’ or ‘screening’ effect and is more or less noticeable depending on the parameters. As the structures that emerge from the simulations are fractal, the proper way to study this effect is by the multifractal formalism, explained in the first section.

In Fig. 7, the generalised dimension $D_{q>1}$ curve is plotted for different morphologies. It can be seen that the probability associated with a ramified colony presents a strong multifractality since $D_{q}$ varies significantly with $q$. It is also included in Fig. 7 the curve for a diffusion-limited aggregate Matsushita1987 for comparison. The standard deviations are presented in Table 3. Unfortunately, this analysis cannot be carried out in very compact colonies since the fractal regime is very short to be reliable or it is not observable. Nor can it be done in the early stages of the growth process for the same reason.

In Fig. 8, the generalised dimension is plotted again, but now including the $q<1$ interval. Only the results obtained by the L.M. can be used in this interval. The value of $D_{q\ll-1}$ is very interesting, because it quantifies the shadowing effect, making evident the differences between different morphologies. It is worth noting that $D_{q=0}$ should not be equal to the fractal dimension $D_{f}$ presented in the previous subsection because $D_{q=0}$ is the fractal dimension of the support of the measure, i.e., the perimeter of the colony, while $D_{f}$ is the fractal dimension of the area. All of these characteristic points are summarised in Table 3, which also includes the values corresponding to diffusion-limited aggregation as a comparison Hayakawa1987 ; Matsushita1987 . The asymptotic values are estimated with $q=25$ and $q=-25$.

Note that the curves of the generalised dimension are always above the case of DLA. In the region of $q>0$, taking into account the calculated standard deviations, the differences are not as noticeable, but they are in the region of $q<0$. In this region, where the measurement best distinguishes each case, they depart notoriously from the case of DLA. Always considering branched cases, it is observed that higher values are associated with higher $C$ and $D$ values, which can be understood if we associate this measure with the screening phenomenon. The larger are $C$ and $D$, the thicker and narrower the branches and the fjords become, respectively, so the screening increases to the interior areas. Nevertheless, note that there is a limit on how large the parameters can be. If $C$ and $D$ are too large, so the branches disappear, the screening effect is barely noticeable. The growth probability becomes almost uniform, and the multifractality should be lost.

4 Conclusions

The goal of this work is the construction of a model based on basic theories of physics, capable of generating a variety of complex patterns observed in bacteria colonies. Under the hypothesis that there is a single collective movement mechanism behind the different morphologies (sliding), the different results are achieved by varying parameters of the environment outside the colony, without changing the behaviour of the agents. Under these precepts, we manage to generate patterns from the most round and compact to extremely ramified, going through different intermediate morphologies. The different approaches used to characterise the structures also allow comparison with the two most studied models that predict patterns of bacterial colonies, the Eden model and DLA.

The fractal dimension analysis of simulations with branched colonies shows a fractal dimension compatible with experiments Matsushita1990 . Although there are distinct differences between the present and the DLA model, their fractal dimensions are in good agreement, which suggests that the most ramified cases considered are close to the diffusion limit.

On the other hand, the characterization using scaling laws show that there are two characteristic growth regimes, one compact and one branched. According to the relation found, the crossover between regimes occurs at a critical number of cells, that depends on the parameters for nutrition concentration $C$ and nutrition diffusion $D$. Thus, branches will always be generated for finite values of these parameters.

The multifractal measurement shows strong multifractality for the ramified cases. In case the order $q$ of the generalised dimension $D_{q}$ is $q>0$, the curves of $D_{q}$ of the simulations of different values of $C$ and $D$ are close to each other, but in the region of $q<0$, their differences become notoriously. All these curves are, however, always above the DLA curve and, as expected, approach to the DLA curve for lower values for $C$ and $D$. Higher values are associated with higher $C$ and $D$ values because fjords into the interior are narrow.

Unfortunately, a more in-depth comparison between simulation and experimental data cannot be carried out due to the lack of quantitative experimental data. To date, there are almost exclusively qualitative characterizations of the morphologies of bacterial colonies, which only in some rare cases provide additional information on the fractal dimension (which is not a complete indicator). Although the methods that we use in this work to calculate the generalised dimensions cannot be used in experiments, there are other methods that might be used, such as the one described by Ohta and Honjo Ohta1988 , based on associating probabilities according to the variation of the area occupied by the colony in a certain section. Only such a deeper experimental analysis would offer a complete characterization of the processes involved in the structure formation of bacterial colonies and would allow contrasting the proposed model with experiments.

Conflict of interest

Authors declare that there are no conflicts of interest.

Acknowledgements

L.A.B. and L.V. thanks UNMdP and CONICET (PIP 00443/2014) for financial support. D.H. thanks UNMdP and CONICET (PIP 0100629/2013) for financial support.

Author contribution statement

L.V. and D.H. designed the model. L.V. implemented the computational framework. All authors contributed to the analysis of the data and in the writing of the manuscript.