Spatiotemporal Self-Organization of Fluctuating Bacterial Colonies

TL;DR

This paper models bacterial colonies' spatiotemporal dynamics, revealing how fluctuations influence colony stability and movement, and connects these behaviors to biofilm life cycles using bifurcation and large deviation analysis.

Contribution

It introduces a coupled model of motility-induced phase separation and population dynamics, highlighting the effects of fluctuations on colony metastability and spatial organization.

Findings

Static colonies become metastable with fluctuations.

Limit cycles lead to random colony formation locations.

The model explains aspects of biofilm and planktonic life cycles.

Abstract

We model an enclosed system of bacteria, whose motility-induced phase separation is coupled to slow population dynamics. Without noise, the system shows both static phase separation and a limit cycle, in which a rising global population causes a dense bacterial colony to form, which then declines by local cell death, before dispersing to re-initiate the cycle. Adding fluctuations, we find that static colonies are now metastable, moving between spatial locations via rare and strongly nonequilibrium pathways, whereas the limit cycle becomes almost periodic such that after each redispersion event the next colony forms in a random location. These results, which hint at some aspects of the biofilm-planktonic life cycle, can be explained by combining tools from large deviation theory with a bifurcation analysis in which the global population density plays the role of control parameter.

Equations

$${\partial }_t\rho =\nabla \cdot ({\mathcal{D}}_e(\rho )\nabla \rho )-{\delta }^2\nabla \cdot \left(\rho D(\rho )\nabla \Delta \rho \right).$$

$$E[\rho ]={\int }_{\Omega }\left(\rho \log \rho -\rho +f(\rho )+\frac{1}{2}{\delta }^2|\nabla \rho {|}^2\right)d\mathbf{x}.$$

$${\partial }_t\rho$$

$$I_T[\rho ]=\frac{1}{2}{\int }_0^{T}dt{\int }_{\Omega }d\mathbf{x}\theta (\mathbf{x},t;\rho ){\partial }_t\rho (\mathbf{x},t)$$

$${\partial }_t\rho$$

$$+\alpha \rho e^{\theta }-\alpha {\rho }^2{e}^{-\theta }/{\rho }_0.$$

Full text

Because they are not bound by the standard laws of equilibrium thermodynamics, active materials such as bird flocks Ballerini et al. (2008), motile bacteria Marchetti et al. (2013), self-organizing bio-polymers Schaller et al. (2010); Sumino et al. (2012), or man-made self-propelled particles Howse et al. (2007) have many more routes towards self-assembly and self-organization than systems whose dynamics satisfy detailed-balance. Motility-induced phase separation (MIPS) is one example Wittkowski et al. (2014); Cates and Tailleur (2015). MIPS arises naturally in systems of self-propelled particles whose locomotive speed decreases monotonically with density, through a feedback in which particles accumulate where they move slowly and vice-versa. This provides a generic path to pattern formation in motile agents, including both micro-organisms such as E. coli Shapiro (1995); Tailleur and Cates (2008), where slowdown can be caused by quorum sensing Parsek and Greenberg (2005); Liu et al. (2011); Fu et al. (2012); Solano et al. (2014), and synthetic colloidal analogues where slowdown is caused by crowding Fily and Marchetti (2012); Palacci et al. (2013); Buttinoni et al. (2013); Cates and Tailleur (2013, 2015). The simplest theories of MIPS describe phase separation through an effective free energy functional Cates and Tailleur (2015), although active gradient terms can alter this structure and the resulting coexistence behavior Wittkowski et al. (2014); Cates and Tailleur (2015).

