Crystallization and melting of bacteria colonies and Brownian Bugs

Spontaneous emergence of clustering is a widespread phenomenon in population biology, ecology, material science, and other fields Murray ; Levin . Either passive or active “particles” (trees, bacteria, plankton, etc.) bunch together forming dense localized clusters embedded in an, otherwise, almost empty surrounding space. The patchy distribution of plankton in the ocean surface plankton , the spatial distribution of vegetation in semi-arid regions vegetation , or the fascinating patterns generated by bacteria grown in Petri dishes Budrene ; BJ are some examples. Some simple mechanisms leading to clustering have been described:

(i) Non-interacting inertial particles moving in a fluctuating environment (as a turbulent flow) may become clustered Coalescence . This path coalescence mechanism has been claimed to contribute to plankton patchiness.

(ii) Branching and annihilating Brownian particles (also called super-Brownian processes) tend to bunch together. This type of clustering has its roots in the fact that offsprings are created within a local neighborhood of their parents and die anywhere, giving rise to an overall tendency to form localized colonies Young .

(iii) In some cases, clustering stems from the interaction with a second “field” influencing the dynamics (water in models for vegetation patchiness vegetation , nutrients or chemicals in models for chemotactic bacteria clustering Levin , etc.) through some type of feedback mechanism. Reaction-diffusion equations have proven to be a convenient tool to model clustering at a mesoscopic level in this case Murray ; Cross ; Walgraef ; Fuentes .

Contrarily to the cases above, there are clustering models in which individuals interact with each other in a direct way. For example, their effective interaction can be such that the reproduction rate (or some other individual dynamical property) is diminished by a factor depending on the relative abundance of neighboring individuals. As a well documented example of this, let us mention the empirically observed Janzen-Connell effect in ecology: owing to various circumstances (presence of specialized pests or predators, competence for resources, etc.), the effective reproduction rate of an individual decreases with the number of conspecific specimens surrounding it JC . Similar effects appear also in bacteria colonies, social phenomena, and in systems exhibiting collective motion collective .

In many of these systems, clusters are distributed in space in a disordered way but, very remarkably, in some other cases they self-organize forming rather regular patterns. As a simple example to bear in mind, let us focus on the strikingly ordered patterns self-generated by Escherichia Coli, Salmonella typhimorium, and other bacteria grown in Petri dishes. When grown in a substrate of nutrients under adequate conditions these bacteria colonies structure themselves into clusters which, in their turn, self-organize in spirals, squares, or crystal-like hexagonal arrays (as a beautiful illustration, see Fig.1 in Budrene or page 529 of BJ ).

How ordered can such two-dimensional patterns be? This question resembles very much the problem of crystal ordering in equilibrium solids and the existence of a melting (freezing) solid/liquid (liquid/solid) transition. For two-dimensional systems in thermodynamic equilibrium the Mermin-Wagner-Hohenberg (MWH) theorem MW rules out continuous symmetries to be spontaneously broken in the presence of fluctuations. This covers the case of translational invariance and, therefore, two dimensional crystals/solids cannot exhibit true long-range translational order (destroyed by low-energy Goldstone modes). Indeed, the most popular equilibrium melting theory (to be described in some detail below KT ; KTHNY ; Strandburg ; Nelson ) predicts the melting transition to occur in two-stages, and to include an intermediate hexatic phase in between a quasi long-range ordered solid phase and the disordered (or isotropic or liquid) phase. Alternatively, a competing theory Chui predicts a unique discontinuous transition between such two phases.

Note, notwithstanding, that all the clustering problems described above occur away from thermodynamic equilibrium. Actually, many of them exhibit absorbing states, a quintessential example of irreversible nonequilibrium dynamics (fluctuations can lead from activity to the absorbing state, but not the other way around AS ).

Does the MWH theorem applies to nonequilibrium problems as self-organized bacteria colonies? Do they exhibit a melting/freezing transition similar to that occurring in equilibrium? Do they have an hexatic phase?

Nonequilibrium problems violating the MWH theorem are well known. A notorious example is provided by the model of self-propelled particles proposed by Vicsek et al. collective which illustrates how communities (flocks, schools) of interacting individuals (birds, fishes) move coherently, with true long-range order, adopting a collective non-vanishing velocity (which breaks the continuous rotational symmetry). Soon after the introduction of such a model it was proved in an elegant way that, indeed, true long-range order emerges owing to inherently nonequilibrium effects TT (see also Racz for a related, though different, nonequilibrium mechanism leading to true long-range two-dimensional ordering). Is the MWH violated in a similar way for the clustering problems described above?

Aimed at answering the previously posed questions and to scrutinize the role of fluctuations in nonequilibrium clustered systems exhibiting self-organized patterns, we analyze here a simple model of individual interacting particles bugs1 ; bugs2 ; framos exhibiting clustering as well as ordered patterns.

Our main conclusions are as follows: while we find a solid phase with quasi long-range translational order and a disordered phase in analogy with equilibrium situations, we do not find any evidence either of the hexatic phase predicted by the most popular theory for equilibrium melting KT ; KTHNY ; Strandburg ; Nelson , nor we find the discontinuous transition predicted by competing theories Chui . Instead, we report on a one-stage continuous transition analogous to what was reported in numerical investigations of other equilibrium JJ and nonequilibrium QG melting problems. An explanation justifying our findings, supported by the analysis of topological defects and likely to be valid for other systems, is proposed.

The paper is structured as follows. First we define the model, describe its basic ingredients and introduce a Langevin equation capturing all its relevant traits at a continuous level. We discuss briefly its main phenomenology: emergence of clusters and ordering into hexagonal patterns. Afterwards, we report on its numerical analysis paying attention to the melting transition (including a careful analysis of topological defects) and compare our results with standard equilibrium melting theories. Finally, we discuss our findings from a general perspective and present the conclusions.

The individual-based model we study here, the “interacting Brownian bug” (IBB) model, was recently introduced by two of us bugs1 ; bugs2 . It consists of branching-annihilating Brownian particles (bugs, bacteria) which interact with each other within a finite distance, $R$ bugs1 . Particles move off-lattice in a $d$-dimensional $[0,L]^{d}$ interval with periodic boundary conditions, obeying a dynamics by which particles can:

• •

  diffuse (at rate $1$) performing Gaussian random jumps of variance $2D$,

• •

  dissappear spontaneously (at rate $\beta_{0}$),

• •

  branch, creating an offspring at their same spatial coordinates with a density-dependent rate $\lambda$:

  $$\lambda(j)=max\{0,\lambda_{0}-{N_{R}(j)}/{N_{s}}\},$$

  (1)

  where $j$ is the particle label, $\lambda_{0}$ (reproduction rate in isolation) and $N_{s}$ (saturation number) are fixed parameters, and $N_{R}(j)$ stands for the number of particles within a radius $R$ from $j$.

The control parameter is $\mu=\lambda_{0}-\beta_{0}$ while the function $max()$ enforces the transition rates positivity. See otros where similar models in which the death or the diffusion rate are density-dependent are studied.

The main phenomenology of the IBB is as follows bugs1 ; bugs2 ; framos . For large values of $\mu$ there is a stationary finite density of bugs (active phase) while for small values the system falls ineluctably into the vacuum (absorbing phase). Separating these two regimes there is a critical point at some value $\mu_{c}$, belonging to the directed percolation (DP) universality class framos characterizing in a robust way transitions into absorbing states AS .

In the active phase, owing to the local-density dependent dynamical rules particles group together forming clusters (see Fig. 1 (a) and (b)) provided that the diffusion constant is small enough (for large values of D, homogeneous distributions are obtained). Such clusters have a well-defined typical diameter (which can be analytically estimated nb ) and a characteristic number of particles within, which depend on the parameters $R$, $N_{s}$, and $\mu$ bugs1 ; bugs2 . Well inside the active phase, when the clusters start filling the available space they self-organize in spatial structures with remarkable hexagonal order (see Fig.1 (b)).

From the theoretical side, an important breakthrough is that the IBB model can be cast into a continuous stochastic equation bugs1 . Indeed, by applying standard Fock-space (Doi-Peliti GF ) techniques, a Langevin equation including the main relevant traits of the problem in a parsimonious way can be derived (see Appendix A in bugs1 and references therein). The Langevin equation for the local density density of bugs $\rho({\bf x},t)$ reads (in the Ito representation)