MIPS, can be arrested by the birth and death of particles. The simplest (logistic) population dynamics has for uniform systems a stable fixed point at some carrying capacity $\rho_{0}$, with decay towards this from higher or lower densities. If $\rho_{0}$ lies within the miscibility gap of a phase separation, the uniform state at $\rho_{0}$ is unstable, but so is a state of coexisting bulk phases (each of which would have nonstationary density). This can result in stationary patterns of finite wavelength, in which bacteria reproduce in dilute regions, migrate by diffusive motility to dense ones, become immotile, and there die off by overcrowding Cates et al. (2010); Cates and Tailleur (2015). In Cates et al. (2010) a field theoretic description of such arrested phase separation was proposed to capture the variety of bacterial colony patterns earlier observed in experiments Murray (2003).

The work of Cates et al. (2010) mainly describes stationary spatial rather than spatiotemporal patterns, and crucially neglects the effects of intrinsic fluctuations in population density and motility of the bacteria. These fluctuations are always present in finite systems and they are relevant to the nature and stability of the self-assembling structures. One purpose of the present Letter is to explore these effects in situations where birth and death processes are slow compared to diffusive exploration times. Another is to show that entirely different physics from that reported in Cates et al. (2010) can arise, involving spatiotemporal rather than static patterning. For simplicity we focus on one-dimensional systems whose extent is smaller than the natural length scale for the spacing between dense patches, controlled by the distance a particle can move during its lifetime. Accordingly, only one dense bacterial domain of low motility (hereafter ‘colony’) is present at a time, with the remaining bacteria in a dilute and motile (i.e., ‘planktonic’) phase. These simplifications allow us to focus on temporal aspects of the pattern formation, and also to give a thorough analysis of the role of noise terms. Our qualitative conclusions, however, apply in higher dimensions.

Concretely, we identify two regimes in which the fluctuations induce nontrivial self-organization pathways. The first arises in situations where the deterministic dynamics possess time-periodic solutions driven by the interplay between MIPS and the competition for resources: in the resulting limit cycles a colony of bacteria periodically appears and disappears at a fixed location in the domain. In this regime, the fluctuations, no matter how small, are shown to have a drastic impact: they allow the colony to explore its environment by randomly jumping from one location to another each time the system revisits its spatially homogeneous planktonic state. With periodic spatial boundary conditions the colony now appears and disappears at spatially random locations; if boundary conditions instead favor localization at container walls, a random choice of wall is made each cycle. It is striking that this behavior, which resembles the biofilm-planktonic lifecycle of several micro-organisms Kostakioti et al. (2013) can arise from the combination of MIPS and logistic growth alone – particularly as the bacterial quorum sensing response, a likely cause of MIPS Tailleur and Cates (2008); Cates et al. (2010) is now thought also to be directly implicated in biofilm dispersion Solano et al. (2014). We are not suggesting that the combination of MIPS and logistic growth directly explains the complex life cycle of real biofilms, but it may nonetheless be a contributory factor in its evolution from simpler beginnings.

In the second regime, we show that the deterministic theory predicts a single static colony with multiple stable locations. In these situations the intrinsic fluctuations induce rare, noise-activated transitions via the uniform state from one such metastable pattern to another. (This again resembles the biofilm lifecycle, but now with random intervals between dispersion events.) As shown below, the rate and mechanism of these transitions can be completely characterized using tools from large deviation theory (LDT) Freidlin and Wentzell (2012). The transition pathways are intrinsically out-of-equilibrium and involve a subtle interplay between phase-separation and reproduction. In particular these transitions do not follow the deterministic relaxation path in reverse, as one would predict from standard equilibrium arguments in systems with microscopic reversibility.

Following Cates et al. (2010), we focus on motile bacteria whose motile diffusivity $D(\rho)=\tfrac{1}{2}v^{2}(\rho)\tau/d$ Tailleur and Cates (2008) depends on their local density $\rho$ via the swim speed $v(\rho)$. Here $\tau$ is a rotational relaxation time and $d$ is dimensionality. Because the swim speed varies in space through $\rho$, the particles have a mean drift velocity $\boldsymbol{V}(\rho)=-\tfrac{1}{2}D^{\prime}(\rho)\nabla\rho$ Tailleur and Cates (2008). Neglecting fluctuations, these effects can be combined into the following equation governing the evolution of the bacterial density $\rho({\boldsymbol{x}},t)$:

[TABLE]

In this equation the diffusivity and drift have been combined into a collective diffusivity $\mathcal{D}_{e}(\rho)=D(\rho)+\tfrac{1}{2}\rho\,D^{\prime}(\rho)$, and a higher-order gradient, or regularizing, term proportional to $\delta^{2}$ has been added to account for the fact that the bacteria only sense each other’s influence over a finite distance $\delta>0~{}$111In general this effect leads to a nonlocal equation, out of which (1) emerges at leading order after expansion in $\delta\ll 1$ Grafke et al. (2017). .

Our choice of regularizer is the simplest form to emerge from explicit coarse-graining of the microscopic dynamics of active particles Tailleur and Cates (2008); Speck et al. (2014); Cates and Tailleur (2015). With this form, and in the absence of birth-death processes, the physics of MIPS, including all noise contributions, maps exactly onto an equilibrium model of phase separation. The regularizer then derives from a square gradient contribution in an underlying free energy that describes passive Brownian particles with attractive interactions Cates and Tailleur (2015). (A more detailed coarse-graining gives further gradient terms that are not representable by any free energy; these shift the MIPS phase boundaries only slightly Cates and Tailleur (2015); Stenhammar et al. (2013).) Thus our model, whose microscopic derivation is in the supplemental material 137137137See supplemental material at [URL will be inserted by publisher] for a derivation from a microscopic model, includes a fluctuating noise term in Eq. (1) such that its dynamics obey detailed balance (DB) with the free energy

[TABLE]

Eq. (1) then takes the form of a generalized gradient flow, $\partial_{t}\rho=-M(\rho)(\delta F/\delta\rho)$, with nonlinear mobility operator $M(\rho)\xi=\nabla\cdot(\rho D(\rho)\nabla\xi)$ Mielke et al. (2016). The DB property of Eq. (1) is inessential since it is violated by the birth and death terms added below 133713371337Detailed balance violation arises because the rates for birth and death are unaffected by the motility induced “attractions” (whose true origin is kinetic) whereas genuine enthalpic attractions would require the birth and death rates to become explicitly density dependent Lefever et al. (1995); Glotzer et al. (1995a, 1994, b). but it simplifies the analysis when these are small. The regularizer used here is an improvement on the one adopted in Cates et al. (2010), which is purely phenomenological and does not emerge from any known coarse-graining of an active-particle model.

Crucially, whenever $d\ln v/d\ln\rho<-1$, we have $\mathcal{D}_{e}(\rho)<0$ in (1). This is the spinodal regime of local instability for MIPS Tailleur and Cates (2008). We assume $v(\rho)=v_{0}e^{-\lambda\rho/2}$ where $v_{0}$ is the speed of an isolated particle and $\lambda>0$ Cates et al. (2010). In a $d$-dimensional box $\Omega=[0,L]^{d}$, after non-dimensionalization via $\tau v_{0}^{2}=\lambda=1$, (1) reduces to

[TABLE]

We study (3) for ${\boldsymbol{x}}\in\Omega=[0,1]^{d}$ with, unless otherwise stated, Neumann boundary conditions: $\hat{\boldsymbol{n}}\cdot\nabla\rho=0\,\forall\,{\boldsymbol{x}}\in\partial\Omega$ for $\hat{\boldsymbol{n}}$ normal to the boundary. These localize the dense phase at a wall, to minimize its interfacial energy.