where the noise amplitude $\sigma$ is a function of the microscopic parameters and $\eta(\textbf{x},t)$ is a normalized Gaussian white noise. This includes only the leading terms in a density expansion; for instance, a higher order noise term appear in the mapping bugs1 but it does not alter the results reported in what follows in any significant way.

Note that, leaving aside the non-local saturation term, this equation coincides with the Reggeon-field theory or Gribov process, describing at a coarse-grained level systems with absorbing states in the directed percolation (DP) class (see AS and references therein). Let us underline the presence of a square-root multiplicative noise. Note also that the deterministic part of Eq.(II), including a non-local saturation term, is identical to the equation proposed by Fuentes et al. Fuentes . Eq.(II) is therefore a simple (the simplest) stochastic generalization of such model.

The main advantage of the continuous Langevin equation above is that it allows for analytical studies and permits to scrutinize the effect of fluctuations (by just tuning the noise amplitude). It also constitutes a more general, elegant, and compact formulation of the problem. For these reasons, we choose Eq.(II) as the starting point of our forthcoming analyses.

Analytical studies of the deterministic part of Eq.(II) (i.e. mean field analyses) have already been performed in bugs1 ; bugs2 ; Fuentes . They permit, for instance, to understand the wave-length instability mechanism leading to pattern formation, as well as many other aspects. In order to analyze the full stochastic Eq.(II) and, in particular, its associated melting/freezing transition, one needs to resort to computational studies.

However, integrating numerically Langevin equations with square-root noise is a highly non-trivial task; owing to the fact that for small density values the square-root term (with multiplies a random number) becomes larger in amplitude than the deterministic terms, standard integration schemes (Euler, Runge-Kutta, etc.) lead ineluctably to unphysical negative values of the density field dickman , and this pathology is not easily cured in any naïve way. Luckily, a very efficient scheme, specifically devised to overcome such a difficulty, has been recently proposed Dornic .

The method is a split-step algorithm in which the system is discretized in space and two evolution operations are performed at each discrete time-step: (i) first, the noise term is treated in an exact way, by sampling the conditional probability distribution coming out of the (exactly solvable) associated Fokker-Planck equation at each site. By sampling in an exact way such a distribution an output is produced at each site. (ii) Afterwards, the remaining deterministic terms are integrated using any standard scheme taking as the input at each site the output of the previous step at each site. More details and applications can be found in Dornic .

To implement the split-step scheme to integrate Eq.(II) in two dimensions, we discretize the space, by introducing a lattice of linear size $L=256$ or $~512$ ($L\approx 1024$ is already at the limit of our present computational power). We fix the discrete time step to $\Delta t=0.25$, $R=10$, $N_{s}=50$, $D=0.25$ (which is small enough to have clustering) and use either $\mu$ or $\sigma$ as a control parameter. For all of the simulations reported here we initialize the system with a homogeneous initial density, $\rho({\bf{x}},t=0)=\rho_{0}$, leave it relax towards its stationary state (reached typically after $10^{5}$ Monte Carlo steps), and perform steady-state measurements (averaging over, at least, $10^{5}$ configurations).

First of all, we verify that Eq.(II) reproduces qualitatively all the basic phenomenology of the microscopic IBB model: (i) an absorbing phase for small values of $\mu$, (ii) an active disordered phase, encountered by increasing $\mu$ for a fixed $\sigma$ and (iii) an active ordered crystal-like phase which is reached by further increasing the value of $\mu$ or alternatively, fixing $\mu$ and reducing the noise amplitude, $\sigma$ (see Fig.(1 (c) and (d))).

Separating (i) and (ii) there is a directed-percolation-like phase transition, while our focus here is on the transition from the disordered active phase (ii) to the ordered self-organized one (iii). In order to study the effect of the noise on such a transition, from now on, we fix $\mu=2.0$, well into the active phase, and use the noise amplitude, $\sigma$, as a tuning parameter. Fig. 1 (plots (c) and (d)) shows two snapshots obtained for $\mu=2.0$ with $\sigma=2.5$ and $\sigma=1.5$ respectively. While in both cases clusters of localized activity ($\rho(\textbf{x})\neq 0$) exist, only in the second one clusters are self-organized into an ordered hexagonal array.

The evidence accumulated so far, both from cluster analyses and structure function measurements, reveals the existence of a melting/freezing transition, somewhere around $\sigma=2.39$. The key question we face now is: does such a nonequilibrium transition exhibit the well-known universal features of standard equilibrium melting/freezing in two-dimensional systems?

Note, in addition, that here the number of clusters (particles) is not constant but fluctuates, and its average value varies as a function of the control parameter, $\sigma$. Is this relevant for the properties at the melting/freezing transition?

Before tackling these problems, for the sake of completeness, and for future reference, in this section we briefly summarize the main results of the celebrated standard theory of two-dimensional melting: the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory KT ; KTHNY ; Strandburg ; Nelson . We also discuss briefly an alternative competing theory.

The KTHNY is based on a statistical physics analysis of topological defects, i.e. particles with a number of nearest neighbors (assuming a Voronoi or Wigner-Seitz construction) other than $6$. Dislocations perturb translational order and disclinations hinder orientational order Nelson ; Strandburg . A detailed inspection of hexagonal ordering in the presence of fluctuations reveals that disclinations correspond to free monopoles, either five-fold or seven-fold, where five and seven refer to the number of nearest neighbors as measured in a Voronoi or Wigner-Seitz construction. Analogously, dislocations can be identified with tight pairs (i.e. dipoles) of a five-fold and a seven-fold disclinations. The disordered (or isotropic or liquid) phase is characterized by the proliferation of defects: both monopoles and dipoles.

The main prediction of the KTHNY theory is that, contrarily to what happens in higher dimensional systems, where the melting occurs discontinuously at a unique transition point, in two-dimensional systems melting occurs in two stages. Translational and orientational order lose their stability at different Kosterlitz-Thouless-like KT critical points where dislocations and disclinations, respectively, unbind.

The theory assumes that in the solid phase there are nor free dipoles nor monopoles, but only quadrupoles (low-energy excitations), that the number of dislocations/dipoles throughout the first (melting) transition is small and that they are generated progressively in a smooth way as the temperature is risen. This allows to treat the system at the melting transition as a weakly interacting gas of dipoles/dislocations. An analogous assumption is made for monopoles/disclinations at the second transition point where monopoles unbind from dipoles.

The three phases put forward by the KTHNY theory are as follows KTHNY ; Strandburg ; Nelson (see Fig. 4 for a graphical illustration):

• •

  Only at zero temperature the lattice ordering can be perfect while, for any non-vanishing temperature, defects appear. Below a first critical point (denoted $\sigma_{m}$ here) defects are tight together in quadrupoles, hindering translational order. As a result, translational correlation functions decay algebraically in space (with continuously varying exponents), as corresponds to quasi long-range order. Orientational correlations (see below for a precise definition) decay at long distances to a non-vanishing constant value, corresponding to true long-range orientational order explanation .

• •

  Above a first critical point, $\sigma_{m}$, dipoles/dislocations unbound from quadrupoles, destroy translational order (i.e. translational correlations decay exponentially fast). They also affect orientational correlations which decay algebraically with a continuously variant exponent $\eta_{6}$ obeying $\frac{1}{4}\leq\eta_{6}\leq\frac{1}{3}$ and a diverging correlation length, i.e. they exhibit quasi-long-range orientational order. This is the, so called, hexatic phase.

• •

  Above a second (isotropic) critical point, $\sigma_{i}$, monopoles/disclinations unbind from dipoles, hindering quasi long range orientational order. Both translational and orientational correlations decay exponentially. The associated correlation lengths diverge as stretched exponentials of the form $a\exp(b\Delta^{-1/2})$ where $\Delta$ is the distance to $\sigma_{i}$. This is the isotropic (also called disordered or liquid) phase.

This scenario has been verified in a number of numerical Hexatic ; Strandburg and experimental Hexatic-exp ; Strandburg studies, while it was not verified in others JJ ; QG .

Despite of its success, the KTHNY is not the only plausible theory of two-dimensional melting. A competing one was proposed by Chui Chui , who argued that some systems should exhibit a unique first-order melting transition mediated by the appearance of “grain boundaries”. In this picture, chains of defects appear limiting ordered grains and separate neighboring mismatching domains. The main difference between this theory and the KTHNY one is that here the transition appears owing to a collective excitation of defects. In some systems, it has been shown that the melting transition can change from KTHNY to first order upon changing some parameter Strandburg , so both theories can be taken as complementary.