The mean bacterial density, $\bar{\rho}=|\Omega|^{-1}\int_{\Omega}\rho(x)\,d{\boldsymbol{x}}$, is conserved by (3) and controls the phase separation. Fig. 1 (left) shows the bifurcation diagram obtained when $d=1$ and $\delta^{2}=2\cdot 10^{-3}$ by projecting the fixed points of (3) onto the $(\bar{\rho},A)$ plane, where the ‘signed asymmetry’ is defined as $A[\rho]=\smash{\int_{0}^{1/2}\rho(x)\,dx}-\smash{\int_{1/2}^{1}\rho(x)\,dx}$. A linear stability analysis finds the spatially uniform solution to be stable if $\bar{\rho}<\rho_{\text{S}}=2/(1-2\delta^{2}\pi^{2})$. (This becomes the bulk spinodal condition in the large system, or small $\delta$, limit.) A subcritical pitchfork bifurcation occurs at $\bar{\rho}=\rho_{\text{S}}$, where two unstable and one stable branches merge into a single unstable branch. The remaining two stable branches, $\rho_{\text{L}}(x)$ and $\rho_{\text{R}}(x)$, correspond to colony formation on either the left or the right wall. Once present, each such phase-separated state remains stable down to a ‘binodal’ density $\bar{\rho}=\rho_{\text{B}}<\rho_{\text{S}}$ 133813381338For the particular choice in (3), $\rho_{\text{B}}\to 0$ in the large system limit but for a finite system the lower stability limit $\rho_{\text{B}}<\rho_{\text{S}}$ remains nonzero, albeit now depending on $\bar{\rho}$ as well as $\delta$., whereas for $\bar{\rho}<\rho_{\text{B}}$ the colony redisperses diffusively.

These transitions in static stability become dynamically significant once the mean bacterial density $\bar{\rho}$ is allowed to change by introducing logistic growth. Focusing again on the deterministic situation first, this adds to (3) a term $\alpha\rho(1-\rho/\rho_{0})$, where $\alpha$ is the birth rate and $\rho_{0}$ is the carrying capacity. Crucially, this breaks detailed balance, allowing not only steady-state fluxes Cates et al. (2010) but also nonstationary patterns (see below). When the population dynamics is much slower than particle motion, $\alpha\ll 1$, the bifurcation diagram shown in Fig. 1 now depicts the projection of a slow manifold $\mathcal{M}$ instead. Diffusive motility leads to convergence of $\rho$ to $\mathcal{M}$ on a fast time-scale $O(1)$ in $\alpha$, during which the global mean $\bar{\rho}$ is almost conserved. On $\mathcal{M}$ itself, the motion is driven solely by changing $\bar{\rho}$ via the slow logistic term on time-scales $O(\alpha^{-1})$.

Three different regimes can then be identified on varying the carrying capacity $\rho_{0}$. For $\rho_{0}<\rho_{\text{S}}$ the system remains in a uniform phase, in which dilute motile bacteria homogeneously cover the domain at $\rho(x)=\rho_{0}$ 133913391339Note that for $\rho_{\text{B}}<\rho_{0}<\rho_{\text{S}}$, the two phase separated states $\rho_{\text{L,R}}$ are destabilized by the $\alpha$ term.. At higher capacities homogeneous solutions cannot exist, and the solution jumps to one of the stable branches, leading to a dense layer of immotile bacteria at one of the boundaries, which forms on the fast time-scale. However, for carrying capacities $\rho_{\text{S}}<\rho_{0}<\rho_{0}^{c}$, the resulting colony is destabilized by the birth/death term; two limit cycles then appear. As soon as a colony emerges at one end of the domain, the bacteria in it start to slowly die out until the global density drops again to $\rho_{\text{B}}$. It then disperses rapidly, reverting to a uniform phase, whose density slowly grows until $\rho_{\text{S}}$ is reached and the cycle begins anew. Symmetry is broken by weak memory of the previous cycle, so that the colony always reforms in the same place. (This also holds with periodic boundary conditions, but the location is then arbitrary.) One of these limit cycles is projected onto the bifurcation diagram in Fig 1 (left), and its spatio-temporal evolution shown in Fig. 1 (right, a).