To determine the plausibility of the KTHNY scenario (or, instead, that of competing theories) for the nonequilibrium transition under study, we analyze in this section translational and orientational order in numerical simulations of Eq.(II).

Inspired by the analyses in QG , in this section we check the KTHNY assumptions rather than its predictions. As briefly explained above, the theory assumes that the number of dipoles/dislocations (monopoles/disclinations) throughout the first (second) transition is small and that they are generated progressively in a smooth way as the temperature is risen. This allows to treat the first (second) transition as a weakly interacting gas of dipoles/dislocations (monopoles/disclinations).

To check if this picture describes properly the behavior of our model, we scrutinize how defects appear and proliferate upon rising the noise amplitude. In Fig. 11 we plot defect maps for $4$ different values of $\sigma$. The lines correspond to the Voronoi construction for a given configuration; five-fold, six-fold, and seven-fold clusters are marked with black circles, central dots, and red stars respectively (higher and lower order defects correspond to blank clusters). As the noise amplitude is increased the total number of defects grows. While below the critical point only quadrupoles are significatively present, nearby the transition point isolated dipoles start unbinding. Even if unbound, they have a clear tendency to bunch together, and indeed, at slightly larger noise-amplitudes, in the isotropic phase, defects bunch together showing a tendency to form string-like structures. Deep into the isotropic phase defects proliferate, and extend through the system keeping, in any case, the propensity to form condensed string-like structures.

In this respect, it is noteworthy that Fisher, et al. Fisher predicted, within the KTHNY theory, that monopoles can appear at correlated locations showing a strong tendency to arrange themselves in small-angle grain boundaries. Contrarily, the condensates we find seem to be dominated by strings of dipoles.

To quantify the phenomenology above, we measure the density of defects (quadrupoles, dipoles and monopoles) as a function of noise amplitude (see Fig.12). The density of defects (of any type) grows from a small value below $\sigma=2.3$ to a large value in the isotropic phase. The increase is very steep around $\sigma=2.38$ for both dipole and monopole densities, suggesting the existence of a unique unbinding transition.

To complement the previous plot, the inset of Fig. 12 shows the density of isolated dipoles/dislocations and monopoles/disclinations versus $\sigma$, confirming the previous observations, that the steepest increase occurs at $\sigma\approx 2.39$. Note that the total number of monopoles grows continuously upon increasing $\sigma$, while the presence of a relative maximum in the density of dipoles reflects the competition between two conflicting tendencies: the unbinding of dipoles from quadrupoles and the liberation of monopoles (and formation of condensed defects) from free dipoles.

Finally, Fig. 13 shows the densities of isolated dipoles and of condensates as a function of the density of monopoles. This plot reveals a strong correlation between monopoles and dipoles and, more importantly, it illustrates how they both vanish at roughly the same point, suggesting again (up to finite size effects) the existence of a unique one-stage transition. Indeed, if the dipoles unbinding occurred before the monopoles one, the line joining the circles should intersect the vertical axis which is not the case. It is also at such a unique transition point that condensed defects appear, being their density (roughly) linearly correlated with the density of monopoles.

In summary, a detailed analysis of defects in Voronoi tessellations supports the interpretation that both dipoles and monopoles unbind at a unique critical point. At such a transition point condensed defects are generated in an apparently continuous way. In a future publication we will analyze in a more detail the correlations between isolated monopoles, dipoles and condensates to further clarify the analogies and differences with the standard KTHNY scenario, and try to develop a theoretical framework for one-stage continuous melting transitions.

Motivated by the observation of remarkably regular arrays of clusters formed by bacteria growing in Petri dishes and related problems, we have revisited the individual-based interacting Brownian bug (IBB) model in which birth rates are local-density dependent bugs1 ; bugs2 . Apart from an absorbing phase transition, this model shows a transition from a disordered active phase, in which particles aggregate in localized clusters, to an ordered or crystallized active phase, in which clusters self-organize forming hexagonal arrays. For the sake of generality and aimed at facilitating analytical studies, rather than studying the discrete model itself we have analyzed its equivalent continuous Langevin representation. This stochastic equation, which involves a non-local saturation term as well as a square-root multiplicative noise, can be numerically integrated by discretizing it and employing a recently introduced efficient integration scheme specifically designed to deal with square-root multiplicative noise Dornic .

The main results we obtain are:

(i) First, we have shown explicitly that the continuous model (truncated to include only the leading terms) reproduces the phenomenology of the original discrete model, including an absorbing phase in which the stationary particle density is zero, a disordered active phase in which the density field is localized in clusters of activity surrounded by empty regions and, finally an ordered phase with hexagonal patterns.

(ii) We have studied the ordering transition by analyzing cluster properties (average velocity, mean-displacement, etc.) as well as by means of analyses of the structure-function. These studies reveal that a melting/freezing transition indeed occurs at some value of the noise amplitude: for small noises, clusters are trapped into hexagonal configurations as a result of a collective effect, while for noise amplitudes above threshold they have much more mobility.

(iii) To better understand and quantify the transition we have studied both translational and orientational correlation functions. In particular, we have verified that, in the solid phase, translational correlation functions exhibit generic power-law decay even if with system-size induced cut-offs. The value of the decay exponent at the melting transition being in excellent agreement with the KTHNY prediction. However, we do not found any evidence of an hexatic phase, contrarily to what predicted by the standard (Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY)) theory. This same conclusion is also borne out by analysis of the global orientational order-parameter. Up to the numerical resolution limit, this suggests the existence of a unique continuous transition point.

(iv) Above such a transition point, i.e. in the disordered or isotropic phase, the correlation functions have been shown to decay as stretched exponentials, in excellent agreement with the KTHNY predictions.

(v) We have performed a detailed analysis of defects and found that both free dipoles/dislocations and monopoles/disclinations seem to unbind from quadrupoles at a unique transition point, above which also defect condensates are formed. This is in contrast with the KTHNY scenario and leads to a continuous one-stage melting transition.

(vi) The unique transition is continuous, so it cannot be explained by the main alternative theory to KTHNY Chui which predicts a first-order melting transition.

Summing up, the phenomenology of this nonequilibrium model can be only partially described by the KTHNY theory. While the melting/freezing transition is indeed characterized by the smooth unbinding of defects this occurs through a unique continuous transition. The reason for this seems to yield in the fact that once dipoles/dislocations unbind, the perturbation they generate around them is large enough as to unbind also monopoles/disclinations: dipoles, monopoles and the string-like condensates they form are strongly correlated.

Note that, the nonequilibrium microscopic dynamics of the interacting Brownian bug model (or its equivalent Langevin representation) is responsible for the generation of mesoscopic clusters. Once such clusters are generated, they interact in an effective way, and the physics at a macroscopic scale does not seem to differ in any essential way from other equilibrium problems JJ ; QG for which a similar one-stage continuous melting transition has been reported. We believe that this scenario is not specific of nonequilibrium systems but is determined by the way in which defects interact among themselves. This might not depend on the equilibrium versus nonequilibrium nature of the process but rather on other structural details influencing the way defects interact among themselves. A more detailed analysis of defects and defect-correlation will be investigated in a future work, where we will try to develop a theoretical framework for one-stage continuous melting transitions.

Let us also emphasize that the models we have analyzed are not the best choice to explore with high numerical resolution the possibility of one-stage melting from a general perspective. Effective models,with a dynamics at the level of clusters (as opposed to having a microscopic dynamics for particles) would be a much better option from the computational point of view.

In summary, despite of interesting and not fully understood differences, the striking patterns produced by the biologically inspired Langevin equation (II) resemble very much the melting/freezing solid/liquid transition in equilibrium systems.

It is our hope that this paper will motivate further studies of (i) the effect of fluctuations on self-organized nonequilibrium patterns and (ii) the analogies and differences between the defect-mediated type of melting transition described here and standard equilibrium melting scenarios.

ACKNOWLEDGMENTS We acknowledge useful comments and discussions with F. de los Santos. This work was supported in part by the Spanish projects FIS2005-00791 and FIS2004-00953 (Ministerio de Educación y Ciencia), FQM-165 (Junta de Andalucía), and the European Commission through the NEST-Complexity project PATRES (043268).