Finally, for $\rho_{0}>\rho_{0}^{c}$, all stable fixed points of the system comprise a dense colony in coexistence with a planktonic ‘vapor’. A steady flux of particles from the vapor balances cell-death within the colony, so that the macroscopic model is stationary, while the microscopic dynamics are not Cates et al. (2010). Such fixed points correspond to points on stable branches of the slow manifold $\mathcal{M}$. For $d=1$ there are two such branches ($\rho_{\text{L,R}}$), with more in $d\geq 2$, corresponding to (say) colonies located in the corners of a square in $d=2$. For these computations and the ones below, we picked $\alpha=10^{-4}$ and $\delta^{2}=2\cdot 10^{-3}$, which leads to $\rho_{\text{S}}\approx 2.08$, $\rho_{\text{B}}\approx 1.64$, and $\rho_{0}^{c}\approx 6.749$. The complete phase diagram in $\delta$ and $\rho_{0}$ is shown in Fig. 2.

So far we have neglected the effect of intrinsic fluctuations, both in the diffusive and the reproductive dynamics. The fluctuations in diffusion can be formally accounted for by adding in (3) a noise-term $\sqrt{N^{-1}}\nabla\cdot(\sqrt{2\rho D(\rho)}{\boldsymbol{\eta}})$, where ${\boldsymbol{\eta}}$ is spatio-temporal white noise and $N$ is the typical number of particles present in the domain $\Omega$. On the other hand, the population dynamics can be modeled by a reversible reaction $A\leftrightarrow A+A$ with forward and backward rates $r_{f}=\alpha$ and $r_{b}=\alpha/(N\rho_{0})$ Note (1337). The combined effect is captured by a Markov jump process. In the limit $N\to\infty$, this gives the logistic growth term $\alpha\rho(1-\rho/\rho_{0})$ considered before, with fluctuations that are Poisson at each location, scaling again with $\sqrt{N^{-1}}$ Shwartz and Weiss (1995).

We consider only the weak fluctuation (large $N$) regime which can be captured by LDT Freidlin and Wentzell (2012). This theory predicts that a noise-driven event occurs, with probability close to $1$, via the path involving the least unlikely fluctuation able to drive this event. The resulting most likely path (MLP, also referred to as the instanton) is the minimizer of the action functional

[TABLE]

which quantifies the likelihood of the ‘fluctuation’ $\theta({\boldsymbol{x}},t;\rho)$. In our context, $\theta({\boldsymbol{x}},t;\rho)$ is related to $\rho({\boldsymbol{x}},t)$ Note (137) via

[TABLE]

The action $I_{T}$ should be minimized over all paths, and all durations $T$, consistent with the event of interest, giving $I^{*}=\inf_{T>0}\inf_{\rho}I_{T}[\rho]$. The minimizer is the MLP for the event, whose rate is then, up to a prefactor, $\exp(-NI^{*})$  134013401340The numerical details of the minimization are discussed in Heymann and Vanden-Eijnden (2008); Grafke et al. (2016)..

In the regime $\rho_{0}<\rho_{\text{S}}$, the spatially homogeneous configuration is stable and LDT predicts that deviations away from it are exponentially rare and transient. This is no longer true for $\rho_{\text{S}}<\rho_{0}<\rho_{0}^{c}$: Since part of the limit cycle tracks $O(\alpha)$ close to the separatrix for times $O(\alpha^{-1})$, even tiny fluctuations can trigger a crossing of the separatrix into the other limit cycle. This is consistent with LDT, for which in the limit $\alpha\to 0$ zero-action minimizers connect these two limit cycles. Therefore, due to the fluctuations, bacterial colonies randomly appear either on the left or right wall, disappearing again after times $O(\alpha^{-1})$. This switching behavior is depicted in Fig. 1 (b). With periodic boundary conditions, the colony instead appears at a random location each cycle, as shown in Fig. 1 (c). Note that the time period of the cycles is not affected significantly by these fluctuations.

For $\rho_{0}\,>\,\rho_{0}^{c}$, on the other hand, $\rho_{\text{L}}$ and $\rho_{\text{R}}$ become metastable; the noise triggers rare and aperiodic transitions between them. A projection of the resulting MLP from $\rho_{\text{L}}$ to $\rho_{\text{R}}$ onto the bifurcation diagram is shown in Fig. 3 (left), whereas its actual shape is shown in Fig. 3 (center); that from $\rho_{\text{R}}$ to $\rho_{\text{L}}$ follows by symmetry. To understand its features, notice that the slow manifold connects the two stable fixed points $\rho_{\text{L}}$ and $\rho_{\text{R}}$ through the bifurcation point, and so can be used as a channel to facilitate the transition. This is indeed correct for $\delta\ll 1$, as the free energy barrier for a jump between stable branches scales like $\delta^{-1}$ Grafke et al. (2017). It is confirmed by a numerical calculation of the MLP as shown in Fig. 3 (center): The colony of bacteria on the left first disperses to form a uniform planktonic phase which then attains the bifurcation point $\rho_{\text{X}}({\boldsymbol{x}})\,=\,\rho_{\text{S}}$. The system then follows the unstable branch (separatrix) of $\mathcal{M}$ – with two symmetric colonies – to the saddle point at $\rho_{\text{Saddle}}({\boldsymbol{x}})$, where it enters the basin of attraction of $\rho_{\text{R}}$. Note that for finite $N$, diffusive noise will push the actual trajectory off the separatrix well before it reaches the transition state at $\rho_{\text{Saddle}}$. The event rate is still found by LDT, since the motion after visiting $\rho_{\text{X}}$ is effectively deterministic and does not contribute to $I^{*}$.

Let us stress the non-equilibrium nature of this transition: If the system were in detailed balance, time-reversal symmetry would require the transition path to follow the deterministic relaxation trajectory in reverse from $\rho_{\text{L}}$ to the transition state $\rho_{\text{Saddle}}$, and then relax deterministically to $\rho_{\text{R}}$ (see Fig. 3). Along this trajectory, the colony initially grows instead of shrinking, bringing the global average density $\bar{\rho}$ up to $\rho_{\text{Saddle}}$. Subsequently, bacteria leave the colony, cross the low-density region in the center, and accumulate on the opposite wall, keeping $\bar{\rho}$ constant. This part of the transition happens diffusively causing the system to leave the slow manifold. After passing through $\rho_{\text{Saddle}}$, the reversible transition necessarily coincides with the true MLP.

In summary, we have analyzed a minimal model for the fluctuating dynamics of self-organization among motile bacteria in a finite domain, in the presence of birth and death processes that are slow compared to diffusive time-scales. For $\rho_{\text{S}}<\rho_{0}<\rho_{0}^{c}$, a bacterial colony is present; fluctuations allow it to explore the domain by random relocation at regular intervals. This exploration remains possible at higher carrying capacities, $\rho_{0}>\rho_{0}^{c}$, where fluctuations now induce exponentially rare, aperiodic relocations. Such transitions are inherently out-of-equilibrium: their most likely path differs significantly from the one followed in systems with time-reversible dynamics. In our model, relocation of a colony is triggered by a slow decline in local population density, followed by rapid dispersal and re-formation elsewhere, rather than by progressive migration between old and new sites. The intermediate state has no colonies in the first case and two in the second – alternatives that should be clearly distinguishable experimentally. Similar considerations may apply to other organisms that alternate between nomadic and cooperative lifestyles. The formalism readily applies to a higher-dimensional setup, e.g. bacterial colonies in the corners of a two-dimensional rectangular domain, and a slow manifold with four stable branches.

Acknowledgments. We thank T. Schäfer for his help with the numerical scheme and A. Donev, O. Hirschberg, and C. Nardini for interesting discussions. MEC is funded by the Royal Society. EVE is supported in part by the Materials Research Science and Engineering Center (MRSEC) program of the National Science Foundation (NSF) under award number DMR-1420073 and by NSF under award number DMS-1522767.